SUBROUTINE SetupCompDomain_test
! By Allan P. Engsig-Karup.
USE GlobalVariables
IMPLICIT NONE

REAL(KIND=long), DIMENSION(:,:), POINTER :: x, y
REAL(KIND=long), DIMENSION(:), POINTER   :: z
INTEGER :: Nx, Ny, Nz
! Local variable...
REAL(KIND=long), DIMENSION(:), POINTER :: R, angle ! Radius and angle...
REAL(KIND=long), DIMENSION(:,:), POINTER :: e, n
INTEGER :: i, j, k, Nxstart, Nxcircle

Nx = FineGrid%Nx
Ny = FineGrid%Ny
Nz = FineGrid%Nz

! Allocate space for physical grids and ghost point layers
ALLOCATE(FineGrid%x(Nx+2*GhostGridX,Ny+2*GhostGridY))
ALLOCATE(FineGrid%y(Nx+2*GhostGridX,Ny+2*GhostGridY))
ALLOCATE(FineGrid%z(Nz+GhostGridZ))

! TEMPORARY POINTERS/VARIABLES
x => FineGrid%x
y => FineGrid%y
z => FineGrid%z

! ALLOCATE BATHYMETRY DATA
ALLOCATE(FineGrid%h(Nx+2*GhostGridX,Ny+2*GhostGridY))
IF (Nx>1) THEN
  ALLOCATE(FineGrid%hx(Nx+2*GhostGridX,Ny+2*GhostGridY))
  ALLOCATE(FineGrid%hxx(Nx+2*GhostGridX,Ny+2*GhostGridY))
ENDIF
IF (Ny>1) THEN
  ALLOCATE(FineGrid%hy(Nx+2*GhostGridX,Ny+2*GhostGridY))
  ALLOCATE(FineGrid%hyy(Nx+2*GhostGridX,Ny+2*GhostGridY))
ENDIF

!$$$$$$ ! Semi-circular channel... (with a zone to generate and a zone to absorb...)
!$$$$$$ ! FIXME it is probably not the correct way to do this... (define grid as before, do this in the initialization...)
!$$$$$$ IF ((Nx>1).AND.(Ny>1)) THEN
!$$$$$$
!$$$$$$     ! Index of the beginning of the circle
!$$$$$$     Nxstart = Nx/4
!$$$$$$     ! Index of the end of the circle
!$$$$$$     !Nxstop = NINT(2*Nx/3)
!$$$$$$     ! Number of points in the circle part
!$$$$$$     Nxcircle = Nx/2
!$$$$$$
!$$$$$$     ! Allocate the temporary radius and angle...
!$$$$$$     ALLOCATE(angle(1:Nxcircle))
!$$$$$$     ALLOCATE(R(-GhostGridY:Ny-1+GhostGridY))
!$$$$$$     !
!$$$$$$     DO i = 1,Nxcircle
!$$$$$$         angle(i) = PI-(REAL(i-1,long)/REAL((Nxcircle-1),long)*PI)
!$$$$$$     ENDDO
!$$$$$$     DO j = -GhostGridY,Ny-1+GhostGridY
!$$$$$$         R(j) = Lx + REAL(j,long)/REAL((Ny-1),long)*(Ly)
!$$$$$$     END DO
!$$$$$$     ! Define fine grid
!$$$$$$     DO i=1,Nxcircle
!$$$$$$         DO j=-GhostGridY,Ny-1+GhostGridY
!$$$$$$             x(i+Nxstart-1,j+1+GhostGridY) = R(j)*COS(angle(i))
!$$$$$$             y(i+Nxstart-1,j+1+GhostGridY) = R(j)*SIN(angle(i))
!$$$$$$         ENDDO
!$$$$$$     ENDDO
!$$$$$$     !
!$$$$$$     DEALLOCATE(angle, R)
!$$$$$$     !
!$$$$$$     DO i=1,Nxstart-1
!$$$$$$       DO j=1,Ny+2*GhostGridY
!$$$$$$         x(i,j)=x(Nxstart,j)
!$$$$$$         y(i,j)=y(Nxstart,j)-REAL(Nxstart-i,long)/REAL(Nxstart-1-GhostGridX,long)*2*SFsol%L!4*Lx !FIXME approx...
!$$$$$$       ENDDO
!$$$$$$     ENDDO
!$$$$$$     !
!$$$$$$     DO i=Nxstart+Nxcircle,Nx+2*GhostGridX
!$$$$$$       DO j=1,Ny+2*GhostGridY
!$$$$$$         x(i,j)= x(Nxstart+Nxcircle-1,j)
!$$$$$$         y(i,j)= y(Nxstart+Nxcircle-1,j)-REAL(i-(Nxstart+Nxcircle)+1,long)/&
!$$$$$$             REAL(Nx+2*GhostGridX-(Nxstart+Nxcircle)+1-GhostGridX,long)*2*SFsol%L!4*Lx !FIXME approx...*4*Lx
!$$$$$$       ENDDO
!$$$$$$     ENDDO
!$$$$$$ ELSE
!$$$$$$   PRINT *,'Error this is not done...'
!$$$$$$   STOP
!$$$$$$ ENDIF
! Cylinder case...
! Define the radius and angle
IF ((Nx>1).AND.(Ny>1)) THEN
    ! Allocate the temporary radius and angle...
    ALLOCATE(angle(-GhostGridX:Nx-1+GhostGridX))
    ALLOCATE(R(-GhostGridY:Ny-1+GhostGridY))
    !
    DO i = -GhostGridX,Nx-1+GhostGridX
        angle(i) = PI-(REAL(i,long)/REAL((Nx-1),long)*PI) !PI-0.5d0/REAL((Nx-1),long)*PI-(REAL(i,long)/REAL((Nx-1),long)*PI) !
!$$$$$$         print*,'Initial',angle(i)
!$$$$$$         angle(i) = (PI+(COS(PI*REAL(i,long)/REAL((Nx-1),long)))*PI)/two
!$$$$$$         angle(i) = (PI+TANH(1.5_long*REAL(Nx-1-2*i,long)/REAL(Nx-1))/TANH(1.5_long)*PI)/two
    ENDDO
    i=0
    DO j = -GhostGridY,Ny-1+GhostGridY
      i=i+1
        R(j) = Lx + REAL(j,long)/REAL((Ny-1),long)*(Ly)
        !angle(j) = REAL(j,long)/REAL((Ny-1),long)*PI
        ! cosine grid (seems too strng)
        !R(j) = Lx + (COS(pi*REAL(j,long)/REAL(2*(Ny-1),long)))*(Ly)
        ! stretch grid
!$$$$$$         R(Ny-1+GhostGridY+1-i) = Lx + (1.0_long-TANH(1.3_long*REAL(j,long)/REAL(Ny-1))/TANH(1.3_long))*(Ly)
    END DO
    ! Define symetrical points...
    angle(-1)=angle(0)-(angle(1)-angle(0))
    angle(Nx-1+GhostGridX)=angle(Nx-1)-(angle(Nx-2)-angle(Nx-1))
    !
    R(-1)=R(0)-(R(1)-R(0))
    R(Ny-1+GhostGridY)=R(Ny-1)-(R(Ny-2)-R(Ny-1))
    ! Define fine grid
    DO i=-GhostGridX,Nx-1+GhostGridX
        DO j=-GhostGridY,Ny-1+GhostGridY
            x(i+1+GhostGridX,j+1+GhostGridY) = R(j)*COS(angle(i))
            y(i+1+GhostGridX,j+1+GhostGridY) = R(j)*SIN(angle(i))
        ENDDO
    ENDDO
    ! Prevent from symetrical ghost points...
    !x(1,:)=x(2,:)
    !y(1,:)=y(2,:)-(y(3,2)-y(2,2))
    !x(Nx+2*GhostGridX,:)=x(Nx+GhostGridX,:)
    !y(Nx+2*GhostGridX,:)=y(Nx+GhostGridX,:)-(y(Nx,2)-y(Nx+GhostGridX,2))
    !
    DEALLOCATE(angle, R)
ELSE
  PRINT *,'Error this is not done...'
  STOP
ENDIF
! Grid with straight boundaries
!$$$$$$ ! Define fine grid
!$$$$$$ IF (Nx==1) THEN
!$$$$$$     x(1,1:Ny+2*GhostGridY) = Lx;
!$$$$$$ ELSE
!$$$$$$     SELECT CASE (GridX)
!$$$$$$         CASE (0) ! Even Grid
!$$$$$$             DO i = -GhostGridX,Nx-1+GhostGridX
!$$$$$$                x(i+1+GhostGridX,1:Ny+2*GhostGridY) = REAL(i,long)/REAL((Nx-1),long)*Lx;
!$$$$$$             END DO
!$$$$$$         CASE DEFAULT
!$$$$$$             PRINT *,'Error: Specified grid for x-direction invalid.'
!$$$$$$             STOP
!$$$$$$     END SELECT
!$$$$$$ ENDIF
!$$$$$$ IF (Ny==1) THEN
!$$$$$$     y(1:Nx+2*GhostGridX,1) = Ly;
!$$$$$$ ELSE
!$$$$$$     SELECT CASE (GridY)
!$$$$$$         CASE (0) ! Even Grid
!$$$$$$             DO i = -GhostGridY,Ny-1+GhostGridY
!$$$$$$                y(1:Nx+2*GhostGridX,i+1+GhostGridY) = REAL(i,long)/REAL((Ny-1),long)*Ly;
!$$$$$$             END DO
!$$$$$$         CASE DEFAULT
!$$$$$$             PRINT *,'Error: Specified grid for y-direction invalid.'
!$$$$$$             STOP
!$$$$$$     END SELECT
!$$$$$$ END IF
! Define fine grid
SELECT CASE (GridZ)
	CASE (0) ! Even Grid
		DO i = 0,Nz-1
		   z(i+1+GhostGridZ) = REAL(i,long)/REAL((Nz-1),long)*Lz;
		END DO
	CASE (1) ! Uneven Grid
		DO i = 0,Nz-1
			z(i+1+GhostGridZ) = SIN(pi*REAL(i,long)/REAL(2*(Nz-1),long))
		END DO
	CASE DEFAULT
		PRINT *,'Error: Specified grid for z-direction invalid.'
		STOP
END SELECT
IF (GhostGridZ==1) THEN
	! include ghost point below bottom. Reflect point just above bottom symmetrically
	! to ghost point by reflection about z=0
	z(1)=-z(3)
ENDIF
! Allocate array for coefficients involving sigma and its derivatives
! to be used in defining the transformed laplace problem.
!ALLOCATE(FineGrid%dsigma((Nx+2*GhostGridX)*(Ny+2*GhostGridY)*(Nz+GhostGridZ),5))
!FineGrid%dsigma=zero
END SUBROUTINE SetupCompDomain_test


!> @brief   functions to cast different unlimited polymorphic variables to different types.
!> @author  Robert schuster
module fplus_converter
    use fplus_fillvalue
    implicit none
    public

contains

    function to_integer(value)
        integer :: to_integer
        class(*) :: value
        select type (value)
            type is (integer)
                to_integer = value
            class default
                to_integer = fplus_fill_int
        end select
    end function

    function to_real(value)
        real :: to_real
        class(*) :: value
        select type (value)
            type is (real (kind=4))
                to_real = value
            class default
                to_real = fplus_fill_real
        end select
    end function

    function to_integer8(value)
        integer (kind=8) :: to_integer8
        class(*) :: value
        select type (value)
            type is (integer (kind=8))
                to_integer8 = value
            class default
                to_integer8 = fplus_fill_intK8
        end select
    end function

    function to_real8(value)
        real (kind=8) :: to_real8
        class(*) :: value
        select type (value)
            type is (real (kind=8))
                to_real8 = value
            class default
                to_real8 = fplus_fill_realK8
        end select
    end function

end module fplus_converter


text:
      subroutine Setup_NodeElems()
!***********************************************************************
! Copyright 2011 Los Alamos National Security, LLC  All rights reserved
! Unless otherwise indicated,  this information has been authored by an
! employee or employees of the Los Alamos National Security, LLC (LANS),
! operator of the  Los  Alamos National  Laboratory  under Contract  No.
! DE-AC52-06NA25396  with  the U. S. Department  of  Energy.  The  U. S.
! Government   has   rights  to  use,  reproduce,  and  distribute  this
! information.  The  public may copy  and  use this  information without
! charge, provided that this  Notice and any statement of authorship are
! reproduced on all copies.  Neither  the  Government nor LANS makes any
! warranty,   express   or   implied,   or   assumes  any  liability  or
! responsibility for the use of this information.
!***********************************************************************
!
! store elements connected to each node
!
! Author : Sai Rapaka

      use comai, only: nei, neq, ns
      use comfem

      implicit none

      integer, dimension(neq)        :: elcount
      integer node, nsize
      integer i_start, i_end
      logical i_done
      integer i, j, k

      ! Currently written only for hex elements in 3-D
      ! Easy to add conditions for other elements
      elcount = 0
      do i=1,nei
        do j=1,ns
          node = elnode(i,j)
          elcount(node) = elcount(node) + 1
        enddo
      enddo

      ! Compute size of NodeElems
      nsize = neq + 1
      do i=1,neq
        nsize = nsize + elcount(i)
      enddo

      allocate(NodeElems(nsize))
      NodeElems = 0

      do i=1,neq+1
        if(i.eq.1) then
          NodeElems(i) = (neq+1)
        else
          NodeElems(i) = NodeElems(i-1) + elcount(i-1)
        endif
      enddo

      do i=1,nei
        do j=1,ns
          node = elnode(i,j)
          i_start = NodeElems(node) + 1
          i_end = NodeElems(node+1)
          i_done = .false.
          do k=i_start,i_end
            if((NodeElems(k).eq.0).and.(.not.i_done)) then
              NodeElems(k) = i
              i_done = .true.
            endif
          enddo
        enddo
      enddo

      end subroutine Setup_NodeElems

ext: f

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


!----------------------------------------------------------------------------------------------
! MODULE: overlap_mo_mod
!> @author Marin Sapunar, Ruđer Bošković Institute
!> @date March, 2019
!
!> @brief Subroutine for calculating overlap matrix between two sets of molecular orbitals..
!----------------------------------------------------------------------------------------------
module overlap_mo_mod
    use global_defs
    implicit none

    private
    public overlap_mo

contains

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: overlap_mo
    !> @brief Calculate overlap matrix between two sets of molecular orbitals.
    !----------------------------------------------------------------------------------------------
    subroutine overlap_mo()
        use overlap_variables
        use matrix_mod, only : mat_ge_mmm
        use read_all_mod
        logical :: uhf1, uhf2
        integer :: i

        if (allocated(s_mo_a)) deallocate(s_mo_a)
        if (allocated(s_mo_b)) deallocate(s_mo_b)

        ! Read MOs.
        time0 =  omp_get_wtime()
        call read_mo(input_format_1, path1, mos1, bs1)
        call read_mo(input_format_2, path2, mos2, bs2)

        ! Check unrestricted/restricted calculations.
        uhf1 = (mos1%n_mo_b > 0)
        uhf2 = (mos2%n_mo_b > 0)
        if (uhf1 .and. (.not. uhf2)) then
            if (print_level >= 1) then
                write(stdout, *)
                write(stdout, '(1x,a)') 'Unrestricted bra and restricted ket MOs...'
            end if
            call mos2%to_unrestricted()
        else if (uhf2 .and. (.not. uhf1)) then
            if (print_level >= 1) then
                write(stdout, *)
                write(stdout, '(1x,a)') 'Restricted bra and unrestricted ket MOs...'
            end if
            call mos1%to_unrestricted()
        end if
        time_in =  time_in + omp_get_wtime() - time0

        ! Calculate overlaps.
        time0 = omp_get_wtime()
        if (print_level >= 2) then
            write(stdout, *)
            write(stdout, '(1x,a)') 'Computing MO overlaps...'
            write(stdout, '(5x,a,2(1x,i0))') 'Number of bra of AOs/MOs:', shape(mos1%ca)
            write(stdout, '(5x,a,2(1x,i0))') 'Number of ket of AOs/MOs:', shape(mos2%ca)
        end if
        allocate(s_mo_a(mos1%n_mo_a, mos2%n_mo_a))
        call mat_ge_mmm(mos1%ca, s_ao, mos2%ca, s_mo_a, transa='T')
        if (uhf1 .or. uhf2) then
            allocate(s_mo_b(mos1%n_mo_b, mos2%n_mo_b))
            call mat_ge_mmm(mos1%cb, s_ao, mos2%cb, s_mo_b, transa='T')
        end if
        if (print_level >= 3) then
            write(stdout, '(5x,a)') 'Diagonal of the alpha MO overlap matrix: '
            write(stdout, '(6x,15f8.4)') [( s_mo_a(i, i), i=1, min(size(s_mo_a, 1), size(s_mo_a, 2)) )]
            if (allocated(s_mo_b)) then
                write(stdout, '(5x,a)') 'Diagonal of the beta MO overlap matrix: '
                write(stdout, '(10x,15(f10.5))') [( s_mo_b(i, i), i=1, min(size(s_mo_b, 1), size(s_mo_b, 2)) )]
            end if
        end if
        time_mo = omp_get_wtime() - time0
    end subroutine overlap_mo

end module overlap_mo_mod


!   DTGEVC Example Program Text
!   Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!   For licence see
!     https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

    Module dtgevc_example_mod

!     DTGEVC Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public :: normalize
!     .. Parameters ..
      Real (Kind=dp), Parameter, Public :: one = 1.0_dp
      Real (Kind=dp), Parameter, Public :: zero = 0.0_dp
      Integer, Parameter, Public :: nin = 5, nout = 6
    Contains
      Subroutine normalize(n, alphai, v, ldv)

!       .. Use Statements ..
        Use blas_interfaces, Only: dnrm2
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Integer, Intent (In) :: ldv, n
!       .. Array Arguments ..
        Real (Kind=dp), Intent (In) :: alphai(n)
        Real (Kind=dp), Intent (Inout) :: v(ldv, *)
!       .. Local Scalars ..
        Real (Kind=dp) :: a, b, r, r1, r2, v1, v2
        Integer :: i, j, k
!       .. Intrinsic Procedures ..
        Intrinsic :: sqrt
!       .. Executable Statements ..

        Do j = 1, n

          If (alphai(j)>=0.0_dp) Then
            If (alphai(j)==0.0_dp) Then
!             Real eigenvalue
!             The 2-norm of Q is calculated using dnrm2.
              r = dnrm2(n, v(1,j), 1)
              v(1:n, j) = v(1:n, j)/r
            Else
!             Complex eigenvalue (positive imaginary part)
!             Make largest element real and positive
              r1 = dnrm2(n, v(1,j), 1)
              r2 = dnrm2(n, v(1,j+1), 1)
              r1 = sqrt(r1**2+r2**2)
              r2 = -1.0_dp
              Do i = 1, n
                r = v(i, j)**2 + v(i, j+1)**2
                If (r>r2) Then
                  r2 = r
                  k = i
                End If
              End Do
              r = r1*sqrt(r2)
              a = v(k, j)/r
              b = v(k, j+1)/r
              Do i = 1, n
                v1 = v(i, j)
                v2 = v(i, j+1)
                v(i, j) = v1*a + v2*b
                v(i, j+1) = v2*a - v1*b
              End Do
            End If
          End If
        End Do
      End Subroutine
    End Module
    Program dtgevc_example

!     DTGEVC Example Main Program

!     .. Use Statements ..
      Use dtgevc_example_mod, Only: nin, normalize, nout, one, zero
      Use lapack_example_aux, Only: nagf_file_print_matrix_real_gen_comp
      Use lapack_interfaces, Only: dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, &
        dlacpy, dlaset, dorgqr, dormqr, dtgevc
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Integer :: i, icols, ifail, ihi, ilo, info, irows, jwork, lda, ldb, &
        ldvl, ldvr, lwork, m, n
      Logical :: ileft, iright
      Character (1) :: compq, compz, howmny, job, side
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), alphai(:), alphar(:), b(:, :), &
        beta(:), lscale(:), rscale(:), tau(:), vl(:, :), vr(:, :), work(:)
      Logical, Allocatable :: select(:)
      Character (0) :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
      Intrinsic :: nint
!     .. Executable Statements ..
      Write (nout, *) 'DTGEVC Example Program Results'
      Flush (nout)

!     ileft  is TRUE if left  eigenvectors are required
!     iright is TRUE if right eigenvectors are required

      ileft = .True.
      iright = .True.

!     Skip heading in data file

      Read (nin, *)
      Read (nin, *) n
      lda = n
      ldb = n
      ldvl = n
      ldvr = n
      lwork = 6*n
      Allocate (a(lda,n), alphai(n), alphar(n), b(ldb,n), beta(n), lscale(n), &
        rscale(n), tau(n), vl(ldvl,ldvl), vr(ldvr,ldvr), work(lwork), &
        select(n))

!     READ matrix A from data file
      Read (nin, *)(a(i,1:n), i=1, n)

!     READ matrix B from data file
      Read (nin, *)(b(i,1:n), i=1, n)

!     Balance matrix pair (A,B)
      job = 'B'
      Call dggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, &
        info)

!     Matrix A after balancing
!     ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call nagf_file_print_matrix_real_gen_comp('General', ' ', n, n, a, lda, &
        'F8.4', 'Matrix A after balancing', 'I', rlabs, 'I', clabs, 80, 0, &
        ifail)
      Write (nout, *)
      Flush (nout)

!     Matrix B after balancing
      ifail = 0
      Call nagf_file_print_matrix_real_gen_comp('General', ' ', n, n, b, ldb, &
        'F8.4', 'Matrix B after balancing', 'I', rlabs, 'I', clabs, 80, 0, &
        ifail)
      Write (nout, *)
      Flush (nout)

!     Reduce B to triangular form using QR
      irows = ihi + 1 - ilo
      icols = n + 1 - ilo
      Call dgeqrf(irows, icols, b(ilo,ilo), ldb, tau, work, lwork, info)

!     Apply the orthogonal transformation to matrix A
      Call dormqr('L', 'T', irows, icols, irows, b(ilo,ilo), ldb, tau, &
        a(ilo,ilo), lda, work, lwork, info)

!     Initialize VL (if left eigenvectors are required)
      If (ileft) Then

        Call dlaset('General', n, n, zero, one, vl, ldvl)
        Call dlacpy('Lower', irows-1, irows-1, b(ilo+1,ilo), ldb, &
          vl(ilo+1,ilo), ldvl)

        Call dorgqr(irows, irows, irows, vl(ilo,ilo), ldvl, tau, work, lwork, &
          info)
      End If

!     Initialize VR (if right eigenvectors are required)
      If (iright) Then
        Call dlaset('General', n, n, zero, one, vr, ldvr)
      End If

!     Compute the generalized Hessenberg form of (A,B)
      compq = 'V'
      compz = 'V'
      Call dgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, vl, ldvl, vr, &
        ldvr, info)

!     Matrix A in generalized Hessenberg form
      ifail = 0
      Call nagf_file_print_matrix_real_gen_comp('General', ' ', n, n, a, lda, &
        'F8.4', 'Matrix A in Hessenberg form', 'I', rlabs, 'I', clabs, 80, 0, &
        ifail)
      Write (nout, *)
      Flush (nout)

!     Matrix B in generalized Hessenberg form
      ifail = 0
      Call nagf_file_print_matrix_real_gen_comp('General', ' ', n, n, b, ldb, &
        'F8.4', 'Matrix B in Hessenberg form', 'I', rlabs, 'I', clabs, 80, 0, &
        ifail)

!     Routine DHGEQZ
!     Workspace query: jwork = -1

      jwork = -1
      job = 'S'

      Call dhgeqz(job, compq, compz, n, ilo, ihi, a, lda, b, ldb, alphar, &
        alphai, beta, vl, ldvl, vr, ldvr, work, jwork, info)

      Write (nout, *)
      Write (nout, 100) nint(work(1))
      Write (nout, 110) lwork
      Write (nout, *)
      Write (nout, 120)
      Write (nout, 130)

!     Compute the generalized Schur form
!     if the workspace lwork is adequate
!     The Schur form also gives parameters
!     required to compute generalized eigenvalues

      If (nint(work(1))<=lwork) Then

        Call dhgeqz(job, compq, compz, n, ilo, ihi, a, lda, b, ldb, alphar, &
          alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)

!       Print the generalized eigenvalues

        Do i = 1, n
          If (beta(i)/=0.0E0_dp) Then
            Write (nout, 140) i, '(', alphar(i)/beta(i), ',', &
              alphai(i)/beta(i), ')'
          Else
            Write (nout, 130) i
          End If
        End Do
        Write (nout, *)
        Flush (nout)

!       Compute left and right generalized eigenvectors
!       of the balanced matrix

        howmny = 'B'
        If (ileft .And. iright) Then
          side = 'B'
        Else If (ileft) Then
          side = 'L'
        Else If (iright) Then
          side = 'R'
        End If

        Call dtgevc(side, howmny, select, n, a, lda, b, ldb, vl, ldvl, vr, &
          ldvr, n, m, work, info)

        If (iright) Then

!         Compute right eigenvectors of the original matrix

          job = 'B'
          side = 'R'

          Call dggbak(job, side, n, ilo, ihi, lscale, rscale, n, vr, ldvr, &
            info)

          Call normalize(n, alphai, vr, ldvr)
!         Print the right eigenvectors

          ifail = 0
          Call nagf_file_print_matrix_real_gen_comp('General', ' ', n, n, vr, &
            ldvr, 'F8.4', 'Right eigenvectors', 'I', rlabs, 'I', clabs, 80, 0, &
            ifail)

          Write (nout, *)
          Flush (nout)
        End If

!       Compute left eigenvectors of the original matrix

        If (ileft) Then
          job = 'B'
          side = 'L'

          Call dggbak(job, side, n, ilo, ihi, lscale, rscale, n, vl, ldvl, &
            info)

          Call normalize(n, alphai, vl, ldvl)
!         Print the left eigenvectors

          ifail = 0
          Call nagf_file_print_matrix_real_gen_comp('General', ' ', n, n, vl, &
            ldvl, 'F8.4', 'Left eigenvectors', 'I', rlabs, 'I', clabs, 80, 0, &
            ifail)

        End If
      Else
        Write (nout, 150)
      End If

100   Format (1X, 'Minimal required LWORK = ', I6)
110   Format (1X, 'Actual value of  LWORK = ', I6)
120   Format (1X, 'Generalized eigenvalues')
130   Format (1X, I4, 5X, 'Infinite eigenvalue')
140   Format (1X, I4, 5X, A, F7.3, A, F7.3, A)
150   Format (1X, 'Insufficient workspace for array WORK', /, ' in DHGEQZ/', &
        'DHGEQZ')
    End Program


text:
c f77 -e extract.f -L$PCMDI/lib -lezget -lcdms -ldrs -L/usr/local/netcdf-2.3.2/libsrc -lnetcdf -o extract
c f77 -e -g extract.f -L/zooks0/ktaylor/pcmdi/util/ezget -lezgetdebug -L$PCMDI/lib -lcdms -ldrs -L/usr/local/netcdf-2.3.2/libsrc -lnetcdf -o extract
c  IBM
c xlf extract.f  -L/u/williams/ezget -lezget -L/u/williams/devel/cdms/lib -lcdms -L/u/williams/drs/lib -ldrs -L/u/williams/netcdf/libsrc -lnetcdf -o extract
c   HP
c fort77  +U77 extract.f -L/afs/nersc.gov/u/ahs/u4160/ezget -lezget -L/afs/nersc.gov/software/pub/drs/lib -lcdms -ldrs -lm -L/afs/nersc.gov/projects/graphics/netcdf/2.3.2pl4/lib -lnetcdf -o extract
c   sun/SOLARIS
c f77 -e extract.f -L/pcmdi/cirrus/ktaylor/ezget/solaris -lezget -L/pcmdi/cirrus/drs_solaris -lcdms -ldrs -L/pcmdi/cirrus/netcdf-2.3.2/libsrc -lnetcdf -o extract
c   SGI
c f77 extract.f -L/pcmdi/ktaylor/pcmdi/util/ezget/sgi -lezget -L/usr/local/cdms/lib -lcdms -L/usr/local/drs/drslib -ldrs -L/usr/local/lib -lnetcdf -o extract

       program extract

c   This program
c     ** retrieves the first 12 months of global data for a variable
c        (named 'tas' in this example) that is stored in a file (named
c        '/scratch/staff/lisa/amip_obs/nasa-amip_t').
c     ** creates an array of "weights" with elements proportional to
c        grid cell area (except for grid cells with missing data
c        where the "weight" is set to 0.0).
c     ** obtains the longitude and latitude coordinates and the
c        length of each dimension retrieved.
c     ** prints out a portion of the retrieved data.
c
c   Concerning the structure of the stored data, we assume the
c   following:
c
c    1.  The field is a function of longitude, latitude, and time, but
c        the dimension order is unknown.
c    2.  The data are stored in a rectangular array.
c    3.  The size of the array retrieved is limited as follows:
c        a. The longitude dimension is no longer than 100 elements.
c        b. The latitude dimension is no longer than 50 elements.
c        c. The time dimension is no longer than 12 elements.
c    4.  The names of the dimensions (as stored in the file) are:
c        'latitude', 'longitude', and 'time' (but not necessarily in
c        that order).
c    5.  The longitude coordinates are evenly spaced.
c    6.  The latitude coordinates are evenly spaced.
c    7.  The units for longitude and latitude are degrees.
c    8.  Only one field in the file has the name, 'tas' (i.e., this
c        variable name is unique in the file).
c    9.  The number, 1.0e20, is stored in the place of the true data
c        anywhere that data are missing (i.e., this is the "missing
c        data" value or indicator).
c
c
c    This program obtains the following information:
c
c    lons = the actual length of the longitude dimension of the data
c             retrieved from storage.
c    lats = the actual length of the latitude dimension of the data
c             retrieved from storage.
c    mons = the actual number of months of data retrieved from storage.
c    alon(100) = a vector containing the longitude coordinates for the
c                 data.
c    alat(50) = a vector containing the latitude coordinates for the
c                 data.
c    adata(100,50,12) = the retrieved array of data, which according to
c       our specifications given below, will be put in the following
c       structure (regardless of how it was originally stored):
c       1. The dimension order will be: longitude, latitude, time
c          (i.e., in the array, "adata", the first index is associated
c          with longitude, the second with latitude, and the third with
c          time).
c       2. The longitudes will be ordered from west to east, starting
c          with the longitude nearest to 180 W.
c       3. The latitudes will be ordered north to south.
c       4. The months will be ordered consecutively.
c    wtsmask(100,50,12) = the created "missing data" mask which will
c       be ordered the same as "adata", with elements set proportional
c       to the grid cell area (except for grid cells with missing data.
c       where the elements will be set to 0.0)
c
c    Note that if lons < nlon and/or lats < nlat, and/or mons < nmon,
c    then the extra elements of the array, "adata" and "wtsmask" will be
c    assigned a value of 0.0 by EzGet.

      parameter (nlon=100, nlat=50, nmon=12, n4=0)

      real adata(nlon,nlat,nmon), wtsmask(nlon,nlat,nmon),
     &      alon(nlon), alat(nlat)
      integer lons, lats, mons, i4

c     -------------------------------------------
c     Initialize EzGet:

      call initget

c     -------------------------------------------
c     Tell EzGet to consider "missing" any data that have
c     values (within a small tolerance) of 1.0e20.

      call defmisc('input missing value', 'real', 1.0e20)

c     -------------------------------------------
c     Define variable 1 as 'tas' and indicate path/filename where
c     it is stored.  In subsequent calls to EzGet (e.g., defdim
c     getdata, and getcoord), this field will be referenced by
c     the index assigned here (i.e., 1).

      call defvar(1, 'tas', '/scratch/staff/lisa/amip_obs/nasa-amip_t')

c     -------------------------------------------
c     Define domain for variable 1 and define order that EzGet will
c     retrieve data.

      call defdim(1, 1, 'longitude', 'width', 'nearest',
     &                                             -180.0, 180.0, 360.0)
      call defdim(1, 2, 'latitude', 'cosine', 'range',
     &                                               90.0, -90.0, 0.0)
      call defdim(1, 3, 'time',    'unit', 'nearest', 1.0, 12.0, 0.0)

c     where in the first call to defdim above the arguments indicate
c     the following:
c
c      1, 1, 'longitude':
c         The first 3 arguments indicate that for variable 1 (as
c         indicated by 1st argument), the data should be retrieved such
c         that the dimension named 'longitude' (as indicated by the
c         3rd argument) is the first dimension  (as indicated by the 2nd
c         argument).
c      'width': this argument controls the creation of the array of
c           "weights" that will be associated with the data. For each
c           array element the weight is equal to the product of the
c           weights defined by each dimension and here 'width' specifies
c           that the "weights" should be proportional to the
c           longitudinal width of each grid cell.
c      'nearest', -180.0, 180.0, 360.0:
c           these 4 arguments define the domain that will be retrieved
c           (and also the order of retrieval).  In this case all data
c           with a longitude coordinate roughly in the range -180 to 180
c           will be extracted, starting near -180.  The value 360.0
c           indicates that this coordinate is cyclical with a period of
c           360, so that EzGet will recognize equivalences such as
c           0. = 360. = -360. and -90. = 270.  Note that with 'nearest'
c           specified and with the range covering a complete cycle,
c           EzGet may shift the domain slightly so as to prevent a grid
c           cell near -180.0 (=180.0) from being split across the
c           boundary.  If, for example, the center of grid cells
c           are located at -180, -170, -160, ... 170, then EzGet will
c           shift the requested domain to -185 to 175 so the grid cell
c           at -180 will not be split.
c
c     and in the second call to defdim above the arguments indicate the
c     following:
c      'cosine': this option for controlling creation of weights
c           specifies that weights be generated equal approximately
c           to the cosine of latitude (i.e.,  abs(sin(bdry1) -
c           sin(bdry2)), where bdry1 and bdry2 are edges of the latitude
c           grid-cell, assumed to be half-way between grid-cell centers)
c           be assigned to all elements (except the weight will be 0.0
c           for grid cells with missing data).  Another option
c           ('gaussian') would be appropriate for spectral models with
c           gaussian grids (as opposed to the evenly spaced grid
c           accessed here).
c      'range', 90.0, -90.0, 0.0:
c           these 4 arguments define the domain that will be retrieved
c           (and also the order of retrieval).  In this case all data
c           with latitude coordinates in the range 90 to -90 will
c           be retrieved, starting near 90 (i.e. data will be retrieved
c           from north to south).  The value 0.0 indicates that this
c           coordinate is not cyclical.
c
c     and in the third call to defdim above the fourth argument
c     indicates the following:
c      'unit': this option for controlling creation of weights specifies
c           that unit weight should be given to each element of this
c           dimension (except the weight will be 0. for grid cells with
c           missing data).
c
c     -------------------------------------------
c     Extract variable 1 from file and create missing data mask.
c     The lengths of the longitude, latitude and time dimensions of
c     the data stored in the file are unknown, so initialize the
c     "expected" dimension lengths to 0.  Return the actual dimensions
c     of the retrieved data.

      lons = 0
      lats = 0
      mons = 12
      i4 = 0
      call getdata(1,  nlon,nlat,nmon,n4,  lons,lats,mons,i4,
     &             wtsmask, adata)

c     where on calling getdata the arguments have been defined as
c     follows:
c
c      1: The first argument has the value 1 and indicates that data
c         will be extracted for variable 1 (defined in defvar above).
c      nlon, nlat, nmon, n4:
c         These are the declared dimensions of the arrays, wtsmask and
c         adata.  Note that n4=0 because these arrays have only 3
c         dimensions.
c      lons, lats, mons, i4:
c         In this example, the user does not know how large the actual
c         longitude and latitude dimensions of the data being retrieved
c         will be, so these are set to 0.  The time dimension is
c         expected to be 12 (months) as specified in the earlier call to
c         defdim.  In general, if the user knows the size of the domain
c         to be retrieved, it is usually prudent to set these arguments
c         to the expected size of the domain, because then EzGet can
c         error exit if the actual size differs from what is expected.
c
c     On return from getdata:
c
c       1, nlon, nlat, nmon, n4:
c          The first 5 arguments will be unchanged.
c       lons, lats, mons, i4:
c          returns the actual dimensions of the data retrieved.
c       wtsmask: returns the weights associated with the data.
c          In this example, the weights will either be 0.0 or
c          proportional to grid cell area, depending on whether or not
c          the data are missing.
c       adata: returns the extracted data.

c     -------------------------------------------
c     Retrieve longitude and latitude coordinates of variable 1 (in
c     the same order as the retrieved data).

      call getcoord(1, 1, alon)
      call getcoord(1, 2, alat)

c     where in the first call to getcoord above, the arguments indicate
c     the following:
c
c      1: The 1st argument has the value 1 and indicates that we want
c           to obtain the coordinates of data extracted for variable 1
c           (as defined in the earlier call to defvar).
c      1: The 2nd argument has the value 1 and indicates that we want
c           to obtain the coordinates for the 1st dimension (as defined
c           in the earlier call to defdim).
c      alon: returns the coordinate values for the 1st dimension of
c           of variable 1 (i.e., 'longitude', as defined by defdim).

c     -------------------------------------------
c     Close all files opened by EzGet.

      call closeget

c     -------------------------------------------
c     Write data retrieved for last month.  Note that any missing data
c     will have been assigned the value 1.e20, so when written with
c     f8.3 format will appear as '********':

      write(*,'("longitudes = " / (10f8.3))') (alon(i), i=1,lons)

      write(*,'(/ "month = ", i3)') mons

      do 100 j=1,lats
        write(*,'(/"latitude = ", f8.3)') alat(j)
        write(*,'(10f8.3)') (adata(i,j,mons), i=1,lons)
  100 continue

      end

max_stars_repo_name: UV-CDAT/EzGet

max_issues_repo_name: UV-CDAT/EzGet

max_forks_repo_name: UV-CDAT/EzGet


! { dg-do compile }
! { dg-options "-O2 -fdump-tree-original" }
  integer :: x(2,2)
  if (any(x(:,:) .ne. reshape ((/ 3, 1, 4, 1 /), (/ 2, 2 /)))) STOP 1
end
! { dg-final { scan-tree-dump-times "atmp" 0 "original" } }


text:
recursive subroutine factorial(i,f)
  implicit none
  integer f,i,j
  intent(in) i
  intent(out) f
  if (i.le.0) then
    stop 'argument must be > 0'
  else if (i.eq.1) then
    f=i
  else
    call factorial(i-1,j)
    f=i*j
  end if
end subroutine factorial

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!    Date:       4 Nov 1992
!    Time:       14:12
!    Program:    MAXGRA.FOR
!    Version:    1.0
!    Programmer: Hans Los
!    Previous version(s):
!    0.0 -- 27 Sep 1989 --  8:42 -- Operating System: DOS
!
!  *********************************************************************
!  *      SUBROUTINE TO CALCULATE MAXIMUM GRAZING RATES AT WHICH       *
!  *                   ECOGROUPS CAN SURVIVE                           *
!  *********************************************************************
!
! 0895 MvdV GRAMX added for output of maximum grazing rate and
!           dimension ZOOD and ZOOPR adapted for more than one grazer
      SUBROUTINE MAXGRA(ROOT,EXTB,DAY,DSURF,EADJ,PMAXJ,CDATE,ZOOD,NUMGR,
     1                  NUMSP,DEP,GRAMX)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'blmdim.inc'
      INCLUDE 'size.inc'
      INCLUDE 'phyt1.inc'
      INCLUDE 'phyt2.inc'
      INCLUDE 'graas.inc'
      INCLUDE 'ioblck.inc'
      DIMENSION ROOT(2),OUTP(2,20),ZOOD(0:MG)
      CHARACTER*8 CDATE
!
!  Stop criteria:
!  1. Extinction of ecogroup j = background ext. +/- 0.05
!     (equivalent to a chlorophyll concentration < 10.0 mg/m3).
!  2. GHI-GLOW < 0.001.
!  3. After 20 iteration steps.
! If the algorithm stops by either criterium 2 or 3,
! Gmax will be preceeded by a "-" sign in the output.
!
!
! Print heading and compute ZMAX.
!
!
      IF (NUMSP .GT. 1) GO TO 50
      IF (NREP .GT. 1) GO TO 30
      CALL FORMFE (IOU(16))
      WRITE (IOU(16),10)
   10 FORMAT('  Maximum possible grazing rate computation.',/,
     1       '  Determine: 1. Gmax = the maximum grazing rate constan',
     2       't at which',/,'                a phytoplankton',
     3       ' type has a non-zero net growth rate.',/,
     4       '             2. Zmax = the corresponding zooplankton',
     5       ' concentration.',/,'  When more zooplankton is present',
     6       ' than Zmax, the growth rate of type j',/,'  MIGHT',
     7       ' (but does not necessarily) become negative.')
      LINEPA = 2
      WRITE (IOU(16),20) (SPNAME(K),K=1,NUSPEC)
   20 FORMAT (///,' Date',2X,'ZOOD',10X,20(A6,1X),/)
      SARG=ZOOK*XMIN
      ZMAX=(ZOOK+XMIN+2*DSQRT(SARG))/ZOOGR
   30 CONTINUE
      IF (NUGRAZ.EQ.0) THEN
        NG = 0
      ELSE
        NG = 1
      ENDIF
      LINEPA = LINEPA + 1
      IF (LINEPA .LE. 27) GO TO 50
      LINEPA = 1
      CALL FORMFE (IOU(16))
      WRITE (IOU(16),20) (SPNAME(K),K=1,NUSPEC)
   50 CONTINUE
      IF (ZOOPR(NUMSP,NG) .GT. 0.0 ) GO TO 70
!
!   Exit for unedible types, putting a large number in OUTP.
!
      DO 60 I=1,2
   60 OUTP(I,NUMSP+1)=1.D7
      IF (NUMSP .EQ. NUSPEC) GO TO 400
      RETURN
   70 IF (ROOT(2) .GT. EXTB) GO TO 90
!
!  Net growth rate is negative, even if G = 0.0
!
      DO 80 I=1,2
   80 OUTP(I,NUMSP+1)=0.0
      IF (NUMSP .EQ. NUSPEC) GO TO 400
      RETURN
   90 CONTINUE
      ITE=0
!
!  Establish initial boundaries for G interval. Start with:
!  Glow=0.0, and GHI is obtained by demanding EMIN=0.3 and R(T) and
!  M(T) of ecogroup j are both 0.1*PGMAX(T).
!
      GLOW=0.0
      GHI=0.03*PMAXJ*DAY-0.5
  100 GHI=GHI+0.5
      EADJST=EADJ+GHI/PMAXJ

!
!  Update nov 4 1992: use total depth rather than partial depth.
!  See also MAXGRO and BLOOM
!
*     CALL CONSTR(DSURF,DMIX(NUMSP),EADJST,ROOT,NUMGR)
      CALL CONSTR(DSURF,DEP,EADJST,ROOT,NUMGR)
!
!  If groupp j still has a positive growth rate (Kmax > Kb), increase
!  GHI with 0.5.
!
      IF (ROOT(2) .LT. EXTB) GO TO 110
      GLOW=GHI
      GO TO 100
!
!  Half G interval.
!
  110 GMID=(GHI+GLOW)/2.
  120 ITE=ITE+1
      EADJST=EADJ+GMID/PMAXJ
!
!  Update nov 4 1992: use total depth rather than partial depth.
!  See also MAXGRO and BLOOM
!
*     CALL CONSTR(DSURF,DMIX(NUMSP),EADJST,ROOT,NUMGR)
      CALL CONSTR(DSURF,DEP,EADJST,ROOT,NUMGR)
!
!  Test whether ROOT(2) (UKmax) is close enough to EXTB:
!  stop criterium 1.
!
      IF (ROOT(2) .GT. EXTB+0.05) GO TO 160
      IF (ROOT(2) .LT. EXTB-0.05) GO TO 130
      ZMAXJ=ZMAX*GMID/ZOOPR(NUMSP,NG)
      OUTP(1,NUMSP+1)=GMID
      OUTP(2,NUMSP+1)=ZMAXJ
      IF (NUMSP .EQ. NUSPEC) GO TO 400
      RETURN
!
!  Establish the next (lower) grazing rate.
!
  130 GHI=GMID
      GMID=(GMID+GLOW)/2.
!
!  Test for stop criterium 2: GHI-GLOW.
!
  140 GSTOP=GHI-GLOW
      IF (GSTOP .LT. 1.D-3) GO TO 150
      IF (ITE .LT. 20) GO TO 120
  150 ZMAXJ=ZMAX*GMID/ZOOPR(NUMSP,NG)
      OUTP(1,NUMSP+1)=GMID
      OUTP(2,NUMSP+1)=ZMAXJ
      IF (NUMSP .EQ. NUSPEC) GO TO 400
      RETURN
!
!  Establish the next (higher) grazing rate.
!
  160 GLOW=GMID
      GMID=(GHI+GMID)/2.
      GO TO 140
!
!   Print output for each ecogoup for this time-step.
!
  400 CONTINUE
      NTOT=NUSPEC+1

      GRAMX = 0.0
      DO 340 J=2,NTOT
        IF (OUTP(1,J).GT.GRAMX) GRAMX = OUTP(1,J)
  340 CONTINUE

      WRITE(IOU(16),180) CDATE,(OUTP(1,K),K=2,NTOT)
      WRITE(IOU(16),190) ZOOD(0),(OUTP(2,K),K=2,NTOT)
  180 FORMAT(1X,A4,10X,'Gmax',2X,20(F6.2,1X))
  190 FORMAT(8X,F6.1,1X,'Zmax',2X,20(F6.1,1X))
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
Drexel Drive starts at J Street and runs toward the east crossing L Street and ending past Chestnut Lane where it meets Snyder Drive at a line of bollards. Mostly residential and a good street for bicyclists traveling east west as it is a part of the Davis Bike Loop.

From west to east:

  J Street
  Auburn Drive
  Holmes Junior High School on the southern side of the street
  Hemlock Lane
  L Street, which at one point wasnt a stop for traffic on L Street going north and south.
  Pacific Drive
  Cypress Lane
  Chestnut Lane
  The bollards to Snyder Drive and the bike path to Pole Line Road and Loyola Drive

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program t
  integer,dimension(3)::i=(/1,2,3/)
  where (i>1)
    i=i*2
  elsewhere
    i=i*3
  endwhere
  print *,i
end program t


text:
	SUBROUTINE GDMMRK  ( marker, kx, ky, ix1, iy1, ix2, iy2, ixinc,
     +			     istag, iyinc, iret )
C************************************************************************
C* GDMMRK								*
C*									*
C* This subroutine plots markers for GDMAP.				*
C*									*
C* GDMMRK  ( MARKER, KX, KY, IX1, IY1, IX2, IY2, IXINC, ISTAG, IYINC, 	*
C*	     IRET )							*
C*									*
C* Input parameters:							*
C*	MARKER		CHAR*		Input for MARKER		*
C*	KX		INTEGER		Number of points in x dir	*
C*	KY		INTEGER		Number of points in y dir	*
C*	IX1		INTEGER		Initial x point			*
C*	IY1		INTEGER		Initial y point			*
C*	IX2		INTEGER		Last x point			*
C*	IY2		INTEGER		Last y point			*
C*	IXINC		INTEGER		X increment			*
C*	ISTAG		INTEGER		Stagger increment		*
C*	IYINC		INTEGER		Y increment			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 4/85						*
C* M. desJardins/GSFC	 8/88	Cleaned up				*
C* G. Huffman/GSC	 1/89	No plotting for color .lt. 1		*
C* S. Jacobs/NMC	 9/94	Added staggered plotting		*
C************************************************************************
	CHARACTER*(*)	marker
C------------------------------------------------------------------------
	iret = 0
C
C*	Translate marker input.
C
	CALL IN_MARK  ( marker, imcolr, ier )
C
C*	If color number .lt. 1, don't plot.
C
	IF  ( imcolr .gt. 0 )  THEN
C
C*	    Set marker color.
C
	    CALL GSCOLR  ( imcolr, ier )
C
C*	    Plot markers.
C
	    ixstrt = ix1
	    DO  j = iy1, iy2, iyinc
		fj = FLOAT (j)
		DO  i = ixstrt, ix2, ixinc
		    fi = FLOAT (i)
		    CALL GMARK  ( 'G', 1, fi, fj, ier )
		END DO
		IF  ( ixstrt .eq. ix1 )  THEN
		    ixstrt = ixstrt + istag
		ELSE
		    ixstrt = ix1
		END IF
	    END DO
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
oh my its me on the Front Page wiki :)

Im interested in evolution, genetics, Town Flora angiosperms, biodiversity, natural history, mycology, invasive species, conservation, and marine biology. I enjoy literature, specifically short stories and science fiction, and intellectual arguments. Accordingly, I studied Evolution, Ecology, and Biodiversity, obtaining my degree in March 2008 alongside a minor in English.

 Wiki Favs

   Dog Park  awww, doggies!
   Good Ideas For Dates  cause you know i cant come up with these things myself
   WiFi Hot Spots  good to check before you leave home
   MU Bus Phone 1  just love this picture

Davis Favs

 UC Davis campus
 Worlds Largest Bicycle Parade
 Sunrise Restaurant
 Nobu Hiro

20050824 16:29:50 nbsp Hey Jasmin, thanks for creating the optometrists page. I cleaned it up a little and put each Dr. on his/her own page. Also, remember to Preview to aviod cluttering up the server with repeated saves. Users/ArlenAbraham

20051011 11:15:02 nbsp RE Allegre Tennis Courts: Okay... so whose tennis courts are they? Is there another complex next door or are they private or what? If you know what the deal is, can you post info on the page? As an outsider, it appears that the courts are part of the complex! Users/AlphaDog

20060314 22:23:26 nbsp The Discovery, EEB Club, and Evolution, Ecology, and Biodiversity pages are confusing are these three different things or three names for the same thing? I have my theory, but I still think they need clarification. IMHO. Users/SteveDavison


text:
program part1
  integer :: n, a, b
  integer, dimension(12) :: bits
  character(len = 12) :: s

  n = 0
  do i = 1, 12
     bits(i) = 0
  end do
  do
     read (5, 12, end=20) s
12   format(a12)
     n = n+1
     do i = 1, 12
        if (s(i:i) == '1') then
           bits(i) = bits(i) + 1
        end if
     end do
  end do
20 a = 0
  b = 0
  do i = 1, 12
     a = a * 2
     b = b * 2
     if (bits(i) > n/2) then
        a = a + 1
     else
        b = b + 1
     end if
  end do
  write(*, '(I0)') a * b
end program part1

max_stars_repo_name: abeaumont/competitive-programming

max_issues_repo_name: abeaumont/competitive-programming

max_forks_repo_name: abeaumont/competitive-programming


text:
program myfortran

    ! --- Put use statements here
    !
    ! use mymodule

    implicit none

    ! --- Declare your variables here
    !
    ! integer :: a

    print*, 'Hello, Fortran!'

contains

    ! --- Add your subroutines and functions here
    !
    ! subroutine mysub
    !
    ! end subroutine mysub

end program myfortran

max_stars_repo_name: jonaslindemann/qtcreator-fortran

max_issues_repo_name: jonaslindemann/qtcreator-fortran

max_forks_repo_name: jonaslindemann/qtcreator-fortran


text:
! RUN: not %flang -fsyntax-only -fopenmp %s 2>&1 | FileCheck %s
! OpenMP Version 4.5
! 2.5 parallel construct.
! A program that branches into or out of a parallel region
! is non-conforming.

program omp_parallel
  integer i, j, k

  !CHECK: invalid branch into an OpenMP structured block
  goto 10

  !$omp parallel
  do i = 1, 10
    do j = 1, 10
      print *, "Hello"
      !CHECK: In the enclosing PARALLEL directive branched into
      !CHECK: STOP statement is not allowed in a PARALLEL construct
      10 stop
    end do
  end do
  !$omp end parallel

end program omp_parallel

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20080412 08:28:58 nbsp  Welcome to the Wiki. Please consider the importance of using your RealName. Also Please read Welcome to the Wiki/Business Owner; it will help explain how you can make the wiki a positive experience for you without clashing with established wiki social norms. Users/JasonAller

20081111 20:07:12 nbsp Your recent addition of your business listing to all those other wikis is pretty far into a gray area. Particularly since the listings didnt even modify the address to mention that the business is in Davis. Users/JasonAller

20081111 20:32:34 nbsp I really dont think anybody will come to Davis from Fresno just for a smog check, so I deleted the page. There are plenty of alternative businesses in Fresno and it would awesome if you could make pages for them. Thanks! Users/NickSchmalenberger


text:
      SUBROUTINE SPCTRM(P,M,K,OVRLAP,W1,W2)
      LOGICAL OVRLAP
      DIMENSION P(M),W1(*),W2(M)
      WINDOW(J)=(1.-ABS(((J-1)-FACM)*FACP))
C     WINDOW(J)=1.
C     WINDOW(J)=(1.-(((J-1)-FACM)*FACP)**2)
      MM=M+M
      M4=MM+MM
      M44=M4+4
      M43=M4+3
      DEN=0.
      FACM=M-0.5
      FACP=1./(M+0.5)
      SUMW=0.
      DO 11 J=1,MM
        SUMW=SUMW+WINDOW(J)**2
11    CONTINUE
      DO 12 J=1,M
        P(J)=0.
12    CONTINUE
      IF(OVRLAP)THEN
        READ (9,*) (W2(J),J=1,M)
      ENDIF
      DO 18 KK=1,K
        DO 15 JOFF=-1,0,1
          IF (OVRLAP) THEN
            DO 13 J=1,M
              W1(JOFF+J+J)=W2(J)
13          CONTINUE
            READ (9,*) (W2(J),J=1,M)
            JOFFN=JOFF+MM
            DO 14 J=1,M
              W1(JOFFN+J+J)=W2(J)
14          CONTINUE
          ELSE
            READ (9,*) (W1(J),J=JOFF+2,M4,2)
          ENDIF
15      CONTINUE
        DO 16 J=1,MM
          J2=J+J
          W=WINDOW(J)
          W1(J2)=W1(J2)*W
          W1(J2-1)=W1(J2-1)*W
16      CONTINUE
        CALL FOUR1(W1,MM,1)
        P(1)=P(1)+W1(1)**2+W1(2)**2
        DO 17 J=2,M
          J2=J+J
          P(J)=P(J)+W1(J2)**2+W1(J2-1)**2
     *        +W1(M44-J2)**2+W1(M43-J2)**2
17      CONTINUE
        DEN=DEN+SUMW
18    CONTINUE
      DEN=M4*DEN
      DO 19 J=1,M
        P(J)=P(J)/DEN
19    CONTINUE
      RETURN
      END

ext: for


text:
program da_vp_bilin

!----------------------------------------------------------------------
! Purpose: Regridding from low to high resolution in control variable space
!                       by using bilinear interpolation
!
! where n  is the grid number in x or y
!       ns is the refinement ratio between two resulotions
!
! Method: follow da_bilin.f90
!
! Compile:
!
!   pgf90 -o da_vp_bilin.exe da_vp_bilin.f90
!
! liuz@ucar.edu , 2016-08, NCAR/MMM
!----------------------------------------------------------------------

  !use netcdf

  implicit none

  !These variables' incremental will be regridded by default
  character (len=6), dimension(1:19) :: vNam

  integer :: ix, jy, kz, k, status
  integer :: ixh, jyh, kzh
  integer :: nLat, nLon, oLat, oLon
  integer :: sLat, eLat, sLon, eLon
  integer :: rLat, rLon

  real, dimension(:,:,:), allocatable :: v1, v2, v3, v4, v5
  real, dimension(:,:,:), allocatable :: v6, v7, v8, v9, v10, v11
  real, dimension(:,:,:), allocatable :: v1h, v2h, v3h, v4h, v5h
  real, dimension(:,:,:), allocatable :: v6h, v7h, v8h, v9h, v10h, v11h
  real, dimension(:,:),   allocatable :: iVar, oVar

  character (len = 255)  :: appname   = ""
  character (len = 255)  :: arg       = ""
  character (len = 19)  :: analysis_date
  character (len = 255)  :: input_file= "vp_output.global"
  character (len = 255)  :: output_file= "vp_output.global_hires"

  integer, parameter :: vp_unit = 8
  integer, parameter :: vp_hires_unit = 9
  integer :: ratio ! resolution ratio
  integer :: cloud_cv_options ! 2 or 3 with cloud cv variables
  integer :: use_cv_w         ! =1 for w control variable
  integer :: io_status
  integer iargc

  LOGICAL :: file_exists

  !These variables' incremental will be regridded by default

  !call getarg(0, appname)
  !n=index(appname, '/', BACK=.true.)
  !appname = trim(appname(n+1:))

  call getarg(1, arg)
  call getarg(2, arg)
  read(arg, '(i3)') ratio

  call getarg(3, arg)
  call getarg(4, arg)
  read(arg, '(i3)') cloud_cv_options

  call getarg(5, arg)
  call getarg(6, arg)
  read(arg, '(i3)') use_cv_w

  write (*, *) 'ratio = ', ratio, 'cloud_cv_options = ', cloud_cv_options, &
               'use_cv_w = ', use_cv_w

!  read vp file
!--------------------
  inquire(FILE=trim(input_file), EXIST=file_exists)

  if ( .not. file_exists ) then
    Write(*,*) "\nError: "//trim(input_file)//" not exists\n"
    call exit(-1)
  else
    Write(*,*) "Found: "//trim(input_file)
  endif

  open(unit=vp_unit,file=trim(input_file),iostat=io_status,form='UNFORMATTED',status='OLD')
  if (io_status /= 0) then
     write(*,*) "Error ",io_status," opening vp file "//trim(input_file)
     call exit(-1)
  end if
  write(*,*) 'Reading vp from : '//trim(input_file)
  !read(vp_unit) analysis_date
  !print *, 'analysis_date = ', analysis_date
  read(vp_unit) ix, jy, kz ! domain dimension (unstagered)
  print *, "input file: ix, jy, kz = ", ix, jy, kz

  allocate ( v1 (1:ix,1:jy,1:kz))
  allocate ( v2 (1:ix,1:jy,1:kz))
  allocate ( v3 (1:ix,1:jy,1:kz))
  allocate ( v4 (1:ix,1:jy,1:kz))
  allocate ( v5 (1:ix,1:jy,1:kz))

  read(vp_unit) v1, v2, v3, v4, v5

  if ( cloud_cv_options >= 2 ) then
    allocate ( v6 (1:ix,1:jy,1:kz))
    allocate ( v7 (1:ix,1:jy,1:kz))
    allocate ( v8 (1:ix,1:jy,1:kz))
    allocate ( v9 (1:ix,1:jy,1:kz))
    allocate ( v10 (1:ix,1:jy,1:kz))
    read(vp_unit) v6, v7, v8, v9, v10
  end if

  if ( use_cv_w == 1 ) then
    allocate ( v11 (1:ix,1:jy,1:kz))
    read(vp_unit) v11
  end if

  write(*,*) 'End Reading vp from : '//trim(input_file)
  close(vp_unit)
!-----------------------------
! end read vp file
!----------------------

    nLon = ix + 2 ! 52
    nLat = jy + 2 ! 52

    rLon = ix * ratio ! 150
    rLat = jy * ratio ! 150

    oLon = ( nLon - 1 ) * ratio + 1 ! 154
    oLat = ( nLat - 1 ) * ratio + 1

    elat = (oLat - rLat) / 2       ! 2
    slat =  oLat - rLat - elat + 1 ! 3

    elon = (oLon - rLon) / 2
    slon =  oLon - rLon - elon + 1

    allocate(iVar(nLon, nLat), stat=status)
    allocate(oVar(oLon, oLat), stat=status)

  ixh = ix*ratio
  jyh = jy*ratio

  allocate ( v1h (1:ixh,1:jyh,1:kz))
  allocate ( v2h (1:ixh,1:jyh,1:kz))
  allocate ( v3h (1:ixh,1:jyh,1:kz))
  allocate ( v4h (1:ixh,1:jyh,1:kz))
  allocate ( v5h (1:ixh,1:jyh,1:kz))

  if ( cloud_cv_options >= 2 ) then
    allocate ( v6h (1:ixh,1:jyh,1:kz))
    allocate ( v7h (1:ixh,1:jyh,1:kz))
    allocate ( v8h (1:ixh,1:jyh,1:kz))
    allocate ( v9h (1:ixh,1:jyh,1:kz))
    allocate ( v10h (1:ixh,1:jyh,1:kz))
  end if

  if ( use_cv_w == 1 ) then
    allocate ( v11h (1:ixh,1:jyh,1:kz))
  end if

  do k = 1, kz
     iVar(2:nlon-1,2:nlat-1) = v1(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v1h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

     iVar(2:nlon-1,2:nlat-1) = v2(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v2h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

     iVar(2:nlon-1,2:nlat-1) = v3(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v3h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

     iVar(2:nlon-1,2:nlat-1) = v4(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v4h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

     iVar(2:nlon-1,2:nlat-1) = v5(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v5h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

   if ( cloud_cv_options >= 2 ) then
     iVar(2:nlon-1,2:nlat-1) = v6(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v6h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

     iVar(2:nlon-1,2:nlat-1) = v7(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v7h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

     iVar(2:nlon-1,2:nlat-1) = v8(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v8h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

     iVar(2:nlon-1,2:nlat-1) = v9(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v9h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)

     iVar(2:nlon-1,2:nlat-1) = v10(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v10h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
   end if

   if ( use_cv_w == 1 ) then
     iVar(2:nlon-1,2:nlat-1) = v11(:,:,k)
     iVar(1,:)               = iVar(2,:)
     iVar(nlon,:)            = iVar(nlon-1,:)
     iVar(:,1)               = iVar(:,2)
     iVar(:,nlat)            = iVar(:,nlat-1)
     call bilin(iVar,nLon,nLat,ratio,oVar,oLon,oLat)
     v11h(:,:,k) = oVar(slon:olon-elon,slat:olat-elat)
   end if
  enddo

     open(unit=vp_hires_unit,file=trim(output_file),iostat=io_status,form='UNFORMATTED',status='UNKNOWN')
     if (io_status /= 0) then
        write(*,*) "Error ",io_status," opening vp file "//trim(output_file)
        call exit(-1)
     end if
     write(*,*) 'Writting vp on hires. to : '//trim(output_file)

    print *, 'output file: ixh, jyh, kz=', ixh, jyh, kz
    write(vp_hires_unit) ixh, jyh, kz
    write(vp_hires_unit) v1h,v2h,v3h,v4h,v5h
    if ( cloud_cv_options >= 2 ) then
      write(vp_hires_unit) v6h,v7h,v8h,v9h,v10h
    end if
    if ( use_cv_w == 1 ) then
      write(vp_hires_unit) v11h
    end if

    deallocate(v1, stat=status)
    deallocate(v2, stat=status)
    deallocate(v3, stat=status)
    deallocate(v4, stat=status)
    deallocate(v5, stat=status)

    deallocate(v1h, stat=status)
    deallocate(v2h, stat=status)
    deallocate(v3h, stat=status)
    deallocate(v4h, stat=status)
    deallocate(v5h, stat=status)

  if ( cloud_cv_options >= 2 ) then
    deallocate(v6, stat=status)
    deallocate(v7, stat=status)
    deallocate(v8, stat=status)
    deallocate(v9, stat=status)
    deallocate(v10, stat=status)

    deallocate(v6h, stat=status)
    deallocate(v7h, stat=status)
    deallocate(v8h, stat=status)
    deallocate(v9h, stat=status)
    deallocate(v10h, stat=status)
  end if

  if ( use_cv_w == 1 ) then
    deallocate(v11, stat=status)
    deallocate(v11h, stat=status)
  end if

  Write(*,*) "Regridding increment completed successfully"

contains
  subroutine show_usage()
     Write(*,*) 'Usage :'//trim(appname)// &
     '[-h] [-fg_lores filename] [-an_lores filename] [-fg_hires filename] [-ns n ] [-o outputfile]'
     Write(*,*) "  -fg_lores    Optional, low resulotion first guess file,                 default - fg"
     Write(*,*) "  -an_lores    Optional, low resulotion analysis file comes from wrfvar,  default - wrfvar_output"
     Write(*,*) "  -fg_hires    Optional, high resultion first guess file,                 default - wrfinput_hires"
     Write(*,*) "  -ns          Optional, the refinement ratio between two resulotions,    default - 3"
     Write(*,*) "  -o           Optional, output high resulotion analysis file,            default - wrfvar_output_hires"
     Write(*,*) "  -h           Show this help"
  end subroutine show_usage

  !subroutine nf90_handle_err(status, errmsg)
  !  integer,          intent(in) :: status
  !  character(len=*), intent(in) :: errmsg
!
!    if(status /= nf90_noerr) then
!      print *, trim(nf90_strerror(status))//" : "//trim(errmsg)
!      Stop
!    end if
!  end subroutine nf90_handle_err

  subroutine bilin(old,xi,yi,ns,new,xo,yo)

! assume: xo = (xi-1)*ns + 1, xi=50, xo=49*3+1=148
!         yo = (yi-1)*ns + 1

    implicit none

    integer,                 intent(in) :: xi,yi,xo,yo
    real, dimension(xi,yi),  intent(in) :: old
    integer,                 intent(in) :: ns
    real, dimension(xo,yo),  intent(out):: new

    real  :: im(1:ns+1,2)
!    real  :: imm(1:ns+3,2)
    integer:: i,j,jm1,im1,ix1,ix2,iy1,iy2

    forall(i=1:ns+1) im(i,2) = real(i-1)/ns
      im(:,1) = 1 - im(:,2)

    do j=2,yi
      jm1 = j - 1
      iy2 = jm1 * ns + 1
      iy1 = iy2 - ns
      do i=2,xi
        im1 = i - 1
        ix2 = im1 * ns + 1
        ix1 = ix2 - ns
        new(ix1:ix2,iy1:iy2) = matmul(im,matmul(old(im1:i,jm1:j),transpose(im)))
      end do
    end do

   !   ns = ns + 2
   ! forall(i=1:ns+1) imm(i,2) = real(i-1)/ns
   !   imm(:,1) = 1 - imm(:,2)
!
!      j=yi
!      jm1 = j - 1
!      iy2 = jm1 * ns + 1
!      iy1 = iy2 - ns
!
!      i=xi
!        im1 = i - 1
!        ix2 = im1 * ns + 1
!        ix1 = ix2 - ns
!        new(ix1:ix2,iy1:iy2) = matmul(imm,matmul(old(im1:i,jm1:j),transpose(imm)))
!      end do
!    end do

  end subroutine bilin

end program da_vp_bilin

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: matzegoebel/WRF-fluxavg

max_forks_repo_name: matzegoebel/WRF-fluxavg


text:
module LatLonOutputModule
!******************************************************************************
!	Peter A. Bosler
!	Department of Mathematics
!	University of Michigan
!	pbosler@umich.edu
!
!******************************************************************************
!
!	Defines output data structures and methods for plotting Lagrangian meshes of the sphere.
!
!	Bosler, P.A., "Particle Methods for Geophysical Flow on the Sphere," PhD Thesis; the University of Michigan, 2013.
!
!----------------
use NumberKindsModule
use SphereGeomModule
use IntegerListModule
use LoggerModule
use ParticlesModule
use EdgesModule
use PanelsModule
use SphereMeshModule
use STRIPACKInterfaceModule
use SSRFPACKInterfaceModule

implicit none
private
public LLSource
public New, Delete
public LLOutputMatlab, UpdateFilename

type LLSource
	character(len = 256) :: filename
	character(len = 128) :: title
	type(STRIPACKData) :: delTri
	type(SSRFPACKData) :: vectorSource
	type(SSRFPACKData), pointer :: scalarSource => null()
	integer(kint) :: nTracer
	integer(kint) :: problemKind
	type(SSRFPACKData), pointer :: tracerSource => null()
	integer(kint) :: nLon
	real(kreal), pointer :: lats(:) => null(),&
				   			lons(:) => null()
end type
!
!----------------
! Logging
!----------------
!
logical(klog), save :: logInit = .FALSE.
type(Logger) :: log
character(len=28), save :: logKey = 'LatLonOutput'
integer(kint), parameter :: logLevel = DEBUG_LOGGING_LEVEL
character(len=28) :: formatString
character(len=128) :: logString

!
!----------------
! Interfaces
!----------------
!
interface New
	module procedure NewPrivate
end interface

interface Delete
	module procedure DeletePrivate
end interface

interface UpdateFilename
	module procedure UpdateFilenameLL
end interface

contains
!
!----------------
! Standard methods : Constructor / Destructor, Copy
!----------------
!

subroutine NewPrivate(self,aMesh,filename,nLon)
	type(LLSource), intent(out) :: self
	type(SphereMesh), intent(in) :: aMesh
	character(len=*), intent(in) :: filename
	integer(kint), intent(in), optional :: nLon
	! local variables
	real(kreal) :: dLambda
	integer(kint) :: j

	if (.not. logInit) call InitLogger(log,procRank)

	self%filename = filename
	self%nLon = nLon

	if ( present(nLon) ) then
		dLambda = 2.0_kreal*PI/real(nLon,kreal)
	else
		dLambda = PI/180.0_kreal
	endif
	allocate(self%lons(nLon))
	allocate(self%lats(nLon/2 + 1))
	do j=1,nLon
		self%lons(j) = (j-1)*dLambda
	enddo
	do j=1,nLon/2+1
		self%lats(j) = -PI/2.0_kreal + (j-1)*dLambda
	enddo
	call New(self%delTri,aMesh)
	call DelaunayTriangulation(self%delTri)

	call New(self%vectorSource,self%delTri,.TRUE.)

	self%problemKind = aMesh%problemKind
	if ( self%problemKind /= ADVECTION_SOLVER ) then
		allocate(self%scalarSource)
		call New(self%scalarSource,self%delTri,.FALSE.)
	endif

	self%nTracer = aMesh%nTracer
	if ( self%nTracer > 0 ) then
		allocate(self%tracerSource)
		call New(self%tracerSource, self%delTri, self%nTracer)
	endif

end subroutine

subroutine DeletePrivate(self)
	type(LLSource), intent(inout) :: self
	deallocate(self%lats)
	deallocate(self%lons)
	call Delete(self%deltri)
	call Delete(self%scalarSource)
	call Delete(self%vectorSource)
	if ( self%nTracer > 0 ) then
		call Delete(self%tracerSource)
		deallocate(self%tracerSource)
	endif
	if ( associated(self%scalarSource) ) then
		call Delete(self%scalarSource)
		deallocate(self%scalarSource)
	endif
end subroutine

!
!----------------
! Public functions
!----------------
!
subroutine UpdateFileNameLL(self,filename)
	type(LLSource), intent(inout) :: self
	character(len=*), intent(in) :: filename
	self%filename = trim(filename)
end subroutine

subroutine LLOutputMatlab(self,aMesh)
	type(LLSource), intent(inout) :: self
	type(SphereMesh), intent(in) :: aMesh
	!
	integer(kint) :: i, j, k, writeStat
	real(kreal) :: interpVector(3), interpScalar, xyzLoc(3), lonUnitVector(3), latUnitVector(3)
	real(kreal), allocatable :: llgrid1(:,:), llgrid2(:,:)

	! open output file
	open(unit = WRITE_UNIT_1, file = self%filename, status = 'REPLACE', action = 'WRITE', iostat = writeStat)
	if ( writeStat /= 0 ) then
		call LogMessage(log,ERROR_LOGGING_LEVEL,'LLOutputMatlab : ', 'ERROR opening output file.')
		return
	endif

	! write lon/lat vectors
	write(WRITE_UNIT_1,'(A)',ADVANCE = 'NO') 'lon = [ '
	do i=1,self%nlon-1
		write(WRITE_UNIT_1,'(F24.15,A)') self%lons(i), '; ...'
	enddo
	write(WRITE_UNIT_1,'(F24.15,A)') self%lons(self%nLon), ' ] ; '

	write(WRITE_UNIT_1,'(A)',advance = 'NO') 'lat = [ '
	do i=1,self%nlon/2
		write(WRITE_UNIT_1,'(F24.15,A)') self%lats(i), '; ...'
	enddo
	write(WRITE_UNIT_1,'(F24.15,A)') self%lats(self%nlon/2+1), ' ] ;'

	allocate(llgrid1(self%nLon/2+1,self%nLon))
	llgrid1 = 0.0_kreal
	allocate(llgrid2(self%nLon/2+1,self%nLon))
	llgrid2 = 0.0_kreal
	!
	! Lagrangian parameter
	!
	call SetSourceLagrangianParameter(self%vectorSource,self%delTri)
	do j=1,self%nLon
		do i=1,self%nLon/2+1
			xyzLoc = EARTH_RADIUS*[ cos(self%lats(i))*cos(self%lons(j)), &
						cos(self%lats(i))*sin(self%lons(j)), sin(self%lats(i)) ]
			interpVector = InterpolateVector(xyzLoc, self%vectorSource, self%delTri)
			llgrid1(i,j) = Longitude(interpVector)
			llgrid2(i,j) = Latitude(interpVector)
		enddo
	enddo

	write(WRITE_UNIT_1,'(A)',advance = 'NO') ' alphaLongitude = [ '
	do i=1,self%nLon/2
		do j=1,self%nLon -1
			write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(i,j)
		enddo
		write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(i,self%nLon), ' ; ...'
	enddo
	do j=1,self%nLon -1
		write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(self%nLon/2+1,j)
	enddo
	write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(self%nLon/2+1,self%nLon), ' ] ; '

	write(WRITE_UNIT_1,'(A)',advance = 'NO') ' alphaLatitude = [ '
	do i=1,self%nLon/2
		do j=1,self%nLon -1
			write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid2(i,j)
		enddo
		write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid2(i,self%nLon), ' ; ...'
	enddo
	do j=1,self%nLon -1
		write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid2(self%nLon/2+1,j)
	enddo
	write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid2(self%nLon/2+1,self%nLon), ' ] ; '

	!
	!	velocity
	!
	call SetSourceVelocity(self%vectorSource,self%delTri)
	llgrid1 = 0.0_kreal
	llgrid2 = 0.0_kreal
	do j=1,self%nLon
		do i=2,self%nLon/2
			xyzLoc = EARTH_RADIUS*[ cos(self%lats(i))*cos(self%lons(j)), &
						cos(self%lats(i))*sin(self%lons(j)), sin(self%lats(i)) ]
			interpVector = InterpolateVector(xyzLoc, self%vectorSource, self%delTri)

			lonUnitVector = [ - xyzLoc(2), &
								xyzLoc(1), 0.0_kreal ]
			lonUnitVector = lonUnitVector/(EARTH_RADIUS*sqrt(EARTH_RADIUS*EARTH_RADIUS - xyzLoc(3)*xyzLoc(3)))

			latUnitVector = [ -xyzLoc(1)*xyzLoc(3), &
							  -xyzLoc(2)*xyzLoc(3), &
							  EARTH_RADIUS*EARTH_RADIUS - xyzLoc(3)*xyzLoc(3) ]
			latUnitVector = latUnitVector/(EARTH_RADIUS*EARTH_RADIUS*sqrt(EARTH_RADIUS*EARTH_RADIUS - xyzLoc(3)*xyzLoc(3)))

			llgrid1(i,j) = sum(interpVector*lonUnitVector)

			llgrid2(i,j) = sum(interpVector*latUnitVector)

		enddo
	enddo

	write(WRITE_UNIT_1,'(A)',advance = 'NO') ' zonalVelocityU = [ '
	do i=1,self%nLon/2
		do j=1,self%nLon -1
			write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(i,j)
		enddo
		write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(i,self%nLon), ' ; ...'
	enddo
	do j=1,self%nLon -1
		write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(self%nLon/2+1,j)
	enddo
	write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(self%nLon/2+1,self%nLon), ' ] ; '

	write(WRITE_UNIT_1,'(A)',advance = 'NO') ' meridionalVelocityV = [ '
	do i=1,self%nLon/2
		do j=1,self%nLon -1
			write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid2(i,j)
		enddo
		write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid2(i,self%nLon), ' ; ...'
	enddo
	do j=1,self%nLon -1
		write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid2(self%nLon/2+1,j)
	enddo
	write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid2(self%nLon/2+1,self%nLon), ' ] ; '

	!
	!	tracers
	!
	if ( self%nTracer > 0 ) then
		call SetSourceTracer(self%tracerSource, self%delTri)
		llgrid1 = 0.0_kreal
		do k = 1, self%nTracer
			do j=1,self%nLon
				do i=1,self%nLon/2 + 1
					xyzLoc = EARTH_RADIUS*[ cos(self%lats(i))*cos(self%lons(j)), &
						cos(self%lats(i))*sin(self%lons(j)), sin(self%lats(i)) ]

					llgrid1(i,j) = InterpolateTracer(xyzLoc, self%tracerSource, self%delTri, k)
				enddo
			enddo

			write(WRITE_UNIT_1,'(A,I1,A)',advance = 'NO') ' tracer', k,' = [ '
			do i=1,self%nLon/2
				do j=1,self%nLon -1
					write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(i,j)
				enddo
				write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(i,self%nLon), ' ; ...'
			enddo
			do j=1,self%nLon -1
				write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(self%nLon/2+1,j)
			enddo
			write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(self%nLon/2+1,self%nLon), ' ] ; '
		enddo
	endif

	!
	!	bve variables
	!
	if ( self%problemKind == BVE_SOLVER) then
		!
		!	absolute vorticity
		!
		llgrid1 = 0.0_kreal
		call SetSourceAbsVort(self%scalarSource,self%delTri)
		do j=1,self%nLon
			do i=1,self%nLon/2 + 1
				xyzLoc = EARTH_RADIUS*[ cos(self%lats(i))*cos(self%lons(j)), &
					cos(self%lats(i))*sin(self%lons(j)), sin(self%lats(i)) ]

				llgrid1(i,j) = InterpolateScalar(xyzLoc, self%scalarSource, self%delTri)
			enddo
		enddo

		write(WRITE_UNIT_1,'(A)',advance = 'NO') ' absVort = [ '
		do i=1,self%nLon/2
			do j=1,self%nLon -1
				write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(i,j)
			enddo
			write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(i,self%nLon), ' ; ...'
		enddo
		do j=1,self%nLon -1
			write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(self%nLon/2+1,j)
		enddo
		write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(self%nLon/2+1,self%nLon), ' ] ; '

		!
		!	relative vorticity
		!
		llgrid1 = 0.0_kreal
		call SetSourceRelVort(self%scalarSource,self%delTri)
		do j=1,self%nLon
			do i=1,self%nLon/2 + 1
				xyzLoc = EARTH_RADIUS*[ cos(self%lats(i))*cos(self%lons(j)), &
					cos(self%lats(i))*sin(self%lons(j)), sin(self%lats(i)) ]

				llgrid1(i,j) = InterpolateScalar(xyzLoc, self%scalarSource, self%delTri)
			enddo
		enddo
		write(WRITE_UNIT_1,'(A)',advance = 'NO') ' relVort = [ '
		do i=1,self%nLon/2
			do j=1,self%nLon -1
				write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(i,j)
			enddo
			write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(i,self%nLon), ' ; ...'
		enddo
		do j=1,self%nLon -1
			write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(self%nLon/2+1,j)
		enddo
		write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(self%nLon/2+1,self%nLon), ' ] ; '

		!
		!	kinetic energy
		!
		llgrid1 = 0.0_kreal
		call SetSourceKineticEnergy(self%scalarSource,self%delTri)
		do j=1,self%nLon
			do i=1,self%nLon/2 + 1
				xyzLoc = EARTH_RADIUS*[ cos(self%lats(i))*cos(self%lons(j)), &
					cos(self%lats(i))*sin(self%lons(j)), sin(self%lats(i)) ]

				llgrid1(i,j) = InterpolateScalar(xyzLoc, self%scalarSource, self%delTri)
			enddo
		enddo
		write(WRITE_UNIT_1,'(A)',advance = 'NO') ' ke = [ '
		do i=1,self%nLon/2
			do j=1,self%nLon -1
				write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(i,j)
			enddo
			write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(i,self%nLon), ' ; ...'
		enddo
		do j=1,self%nLon -1
			write(WRITE_UNIT_1,'(F24.10)',advance='NO') llgrid1(self%nLon/2+1,j)
		enddo
		write(WRITE_UNIT_1,'(F24.10, A)',advance='YES') llgrid1(self%nLon/2+1,self%nLon), ' ] ; '
	endif

	!
	!	TO DO : SWE VARIABLES
	!

	deallocate(llgrid1)
	deallocate(llgrid2)
	close(WRITE_UNIT_1)
end subroutine

!
!----------------
! Module methods : type-specific functions
!----------------
!
subroutine InitLogger(aLog,rank)
	type(Logger), intent(inout) :: aLog
	integer(kint), intent(in) :: rank
	write(logKey,'(A,A,I0.2,A)') trim(logKey),'_',rank,' : '
	call New(aLog,logLevel)
	logInit = .TRUE.
end subroutine

end module

max_stars_repo_name: pbosler/LPPM

max_issues_repo_name: pbosler/LPPM

max_forks_repo_name: pbosler/LPPM