text:
MODULE m_greensfTorgue

   USE m_types
   USE m_juDFT
   USE m_constants
   USE m_intgr
   USE m_gaunt
   USE m_xmlOutput

   IMPLICIT NONE

   CONTAINS

   SUBROUTINE greensfTorgue(greensFunction,sphhar,atoms,sym,noco,nococonv,input,f,g,flo,atomType,torgue,vTot)

      !--------------------------------------------------------------------------
      ! This Subroutine implements the formula:
      !   alpha     1                  ->    Ef        i ->    alpha     ->->
      !  J      = - -  Im Tr      int dx  int   dE    B (x) sigma    G  (x,x,E)
      !   i         pi      sigma           -infinity  xc             ii
      !
      ! For the evaluation of the torgue at site i (argument atomType)
      !--------------------------------------------------------------------------

      TYPE(t_greensf),        INTENT(IN)     :: greensFunction(:)
      TYPE(t_sphhar),         INTENT(IN)     :: sphhar
      TYPE(t_atoms),          INTENT(IN)     :: atoms
      TYPE(t_sym),            INTENT(IN)     :: sym
      TYPE(t_noco),           INTENT(IN)     :: noco
      TYPE(t_nococonv),       INTENT(IN)     :: nococonv
      TYPE(t_input),          INTENT(IN)     :: input
      REAL,                   INTENT(IN)     :: f(:,:,0:,:,:)
      REAL,                   INTENT(IN)     :: g(:,:,0:,:,:)
      REAL,                   INTENT(IN)     :: flo(:,:,:,:,:)
      INTEGER,                INTENT(IN)     :: atomType
      REAL,                   INTENT(INOUT)  :: torgue(:)
      TYPE(t_potden),         INTENT(IN)     :: vTot

      INTEGER :: na,nsym,nh,i_gf,l,lp,iContour
      INTEGER :: lh,mem,mu,m,mp,iz,ipm,lamda,jr,alpha
      COMPLEX :: phaseFactor
      REAL    :: realIntegral, imagIntegral
      COMPLEX :: sigma(2,2,3),torgue_cmplx(3),g_Spin(2,2)
      CHARACTER(LEN=20) :: attributes(5)

      REAL,    ALLOCATABLE :: bxc(:,:)
      COMPLEX, ALLOCATABLE :: g_ii(:,:),g_iiSpin(:,:,:,:)

      IF(input%jspins.NE.2) CALL juDFT_error("Torgue calculation only for magnetic systems", calledby="greensFunctionTorgue")
      IF(sym%nop>1) CALL juDFT_warn("Torgue calculation only without symmetries", calledby="greensFunctionTorgue")

      CALL timestart("Green's Function Torgue: init")
      !Get Bxc from the total potential (local frame)
      ALLOCATE(bxc(SIZE(vTot%mt,1),0:SIZE(vTot%mt,2)-1))
      bxc = vTot%mt(:,:,atomType,1) - vTot%mt(:,:,atomType,2)

      na=SUM(atoms%neq(:atomType-1))+1
      nsym = sym%ntypsy(na)
      nh = sphhar%nlh(nsym)

      ! sigma are the Pauli matrices
      sigma=cmplx_0
      sigma(1,2,1)=CMPLX(1.0,0.0)
      sigma(2,1,1)=CMPLX(1.0,0.0)
      sigma(1,2,2)=CMPLX(0.0,-1.0)
      sigma(2,1,2)=CMPLX(0.0,1.0)
      sigma(1,1,3)=CMPLX(1.0,0.0)
      sigma(2,2,3)=CMPLX(-1.0,0.0)

      CALL timestop("Green's Function Torgue: init")
      CALL timestart("Green's Function Torgue: Integration")
      torgue_cmplx = cmplx_0
      iContour = -1
      DO i_gf = 1, SIZE(greensFunction)

         IF(greensFunction(i_gf)%elem%atomType.NE.atomType.OR.&
            greensFunction(i_gf)%elem%atomTypep.NE.atomType) THEN
            CALL juDFT_error("Provided greensFunction for wrong atomType", calledby="greensFunctionTorgue")
         ENDIF

         l  = greensFunction(i_gf)%elem%l
         lp = greensFunction(i_gf)%elem%lp

         IF(iContour == -1) THEN
            iContour = greensFunction(i_gf)%elem%iContour
         ELSE IF(greensFunction(i_gf)%elem%iContour/=iContour) THEN
            CALL juDFT_error("Provided different energy contours", calledby="greensFunctionTorgue")
         ENDIF

!         !$OMP parallel default(none) &
!         !$OMP shared(sphhar,atoms,greensFunction,f,g,flo,sigma,bxc) &
!         !$OMP shared(nh,nsym,l,lp,i_gf,atomType,torgue_cmplx) &
!         !$OMP private(lh,m,lamda,mem,mu,mp,phaseFactor,ipm,iz,alpha,jr) &
!         !$OMP private(realIntegral,imagIntegral,g_ii,g_iiSpin,g_Spin)
         ALLOCATE(g_ii(atoms%jmtd,greensFunction(i_gf)%contour%nz),source=cmplx_0)
         ALLOCATE(g_iiSpin(2,2,atoms%jmtd,greensFunction(i_gf)%contour%nz),source=cmplx_0)
!         !$OMP do collapse(2) reduction(+:torgue_cmplx)
         DO lh = 0, nh
            DO m = -l, l
               lamda = sphhar%llh(lh,nsym)
               IF(MOD(lamda+l+lp,2) .NE. 0) CYCLE
               IF(lamda.GT.l+lp) CYCLE
               IF(lamda.LT.abs(l-lp)) CYCLE
               DO mem = 1,sphhar%nmem(lh,nsym)
                  mu = sphhar%mlh(mem,lh,nsym)
                  mp = m - mu
                  IF(ABS(mp).GT.lp) CYCLE
                  phaseFactor = (sphhar%clnu(mem,lh,nsym))*gaunt1(lp,lamda,l,mp,mu,m,atoms%lmaxd)
                  IF(ABS(phaseFactor).LT.1e-12) CYCLE !Naive approach just skip all elements with zero gaunt coefficient
                  DO ipm = 1, 2
                     CALL greensFunction(i_gf)%getRadialSpin(atoms,m,mp,ipm==2,f(:,:,:,:,atomType),g(:,:,:,:,atomType),&
                                                             flo(:,:,:,:,atomType),g_iiSpin)
                     DO iz = 1, SIZE(g_ii,2)
                        DO alpha = 1, 3 !(x,y,z)
                           DO jr = 1, atoms%jri(atomType)
                              IF(ipm==1) THEN
                                 g_Spin = matmul(sigma(:,:,alpha),g_iiSpin(:,:,jr,iz))
                              ELSE
                                 g_Spin = matmul(conjg(sigma(:,:,alpha)),g_iiSpin(:,:,jr,iz))
                              ENDIF
                              g_ii(jr,iz) = g_Spin(1,1) + g_Spin(2,2)
                           ENDDO
                           CALL intgr3(REAL(g_ii(:,iz)*bxc(:,lh)),atoms%rmsh(:,atomType),atoms%dx(atomType),atoms%jri(atomType),realIntegral)
                           CALL intgr3(AIMAG(g_ii(:,iz)*bxc(:,lh)),atoms%rmsh(:,atomType),atoms%dx(atomType),atoms%jri(atomType),imagIntegral)
                           torgue_cmplx(alpha) = torgue_cmplx(alpha) - 1/(2*ImagUnit*pi_const) * (-1)**(ipm-1) * (realIntegral+ImagUnit*imagIntegral) &
                                                * MERGE(phaseFactor*greensFunction(i_gf)%contour%de(iz),conjg(phaseFactor*greensFunction(i_gf)%contour%de(iz)),ipm.EQ.1)
                        ENDDO
                     ENDDO
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
!         !$OMP end do
         DEALLOCATE(g_ii,g_iiSpin)
!         !$OMP end parallel

      ENDDO
      torgue = REAL(torgue_cmplx)
      CALL timestop("Green's Function Torgue: Integration")

      WRITE(oUnit,'(A,I4,A,3f14.8,A)') '  atom: ', atomType, '   torgue: ', torgue * hartree_to_ev_const * 1000, ' meV'

      attributes = ''
      WRITE(attributes(1),'(i0)') atomType
      WRITE(attributes(2),'(f14.8)') torgue(1) * hartree_to_ev_const * 1000
      WRITE(attributes(3),'(f14.8)') torgue(2) * hartree_to_ev_const * 1000
      WRITE(attributes(4),'(f14.8)') torgue(3) * hartree_to_ev_const * 1000
      WRITE(attributes(5),'(a3)') 'meV'
      CALL writeXMLElementForm('torgue',['atomType','sigma_x ','sigma_y ','sigma_z ','unit    '],&
                               attributes,reshape([8,7,7,7,4,6,14,14,14,3],[5,2]))

   END SUBROUTINE greensfTorgue

   SUBROUTINE greensfSOTorgue(greensFunction,sphhar,atoms,sym,noco,nococonv,input,f,g,flo,atomType,torgue,vso)

      !--------------------------------------------------------------------------
      ! This Subroutine implements the formula:
      !   alpha     1                  ->    Ef        i ->    alpha     ->->
      !  J      = - -  Im Tr      int dx  int   dE    V (x) sigma    G  (x,x,E)
      !   i         pi      sigma           -infinity  SO             ii
      !
      ! For the evaluation of the torgue at site i (argument atomType) caused by Spin-orbit-coupling
      !--------------------------------------------------------------------------

      TYPE(t_greensf),        INTENT(IN)     :: greensFunction(:)
      TYPE(t_atoms),          INTENT(IN)     :: atoms
      TYPE(t_sphhar),         INTENT(IN)     :: sphhar
      TYPE(t_sym),            INTENT(IN)     :: sym
      TYPE(t_noco),           INTENT(IN)     :: noco
      TYPE(t_nococonv),       INTENT(IN)     :: nococonv
      TYPE(t_input),          INTENT(IN)     :: input
      REAL,                   INTENT(IN)     :: f(:,:,0:,:,:)
      REAL,                   INTENT(IN)     :: g(:,:,0:,:,:)
      REAL,                   INTENT(IN)     :: flo(:,:,:,:,:)
      INTEGER,                INTENT(IN)     :: atomType
      REAL,                   INTENT(INOUT)  :: torgue(:)
      REAL,                   INTENT(IN)     :: vso(:,0:)

      INTEGER :: jspin,na,nsym,nh,i_gf,l,lp,spin,iContour
      INTEGER :: lh,mems,mem,mh,m,mp,iz,ipm,lamda,jr,alpha
      COMPLEX :: phaseFactor
      REAL    :: realIntegral, imagIntegral
      COMPLEX :: sigma(2,2,3),chi(2,2),torgue_cmplx(3,2),g_Spin(2,2)
      CHARACTER(LEN=20) :: attributes(5)

      COMPLEX, ALLOCATABLE :: g_ii(:,:),g_iiSpin(:,:,:,:)

      TYPE(t_usdus) :: usdus

      IF(input%jspins.NE.2) CALL juDFT_error("Torgue calculation only for magnetic systems", calledby="greensFunctionTorgue")
      IF(sym%nop>1) CALL juDFT_warn("Torgue calculation only without symmetries", calledby="greensFunctionTorgue")

      na=SUM(atoms%neq(:atomType-1))+1

      ! sigma are the Pauli matrices
      sigma=cmplx_0
      sigma(1,2,1)=CMPLX(1.0,0.0)
      sigma(2,1,1)=CMPLX(1.0,0.0)
      sigma(1,2,2)=CMPLX(0.0,-1.0)
      sigma(2,1,2)=CMPLX(0.0,1.0)
      sigma(1,1,3)=CMPLX(1.0,0.0)
      sigma(2,2,3)=CMPLX(-1.0,0.0)

      CALL timestart("Green's Function SOTorgue: Integration")
      torgue_cmplx = cmplx_0
      iContour = -1
      ALLOCATE(g_ii(atoms%jmtd,greensFunction(1)%contour%nz),source=cmplx_0)
      DO i_gf = 1, SIZE(greensFunction)

         IF(greensFunction(i_gf)%elem%atomType.NE.atomType.OR.&
            greensFunction(i_gf)%elem%atomTypep.NE.atomType) THEN
            CALL juDFT_error("Provided greensFunction for wrong atomType", calledby="greensFunctionTorgue")
         ENDIF

         l  = greensFunction(i_gf)%elem%l
         lp = greensFunction(i_gf)%elem%lp

         IF(iContour == -1) THEN
            iContour = greensFunction(i_gf)%elem%iContour
         ELSE IF(greensFunction(i_gf)%elem%iContour/=iContour) THEN
            CALL juDFT_error("Provided different energy contours", calledby="greensFunctionTorgue")
         ENDIF

         DO lh = 0, SIZE(vso,2)-1
            DO mh = -lh,lh
               DO m = -l, l
                  DO mp = -lp, lp
                     phaseFactor = gaunt1(lp,lh,l,mp,mh,m,atoms%lmaxd)
                     IF(ABS(phaseFactor).LT.1e-12) CYCLE !Naive approach just skip all elements with zero gaunt coefficient
                     DO ipm = 1, 2
                        CALL greensFunction(i_gf)%getRadialSpin(atoms,m,mp,ipm==2,f(:,:,:,:,atomType),g(:,:,:,:,atomType),&
                                                                flo(:,:,:,:,atomType),g_iiSpin)
                        DO iz = 1, SIZE(g_ii,2)
                           DO alpha = 1, 3 !(x,y,z)
                              DO jr = 1, atoms%jri(atomType)
                                 IF(ipm==1) THEN
                                    g_Spin = matmul(sigma(:,:,alpha),g_iiSpin(:,:,jr,iz))
                                 ELSE
                                    g_Spin = matmul(conjg(sigma(:,:,alpha)),g_iiSpin(:,:,jr,iz))
                                 ENDIF
                                 g_ii(jr,iz) = g_Spin(1,1) + g_Spin(2,2)
                              ENDDO
                              CALL intgr3(REAL(g_ii(:,iz)*vso(:,lh)),atoms%rmsh(:,atomType),atoms%dx(atomType),atoms%jri(atomType),realIntegral)
                              CALL intgr3(AIMAG(g_ii(:,iz)*vso(:,lh)),atoms%rmsh(:,atomType),atoms%dx(atomType),atoms%jri(atomType),imagIntegral)
                              torgue_cmplx(alpha,ipm) = torgue_cmplx(alpha,ipm) - 1/(2*ImagUnit*pi_const) * (-1)**(ipm-1) * (realIntegral+ImagUnit*imagIntegral) &
                                                   * MERGE(phaseFactor*greensFunction(i_gf)%contour%de(iz),conjg(phaseFactor*greensFunction(i_gf)%contour%de(iz)),ipm.EQ.1)
                           ENDDO
                        ENDDO
                     ENDDO
                  ENDDO
               ENDDO
            ENDDO
         ENDDO

      ENDDO
      WRITE(*,*) torgue_cmplx(:,1) * hartree_to_ev_const * 1000
      WRITE(*,*) torgue_cmplx(:,2) * hartree_to_ev_const * 1000
      torgue = REAL(torgue_cmplx(:,1)+torgue_cmplx(:,2))
      CALL timestop("Green's Function SOTorgue: Integration")

      WRITE(oUnit,'(A,I4,A,3f14.8,A)') '  atom: ', atomType, '   torgue: ', torgue * hartree_to_ev_const * 1000, ' meV'

      attributes = ''
      WRITE(attributes(1),'(i0)') atomType
      WRITE(attributes(2),'(f14.8)') torgue(1) * hartree_to_ev_const * 1000
      WRITE(attributes(3),'(f14.8)') torgue(2) * hartree_to_ev_const * 1000
      WRITE(attributes(4),'(f14.8)') torgue(3) * hartree_to_ev_const * 1000
      WRITE(attributes(5),'(a3)') 'meV'
      CALL writeXMLElementForm('torgue',['atomType','sigma_x ','sigma_y ','sigma_z ','unit    '],&
                               attributes,reshape([8,7,7,7,4,6,14,14,14,3],[5,2]))

   END SUBROUTINE greensfSOTorgue

END MODULE m_greensfTorgue

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
!-----subroutine---------------------------------------------g.guidati--
!
     subroutine RHSINT(value,dval1,dval2,dval11,dval12,dval22)
!
!-----------------------------------------------------------------------
!
!     franco guidati    IAG
!                       Uni Stuttgart
!
!
!     scope             determine incident wave-field
!
!.......................................................................
!               declarations
!.......................................................................
      USE TINoiseGeneric
      USE TINoiseInput
      USE TINoiseDDD
      IMPLICIT                        NONE

   ! Local variables.
INTEGER (4)                  :: ivar,nvar,nok,nbad

REAL(DbKi)                   :: eps,h1, hmin
REAL(DbKi)                   :: phif (10)
EXTERNAL          :: RHSI,rkqs

COMPLEX(DbKi)                :: value, F2, arg1, arg2, ampli, dval1, dval2,dval11,dval12,dval22

      nvar    =  10
      eps     =  1.0D-5
      h1      =  0.1d0
      hmin    =  0.0d0
      s1      = tim(1)
      s2      = tim(nstr)

!.......................................................................
!     determine right-hand-side
!.......................................................................

      do ivar=1,nvar
        phif(ivar) = 0.0d0
      end do

      call odeint(phif,nvar,s1,s2,eps,h1,hmin,nok,nbad,RHSI,rkqs)
      value = dcmplx(phif(1),phif(2))
      dval1 = dcmplx(phif(3),phif(4))
      dval2 = dcmplx(phif(5),phif(6))
      dval11 = dcmplx(phif(7),phif(8))
      dval12 = dcmplx(phif(9),phif(10))
      dval22 = -dval11 - kwave2 * value

!.......................................................................
!               end of subroutine
!.......................................................................
      return
!***********************************************************************
      end

!***********************************************************************

      subroutine RHSI(s,y,dydx)
      USE TINoiseGeneric
      USE TINoiseInput
      USE TINoiseDDD
      USE TINoiseRHSin
      IMPLICIT                        NONE

   ! Local variables.
INTEGER (4)                  :: khi,klo

REAL(DbKi)                   :: s,y1,y2,d1y1,d1y2,d1,d2,d1s,d2s,r2
REAL(DbKi)                   :: d1sdr2,d2sdr2,d1dr2,d2dr2,arg
REAL(DbKi)                   :: y(10), dydx(10)

COMPLEX(DbKi)                :: HANK1, HANK0, HA0, HA1,n1amp,n2amp
COMPLEX(DbKi)                :: dipol,dummy,d1ss,d2ss
COMPLEX(DbKi)                :: gr,gr1,gr2,gr11,gr12,gr22,gr111,gr112,gr122,gr222
COMPLEX(DbKi)                :: dkern1,dkern2,dkern3,dkern4,dkern5

      call SPL_EX  (tim,nstr,s,khi,klo)
      call SPL_EX1 (tim,pstr1(1,ipath),d2pstr1(1,ipath),nstr,s,y1,d1y1,khi,klo)
      call SPL_EX1 (tim,pstr2(1,ipath),d2pstr2(1,ipath),nstr,s,y2,d1y2,khi,klo)
      call SPL_E1A (tim,dipole_strength(1,ipath),d2dipole_strength(1,ipath),nstr,s,dipol,dummy,khi,klo)

      d1 = y1-x1
      d2 = y2-x2
      d1s = d1**2
      d2s = d2**2
      r2 = d1s + d2s
      d1sdr2 = d1s / r2
      d2sdr2 = d2s / r2
      d1dr2 = d1/r2
      d2dr2 = d2/r2
      n1amp = - d1y2 * dipol * abb1
      n2amp =   d1y1 * dipol * abb1

      arg = r2 * kwave2
      HA1 = HANK1(arg)
      HA0 = HANK0(arg)
      gr = (0.0,-0.25) * HA0
      gr1   = d1*HA1
      gr2   = d2*HA1
      dkern1 = (n1amp * gr1 + n2amp * gr2 )

      dydx( 1) =  dble(dkern1)
      dydx( 2) =  dble(-imag*dkern1)

      if(lderiv) then
        gr11  = -(HA0*d1sdr2+(1.0d0-2.0d0*d1sdr2)*HA1)
        gr12  = -d1*d2/r2*(HA0-2.0d0*HA1)
        gr22  = -(HA0*d2sdr2+(1.0d0-2.0d0*d2sdr2)*HA1)
        gr111 = d1dr2*(-kwave2*d1s+8.0d0*d1sdr2-6.0d0)*HA1+d1dr2*(-4.0d0*d1sdr2+3.0d0)*HA0
        gr222 = d2dr2*(-kwave2*d2s+8.0d0*d2sdr2-6.0d0)*HA1+d2dr2*(-4.0d0*d2sdr2+3.0d0)*HA0
        gr112 = -kwave2 * gr2 - gr222
        gr122 = -kwave2 * gr1 - gr111

        dkern2 = (n1amp * gr11   + n2amp * gr12  )
        dkern3 = (n1amp * gr12   + n2amp * gr22  )
        dkern4 = (n1amp * gr111  + n2amp * gr112 )
        dkern5 = (n1amp * gr112  + n2amp * gr122 )

        dydx( 3) =  dble(dkern2)
        dydx( 4) =  dble(-imag*dkern2)
        dydx( 5) =  dble(dkern3)
        dydx( 6) =  dble(-imag*dkern3)
        dydx( 7) =  dble(dkern4)
        dydx( 8) =  dble(-imag*dkern4)
        dydx( 9) =  dble(dkern5)
        dydx(10) =  dble(-imag*dkern5)
      else
        dydx( 3) = 0.0
        dydx( 4) = 0.0
        dydx( 5) = 0.0
        dydx( 6) = 0.0
	dydx( 7) = 0.0
        dydx( 8) = 0.0
      	dydx( 9) = 0.0
      	dydx(10) = 0.0
      endif
      return
      end

!***********************************************************************

max_stars_repo_name: CarloSucameli/Openfast-noise

max_issues_repo_name: CarloSucameli/Openfast-noise

max_forks_repo_name: CarloSucameli/Openfast-noise


text:
C PR rtl-optimization/42246
C { dg-do compile { target powerpc*-*-* ia64-*-* x86_64-*-* } }
C { dg-options "-O2 -fselective-scheduling -fsel-sched-pipelining -fsel-sched-pipelining-outer-loops" }

      subroutine distance(x,clo)
      implicit real*8 (a-h,o-z)
      dimension x(2,6),x1(2,6),clo(6)
      do 60 i=1,2
        do 20 j=1,6
          x(i,j)=clo(j)
   20   continue
        do 40 iq=1,6
          x1(i,iq)=0.0d0
   40   continue
        do 50 j=1,6
          x(i,j)=x1(i,j)
   50   continue
   60 continue
      return
      end

ext: f

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!> 通过应用密度求和法更新密度的子程序。详见第 4 章中的论述。
!>   subroutine to calculate the density with sph summation algorithm.
!>   see equ.(4.35)
!>
!>     ntotal : number of particles                                  [in]
!>     hsml   : smoothing length                                     [in]
!>     mass   : particle masses                                      [in]
!>     niac   : number of interaction pairs                          [in]
!>     pair_i : list of first partner of interaction pair            [in]
!>     pair_j : list of second partner of interaction pair           [in]
!>     w      : kernel for all interaction pairs                     [in]
!>     itype   : type of particles                                   [in]
!>     x       : coordinates of all particles                        [in]
!>     rho    : density                                             [out]

subroutine sum_density(ntotal, hsml, mass, niac, pair_i, pair_j, w, itype, rho)

    use sph_kind, only: rk
    use parameter
    implicit none

    !> 在模拟中所使用的粒子总数
    !> number of particles in simulation
    integer, intent(in) :: ntotal
    !> 相互作用对的数目
    integer, intent(in) :: niac
    integer :: pair_i(max_interaction), pair_j(max_interaction), itype(maxn)
    real(rk) :: hsml(maxn), mass(maxn), w(max_interaction), rho(maxn)
    integer :: i, j, k, d
    real(rk) :: selfdens, hv(dim), r, wi(maxn)

    !     wi(maxn)---integration of the kernel itself

    do d = 1, dim
        hv(d) = 0._rk
    end do

    !     self density of each particle: wii (kernel for distance 0)
    !     and take contribution of particle itself:

    r = 0._rk

    !     firstly calculate the integration of the kernel over the space

    do i = 1, ntotal
        call kernel(r, hv, hsml(i), selfdens, hv)
        wi(i) = selfdens*mass(i)/rho(i)
    end do

    do k = 1, niac
        i = pair_i(k)
        j = pair_j(k)
        wi(i) = wi(i) + mass(j)/rho(j)*w(k)
        wi(j) = wi(j) + mass(i)/rho(i)*w(k)
    end do

    !     secondly calculate the rho integration over the space

    do i = 1, ntotal
        call kernel(r, hv, hsml(i), selfdens, hv)
        rho(i) = selfdens*mass(i)
    end do

    !     calculate sph sum for rho:
    do k = 1, niac
        i = pair_i(k)
        j = pair_j(k)
        rho(i) = rho(i) + mass(j)*w(k)
        rho(j) = rho(j) + mass(i)*w(k)
    end do

    !     thirdly, calculate the normalized rho, rho=sum(rho)/sum(w)

    if (nor_density) then
        do i = 1, ntotal
            rho(i) = rho(i)/wi(i)
        end do
    end if

end subroutine sum_density

!> 通过应用连续密度法更新密度的子程序。详见第 4 章中的论述。
!>     subroutine to calculate the density with sph continuity approach.
!>     see equ.(4.34)
!>
!>     ntotal : number of particles                                  [in]
!>     mass   : particle masses                                      [in]
!>     niac   : number of interaction pairs                          [in]
!>     pair_i : list of first partner of interaction pair            [in]
!>     pair_j : list of second partner of interaction pair           [in]
!>     dwdx   : derivation of kernel for all interaction pairs       [in]
!>     vx     : velocities of all particles                          [in]
!>     itype   : type of particles                                   [in]
!>     x      : coordinates of all particles                         [in]
!>     rho    : density                                              [in]
!>     drhodt : density change rate of each particle                [out]

subroutine con_density(ntotal, mass, niac, pair_i, pair_j, dwdx, vx, itype, x, rho, drhodt)

    use sph_kind, only: rk
    use parameter
    implicit none

    !> 在模拟中所使用的粒子总数
    !> number of particles in simulation
    integer, intent(in) :: ntotal
    !> 相互作用对的数目
    integer, intent(in) :: niac
    integer :: pair_i(max_interaction), pair_j(max_interaction), itype(maxn)
    real(rk) :: mass(maxn), dwdx(dim, max_interaction), vx(dim, maxn), x(dim, maxn), rho(maxn), drhodt(maxn)
    integer :: i, j, k, d
    real(rk) :: vcc, dvx(dim)

    do i = 1, ntotal
        drhodt(i) = 0._rk
    end do

    do k = 1, niac

        i = pair_i(k)
        j = pair_j(k)
        do d = 1, dim
            dvx(d) = vx(d, i) - vx(d, j)
        end do
        vcc = dvx(1)*dwdx(1, k)
        do d = 2, dim
            vcc = vcc + dvx(d)*dwdx(d, k)
        end do
        drhodt(i) = drhodt(i) + mass(j)*vcc
        drhodt(j) = drhodt(j) + mass(i)*vcc

    end do

end subroutine con_density

max_stars_repo_name: zoziha/SPH

max_issues_repo_name: zoziha/SPH

max_forks_repo_name: zoziha/SPH


text:
program main
  use structure_c
  use qsort_struc_m
  use rmsd_m
  implicit none
  type tag_list
    integer, dimension(:), allocatable :: i
  end type

  integer :: num_files, i, j, cand, k, cnt, curnum, tmpsize, tmpratio, curratio
  character(len=50) :: filename
  integer, dimension(:), allocatable :: tmp_array
  integer, dimension(0:19) :: ref_array, com_array
  type(structure), dimension(:), allocatable :: struc_rep
  type(tag_list), dimension(:), allocatable :: tl, tl_tmp
  double precision :: rmsd, rmsd_tmp, ediff, ediff_tmp, rmsd_lim = .5d0, ediff_lim = 5.d0, e_ub = 0.d0, anceps = 1.d0
  double precision, dimension(:), allocatable :: maxrmsd, maxediff, maxrmsd_tmp, maxediff_tmp
  double precision :: econv = 627.509469d0

  open(123,file='tmp.txt',status='old')
  read(123,*) num_files
  read(123,*) rmsd_lim, ediff_lim, econv
  read(123,*) anceps

  write(6,'(A)') 'Reading structures ...'

  allocate(struc_rep(num_files))

  read(123,*) curnum, filename
  call struc_read_file(struc_rep(1),trim(filename),'xyz')
  call struc_set(struc_rep(1),t=curnum)
  call struc_ctrans(struc_rep(1))
  do j = 0, 19
    ref_array(j) = struc_get_natoms(struc_rep(1),j)
  end do

  do i = 2, num_files
    read(123,*) curnum, filename
    call struc_read_file(struc_rep(i),trim(filename),'xyz')
    call struc_set(struc_rep(i),t=curnum)

    do j = 0, 19
      com_array(j) = struc_get_natoms(struc_rep(i),j)
    end do

    if (any(com_array .ne. ref_array)) then
      write(6,'(A)') 'Structure '//trim(filename)//' does not match the reference (atm fatal)'
      stop
    end if
    call struc_ctrans(struc_rep(i))
  end do

  close(123)

  call qsort(struc_rep)
  call rmsd_set(rmsd_lim,.true.,anceps)
  call rmsd_init(struc_rep(1))
  e_ub = struc_get_energy(struc_rep(1))

  open(95,file='rmsd.txt')

  allocate(tl(1),maxrmsd(1),maxediff(1))
  allocate(tl(1)%i(1))
  tl(1)%i(1) = 1
  maxrmsd(1) = 0.d0
  maxediff(1) = 0.d0

  write(6,'(A)') 'Starting analysis ...'

  tmpratio=0
  write(6,'(A)') '---|---|---|---|---|---|---|---|---|---|'

  do i=2,num_files
    rmsd = 1111.d0
    cand = 0
    call rmsd_set_anc_s1(struc_rep(i))
    !~write(*,*)
	curratio = i*40/num_files
	call progressbar_add(tmpratio,curratio)
	tmpratio = curratio
    do j=1,size(tl)
      ediff_tmp = abs(struc_get_energy(struc_rep(i)) - struc_get_energy(struc_rep(tl(j)%i(1))))*econv
      if (ediff_tmp .gt. ediff_lim) cycle

      call rmsd_set_anc_s2(struc_rep(tl(j)%i(1)))
      call rmsd_calc(rmsd_tmp,cnt)

      write(95,'(A,I4.1,A,I4.1,A,F10.4,A,F10.4,A,I4)') 'compare ',struc_rep(i)%tag,' and ',struc_rep(tl(j)%i(1))%tag,&
      &'; RMSD [A]: ',rmsd_tmp,';    Delta E [kcal/mol]: ',(struc_rep(tl(j)%i(1))%energy-struc_rep(i)%energy)*econv,&
      &';  Match count:', cnt

      if (rmsd_tmp .ge. rmsd) cycle
      cand = j
      rmsd = rmsd_tmp
      ediff = ediff_tmp
    end do
    if (rmsd .le. rmsd_lim) then
      !~write(*,'(I4.1,A,I4.1)') i, ' is higher in energy than ', tl(cand)%i(1)
      allocate(tmp_array(size(tl(cand)%i)))
      tmp_array = tl(cand)%i
      deallocate(tl(cand)%i)
      allocate(tl(cand)%i(size(tmp_array)+1))
      tl(cand)%i(:size(tmp_array)) = tmp_array
      tl(cand)%i(size(tmp_array)+1) = i
      deallocate(tmp_array)

      maxrmsd(cand) = max(maxrmsd(cand),rmsd)
      maxediff(cand) = max(maxediff(cand),ediff)
    else
      !~write(*,'(I4.1,A)') i, ' is a new structure'
      tmpsize = size(tl)
      allocate(tl_tmp(tmpsize),maxediff_tmp(tmpsize),maxrmsd_tmp(tmpsize))
      maxediff_tmp = maxediff
      maxrmsd_tmp = maxrmsd

      do k=1,tmpsize
        allocate(tl_tmp(k)%i(size(tl(k)%i)))
        tl_tmp(k)%i = tl(k)%i
      end do

      deallocate(tl,maxediff,maxrmsd)
      allocate(tl(tmpsize+1),maxediff(tmpsize+1),maxrmsd(tmpsize+1))
      maxediff(:tmpsize) = maxediff_tmp
      maxrmsd(:tmpsize) = maxrmsd_tmp

      do k=1,tmpsize
        allocate(tl(k)%i(size(tl_tmp(k)%i)))
        tl(k)%i = tl_tmp(k)%i
      end do

      allocate(tl(tmpsize+1)%i(1))
      tl(tmpsize+1)%i(1) = i
      maxediff(tmpsize+1) = 0.d0
      maxrmsd(tmpsize+1) = 0.d0

      deallocate(tl_tmp,maxediff_tmp,maxrmsd_tmp)
    end if
  end do
  close(95)
  write(6,'(A)')

  open(111,file='out.txt',status='replace')
  !~write(111,*) '#energy,', 'occurrence,', 'matching structures with increasing energy'
  do i=1,size(tl)
    write(111,'(F16.5,I4,F10.5,F10.5,9999I5)') (struc_rep(tl(i)%i(1))%energy-e_ub)*econv, size(tl(i)%i), maxrmsd(i), maxediff(i),&
    & struc_rep(tl(i)%i(:))%tag
  end do
  close(111)
  deallocate(struc_rep)
  stop
666 write(*,*) 'bad input'

contains
  subroutine progressbar_add(old,current)
    integer :: current, old, loopcnt

	do loopcnt=old+1,current
	  write(6,'(A)',advance='no') '*'
	end do

  end subroutine

end program

max_stars_repo_name: frostbaum/RMSDysator

max_issues_repo_name: frostbaum/RMSDysator

max_forks_repo_name: frostbaum/RMSDysator


text:
The Mac Paramedic is now called White Wave Digital

The Mac Paramedic serves both homes and businesses.  They offer onsite service and support (including new setups, purchase advice, printer problems, computer repair, etc.), primarily for Apple Macs and Linux Unix machines, and reluctantly for Windows machines.  They also offer networking services, such as networking setup and troubleshooting.

The owner is Certified by Apple Computer as a support technician and technical coordinator.

  See the Computers and Electronics entry for businesses that sell computers, or the Computer Repair entry for help when they break.

20090624 22:04:11 nbsp Very professional work done to get my macbook back up and running, not to mention in the span of about 2 days.  He points out any problems in your presence associated with the computer in question and gives you an idea of what he will do, then gets to work.  Upfront and honest, no bull, especially with logistics.  Pricing is very reasonable and he works efficiently.

I would easily recommend going to the Mac Paramedic over an Apple Store due to his availability and extremely flexible work hours.
In short, this IS the person to take your mac to when you have problems or perhaps questions.

Thanks Larry!! Users/ecj

20091008 11:06:23 nbsp Ive met Larry, and I think he is awesome. I would send any problem with an Apple I couldnt fix myself to him without hesitation. Users/NickSchmalenberger

20091011 19:49:10 nbsp Larry has been helping me setup, upgrade, troubleshoot, connect to DSL, fix my MacIntosh computers since 2001.  He promptly responds to emails/phone calls and is usually available within 2448 hours to work on the computer.  He is EXTREMELY competent and wellversed in his computer knowledge.  He is my one stop shop answer to all MacIntosh related matters!  Its great to have a good doctor . . . a good auto mechanic . . . and a good computer technician!   Users/cregross

20091023 20:40:02 nbsp Efficient, friendly, reasonably priced, and best of all, effective. It really was a pleasure to work with Larry, who returned my call right away, met me at home the next morning despite a busy schedule, and had my Macbook up and running within 24 hours. And the paramedic title is completely apt; he assessed my computers vitals with the mark of a pro. If I have any problems again with my Mac, my consolation will be working with Larry again! Users/cwbeebe

20100903 01:08:15 nbsp Mac Paramedic, is the place to go for apple assistance in Davis.  Larry understands all the new stuff coming out, and can explain how they work, not just with apple also with PC as well.  As more and more of us pc users purchase ipods and ipads it is complex how they connect, and when problems happen it is wonderful to be able to get assistance from a respected business in Davis. I appreciate all your help with my ipod/ pc interactions, I would be lost without your assistance.  jamie    Users/Jamie523

20101117 10:48:32 nbsp Larry Dietrich has helped me several times with my Macbook and has been able to do this very soon after I called for help. He can fix anything! and he manages to do it without making you feel stupid. This is a big plus. Larry is extremely honest and helpful. I would never go anywhere else with my Mac problems. Thank you Larry for helping me enjoy and learn about my computer! Users/AnneSyer

20101217 13:33:31 nbsp I would highly recommend Larry for any of your Mac repair needs. Right from the getgo, when I spoke to him on the phone he was very accommodating, offering to either make a house call or for me to come to his workshop. He was able to calmly diagnose the problem with my laptop, and his pricing is very reasonable.

The nearest Mac store is in Arden Fair, Sacramento and its usually overrun with people. Larry was able to see me in Davis the same day I contacted him, and so I avoided the inconvenience of the commute to Sac and waiting in a crowded store. I will definitely be calling him the next time I have a problem with my Mac, and you should do the same! Users/SoniaFerrandiz

20110303 13:41:53 nbsp Larry came over to our house and then found and fixed our Mac problem in about 2 minutes.  I highly recommend him! Users/TeresaKing

20110405 14:17:09 nbsp Recently my Mac Workbook started clicking at start up and showing only a gray screen with a folder and question mark blinking in the middle.  Uh oh.  I consulted Davis Wiki and found the Mac Paramedic.  I called him late on a Saturday morning.  He made an appointment to come to my house and diagnose the problem for the middle of that afternoon.  The news was bad  new hard drive needed.  He took the computer, replaced the hard drive, and recovered my data, and had the computer back to me by that same evening.  Wow!  Great service, reasonable rates.  He has a lot of parts at his shop so in many cases you can get even complex jobs back very quickly.  Most highly recommended. Users/Dora5

20110504 09:20:13 nbsp I highly recommend The Mac Paramedic!
After accidentally downloading a fake antivirus software on my MacBook Pro I called Larry. He had heard of the virus, and was able to guide me through each step to get rid of it right over the phone. I expected to spend the day dealing with this problem but Larry was able to help me in just a few minutes without even leaving my home.
Thank you so much Larry!

 Users/lisahicklin

20130109 13:32:14 nbsp I highly recommend Larry. He is active on the listserv for the Davis Mac Users Group. Last weekend he replied within two hours to my voicemail, and he solved the problem with my iMac over the phone in under 15 minutes. I would use his services again. Users/SylviaWright

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! 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.

!>    @brief return normalisation-vector for initial normalisation of the linear-system
!>    @param[in] mw right side, solution, diagonal (system matrix)
!>    @param[in] mbc_mask boundary mask
!>    @param[in] N number of elements, vector length
!>    @param[in] md diagonal of the system matrx
!>    @param[in] mx currentr solution vector (used for step by step updates)
!>    @param[out] dnormalise updated (new) normalisation vector
!>    @param[in] soddmx second offset (dominant) diagonal*mx
      SUBROUTINE get_norm2(n,mw,mbc_mask,md,mx,dnormalise,soddmx)
        IMPLICIT NONE
        INTEGER n, i
        DOUBLE PRECISION dnormalise(n), dcrit_nrm
!     right side, solution, diagonal (system matrix)
        DOUBLE PRECISION mw(n), mx(n), md(n), soddmx(n)
!     boundary mask
        CHARACTER mbc_mask(n)
        INTRINSIC dsqrt, dabs, dble
        LOGICAL test_null
        EXTERNAL test_null

!     correction for 2-Norm
        dcrit_nrm = dsqrt(dble(n))

        DO i = 1, n
          dnormalise(i) = 1.D0
!        add, since it is no boundary condition element
          IF (mbc_mask(i)=='+') THEN
!           important: diagonal dominance
            IF ( .NOT. test_null(mx(i))) THEN
              dnormalise(i) = 1.D0 /(max(dabs(md(i)*mx(i)), soddmx(i))*dcrit_nrm)
            ELSE
              dnormalise(i) = 1.D0 /(max(dabs(md(i)), dabs(mw(i)), soddmx(i))*dcrit_nrm)
            END IF
          END IF
        END DO
!
        RETURN
      END


!
! SDF (Self-Describing Format) Fortran Library
! Copyright (c) 2013-2016, SDF Development Team
!
! Distributed under the terms of the BSD 3-clause License.
! See the LICENSE file for details.
!

MODULE sdf_input_station

  USE sdf_input_station_ru
  USE sdf_input_station_r4
  USE sdf_input_station_r8

  IMPLICIT NONE

  INTERFACE sdf_read_station_array
    MODULE PROCEDURE &
        read_station_array_r4_r4, &
        read_station_array_r8_r4, &
        read_station_array_r4_r8, &
        read_station_array_r8_r8
  END INTERFACE sdf_read_station_array

END MODULE sdf_input_station


text:
      SUBROUTINE BAL2T4
C
C CHANGE RECORD
C  SUBROUTINE ADDED FOR 2 TIME-LEVEL BALANCES INCLUDING SED,SND,TOX
C **  SUBROUTINES CALBAL CALCULATE GLOBAL VOLUME, MASS, MOMENTUM,
C **  AND ENERGY BALANCES
C
      USE GLOBAL
	IMPLICIT NONE
	INTEGER::L,LN,K
	REAL::DUTMP,DVTMP
      IF(ISDYNSTP.EQ.0)THEN
        DELT=DT
      ELSE
        DELT=DTDYN
      END IF
C
C **  CALCULATE MOMENTUM AND ENERGY DISSIPATION
C
      DO L=2,LA
        LN=LNC(L)
        UUEOUT=UUEOUT+0.5*DELT*SPB(L)*DXYP(L)*(U(L,1)*TBX(L)
     &      +U(LEAST(L),1)*TBX(LEAST(L))-U(L,KC)*TSX(L)-U(LEAST(L),KC)*TSX(LEAST(L)))
        VVEOUT=VVEOUT+0.5*DELT*SPB(L)*DXYP(L)*(V(L,1)*TBY(L)
     &      +V(LN,1)*TBX(LN)-V(L,KC)*TSY(L)-V(LN,KC)*TSX(LN))
      ENDDO
      DO K=1,KS
        DO L=2,LA
          LN=LNC(L)
          DUTMP=0.5*( U(L,K+1)+U(LEAST(L),K+1)-U(L,K)-U(LEAST(L),K) )
          DVTMP=0.5*( V(L,K+1)+V(LN,K+1)-V(L,K)-V(LN,K) )
          UUEOUT=UUEOUT+DELT*SPB(L)*2.0*DXYP(L)*AV(L,K)
     &        *( DUTMP*DUTMP )/(DZC(K+1)+DZC(K))
          VVEOUT=VVEOUT+DELT*SPB(L)*2.0*DXYP(L)*AV(L,K)
     &        *( DVTMP*DVTMP )/(DZC(K+1)+DZC(K))
          BBEOUT=BBEOUT+DELT*SCB(L)*DXYP(L)*HP(L)
     &        *GP*AB(L,K)*(B(L,K+1)-B(L,K))
        ENDDO
      ENDDO
      RETURN
      END


text:
!! https://github.com/flang-compiler/flang/issues/713
!!
!! BAD SOURCED ALLOCATION
!!
!! The following example allocates an unlimited polymorphic array with a SOURCE=
!! argument that is an unlimited polymorphic array whose dynamic type is CHARACTER.
!! With the March 2019 binary release the value of the allocated array is not
!! correctly defined, and the code exits at the STOP 3 statement. (The value is
!! the character 'f' followed by 4 null characters.)
!!

character(:), allocatable :: ref(:,:)
class(*), allocatable :: b(:,:)
ref = reshape(['foo','bar'], shape=[1,2])
call foo(ref, b)
contains
  subroutine foo(a, b)
    class(*), intent(in)  :: a(:,:)
    class(*), allocatable :: b(:,:)
    allocate(b(size(a,1),size(a,2)), source=a) ! <== VALUE OF B NOT CORRECTLY DEFINED
    !allocate(b, source=a) ! unspecified bounds allowed by F2008
    select type (b)
    type is (character(*))
      if (len(b) /= len(ref)) stop 1
      if (any(shape(b) /= shape(ref))) stop 2
      if (any(b /= ref)) stop 3
    class default
      stop 4
    end select
  end subroutine
end

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


        !COMPILER-GENERATED INTERFACE MODULE: Thu Oct 25 10:21:43 2018
        MODULE VERT_ANOM__genmod
          INTERFACE
            FUNCTION VERT_ANOM(XX,YY,ZZZ,IPS)
              REAL(KIND=4) :: XX
              REAL(KIND=4) :: YY
              REAL(KIND=4) :: ZZZ
              INTEGER(KIND=4) :: IPS
              REAL(KIND=4) :: VERT_ANOM
            END FUNCTION VERT_ANOM
          END INTERFACE
        END MODULE VERT_ANOM__genmod


text:
      program main
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C                                                                      %
C Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
C Junior University.  All rights reserved.                             %
C                                                                      %
C The programs in GSLIB are distributed in the hope that they will be  %
C useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
C responsibility to anyone for the consequences of using them or for   %
C whether they serve any particular purpose or work at all, unless he  %
C says so in writing.  Everyone is granted permission to copy, modify  %
C and redistribute the programs in GSLIB, but only under the condition %
C that this notice and the above copyright notice remain intact.       %
C                                                                      %
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c-----------------------------------------------------------------------
c
c                  Normal/Lognormal Probability Plot
c                  *********************************
c
c Displays a set of data values on a probability plot with either an
c arithmetic or logarithmic scaling.
c
c INPUT/OUTPUT Parameters:
c
c   datafl         the data file
c   ivr,iwt        columns for the variable and the weight
c   tmin,tmax      trimming limits
c   outfl          output file for PostScript probability plot
c   npts           number of points to plot (negative means plot all)
c   ilog           1=logarithmic scale, 0=arithmetic
c   min, max, inc  min value, max value and increment for labeling
c   title          the title
c
c
c
c PROGRAM NOTES:
c
c 1. The program is executed with no command line arguments.  The user
c    will be prompted for the name of a parameter file.  The parameter
c    file is described in the documentation (see example scatplt.par)
c
c 2. The logarithmic scaling is base 10
c
c 3. All acceptable values outside the plotting limits (including
c    zero or negative values when a logarithmic scale is used) are
c    used to establish the cdf but are not shown.
c
c
c
c AUTHOR: Clayton V. Deutsch                             DATE: 1989-1999
c-----------------------------------------------------------------------
      use       msflib
      parameter(MV=500, NPVAL = 21, ZERO=0.0, BIGNUM=1.0e21,
     +          EPSLON=1.0e-21, VERSION=2.905)

      real      var(MV),xx(MV),yy(MV),pval(NPVAL)
      integer   test
      real*8    pdp
      character datafl*512,outfl*512,varlab*24,title*40,str*512
      logical   testfl,plotall
c
c Dynamic allocation of arrays that depend on the number of data.
c
      real, allocatable :: ar1(:), ar2(:)

      data lin/1/,lpsout/1/,pscl/0.24/,pxmin/0.0/,pxmax/288.0/
     +     pymin/0.0/,pymax/216.0/,xmin/-10.0/,xmax/60.0/,
     +     ymin/-10.0/,ymax/60.0/,hpxmin/1.0/,hpxmax/59.5/,
     +     hpymin/0.0/,hpymax/58.0/

      common /psdata/ lpsout,pscl,pxmin,pxmax,pymin,pymax,xmin,
     +                xmax,ymin,ymax
c
c The probability scaling and labelling:
c
      data pval/0.0001,0.001,0.002,0.01,0.02,0.05,0.10,0.20,0.30,
     +          0.40,0.50,0.60,0.70,0.80,0.90,0.95,0.98,0.99,0.998,
     +          0.999,0.9999/
c
c
c Note VERSION number:
c
      write(*,9999) VERSION
 9999 format(/' PROBPLT Version: ',f5.3/)
c
c Get the name of the parameter file - try the default name if no input:
c
      do i=1,512
            str(i:i) = ' '
      end do
      call getarg(1,str)
      if(str(1:1).eq.' ')then
            write(*,*) 'Which parameter file do you want to use?'
            read (*,'(a)') str
      end if
      if(str(1:1).eq.' ') str(1:20) = 'probplt.par         '
      inquire(file=str,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR - the parameter file does not exist,'
            write(*,*) '        check for the file and try again  '
            write(*,*)
            if(str(1:20).eq.'probplt.par         ') then
                  write(*,*) '        creating a blank parameter file'
                  call makepar
                  write(*,*)
            end if
            stop
      endif
      open(lin,file=str,status='OLD')
c
c Find Start of Parameters:
c
 1    read(lin,'(a4)',end=97) str(1:4)
      if(str(1:4).ne.'STAR') go to 1
c
c Read Input Parameters:
c
      read(lin,'(a512)',err=97) datafl
      call chknam(datafl,512)
      write(*,*) ' data file = ',datafl(1:40)

      read(lin,*,err=97) ivr,iwt
      write(*,*) ' columns = ',ivr,iwt

      read(lin,*,err=97) tmin,tmax
      write(*,*) ' tmin and tmax = ',tmin,tmax

      read(lin,'(a512)',err=97) outfl
      call chknam(outfl,512)
      write(*,*) ' outfl = ',outfl(1:40)

      read(lin,*,err=97) npts
      write(*,*) ' number of points to plot = ',npts

      read(lin,*,err=97) ilog
      write(*,*) ' logarithm scale option = ',ilog

      read(lin,*,err=97) pmin,pmax,pinc
      write(*,*) ' min, max, increment = ',pmin,pmax,pinc

      read(lin,'(a40)',err=97) title
      call chktitle(title,40)
      write(*,*) ' title = ',title
      write(*,*)

      close(lin)
c
c Check to make sure the data file exists, then either read in the
c data or write an error message and stop:
c
      inquire(file=datafl,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR - the data file does not exist,'
            write(*,*) '        check for the file and try again  '
            stop
      endif
c
c The data file exists so open the file and read in the header
c information. Initialize the storage that will be used to summarize
c the data found in the file:
c
      open(lin,file=datafl,status='OLD')
      read(lin,*,err=99)
      read(lin,*,err=99)     nvari
      do i=1,nvari
            read(lin,*)
      end do
      maxdat = 0
 20   read(lin,*,end=40,err=99)(var(j),j=1,nvari)
      if(var(ivr).lt.tmin.or.var(ivr).ge.tmax)go to 20
      maxdat = maxdat + 1
      go to 20
 40   continue
c
c Allocate the needed memory.
c
      allocate (ar1(maxdat),stat = test)
      if (test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
c
      allocate (ar2(maxdat),stat = test)
      if (test.ne.0) then
            write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
            stop
      end if
c
c Now read the data for real.
c
      rewind(lin)
      read(lin,'(a)',err=99) str
      read(lin,*,err=99)     nvari
      do i=1,nvari
            read(lin,'(a24)',err=99) str(1:24)
            if(i.eq.ivr) varlab = str(1:24)
      end do
      if(ivr.gt.nvari) then
            write(*,*) ' ivr is greater than the number in data file!'
            stop
      end if
c
c Read as much data as possible:
c
      nd = 0
      nt = 0
      vrmin =  BIGNUM
      vrmax = -BIGNUM
 7    read(lin,*,end=8,err=99) (var(j),j=1,nvari)
c
c Trim this data?
c
      if(var(ivr).lt.tmin.or.var(ivr).ge.tmax) then
            nt = nt + 1
            go to 7
      endif
      if(iwt.ge.1) then
            if(var(iwt).le.EPSLON) then
                  nt = nt + 1
                  go to 7
            endif
      endif
c
c Accept this data:
c
      nd = nd + 1
      ar1(nd) = var(ivr)
      vrmin = min(ar1(nd),vrmin)
      vrmax = max(ar1(nd),vrmax)
      if(iwt.ge.1) then
            ar2(nd) = var(iwt)
      else
            ar2(nd) = 1.0
      endif
c
c Go back for another data:
c
      go to 7
 8    close(lin)
      if(nd.lt.1)then
            write(*,*)'Not reading any data.'
            stop
      endif
c
c Get mean and total weight:
c
      xtwt = 0.0
      xmen = 0.0
      do i=1,nd
            xmen = xmen + ar1(i)*ar2(i)
            xtwt = xtwt + ar2(i)
      end do
      if(xtwt.lt.EPSLON) stop 'Cumulative Probability too LOW'
      xtwti = 1.0  / xtwt
      xmen  = xmen * xtwti
c
c Get the variance:
c
      xvar = 0.0
      do i=1,nd
            xvar = xvar + (ar1(i)-xmen) * (ar1(i)-xmen) * ar2(i)
      end do
      xvar  = xvar * xtwti
c
c Sort the Data in Ascending Order:
c
      call sortem(1,nd,ar1,1,ar2,c,d,e,f,g,h)
c
c Get cumulative probability and normalize:
c
      oldcp = 0.0
      cp    = 0.0
      do i=1,nd
            cp     = cp + ar2(i)*xtwti
            ar2(i) = 0.5*(cp+oldcp)
            oldcp  = cp
      end do
      call locate(ar2,nd,1,nd,0.50,i)
      if(i.eq.0) then
            xmed = ar1(1)
      else if(i.eq.nd) then
            xmed = ar1(nd)
      else
            xmed = ar1(i) +      (ar1(i+1)-ar1(i)) *
     +                      (0.50-ar2(i))/(ar2(i+1)-ar2(i))
      endif
c
c Write Some of the Statistics to the screen:
c
      write(*,900) nd,xmen,xmed,sqrt(max(xvar,0.0)),vrmin,vrmax,nt
 900  format(/' There are ',i8,' data with:',/,
     +        '   mean value          = ',f12.5,/,
     +        '   median value        = ',f12.5,/,
     +        '   standard deviation  = ',f12.5,/,
     +        '   minimum and maximum = ',2f12.5,/,
     +        '   number trimmed      = ',i8,/)
c
c Plot the probability plot:
c
      open(lpsout,file=outfl,status='UNKNOWN')
c
c Add a header:
c
      write(lpsout,998) title(1:20)
 998  format('%!PS                                 %    Remove     ',
     +    /, '90 234 translate 1.5 1.5 scale       %  these lines  ',
     +    /, '                                     % for EPSF file ',
     +    /, '%!PS-Adobe-3.0 EPSF-3.0',
     +    /, '%%BoundingBox: 0 0 288 216',
     +    /, '%%Creator: GSLIB',
     +    /, '%%Title:   ',a20,
     +    /, '%%CreationDate: ',
     +    /, '%%EndComments',/,/,/,'%',/,'%',/,'%',/,
     +    /, '/m {moveto} def /l {lineto} def /r {rlineto} def',
     +    /, '/s {stroke} def /n {newpath} def /c {closepath} def',
     +    /, '/rtext{ dup stringwidth pop -1 div 0 rmoveto show } def',
     +    /, '/ctext{ dup stringwidth pop -2 div 0 rmoveto show } def',
     +    /, '/ltext{show} def /gr{grestore} def /gs{gsave} def',
     +    /, '/tr{translate} def /setc{setrgbcolor} def',
     +    /, '/bullet{ 6 0 360 arc c fill } def',/,/,
     +    /, '%72 72 translate',/,/,
     +    /, '0.240000 0.240000 scale')
c
c Adjust the scaling parameters if undefined or if log scaling is used:
c
      if(pmax.le.pmin) then
            pmin = ar1(1)
            pmax = ar1(nd)
            pinc = (pmax-pmin)/10.0
            write(*,*) 'Estimating limits for the plot:',pmin,pmax
      endif
      if(ilog.eq.1) then
c           pmin   = max(ar1(1), 0.0001)
c           pmax   = max(ar1(nd),0.0001)
            istart = int(alog10(pmin)-0.999)
            start  = 10.0**istart
            iend   = int(alog10(pmax)+0.999)
            ncyc   = iend - istart
            pmin   = real(istart)
            pmax   = real(iend)
      endif
c
c Initialize:
c
      xhsmin = pmin
      xhsmax = pmax
      pdp = dble(pval(1))
      call gauinv(pdp,yhsmin,ierr)
      pdp = dble(pval(npval))
      call gauinv(pdp,yhsmax,ierr)
      xrange = hpxmax - hpxmin
      yrange = hpymax - hpymin
c
c Write the title and the labels:
c
      ts   = 7.5
      xloc = hpxmin - 0.15*xrange
      yloc = hpymin + 0.50*yrange
      call pstext(xloc,yloc,22,'Cumulative Probability',ts,
     +                      1,90.0,1)
      xloc = hpxmin + 0.50*xrange
      yloc = hpymin - 0.15*yrange
      call pstext(xloc,yloc,24,varlab,ts,1,0.0,1)
      yloc = hpymax + 0.01*(hpxmax-hpxmin)
      call pstext(hpxmin,yloc,40,title,9.,3,0.,0)
c
c Draw a border box:
c
      xx(1) = hpxmin
      yy(1) = hpymin
      xx(2) = hpxmin
      yy(2) = hpymax
      xx(3) = hpxmax
      yy(3) = hpymax
      xx(4) = hpxmax
      yy(4) = hpymin
      xx(5) = hpxmin
      yy(5) = hpymin
      np    = 5
      siz   = 0.7
      idsh  = 0
      call psline(np,xx,yy,siz,idsh)
c
c Scale and Draw the Probability Axis:
c
      np    = 2
      siz   = 0.4
      xx(1) = hpxmin
      xx(2) = hpxmax
      xloc  = hpxmin - 0.01*xrange
      ts    = 7.0
      do i=1,npval
            pdp = dble(pval(i))
            call gauinv(pdp,yloc,ierr)
            yy(1) = resc(yhsmin,yhsmax,hpymin,hpymax,yloc)
            yy(2) = yy(1)
            call psline(np,xx,yy,siz,idsh)
            yloc = yy(1) - 0.01*yrange
            cp   = pval(i) * 100.0
            write(title(1:5),'(i5)') int(cp)
            if(i.eq.1.or.i.eq.npval) write(title(1:5),'(f5.2)') cp
            if(i.eq.2.or.i.eq.3.or.i.eq.19.or.i.eq.20)
     +          write(title(1:5),'(f5.1)') cp
            call pstext(xloc,yloc,5,title,ts,1,0.0,2)
      end do
c
c Scale and Draw the Variable Axis:
c
      yy(1) = hpymin
      yy(2) = hpymax
      yloc  = hpymin - 0.05*xrange
      if(ilog.eq.1) then
c
c      LOGARITHMIC SCALE:
c
            start  = start / 10.0
            do i=1,ncyc+1
                  start = start * 10.0
                  xloc = alog10(start)
                  xloc = resc(xhsmin,xhsmax,hpxmin,hpxmax,xloc)
                  a = abs(start)
                  s = start
                           write(title(1:8),'(f8.0)')s
                  if(a.lt.10) write(title(1:8),'(f8.1)')s
                  if(a.lt.1.) write(title(1:8),'(f8.2)')s
                  if(a.lt.0.1) write(title(1:8),'(f8.3)')s
                  if(a.lt.0.01) write(title(1:8),'(f8.4)')s
                  if(a.lt.0.001) write(title(1:8),'(f8.5)')s
                  if(a.lt.0.0001) write(title(1:8),'(f8.6)')s
                  if(a.lt.0.00001) write(title(1:8),'(f8.7)')s
                  call pstext(xloc,yloc,8,title,ts,1,0.0,1)
                  if(i.le.ncyc) then
                        do j=1,9
                              xloc = alog10(start*j)
                              xx(1) = resc(xhsmin,xhsmax,
     +                                          hpxmin,hpxmax,xloc)
                              xx(2) = xx(1)
                              call psline(np,xx,yy,siz,idsh)
                        end do
                  endif
            end do
      else
c
c      ARITHMETIC SCALE:
c
            nloop = int((pmax-pmin)/pinc) + 1
            temp  = pmin-pinc
            do i=1,nloop
                  temp = temp + pinc
                  xloc = resc(xhsmin,xhsmax,hpxmin,hpxmax,temp)
                  s = temp
                  a = abs(temp)
                           write(title(1:8),'(f8.0)')s
                  if(a.lt.10) write(title(1:8),'(f8.1)')s
                  if(a.lt.1.) write(title(1:8),'(f8.2)')s
                  if(a.lt.0.1) write(title(1:8),'(f8.3)')s
                  if(a.lt.0.01) write(title(1:8),'(f8.4)')s
                  if(a.lt.0.001) write(title(1:8),'(f8.5)')s
                  if(a.lt.0.0001) write(title(1:8),'(f8.6)')s
                  if(a.lt.0.00001) write(title(1:8),'(f8.7)')s
                  if(s.eq.0.0) write(title(1:8),'(f8.0)')s
                  call pstext(xloc,yloc,8,title,ts,1,0.0,1)
                  xx(1) = xloc
                  xx(2) = xx(1)
                  call psline(np,xx,yy,siz,idsh)
            end do
      endif
c
c Scale and Draw the points on the plot:
c
      nloop   =  nd
      plotall = .true.
      if(npts.le.0) npts = nd
      if(npts.lt.nd) then
            nloop   =  npts
            plotall = .false.
            write(*,901) nloop
 901        format(' Plotting  ',i8,' equally spaced quantiles',/)
      end if
      do i=1,nloop
            if(plotall) then
                  xloc =      ar1(i)
                  pdp  = dble(ar2(i))
            else
                  pp   = (real(i)-0.5)/real(nloop)
                  pdp  = dble(pp)
                  call locate(ar2,nd,1,nd,pp,j)
                  if(j.eq.0) then
                        xloc = ar1(1)
                  else if(j.eq.nd) then
                        xloc = ar1(nd)
                  else
                        xloc = ar1(j) +      (ar1(j+1)-ar1(j)) *
     +                                  (pp-ar2(j))/(ar2(j+1)-ar2(j))
                  end if
            end if
            if(ilog.ne.0) then
                  if(xloc.le.EPSLON) go to 10
                  xloc  = alog10(xloc)
            endif
            call gauinv(pdp,yloc,ierr)
            if(xloc.lt.xhsmin.or.xloc.gt.xhsmax.or.
     +         yloc.lt.yhsmin.or.yloc.gt.yhsmax) go to 10
            xx(1) = resc(xhsmin,xhsmax,hpxmin,hpxmax,xloc)
            yy(1) = resc(yhsmin,yhsmax,hpymin,hpymax,yloc)
            ix = int((resc(xmin,xmax,pxmin,pxmax,xx(1)))/ pscl)
            iy = int((resc(ymin,ymax,pymin,pymax,yy(1)))/ pscl)
            write(lpsout,101) ix,iy
 101        format('n ',i4.4,1x,i4.4,' bullet')
 10         continue
      end do
c
c Add a footer to the Postscript plot file:
c
      write(lpsout,999)
 999  format('%END OF POSTSCRIPT FILE',/,'4.166667 4.166667 scale',/,/,
     +       '%%EOF',/,'showpage')
c
c Finished:
c
      close(lpsout)
      write(*,9998) VERSION
 9998 format(/' PROBPLT Version: ',f5.3, ' Finished'/)
      stop
 97   stop 'ERROR in parameters somewhere'
 99   stop 'ERROR in data somewhere'
      end

      subroutine makepar
c-----------------------------------------------------------------------
c
c                      Write a Parameter File
c                      **********************
c
c
c
c-----------------------------------------------------------------------
      lun = 99
      open(lun,file='probplt.par',status='UNKNOWN')
      write(lun,10)
 10   format('                  Parameters for PROBPLT',/,
     +       '                  **********************',/,/,
     +       'START OF PARAMETERS:')

      write(lun,11)
 11   format('../data/cluster.dat            ',
     +       '-file with data')
      write(lun,12)
 12   format('3   5                          ',
     +       '-  columns for variable and weight')
      write(lun,13)
 13   format('-1.0  1.0e21                   ',
     +       '-  trimming limits')
      write(lun,14)
 14   format('probplt.ps                     ',
     +       '-file for PostScript output')
      write(lun,15)
 15   format('0                              ',
     +       '-number of points to plot (<0 for all)')
      write(lun,16)
 16   format('0                              ',
     +       '-0=arithmetic, 1=log scaling')
      write(lun,17)
 17   format('0.0   30.0   5.0               ',
     +       '-min,max,increment for labeling')
      write(lun,18)
 18   format('Clustered Data                 ',
     +       '-title')

      close(lun)
      return
      end

ext: for


text:
      SUBROUTINE DSWAP(N,X,INCX,Y,INCY)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 2006-06-07 DSWAP  Krogh   Removed arithmetic ifs
C>> 1994-11-11 DSWAP  Krogh   Declared all vars.
C>> 1994-10-20 DSWAP  Krogh   Changes to use M77CON
C>> 1985-08-02 DSWAP  Lawson  Initial code.
c--D replaces "?": ?SWAP
C
C     INTERCHANGE X and Y.
C     FOR I = 0 TO N-1, INTERCHANGE  X(LX+I*INCX) AND Y(LY+I*INCY),
C     WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS
C     DEFINED IN A SIMILAR WAY USING INCY.
C
      INTEGER N, INCX, INCY, IX, IY, I, M, MP1, NS
      DOUBLE PRECISION X(*),Y(*),TEMP1,TEMP2,TEMP3
      IF(N.LE.0)RETURN
      IF ((INCX .NE. INCY) .OR. (INCX .LT. 0)) THEN
C       CODE FOR UNEQUAL OR NONPOSITIVE INCREMENTS.
        IX = 1
        IY = 1
        IF(INCX.LT.0)IX = (-N+1)*INCX + 1
        IF(INCY.LT.0)IY = (-N+1)*INCY + 1
        DO 10 I = 1,N
          TEMP1 = X(IX)
          X(IX) = Y(IY)
          Y(IY) = TEMP1
          IX = IX + INCX
          IY = IY + INCY
 10     CONTINUE
      ELSE IF (INCX .EQ. 1) THEN
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C       CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 3.
        M = MOD(N,3)
        IF( M .NE. 0 ) THEN
          DO 30 I = 1,M
            TEMP1 = X(I)
            X(I) = Y(I)
            Y(I) = TEMP1
 30       CONTINUE
          IF( N .LT. 3 ) RETURN
        end if
        MP1 = M + 1
        DO 50 I = MP1,N,3
          TEMP1 = X(I)
          TEMP2 = X(I+1)
          TEMP3 = X(I+2)
          X(I) = Y(I)
          X(I+1) = Y(I+1)
          X(I+2) = Y(I+2)
          Y(I) = TEMP1
          Y(I+1) = TEMP2
          Y(I+2) = TEMP3
 50     CONTINUE
      ELSE
C     CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS.
        NS = N*INCX
        DO 70 I=1,NS,INCX
          TEMP1 = X(I)
          X(I) = Y(I)
          Y(I) = TEMP1
 70     CONTINUE
      END IF
      RETURN
      END


text:
module types

    implicit none

    ! Nonrelativistic operator coefficients
    type coeff
        real(kind=8), dimension(:), allocatable :: c
    end type coeff

    type eftheory
        type(coeff) :: isoc(0:1) !c0 = (cp+cn)/2   c1 = (cp-cn)/2
        type(coeff) :: xpnc(0:1) !0 = protons, 1 = neutrons
    end type eftheory

    ! Dark matter particles
    type particle
        real(kind=8) :: j = 0.5d0! darkmatter spin
        real(kind=8) :: mass = 50.0d0
        real(kind=8) :: localdensity = 0.3d0
    end type particle

    type onebdenmat
        logical good
        real(kind=8), allocatable :: rho(:,:,:,:)
    end type onebdenmat

    type state
        integer Jx2 ! spin times 2
        integer Tx2 ! isospin times 2
    end type state

    type nucleus
        integer :: Z
        integer :: N
        integer :: A
        integer :: Mt
        real(kind=8) :: mass
        real(kind=8) :: Nt ! mass density of target nuclei
        type(state) :: groundstate

        type (onebdenmat) :: densitymats
        logical :: evenA
    end type nucleus

end module types

max_stars_repo_name: ogorton/dmfortfactor

max_issues_repo_name: ogorton/dmfortfactor

max_forks_repo_name: ogorton/dmfortfactor


MODULE module_sf_gfs

CONTAINS

   SUBROUTINE SF_GFS(U3D,V3D,T3D,QV3D,P3D, &
   CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
                     ZNT,UST,PSIM,PSIH, &
   XLAND,HFX,QFX,LH,TSK,FLHC,FLQC, &
                     QGH,QSFC,U10,V10, &
                     GZ1OZ0,WSPD,BR,ISFFLX, &
                     EP1,EP2,KARMAN,itimestep, &
                     ids,ide, jds,jde, kds,kde, &
                     ims,ime, jms,jme, kms,kme, &
                     its,ite, jts,jte, kts,kte )

      USE MODULE_GFS_MACHINE, ONLY : kind_phys
      USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, &
                                        ISFFLX,itimestep
      REAL, INTENT(IN) :: &
                                        CP, &
                                        EP1, &
                                        EP2, &
                                        KARMAN, &
                                        R, &
                                        ROVCP, &
                                        XLV
      REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
                                        P3D, &
                                        QV3D, &
                                        T3D, &
                                        U3D, &
                                        V3D
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
                                        TSK, &
                                        PSFC, &
                                        XLAND
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
                                        UST, &
                                        ZNT
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
                                        BR, &
                                        CHS, &
                                        CHS2, &
                                        CPM, &
                                        CQS2, &
                                        FLHC, &
                                        FLQC, &
                                        GZ1OZ0, &
                                        HFX, &
                                        LH, &
                                        PSIM, &
                                        PSIH, &
                                        QFX, &
                                        QGH, &
                                        QSFC, &
                                        U10, &
                                        V10, &
                                        WSPD
      REAL :: ESAT
      REAL (kind=kind_phys) :: &
                                        RHOX
      REAL (kind=kind_phys), DIMENSION(its:ite) :: &
                                        CH, &
                                        CM, &
                                        DDVEL, &
                                        DRAIN, &
                                        EP, &
                                        EVAP, &
                                        FH, &
                                        FH2, &
                                        FM, &
                                        HFLX, &
                                        PH, &
                                        PM, &
                                        PRSL1, &
                                        PRSLKI, &
                                        PS, &
                                        Q1, &
                                        Q2M, &
                                        QSS, &
                                        QSURF, &
                                        RB, &
                                        RCL, &
                                        RHO1, &
                                        SLIMSK, &
                                        STRESS, &
                                        T1, &
                                        T2M, &
                                        THGB, &
                                        THX, &
                                        TSKIN, &
                                        SHELEG, &
                                        U1, &
                                        U10M, &
                                        USTAR, &
                                        V1, &
                                        V10M, &
                                        WIND, &
                                        Z0RL, &
                                        Z1
      INTEGER :: &
                                        I, &
                                        IM, &
                                        J, &
                                        K, &
                                        KM
   if(itimestep.eq.0) then
     CALL GFUNCPHYS
   endif
   IM=ITE-ITS+1
   KM=KTE-KTS+1
   DO J=jts,jte
      DO i=its,ite
        DDVEL(I)=0.
        RCL(i)=1.
        PRSL1(i)=P3D(i,kts,j)*.001
        PS(i)=PSFC(i,j)*.001
        Q1(I) = QV3D(i,kts,j)
        QSURF(I)=0.
        SHELEG(I)=0.
        SLIMSK(i)=ABS(XLAND(i,j)-2.)
        TSKIN(i)=TSK(i,j)
        T1(I) = T3D(i,kts,j)
        U1(I) = U3D(i,kts,j)
        USTAR(I) = UST(i,j)
        V1(I) = V3D(i,kts,j)
        Z0RL(I) = ZNT(i,j)*100.
      ENDDO
      DO i=its,ite
         PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
         THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP
         THX(I)=T1(i)*(100./PRSL1(I))**ROVCP
         RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I)))
         Q1(I)=Q1(I)/(1.+Q1(I))
      ENDDO
      CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
                  SHELEG,TSKIN,QSURF, &
                  Z0RL, &
                  U10M,V10M,T2M,Q2M, &
                  CM,CH,RB, &
                  RCL,PRSL1,PRSLKI,SLIMSK, &
                  DRAIN,EVAP,HFLX,STRESS,EP, &
                  FM,FH,USTAR,WIND,DDVEL, &
                  PM,PH,FH2,QSS,Z1 )
      DO i=its,ite
        U10(i,j)=U10M(i)
        V10(i,j)=V10M(i)
        BR(i,j)=RB(i)
        CHS(I,J)=CH(I)*WIND(I)
        CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
        CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
        esat = fpvs(t1(i))
        QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
        QSFC(I,J)=qss(i)
        PSIH(i,j)=PH(i)
        PSIM(i,j)=PM(i)
        UST(i,j)=ustar(i)
        WSPD(i,j)=WIND(i)
        ZNT(i,j)=Z0RL(i)*.01
      ENDDO
      DO i=its,ite
        FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J)
        FLQC(i,j)=RHO1(I)*CHS(I,J)
        GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01))
        CQS2(i,j)=CHS2(I,J)
      ENDDO
      IF (ISFFLX.EQ.0) THEN
        DO i=its,ite
          HFX(i,j)=0.
          LH(i,j)=0.
          QFX(i,j)=0.
        ENDDO
      ELSE
        DO i=its,ite
          IF(XLAND(I,J)-1.5.GT.0.)THEN
            HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
          ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
            HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
            HFX(I,J)=AMAX1(HFX(I,J),-250.)
          ENDIF
          QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
          QFX(I,J)=AMAX1(QFX(I,J),0.)
          LH(I,J)=XLV*QFX(I,J)
        ENDDO
      ENDIF
    ENDDO
   END SUBROUTINE SF_GFS
      SUBROUTINE PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
     & SHELEG,TSKIN,QSURF, &
     & Z0RL, &
     & U10M,V10M,T2M,Q2M, &
     & CM, CH, RB, &
     & RCL,PRSL1,PRSLKI,SLIMSK, &
     & DRAIN,EVAP,HFLX,STRESS,EP, &
     & FM,FH,USTAR,WIND,DDVEL, &
     & PM,PH,FH2,QSS,Z1 )
      USE MODULE_GFS_MACHINE, ONLY : kind_phys
      USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
      USE MODULE_GFS_PHYSCONS, grav => con_g, SBC => con_sbc, HVAP => con_HVAP &
     &, CP => con_CP, HFUS => con_HFUS, JCAL => con_JCAL &
     &, EPS => con_eps, EPSM1 => con_epsm1, t0c => con_t0c &
     &, RVRDM1 => con_FVirt, RD => con_RD
      implicit none
      integer IM, km
      real(kind=kind_phys), parameter :: cpinv=1.0/cp, HVAPI=1.0/HVAP
      real(kind=kind_phys) DELT
      INTEGER SOILTYP(IM), VEGTYPE(IM)
      real(kind=kind_phys) PS(IM), U1(IM), V1(IM), &
     & T1(IM), Q1(IM), SHELEG(IM), &
     & TSKIN(IM), QSURF(IM), SMC(IM,KM), &
     & STC(IM,KM), DM(IM), SIGMAF(IM), &
     & CANOPY(IM), DLWFLX(IM), SLRAD(IM), &
     & SNOWMT(IM), Z0RL(IM), TG3(IM), &
     & GFLUX(IM), F10M(IM), U10M(IM), &
     & V10M(IM), T2M(IM), Q2M(IM), &
     & ZSOIL(IM,KM), CM(IM), CH(IM), &
     & RB(IM), RHSCNPY(IM), RHSMC(IM,KM), &
     & AIM(IM,KM), BIM(IM,KM), CIM(IM,KM), &
     & RCL(IM), PRSL1(IM), PRSLKI(IM), &
     & SLIMSK(IM), DRAIN(IM), EVAP(IM), &
     & HFLX(IM), RNET(IM), EP(IM), &
     & FM(IM), FH(IM), USTAR(IM), &
     & WIND(IM), DDVEL(IM), STRESS(IM)
      integer k,i
      real(kind=kind_phys) CANFAC(IM), &
     & DDZ(IM), DDZ2(IM), DELTA(IM), &
     & DEW(IM), DF1(IM), DFT0(IM), &
     & DFT2(IM), DFT1(IM), &
     & DMDZ(IM), DMDZ2(IM), DTDZ1(IM), &
     & DTDZ2(IM), DTV(IM), EC(IM), &
     & EDIR(IM), ETPFAC(IM), &
     & FACTSNW(IM), FH2(IM), FM10(IM), &
     & FX(IM), GX(IM), &
     & HCPCT(IM), HL1(IM), HL12(IM), &
     & HLINF(IM), PARTLND(IM), PH(IM), &
     & PH2(IM), PM(IM), PM10(IM), &
     & PSURF(IM), Q0(IM), QS1(IM), &
     & QSS(IM), RAT(IM), RCAP(IM), &
     & RCH(IM), RHO(IM), RS(IM), &
     & RSMALL(IM), SLWD(IM), SMCZ(IM), &
     & SNET(IM), SNOEVP(IM), SNOWD(IM), &
     & T1O(IM), T2MO(IM), TERM1(IM), &
     & TERM2(IM), THETA1(IM), THV1(IM), &
     & TREF(IM), TSURF(IM), TV1(IM), &
     & TVS(IM), TSURFO(IM), TWILT(IM), &
     & XX(IM), XRCL(IM), YY(IM), &
     & Z0(IM), Z0MAX(IM), Z1(IM), &
     & ZTMAX(IM), ZZ(IM), PS1(IM)
      real(kind=kind_phys) a0, a0p, a1, a1p, aa, aa0, &
     & aa1, adtv, alpha, arnu, b1, b1p, &
     & b2, b2p, bb, bb0, bb1, bb2, &
     & bfact, ca, cc, cc1, cc2, cfactr, &
     & ch2o, charnock, cice, convrad, cq, csoil, &
     & ctfil1,ctfil2, delt2, df2, dfsnow, &
     & elocp, eth, ff, FMS, &
     & fhs, g, hl0, hl0inf, &
     & hl110, hlt, hltinf,OLINF, rcq, rcs, &
     & rct, restar, rhoh2o,rnu, RSI, &
     & rss, scanop, sig2k, sigma, smcdry, &
     & t12, t14, tflx, tgice, topt, &
     & val, vis, zbot, snomin, tem
      PARAMETER (CHARNOCK=.014,CA=.4)
      PARAMETER (G=grav,sigma=sbc)
      PARAMETER (ALPHA=5.,A0=-3.975,A1=12.32,B1=-7.755,B2=6.041)
      PARAMETER (A0P=-7.941,A1P=24.75,B1P=-8.705,B2P=7.899,VIS=1.4E-5)
      PARAMETER (AA1=-1.076,BB1=.7045,CC1=-.05808)
      PARAMETER (BB2=-.1954,CC2=.009999)
      PARAMETER (ELOCP=HVAP/CP,DFSNOW=.31,CH2O=4.2E6,CSOIL=1.26E6)
      PARAMETER (SCANOP=.5,CFACTR=.5,ZBOT=-3.,TGICE=271.2)
      PARAMETER (CICE=1880.*917.,topt=298.)
      PARAMETER (RHOH2O=1000.,CONVRAD=JCAL*1.E4/60.)
      PARAMETER (CTFIL1=.5,CTFIL2=1.-CTFIL1)
      PARAMETER (RNU=1.51E-5,ARNU=.135*RNU)
      parameter (snomin=1.0e-9)
      LOGICAL FLAG(IM), FLAGSNW(IM)
      real(kind=kind_phys) KT1(IM), KT2(IM), &
     & ET(IM,KM), &
     & STSOIL(IM,KM), AI(IM,KM), BI(IM,KM), &
     & CI(IM,KM), RHSTC(IM,KM)
      real(kind=kind_phys) rsmax(13), rgl(13), rsmin(13), hs(13), &
     & smmax(9), smdry(9), smref(9), smwlt(9)
      data rsmax/13*5000./
      data rsmin/150.,100.,125.,150.,100.,70.,40., &
     & 300.,400.,150.,999.,40.,999./
      data rgl/5*30.,65.,4*100.,999.,100.,999./
      data hs/41.69,54.53,51.93,47.35,47.35,54.53,36.35, &
     & 3*42.00,999.,36.35,999./
      data smmax/.421,.464,.468,.434,.406,.465,.404,.439,.421/
      data smdry/.07,.14,.22,.08,.18,.16,.12,.10,.07/
      data smref/.283,.387,.412,.312,.338,.382,.315,.329,.283/
      data smwlt/.029,.119,.139,.047,.010,.103,.069,.066,.029/
      SIG2K = 1. - 4. * G * 2. / (CP * 280.)
      DO I=1,IM
        XRCL(I) = SQRT(RCL(I))
        PSURF(I) = 1000. * PS(I)
        PS1(I) = 1000. * PRSL1(I)
        WIND(I) = XRCL(I) * SQRT(U1(I) * U1(I) + V1(I) * V1(I)) &
     & + MAX(0.0_kind_phys, MIN(DDVEL(I), 30.0_kind_phys))
        WIND(I) = MAX(WIND(I),1._kind_phys)
        Q0(I) = MAX(Q1(I),1.E-8_kind_phys)
        TSURF(I) = TSKIN(I)
        THETA1(I) = T1(I) * PRSLKI(I)
        TV1(I) = T1(I) * (1. + RVRDM1 * Q0(I))
        THV1(I) = THETA1(I) * (1. + RVRDM1 * Q0(I))
        TVS(I) = TSURF(I) * (1. + RVRDM1 * Q0(I))
        RHO(I) = PS1(I) / (RD * TV1(I))
        qs1(i) = fpvs(t1(i))
        QS1(I) = EPS * QS1(I) / (PS1(I) + EPSM1 * QS1(I))
        QS1(I) = MAX(QS1(I), 1.E-8_kind_phys)
        Q0(I) = min(QS1(I),Q0(I))
        qss(i) = fpvs(tskin(i))
        QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
        RS(I) = 0.
        Z0(I) = .01 * Z0RL(i)
        DM(I) = 1.
      GOTO 1111
        FACTSNW(I) = 10.
        IF(SLIMSK(I).EQ.2.) FACTSNW(I) = 3.
        SNOWD(I) = SHELEG(I) / 1000.
        FLAGSNW(I) = .FALSE.
        IF(SNOWD(I).GT..001.OR.SLIMSK(I).EQ.2.) RS(I) = 0.
        IF(SNOWD(I).GT..001) FLAGSNW(I) = .TRUE.
        IF(SLIMSK(I).EQ.0.) THEN
          ZSOIL(I,1) = 0.
        ELSEIF(SLIMSK(I).EQ.1.) THEN
          ZSOIL(I,1) = -.10
        ELSE
          ZSOIL(I,1) = -3. / KM
        ENDIF
1111 CONTINUE
      ENDDO
      GOTO 2222
      DO K = 2, KM
        DO I=1,IM
          IF(SLIMSK(I).EQ.0.) THEN
            ZSOIL(I,K) = 0.
          ELSEIF(SLIMSK(I).EQ.1.) THEN
            ZSOIL(I,K) = ZSOIL(I,K-1) &
     & + (-2. - ZSOIL(I,1)) / (KM - 1)
          ELSE
            ZSOIL(I,K) = - 3. * FLOAT(K) / FLOAT(KM)
          ENDIF
        ENDDO
      ENDDO
2222 CONTINUE
      DO I=1,IM
        Z1(I) = -RD * TV1(I) * LOG(PS1(I)/PSURF(I)) / G
        DRAIN(I) = 0.
      ENDDO
      DO K = 1, KM
        DO I=1,IM
          ET(I,K) = 0.
          RHSMC(I,K) = 0.
          AIM(I,K) = 0.
          BIM(I,K) = 1.
          CIM(I,K) = 0.
          STSOIL(I,K) = STC(I,K)
        ENDDO
      ENDDO
      DO I=1,IM
        EDIR(I) = 0.
        EC(I) = 0.
        EVAP(I) = 0.
        EP(I) = 0.
        SNOWMT(I) = 0.
        GFLUX(I) = 0.
        RHSCNPY(I) = 0.
        FX(I) = 0.
        ETPFAC(I) = 0.
        CANFAC(I) = 0.
      ENDDO
      DO I=1,IM
        IF(SLIMSK(I).EQ.0.) THEN
          USTAR(I) = SQRT(G * Z0(I) / CHARNOCK)
        ENDIF
        Z0MAX(I) = MIN(Z0(I),0.1 * Z1(I))
        ZTMAX(I) = Z0MAX(I)
        IF(SLIMSK(I).EQ.0.) THEN
          RESTAR = USTAR(I) * Z0MAX(I) / VIS
          RESTAR = MAX(RESTAR,.000001_kind_phys)
          RAT(I) = 2.67 * restar ** .25 - 2.57
          RAT(I) = min(RAT(I),7._kind_phys)
          ZTMAX(I) = Z0MAX(I) * EXP(-RAT(I))
        ENDIF
      ENDDO
      DO I = 1, IM
        DTV(I) = THV1(I) - TVS(I)
        ADTV = ABS(DTV(I))
        ADTV = MAX(ADTV,.001_kind_phys)
        DTV(I) = SIGN(1._kind_phys,DTV(I)) * ADTV
        RB(I) = G * DTV(I) * Z1(I) / (.5 * (THV1(I) + TVS(I)) &
     & * WIND(I) * WIND(I))
        RB(I) = MAX(RB(I),-5000._kind_phys)
        FM(I) = LOG((Z1(I)) / Z0MAX(I))
        FH(I) = LOG((Z1(I)) / ZTMAX(I))
        HLINF(I) = RB(I) * FM(I) * FM(I) / FH(I)
        FM10(I) = LOG((Z0MAX(I)+10.) / Z0MAX(I))
        FH2(I) = LOG((ZTMAX(I)+2.) / ZTMAX(I))
      ENDDO
      DO I = 1, IM
        IF(DTV(I).GE.0.) THEN
          HL1(I) = HLINF(I)
        ENDIF
        IF(DTV(I).GE.0..AND.HLINF(I).GT..25) THEN
          HL0INF = Z0MAX(I) * HLINF(I) / Z1(I)
          HLTINF = ZTMAX(I) * HLINF(I) / Z1(I)
          AA = SQRT(1. + 4. * ALPHA * HLINF(I))
          AA0 = SQRT(1. + 4. * ALPHA * HL0INF)
          BB = AA
          BB0 = SQRT(1. + 4. * ALPHA * HLTINF)
          PM(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
          PH(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
          FMS = FM(I) - PM(I)
          FHS = FH(I) - PH(I)
          HL1(I) = FMS * FMS * RB(I) / FHS
        ENDIF
      ENDDO
      DO I = 1, IM
        IF(DTV(I).GE.0.) THEN
          HL0 = Z0MAX(I) * HL1(I) / Z1(I)
          HLT = ZTMAX(I) * HL1(I) / Z1(I)
          AA = SQRT(1. + 4. * ALPHA * HL1(I))
          AA0 = SQRT(1. + 4. * ALPHA * HL0)
          BB = AA
          BB0 = SQRT(1. + 4. * ALPHA * HLT)
          PM(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
          PH(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
          HL110 = HL1(I) * 10. / Z1(I)
          AA = SQRT(1. + 4. * ALPHA * HL110)
          PM10(I) = AA0 - AA + LOG((AA + 1.) / (AA0 + 1.))
          HL12(I) = HL1(I) * 2. / Z1(I)
          BB = SQRT(1. + 4. * ALPHA * HL12(I))
          PH2(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
        ENDIF
      ENDDO
      DO I=1,IM
        IF(DTV(I).LT.0.) THEN
          OLINF = Z1(I) / HLINF(I)
          IF(ABS(OLINF).LE.50. * Z0MAX(I)) THEN
            HLINF(I) = -Z1(I) / (50. * Z0MAX(I))
          ENDIF
        ENDIF
      ENDDO
      DO I = 1, IM
        IF(DTV(I).LT.0..AND.HLINF(I).GE.-.5) THEN
          HL1(I) = HLINF(I)
          PM(I) = (A0 + A1 * HL1(I)) * HL1(I) &
     & / (1. + B1 * HL1(I) + B2 * HL1(I) * HL1(I))
          PH(I) = (A0P + A1P * HL1(I)) * HL1(I) &
     & / (1. + B1P * HL1(I) + B2P * HL1(I) * HL1(I))
          HL110 = HL1(I) * 10. / Z1(I)
          PM10(I) = (A0 + A1 * HL110) * HL110 &
     & / (1. + B1 * HL110 + B2 * HL110 * HL110)
          HL12(I) = HL1(I) * 2. / Z1(I)
          PH2(I) = (A0P + A1P * HL12(I)) * HL12(I) &
     & / (1. + B1P * HL12(I) + B2P * HL12(I) * HL12(I))
        ENDIF
        IF(DTV(I).LT.0.AND.HLINF(I).LT.-.5) THEN
          HL1(I) = -HLINF(I)
          PM(I) = LOG(HL1(I)) + 2. * HL1(I) ** (-.25) - .8776
          PH(I) = LOG(HL1(I)) + .5 * HL1(I) ** (-.5) + 1.386
          HL110 = HL1(I) * 10. / Z1(I)
          PM10(I) = LOG(HL110) + 2. * HL110 ** (-.25) - .8776
          HL12(I) = HL1(I) * 2. / Z1(I)
          PH2(I) = LOG(HL12(I)) + .5 * HL12(I) ** (-.5) + 1.386
        ENDIF
      ENDDO
      DO I = 1, IM
        FM(I) = FM(I) - PM(I)
        FH(I) = FH(I) - PH(I)
        FM10(I) = FM10(I) - PM10(I)
        FH2(I) = FH2(I) - PH2(I)
        CM(I) = CA * CA / (FM(I) * FM(I))
        CH(I) = CA * CA / (FM(I) * FH(I))
        CQ = CH(I)
        STRESS(I) = CM(I) * WIND(I) * WIND(I)
        USTAR(I) = SQRT(STRESS(I))
      ENDDO
      DO I = 1, IM
        IF(SLIMSK(I).EQ.0.) THEN
          Z0(I) = (CHARNOCK / G) * USTAR(I) ** 2
          Z0(I) = MIN(Z0(I),.1_kind_phys)
          Z0(I) = MAX(Z0(I),1.E-7_kind_phys)
          Z0RL(I) = 100. * Z0(I)
        ENDIF
      ENDDO
   GOTO 5555
      DO I = 1, IM
        RCH(I) = RHO(I) * CP * CH(I) * WIND(I)
      ENDDO
      DO I = 1, IM
        IF(SLIMSK(I).EQ.0.) THEN
          EVAP(I) = ELOCP * RCH(I) * (QSS(I) - Q1(I))
          DM(I) = 1.
          QSURF(I) = QSS(I)
        ENDIF
      ENDDO
     DO I = 1, IM
       GFLUX(I) = 0.
       IF(SLIMSK(I).EQ.1.) THEN
         SMCZ(I) = .5 * (SMC(I,1) + .20)
         DFT0(I) = KTSOIL(SMCZ(I),SOILTYP(I))
       ELSEIF(SLIMSK(I).EQ.2.) THEN
         DFT0(I) = 2.2
       ENDIF
     ENDDO
     DO I=1,IM
       IF(SLIMSK(I).NE.0.) THEN
         IF(FLAGSNW(I)) THEN
  TFLX = MIN(T1(I), TSURF(I))
  GFLUX(I) = -DFSNOW * (TFLX - STSOIL(I,1)) &
    & / (FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
         ELSE
  GFLUX(I) = DFT0(I) * (STSOIL(I,1) - TSURF(I)) &
    & / (-.5 * ZSOIL(I,1))
         ENDIF
         GFLUX(I) = MAX(GFLUX(I),-200._kind_phys)
         GFLUX(I) = MIN(GFLUX(I),+200._kind_phys)
       ENDIF
     ENDDO
     DO I = 1, IM
       FLAG(I) = SLIMSK(I).NE.0.
       PARTLND(I) = 1.
       IF(SNOWD(I).GT.0..AND.SNOWD(I).LE..001) THEN
         PARTLND(I) = 1. - SNOWD(I) / .001
       ENDIF
     ENDDO
     DO I = 1, IM
       SNOEVP(I) = 0.
       if(SNOWD(I).gt..001) PARTLND(I) = 0.
     ENDDO
     DO I = 1, IM
       IF(FLAG(I)) THEN
         T12 = T1(I) * T1(I)
         T14 = T12 * T12
         RCAP(I) = -SLWD(I) - SIGMA * T14 + GFLUX(I) &
    & - RCH(I) * (T1(I) - THETA1(I))
         RSMALL(I) = 4. * SIGMA * T1(I) * T12 / RCH(I) + 1.
         DELTA(I) = ELOCP * EPS * HVAP * QS1(I) / (RD * T12)
         TERM1(I) = ELOCP * RSMALL(I) * RCH(I)*(QS1(I)-Q0(I))
         TERM2(I) = RCAP(I) * DELTA(I)
         EP(I) = (ELOCP * RSMALL(I) * RCH(I) * (QS1(I) - Q0(I)) &
    & + RCAP(I) * DELTA(I))
         EP(I) = EP(I) / (RSMALL(I) + DELTA(I))
       ENDIF
     ENDDO
     DO I = 1, IM
       FLAG(I) = SLIMSK(I).EQ.1..AND.EP(I).GT.0.
     ENDDO
     DO I = 1, IM
       IF(FLAG(I)) THEN
         DF1(I) = FUNCDF(SMC(I,1),SOILTYP(I))
         KT1(I) = FUNCKT(SMC(I,1),SOILTYP(I))
       endif
       if(FLAG(I).and.STC(I,1).lt.t0c) then
         DF1(I) = 0.
         KT1(I) = 0.
       endif
       IF(FLAG(I)) THEN
         TREF(I) = smref(SOILTYP(I))
         TWILT(I) = smwlt(SOILTYP(I))
         smcdry = smdry(SOILTYP(I))
         FX(I) = -2. * DF1(I) * (SMC(I,1) - smcdry) / ZSOIL(I,1) &
    & - KT1(I)
         FX(I) = MIN(FX(I), EP(I)/HVAP)
         FX(I) = MAX(FX(I),0._kind_phys)
         EDIR(I) = FX(I) * (1. - SIGMAF(I)) * PARTLND(I)
       ENDIF
     ENDDO
     DO I = 1, IM
       if(FLAG(I)) then
         ff = .55 * 2. * SNET(I) / rgl(VEGTYPE(I))
         rcs = (ff + RS(I)/rsmax(VEGTYPE(I))) / (1. + ff)
         rcs = max(rcs,.0001_kind_phys)
         rct = 1.
         rcq = 1.
         RS(I) = RS(I) / (rcs * rct * rcq)
       endif
     ENDDO
     DO I = 1, IM
       IF(FLAG(I)) THEN
         CANFAC(I) = (CANOPY(I) / SCANOP) ** CFACTR
       endif
       IF(FLAG(I)) THEN
         ETPFAC(I) = SIGMAF(I) &
    & * (1. - CANFAC(I)) / HVAP
         GX(I) = (SMC(I,1) - TWILT(I)) / (TREF(I) - TWILT(I))
         GX(I) = MAX(GX(I),0._kind_phys)
         GX(I) = MIN(GX(I),1._kind_phys)
         rss = GX(I) * (ZSOIL(I,1) / ZSOIL(I,km))
         rss = max(rss,.0001_kind_phys)
         RSI = RS(I) / rss
         eth = (TERM1(I) + TERM2(I)) / &
    & (DELTA(I) + RSMALL(I) * (1. + RSI * CH(I) * WIND(I)))
         ET(I,1) = ETPFAC(I) * eth &
    & * PARTLND(I)
       ENDIF
     ENDDO
     DO K = 2, KM
       DO I=1,IM
         IF(FLAG(I)) THEN
  GX(I) = (SMC(I,K) - TWILT(I)) / (TREF(I) - TWILT(I))
  GX(I) = MAX(GX(I),0._kind_phys)
  GX(I) = MIN(GX(I),1._kind_phys)
         rss = GX(I) * ((ZSOIL(I,k) - ZSOIL(I,k-1))/ZSOIL(I,km))
         rss = max(rss,1.e-6_kind_phys)
         RSI = RS(I) / rss
         eth = (TERM1(I) + TERM2(I)) / &
    & (DELTA(I) + RSMALL(I) * (1. + RSI * CH(I) * WIND(I)))
  ET(I,K) = eth &
    & * ETPFAC(I) * PARTLND(I)
         ENDIF
       ENDDO
     ENDDO
         400 CONTINUE
     DO I=1,IM
       IF(FLAG(I)) THEN
         EC(I) = SIGMAF(I) * CANFAC(I) * EP(I) / HVAP
         EC(I) = EC(I) * PARTLND(I)
         EC(I) = min(EC(I),CANOPY(I)/delt)
       ENDIF
     ENDDO
     DO I = 1, IM
       IF(FLAG(I)) THEN
        EVAP(I) = EDIR(I) + EC(I)
       ENDIF
     ENDDO
     DO K = 1, KM
       DO I=1,IM
         IF(FLAG(I)) THEN
  EVAP(I) = EVAP(I) + ET(I,K)
         ENDIF
       ENDDO
     ENDDO
     DO I=1,IM
       IF(FLAG(I)) THEN
         EVAP(I) = MIN(EVAP(I)*HVAP,EP(I))
       ENDIF
     ENDDO
     DO I = 1, IM
       IF(SLIMSK(I).EQ.2.) THEN
         EVAP(I) = PARTLND(I) * EP(I)
       ENDIF
     ENDDO
     DO I = 1, IM
       FLAG(I) = SLIMSK(I).NE.0..AND.EP(I).LE.0.
       DEW(I) = 0.
     ENDDO
     DO I = 1, IM
       IF(FLAG(I)) THEN
         DEW(I) = -EP(I) * DELT / (HVAP * RHOH2O)
         EVAP(I) = EP(I)
         DEW(I) = DEW(I) * PARTLND(I)
         EVAP(I) = EVAP(I) * PARTLND(I)
         DM(I) = 1.
       ENDIF
     ENDDO
     DO I = 1, IM
       FLAG(I) = SLIMSK(I).NE.0..AND.SNOWD(I).GT.0.
     ENDDO
     DO I = 1, IM
       IF(FLAG(I)) THEN
         BFACT = SNOWD(I) / (DELT * EP(I) / (HVAP * RHOH2O))
         BFACT = MIN(BFACT,1._kind_phys)
         IF(EP(I).LE.0.) BFACT = 1.
         IF(SNOWD(I).LE..001) THEN
  SNOEVP(I) = bfact * EP(I)
  EVAP(I)=EVAP(I)+SNOEVP(I)*(1.-PARTLND(I))
         ELSE
  SNOEVP(I) = bfact * EP(I)
  EVAP(I) = SNOEVP(I)
         ENDIF
         TSURF(I) = T1(I) + &
    & (RCAP(I) - GFLUX(I) - DFSNOW * (T1(I) - STSOIL(I,1)) &
    & /(FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys)) &
    & - SNOEVP(I)) / (RSMALL(I) * RCH(I) &
    & + DFSNOW / (FACTSNW(I)* MAX(SNOWD(I),.001_kind_phys)))
         SNOWD(I) = SNOWD(I) - SNOEVP(I) * delt / (rhoh2o * hvap)
         SNOWD(I) = MAX(SNOWD(I),0._kind_phys)
       ENDIF
     ENDDO
         500 CONTINUE
     DO I = 1, IM
       FLAG(I) = SLIMSK(I).NE.0. &
    & .AND.SNOWD(I).GT..0
     ENDDO
     DO I = 1, IM
       IF(FLAG(I).AND.TSURF(I).GT.T0C) THEN
         SNOWMT(I) = RCH(I) * RSMALL(I) * DELT &
    & * (TSURF(I) - T0C) / (RHOH2O * HFUS)
         SNOWMT(I) = min(SNOWMT(I),SNOWD(I))
         SNOWD(I) = SNOWD(I) - SNOWMT(I)
         SNOWD(I) = MAX(SNOWD(I),0._kind_phys)
         TSURF(I) = MAX(T0C,TSURF(I) &
    & -HFUS*SNOWMT(I)*RHOH2O/(RCH(I)*RSMALL(I)*DELT))
       ENDIF
     ENDDO
     DO I = 1, IM
       IF (SNOWD(I) .GT. snomin) THEN
         qss(i) = fpvs(tsurf(i))
         QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
         EVAP(I) = elocp * RCH(I) * (QSS(I) - Q0(I))
       ENDIF
     ENDDO
     DO I = 1, IM
       FLAG(I) = SLIMSK(I).EQ.1.
       if(FLAG(I)) then
         DF1(I) = FUNCDF(SMCZ(I),SOILTYP(I))
         KT1(I) = FUNCKT(SMCZ(I),SOILTYP(I))
       endif
       if(FLAG(I).and.STC(I,1).lt.t0c) then
         DF1(I) = 0.
         KT1(I) = 0.
       endif
       IF(FLAG(I)) THEN
         RHSCNPY(I) = -EC(I) + SIGMAF(I) * RHOH2O * DEW(I) / DELT
         SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
         DMDZ(I) = (SMC(I,1) - SMC(I,2)) / (-.5 * ZSOIL(I,2))
         RHSMC(I,1) = (DF1(I) * DMDZ(I) + KT1(I) &
    & + (EDIR(I) + ET(I,1))) / (ZSOIL(I,1) * RHOH2O)
         RHSMC(I,1) = RHSMC(I,1) - (1. - SIGMAF(I)) * DEW(I) / &
    & ( ZSOIL(I,1) * delt)
         DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
         AIM(I,1) = 0.
         BIM(I,1) = DF1(I) * DDZ(I) / (-ZSOIL(I,1) * RHOH2O)
         CIM(I,1) = -BIM(I,1)
       ENDIF
     ENDDO
     DO K = 2, KM
       IF(K.LT.KM) THEN
         DO I=1,IM
  IF(FLAG(I)) THEN
    DF2 = FUNCDF(SMCZ(I),SOILTYP(I))
    KT2(I) = FUNCKT(SMCZ(I),SOILTYP(I))
  ENDIF
  IF(FLAG(I).and.STC(I,k).lt.t0c) THEN
    df2 = 0.
    KT2(I) = 0.
  ENDIF
  IF(FLAG(I)) THEN
    DMDZ2(I) = (SMC(I,K) - SMC(I,K+1)) &
    & / (.5 * (ZSOIL(I,K-1) - ZSOIL(I,K+1)))
    SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
    RHSMC(I,K) = (DF2 * DMDZ2(I) + KT2(I) &
    & - DF1(I) * DMDZ(I) - KT1(I) + ET(I,K)) &
    & / (RHOH2O*(ZSOIL(I,K) - ZSOIL(I,K-1)))
    DDZ2(I) = 2. / (ZSOIL(I,K-1) - ZSOIL(I,K+1))
    CIM(I,K) = -DF2 * DDZ2(I) &
    & / ((ZSOIL(I,K-1) - ZSOIL(I,K))*RHOH2O)
  ENDIF
         ENDDO
       ELSE
         DO I = 1, IM
  IF(FLAG(I)) THEN
    KT2(I) = FUNCKT(SMC(I,K),SOILTYP(I))
  ENDIF
  if(FLAG(I).and.STC(I,k).lt.t0c) KT2(I) = 0.
  IF(FLAG(I)) THEN
    RHSMC(I,K) = (KT2(I) &
    & - DF1(I) * DMDZ(I) - KT1(I) + ET(I,K)) &
    & / (RHOH2O*(ZSOIL(I,K) - ZSOIL(I,K-1)))
    DRAIN(I) = KT2(I)
    CIM(I,K) = 0.
  ENDIF
         ENDDO
       ENDIF
       DO I = 1, IM
         IF(FLAG(I)) THEN
  AIM(I,K) = -DF1(I) * DDZ(I) &
    & / ((ZSOIL(I,K-1) - ZSOIL(I,K))*RHOH2O)
  BIM(I,K) = -(AIM(I,K) + CIM(I,K))
  DF1(I) = DF2
  KT1(I) = KT2(I)
  DMDZ(I) = DMDZ2(I)
  DDZ(I) = DDZ2(I)
         ENDIF
       ENDDO
     ENDDO
         600 CONTINUE
     DO I=1,IM
       FLAG(I) = SLIMSK(I).NE.0.
     ENDDO
     DO I=1,IM
       IF(FLAG(I).AND..NOT.FLAGSNW(I)) THEN
         YY(I) = T1(I) + &
    & (RCAP(I)-GFLUX(I) &
    & - EVAP(I)) / (RSMALL(I) * RCH(I))
         ZZ(I) = 1. + DFT0(I) / (-.5 * ZSOIL(I,1) * RCH(I) * RSMALL(I))
         XX(I) = DFT0(I) * (STSOIL(I,1) - YY(I)) / &
    & (.5 * ZSOIL(I,1) * ZZ(I))
       ENDIF
       IF(FLAG(I).AND.FLAGSNW(I)) THEN
         YY(I) = STSOIL(I,1)
         ZZ(I) = 1.
         XX(I) = DFSNOW * (STSOIL(I,1) - TSURF(I)) &
    & / (-FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
       ENDIF
     ENDDO
     DO I = 1, IM
       IF(FLAG(I)) THEN
         SMCZ(I) = MAX(SMC(I,1), SMC(I,2))
         DTDZ1(I) = (STSOIL(I,1) - STSOIL(I,2)) / (-.5 * ZSOIL(I,2))
         IF(SLIMSK(I).EQ.1.) THEN
  DFT1(I) = KTSOIL(SMCZ(I),SOILTYP(I))
  HCPCT(I) = SMC(I,1) * CH2O + (1. - SMC(I,1)) * CSOIL
         ELSE
  DFT1(I) = DFT0(I)
  HCPCT(I) = CICE
         ENDIF
         DFT2(I) = DFT1(I)
         DDZ(I) = 1. / (-.5 * ZSOIL(I,2))
         AI(I,1) = 0.
         BI(I,1) = DFT1(I) * DDZ(I) / (-ZSOIL(I,1) * HCPCT(I))
         CI(I,1) = -BI(I,1)
         BI(I,1) = BI(I,1) &
    & + DFT0(I) / (.5 * ZSOIL(I,1) **2 * HCPCT(I) * ZZ(I))
         RHSTC(I,1) = (DFT1(I) * DTDZ1(I) - XX(I)) &
    & / (ZSOIL(I,1) * HCPCT(I))
       ENDIF
     ENDDO
     DO K = 2, KM
       DO I=1,IM
         IF(SLIMSK(I).EQ.1.) THEN
  HCPCT(I) = SMC(I,K) * CH2O + (1. - SMC(I,K)) * CSOIL
         ELSEIF(SLIMSK(I).EQ.2.) THEN
  HCPCT(I) = CICE
         ENDIF
       ENDDO
       IF(K.LT.KM) THEN
         DO I = 1, IM
  IF(FLAG(I)) THEN
    DTDZ2(I) = (STSOIL(I,K) - STSOIL(I,K+1)) &
    & / (.5 * (ZSOIL(I,K-1) - ZSOIL(I,K+1)))
    SMCZ(I) = MAX(SMC(I,K), SMC(I,K+1))
    IF(SLIMSK(I).EQ.1.) THEN
      DFT2(I) = KTSOIL(SMCZ(I),SOILTYP(I))
    ENDIF
    DDZ2(I) = 2. / (ZSOIL(I,K-1) - ZSOIL(I,K+1))
    CI(I,K) = -DFT2(I) * DDZ2(I) &
    & / ((ZSOIL(I,K-1) - ZSOIL(I,K)) * HCPCT(I))
  ENDIF
         ENDDO
       ELSE
         DO I = 1, IM
  IF(SLIMSK(I).EQ.1.) THEN
    DTDZ2(I) = (STSOIL(I,K) - TG3(I)) &
    & / (.5 * (ZSOIL(I,K-1) + ZSOIL(I,K)) - ZBOT)
    DFT2(I) = KTSOIL(SMC(I,K),SOILTYP(I))
    CI(I,K) = 0.
  ENDIF
  IF(SLIMSK(I).EQ.2.) THEN
    DTDZ2(I) = (STSOIL(I,K) - TGICE) &
    & / (.5 * ZSOIL(I,K-1) - .5 * ZSOIL(I,K))
    DFT2(I) = DFT1(I)
    CI(I,K) = 0.
  ENDIF
         ENDDO
       ENDIF
       DO I = 1, IM
         IF(FLAG(I)) THEN
  RHSTC(I,K) = (DFT2(I) * DTDZ2(I) - DFT1(I) * DTDZ1(I)) &
    & / ((ZSOIL(I,K) - ZSOIL(I,K-1)) * HCPCT(I))
  AI(I,K) = -DFT1(I) * DDZ(I) &
    & / ((ZSOIL(I,K-1) - ZSOIL(I,K)) * HCPCT(I))
  BI(I,K) = -(AI(I,K) + CI(I,K))
  DFT1(I) = DFT2(I)
  DTDZ1(I) = DTDZ2(I)
  DDZ(I) = DDZ2(I)
         ENDIF
       ENDDO
     ENDDO
         700 CONTINUE
     DO K = 1, KM
       DO I=1,IM
         IF(FLAG(I)) THEN
  RHSTC(I,K) = RHSTC(I,K) * DELT2
  AI(I,K) = AI(I,K) * DELT2
  BI(I,K) = 1. + BI(I,K) * DELT2
  CI(I,K) = CI(I,K) * DELT2
         ENDIF
       ENDDO
     ENDDO
     DO I=1,IM
       IF(FLAG(I)) THEN
         CI(I,1) = -CI(I,1) / BI(I,1)
         RHSTC(I,1) = RHSTC(I,1) / BI(I,1)
       ENDIF
     ENDDO
     DO K = 2, KM
       DO I=1,IM
         IF(FLAG(I)) THEN
  CC = 1. / (BI(I,K) + AI(I,K) * CI(I,K-1))
  CI(I,K) = -CI(I,K) * CC
  RHSTC(I,K) = (RHSTC(I,K) - AI(I,K) * RHSTC(I,K-1)) * CC
         ENDIF
       ENDDO
     ENDDO
     DO I=1,IM
       IF(FLAG(I)) THEN
         CI(I,KM) = RHSTC(I,KM)
       ENDIF
     ENDDO
     DO K = KM-1, 1
       DO I=1,IM
         IF(FLAG(I)) THEN
  CI(I,K) = CI(I,K) * CI(I,K+1) + RHSTC(I,K)
         ENDIF
       ENDDO
     ENDDO
     DO K = 1, KM
       DO I=1,IM
         IF(FLAG(I)) THEN
  STSOIL(I,K) = STSOIL(I,K) + CI(I,K)
         ENDIF
       ENDDO
     ENDDO
     DO I=1,IM
       IF(SLIMSK(I).NE.0..AND..NOT.FLAGSNW(I)) THEN
         TSURF(I) = (YY(I) + (ZZ(I) - 1.) * STSOIL(I,1)) / ZZ(I)
       ENDIF
       IF(SLIMSK(I).EQ.2..AND..NOT.FLAGSNW(I)) THEN
         TSURF(I) = MIN(TSURF(I),T0C)
       ENDIF
     ENDDO
     DO K = 1, KM
       DO I=1,IM
         IF(SLIMSK(I).EQ.2) THEN
  STSOIL(I,K) = MIN(STSOIL(I,K),T0C)
         ENDIF
       ENDDO
     ENDDO
     DO I=1,IM
       IF(SLIMSK(I).NE.0.) THEN
         TSKIN(I) = CTFIL1 * TSURF(I) + CTFIL2 * TSKIN(I)
       ENDIF
     ENDDO
     DO K = 1, KM
       DO I=1,IM
         IF(SLIMSK(I).NE.0.) THEN
  STC(I,K) = CTFIL1 * STSOIL(I,K) + CTFIL2 * STC(I,K)
         ENDIF
       ENDDO
     ENDDO
     DO I=1,IM
       FLAG(I) = SLIMSK(I).NE.0. &
    & .AND.FLAGSNW(I)
     ENDDO
     DO I = 1, IM
       IF(FLAG(I)) THEN
         GFLUX(I) = -DFSNOW * (TSKIN(I) - STC(I,1)) &
    & / (FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
       ENDIF
     ENDDO
     DO I = 1, IM
       IF( SLIMSK(I).NE.0..AND..NOT.FLAGSNW(I)) THEN
         GFLUX(I) = DFT0(I) * (STC(I,1) - TSKIN(I)) &
    & / (-.5 * ZSOIL(I,1))
       ENDIF
     ENDDO
5555 CONTINUE
      DO I = 1, IM
        F10M(I) = FM10(I) / FM(I)
        F10M(I) = min(F10M(I),1._kind_phys)
        U10M(I) = F10M(I) * XRCL(I) * U1(I)
        V10M(I) = F10M(I) * XRCL(I) * V1(I)
      ENDDO
 6000 FORMAT(8(F8.2,','))
      RETURN
      END SUBROUTINE PROGTM
      SUBROUTINE PROGT2(IM,KM,RHSCNPY, &
     & RHSMC,AI,BI,CI,SMC,SLIMSK, &
     & CANOPY,PRECIP,RUNOFF,SNOWMT, &
     & ZSOIL,SOILTYP,SIGMAF,DELT,me)
      USE MODULE_GFS_MACHINE , ONLY : kind_phys
      implicit none
      integer km, IM, me
      real(kind=kind_phys) delt
      real(kind=kind_phys) RHSCNPY(IM), RHSMC(IM,KM), AI(IM,KM), &
     & BI(IM,KM), CI(IM,KM), SMC(IM,KM), &
     & SLIMSK(IM), CANOPY(IM), PRECIP(IM), &
     & RUNOFF(IM), SNOWMT(IM), ZSOIL(IM,KM), &
     & SIGMAF(IM)
      INTEGER SOILTYP(IM)
      integer k, lond, i
      real(kind=kind_phys) CNPY(IM), PRCP(IM), TSAT(IM), &
     & INF(IM), INFMAX(IM), SMSOIL(IM,KM)
      real(kind=kind_phys) cc, ctfil1, ctfil2, delt2, &
     & drip, rffact, rhoh2o, &
     & rzero, scanop, tdif, KSAT
      LOGICAL FLAG(IM)
      PARAMETER (SCANOP=.5, RHOH2O=1000.)
      PARAMETER (CTFIL1=.5, CTFIL2=1.-CTFIL1)
      PARAMETER (RFFACT=.15)
      LOND = 353
      DELT2 = DELT * 2.
      DO I=1,IM
        PRCP(I) = RHOH2O * (PRECIP(I)+SNOWMT(I)) / DELT
        RUNOFF(I) = 0.
        CNPY(I) = CANOPY(I)
      ENDDO
      DO I=1,IM
        IF(SLIMSK(I).EQ.1.) THEN
          RHSCNPY(I) = RHSCNPY(I) + SIGMAF(I) * PRCP(I)
          CANOPY(I) = CANOPY(I) + DELT * RHSCNPY(I)
          CANOPY(I) = MAX(CANOPY(I),0._kind_phys)
          PRCP(I) = PRCP(I) * (1. - SIGMAF(I))
          IF(CANOPY(I).GT.SCANOP) THEN
            DRIP = CANOPY(I) - SCANOP
            CANOPY(I) = SCANOP
            PRCP(I) = PRCP(I) + DRIP / DELT
          ENDIF
          INF(I) = PRCP(I)
          TSAT(I) = THSAT(SOILTYP(I))
          INFMAX(I) = (-ZSOIL(I,1)) * &
     & ((TSAT(I) - SMC(I,1)) / DELT - RHSMC(I,1)) &
     & * RHOH2O
          INFMAX(I) = MAX(RFFACT*INFMAX(I),0._kind_phys)
          IF(INF(I).GT.INFMAX(I)) THEN
            RUNOFF(I) = INF(I) - INFMAX(I)
            INF(I) = INFMAX(I)
          ENDIF
          INF(I) = INF(I) / RHOH2O
          RHSMC(I,1) = RHSMC(I,1) - INF(I) / ZSOIL(I,1)
        ENDIF
      ENDDO
      DO I=1,IM
        FLAG(I) = SLIMSK(I).EQ.1.
      ENDDO
      DO K = 1, KM
        DO I=1,IM
          IF(FLAG(I)) THEN
            RHSMC(I,K) = RHSMC(I,K) * DELT2
            AI(I,K) = AI(I,K) * DELT2
            BI(I,K) = 1. + BI(I,K) * DELT2
            CI(I,K) = CI(I,K) * DELT2
          ENDIF
        ENDDO
      ENDDO
      DO I=1,IM
        IF(FLAG(I)) THEN
          CI(I,1) = -CI(I,1) / BI(I,1)
          RHSMC(I,1) = RHSMC(I,1) / BI(I,1)
        ENDIF
      ENDDO
      DO K = 2, KM
        DO I=1,IM
          IF(FLAG(I)) THEN
            CC = 1. / (BI(I,K) + AI(I,K) * CI(I,K-1))
            CI(I,K) = -CI(I,K) * CC
            RHSMC(I,K)=(RHSMC(I,K)-AI(I,K)*RHSMC(I,K-1))*CC
          ENDIF
        ENDDO
      ENDDO
      DO I=1,IM
        IF(FLAG(I)) THEN
          CI(I,KM) = RHSMC(I,KM)
        ENDIF
      ENDDO
      DO K = KM-1, 1
        DO I=1,IM
          IF(FLAG(I)) THEN
            CI(I,K) = CI(I,K) * CI(I,K+1) + RHSMC(I,K)
          ENDIF
        ENDDO
      ENDDO
 100 CONTINUE
      DO K = 1, KM
        DO I=1,IM
          IF(FLAG(I)) THEN
            SMSOIL(I,K) = SMC(I,K) + CI(I,K)
            SMSOIL(I,K) = MAX(SMSOIL(I,K),0._kind_phys)
            TDIF = MAX(SMSOIL(I,K) - TSAT(I),0._kind_phys)
            RUNOFF(I) = RUNOFF(I) - &
     & RHOH2O * TDIF * ZSOIL(I,K) / DELT
            SMSOIL(I,K) = SMSOIL(I,K) - TDIF
          ENDIF
        ENDDO
      ENDDO
      DO K = 1, KM
        DO I=1,IM
          IF(FLAG(I)) THEN
            SMC(I,K) = CTFIL1 * SMSOIL(I,K) + CTFIL2 * SMC(I,K)
          ENDIF
        ENDDO
      ENDDO
      RETURN
      END SUBROUTINE PROGT2
      FUNCTION KTSOIL(THETA,KTYPE)
      USE MODULE_GFS_MACHINE , ONLY : kind_phys
      USE module_progtm , ONLY : TSAT, DFKT
      implicit none
      integer ktype,kw
      real(kind=kind_phys) ktsoil, theta, w
      W = (THETA / TSAT(KTYPE)) * 20. + 1.
      KW = W
      KW = MIN(KW,21)
      KW = MAX(KW,1)
      KTSOIL = DFKT(KW,KTYPE) &
     & + (W - KW) * (DFKT(KW+1,KTYPE) - DFKT(KW,KTYPE))
      RETURN
      END FUNCTION KTSOIL
      FUNCTION FUNCDF(THETA,KTYPE)
      USE MODULE_GFS_MACHINE , ONLY : kind_phys
      USE module_progtm , ONLY : TSAT, DFK
      implicit none
      integer ktype,kw
      real(kind=kind_phys) funcdf,theta,w
      W = (THETA / TSAT(KTYPE)) * 20. + 1.
      KW = W
      KW = MIN(KW,21)
      KW = MAX(KW,1)
      FUNCDF = DFK(KW,KTYPE) &
     & + (W - KW) * (DFK(KW+1,KTYPE) - DFK(KW,KTYPE))
      RETURN
      END FUNCTION FUNCDF
      FUNCTION FUNCKT(THETA,KTYPE)
      USE MODULE_GFS_MACHINE , ONLY : kind_phys
      USE module_progtm , ONLY : TSAT, KTK
      implicit none
      integer ktype,kw
      real(kind=kind_phys) funckt,theta,w
      W = (THETA / TSAT(KTYPE)) * 20. + 1.
      KW = W
      KW = MIN(KW,21)
      KW = MAX(KW,1)
      FUNCKT = KTK(KW,KTYPE) &
     & + (W - KW) * (KTK(KW+1,KTYPE) - KTK(KW,KTYPE))
      RETURN
      END FUNCTION FUNCKT
      FUNCTION THSAT(KTYPE)
      USE MODULE_GFS_MACHINE , ONLY : kind_phys
      USE module_progtm , ONLY : TSAT
      implicit none
      integer ktype
      real(kind=kind_phys) thsat
      THSAT = TSAT(KTYPE)
      RETURN
      END FUNCTION THSAT
      FUNCTION TWLT(KTYPE)
      USE MODULE_GFS_MACHINE , ONLY : kind_phys
      implicit none
      integer ktype
      real(kind=kind_phys) twlt
      TWLT = .1
      RETURN
      END FUNCTION TWLT
 END MODULE module_sf_gfs