text:
      module  m_GFIO_GetFld

! !USES:

      implicit none

      PRIVATE
      PUBLIC GFIO_GetFld

      interface GFIO_GetFld; module procedure
     &          GFIO_Get_Fall_,  ! GEOS-3
     &          GFIO_Get_2d_,    ! GEOS-3
     &          GFIO_Get_3d_     ! GEOS-3
      end interface
      interface Fix_Undef_In; module procedure
     &          Fix_Undef_In_
      end interface

!  Internal Interfaces
!  -------------------
      interface Hflip_
        module procedure HFlip2_    ! Horizontal flip of 2d fields
        module procedure HFlip3_    ! Horizontal flip of 3d fields
      end interface

      character(len=*), parameter :: myname = 'GFIO_GetFld'

      logical, parameter :: verb   = .true.
      integer, parameter :: stdout = 6
      integer, parameter :: stderr = 6
      real   , parameter :: UNDEF  = 1.e15
      real   , parameter :: rtol   = 1.e-3  ! 25Sep2007 Todling: Adjusted tol

      CONTAINS

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE:  GFIO_Get_Fall_ ---
!
! !INTERFACE:
!
      subroutine GFIO_Get_Fall_ ( gfiofile, nymd, nhms, im, jm, mlev, plevs,
     &                            nvars2d, vname2d, fld2d,
     &                            nvars3d, vname3d, fld3d,
     &                            stat )

! !USES:

      Implicit NONE

!
! !INPUT PARAMETERS:
!

      character*(*) gfiofile       ! Input GFIO file name
                                   ! Synoptic time:
      integer       nymd           !  Year-month-day, e.g.,  19971012
      integer       nhms           !  hour-minute-sec, e.g.,   120000
      integer       im
      integer       jm
      integer       mlev
      integer       nvars2d
      integer       nvars3d
      character*(*) vname2d(nvars2d)
      character*(*) vname3d(nvars3d)
      real          plevs(mlev)

! !OUTPUT PARAMETERS:

      real          fld2d(im,jm,nvars2d)
      real          fld3d(im,jm,mlev,nvars3d)

      integer, intent(out), optional :: stat

!
! !DESCRIPTION: Read grided analysis fields using GFIO pre- and post-
!               analysis files.
!
! !TO DO:
!   replace all exits by output return code
!
! !REMARKS:
!
! !REVISION HISTORY:
!
!  21Nov2000 Todling   Initial code.
!  24Feb2006 Todling   Adapted to handle flipped fields from GEOS-5.
!  11Jul2007 Stassi/RT Fix for lon precision issue between G5 and G4 files.
!  16Jul2008 Todling   Lon check now using 64-bin array
!
!EOP
!-------------------------------------------------------------------------
!_BOC

      character(len=*), parameter :: myname_ = myname//'*'//'GFIO_Get_Fall_'

                                   ! Grid Information
                                   !    direction
      integer       ier, year, month, day, hour
      integer       nobs_oi, i, j, k, kk, ia, ios

      integer       lm_max_max          ! Maximum number of time stamps in file

      integer       ngatts              ! number of attributes on file
      integer       nvars, nvars_in
      integer       im_max, jm_max, mlev_max, lm_max
      integer       n2dim, n3dim
      integer       iv, icount
      integer       ncid
      integer       dimSize, dimId
      integer       ncvid
      character*255 dimName
      integer       kloc(1), levcnt
      logical       usrlevs
      logical       fliplon

      real*8  mxlon, mnlon
      real*4, allocatable :: lon32(:)
      real*8, allocatable :: lon64(:)
      real, allocatable :: lev(:)
      real, allocatable :: aux(:,:,:,:)

!     Define constants
!     ----------------
      nvars = nvars2d + nvars3d
      usrlevs = .true.
      if ( present(stat) ) stat = 0

!     Open GFIO file
!     --------------
      im_max   = im
      jm_max   = jm
      mlev_max = mlev
      call GFIO_Open ( trim(gfiofile), 0, ncid, ier )

      if ( ier .eq. 0 ) then

         call GFIO_DimInquire ( ncid, im_max, jm_max, mlev_max, lm_max,
     &                          nvars_in, ngatts, ier)

         if ( nvars .gt. nvars_in ) then
             write(stdout,'(a,2(a,i5))') myname_,
     &                    ': Warning, no. of variable in file: ', nvars_in,
     &                    ' smaller than requested no. of vars: ', nvars
         end if

!        Check consistency of Gfio file
!        ------------------------------
         if (im_max .ne. im) then
            write(stderr,'(2a,i5)') myname_,
     &           ': trying to read unallowed longitudinal dim'
            call exit(7)
         endif

         if (jm_max .ne. jm) then
            write(stderr,'(2a,i5)') myname_,
     &           ': trying to read unallowed latitudinal dim'
            call exit(7)
         endif

         if (mlev_max .ne. mlev) then
            write(stdout,'(2a,i5)') myname_,
     &           '-WARNING: trying to read unallowed number of levels'
            usrlevs = .false.
         endif

      else

         write(stdout,'(2a)') myname_,
     .                ': No GFIO file for current date '
         call flush (stdout)
         if ( present(stat) ) then
           stat = 99
           return
         else
           call exit(7)
         end if

      endif

!     OK, let's read data from file
!     =============================
      icount = 0
      n2dim  = im * jm
      n3dim  = n2dim * mlev

!     Get longitude array to check on GEOS-5 flipped fields
!     -----------------------------------------------------
      allocate ( lon64(im), stat=ier )
        if (ier/=0) then
            write(stderr,'(2a,i10,2a)') myname_, 'Alloc(lon) error '
            call flush (stderr); call exit(7)
        endif

      call ncdinq (ncid, 1, dimName, dimSize, ier)
         if (ier/=0)then
             write(stderr,'(2a)') myname_, 'Error(ncdinq)'
             call flush (stderr); call exit(7)
         endif
      dimId = ncvid (ncid, dimName, ier)
         if (ier/=0)then
             write(stderr,'(2a)') myname_, 'Error(ncvid)'
             call flush (stderr); call exit(7)
         endif

!     Caution: will try to read as 64-bit first ...
!     ---------------------------------------------
      call ncvgt ( ncid, dimId, 1, im, lon64, ier )
           if ( ier .ne. 0 ) then
             write(stderr,'(2a,i10,2a)') myname_,
     &                    ': Gfio_GetVar, ier = ', ier,
     &                    ' for variable ', 'lon'
             call flush(stderr); call exit(7)
            endif

!     ... if absurd lons, read as 32-bit
!     ----------------------------------
      mnlon = minval(lon64)
      mxlon = maxval(lon64)
      if ( mnlon < -200.0 .or. mxlon > 400.0 ) then
         allocate ( lon32(im), stat=ier )
         call ncvgt ( ncid, dimId, 1, im, lon32, ier )
         lon64=lon32
         deallocate ( lon32 )
      endif

!     Is orientation of fields shifted by 180 degs from internal def?
!     ---------------------------------------------------------------
      fliplon = .false.
      if ( abs(lon64(1)+180.0)<rtol ) then
        fliplon = .true.  ! orientation of fields shifted by 180 degs from internal def
        print *, ' Flipping Fields Longitudinally ... '
      endif
      deallocate ( lon64 )

!     First read surface fields
!     -------------------------
      do iv = 1, nvars2d
      call Gfio_GetVar ( ncid, vname2d(iv), nymd, nhms,
     &                   im, jm, 0, 1, fld2d(1,1,iv), ier )
           if ( ier .ne. 0 ) then
             write(stderr,'(2a,i10,2a)') myname_,
     &                    ': Gfio_GetVar, ier = ', ier,
     &                    ' for variable ', trim(vname2d(iv))
             call exit(7)
           else
             call fix_undef_in_ ( n2dim, fld2d(:,:,iv:iv), vname2d(iv) )
           end if
           if( fliplon ) call hflip_ ( fld2d(1:im,1:jm,iv), im, jm )
           icount = icount + 1
      end do

!     Now read upper-air fields
!     -------------------------
      if ( usrlevs ) then
        do iv = 1, nvars3d
           call Gfio_GetVar ( ncid, vname3d(iv), nymd, nhms,
     &                        im, jm, 1, mlev, fld3d(1,1,1,iv), ier )
           if ( ier .ne. 0 ) then
             write(stderr,'(2a,i10)') myname_,
     &                    ': Gfio_GetVar, ier = ', ier,
     &                    ' for variable ', trim(vname3d(iv))
             call exit(7)
           else
             call fix_undef_in_ ( n3dim, fld3d(:,:,:,iv:iv), vname3d(iv) )
           end if
           if( fliplon ) call hflip_ ( fld3d(1:im,1:jm,1:mlev,iv), im, jm, mlev )
           icount = icount + 1
        end do

      else

        allocate ( aux(im,jm,mlev_max,1), lev(mlev_max), stat=ier )
             if ( ier .ne. 0 ) then
                  write(stderr,'(2a,i10)') myname_,
     .                         ': Alloc(aux), ier = ', ier
                  call exit(7)
             end if
        call Gfio_GetVar ( ncid, 'Height', nymd, nhms,
     &                     mlev_max, 1, 0, 1, lev, ier )
           if ( ier .ne. 0 ) then
             write(stderr,'(2a,i10)') myname_,
     &                    ': Gfio_GetVar (lev), ier = ', ier
             call exit(7)
           end if

        do iv = 1, nvars3d
           call Gfio_GetVar ( ncid, vname3d(iv), nymd, nhms,
     &                        im, jm, 1, mlev_max, aux, ier )
           if ( ier .ne. 0 ) then
             write(stderr,'(2a,i10,2a)') myname_,
     &                    ': Gfio_GetVar, ier = ', ier,
     &                    ' for variable ', trim(vname3d(iv))
             call exit(7)
           else
             call fix_undef_in_ ( im*jm*mlev_max, aux, vname3d(iv) )
           end if
           if( fliplon ) call hflip_ ( fld3d(1:im,1:jm,1:mlev_max,iv), im, jm, mlev_max )

           levcnt = 0
           fld3d(:,:,:,iv:iv) = undef
           do k = 1, mlev
              kloc = maxloc(lev(1:mlev_max),
     .               mask=abs(log(plevs(k))-log(lev(1:mlev_max))).lt.1e-5)
              j    = kloc(1)
              if (j>0) then
                 fld3d(:,:,k:k,iv:iv) = aux(:,:,j:j,1:1)
                 levcnt = levcnt + 1
              end if
           end do
           if ( levcnt.ne.mlev ) then
             write(stdout,'(a,2(i10))') myname_,
     &                    'WARNING: Requested levs not found, levcnt = ',
     &                     levcnt, ' levs requested ', mlev
           end if
           icount = icount + 1
        end do

        deallocate ( aux, lev, stat=ier )

      end if

!     Close GFIO file now
!     -------------------
      call Gfio_Close ( ncid, ier )
      if ( ier .ne. 0 ) then
           write(stdout,'(2a,i7)') myname_,
     &                  ': WARNING, problems closing GFIO file, ier = ', ier
      end if

      write(stdout,'(2a,i12,2a)') myname_,
     &            ':  read ', icount, ' variables from ', trim(gfiofile)

!     All done
!     --------
      return
      end subroutine GFIO_Get_Fall_
!EOC

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE:  GFIO_Get_2d_ ---
!
! !INTERFACE:
!
      subroutine GFIO_Get_2d_ ( gfiofile, nymd, nhms, im, jm,
     &                          nvars, vname, fld )

! !USES:

      Implicit NONE

!
! !INPUT PARAMETERS:
!

      character*(*) gfiofile       ! Input GFIO file name
                                   ! Synoptic time:
      integer       nymd           !  Year-month-day, e.g.,  19971012
      integer       nhms           !  hour-minute-sec, e.g.,   120000
      integer       im
      integer       jm
      integer       nvars
      character*(*) vname(nvars)

! !OUTPUT PARAMETERS:

      real          fld(im,jm,nvars)

!
! !DESCRIPTION: Read grided analysis fields using GFIO pre- and post-
!               analysis files.
!
! !TO DO:
!
! !REMARKS:
!
! !REVISION HISTORY:
!
!  21Nov2000 Todling   Initial code.
!  24Feb2006 Todling   Adapted to handle flipped fields from GEOS-5.
!  11Jul2007 Stassi/RT Fix for lon precision issue between G5 and G4 files.
!  16Jul2008 Todling   Lon check now using 64-bin array
!
!EOP
!-------------------------------------------------------------------------
!_BOC

      character(len=*), parameter :: myname_ = myname//'*'//'GFIO_Get_2d_'

                                   ! Grid Information
                                   !    direction
      integer       ier, year, month, day, hour
      integer       nobs_oi, i, j, k, kk, ia, ios

      integer       lm_max_max          ! Maximum number of time stamps in file

      integer       ngatts              ! number of attributes on file
      integer       nvars_in
      integer       im_max, jm_max, mlev_max, lev_max, lm_max
      integer       n2dim
      integer       iv, icount
      integer       ncid
      integer       ncvid
      integer       dimId, dimSize
      character*255 dimName

      real*8  mxlon, mnlon
      real*4, allocatable :: lon32(:)
      real*8, allocatable :: lon64(:)

      logical       fliplon

!     Open GFIO file
!     --------------
      im_max   = im
      jm_max   = jm
      call GFIO_Open ( trim(gfiofile), 0, ncid, ier )

      if ( ier .eq. 0 ) then

         call GFIO_DimInquire ( ncid, im_max, jm_max, mlev_max, lm_max,
     &                          nvars_in, ngatts, ier)

         if ( nvars .gt. nvars_in ) then
             write(stdout,'(a,2(a,i5))') myname_,
     &                    ': Warning, no. of variable in file: ', nvars_in,
     &                    ' smaller than requested no. of vars: ', nvars
         end if

!        Check consistency of Gfio file
!        ------------------------------
         if (im_max .ne. im) then
            write(stderr,'(2a,i5)') myname_,
     &           ': trying to read unallowed longitudinal dim'
            call exit(7)
         endif

         if (jm_max .ne. jm) then
            write(stderr,'(2a,i5)') myname_,
     &           ': trying to read unallowed latitudinal dim'
            call exit(7)
         endif

      else

         write(stdout,'(2a)') myname_,
     .                ': No GFIO file for current date '
         call flush (stdout)
         call exit(7)

      endif

!     OK, let's read data from file
!     =============================
      icount = 0
      n2dim  = im * jm

!     Get longitude array to check on GEOS-5 flipped fields
!     -----------------------------------------------------
      allocate ( lon64(im), stat=ier )
        if (ier/=0) then
            write(stderr,'(2a,i10,2a)') myname_, 'Alloc(lon) error '
            call flush (stdout); call exit(7)
        endif

      call ncdinq (ncid, 1, dimName, dimSize, ier)
         if (ier/=0)then
             write(stderr,'(2a)') myname_, 'Error(ncdinq)'
             call flush (stderr); call exit(7)
         endif
      dimId = ncvid (ncid, dimName, ier)
         if (ier/=0)then
             write(stderr,'(2a)') myname_, 'Error(ncvid)'
             call flush (stderr); call exit(7)
         endif

!     Caution: will try to read as 64-bit first ...
!     ---------------------------------------------
      call ncvgt ( ncid, dimId, 1, im, lon64, ier )
           if ( ier .ne. 0 ) then
             write(stderr,'(2a,i10,2a)') myname_,
     &                    ': Gfio_GetVar, ier = ', ier,
     &                    ' for variable ', 'lon'
             call flush(stderr); call exit(7)
            endif

!     ... if absurd lons, read as 32-bit
!     ----------------------------------
      mnlon = minval(lon64)
      mxlon = maxval(lon64)
      if ( mnlon < -200.0 .or. mxlon > 400.0 ) then
         allocate ( lon32(im), stat=ier )
         call ncvgt ( ncid, dimId, 1, im, lon32, ier )
         lon64=lon32
         deallocate ( lon32 )
      endif

!     Is orientation of fields shifted by 180 degs from internal def?
!     ---------------------------------------------------------------
      fliplon = .false.
      if ( abs(lon64(1)+180.0)<rtol ) then
        fliplon = .true.  ! orientation of fields shifted by 180 degs from internal def
        print *, ' Flipping Fields Longitudinally ... '
      endif
      deallocate ( lon64 )

!     First read surface fields
!     -------------------------
      do iv = 1, nvars
      call Gfio_GetVar ( ncid, vname(iv), nymd, nhms,
     &                   im, jm, 0, 1, fld(1,1,iv), ier )
           if ( ier .ne. 0 ) then
             write(stderr,'(2a,i10)') myname_,
     &                    ': Gfio_GetVar, ier = ', ier,
     &                    ' for variable ', trim(vname(iv))
             call exit(7)
           else
             call fix_undef_in_ ( n2dim, fld(:,:,iv:iv), vname(iv) )
           end if
           if( fliplon ) call hflip_ ( fld(1:im,1:jm,iv), im, jm )
           icount = icount + 1
      end do

!     Close GFIO file now
!     -------------------
      call Gfio_Close ( ncid, ier )
      if ( ier .ne. 0 ) then
           write(stdout,'(2a,i7)') myname_,
     &                  ': WARNING, problems closing GFIO file, ier = ', ier
      end if

      write(stdout,'(2a,i12,2a)') myname_,
     &            ':  read ', icount, ' variables from ', trim(gfiofile)

!     All done
!     --------
      return
      end subroutine GFIO_Get_2d_
!_EOC

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE:  GFIO_Get_3d_ ---
!
! !INTERFACE:
!
      subroutine GFIO_Get_3d_ ( gfiofile, nymd, nhms, im, jm, mlev,
     &                          nvars, vname, fld,
     &                          plevs, verbose, stat )   ! optional

! !USES:

      Implicit NONE

!
! !INPUT PARAMETERS:
!

      character*(*) gfiofile       ! Input GFIO file name
                                   ! Synoptic time:
      integer       nymd           !  Year-month-day, e.g.,  19971012
      integer       nhms           !  hour-minute-sec, e.g.,   120000
      integer       im
      integer       jm
      integer       mlev
      integer       nvars
      character*(*) vname(nvars)
      real, intent(in), optional ::  plevs(mlev)
      logical, intent(in), optional :: verbose

! !OUTPUT PARAMETERS:

      real          fld(im,jm,mlev,nvars)

! !INPUT/OUTPUT PARAMETERS:

      integer, intent(out), optional :: stat   ! stat=1 incompatible dims
                                               ! stat=2 date/time not found
                                               ! stat=3 var not found
                                               ! stat=4 cannot alloc mem
!
! !DESCRIPTION: Read grided analysis fields using GFIO pre- and post-
!               analysis files.
!
! !TO DO:
!
! !REMARKS:
!
! !REVISION HISTORY:
!
!  21Nov2000 Todling   Initial code.
!  05Oct2004 Todling   Generalized interface in an attempt to remove redundant
!                      interface added by adjoint testing group.
!  24Feb2006 Todling   Adapted to handle flipped fields from GEOS-5.
!  01Feb2007 Todling   Add check for ROOT pe
!  11Jul2007 Stassi/RT Fix for lon precision issue between G5 and G4 files.
!  27Sep2007 Todling   Replace ROOT check by verbose
!  16Jul2008 Todling   Lon check now using 64-bin array
!
!EOP
!-------------------------------------------------------------------------
!_BOC

      character(len=*), parameter :: myname_ = myname//'*'//'GFIO_Get_3d_'

                                   ! Grid Information
                                   !    direction
      integer       ier, year, month, day, hour
      integer       nobs_oi, i, j, k, kk, ia, ios

      integer       lm_max_max          ! Maximum number of time stamps in file

      integer       ngatts              ! number of attributes on file
      integer       nvars_in
      integer       im_max, jm_max, mlev_max, lm_max
      integer       n2dim, n3dim
      integer       iv, icount
      integer       ncid
      integer       dimId, dimSize
      integer       ncvid
      integer       kloc(1), levcnt
      character*255 dimName
      logical       usrlevs
      logical       fliplon
      logical       verbose_

      real*8  mxlon, mnlon
      real*4, allocatable :: lon32(:)
      real*8, allocatable :: lon64(:)
      real, allocatable :: lev(:)
      real, allocatable :: aux(:,:,:,:)

      verbose_ = verb
      if(present(verbose)) verbose_ = verbose

!     Open GFIO file
!     --------------
      im_max   = im
      jm_max   = jm
      mlev_max = mlev
      usrlevs = .true.
      if(present(stat)) stat = 0
      call GFIO_Open ( trim(gfiofile), 0, ncid, ier )

      if ( ier .eq. 0 ) then

         call GFIO_DimInquire ( ncid, im_max, jm_max, mlev_max, lm_max,
     &                          nvars_in, ngatts, ier)

         if ( nvars .gt. nvars_in ) then
             if(verbose_) write(stdout,'(a,2(a,i5))') myname_,
     &                    ': Warning, no. of variable in file: ', nvars_in,
     &                    ' smaller than requested no. of vars: ', nvars
         end if

!        Check consistency of Gfio file
!        ------------------------------
         if (im_max .ne. im) then
            if(verbose_) write(stderr,'(2a,i5)') myname_,
     &           ': trying to read unallowed longitudinal dim'
            if ( present(stat) ) then
                 stat = 1
                 return
            else
                 call exit(7)
            endif
         endif

         if (jm_max .ne. jm) then
            if(verbose_) write(stderr,'(2a,i5)') myname_,
     &           ': trying to read unallowed latitudinal dim'
            if ( present(stat) ) then
                 stat = 1
                 return
            else
                 call exit(7)
            endif
         endif

         if (mlev_max .ne. mlev) then
            if(verbose_) write(stdout,'(2a,i5)') myname_,
     &           '-WARNING: trying to read unallowed number of levels'
            usrlevs = .false.
         endif

      else

         if(verbose_) write(stdout,'(2a)') myname_,
     .                ': No GFIO file for current date '
         call flush (stdout)
         if ( present(stat) ) then
              stat = 2
              return
         else
              call exit(7)
         endif
      endif

!     OK, let's read data from file
!     =============================
      icount = 0
      n2dim  = im * jm
      n3dim  = n2dim * mlev

!     Get longitude array to check on GEOS-5 flipped fields
!     -----------------------------------------------------
      allocate ( lon64(im), stat=ier )
        if (ier/=0) then
            if(verbose_) write(stderr,'(2a,i10,2a)') myname_, 'Alloc(lon) error '
            call flush (stdout)
            call exit(7)
        endif

      call ncdinq (ncid, 1, dimName, dimSize, ier)
         if (ier/=0)then
             if(verbose_) write(stderr,'(2a)') myname_, 'Error(ncdinq)'
             call flush (stderr); call exit(7)
         endif
      dimId = ncvid (ncid, dimName, ier)
         if (ier/=0)then
             if(verbose_) write(stderr,'(2a)') myname_, 'Error(ncvid)'
             call flush (stderr); call exit(7)
         endif

!     Caution: will try to read as 64-bit first ...
!     ---------------------------------------------
      call ncvgt ( ncid, dimId, 1, im, lon64, ier )
           if ( ier .ne. 0 ) then
             if(verbose_) write(stderr,'(2a,i10,2a)') myname_,
     &                    ': Gfio_GetVar, ier = ', ier,
     &                    ' for variable ', 'lon'
             call flush(stderr); call exit(7)
            endif

!     ... if absurd lons, read as 32-bit
!     ----------------------------------
      mnlon = minval(lon64)
      mxlon = maxval(lon64)
      if ( mnlon < -200.0 .or. mxlon > 400.0 ) then
         allocate ( lon32(im), stat=ier )
         call ncvgt ( ncid, dimId, 1, im, lon32, ier )
         lon64=lon32
         deallocate ( lon32 )
      endif

!     Is orientation of fields shifted by 180 degs from internal def
!     -------------------------------------------------------------
      fliplon = .false.
      if ( abs(lon64(1)+180.0)<rtol ) then
        fliplon = .true.  ! orientation of fields shifted by 180 degs from internal def
        print *, ' Flipping Fields Longitudinally ... '
      endif
      deallocate ( lon64 )

!     Now read upper-air fields
!     -------------------------
      if ( usrlevs ) then

        do iv = 1, nvars
          call Gfio_GetVar ( ncid, vname(iv), nymd, nhms,
     &                     im, jm, 1, mlev, fld(1,1,1,iv), ier )
             if ( ier .ne. 0 ) then
               if(verbose_) write(stderr,'(2a,i10,2a)') myname_,
     &                      ': Gfio_GetVar  ier = ', ier,
     &                      ' for variable ', trim(vname(iv))
               if ( present(stat) ) then
                    stat = 3
                    return
               else
                    call exit(7)
               endif
             else
               call fix_undef_in_ ( n3dim, fld(:,:,:,iv:iv), vname(iv) )
             end if
             if( fliplon ) call hflip_ ( fld(1:im,1:jm,1:mlev,iv), im, jm, mlev )
             icount = icount + 1
        end do

      else

        allocate ( aux(im,jm,mlev_max,1), lev(mlev_max), stat=ier )
             if ( ier .ne. 0 ) then
                  if(verbose_) write(stderr,'(2a,i10)') myname_,
     .                         ': Alloc(aux), ier = ', ier
                  if ( present(stat) ) then
                       stat = 4
                       return
                  else
                       call exit(7)
                  endif
             end if
        call Gfio_GetVar ( ncid, 'Height', nymd, nhms,
     &                     mlev_max, 1, 0, 1, lev, ier )
           if ( ier .ne. 0 ) then
             if(verbose_) write(stderr,'(2a,i10)') myname_,
     &                    ': Gfio_GetVar (lev), ier = ', ier
             if ( present(stat) ) then
                  stat = 3
                  return
             else
                  call exit(7)
             endif
           end if

        do iv = 1, nvars
           call Gfio_GetVar ( ncid, vname(iv), nymd, nhms,
     &                        im, jm, 1, mlev_max, aux, ier )
           if ( ier .ne. 0 ) then
             if(verbose_) write(stderr,'(2a,i10,2a)') myname_,
     &                    ': Gfio_GetVar, ier = ', ier,
     &                    ' for variable ', trim(vname(iv))
             if ( present(stat) ) then
                  stat = 3
                  return
             else
                  call exit(7)
             endif
           else
             call fix_undef_in_ ( im*jm*mlev_max, aux, vname(iv) )
           end if
           if( fliplon ) call hflip_ ( fld(1:im,1:jm,1:mlev_max,iv), im, jm, mlev_max )

           levcnt = 0
           fld(:,:,:,iv:iv) = undef
           if (present(plevs) ) then
               do k = 1, mlev
                  kloc = maxloc(lev(1:mlev_max),
     .                   mask=abs(log(plevs(k))-log(lev(1:mlev_max))).lt.1e-5)
                  j    = kloc(1)
                  if (j>0) then
                     fld(:,:,k:k,iv:iv) = aux(:,:,j:j,1:1)
                     levcnt = levcnt + 1
                  end if
               end do
               if ( levcnt.ne.mlev ) then
                    if(verbose_) write(stderr,'(a,2(i10))') myname_,
     &                      'WARNING: Requested levs not found, levcnt = ',
     &                      levcnt, ' levs requested ', mlev
               end if
          else
               fld(:,:,1:mlev,iv:iv) = aux(:,:,1:mlev,1:1)
          endif  ! < present(plevs) >
          icount = icount + 1
        end do

        deallocate ( aux, lev, stat=ier )

      end if

!     Close GFIO file now
!     -------------------
      call Gfio_Close ( ncid, ier )
      if ( ier .ne. 0 ) then
           if(verbose_) write(stdout,'(2a,i7)') myname_,
     &                  ': WARNING, problems closing GFIO file, ier = ', ier
      end if

      if(verbose_) write(stdout,'(2a,i12,2a)') myname_,
     &            ':  read ', icount, ' variables from ', trim(gfiofile)
      if(verbose_) write(stdout,'(2a,i8.8,a,i6.6)') myname_,
     &            ': nymd, nhms = ', nymd,' ', nhms

!     All done
!     --------
      return
      end subroutine GFIO_Get_3d_
!EOC

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE:  Fix_Undef_In_ - Fixes undef values when input is 32 bit
!
! !INTERFACE:
!
      subroutine Fix_Undef_In_ ( ndim, var, vname, verbose )
!
! !USES:

      use m_fpe, only: isnan

      implicit none

! !INPUT PARAMETERS:

      integer  ndim
      character*(*) vname
      logical, intent(in), optional ::  verbose

! !INPUT/OUTPUT PARAMETERS:

      real     var(ndim)

! !DESCRIPTION:
!
! !REVISION HISTORY:
!
!  13Apr00   Todling    - Initial code.
!  23Jun2010 Kokron     - eliminate previous calculation of val_max_in
!                         because it would fail if the input field
!                         contained NaNs
!
!EOP
!-------------------------------------------------------------------------

      character(len=*), parameter :: myname_ = myname//'*'//'Fix_Undef_In_'

      integer  i
      integer  nonmiss, numnan
      real     tol
      real     val_max_in, val_max_out
      logical  verbose_

      tol   = 0.001 * undef

      verbose_ = verb
      if ( present(verbose) ) verbose_ = verbose

!     Make sure there is no precision problem with undef's
!     ----------------------------------------------------
      nonmiss = 0; numnan = 0.; val_max_in=0.
      do i= 1,ndim
         if (isnan(var(i))) then
            numnan=numnan+1
            var(i) = undef
         elseif ( abs(var(i)-undef) .lt. tol ) then
            nonmiss = nonmiss + 1
            var(i) = undef
         end if
      end do
      val_max_out = maxval(var)
      if (nonmiss.ne.0) then
         if(verbose_)
     &   write(stdout,'(2a,i10,3a)') myname_, ': Fixed ', nonmiss,
     .                             ' values from input with strange undef',
     .                             ' for variable ', trim(vname)
         call flush(stdout)
         if ( val_max_out .ne. undef ) then
            write(stdout,*) ' Largest  value on  input: ',  val_max_in
            write(stdout,*) ' Largest  value on output: ',  val_max_out
            write(stdout,*) ' Undef        value spec.: ',  undef
            write(stdout,*) ' Correction not done. Aborting ... '
            call exit (7)
         end if
      end if
      if (numnan .ne. 0) then
         write(stdout,'(2a,i10,3a)') myname_, ': There were ',numnan,
     .                             ' NaN values from input',
     .                             ' for variable ', trim(vname)
         call flush(stdout)
      endif

      return
      end subroutine Fix_Undef_In_

      subroutine hflip3_ ( q,im,jm,km )
      implicit none
      integer  im,jm,km,i,j,k
      real, intent(inout) :: q(im,jm,km)
      real, allocatable   :: dum(:)
      allocate ( dum(im) )
      do k=1,km
      do j=1,jm
      do i=1,im/2
         dum(i) = q(i+im/2,j,k)
         dum(i+im/2) = q(i,j,k)
      enddo
         q(:,j,k) = dum(:)
      enddo
      enddo
      deallocate ( dum )
      end subroutine hflip3_

      subroutine hflip2_ ( q,im,jm )
      implicit none
      integer  im,jm,i,j
      real, intent(inout) :: q(im,jm)
      real, allocatable   :: dum(:)
      allocate ( dum(im) )
      do j=1,jm
      do i=1,im/2
         dum(i) = q(i+im/2,j)
         dum(i+im/2) = q(i,j)
      enddo
         q(:,j) = dum(:)
      enddo
      deallocate ( dum )
      end subroutine hflip2_

      end module m_GFIO_GetFld


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.

      subroutine heatfl ( pmsa   , fl     , ipoint , increm, noseg ,
     &                    noflux , iexpnt , iknmrk , noq1  , noq2  ,
     &                    noq3   , noq4   )
!>\file
!>       Total heat flux for surface water absolute temperature model

!
!     Description of the module :
!
! Name    T   L I/O   Description                                    Units
! ----    --- -  -    -------------------                            -----
! beta    R*4 1 L Bowens constant                                [mbar/°C]
! cloud   R*4 1 L Cloud coverage (fraction)                            [-]
! cpa     R*4 1 I Specific heat capacity of air                   [J/kg,K]
! cp      R*4 1 I Specific heat capacity of water                [J/kg,°C]
! CwindA  R*4 1 I Coefficient in wind function a                       [-]
! CwindB  R*4 1 I Coefficient in wind function b                       [-]
! CwindC  R*4 1 I Coefficient in wind function c                       [-]
! DeltaT  R*4 1 I Correction factor for TempWa near discharge          [-]
! Depth   R*4 1 I Water depth of segment                               [m]
! dTemp   R*4 1 O Flux on water temperature                         [°C/d]
! Emiss   R*4 1 L Emissivity (colour factor) for the atmosphere        [-]
! EWater  R*4 1 I Emissivity (colour factor) of surface water          [-]
! Fa      R*4 1 I Reflection factor for Qa                             [-]
! Fsw     R*4 1 I Reflection factor for Qsw                            [-]
! Fwind   R*4 1 L Wind function for evaporation                 [m/s,mbar]
! Fwind2  R*4 1 L Wind function for evaporation, incl corr.    [W/m2,mbar]
! Ha      R*4 1 I Height of wind speed needed for wind function        [m]
! Hm      R*4 1 I Height at which wind speed is measured               [m]
! HTVap   R*4 1 L Heat of evaporation                               [J/kg]
! HTVap   R*4 1 L Heat of evaporation at reference temperature      [J/kg]
! K       R*4 1 I Roughness coef for calculation vwind                 [-]
! Patm    R*4 1 I Air pressure                                 [mb of HPa]
! Pvap    R*4 1 L Actual vapour pressure in air                     [mbar]
! Psvap   R*4 1 L Saturated vapour pressure in air                  [mbar]
! PvapWa  R*4 1 L Saturated vapour pressure at water-air interface  [mbar]
! Qa      R*4 1 L Long wave atmospheric radiation reaching water    [W/m2]
! Qan     R*4 1 L Non refl. (net) atm. long wave rad. reaching water[W/m2]
! Qbr     R*4 1 L Long wave back radiation from water               [W/m2]
! Ql      R*4 1 L Latent heat flux by evaporation or condensation   [W/m2]
! Qrb     R*4 1 I Heat exchange between water and riverbed          [W/m2]
! Qsg     R*4 1 L Sensible heat of conduction                       [W/m2]
! Qsn     R*4 1 L Non refl. (net) short wave rad. reaching water    [W/m2]
! Qsw     R*4 1 L Short wave radiation reaching water               [W/m2]
! Qt      R*4 1 L Total heat flux                                   [W/m2]
! RatBow  R*4 1 L Bowens ratio                                         [-]
! RelHum  R*4 1 I Relative air humidity                                [%]
! Rho0    R*4 1 I Density of surface water at 4°C                  [kg/m3]
! RhoRef  R*4 1 L Density of surface water at reference temperature[kg/m3]
! SBC     R*4 1 I StefanBolzman Constant                        [W/m2/K-4]
! SunFac  R*4 1 I Percentage sunshine                                  [%]
! SWEmis  R*4 1 I Switch for calculation of emissivity                 [-]
! TempAt  R*4 1 I Air temperature                                     [°C]
! TempWa  R*4 1 I Surface water temperature                           [°C]
! Tref    R*4 1 I Reference water temperature for wind function       [°C]
! Vevap   R*4 1 L Evaporation rate                                   [m/s]
! VWinda  R*4 1 L Calculated wind speed                              [m/s]
! VWindm  R*4 1 I Measured wind speed                                [m/s]
!
!     Logical Units : -

!     Modules called : -

!     Name     Type   Library
!     ------   -----  ------------

      IMPLICIT NONE
      REAL     PMSA  ( * ) , FL  (*)
      INTEGER  IPOINT( * ) , INCREM(*) , NOSEG , NOFLUX,
     +         IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4,
     +         ISEG, IFLUX, IKMRK1, IKMRK2
      INTEGER  IP1, IP2, IP3, IP4, IP5, IP6, IP7, IP8, IP9, IP10,
     +         IP11, IP12, IP13, IP14, IP15, IP16, IP17, IP18, IP19,
     +         IP20, IP21, IP22, IP23, IP24, IP25, IP26, IP27, IP28,
     +         IP29, IP30, IP31, IP32, IP33, IP34, IP35, IP36, IP37,
     +         IP38, IP39, IP40, IP41
      INTEGER  ISWTEMP, ISWEMIS
      REAL     Qsw, Fsw, SunFac, Pvap, TempAt, SBC, Fa, EWater,
     j         TempWa, RhoWat, CWindA, CWindB, CWindC, VWindm, K,
     j         PvapWa, Cpa, Patm, Qsn, Cloud, Emiss, Qa, Qan, Qbr, C,
     j         HTVap, Fwind, Ql, Beta, RatBow, Qsg, Qt, Qrb, Tref,
     j         dTemp, cp, DEPTH, RelHum, Psvap, Hm, Ha, Vevap, Fwind2,
     j         VWinda, DeltaT, HtVRef, RhoRef, Rho0, MODTEMP, mindeptht

      IP1   = IPOINT( 1)
      IP2   = IPOINT( 2)
      IP3   = IPOINT( 3)
      IP4   = IPOINT( 4)
      IP5   = IPOINT( 5)
      IP6   = IPOINT( 6)
      IP7   = IPOINT( 7)
      IP8   = IPOINT( 8)
      IP9   = IPOINT( 9)
      IP10  = IPOINT(10)
      IP11  = IPOINT(11)
      IP12  = IPOINT(12)
      IP13  = IPOINT(13)
      IP14  = IPOINT(14)
      IP15  = IPOINT(15)
      IP16  = IPOINT(16)
      IP17  = IPOINT(17)
      IP18  = IPOINT(18)
      IP19  = IPOINT(19)
      IP20  = IPOINT(20)
      IP21  = IPOINT(21)
      IP22  = IPOINT(22)
      IP23  = IPOINT(23)
      IP24  = IPOINT(24)
      IP25  = IPOINT(25)
      IP26  = IPOINT(26)
      IP27  = IPOINT(27)
      IP28  = IPOINT(28)
      IP29  = IPOINT(29)
      IP30  = IPOINT(30)
      IP31  = IPOINT(31)
      IP32  = IPOINT(32)
      IP33  = IPOINT(33)
      IP34  = IPOINT(34)
      IP35  = IPOINT(35)
      IP36  = IPOINT(36)
      IP37  = IPOINT(37)
      IP38  = IPOINT(38)
      IP39  = IPOINT(39)
      IP40  = IPOINT(40)
      IP41  = IPOINT(41)
!
      IFLUX = 0
      DO 9000 ISEG = 1 , NOSEG

!     Heat exchange only for active water segments
!
!!    CALL DHKMRK(1,IKNMRK(ISEG),IKMRK1)
!!    IF (IKMRK1.EQ.1) THEN
      IF (BTEST(IKNMRK(ISEG),0)) THEN
!
!     Heat exchange only for top layer segments
!
         CALL DHKMRK(2,IKNMRK(ISEG),IKMRK2)
         IF (IKMRK2.EQ.0 .OR. IKMRK2.EQ.1) THEN
!
            Qsw     = PMSA(IP1 )
            Fsw     = PMSA(IP2 )
            SunFac  = PMSA(IP3 )
            ISWEmis = NINT(PMSA(IP4 ))
            TempAt  = PMSA(IP5 )
            SBC     = PMSA(IP6 )
            Fa      = PMSA(IP7 )
            Ewater  = PMSA(IP8 )
            TempWa  = PMSA(IP9 )
            Rho0    = PMSA(IP10)
            CWinda  = PMSA(IP11)
            CWindb  = PMSA(IP12)
            CWindc  = PMSA(IP13)
            VWindm  = PMSA(IP14)
            hm      = PMSA(IP15)
            ha      = PMSA(IP16)
            k       = PMSA(IP17)
            cpa     = PMSA(IP18)
            Patm    = PMSA(IP19)
            Qrb     = PMSA(IP20)
            Cp      = PMSA(IP21)
            DEPTH   = PMSA(IP22)
            RelHum  = PMSA(IP23) / 100.0
            DeltaT  = PMSA(IP24)
            Tref    = PMSA(IP25)
            ISWTEMP = NINT(PMSA(IP26))
            mindeptht = PMSA(IP27)

!
!     ------Short wave solar radiation-----------
!
            Qsn = Qsw * (1- Fsw)
!
!     ------Long wave atmospheric radiation------
!
            cloud = (100.0 - SunFac) / 100.0
!
            Psvap = 6.131 + 0.467 * TempAt + 0.0089 *
     j               TempAt ** 2.0 + 0.000527 * TempAt ** 3.0
!
            Pvap = RelHum * Psvap
!
!           ----1.  Option 1 = Brunt, 1932
!
            IF (ISWEmis.EQ.1) THEN
!
               Emiss = (0.51 + 0.066 * SQRT(Pvap)) *
     j                 (1.0 + 0.17 * cloud ** 2.0)
!
               Qa = Emiss * SBC * (TempAt + 273.15) ** 4.0
!
!           ----2.  Option 2 = Edinger, 1965 and Koberg, 1962
!
            ELSEIF (ISWEmis.EQ.2) THEN
!
!           Pascal heeft formulering voor C (Brunt coefficient)
!           aangevraagd, wordt later ingevuld...
!
               C = 1
!
               Emiss = 1.1 * C + 0.030 * SQRT(Pvap)
!
               Qa = Emiss * SBC * (TempAt + 273.15) ** 4.0
!
!           ----3.  Option 3 = Edinger, 1965
!
            ELSEIF (ISWEmis.EQ.3) THEN
!
               Emiss = 0.74 * (1 + 0.17 * cloud) +
     j                 0.0045 * (1 - 0.4 * cloud) *
     j                 Pvap
!
               Qa = Emiss * SBC * (TempAt + 273.15) ** 4.0
!
!           ----4.  Option 4 = Ludikhuize, 1994 as in WAQUA
!
            ELSEIF (ISWEmis.EQ.4) THEN
!
               Qa = (218.0 + 6.3 * TempAt) *
     j                 (1.0 + 0.17 * cloud ** 2.0)
!
            ELSE
               WRITE (*,*) ' Illegal option for emissivity formula'
               CALL SRSTOP(1)
            ENDIF
!
            Qan = Qa * (1 - Fa)
!
!     ------Long wave back radiation from water------
!
            Qbr = Ewater * SBC * ((TempWa + DeltaT + 273.15) ** 4.0)
!
!     ------Latent heat (evaporation or condensation)-----
!
            RhoWat = Rho0 * (1.0 - 7.17e-6 *
     j               (TempWa + DeltaT - 4.0) ** 2.0)
!
            RhoRef = Rho0 * (1.0 - 7.17e-6 *
     j               (Tref - 4.0) ** 2.0)
!
            HtVap = 2.5e+06 - 2300.0 * (TempWa + DeltaT)
!
            HtVRef = 2.5e+06 - 2300.0 * Tref
!
            PvapWa = 6.131 + 0.467 * (TempWa + DeltaT) + 0.0089 *
     j               (TempWa + DeltaT) ** 2.0 + 0.000527 *
     j               (TempWa + DeltaT) ** 3.0
!
            IF (ABS(Hm - Ha) .LT. 0.0001) THEN
               VWinda = VWindm
            ELSE
               IF  ((Hm .GT. (10 * k)) .AND. (Ha .GT. (10 * k))) THEN
                   VWinda = VWindm * LOG(Hm / k) / LOG(Ha / k)
               ELSE
                   VWinda = VWindm
               ENDIF
            ENDIF
!
            FWind = CwindA + CwindB * VWinda + CwindC * VWinda ** 2.0
!
            FWind2 = Fwind / (RhoRef * HtVRef)
!
            Vevap = FWind2 * (PvapWa - Pvap)
!
            Ql = RhoWat * HtVap * Vevap
!
!     ------Convective heat-----------------------------
!
            beta = (cpa * Patm) / (0.62 * HtVap)
!
            Qsg = beta * (TempWa + DeltaT - TempAt) * RhoWat * HtVap * Fwind2
!
!     ------Convective heat to river bed----------------
!           formulations can be included later...
!
!     ------Total heat flux-----------------------------
!
            Qt = Qsn + Qan - Qbr - Ql - Qsg + Qrb
!
!     ------Change of water temperature-----------------
!
            dTemp = (Qt * 86400.0) / (RhoWat * Cp * Depth)
!
!     ------For output to exces temperature moddeling
!
            IF ( ISWTEMP .EQ. 0 ) THEN
               MODTEMP = TEMPWA
            ELSE
               MODTEMP = 0.0
            ENDIF

            PMSA(IP28 ) = Qsn
            PMSA(IP29 ) = cloud
            PMSA(IP30 ) = Psvap
            PMSA(IP31 ) = Emiss
            PMSA(IP32 ) = Qan
            PMSA(IP33 ) = Qbr
            PMSA(IP34 ) = Ql
            PMSA(IP35 ) = Qsg
            PMSA(IP36 ) = Qt
            PMSA(IP37 ) = PvapWa
            PMSA(IP38 ) = VWinda
            PMSA(IP39 ) = Vevap
            PMSA(IP40 ) = RhoWat
            PMSA(IP41 ) = MODTEMP

            IF ( depth .lt. mindeptht ) THEN

               ! depth less then mindeptht no fluxes

               dTemp = 0.0

            ENDIF
!
         ELSE
!
!           For non top layer segments no heat flux is calculated
!
            dTemp = 0.0
!
         ENDIF
!
      ENDIF
!
      FL( 1+IFLUX) = dTemp

      IFLUX = IFLUX + NOFLUX
      IP1   = IP1   + INCREM (  1 )
      IP2   = IP2   + INCREM (  2 )
      IP3   = IP3   + INCREM (  3 )
      IP4   = IP4   + INCREM (  4 )
      IP5   = IP5   + INCREM (  5 )
      IP6   = IP6   + INCREM (  6 )
      IP7   = IP7   + INCREM (  7 )
      IP8   = IP8   + INCREM (  8 )
      IP9   = IP9   + INCREM (  9 )
      IP10  = IP10  + INCREM ( 10 )
      IP11  = IP11  + INCREM ( 11 )
      IP12  = IP12  + INCREM ( 12 )
      IP13  = IP13  + INCREM ( 13 )
      IP14  = IP14  + INCREM ( 14 )
      IP15  = IP15  + INCREM ( 15 )
      IP16  = IP16  + INCREM ( 16 )
      IP17  = IP17  + INCREM ( 17 )
      IP18  = IP18  + INCREM ( 18 )
      IP19  = IP19  + INCREM ( 19 )
      IP20  = IP20  + INCREM ( 20 )
      IP21  = IP21  + INCREM ( 21 )
      IP22  = IP22  + INCREM ( 22 )
      IP23  = IP23  + INCREM ( 23 )
      IP24  = IP24  + INCREM ( 24 )
      IP25  = IP25  + INCREM ( 25 )
      IP26  = IP26  + INCREM ( 26 )
      IP27  = IP27  + INCREM ( 27 )
      IP28  = IP28  + INCREM ( 28 )
      IP29  = IP29  + INCREM ( 29 )
      IP30  = IP30  + INCREM ( 30 )
      IP31  = IP31  + INCREM ( 31 )
      IP32  = IP32  + INCREM ( 32 )
      IP33  = IP33  + INCREM ( 33 )
      IP34  = IP34  + INCREM ( 34 )
      IP35  = IP35  + INCREM ( 35 )
      IP36  = IP36  + INCREM ( 36 )
      IP37  = IP37  + INCREM ( 37 )
      IP38  = IP38  + INCREM ( 38 )
      IP39  = IP39  + INCREM ( 39 )
      IP40  = IP40  + INCREM ( 40 )
      IP41  = IP41  + INCREM ( 41 )
!
 9000 CONTINUE
!

      RETURN
      END

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine eabxag: calculate energy/gradient of single halogen bond
!          (numerical gradient)
!
!     part of QMDFF
!
subroutine eabxag(n,A,B,H,xyz,ca,energy,gdr)
use debug
use qmdff
implicit none
integer::A,B,H,n
real(kind=8)::xyz(3,n),ca,energy,gdr(3,n)
real(kind=8)::cosabh,aterm,rdampl
real(kind=8)::rab,rah,rbh,rab2,rah2,rbh2,rah4,rbh4
real(kind=8)::drah(3),drbh(3),drab(3)
real(kind=8)::ga(3),gb(3),gh(3),dg(3),dg2(3)
real(kind=8)::gi,denom,ratio,apref,aprod
real(kind=8)::eabh
REAL(kind=8)::longcut=120.d0
REAL(kind=8)::alp2= 6.d0
REAL(kind=8)::alp3=6.d0
real(kind=8)::step,er,el,dum
integer::i,j,i1,i2
!   The newly added gradient components
real(kind=8)::g_local_a(3),g_local_b(3),g_local_c(3)
real(kind=8)::vec_ah(3),vec_bh(3)

i1=A
i2=B
i =H

call eabx(n,i1,i2,i,xyz,ca,er)
energy=energy+er

!
!     Calculate components of single distance vectors
!
vec_ah=xyz(:,i1)-xyz(:,i)
vec_bh=xyz(:,i2)-xyz(:,i)
!
!     If debugging is activated: store single energy parts
!
!if (do_debug) then
!   if (qmdff_num.eq.1) then
!      comp_no=comp_no+1
!      act_part=er
!      qmdff_parts(1,struc_no,comp_no)=act_part
!      write(parts_labels(1,comp_no),'(a,i4,i4,i4,i4)') "x-bond",A,B,H
!  else if (qmdff_num.eq.2) then
!      comp_no2=comp_no2+1
!      act_part=er
!      qmdff_parts(2,struc_no,comp_no2)=act_part
!      write(parts_labels(2,comp_no2),'(a,i4,i4,i4,i4)') "x-bond",A,B,H
!   end if
!end if

!
!     calculate numerical gradient: elogante coordinates!
!
step=1.d-6
dum=1./(2.0d0*step)
do j=1,3
   xyz(j,i1)=xyz(j,i1)+step
   call eabx(n,i1,i2,i,xyz,ca,er)
   xyz(j,i1)=xyz(j,i1)-step*2.0d0
   call eabx(n,i1,i2,i,xyz,ca,el)
   xyz(j,i1)=xyz(j,i1)+step
   g_local_a=(er-el)*dum
end do
gdr(:,i1)=gdr(:,i1)+g_local_a
do j=1,3
   xyz(j,i2)=xyz(j,i2)+step
   call eabx(n,i1,i2,i,xyz,ca,er)
   xyz(j,i2)=xyz(j,i2)-step*2.0d0
   call eabx(n,i1,i2,i,xyz,ca,el)
   xyz(j,i2)=xyz(j,i2)+step
   g_local_c=(er-el)*dum
end do
gdr(:,i2)=gdr(:,i2)+g_local_c

do j=1,3
   xyz(j,i )=xyz(j,i )+step
   call eabx(n,i1,i2,i,xyz,ca,er)
   xyz(j,i )=xyz(j,i )-step*2.0d0
   call eabx(n,i1,i2,i,xyz,ca,el)
   xyz(j,i )=xyz(j,i )+step
   g_local_b=(er-el)*dum
end do
gdr(:,i)=gdr(:,i)+g_local_b

!
!     Calculate the virial tensor components, if needed!
!

if (calc_vir) then
  vir_ten(1,1)=vir_ten(1,1)+vec_ah(1)*g_local_a(1)+vec_bh(1)*g_local_b(1)
  vir_ten(2,1)=vir_ten(2,1)+vec_ah(2)*g_local_a(1)+vec_bh(2)*g_local_b(1)
  vir_ten(3,1)=vir_ten(3,1)+vec_ah(3)*g_local_a(1)+vec_bh(3)*g_local_b(1)
  vir_ten(1,2)=vir_ten(1,2)+vec_ah(1)*g_local_a(2)+vec_bh(1)*g_local_b(2)
  vir_ten(2,2)=vir_ten(2,2)+vec_ah(2)*g_local_a(2)+vec_bh(2)*g_local_b(2)
  vir_ten(3,2)=vir_ten(3,2)+vec_ah(3)*g_local_a(2)+vec_bh(3)*g_local_b(2)
  vir_ten(1,3)=vir_ten(1,3)+vec_ah(1)*g_local_a(3)+vec_bh(1)*g_local_b(3)
  vir_ten(2,3)=vir_ten(2,3)+vec_ah(2)*g_local_a(3)+vec_bh(2)*g_local_b(3)
  vir_ten(3,3)=vir_ten(3,3)+vec_ah(3)*g_local_a(3)+vec_bh(3)*g_local_b(3)
end if

return
end subroutine eabxag

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
	SUBROUTINE WN_DCOD ( curtim, gemfil, stntbl, prmfil,
     +			      iadstn, maxtim, nhours, iret )
C************************************************************************
C* WN_DCOD								*
C*									*
C* This subroutine decodes tornado, thunderstorm, snow squall, and      *
C* flash flood warning bulletins and writes the data to an ASCII file.  *
C*									*
C* WN_DCOD ( CURTIM, GEMFIL, STNTBL, PRMFIL, IADSTN, MAXTIM, NHOURS,	*
C*		IRET )							*
C*									*
C* Input parameters:							*
C*	CURTIM		CHAR*		Current time for input data	*
C*	GEMFIL		CHAR*		Output file name template	*
C*	STNTBL		CHAR*		Station table			*
C*	PRMFIL		CHAR*		Parameter table file name	*
C*	IADSTN		INTEGER		Number of additional stations	*
C*	MAXTIM		INTEGER		Number of times allowed		*
C*	NHOURS		INTEGER		Number of hours prior to CURTIM	*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					 -4 = Unable to get bultim	*
C*									*
C**									*
C* Log:									*
C* A. Hardy/GSC          5/99						*
C* A. Hardy/GSC          5/99	Added error check for dattim		*
C* A. Hardy/GSC          5/99	Added wtype to if statements		*
C* A. Hardy/GSC          6/99	Reworked and cleaned up			*
C* A. Hardy/GSC          2/00	Increased cnties 65->450		*
C* A. Hardy/GSC          2/00	Added check for header error		*
C* D. Kidwell/NCEP	 3/00	Added check for TEST bulletins		*
C* A. Hardy/GSC          3/00	Using bultim to create filename		*
C* A. Hardy/GSC          5/00	Added DC_GHDR call			*
C* A. Hardy/GSC         11/00   Changed FFLD,TORN,STSM -> STFS  	*
C* A. Hardy/GSC         12/00   Cleaned up ierr values			*
C* F. J. Yen/NCEP	 1/01	Replaced WN_RDTB with DC_STNS		*
C* F. J. Yen/NCEP	 3/02	Checked issue time & added DC_WLOG call.*
C* B. Yin/SAIC           3/04   Changed SS_GTIM to CSS_GTIM		*
C* D. Kidwell/NCEP	11/05	Improved error checks			*
C* T. Piper/SAIC	06/07	Extract warning polygon from bulletin	*
C* T. Piper/SAIC	08/07	Format change; ignore TIME...MOT...LOC	*
C* F. J. Yen/NCEP	 3/08	Added ETN; made polygon search robust;  *
C*				added correction letter to cor/test flag*
C* S. Guan/NCEP         04/15   Increased sections (starry) from 5 to 15*
C* S. Guan/NCEP         11/17   Modified to add snow squall warn (SQW)  *
C* S. Guan/NCEP         05/18   Initialized bultin                      *
C* B. Hebbard/NCEP      04/20   NAWIPS-128: Ignore new IBW tags after   *
C*                              warning polygon point list (SCN19-60)   *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'BRIDGE.PRM'
C*
	CHARACTER*(*)	curtim, gemfil, prmfil, stntbl
C*
	CHARACTER	bultin*(DCMXBF), parms(MMPARM)*4, errstr*80,
     +			starry(15)*2800, filnam*132, wtype*3,
     +			wfrom*3, dattim*11, bultim*12, wmohdr*8,
     +			oristn*8, seqnum*4, cnties*700, edttim*6,
     +			sysdt*12, dattmp*12, county(LLSTFL)*6,
     +			strtim*11, stptim*11, bbb*4, poly*512, etn*4
	CHARACTER	adstn(LLSTFL)*8, stnnam(LLSTFL)*32,
     +			stat(LLSTFL)*2, coun(LLSTFL)*2,
     +			tbchrs(LLSTFL)*20, tstchr*1
	INTEGER	istarr (15), irdtar (15)
	INTEGER	istnm(LLSTFL), ispri(LLSTFL), itype
	REAL		adlat(LLSTFL), adlon(LLSTFL), selv(LLSTFL)
	LOGICAL		more, done
C------------------------------------------------------------------------
	iret = 0
C
C*  Initialize open file lists. Set the max number of open files.
C*  Set the type of output file and file source.
C
	maxfil = 2
	iftype = 6
	CALL DC_FINT ( maxfil, iftype, prmfil, ier )
	iflsrc = 2
C
C*  Read the county table.
C
	CALL DC_STNS  ( stntbl, adstn, istnm, stnnam, stat, coun,
     +			adlat, adlon, selv, ispri, tbchrs, nade, ier )
	IF  ( ier .ne. 0 )  THEN
	    RETURN
	END IF
C
C*  Loop until a timeout occurs.
C
	iperr = 0
	DO WHILE ( iperr .eq. 0 )
	    icor = 0
C
C* Initialize bultin
C
            DO I = 1, DCMXBF
               bultin(I:I) =' '
            END DO
C
C*  Get the bulletin.
C
	    CALL DC_GBUL ( bultin, lenbul, ifdtyp, iperr )
	    IF ( iperr .eq. 0 ) THEN
		more = .true.
		CALL DC_TEST ( bultin, lenbul, itest, ier )
C
C*  Break bulletin into 15 strings.
C
		CALL ST_CLSL ( bultin, '*', ' ', 15, starry, inumb,
     +								ierr )
C
C*  Parse the header info from the bulletin.
C
		CALL DC_GHDR ( bultin, lenbul, seqnum, wmohdr,
     +			  oristn, bultim, bbb, nchar, ierr )
		IF ( ierr .ne. 0 ) THEN
		    CALL DC_WLOG ( 2, 'DC', ierr, ' ', ier )
		    more = .false.
		END IF
C
C*              To handle partial cancellation with two VTEC lines
C*              nchar indicates whether there is more than one vtec exist
C
                nchar = 1
                iierr = 0
                DO WHILE ( ( nchar .gt. 0 ) .and. ( iierr .eq. 0 ) )
                  IF ( nchar .gt. 1) THEN
                     islash = INDEX (bultin, '$$')
                     bultin = bultin(islash:lenbul)
                     CALL ST_CLSL ( bultin, '*', ' ', 15, starry,
     +                                               inumb, ierr )
                  END IF
		  IF ( more ) THEN
		    CALL WN_GHDR ( starry(1), icor, wtype, wfrom,
     +			           cnties, edttim, nchar, etn, iierr )
		    IF (  bbb .eq. 'COR'  ) THEN
			icor = 1
		      ELSE IF ( bbb(1:2) .eq. 'CC' ) THEN
C
C*			bbb(3:3) can be from 'A' through 'Z'.
C*			So icor will be from 65 through 90.
C
			icor = ichar (bbb(3:3))
		    END IF
		    IF ( ierr .lt. 0 ) THEN
			CALL ST_UNPR ( bultin(:80), 80, errstr, len1,
     +				       ier )
      			CALL DC_WLOG ( 2, 'DCWARN', ierr,
     +					errstr(:len1), ier )
			more = .false.
		    END IF
		    IF ( inumb .ge. 1 .and. inumb .le. 15 ) THEN
C*
C*			The polygon should be in the fifth (last),
C*			section.  If there is no fifth section, then
C*			search for the polygon in the last section.
C*
		        CALL ST_LSTR(starry(inumb), lens, ier)
		        CALL ST_UNPR(starry(inumb), lens, starry(inumb),
     +				     lenout, ier)
		        ibeg = INDEX(starry(inumb), 'LAT...LON')
		        IF ( ibeg .gt. 0 )  THEN
			    ibeg = ibeg + 10
C*                  "iend" needs to be index of end of polygon string.  Set it to point to character
C*                  *before* the first non-numeric, non-blank character encountered
                            iend = ibeg
                            done = .false.
                            DO WHILE ( .not. done )
                                tstchr = starry(inumb)(iend:iend)
C*                              If this character is numeric or a blank, go to next character
                                IF ( ((tstchr .ge. '0') .and.
     +                                (tstchr .le. '9')) .or.
     +                                (tstchr .eq. ' ') ) THEN
                                    iend = iend + 1
C*                              Otherwise we've found the first non-numeric, non-blank character,
C*                              so back up iend one to previous character (last numeric or blank)
                                ELSE
                                    iend = iend -1
                                    done = .true.
                                END IF
                            END DO

			    IF ( iend .gt. 0 )  THEN
			        poly = starry(inumb)(ibeg:iend)
			    ELSE
				ierr = 4
				mnln = MIN (lens,48)
			        CALL DC_WLOG ( 2, 'DCWARN', ierr,
     +					starry(inumb)(:mnln), ier )
			        more = .false.
			    END IF
		        ELSE
			    ierr = 4
			    mnln = MIN (lens,48)
			    CALL DC_WLOG ( 2, 'DCWARN', ierr,
     +					starry(inumb)(:mnln), ier )
			    more = .false.
		        END IF
		    END IF
		  END IF
C
		IF ( more ) THEN
C
C*		   For icor values .GE. 65, then the correction letter in
C*		   the bbb field was used in icor, so add 256 to it if it
C*		   is a test, else add 2 to make it compatible with the
C*		   old way that didn't use the bbb field.
C
		    IF ( itest .eq. 2 ) THEN
			IF ( icor .ge. 65 ) THEN
		            itest = 256
			END IF
		    END IF
		    icor = icor + itest
C
C*  Get the system time, and make a standard GEMPAK time from the
C*  "current" time.
C
		    itype = 1
		    CALL CSS_GTIM ( itype, sysdt, ier )
		    IF ( curtim .eq. 'SYSTEM' ) THEN
			dattmp = sysdt
		    ELSE
			CALL TI_STAN ( curtim, sysdt, dattmp, ier )
		    END IF
		    CALL TI_CTOI ( dattmp, istarr, ier )
C
C*  Set the GEMPAK date/time.  Use bultim to set the output file name.
C
		    CALL ST_LSTR ( bultim, len, ier )
		    IF ( len .ne. 6 ) THEN
			IF ( curtim .eq. 'SYSTEM' ) THEN
			    bultim = dattmp( 5:6 ) // dattmp( 8:11 )
			ELSE
			    iret = -4
			    CALL ST_UNPR ( bultin (:80), 80, errstr,
     +					   len1, ier )
			    CALL DC_WLOG ( 2, 'DCWARN', iret,
     +					   errstr(:len1), ier )
			    more = .false.
			END IF
		    END IF
		    IF ( more ) THEN
			CALL ST_NUMB ( bultim, issue, ier )
			irday  = issue / 10000
			irhour = mod ( issue, 10000 ) / 100
			irmin  = mod ( mod ( issue, 10000 ), 100 )
			CALL DC_ITIM ( istarr, irday, irhour, irmin,
     +				   irdtar, ier )
C
C*  If invalid values in issue time, use system or current time.
C
		        CALL TI_DAYM ( irdtar(1), irdtar(2), ndays,
     +				ier )
			IF ( irday  .lt. 1 .or. irday  .gt. ndays .or.
     +			     irhour .lt. 0 .or. irhour .gt. 23    .or.
     +			     irmin  .lt. 0 .or. irmin  .gt. 59 )  THEN
			    irday  = istarr(3)
			    irhour = istarr(4)
			    irmin  = istarr(5)
			END IF
			CALL DC_ITIM ( istarr, irday, irhour, irmin,
     +                             irdtar, ier )
			CALL TI_ITOC ( irdtar, dattim, ier )
			IF ( ier .lt. 0 ) more = .false.
C
C*  Remove control characters from entire bulletin.
C
		        IF ( more ) THEN
C
C*  Set bultim from dattim, in case original value was invalid.
C
			    bultim = dattim( 5:6 ) // dattim( 8:11 )
			    CALL ST_UNPR ( bultin, lenbul, bultin,
     +					lenb, ier )
C
C*  Make a file name from the template and the time.  Open the file as
C*  part of the open file list.
C
			    CALL FL_MNAM ( dattim, gemfil, filnam, ier )
			    CALL DC_FCYL ( filnam, iflsrc, stntbl,
     +				    iadstn, maxtim, lunf, nparm,
     +				    parms, ierr )
C
C*  Check that the file was opened properly.
C
			    IF ( ierr .ne. 0 ) THEN
C
C*  If not, write an error to the decoder log file.
C
				CALL DC_WLOG ( 0, 'SF', ierr, filnam,
     +					ier )
				more = .false.
			    END IF
C
			    IF ( more ) THEN
				IF ( ( wtype .eq. 'SVR' ) .or.
     +					( wtype .eq. 'TOR' ) .or.
     +                                  ( wtype .eq. 'SQW' ) .or.
     +					( wtype .eq. 'FFW' ) ) THEN
				    IF ( wtype .eq. 'FFW' ) THEN
					minalw = 480
					mindif = 360
				    ELSE
					minalw = 240
					mindif = 120
				    END IF
				    CALL WN_STFS ( bultim, istarr,
     +					      edttim, cnties, minalw,
     +					      mindif, bultin, 'DCWARN',
     +					      strtim, stptim, county,
     +					      numcty, iret )
				ELSE
				    more = .false.
				END IF
				IF ( ( wtype  .eq. ' ' ) .or.
     +				     ( strtim .eq. ' ' ) .or.
     +				( stptim .eq. ' ' ) )
     +				more = .false.
				IF ( more ) THEN
				   CALL WN_OUT (lunf, wtype, strtim,
     +					  stptim, oristn, icor, etn,
     +					  poly, county, numcty, adstn,
     +					  istnm, stnnam, stat, coun,
     +					  adlat, adlon, selv, ispri,
     +					  tbchrs, nade, iret )
				END IF
			    END IF
			END IF
		    END IF
		END IF
              END DO
	      CALL DC_FCLS ( ier )
	    END IF
	END DO
C
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program main
    use iso_fortran_env, only : output_unit
    write(output_unit,'(i2)') __PGIC__
end program

max_stars_repo_name: souchtal/fti

max_issues_repo_name: souchtal/fti

max_forks_repo_name: souchtal/fti


text:
!======================================================================!
  SUBROUTINE IniDiv
!----------------------------------------------------------------------!
!   Initialize processor numbers to 1                                  !
!----------------------------------------------------------------------!
!------------------------------[Modules]-------------------------------!
  USE all_mod
  USE par_mod
  USE div_mod
!----------------------------------------------------------------------!
  IMPLICIT NONE
!-------------------------------[Locals]-------------------------------!
  INTEGER :: c
!--------------------------------[CVS]---------------------------------!
!  $Id: IniDiv.f90,v 1.6 2002/10/30 16:29:18 niceno Exp $
!  $Source: /home/muhamed/.CVSROOT/T-Rex/Divide/IniDiv.f90,v $
!======================================================================!

!------------------------------!
!     Initialize parameters    !
!------------------------------!
  COORDINATE = 1
  INERTIAL   = 2

!------------------------------!
!     Initialize processors    !
!------------------------------!
  do c=1,NC
    proces(c)=1
  end do

  END SUBROUTINE IniDiv

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


module deq

  implicit none

  private

  public :: dfit_init, dfitin, gradient, eqitem, bgradient, dfit_finish
  public :: invR
  public :: Rpos
  public :: Zpos
  public :: btori,    initialize_btori
  public :: dbtori,   initialize_dbtori
  public :: qfun,     initialize_q
  public :: pfun,     initialize_pressure
  public :: dpfun,    initialize_dpressure
  public :: betafun,  initialize_beta
  public :: bound,    initialize_bound
  public :: psi,      initialize_psi
  public :: rhofun,   initialize_rho

  integer :: nw, nh, nbbbs, nrho

  real, allocatable, dimension (:) :: psi_bar, fp, qsf, pressure, beta, spsi_bar
  real, allocatable, dimension (:) :: dummy, dfit_R, dfit_Z, sdfit_R, sdfit_Z
  real, allocatable, dimension (:) :: rho_mid, psi_mid
  real, allocatable, dimension (:,:) :: dfit_psi, sdfit_psi
  real, allocatable, dimension (:,:,:) :: dpm, dtm
  real, allocatable, dimension (:) :: rbbbs, zbbbs, thetab, r_bound !boundary of plasma

  real :: dfit_dR, dfit_dZ, psi_N

  real :: R_mag, Z_mag, B_T0, aminor, beta_0

  logical :: init_bound = .true.
  logical :: init_btori = .true.
  logical :: init_dbtori = .true.
  logical :: init_q = .true.
  logical :: init_pressure = .true.
  logical :: init_dpressure = .true.
  logical :: init_beta = .true.
  logical :: init_rho = .true.
  logical :: init_psi = .true.
  logical :: init_R = .true.
  logical :: init_Z = .true.

contains
  subroutine dfitin(eqfile, psi_0, psi_a, rmaj, B_T, amin, initeq)
!
!     This subroutine reads an DFIT output file containing
!     the axisymmetric magnetic field geometry on a rectangular
!     domain defined by the coordinates (R,Z).
!
!     dfit_R     R grid
!     dfit_Z     Z grid
!     fp    F on psibar grid
!     dfit_psi   psi on (R,Z) grid
!     R_mag        major radius of the magnetic axis
!     Z_mag        elevation of the magnetic axis
!     rwid       total width of the domain
!     rleft      position of leftmost point of domain
!     zhei       total height of domain
!
    use mp, only: mp_abort
    use constants, only: pi, twopi
    use splines, only: inter_cspl
    implicit none

    character(len=80), intent(in) :: eqfile
    real, intent(out) :: psi_0, psi_a, B_T, amin, rmaj
    integer, intent(in) :: initeq
    real :: p_0
    real :: rwid, zhei
    real :: bcentr
    real :: delta_R, delta_Z
!    real :: fitp_surf2, ierr, rleft, zxy11, zxy1n, zxym1, zxymn
!    real, dimension(:), allocatable :: temp, zp, zx1, zxm, zy1, zyn
    character(80) :: filename

    integer :: i, j, init, nhb, nwb
    integer :: jmin, jmax
!What is the following for?
    data init /1/
    save init

! Need to generalize initialization condition if equilibrium changes

    if(initeq == 0) return
    init=0

    i=index(eqfile,' ')-1
    filename = eqfile(1:i)
    open(unit=5,file=filename,status='old',form='formatted')

! Read the data

    read(5,*) nw, nh, nrho

!    nwb = nw * big
!    nhb = nh * big

    nwb = nw
    nhb = nh

    call alloc_module_arrays(nwb, nwb, nhb, nw, nh, nrho)

    read(5,2020) rwid, zhei
    read(5,2020) R_mag, Z_mag
    read(5,2020) psi_0, psi_a, bcentr
!
! pbar is defined by
! pbar = (psi-psi_0)/(psi_a-psi_0)
! fp and q are functions of pbar
!
    do i=1,nw
       spsi_bar(i) = float(i) / float(nw)
       sdfit_R(i) = rwid * float(i) / float(nw)
    enddo

    do j=1,nh
       sdfit_Z(j) = zhei*(float(j-1)/float(nh-1)-0.5)
    enddo

    do i=1,nwb
       psi_bar(i) = float(i) / float(nwb)
       dfit_R(i) = rwid * float(i) / float(nwb)
    enddo

    do j=1,nhb
       dfit_Z(j) = zhei*(float(j-1)/float(nhb-1)-0.5)
    enddo

    read(5,2020) (dummy(j), j = 1, nw)  ! pressure read
    call inter_cspl(nw, spsi_bar, dummy, nwb, psi_bar, pressure)
    read(5,2020) ((dfit_psi(i,j) , i = 1, nw) , j = 1, nh)
!    read(5,2020) ((sdfit_psi(i,j) , i = 1, nw) , j = 1, nh)

!    allocate(zp(3*nw*nh), temp(nw+2*nh))
!    allocate(zx1(nh), zxm(nh), zy1(nw), zyn(nw))

!    call fitp_surf1(nw, nh, sdfit_R, sdfit_Z, sdfit_psi, &
!         nw, zx1, zxm, zy1, zyn, zxy11, zxym1, zxy1n, zxymn, &
!         255, zp, temp, 1., ierr)

    do j = 1, nhb
       do i = 1, nwb
!          dfit_psi(i,j) = fitp_surf2(dfit_R(i), dfit_Z(j), nw, nh, &
!               sdfit_R, sdfit_Z, sdfit_psi, nw, zp, 1.)
       enddo
    enddo

!    deallocate(zp, temp)

    read (5,2020) (rho_mid(i), i=1,nrho)
    read (5,2020) (psi_mid(i), i=1,nrho)

    fp = 0.
    qsf = 0.

    nw = nwb
    nh = nhb

    nbbbs = 2*(nw+nh)-5

! stay away from edge of box
    delta_R = 1.02 ; delta_Z = 1.02

    allocate(rbbbs(nbbbs), zbbbs(nbbbs), thetab(nbbbs), r_bound(nbbbs))
    jmin = 1
    jmax = nh/2
    do j=jmin, jmax
       rbbbs(j) = (dfit_R(1)-0.5*rwid)/delta_R+0.5*rwid
       zbbbs(j) = dfit_Z(j+nh/2)/delta_Z
    end do

    jmin = jmax + 1
    jmax = jmax + nw - 2
    do j=jmin, jmax
       rbbbs(j) = (dfit_R(j-jmin+2)-0.5*rwid)/delta_R+0.5*rwid
       zbbbs(j) = dfit_Z(nh)/delta_Z
    end do

    jmin = jmax + 1
    jmax = jmax + nh - 1
    do j=jmin, jmax
       rbbbs(j) = (dfit_R(nw)-0.5*rwid)/delta_R+0.5*rwid
       zbbbs(j) = dfit_Z(jmax-j+2)/delta_Z
    end do

    jmin = jmax + 1
    jmax = jmax + nw - 1
    do j=jmin, jmax
       rbbbs(j) = (dfit_R(jmax-j+2)-0.5*rwid)/delta_R+0.5*rwid
       zbbbs(j) = dfit_Z(1)/delta_Z
    end do

    jmin = jmax + 1
    jmax = jmax + nh/2
    do j=jmin, jmax
       rbbbs(j) = (dfit_R(1)-0.5*rwid)/delta_R+0.5*rwid
       zbbbs(j) = dfit_Z(j-jmin+2)/delta_Z
    end do

    do j=1,nbbbs-1
       if ((rbbbs(j) == rbbbs(j+1)) .and. (zbbbs(j) == zbbbs(j+1))) then
          write(*,*) 'duplicates: ',rbbbs(j),zbbbs(j)
       end if
    end do

! get r_boundary(theta)

    thetab = atan2 ((zbbbs-Z_mag), (rbbbs-R_mag))
    r_bound = sqrt( (rbbbs - R_mag)**2 + (zbbbs - Z_mag)**2 )

    call sort(thetab, r_bound, zbbbs, rbbbs)

! Allow for duplicated points near +- pi:

    if(thetab(1) == thetab(2)) then
       thetab(1) = thetab(1) + twopi
       call sort(thetab, r_bound, zbbbs, rbbbs)
    endif

    if(thetab(nbbbs-1) == thetab(nbbbs)) then
       thetab(nbbbs) = thetab(nbbbs) - twopi
       call sort(thetab, r_bound, zbbbs, rbbbs)
    endif

! It isn't likely that a duplicate point would exist near theta = 0,
! so I am not allowing this possibility for now.

    do i=1,nbbbs-1
       if(thetab(i) == thetab(i+1)) then
          write(*,*) 'Duplicates near theta = 0 not allowed.'
          write(*,*) i, i+1, ' Stopping.'
          call mp_abort('Duplicates near theta = 0 not allowed.')
       endif
    enddo
    deallocate (rbbbs, zbbbs)

    aminor=R_mag
    amin=aminor

    close(5)

    r_bound = r_bound / aminor

    R_mag = R_mag / aminor
    Z_mag = Z_mag / aminor
!    rleft = rleft / aminor
!    rwid = rwid / aminor
!    zhei = zhei / aminor
    dfit_R = dfit_R / aminor
    dfit_Z = dfit_Z / aminor

! should rmaj be R_mag? use R_mag for now.

    rmaj = R_mag
    B_T0 = abs(bcentr)
    B_T = B_T0
    psi_a = psi_a / (B_T0*aminor**2)
    psi_0 = psi_0 / (B_T0*aminor**2)
    psi_N = psi_a - psi_0

    p_0=pressure(1)

    fp = fp / (B_T0*aminor)

! MKS: beta = 2 mu_0 p / B**2

    beta = 8. * pi * pressure * 1.e-7 / B_T0**2
    beta_0 = beta(1)

    pressure = pressure / p_0

    dfit_psi = dfit_psi / (B_T0*aminor**2)

    dfit_dR = dfit_R(2) - dfit_R(1)
    dfit_dZ = dfit_Z(2) - dfit_Z(1)

    2020 format (5e16.8)

  end subroutine dfitin

  subroutine dfit_init

    real, dimension(nw, nh) :: eqth
    integer :: i, j

    do i = 1, nw
       do j = 1,nh
          if(dfit_Z(j) == Z_mag .and. dfit_R(i) == R_mag) then
             eqth(i,j) = 0.  ! value should not matter
          else
             eqth(i,j) = atan2( (dfit_Z(j)-Z_mag), (dfit_R(i)-R_mag))
          endif
       enddo
    enddo

    call derm(dfit_psi, dpm)
    call tderm(eqth, dtm)

  end subroutine dfit_init

  subroutine tderm(f, dfm)
    use constants, only: pi
    implicit none
    real, dimension(:,:), intent(in) :: f
    real, dimension(:,:,:), intent(out) :: dfm
    integer :: i, j

! DFIT grid is equally spaced in R, Z -- this routine uses that fact and
! is therefore not completely general.  It is fine for DFIT output.

    i=1
    dfm(i,:,1) = -0.5*(3*f(i,:)-4*f(i+1,:)+f(i+2,:))/dfit_dR

    i=nw
    dfm(i,:,1) = 0.5*(3*f(i,:)-4*f(i-1,:)+f(i-2,:))/dfit_dR

    j=1
    dfm(:,j,2) = -0.5*(3*f(:,j)-4*f(:,j+1)+f(:,j+2))/dfit_dZ

    j=nh
    dfm(:,j,2) = 0.5*(3*f(:,j)-4*f(:,j-1)+f(:,j-2))/dfit_dZ

    do i=2,nw-1
       dfm(i,:,1)=0.5*(f(i+1,:)-f(i-1,:))/dfit_dR
    enddo

    do j=2,nh-1
       do i = 1,nw
          if(f(i,j+1)-f(i,j-1) > pi) then
             dfm(i,j,2)=0.5*(f(i,j+1)-f(i,j-1)-2.*pi)/dfit_dZ
          else
             dfm(i,j,2)=0.5*(f(i,j+1)-f(i,j-1))/dfit_dZ
          endif
       enddo
    enddo

  end subroutine tderm

  subroutine derm(f, dfm)

    implicit none
    real, dimension(:,:), intent(in) :: f
    real, dimension(:,:,:), intent(out) :: dfm
    integer :: i, j

! DFIT grid is equally spaced in R, Z -- this routine uses that fact and
! is therefore not completely general.  It is fine for DFIT output.

    i=1
    dfm(i,:,1) = -0.5*(3*f(i,:)-4*f(i+1,:)+f(i+2,:))/dfit_dR

    i=nw
    dfm(i,:,1) = 0.5*(3*f(i,:)-4*f(i-1,:)+f(i-2,:))/dfit_dR

    j=1
    dfm(:,j,2) = -0.5*(3*f(:,j)-4*f(:,j+1)+f(:,j+2))/dfit_dZ

    j=nh
    dfm(:,j,2) = 0.5*(3*f(:,j)-4*f(:,j-1)+f(:,j-2))/dfit_dZ

    do i=2,nw-1
       dfm(i,:,1)=0.5*(f(i+1,:)-f(i-1,:))/dfit_dR
    enddo

    do j=2,nh-1
       dfm(:,j,2)=0.5*(f(:,j+1)-f(:,j-1))/dfit_dZ
    enddo

  end subroutine derm

  subroutine gradient(rgrid, theta, grad, char, rp, nth, ntm)
    use splines, only: inter_d_cspl
    implicit none
    real, dimension(-ntm:), intent(in) :: rgrid, theta
    real, dimension(-ntm:,:), intent(out) :: grad
    character(len=1), intent(in) :: char
    real, intent(in) :: rp
    integer, intent(in) :: nth, ntm
    real :: aa(1), daa(1), rpt(1)
    integer :: i

    grad = 0.
    do i=-nth, nth
       call eqitem(rgrid(i), theta(i), dpm(:,:,1), grad(i,1))
       call eqitem(rgrid(i), theta(i), dpm(:,:,2), grad(i,2))
    enddo

!     to get grad(pressure), multiply grad(psi) by dpressure/dpsi

    if(char == 'R') then
       rpt(1) = rp
       call inter_d_cspl(nw, psi_bar, pressure, 1, rpt, aa, daa)
       grad = grad*daa(1)*0.5* beta_0/psi_N
    endif

  end subroutine gradient

  subroutine bgradient(rgrid, theta, grad, char, rp, nth_used, ntm)
    use splines, only: inter_d_cspl
    implicit none

    real, dimension(-ntm:), intent(in) :: rgrid, theta
    real, dimension(-ntm:,:), intent(out) :: grad
    character(len=1), intent(in) :: char
    real, intent(in) :: rp
    integer, intent(in) :: nth_used, ntm
    real :: tmp(2), aa(1), daa(1), rpt(1)
    real, dimension(nw, nh, 2) ::  dbish
    integer :: i

    dbish(:,:,1) = sqrt(dpm(:,:,1)**2 + dpm(:,:,2)**2)
    dbish(:,:,2) = 0.

    if(char == 'T') then  ! the order of the next two statements matters
       dbish(:,:,2) = (dtm(:,:,2)*dpm(:,:,1)-dtm(:,:,1)*dpm(:,:,2))/dbish(:,:,1)
       dbish(:,:,1) = (dtm(:,:,1)*dpm(:,:,1)+dtm(:,:,2)*dpm(:,:,2))/dbish(:,:,1)
    endif

    do i=-nth_used,-1
       call eqitem(rgrid(i), theta(i), dbish(:,:,1), tmp(1))
       call eqitem(rgrid(i), theta(i), dbish(:,:,2), tmp(2))
       grad(i,1) = tmp(1)
       grad(i,2) = tmp(2)
    enddo

    do i=0,nth_used
       call eqitem(rgrid(i), theta(i), dbish(:,:,1), tmp(1))
       call eqitem(rgrid(i), theta(i), dbish(:,:,2), tmp(2))
       grad(i,1)=tmp(1)
       grad(i,2)=tmp(2)
    enddo

!     to get grad(pressure), multiply grad(psi) by dpressure/dpsi

    if(char == 'R') then
       rpt(1) = rp
       call inter_d_cspl(nw, psi_bar, pressure, 1, rpt, aa, daa)
       do i=-nth_used, nth_used
          grad(i,1)=grad(i,1)*daa(1) * 0.5*beta_0/psi_N
          grad(i,2)=grad(i,2)*daa(1) * 0.5*beta_0/psi_N
       enddo
    endif

  end subroutine bgradient

  subroutine eqitem(r, thetin, f, fstar)
    use mp, only: mp_abort
    implicit none
    real, intent (in) :: r, thetin, f(:,:)
    real, intent (out) :: fstar
    integer :: i, j, istar, jstar
    real :: st, dt, sr, dr
    real :: r_pos, z_pos

    r_pos = Rpos(r, thetin)
    z_pos = Zpos(r, thetin)

! find point on R mesh

    if(r_pos >= dfit_R(nw) .or. r_pos <= dfit_R(1)) then
       write(*,*) 'No evaluation of eqitem allowed outside'
       write(*,*) 'or on edge of R domain'
       write(*,*) r, thetin, dfit_R(nw), r_pos
       call mp_abort('No evaluation of eqitem allowed outside or on edge of R domain')
    endif

! ensure point is on Z mesh

    if(z_pos >= dfit_Z(nh) .or. z_pos <= dfit_Z(1)) then
       write(*,*) 'No evaluation of eqitem allowed outside'
       write(*,*) 'or on edge of Z domain'
       write(*,*) r, thetin, dfit_Z(1), dfit_Z(nh), z_pos
       call mp_abort('No evaluation of eqitem allowed outside or on edge of Z domain')
    endif

    istar=0
    do i=2,nw
       if(istar /= 0) cycle
       if(r_pos < dfit_R(i)) then
          dr = r_pos - dfit_R(i-1)
          sr = dfit_R(i) - r_pos
          istar=i-1
       endif
    enddo

! Now do Z direction

    jstar=0
    do j=1,nh
       if(jstar /= 0) cycle
       if(z_pos < dfit_Z(j)) then
          dt = z_pos - dfit_Z(j-1)
          st = dfit_Z(j) - z_pos
          jstar=j-1
       endif
    enddo

! use opposite area stencil to interpolate

    fstar=f(istar    , jstar    ) * sr * st &
         +f(istar + 1, jstar    ) * dr * st &
         +f(istar    , jstar + 1) * sr * dt &
         +f(istar + 1, jstar + 1) * dr * dt
    fstar = fstar &
         /abs(dfit_R(istar+1)-dfit_R(istar)) &
         /(dfit_Z(jstar+1)-dfit_Z(jstar))
!     write(*,*) i, dr, dt, sr, st
!     write(*,*) f(istar,jstar+1),f(istar+1,jstar+1)
!     write(*,*) f(istar,jstar),f(istar+1,jstar)
!     write(*,*) dfit_R(istar),dfit_R(istar+1)
!     write(*,*) dfit_Z(jstar),dfit_Z(jstar+1)

  end subroutine eqitem

  function Zpos (r, theta)
    real, intent (in) :: r, theta
    real :: Zpos

    Zpos = Z_mag + r * sin(theta)
  end function Zpos

  function Rpos (r, theta)
    real, intent (in) :: r, theta
    real :: Rpos

    Rpos = R_mag + r * cos(theta)
  end function Rpos

  function invR (r, theta)
    real, intent (in) :: r, theta
    real :: invR

    invR = 1/(R_mag + r*cos(theta))
  end function invR

  function initialize_psi (init)
    integer, intent(in) :: init
    integer :: initialize_psi

    init_psi = .false.
    if(init == 1) init_psi = .true.
    initialize_psi = 1
  end function initialize_psi

  function psi (r, theta)
    real, intent (in) :: r, theta
    real :: f, psi

    call eqitem(r, theta, dfit_psi, f)
    psi = f
  end function psi

  function initialize_btori (init)
    integer, intent(in) :: init
    integer :: initialize_btori

    init_btori = .false.
    if(init == 1) init_btori = .true.
    initialize_btori = 1
  end function initialize_btori

  function btori (pbar)
    use splines, only: new_spline, splint, spline
    implicit none
    real, intent(in) :: pbar
    real :: btori
    type (spline), save :: spl

    if(init_btori) call new_spline(nw, psi_bar, fp, spl)

    btori = splint(pbar, spl)
  end function btori

  function initialize_dbtori (init)
    integer, intent(in) :: init
    integer :: initialize_dbtori

    init_dbtori = .false.
    if(init == 1) init_dbtori = .true.
    initialize_dbtori = 1
  end function initialize_dbtori

  function dbtori (pbar)
    use splines, only: new_spline, dsplint, spline
    implicit none
    real, intent(in) :: pbar
    real :: dbtori
    type (spline), save :: spl

    if(init_dbtori) call new_spline(nw, psi_bar, fp, spl)

    dbtori = dsplint(pbar, spl)/psi_N
  end function dbtori

  function initialize_q (init)
    integer, intent(in) :: init
    integer :: initialize_q

    init_q = .false.
    if(init == 1) init_q = .true.
    initialize_q = 1
  end function initialize_q

  function qfun (pbar)
    use splines, only: new_spline, splint, spline
    implicit none
    real, intent(in) :: pbar
    real :: qfun
    type (spline), save :: spl

    if(init_q) then
       call new_spline(nw, psi_bar, qsf, spl)
    endif

    qfun = splint(pbar, spl)
  end function qfun

  function initialize_pressure (init)
    integer, intent(in) :: init
    integer :: initialize_pressure

    init_pressure = .false.
    if(init == 1) init_pressure = .true.
    initialize_pressure = 1
  end function initialize_pressure

  function pfun (pbar)
    use splines, only: new_spline, splint, spline
    implicit none
    real, intent(in) :: pbar
    real :: pfun
    type (spline), save :: spl

    if(init_pressure) call new_spline(nw, psi_bar, pressure, spl)

    pfun = splint(pbar, spl) * beta_0/2.
  end function pfun

  function initialize_dpressure (init)
    integer, intent(in) :: init
    integer :: initialize_dpressure

    init_dpressure = .false.
    if(init == 1) init_dpressure = .true.
    initialize_dpressure = 1
  end function initialize_dpressure

  function dpfun (pbar)
    use splines, only: new_spline, dsplint, spline
    implicit none
    real, intent(in) :: pbar
    real :: dpfun
    type (spline), save :: spl

    if(init_dpressure) call new_spline(nw, psi_bar, pressure, spl)

    dpfun = dsplint(pbar, spl)/psi_N * beta_0/2.
  end function dpfun

  function initialize_beta (init)
    integer, intent(in) :: init
    integer :: initialize_beta

    init_beta = .false.
    if(init == 1) init_beta = .true.
    initialize_beta = 1
  end function initialize_beta

  function betafun (pbar)
    use splines, only: new_spline, splint, spline
    implicit none
    real, intent(in) :: pbar
    real :: betafun
    type (spline), save :: spl

    if(init_beta) call new_spline(nw, psi_bar, beta, spl)

    betafun = splint(pbar, spl)
  end function betafun

  function initialize_rho (init)
    integer, intent(in) :: init
    integer :: initialize_rho

    init_rho = .false.
    if(init == 1) init_rho = .true.
    initialize_rho = 1
  end function initialize_rho

  function rhofun (pbar)
    use splines, only: new_spline, splint, spline
    implicit none
    real, intent(in) :: pbar
    real :: rhofun
    type (spline), save :: spl

    if(init_rho) call new_spline(nrho, psi_mid, rho_mid, spl)

    rhofun = splint(pbar, spl)
  end function rhofun

  function initialize_bound (init)
    integer, intent(in) :: init
    integer :: initialize_bound

    init_bound = .false.
    if(init == 1) init_bound = .true.
    initialize_bound = 1
  end function initialize_bound

  function bound(theta)
    use splines, only: new_spline, splint, spline
    implicit none
    real, intent(in) :: theta
    real :: bound
    type (spline), save :: spl

    if(init_bound) call new_spline(nbbbs, thetab, r_bound, spl)
    init_bound = .false.

    bound = splint(theta, spl)
  end function bound

  subroutine alloc_module_arrays(np, nw, nh, nws, nhs, nrho)
    integer, intent(in) :: np, nw, nh, nws, nhs, nrho

    allocate (rho_mid(nrho), psi_mid(nrho))
    allocate (psi_bar(np), fp(np), qsf(np), pressure(np), beta(np))
    allocate (dummy(nws), dfit_R(nw), dfit_Z(nh))
    allocate (spsi_bar(nws), sdfit_R(nws), sdfit_Z(nhs))
!  allocate (dfit_psi(nw, nh), sdfit_psi(nws, nhs))
    allocate (dfit_psi(nw, nh))
    allocate (dpm(nw, nh, 2), dtm(nw, nh, 2))

  end subroutine alloc_module_arrays

  subroutine dealloc_module_arrays
    implicit none
    if(allocated(rho_mid)) deallocate(rho_mid,psi_mid)
    if(allocated(psi_bar)) deallocate(psi_bar,fp,qsf,pressure,beta)
    if(allocated(dummy)) deallocate(dummy,dfit_R,dfit_Z)
    if(allocated(spsi_bar)) deallocate(spsi_bar,sdfit_R,sdfit_Z)
    if(allocated(dfit_psi)) deallocate(dfit_psi)
    if(allocated(dpm)) deallocate(dpm,dtm)
  end subroutine dealloc_module_arrays

  subroutine dfit_finish
    implicit none
    call dealloc_module_arrays
  end subroutine dfit_finish

!Could we use routines from the sorting module instead?
  subroutine sort(a, b, c, d)
    real, dimension(:), intent(in out) :: a, b, c, d
    real :: tmp
    integer :: i, j, jmax

    jmax = size(a)

    do j=1,jmax
       do i=1,jmax-j
          if(a(i+1) < a(i)) then
             tmp=a(i); a(i)=a(i+1); a(i+1)=tmp
             tmp=b(i); b(i)=b(i+1); b(i+1)=tmp
             tmp=c(i); c(i)=c(i+1); c(i+1)=tmp
             tmp=d(i); d(i)=d(i+1); d(i+1)=tmp
          endif
       enddo
    enddo
  end subroutine sort
end module deq


text:
Subroutine ppxs(lb1, lb2, srt, ppsig, spprho, ipp)
  Parameter (amp=0.14, pi=3.1415926)
  Save
  ppsig = 0.0
  spprho = 0.0
  ipp = 0
  If (srt<=0.3) Return
  q = sqrt((srt/2)**2-amp**2)
  esigma = 5.8*amp
  tsigma = 2.06*q
  erho = 0.77
  trho = 0.095*q*(q/amp/(1.+(q/erho)**2))**2
  esi = esigma - srt
  If (esi==0) Then
     d00 = pi/2.
     Goto 10
  End If
  d00 = atan(tsigma/2./esi)
10 erh = erho - srt
  If (erh==0.) Then
     d11 = pi/2.
     Goto 20
  End If
  d11 = atan(trho/2./erh)
20 d20 = -0.12*q/amp
  s0 = 8.*pi*sin(d00)**2/q**2
  s1 = 8*pi*3*sin(d11)**2/q**2
  s2 = 8*pi*5*sin(d20)**2/q**2
  s0 = s0*0.197**2*10.
  s1 = s1*0.197**2*10.
  s2 = s2*0.197**2*10.
  spprho = s1/2.
  If (lb1==5 .And. lb2==5) Then
     ipp = 1
     ppsig = s2
     Return
  End If
  If ((lb1==5 .And. lb2==4) .Or. (lb1==4 .And. lb2==5)) Then
     ipp = 2
     ppsig = s2/2. + s1/2.
     Return
  End If
  If ((lb1==5 .And. lb2==3) .Or. (lb1==3 .And. lb2==5)) Then
     ipp = 3
     ppsig = s2/6. + s1/2. + s0/3.
     Return
  End If
  If (lb1==4 .And. lb2==4) Then
     ipp = 4
     ppsig = 2*s2/3. + s0/3.
     Return
  End If
  If ((lb1==4 .And. lb2==3) .Or. (lb1==3 .And. lb2==4)) Then
     ipp = 5
     ppsig = s2/2. + s1/2.
     Return
  End If
  If (lb1==3 .And. lb2==3) Then
     ipp = 6
     ppsig = s2
  End If
  Return
End Subroutine ppxs

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
create ascii-table
	BL ,
	chars `'.-,_:^~!";r/\+()><|=?lci[]vtzjL7fx*s}{YTJ1unCyIFo2ewVh3ka%Z45SX$PmGApqbdUE&K69OHg#D8RQW0MB@N

here ascii-table - const MAX_CHARS ( number of chars in ascii-table )

-2.00 const X1
 0.60 const X2
-1.12 const Y1
 1.12 const Y2

0 var x
0 var y

: sqr ( n -- n^2 ) dup * ;

: mandel-at ( x0 y0 -- iter )
	0 x !  ( clear x & y )
	0 y !
	0 ( s: x0 y0 iter )
	begin
		( check termination )
		dup [ MAX_CHARS 1- lit ] =   x @ sqr   y @ sqr +   4 >   or  if
			nip nip exit
		then
		( step )
		x @ sqr  y @ sqr -   4 pick + ( x' = x*x - y*y + x0 )
		x @ y @ * 2 *   4 pick + ( y' = 2*x*y + y0 )
		y ! ( y = y' )
		x ! ( x = x' )
		1+ ( increment iteration )
	again ;

: mandel ( width height -- )
	dup 0 do
		over 0 do
			( map to range )
			over i swap 1- / [ X2 X1 - lit ] * X1 + ( x )
			over j swap 1- / [ Y2 Y1 - lit ] * Y1 + ( y )
			mandel-at
			ascii-table + @ emit
		loop
		cr
	loop 2drop ;

50 24 mandel

max_stars_repo_name: petrihakkinen/f

max_issues_repo_name: petrihakkinen/f

max_forks_repo_name: petrihakkinen/f


text:
! $UWHPSC/codes/fortran/multifile1/fullcode.f90

program demo
    print *, "In main program"
    call sub1()
    call sub2()
end program demo

subroutine sub1()
    print *, "In sub1"
end subroutine sub1

subroutine sub2()
    print *, "In sub2"
end subroutine sub2

max_stars_repo_name: philipwangdk/HPC

max_issues_repo_name: philipwangdk/HPC

max_forks_repo_name: philipwangdk/HPC


text:
!==============================================================================!
  subroutine User_Mod_End_Of_Correct_Velocity(flow, dt, ini)
!------------------------------------------------------------------------------!
!   This function is called at the end of Correct_Velocity function.           !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Field_Type) :: flow
  real             :: dt    ! time step
  integer          :: ini   ! iner itteration
!==============================================================================!

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows