text:
!*******************************************************************************
!>  @file bivariate.f
!>  @brief Contains module @ref bivariate.
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!>  This modules contains routines for interpolating points inside a grid. This
!>  was originally written by S. P. Hirshman. It has been modified for thread
!>  safety converting it to be object oriented by M. Cianciosa.
!*******************************************************************************
      MODULE bivariate

      USE stel_kinds
      USE profiler

      IMPLICIT NONE

!*******************************************************************************
!  DERIVED-TYPE DECLARATIONS
!  1) bivariate_type
!
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  An object containing persistent data for the bivariate interpolation.
!-------------------------------------------------------------------------------
      TYPE :: bivariate_class
!>  Number of interpolation points.
         INTEGER                            :: nsu
!>  Number of interpolation points.
         INTEGER                            :: nrz
!>  Lower left interpolation point.
         INTEGER, DIMENSION(:), POINTER     :: irz11_bi => null()
!>  Lower right interpolation point.
         INTEGER, DIMENSION(:), POINTER     :: irz12_bi => null()
!>  Upper left interpolation point.
         INTEGER, DIMENSION(:), POINTER     :: irz21_bi => null()
!>  Upper right interpolation point.
         INTEGER, DIMENSION(:), POINTER     :: irz22_bi => null()
!>  Lower left interpolation weight.
         REAL(rprec), DIMENSION(:), POINTER :: w11_bi => null()
!>  Lower right interpolation weight.
         REAL(rprec), DIMENSION(:), POINTER :: w12_bi => null()
!>  Upper left interpolation weight.
         REAL(rprec), DIMENSION(:), POINTER :: w21_bi => null()
!>  Upper right interpolation weight.
         REAL(rprec), DIMENSION(:), POINTER :: w22_bi => null()
      CONTAINS
         PROCEDURE :: set_grids_1d => bivariate_set_grids_1d
         PROCEDURE :: set_grids_2d => bivariate_set_grids_2d
         GENERIC   :: set_grids => set_grids_1d, set_grids_2d
         PROCEDURE :: get_4pt => bivariate_get_4pt
         FINAL     :: bivariate_destruct
      END TYPE

*******************************************************************************
!  INTERFACE BLOCKS
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  Interface for bivariate constructor.
!-------------------------------------------------------------------------------
      INTERFACE bivariate_class
         MODULE PROCEDURE bivariate_construct
      END INTERFACE

      CONTAINS
!*******************************************************************************
!  CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Construct a @ref bivariate_type object.
!>
!>  Allocates memory and initializes a @ref bivariate_type object.
!>
!>  @param[in] ns Number of radial grid points.
!>  @param[in] nu Number of poloidal grid points.
!>  @returns A pointer to a constructed @ref bivariate_type object.
!-------------------------------------------------------------------------------
      FUNCTION bivariate_construct(ns, nu)

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bivariate_class), POINTER :: bivariate_construct
      INTEGER, INTENT(in)              :: ns
      INTEGER, INTENT(in)              :: nu

!  local variables
      REAL (rprec)                     :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      ALLOCATE(bivariate_construct)

      bivariate_construct%nsu = ns*nu

      ALLOCATE (bivariate_construct%irz11_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%irz12_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%irz21_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%irz22_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%w11_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%w12_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%w21_bi(bivariate_construct%nsu))
      ALLOCATE (bivariate_construct%w22_bi(bivariate_construct%nsu))

      CALL profiler_set_stop_time('bivariate_construct', start_time)

      END FUNCTION

!*******************************************************************************
!  DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Deconstruct a @ref bivariate_type object.
!>
!>  Deallocates memory and uninitializes a @ref bivariate_type object.
!>
!>  @param[inout] this A @ref bivariate_type instance.
!-------------------------------------------------------------------------------
      SUBROUTINE bivariate_destruct(this)

      IMPLICIT NONE

!  Declare Arguments
      TYPE (bivariate_class), INTENT(inout) :: this

!  Start of executable code
      IF (ASSOCIATED(this%irz11_bi)) THEN
         DEALLOCATE (this%irz11_bi)
         this%irz11_bi => null()
      END IF

      IF (ASSOCIATED(this%irz12_bi)) THEN
         DEALLOCATE (this%irz12_bi)
         this%irz12_bi => null()
      END IF

      IF (ASSOCIATED(this%irz21_bi)) THEN
         DEALLOCATE (this%irz21_bi)
         this%irz21_bi => null()
      END IF

      IF (ASSOCIATED(this%irz22_bi)) THEN
         DEALLOCATE (this%irz22_bi)
         this%irz22_bi => null()
      END IF

      IF (ASSOCIATED(this%w11_bi)) THEN
         DEALLOCATE (this%w11_bi)
         this%w11_bi => null()
      END IF

      IF (ASSOCIATED(this%w12_bi)) THEN
         DEALLOCATE (this%w12_bi)
         this%w12_bi => null()
      END IF

      IF (ASSOCIATED(this%w21_bi)) THEN
         DEALLOCATE (this%w21_bi)
         this%w21_bi => null()
      END IF

      IF (ASSOCIATED(this%w22_bi)) THEN
         DEALLOCATE (this%w22_bi)
         this%w22_bi => null()
      END IF

      END SUBROUTINE

!*******************************************************************************
!  SETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Set up the interpolation grid for a phi plane.
!>
!>  Compute the indices (ir_bi, jz_bi) and weight factors (wij_bi) for
!>  performing bivariate (4 pt) interpolation FROM a rectangular (R X Z) grid TO
!>  a general (non-orthogonal) grid (s, u). R and Z grids may have different
!>  mesh sizes.
!>
!>  @param[inout] this  A @ref bivariate_type instance.
!>  @param[in]    rsu   r coordinate at s,u
!>  @param[in]    zsu   z coordinate at s,u
!>  @param[in]    rgrid R coordinate on a fixed, equally spaced grid
!>  @param[in]    zgrid Z coordinate on a fixed, equally spaced grid
!-------------------------------------------------------------------------------
      SUBROUTINE bivariate_set_grids_2d(this, rsu, zsu, rgrid, zgrid)

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bivariate_class), INTENT(inout)   :: this
      REAL (rprec), DIMENSION(:,:), INTENT(in) :: rsu
      REAL (rprec), DIMENSION(:,:), INTENT(in) :: zsu
      REAL (rprec), DIMENSION(:), INTENT(in)   :: rgrid
      REAL (rprec), DIMENSION(:), INTENT(in)   :: zgrid

!  local variables
      INTEGER                                  :: ku
      INTEGER                                  :: js
      INTEGER                                  :: ns
      INTEGER                                  :: nu
      INTEGER                                  :: nr
      INTEGER                                  :: nz
      INTEGER                                  :: index1d
      INTEGER                                  :: ir
      INTEGER                                  :: jz
      INTEGER                                  :: ir1
      INTEGER                                  :: jz1
      REAL (rprec)                             :: rad0
      REAL (rprec)                             :: zee0
      REAL (rprec)                             :: ri
      REAL (rprec)                             :: zj
      REAL (rprec)                             :: pr
      REAL (rprec)                             :: qz
      REAL (rprec)                             :: temp
      REAL (rprec)                             :: delr
      REAL (rprec)                             :: delz
      REAL (rprec)                             :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      nr = SIZE(rgrid, 1)
      nz = SIZE(zgrid, 1)
      this%nrz = nr*nz
      ns = SIZE(rsu, 1)
      nu = SIZE(rsu, 2)

      delr = rgrid(2) - rgrid(1)
      delz = zgrid(2) - zgrid(1)

      index1d = 0

      DO ku = 1, nu
         DO js = 1, ns

            index1d = index1d + 1

!  Check that boundary points are inside grid. Stop if not.
            rad0 = rsu(js, ku)
            zee0 = zsu(js, ku)

            IF (rad0.lt.rgrid(1) .or. rad0.gt.rgrid(nr) .or.                   &
     &          zee0.lt.zgrid(1) .or. zee0.gt.zgrid(nz)) THEN
               STOP 'Plasma point outside response function grid!'
            END IF

!  Detremine integer indices (ir, jz) for the lower left R, Z corner grid point.
            ir = INT((rad0 - rgrid(1))/delr) + 1
            jz = INT((zee0 - zgrid(1))/delz) + 1
            ir1 = MIN(nr, ir + 1)
            jz1 = MIN(nz, jz + 1)

!  Store indices in 1d attays.
            this%irz11_bi(index1d) = ir  + nr*(jz - 1)
            this%irz22_bi(index1d) = ir1 + nr*(jz1 - 1)
            this%irz12_bi(index1d) = ir  + nr*(jz1 - 1)
            this%irz21_bi(index1d) = ir1 + nr*(jz - 1)

!  Compute (ri, zj) and (pr, qz) at grid point (ir, jz) also compute weights
!  for the 4 corner grid points.
            ri = rgrid(ir)
            zj = zgrid(jz)
            pr = (rad0 - ri)/delr
            qz = (zee0 - zj)/delz
            temp = pr*qz
            this%w22_bi(index1d) = temp                   !p*q
            this%w21_bi(index1d) = pr - temp              !p*(1-q)
            this%w12_bi(index1d) = qz - temp              !q*(1-p)
            this%w11_bi(index1d) = 1 + temp - (pr + qz)   !(1-p)*(1-q)

         END DO
      END DO

      CALL profiler_set_stop_time('bivariate_set_grids_2d', start_time)

      END SUBROUTINE

!-------------------------------------------------------------------------------
!>  @brief Set up the interpolation grid for a phi plane.
!>
!>  Compute the indices (ir_bi, jz_bi) and weight factors (wij_bi) for
!>  performing bivariate (4 pt) interpolation FROM a rectangular (R X Z) grid TO
!>  a general (non-orthogonal) grid (s, u). R and Z grids may have different
!>  mesh sizes.
!>
!>  @param[inout] this  A @ref bivariate_type instance.
!>  @param[in]    rsu   r coordinate at s,u
!>  @param[in]    zsu   z coordinate at s,u
!>  @param[in]    rgrid R coordinate on a fixed, equally spaced grid
!>  @param[in]    zgrid Z coordinate on a fixed, equally spaced grid
!-------------------------------------------------------------------------------
      SUBROUTINE bivariate_set_grids_1d(this, rsu, zsu, rgrid, zgrid)

      IMPLICIT NONE

!  Declare Arguments
      CLASS (bivariate_class), INTENT(inout) :: this
      REAL (rprec), DIMENSION(:), INTENT(in) :: rsu
      REAL (rprec), DIMENSION(:), INTENT(in) :: zsu
      REAL (rprec), DIMENSION(:), INTENT(in) :: rgrid
      REAL (rprec), DIMENSION(:), INTENT(in) :: zgrid

!  local variables
      INTEGER                                :: ku
      INTEGER                                :: ns
      INTEGER                                :: nu
      INTEGER                                :: nr
      INTEGER                                :: nz
      INTEGER                                :: index1d
      INTEGER                                :: ir
      INTEGER                                :: jz
      INTEGER                                :: ir1
      INTEGER                                :: jz1
      REAL (rprec)                           :: rad0
      REAL (rprec)                           :: zee0
      REAL (rprec)                           :: ri
      REAL (rprec)                           :: zj
      REAL (rprec)                           :: pr
      REAL (rprec)                           :: qz
      REAL (rprec)                           :: temp
      REAL (rprec)                           :: delr
      REAL (rprec)                           :: delz
      REAL (rprec)                           :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      nr = SIZE(rgrid, 1)
      nz = SIZE(zgrid, 1)
      this%nrz = nr*nz
      ns = 1
      nu = SIZE(rsu)

      delr = rgrid(2) - rgrid(1)
      delz = zgrid(2) - zgrid(1)

      index1d = 0

      DO ku = 1, nu
         index1d = index1d + 1
!
!       CHECK THAT BOUNDARY POINTS ARE INSIDE GRID.  IF NOT, STOP!
!
         rad0 = rsu(ku)
         zee0 = zsu(ku)

         IF (rad0.lt.rgrid(1) .or. rad0.gt.rgrid(nr) .or.                      &
     &       zee0.lt.zgrid(1) .or. zee0.gt.zgrid(nz)) THEN
            STOP 'Plasma point outside response function grid!'
         END IF
!
!       DETERMINE INTEGER INDICES (IR,JZ) FOR LOWER LEFT R, Z CORNER GRID POINT
!
         ir = INT((rad0 - rgrid(1))/delr) + 1
         jz = INT((zee0 - zgrid(1))/delz) + 1
         ir1 = MIN(nr, ir + 1)
         jz1 = MIN(nz, jz + 1)

!
!        STORE INDICES IN 1D ARRAYS
!
         this%irz11_bi(index1d) = ir  + nr*(jz - 1)
         this%irz22_bi(index1d) = ir1 + nr*(jz1 - 1)
         this%irz12_bi(index1d) = ir  + nr*(jz1 - 1)
         this%irz21_bi(index1d) = ir1 + nr*(jz - 1)
!
!       COMPUTE RI, ZJ AND PR , QZ AT GRID POINT (IR , JZ)
!       ALSO, COMPUTE WEIGHTS WIJ FOR 4 CORNER GRID POINTS
!
         ri = rgrid(ir)
         zj = zgrid(jz)
         pr = (rad0 - ri)/delr
         qz = (zee0 - zj)/delz
         temp = pr*qz
         this%w22_bi(index1d) = temp                   !p*q
         this%w21_bi(index1d) = pr - temp              !p*(1-q)
         this%w12_bi(index1d) = qz - temp              !q*(1-p)
         this%w11_bi(index1d) = 1 + temp - (pr + qz)   !(1-p)*(1-q)
      END DO

      CALL profiler_set_stop_time('bivariate_set_grids_2d', start_time)

      END SUBROUTINE

!*******************************************************************************
!  GETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Interpolate points on to response function grid.
!>
!>  Determin response function on the s,u grid by 4 point bivariate
!>  interpolation of response function on the r,z mesh at a fixed toroidal
!>  plane. This uses a 2-D interpolation based on the four point formula in
!>  Abramowitz and Stegun, Eq. 25.2.66. Note that the arrays need an explicit
!>  size in dimension to allow the 2D arrays to be passed to the 1D arguments.
!>
!>  @param[inout] this    A @ref bivariate_type instance.
!>  @param[in]    resp_rz Response function on the rz grid.
!>  @param[out]   resp_su Response function on the su grid.
!-------------------------------------------------------------------------------
      SUBROUTINE bivariate_get_4pt(this, resp_rz, resp_su)
      IMPLICIT NONE
      CLASS (bivariate_class), INTENT(inout)         :: this
      REAL (rprec), DIMENSION(this%nrz), INTENT(in)  :: resp_rz
      REAL (rprec), DIMENSION(this%nsu), INTENT(out) :: resp_su

!  local variables
      INTEGER                                        :: jsu
      REAL (rprec)                                   :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      DO jsu = 1, this%nsu

!  Compute repsonse at R, Phi (fixed kin), Z by 4 point interpolation.
         resp_su(jsu) = this%w11_bi(jsu)*resp_rz(this%irz11_bi(jsu))           &
     &                + this%w22_bi(jsu)*resp_rz(this%irz22_bi(jsu))           &
     &                + this%w21_bi(jsu)*resp_rz(this%irz21_bi(jsu))           &
     &                + this%w12_bi(jsu)*resp_rz(this%irz12_bi(jsu))

      END DO

      CALL profiler_set_stop_time('bivariate_get_4pt', start_time)

      END SUBROUTINE

      END MODULE


submodule (mod_integration) smod_integration
    implicit none

contains

    module subroutine gaussquad(N, x, w)
        integer,  intent(in)                :: N
        real(wp), intent(out), dimension(N) :: x, w

        select case (N)
        case (1)
            x = [ 0._wp ]
            w = [ 2._wp ]

        case (2)
            x = [ -0.5773502691896257_wp, &
                0.5773502691896257_wp ]
            w = [ 1._wp, &
                1._wp ]

        case (3)
            x = [ -0.7745966692414834_wp, &
                0._wp, &
                0.7745966692414834_wp ]
            w = [ 0.5555555555555556_wp, &
                0.8888888888888888_wp, &
                0.5555555555555556_wp ]

        case (4)
            x = [ -0.8611363115940526_wp, &
                -0.3399810435848563_wp, &
                0.3399810435848563_wp, &
                0.8611363115940526_wp ]
            w = [ 0.3478548451374538_wp, &
                0.6521451548625461_wp, &
                0.6521451548625461_wp, &
                0.3478548451374538_wp ]

        case (5)
            x = [ -0.9061798459386640_wp, &
                -0.5384693101056831_wp, &
                0._wp, &
                0.5384693101056831_wp, &
                0.9061798459386640_wp ]
            w = [ 0.2369268850561891_wp, &
                0.4786286704993665_wp, &
                0.5688888888888889_wp, &
                0.4786286704993665_wp, &
                0.2369268850561891_wp ]

        case default

            ! call lgwt(-1._wp, 1._wp, N, x, w)
            ! call cgwt(N, x, w)
            call gaussquad_rosetta(N, x, w)

        end select

        return
    end subroutine gaussquad

    subroutine lgwt(a, b, num_pts, x, w)

        !*  This function is a fortran90 port of the matlab function,
        !   lgwt.m The source code of lgwt.m was originally found at:
        !   http://www.mathworks.com/matlabcentral/fileexchange/4540

        ! Variables in/out
        integer, intent(in)   :: num_pts
        real(wp), intent(in)  :: a, b
        real(wp), intent(out) :: x(:), w(:)

        ! Local variables
        integer:: ii, jj, N, N1, N2
        ! real(wp), parameter:: eps=sqrt(epsilon(1.0_wp))
        real(wp), parameter                       :: eps=1d-10
        real(wp), dimension(:),     allocatable   :: xu, array1, y, y0, Lpp
        real(wp), dimension(:, :),  allocatable   :: L, Lp
        real(wp), parameter                       :: pi = 4.0_wp*datan(1.0_wp)

        N = num_pts - 1
        N1 = N + 1
        N2 = N + 2

        ! Allocate and initialize arrays
        allocate(xu(N1), array1(N1), y(N1), y0(N1))
        allocate(L(N1, N2), Lp(N1, N2), Lpp(N1))

        L = 0.0_wp
        Lp = 0.0_wp

        call linspace(-1.0_wp, 1.0_wp, xu)

        array1 = [ (ii, ii = 0, N) ]
        ! Initial guess of the roots of the Legendre polynomial of order N
        y = cos((2.0_wp * dble(array1) + 1.0_wp) * &
            pi / (2.0_wp * dble(N) + 2.0_wp)) + &
            (0.27_wp/dble(N1)) * sin(pi*xu*dble(N)/dble(N2))

        y0 = 2.0_wp

        outer: do
            L(:, 1) = 1.0_wp
            Lp(:, 1) = 0.0_wp

            L(:, 2) = y
            Lp(:, 2) = 1.0_wp

            inner: do jj = 2, N1
                L(:, jj+1) = ( &
                    (2.0_wp*dble(jj)-1.0_wp) * y * L(:, jj) - &
                    dble(jj-1) * L(:, jj-1)) &
                    / dble(jj)
            enddo inner

            Lpp = dble(N2) * (L(:,N1) - y * L(:,N2)) &
                / (1.0_wp - y**2.0_wp)

            y0 = y
            y = y0 - L(:,N2)/Lpp

            if ( norm2(y-y0) < eps ) then
                exit outer
            endif
        enddo outer

        x = ( a*(1.0_wp-y) + b*(1.0_wp+y) ) / 2.0_wp
        w = ( b-a ) / ((1.0_wp - y**2.0_wp)*Lpp**2.0_wp) * &
            (dble(N2) / dble(N1))**2.0_wp

        return
    end subroutine lgwt

    subroutine cgwt(num_pts, x, w)

        ! This function determines the points and weights associated
        ! with Chebyshev-Gauss quadrature

        ! Variables in/out
        integer, intent(in) :: num_pts
        real(wp), intent(out), dimension(:) :: x, w

        ! Local variables
        integer:: ii
        real(wp), parameter:: pi = 4._wp*datan(1._wp)

        x = cos( (dble(2*[( ii, ii=1,num_pts )]-1) )/dble(2*num_pts) * pi)

        w = pi/dble(num_pts) / ((1.0_wp - x**2._wp)**(-0.5_wp))

        ! print*, x
        ! print*, w
        ! stop
        return
    end subroutine cgwt

    subroutine gaussquad_rosetta(n, r1, r2)

        ! This code was originally found at the following website:
        !  http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature#Fortran

        integer,  intent(in)    :: n
        real(wp), intent(out)   :: r1(n), r2(n)

        integer                 :: k
        real(wp), parameter     :: pi = 4._wp*atan(1._wp)
        real(wp)                :: x, f, df, dx
        integer                 :: i,  iter
        real(wp), allocatable   :: p0(:), p1(:), tmp(:)

        p0 = [1._wp]
        p1 = [1._wp, 0._wp]

        do k = 2, n
            tmp = ((2*k-1)*[p1,0._wp]-(k-1)*[0._wp, 0._wp,p0])/k
            p0 = p1; p1 = tmp
        enddo

        outer: do i = 1, n
            x = cos(pi*(i-0.25_wp)/(n+0.5_wp))
            inner: do iter = 1, 10
                f = p1(1); df = 0._wp
                deep: do k = 2, size(p1)
                    df = f + x*df
                    f  = p1(k) + x * f
                enddo deep
                dx =  f / df
                x = x - dx
                if (abs(dx)<10*epsilon(dx)) exit
            enddo inner
            r1(i) = x
            r2(i) = 2/((1-x**2)*df**2)
        enddo outer
        return
    end subroutine gaussquad_rosetta

end submodule


text:
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
      SUBROUTINE RESSYM ( INPUT, OUTPUT )

      IMPLICIT NONE

      LOGICAL               FAILED
      INTEGER               LASTNB

      CHARACTER*(*)         INPUT
      CHARACTER*(*)         OUTPUT

      INTEGER               LNSIZE
      PARAMETER           ( LNSIZE = 255 )

      CHARACTER*(LNSIZE)    SUBSTR

      CHARACTER*(1)         TAB
      CHARACTER*(1)         SPACE
      CHARACTER*(1)         EQUOTE

      INTEGER               I
      INTEGER               LOC
      INTEGER               R
      INTEGER               E

      LOGICAL               CHANGE
      LOGICAL               TRAN1
      LOGICAL               TRAN2

      CALL CHKIN ( 'RESSYM' )

      TAB    = CHAR(9)
      SPACE  = ' '

      CALL GETEQ ( EQUOTE )

      CALL REPLCH ( INPUT, TAB, SPACE, OUTPUT )

      CALL PRTRAP ( OUTPUT, CHANGE )

C
C     Now we just loop until all translations have
C     been performed.  We do:
C
C        1) symbol resolution
C        2) query resolution
C        3) tab removal
C
      DO WHILE ( CHANGE )

         CHANGE = .FALSE.
         TRAN1  = .TRUE.
         TRAN2  = .TRUE.
C
C        First we resolve all symbols.  After each pass we check
C        that we have not created a command that must be trapped.
C
         DO WHILE ( TRAN1 .AND. TRAN2 )

            CALL STRAN  ( OUTPUT, OUTPUT, TRAN1 )
            CALL PRTRAP ( OUTPUT,         TRAN2        )
C
C           Determine whether or not more changes are possible
C           at this point.
C
            CHANGE =          ( CHANGE .OR. TRAN1 )
     .               .AND. TRAN2
     .               .AND. .NOT. FAILED()

         END DO
C
C        If we don't have any errors we take a stab at replacing
C        all queries.  Note that queries can not result in changing
C        anything that isn't a query so we don't have to trap
C        inside the loop.  Note that this means you can't have
C        a command like DEFINE? SYMBOL? VALUE? and just replace
C        the first two queries.  You've got to do them all.  If
C        you want a symbol to have a query you must do it this
C        way:  DEFINE SYMBOL  QUERY?  That way the queries won't
C        get resolve too soon.
C
C        Note:  This can easily be changed so that if a query
C        introduces a symbol, we immediately loop back to the
C        symbol resolution branch.  Simply change the DO WHILE
C        loop below to an IF.  The "loop" will then terminate
C        after one execution leaving any remaining queries
C        untouched until the next pass through the loop.
C
         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'RESSYM' )
            RETURN
         END IF

         TRAN1 = .NOT. FAILED()

         DO WHILE ( TRAN1 )

            CALL QTRAN  ( OUTPUT, OUTPUT, TRAN1 )
            CALL REPLCH ( OUTPUT, TAB, SPACE, OUTPUT )

            CHANGE = CHANGE .OR. TRAN1

         END DO

         CALL PRTRAP ( OUTPUT, TRAN2 )
         CHANGE = CHANGE .AND. TRAN2

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'RESSYM' )
            RETURN
         END IF

      END DO

      IF ( TRAN2 ) THEN
C
C        We remove the special markers that may have been present to
C        protect symbol or query resolution.
C
         I = 1
         CALL NTHUQW ( OUTPUT, I, ' ', SUBSTR, LOC )

         DO WHILE ( LOC .GT. 0 )

            R = LASTNB ( SUBSTR ) - 1
            E = LOC + R

            CALL REPLCH( OUTPUT(LOC:E), EQUOTE, SPACE, OUTPUT(LOC:E) )

            I = I+1
            CALL NTHUQW ( OUTPUT, I, ' ', SUBSTR, LOC )

         END DO

      END IF
C
C     Finally, left justify the commmand.
C
      CALL LJUST ( OUTPUT, OUTPUT )

      CALL CHKOUT ( 'RESSYM' )
      RETURN
      END

ext: f

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
module type_pseudo_timestep
    use mod_kinds,          only: ik, rk
    use type_mesh,      only: mesh_t
    use type_solverdata,    only: solverdata_t
    use type_properties,    only: properties_t
    implicit none

    !>  A class that handles computing a pseudo timestep for Quasi-Newton
    !!  methods.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   11/17/2016
    !!
    !------------------------------------------------------------------------
    type, abstract, public :: pseudo_timestep_t

    contains

        procedure(compute_interface), deferred :: compute

    end type pseudo_timestep_t
    !************************************************************************

    abstract interface
        subroutine compute_interface(self,idomain,mesh,prop,sdata,cfl,itime)
            import ik, rk
            import pseudo_timestep_t
            import mesh_t
            import properties_t
            import solverdata_t
            class(pseudo_timestep_t),   intent(in)      :: self
            integer(ik),                intent(in)      :: idomain
            type(mesh_t),           intent(inout)   :: mesh
            type(properties_t),         intent(in)      :: prop
            type(solverdata_t),         intent(inout)   :: sdata
            real(rk),                   intent(in)      :: cfl(:)
            integer(ik),                intent(in)      :: itime
        end subroutine
    end interface

    !>  Default pseudo_timestep implementation.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   11/17/2016
    !!
    !!
    !-------------------------------------------------------------------------
    type, extends(pseudo_timestep_t), public :: default_pseudo_timestep_t

    contains

        procedure   :: compute

    end type default_pseudo_timestep_t
    !*************************************************************************

contains

    !>  Default procedure for computing a pseudo timestep.
    !!
    !!  Default is based strictly
    !!
    !!  @author Nathan A. Wukie
    !!
    !------------------------------------------------------------------------
    subroutine compute(self,idomain,mesh,prop,sdata,cfl,itime)
        class(default_pseudo_timestep_t),   intent(in)      :: self
        integer(ik),                        intent(in)      :: idomain
        type(mesh_t),                       intent(inout)   :: mesh
        type(properties_t),                 intent(in)      :: prop
        type(solverdata_t),                 intent(inout)   :: sdata
        real(rk),                           intent(in)      :: cfl(:)
        integer(ik),                        intent(in)      :: itime

        integer(ik) :: ielem, ieqn
        real(rk)    ::  h

        !
        ! Loop through elements and compute time-step function
        !
       do ielem = 1,mesh%domain(idomain)%nelem

            !
            ! Compute element spacing parameter
            !
            h = mesh%domain(idomain)%elems(ielem)%vol**(1._rk / 3._rk)

            !
            ! Compute elemen-local timestep
            !
            if (size(cfl) == 1) then
                mesh%domain(idomain)%elems(ielem)%dtau = cfl*h
            else
                do ieqn = 1,size(cfl)
                    mesh%domain(idomain)%elems(ielem)%dtau(ieqn) = cfl(ieqn)*h
                end do
            end if

        end do  ! ielem

    end subroutine compute
    !***********************************************************************

end module type_pseudo_timestep

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
+incdir+../../../rtl/TOP
+incdir+../../../rtl/CPU
+incdir+../../../rtl/APB_DEV
+incdir+../../../rtl/APB_DEV/NAND
+incdir+../../../rtl/APB_DEV/URT
+incdir+../../../rtl/SPI
+incdir+../../../rtl/MIG
+incdir+../../../rtl/MAC
+incdir+../../../rtl/AMBA
+incdir+../../../rtl/CONFREG
+incdir+../../testbench
+incdir+../../testbench/vmac
+nclog_ncsim+ncsim.log
+access+rw
-sv
-64bit
-nowarn MRSTAR
+pathpulse
+profile
+define+x1Gb
+define+sg15E
+define+x16
+define+VIRTUAL_MAC
+define+_DUMPDUMP
-f sys.list
-f rtl.list
-f lib.list

max_stars_repo_name: Ewenwan/nontrivial-mips


text:
!
! Copyright (c) 2015, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!

          subroutine initialize(a)
          character*5 a(10,10)
          a(1,:) = "jbc"
          a(2,:) = "def"
          a(3,:) = "gh"
          a(4,:) = "ij"
          a(5,:) = "kj"
          a(6,:) = "abc"
          a(7,:) = "op"
          a(8,:) = "qr"
          a(9,:) = "st"
          a(10,:) = "ug"

          end subroutine

          program mymax
          parameter(N=86)
          integer result(86), expect(N)
          character*5 char_result(86)
          character*5 a(10,10)
          data expect /-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,&
                       &-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,&
                       &-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,&
                       &-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,&
                       &-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,&
                       &-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,&
                       &-1,0,-1,-1,-1,-1,-1,-1,-1,-1,&
                       &-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,&
                       &-1,0,-1,-1,-1,-1 /

          call initialize(a)

! test basic min/maxloc
! output "a", "mf"
          char_result(1:1) = minval((/'j','d','a'/))
          char_result(2:2) = maxval((/' dk','mf','ah'/))
          result(1) = (char_result(1) == 'a')
          result(2) = (char_result(2) == 'mf')

! test dim
          char_result(3:12) = minval(a, dim=1)
          char_result(13:22) = minval(a, dim=2)
          char_result(23:32) = maxval(a, dim=1)
          char_result(33:42) = maxval(a, dim=2)
          result(3:12) = (char_result(3:12) == a(6,:))
          result(13:22) = (char_result(13:22) == a(:,1))
          result(23:32) = (char_result(23:32) == a(10,:))
          result(33:42) = (char_result(33:42) == a(:,1))

! test mask
          char_result(43:44) = minval(a, mask = a > 'bc')
          char_result(45:46) = maxval(a, mask = a > 'br')
          result(43:44) = (char_result(43:44) == a(2,1:2) )
          result(45:46) = (char_result(45:46) == a(10,1:2))

! test mask with dim
          char_result(47:56) = minval(a, mask = a > 'bc', dim=1)
          char_result(57:66) = minval(a, mask = a > 'bc', dim=2)
          char_result(67:76) = maxval(a, mask = a > 'br', dim=1)
          char_result(77:86) = maxval(a, mask = a > 'br', dim=2)
          result(47:56) = (char_result(47:56) == a(2,:) )
          result(57:66) = (char_result(57:66) == a(:,1))
          result(67:76) = (char_result(67:76) == a(10,:) )
          result(77:86) = (char_result(77:86) == a(:,1))

          call check(result, expect, N)

          end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


! { dg-do compile }
!
! PR fortran/47569
!
! Contributed by Jos de Kloe
!
module teststr
  implicit none
  integer, parameter :: GRH_SIZE = 20, NMAX = 41624
  type strtype
    integer   :: size
    character :: mdr(NMAX)
  end type strtype
contains
  subroutine sub2(string,str_size)
    integer,intent(in)    :: str_size
    character,intent(out) :: string(str_size)
    string(:) = 'a'
  end subroutine sub2
  subroutine sub1(a)
    type(strtype),intent(inout) :: a
    call sub2(a%mdr(GRH_SIZE+1),a%size-GRH_SIZE)
  end subroutine sub1
end module teststr


text:
! Copyright 2011 Free Software Foundation, Inc.
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! 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/>.

program test
  integer :: foo (2, 3, 4)
  integer :: singledim (4)
  foo (:, :, :) = 10
  foo (2, 3, 4) = 20
  foo (2, 3, 4) = 20            ! break-static
  singledim (:) = 1
  singledim (4) = 2
  call sub (singledim, 4, singledim)
end
subroutine sub (varbound, n, unbound)
  integer :: n, varbound (n), unbound (*)
  varbound (4) = unbound (4)    ! break-variable
end

max_stars_repo_name: vidkidz/crossbridge

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
      program cadena_de_caracteres_asignable
        implicit none

        character(:), allocatable :: nombre
        character(:), allocatable :: apellido

        allocate(character(8) :: nombre)
        nombre='Diogenes'
        !Declaración de asignación explícita

        apellido='De Sinope'
        !Asignación en la asignación

        print *, nombre//' '//apellido
      end program cadena_de_caracteres_asignable

max_stars_repo_name: Marc-xyz/InicioRapidoEnFortran

max_issues_repo_name: Marc-xyz/InicioRapidoEnFortran

max_forks_repo_name: Marc-xyz/InicioRapidoEnFortran


text:
      SUBROUTINE BARTYP(KIND)
C
C          --------------------------------------------------
C          ROUTINE NO. (  89)    VERSION (A8.1)     10:DEC:85
C          --------------------------------------------------
C
C          THIS SELECTS HORIZONTAL OR VERTICAL BARS FOR BARCHARTS
C          AND HISTOGRAMS.
C
C
C          <KIND> = ZERO     FOR VERTICAL   BARS.
C          <KIND> = NON-ZERO FOR HORIZONTAL BARS.
C
C
      COMMON /T0BRTY/ IBRTYP
      COMMON /T0TRAC/ IPRINT
      COMMON /T0TRAI/ ITRAC1,ITRAC2,ITRAC3,ITRAC4
C
C
      CALL G3INIT(2)
C
      ITRAC1= KIND
      IF (IPRINT.EQ.1) CALL G0MESG(159,5)
C
      IBRTYP= KIND
      IF (IBRTYP.NE.0) IBRTYP= 1
C
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
    MODULE EFDCPROMOD
    ! ** AUTHOR: DH CHUNG

    USE GLOBALVARS
    USE INFOMOD

    IMPLICIT NONE

    CONTAINS

    SUBROUTINE EFDCPRO
    INTEGER(4)::NS,LL,I,J,L,K,K1
    REAL(4)::ANG1,ANG2,ANG,XLNUTME,YLTUTMN,CCUE,CCVE,CCUN,CCVN,DZPC,MAXTHICK,DEPTHMIN,TMP,TMPCOR
    REAL(4),ALLOCATABLE:: THICK(:),DZC(:,:)

    HFREOUT = 0

    CALL SCANEFDC
    CALL VARALLOC

    LIJ = 1
    XCC = BLK
    YCC = BLK

    CALL INPUT      !GET ZB(IC,JC)

    LAYK = MIN(KC,LAYK)

    PRINT *,'*** READING LXLY.INP'
    OPEN(1,FILE=LXLYFILE,STATUS='OLD',ERR=998)
    CALL SKIPCOM(1, 'C')
    CUE=0
    CVE=0
    CUN=0
    CVN=0
    XMIN = 1E20
    YMIN = 1E20
    DO LL=1,LVC
        READ(1,*,ERR=999)I,J,XCC(LIJ(I,J)),YCC(LIJ(I,J)),CCUE,CCVE,CCUN,CCVN !,TMPVAL,TMPCOR
        L=LIJ(I,J)
        ANG1=ATAN2(CCUN,CCUE)
        ANG2=ATAN2(-CCVE,CCVN)
        ANG=0.5*(ANG1+ANG2)
        IF(SIGN(1.,ANG1).NE.SIGN(1.,ANG2))THEN
            IF(ABS(ANG1).GT.(1.57).OR.ABS(ANG2).GT.(1.57)) THEN  !(>pi/2)
                ! *** HANDLE THE DISCONTINUITY AT THE 180 DEGREES ANGLE
                ANG = ANG + ACOS(-1.0)
            ENDIF
        END IF
        CUE(L)= COS(ANG)
        CVE(L)=-SIN(ANG)
        CUN(L)= SIN(ANG)
        CVN(L)= COS(ANG)
        XMIN = MIN(XMIN,XCC(LIJ(I,J)))
        YMIN = MIN(YMIN,YCC(LIJ(I,J)))
    ENDDO
    XMAX = MAXVAL(XCC)+1
    YMAX = MAXVAL(YCC)+1
    CLOSE(1)

    KSZ = 1

    IF( IGRIDV > 0 )THEN
        IF (IGRIDV == 2) SGZFILE = trim(OUTFOLDER)//'SGZLAYER.OUT'
        ! *** READ SGZ BOTTOM ACTIVE LAYER FOR SIGMA-ZED CELLS
        OPEN(1,FILE=SGZFILE,ACTION='READ',ERR=997)
        CALL SKIPCOM(1, 'C')
        DO LL=2,LA
            READ(1,*,END=1000)I,J,K
            L = LIJ(I,J)
            KSZ(L)=K  ! BED LAYER
        ENDDO
1000    CLOSE(1)
    ENDIF

    ! **
    DZCS = 0
    DO L=2,LA
        DZPC = SUM(DZCK(KSZ(L):KC))
        DO K=KSZ(L),KC
            DZCS(L,K) = DZCK(K)/DZPC
        ENDDO
    ENDDO

    RETURN
997 PRINT*, ' **** OPEN ERROR: SGZFILE '
    PAUSE(1)
    STOP
998 PRINT*, ' **** OPEN ERROR: LXLY.INP '
    PAUSE(1)
    STOP
999 PRINT*, ' **** READ ERROR: LXLY.INP '
    PAUSE(1)
    STOP
    END SUBROUTINE

    SUBROUTINE SCANEFDC
    INTEGER(4)::N,NS,IMDXDY,KDUM,K,STAT,ITMP,I,ISO,J,LINE,LL,L
    REAL(4)::DX,DY,DXYCVT,ZBRADJ,ZBRCVRT,HMIN,HADADJ,HCVRT,HDRY,HWET,RTMP
    CHARACTER(4)::NCRD

    PRINT*,'*** READING EFDC.INP'
    OPEN(1,FILE=EFDCFILE,STATUS='OLD',ERR=998)

    NCRD='1A'
    CALL SEEK('C1A')
    READ(1,*,ERR=10,END=30)ITMP,ITMP,IGRIDV !,IS2TIM,ISHOUSATONIC

    NCRD = 'C6'
    CALL SEEK('C6')
    DO N=0,8
        READ(1,*,ERR=10,END=30) ISTRAN(N)
    ENDDO

    NCRD = 'C9'
    CALL SEEK('C9')
    READ(1,*,ERR=10,END=30)IC,JC,LC,LVC,ISCLO,NDM,LDM,ISMASK,ISCONNECT
    IF(IC.GE.3)THEN
        ICM=IC+1
    ELSE
        PRINT*, 'IC MUST BE AT LEAST 3'
        PAUSE(1)
        STOP
    ENDIF
    IF(JC.GE.3)THEN
        JCM=JC+1
    ELSE
        PRINT*, 'IJ MUST BE AT LEAST 3'
        PAUSE(1)
        STOP
    ENDIF
    IF(LC.GE.3)THEN
        LCM=LC+1
    ELSE
        PRINT*, 'LC MUST BE AT LEAST 3'
        PAUSE(1)
        STOP
    ENDIF

    NCRD='C9A'
    CALL SEEK('C9A')
    READ(1,*,ERR=10,END=30)KC !,KSIG
    KC=ABS(KC)
    IF(KC.GE.1)THEN
        KCM=KC+1
    ELSE
        PRINT*,'KC MUST BE AT LEAST 1'
        PAUSE(1)
        STOP
    ENDIF

    KSM=KCM
    KS=MAX(KC-1,1)

    ALLOCATE(DZCK(KC))
    NCRD='10'
    CALL SEEK('C10')
    DO K=1,KC
        READ(1,*,ERR=10,END=30)KDUM,DZCK(K)
    ENDDO

    NCRD = 'C11'
    CALL SEEK('C11')
    READ(1,*,ERR=10,END=30) DX,DY,DXYCVT,IMDXDY,ZBRADJ,ZBRCVRT,HMIN,HADADJ,HCVRT,HDRY,HWET,BELADJ,BELCVRT

    NCRD = 'C14'
    CALL SEEK('C14')
    READ(1,*,ERR=10,END=30) MTIDE,NWSER,NASER,ISGWIT,ISCHAN,ISWAVE

    NCRD='16'
    CALL SEEK('C16')
    READ(1,*,ERR=10,END=30) NPBS,NPBW,NPBE,NPBN !,NPFOR,NPFORT,NPSER,PDGINIT

    NCRD = 'C22'
    N= SEEK_SCANLINE('C22','NDYE')
    IF(N.GT.0) THEN
        READ(1,*,ERR=10,END=30)NDYE,NTOX,NSED,NSND
    ELSE
        READ(1,*,ERR=10,END=30)NTOX,NSED,NSND
    ENDIF
    NTXM=MAX(1,NTOX)
    NSCM=MAX(1,NSED)
    NSNM=MAX(1,NSND)
    NSTM=MAX(3,NSCM+NSNM+NTXM)

    IF(ISTRAN(8).GT.0)THEN
        LCMWQ=LCM
    ELSE
        LCMWQ=1
    ENDIF

    !C23*  READ VELOCITY, VOL SOUR/SINK, FLOW CONTROL, & WITHDRAW/RETURN DATA
    NCRD='23'
    CALL SEEK('C23')
    !READ(1,*,ERR=10,END=30) RTMP,RTMP,RTMP,NQCTL
    READ(1,*,ERR=10,END=30) ITMP,ITMP,ITMP,NQCTL,ITMP,ITMP,NQWR,ITMP,ITMP,NQCTLSER,NQCRULES
    !NQSIJ NQJPIJ   NQSER   NQCTL  NQCTLT  NHYDST    NQWR  NQWRSR   ISDIQ   NQCTLSER  NQCRULES

    ! *** SET KB AND CHECK FOR SEDZLJ USEAGE
    NSEDFLUME=0
    LSEDZLJ = .FALSE.
    IF(NSED.GT.0.OR.NSND.GT.0)THEN
        CALL SEEK('C36')
        READ(1,*,ERR=10,END=30)ISEDINT,ISEDBINT,ISEDWC,ISMUD,ISNDWC,ISEDVW,ISNDVW,KB,ISDTXBUG
        IF(KB.GE.1)THEN
            KBM=KB+1
        ELSE
            STOP 'KB MUST BE AT LEAST 1'
        ENDIF
        IF(ISTRAN(6)>0)THEN
            IF(ISEDWC==98)THEN
                LSEDZLJ = .TRUE.
                NSEDFLUME=1
            ELSEIF(ISEDWC==99)THEN
                LSEDZLJ = .TRUE.
                NSEDFLUME=2
            ENDIF
        ENDIF

        IF (ISTRAN(6)>0.OR.ISTRAN(7)>0) THEN
            NCRD = 'C36A'
            CALL SEEK('C36A')
            READ(1,*,ERR=10,END=30)ISBEDSTR
        ENDIF

        ! *** LOOK FOR LEGACY APPROACH TO SET SEDZLJ, I.E. IWRSP(1)=98
        IF(.NOT.LSEDZLJ .AND. ISTRAN(6)>0)THEN
            ISDTXBUG=0
            NCRD = 'C40'
            CALL SEEK('C40')
            READ(1,*,ERR=10,END=30)ISDTXBUG
            IF(ISDTXBUG==98)THEN
                LSEDZLJ = .TRUE.
                NSEDFLUME=1
            ENDIF
        ENDIF
    ELSE
        KBM=1
    ENDIF

    NCRD = 'C42A'
    CALL SEEK('C42A')
    DO NS=1,NSND
        READ(1,*,ERR=10,END=30) ISBDLDBC
    ENDDO

    NCRD='46A'
    CALL SEEK('C46A')
    READ(1,*,ERR=10,END=30)ISICE !,ISICECOV,ISICETHK,NISER,ICETHKFUN,DYICEBEG,DYICEEND,DYICEM1,DYICEM2,RICETHKMAX,TEMPICE

    K=0
    EFDCVER = 730
    TBEDIT = 0
    DO N=1,1000
        READ(1,*,ERR=10,END=30) NCRD
        IF (NCRD(1:4)=='C46D'.AND.K==0) THEN
            K=1
        ELSEIF (NCRD(1:4)=='C46D'.AND.K==1) THEN
            READ(1,*,ERR=10,END=30) ITMP,(RTMP,I=1,6),TBEDIT
            EXIT
        ELSEIF(NCRD(1:3)=='C47') THEN
            EFDCVER = 720
            EXIT
        ENDIF
    ENDDO

    NCRD='67'
    CALL SEEK('C67')
    READ(1,*,ERR=10,END=30) ISPD  !,NPD,NPDRT,NWPD,ISLRPD,ILRPD1,ILRPD2,JLRPD1, JLRPD2, MLRPDRT,IPLRPD
    NPD = 0

    NCRD='71A'
    CALL SEEK('C71A')
    READ(1,*,ERR=10,END=30) ITMP,ISBEXP !,NPBPH

    NCRD = 'C88'
    CALL SEEK('C88')
    READ(1,*,ERR=10,END=30) HFREOUT

    CLOSE(1)

    ! ** MAPPING
    !IF(ISCONNECT.GE.1)THEN
    !  !PRINT *,'***   READING MAPPGNS.INP'
    !  OPEN(1,FILE=MAPGFILE,STATUS='UNKNOWN',ERR=999)
    !  DO I=1,8
    !    READ(1,*)
    !  ENDDO
    !  READ(1,*,IOSTAT=ISO) NPNSBP
    !  IF(ISO.GT.0) STOP 'READ ERROR FOR FILE MAPPGNS.INP'
    !  !DO NPNS=1,NPNSBP
    !  !  READ(1,*,IOSTAT=ISO) ISPNS(NPNS),JSPNS(NPNS),INPNS(NPNS),JNPNS(NPNS)
    !  !  IF(ISO.GT.0) STOP 'READ ERROR FOR FILE MAPPGNS.INP'
    !  !ENDDO
    !  CLOSE(1)
    !ENDIF

    ! **
    IF( ISGWIT == 2 )THEN
        WRITE(*,'(A)')'READING GWSER.INP'
        OPEN(1,FILE='gwser.inp',STATUS='UNKNOWN')
        ! *** SKIP OVER TITLE AND AND HEADER LINES
        CALL SKIPCOM(1, '*')
        READ(1,*) NGWSER
        CLOSE(1)
    ENDIF

    RETURN
10  WRITE(*,'(A)') 'EFDC.INP: ERROR READING '//NCRD
    PAUSE(1)
    STOP
30  WRITE(*,'(A)') 'UNEXPECTED END OF INPUT FILE'
    PAUSE(1)
    STOP
998 PRINT*, 'OPEN ERROR: EFDC.INP'
    PAUSE(1)
    STOP
999 PRINT*, 'OPEN ERROR: MAPGFILE.INP'
    PAUSE(1)
    STOP
    END SUBROUTINE

    SUBROUTINE VARALLOC
    ALLOCATE(IJCT(ICM,JCM))
    ALLOCATE(IL(LCM))
    ALLOCATE(JL(LCM))
    ALLOCATE(LNC(LCM),LSC(LCM),LEC(LCM),LWC(LCM))
    ALLOCATE(LNEC(LCM),LNWC(LCM),LSEC(LCM),LSWC(LCM))

    ALLOCATE(LIJ(ICM,JCM))
    ALLOCATE(LCT(LCM))

    ALLOCATE(CUE(LCM))
    ALLOCATE(CUN(LCM))
    ALLOCATE(CVE(LCM))
    ALLOCATE(CVN(LCM))
    ALLOCATE(XCC(LCM))
    ALLOCATE(YCC(LCM))
    ALLOCATE(BELV(LCM))
    !ALLOCATE(HP(LCM))
    ALLOCATE(HP0(LCM))

    ALLOCATE(RSSBCE(LCM))
    ALLOCATE(RSSBCN(LCM))
    ALLOCATE(RSSBCS(LCM))
    ALLOCATE(RSSBCW(LCM))

    !ALLOCATE(ISPNS(NPNSBP))
    !ALLOCATE(JSPNS(NPNSBP))
    !ALLOCATE(INPNS(NPNSBP))
    !ALLOCATE(JNPNS(NPNSBP))

    ALLOCATE(U(LCM,KCM))
    ALLOCATE(V(LCM,KCM))
    ALLOCATE(W(LCM,0:KCM))
    ALLOCATE(UK(LCM,KCM))
    ALLOCATE(VK(LCM,KCM))

    ALLOCATE(UA(LCM))
    ALLOCATE(VA(LCM))
    ALLOCATE(WA(LCM))

    !ALLOCATE(ZBOT(JC,IC),WL(JC,IC))
    ALLOCATE(WL(JC,IC))
    ALLOCATE(XCOR(LCM,5),YCOR(LCM,5),AREA(LCM))
    ALLOCATE(KSZ(LCM))
    ALLOCATE(DZCS(LCM,KCM))
    ALLOCATE(ARR(LCM))

    END SUBROUTINE

    SUBROUTINE INPUT
    INTEGER::L,IS,JCTMP,IACROSS,JACROSS,JT,JF,I,J,JLAST
    INTEGER(4)::ISO,K,IFIRST,ILAST,IT,LT,N
    INTEGER(4)::MASER(100),TCASER,IRELH
    INTEGER(4)::IASWRAD,M,NPNS,MVEGIJT
    REAL(4)   ::TASER,PATM,TDRY,TWET,RAIN,EVAP,SOLSWR,CLOUD,TAASER
    REAL(4)   ::RAINCVT,EVAPCVT,SOLRCVT,CLDCVT,REVC,RCHC,SWRATNF,DABEDT,HTBED1,HTBED2,SWRATNS,FSWRATF
    REAL(4)   ::DXIJ,DYIJ,HIJ,BELVIJ

    CHARACTER ::ADUMMY*5
    CHARACTER*80 TEXT, STR*200

    ! *** READ CELL TYPES FROM FILES CELL.INP
    PRINT *,'*** READING CELL.INP'
    OPEN(1,FILE=CELLFILE,STATUS='UNKNOWN',ERR=995)
    CALL SKIPCOM(1, 'C')
    READ(1,'(A5,640I1)')ADUMMY
    READ(ADUMMY,*)JCTMP
    CLOSE(1)

    OPEN(1,FILE=CELLFILE,STATUS='UNKNOWN')
    CALL SKIPCOM(1, 'C')
    IF(JCTMP.NE.JC)THEN
        ! ***   READ OLD FILE FORMAT
        JACROSS=JC
        IF(JC.GT.640)JACROSS=640
        DO JT=1,JC,JACROSS
            JF=JT
            JLAST=JT+JACROSS-1
            IF(JLAST.GT.JC) JLAST=JC
            DO I=1,IC
                READ(1,'(640I1)',IOSTAT=ISO) (IJCT(I,J),J=JF,JLAST)
                IF(ISO.GT.0) STOP ' **** ERROR IN READING CELL.INP !'
            ENDDO
        ENDDO
    ELSE
        IF(IC.GT.640)THEN
            IACROSS=640
            DO IT=1,IC,IACROSS
                IFIRST=IT
                ILAST=IT+IACROSS-1
                IF(ILAST.GT.IC) ILAST=IC
                DO J=JC,1,-1
                    READ(1,'(A5,640I1)',IOSTAT=ISO)ADUMMY,(IJCT(I,J),I=IFIRST,ILAST)
                    IF(ISO.GT.0) STOP ' **** ERROR IN READING CELL.INP !'
                ENDDO
            ENDDO
        ELSE
            IFIRST=1
            ILAST=IC
            DO J=JC,1,-1
                READ(1,'(A5,640I1)',IOSTAT=ISO)ADUMMY,(IJCT(I,J),I=IFIRST,ILAST)
                IF(ISO.GT.0) STOP ' **** ERROR IN READING CELL.INP !'
            ENDDO
        ENDIF
    ENDIF
    CLOSE(1)

    IF (EFDCVER<730) THEN
        IF (ISTRAN(2)==1 .AND. NASER>0)THEN
            CALL GETASER
        ELSE
            TBEDIT = 0
        ENDIF
    ENDIF

    ! *** IF ISCONNECT GE 1, READ IN NORTH-SOUTH BOUNDARY CELLS FROM
    ! *** FILE MAPPGNS.INP TO SPECIFY A PERIODIC DOMAIN IN THE NORTH-SOUTH
    !IF(ISCONNECT.GE.1)THEN
    !  PRINT *,'***   READING MAPPGNS.INP'
    !  OPEN(1,FILE=MAPGFILE,STATUS='UNKNOWN',ERR=996)
    !  DO IS=1,8
    !    READ(1,*)
    !  ENDDO
    !  READ(1,*,IOSTAT=ISO) NPNSBP
    !  IF(ISO.GT.0) STOP 'READ ERROR FOR FILE MAPPGNS.INP'
    !  DO NPNS=1,NPNSBP
    !    READ(1,*,IOSTAT=ISO) ISPNS(NPNS),JSPNS(NPNS),INPNS(NPNS),JNPNS(NPNS)
    !    IF(ISO.GT.0) STOP 'READ ERROR FOR FILE MAPPGNS.INP'
    !  ENDDO
    !  CLOSE(1)
    !ENDIF

    ! *** IF ISCONNECT GE 2, READ IN EAST-WEST BOUNDARY CELLS FROM
    ! *** FILE MAPPGEW.INP TO SPECIFY A PERIODIC DOMAIN IN THE EAST-WEST DIRECTION
    IF( ISCONNECT >= 2 )THEN
        WRITE(*,'(A)')' *** READING MAPPGEW.INP'
        OPEN(1,FILE=MAPEWFILE,STATUS='UNKNOWN')

        ! *** SKIP OVER TITLE AND AND HEADER LINES
        DO IS=1,8
            READ(1,*)
        ENDDO
        READ(1,*,IOSTAT=ISO) NPEWBP
        IF( ISO > 0 ) STOP ' **** READ ERROR FOR FILE MAPPGEW.INP !'

        ALLOCATE(IWPEW(NPEWBP))
        ALLOCATE(JWPEW(NPEWBP))
        ALLOCATE(IEPEW(NPEWBP))
        ALLOCATE(JEPEW(NPEWBP))

        DO NP=1,NPEWBP
            READ(1,*,IOSTAT=ISO) IWPEW(NP),JWPEW(NP),IEPEW(NP),JEPEW(NP)
            IF( ISO > 0 ) STOP ' **** READ ERROR FOR FILE MAPPGEW.INP !'
        ENDDO
        CLOSE(1)
    ENDIF

    ! *** IF ISCONNECT GE 1, READ IN NORTH-SOUTH BOUNDARY CELLS FROM
    ! *** FILE MAPPGNS.INP TO SPECIFY A PERIODIC DOMAIN IN THE NORTH-SOUTH DIRECTION
    IF( ISCONNECT == 1 .OR. ISCONNECT == 3 )THEN
        WRITE(*,'(A)')' *** READING MAPPGNS.INP'
        OPEN(1,FILE=MAPNSFILE,STATUS='UNKNOWN')

        ! *** SKIP OVER TITLE AND AND HEADER LINES
        DO IS=1,8
            READ(1,*)
        ENDDO
        READ(1,*,IOSTAT=ISO) NPNSBP
        IF( ISO > 0 ) STOP ' **** READ ERROR FOR FILE MAPPGNS.INP !'

        ALLOCATE(ISPNS(NPNSBP))
        ALLOCATE(JSPNS(NPNSBP))
        ALLOCATE(INPNS(NPNSBP))
        ALLOCATE(JNPNS(NPNSBP))

        DO NP=1,NPNSBP
            READ(1,*,IOSTAT=ISO) ISPNS(NP),JSPNS(NP),INPNS(NP),JNPNS(NP)
            IF( ISO > 0 ) STOP ' **** READ ERROR FOR FILE MAPPGNS.INP !'
        ENDDO
        CLOSE(1)
    ENDIF

    CALL CELLMAP

    !IMN = MINVAL(IL(1:LA),IL>0)
    !JMN = MINVAL(JL(1:LA),JL>0)
    !IMX = MAXVAL(IL(1:LA))
    !JMX = MAXVAL(JL(1:LA))

22  FORMAT (A80)

    IF(LVC > 0)THEN
        PRINT *,'*** READING DXDY.INP'
        OPEN(1,FILE=DXDYFILE,STATUS='UNKNOWN',ERR=997)
        CALL SKIPCOM(1, 'C')
        DX =0
        DY =0

        DO LT=1,LVC
            READ(1,*,IOSTAT=ISO)I,J,DXIJ,DYIJ,HIJ,BELVIJ
            IF(ISO.GT.0) STOP ' **** READ ERROR FOR FILE DXDY.INP !'
            L=LIJ(I,J)
            BELV(L)=BELADJ + BELCVRT*BELVIJ
            HP0(L) = HIJ
            DX = DX+DXIJ
            DY = DY+DYIJ
        ENDDO
        DX = DX/LVC
        DY = DY/LVC
    ENDIF

    ISGWIE=0
    IF(ISGWIT==1)THEN
        OPEN(1,FILE='gwater.inp',STATUS='UNKNOWN',ERR=998)
        CALL SKIPCOM(1,'*')
        READ(1,*) ISGWIE
        CLOSE(1)
    ENDIF

    IF(ISPD >= 2) CALL DRIFTERINP

    RETURN
995 PRINT*,' **** OPEN ERROR: CELL.INP !'
    PAUSE(1)
    STOP
996 PRINT*,' **** OPEN ERROR: MAPPGNS.INP !'
    PAUSE(1)
    STOP
997 PRINT*,' **** OPEN ERROR: DXDY.INP !'
    PAUSE(1)
    STOP
998 PRINT*,' **** OPEN ERROR: GWATER.INP !'
    PAUSE(1)
    STOP
    END SUBROUTINE

    SUBROUTINE CELLMAP
    INTEGER :: L,I,J,LCTT,NPN,LE,LW,LS,LN,NW,IGTMP,JGTMP,ICOMP,JCOMP,LTMP
    CHARACTER(200) :: STR

    ! **  SET 1D CELL INDEX SEQUENCE AND MAPPINGS
    L=2
    IL=0
    JL=0

    DO J=1,JC
        DO I=1,IC
            IF( IJCT(I,J) > 0 .AND. IJCT(I,J) < 9 )THEN
                IL(L)=I
                JL(L)=J
                LCT(L)=IJCT(I,J)
                LIJ(I,J)=L
                L=L+1
            ENDIF
        ENDDO
    ENDDO

    LA=L-1

    ! **  SET NORTH AND SOUTH CELL IDENTIFIER ARRAYS
    LNC(1)=LC
    LSC(1)=LC
    LNEC(1)=LC
    LNWC(1)=LC
    LSEC(1)=LC
    LSWC(1)=LC
    LNC(LC)=1
    LSC(LC)=1
    LNEC(LC)=1
    LNWC(LC)=1
    LSEC(LC)=1
    LSWC(LC)=1
    DO L=2,LA
        I=IL(L)
        J=JL(L)
        IF( IJCT(I,J+1) == 9 )THEN
            LNC(L)=LC
        ELSE
            LNC(L)=LIJ(I,J+1)
        ENDIF
        IF( IJCT(I,J-1) == 9 )THEN
            LSC(L)=LC
        ELSE
            LSC(L)=LIJ(I,J-1)
        ENDIF
        IF( IJCT(I+1,J+1) == 9 )THEN
            LNEC(L)=LC
        ELSE
            LNEC(L)=LIJ(I+1,J+1)
        ENDIF
        IF( IJCT(I-1,J+1) == 9 )THEN
            LNWC(L)=LC
        ELSE
            LNWC(L)=LIJ(I-1,J+1)
        ENDIF
        IF( IJCT(I+1,J-1) == 9 )THEN
            LSEC(L)=LC
        ELSE
            LSEC(L)=LIJ(I+1,J-1)
        ENDIF
        IF( IJCT(I-1,J-1) == 9 )THEN
            LSWC(L)=LC
        ELSE
            LSWC(L)=LIJ(I-1,J-1)
        ENDIF
    ENDDO

    ! *** CONNECT EAST-WEST CELLS
    LEC(1)=LC
    LWC(1)=LC
    LEC(LC)=1
    LWC(LC)=1
    DO L=2,LA
        I=IL(L)
        J=JL(L)
        IF( IJCT(I+1,J) == 9 )THEN
            LEC(L)=LC
        ELSE
            LEC(L)=LIJ(I+1,J)
        ENDIF
        IF( IJCT(I-1,J) == 9 )THEN
            LWC(L)=LC
        ELSE
            LWC(L)=LIJ(I-1,J)
        ENDIF
    ENDDO

    ! **  MODIFY EAST-WEST CELL MAPPING FOR PERIOD GRID IN E-W DIRECTION
    IF( ISCONNECT >= 2 )THEN
        DO NPN=1,NPEWBP
            L=LIJ(IWPEW(NPN),JWPEW(NPN))
            LWC(L)=LIJ(IEPEW(NPN),JEPEW(NPN))
            IF( IJCT(IEPEW(NPN),JEPEW(NPN)-1) == 9 )THEN
                LSWC(L)=LC
            ELSE
                LSWC(L)=LIJ(IEPEW(NPN),JEPEW(NPN)-1)
            ENDIF
            IF( IJCT(IEPEW(NPN),JEPEW(NPN)+1) == 9 )THEN
                LNWC(L)=LC
            ELSE
                LNWC(L)=LIJ(IEPEW(NPN),JEPEW(NPN)+1)
            ENDIF

            L=LIJ(IEPEW(NPN),JEPEW(NPN))
            LEC(L)=LIJ(IWPEW(NPN),JWPEW(NPN))
            IF( IJCT(IWPEW(NPN),JWPEW(NPN)-1) == 9 )THEN
                LSEC(L)=LC
            ELSE
                LSEC(L)=LIJ(IWPEW(NPN),JWPEW(NPN)-1)
            ENDIF
            IF( IJCT(IWPEW(NPN),JWPEW(NPN)+1) == 9 )THEN
                LNEC(L)=LC
            ELSE
                LNEC(L)=LIJ(IWPEW(NPN),JWPEW(NPN)+1)
            ENDIF
        ENDDO
    ENDIF

    ! **  MODIFY NORTH-SOUTH CELL MAPPING FOR PERIOD GRID IN N-S DIRECTION
    IF( ISCONNECT == 1 .OR. ISCONNECT == 3 )THEN
        DO NPN=1,NPNSBP
            LS=LIJ(ISPNS(NPN),JSPNS(NPN))
            LSC(LS)=LIJ(INPNS(NPN),JNPNS(NPN))
            IF( IJCT(INPNS(NPN)+1,JNPNS(NPN)) == 9 )THEN
                LSEC(LS)=LC
            ELSE
                LSEC(LS)=LIJ(INPNS(NPN)+1,JNPNS(NPN))
            ENDIF
            IF( IJCT(INPNS(NPN)-1,JNPNS(NPN)) == 9 )THEN
                LSWC(LS)=LC
            ELSE
                LSWC(LS)=LIJ(INPNS(NPN)-1,JNPNS(NPN))
            ENDIF
            LN=LIJ(INPNS(NPN),JNPNS(NPN))
            LNC(LN)=LIJ(ISPNS(NPN),JSPNS(NPN))
            IF( IJCT(ISPNS(NPN)+1,JSPNS(NPN)) == 9 )THEN
                LNEC(LN)=LC
            ELSE
                LNEC(LN)=LIJ(ISPNS(NPN)+1,JSPNS(NPN))
            ENDIF
            IF( IJCT(ISPNS(NPN)-1,JSPNS(NPN)) == 9 )THEN
                LNWC(LN)=LC
            ELSE
                LNWC(LN)=LIJ(ISPNS(NPN)-1,JSPNS(NPN))
            ENDIF
        ENDDO
    ENDIF

    END SUBROUTINE

    SUBROUTINE SEEK(TAG,STAT)
    INTEGER(4),INTENT(OUT),OPTIONAL::STAT
    CHARACTER TAG*(*)
    CHARACTER*80 TEXT
    INTEGER(4)::I,J,K,L,M

    IF (PRESENT(STAT)) STAT=1
    L=LEN(TAG)
    DO I=1,L
        J=ICHAR(TAG(I:I))
        IF(97.LE.J.AND.J.LE.122)THEN
            TAG(I:I)=CHAR(J-32)
        ENDIF
    ENDDO

    DO K=1,2
10      READ(1,'(A)',END=20)TEXT
        M=MAX(1,LEN_TRIM(TEXT))
        DO WHILE(M.GT.L.AND.TEXT(1:1).EQ.'')
            TEXT(1:M-1)=TEXT(2:M)
            TEXT(M:M)=' '
            M=M-1
        ENDDO
        IF(M.LT.L)GO TO 10
        DO I=1,M
            J=ICHAR(TEXT(I:I))
            IF(97.LE.J.AND.J.LE.122)THEN
                TEXT(I:I)=CHAR(J-32)
            ENDIF
        ENDDO
        IF(TEXT(1:L).NE.TAG)GO TO 10
        IF(TEXT(L+1:L+1).NE.' ')GO TO 10
    ENDDO
    RETURN
20  WRITE(*,'(A,A,A)')'GROUP: ',TAG,' NOT FOUND BEFORE END OF FILE'
    IF (PRESENT(STAT)) STAT=0
    END SUBROUTINE

    ! *** NDLuan Added-2018-12-04
    FUNCTION SEEK_SCANLINE(TAG1, TAG2)  RESULT(FOUND)
    CHARACTER TAG1*(*)
    CHARACTER TAG2*(*)
    CHARACTER*80 TEXT
    INTEGER(4)::I,J,K,L1,L2,M,FOUND

    L1=LEN(TAG1)
    DO I=1,L1
        J=ICHAR(TAG1(I:I))
        IF(97.LE.J.AND.J.LE.122)THEN
            TAG1(I:I)=CHAR(J-32)
        ENDIF
    ENDDO

    DO K=1,2
10      READ(1,'(A)',END=20)TEXT
        M=MAX(1,LEN_TRIM(TEXT))
        DO WHILE(M.GT.L1.AND.TEXT(1:1).EQ.'')
            TEXT(1:M-1)=TEXT(2:M)
            TEXT(M:M)=' '
            M=M-1
        ENDDO
        IF(M.LT.L1)GO TO 10
        DO I=1,M
            J=ICHAR(TEXT(I:I))
            IF(97.LE.J.AND.J.LE.122)THEN
                TEXT(I:I)=CHAR(J-32)
            ENDIF
        ENDDO
        IF(TEXT(1:L1).NE.TAG1)GO TO 10
        IF(TEXT(L1+1:L1+1).NE.' ')GO TO 10
    ENDDO
    FOUND=0
    L2=LEN(TAG2)

    DO I=1,L2
        J=ICHAR(TAG2(I:I))
        IF(97.LE.J.AND.J.LE.122)THEN
            TAG2(I:I)=CHAR(J-32)
        ENDIF
    ENDDO
    M=MAX(1,LEN_TRIM(TEXT))
    I = SCAN(TEXT,TAG2)
    IF (I.GT.0.OR.I.LE.M) THEN
        IF(TEXT(I:L2+I-1).EQ.TAG2) THEN
            FOUND=1
            RETURN
        ENDIF
    ENDIF

    RETURN
20  WRITE(*,'(A,A,A)')'GROUP: ',TAG1,' NOT FOUND BEFORE END OF FILE'
    END FUNCTION

    SUBROUTINE GETASER
    INTEGER(4):: N,POS,LS,M
    INTEGER(4):: MASER,IRELH,IASERVER
    REAL(4)   :: VAL(10),TCASER,TAASER,RAINCVT,EVAPCVT,SOLRCVT,CLDCVT
    REAL(4)   :: DS_LONG,DS_LAT,WINDFA,WINDFB,WINDFC
    LOGICAL(4):: COMPUTESOLRAD,USESHADE
    CHARACTER(200)::STR

    OPEN(UNAS,FILE=ASERFILE,ACTION='READ',ERR=998)
    TBEDIT = 0

    ! ** ONLY FOR EFDC_70-71-72:
    CALL SKIPCOM(UNAS,'C')
    DO N=1,4
        READ(UNAS,'(A)') STR  ! READING EE LINES
    ENDDO
    CALL SKIPCOM(UNAS,'C')

    READ(UNAS,*,END=100,ERR=999) MASER,TCASER,TAASER,IRELH,RAINCVT,EVAPCVT,SOLRCVT,CLDCVT
    READ(UNAS,*,END=100,ERR=999) VAL

    TBEDIT = VAL(8)

100 CLOSE(UNAS)
    RETURN
998 PRINT*,' **** OPEN ERROR: ASERFILE !'
    PAUSE(1)
    STOP
999 STOP ' **** READING ASER.INP ERROR !'
    PAUSE(1)
    STOP
    END SUBROUTINE

    SUBROUTINE FREEMEM
    DEALLOCATE(IJCT )
    DEALLOCATE(IL )
    DEALLOCATE(JL )
    DEALLOCATE(LNC )
    DEALLOCATE(LSC )
    DEALLOCATE(LIJ )

    DEALLOCATE(CUE )
    DEALLOCATE(CUN )
    DEALLOCATE(CVE )
    DEALLOCATE(CVN )
    DEALLOCATE(BELV )
    IF (ALLOCATED(HP))  DEALLOCATE(HP)
    IF (ALLOCATED(HP0)) DEALLOCATE(HP0)

    DEALLOCATE(RSSBCE )
    DEALLOCATE(RSSBCN )
    DEALLOCATE(RSSBCS )
    DEALLOCATE(RSSBCW )

    DEALLOCATE(ISPNS )
    DEALLOCATE(JSPNS )
    DEALLOCATE(INPNS )
    DEALLOCATE(JNPNS )

    DEALLOCATE(U)
    DEALLOCATE(V)
    DEALLOCATE(W)
    DEALLOCATE(UK)
    DEALLOCATE(VK)

    DEALLOCATE(XCC)
    DEALLOCATE(YCC)

    IF (ALLOCATED(HSIG)) DEALLOCATE(HSIG,DIR,RTP,RLS,EDIS)
    IF (ALLOCATED(SXX))  DEALLOCATE(SXX,SXY,SYY,WVENE)
    IF (ALLOCATED(XP))   DEALLOCATE(XP,YP,SXXP,SXYP,SYYP)
    IF (ALLOCATED(HSP))  DEALLOCATE(HSP,DIRP,RTPP,RLSP,EDISP)

    END SUBROUTINE

    SUBROUTINE DRIFTERINP
    ! ********************************************************************
    !READING INPUT DATA OF INITIAL LOCATIONS OF DRIFTERS
    !OUTPUT: NPD,XLA,YLA,ZLA,NP=1:NPD
    !        LA_BEGTI, LA_ENDTI, LA_FREQ,LANDT
    INTEGER    :: NP,N,K,NPN,IN,JN,IS,JS,ITMP
    REAL(RKD)  :: RANVAL,RTMP
    CHARACTER(200)  :: STR

    OPEN(ULOC,FILE=DRIFILE,ACTION='READ')
    CALL SKIPCOM(ULOC,'*')
    READ(ULOC,*) ! LA_ZCAL,LA_PRAN,LA_DIFOP,LA_HORDIF,LA_VERDIF,DEPOP,ADJVERP    ! AUG 2016: NEW STRUCTURE
    CALL SKIPCOM(ULOC,'*')
    READ(ULOC,*) LA_BEGTI0,LA_ENDTI0,LA_FREQ
    CALL SKIPCOM(ULOC,'*')
    READ(ULOC,*) NPD,NGRP
    !LA_FREQ = LA_FREQ/1440.                             !Output Frequency

    ! *** ALLOCATE THE DRIFTER ARRAYS
    ALLOCATE(XLA(NPD),YLA(NPD),ZLA(NPD),LLA(NPD))
    ALLOCATE(LA_BEGTI(NPD),LA_ENDTI(NPD),LA_GRP(NPD))
    ALLOCATE(ISOILSPI(NGRP),DDEN(NGRP),MOC(LA))
    ALLOCATE(DVOL(NPD))

    ! *** INIITALIZE THE ARRAYS (NPD)
    XLA=0.0
    YLA=0.0
    ZLA=0.0
    LLA = 0
    ZLA = 0
    MOC = 0

    CALL SKIPCOM(ULOC,'*')
    DO N=1,NGRP
        READ(ULOC,'(A)') STR
        READ(STR,*) K,ISOILSPI(K) !,GRPWS(K),BEDFIX(K)
        IF( ISOILSPI(K) > 0 )THEN
            READ(STR,*) K,ISOILSPI(K),RTMP,ITMP,RTMP,DDEN(K) !,DRAT(K),DTEM(K),DVAP(K),DVMO(K)
            !DRAT(K)=DRAT(K)/86400.
        ENDIF
    ENDDO

    CALL SKIPCOM(ULOC,'*')
    DO NP=1,NPD
        ! *** Read Elevations
        READ(ULOC,*,ERR=999) XLA(NP),YLA(NP),ZLA(NP),LA_BEGTI(NP),LA_ENDTI(NP),LA_GRP(NP)
        IF( LA_GRP(NP) < 1)LA_GRP(NP)=1
    ENDDO
    CLOSE(ULOC)

    LPTSNP = NINT((LA_ENDTI0-LA_BEGTI0)*1440._8/LA_FREQ) + 1

    RETURN
999 STOP 'DRIFTER.INP READING ERROR!'
    END SUBROUTINE

    END MODULE

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus