!!WRF:MODEL_LAYER:PHYSICS
!
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
!-------------------------------------------------------------------
!-- U3D    	3D u-velocity interpolated to theta points (m/s)
!-- V3D    	3D v-velocity interpolated to theta points (m/s)
!-- T3D     	temperature (K)
!-- QV3D        3D water vapor mixing ratio (Kg/Kg)
!-- P3D         3D pressure (Pa)
!-- CP		heat capacity at constant pressure for dry air (J/kg/K)
!-- ROVCP       R/CP
!-- R           gas constant for dry air (J/kg/K)
!-- XLV         latent heat of vaporization for water (J/kg)
!-- PSFC	surface pressure (Pa)
!-- ZNT		roughness length (m)
!-- UST		u* in similarity theory (m/s)
!-- PSIM        similarity stability function for momentum
!-- PSIH        similarity stability function for heat
!-- XLAND	land mask (1 for land, 2 for water)
!-- HFX		upward heat flux at the surface (W/m^2)
!-- QFX   	upward moisture flux at the surface (kg/m^2/s)
!-- LH          net upward latent heat flux at surface (W/m^2)
!-- TSK		surface temperature (K)
!-- FLHC        exchange coefficient for heat (m/s)
!-- FLQC        exchange coefficient for moisture (m/s)
!-- QGH         lowest-level saturated mixing ratio
!-- U10         diagnostic 10m u wind
!-- V10         diagnostic 10m v wind
!-- GZ1OZ0      log(z/z0) where z0 is roughness length
!-- WSPD        wind speed at lowest model level (m/s)
!-- BR          bulk Richardson number in surface layer
!-- ISFFLX      isfflx=1 for surface heat and moisture fluxes
!-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
!-- KARMAN      Von Karman constant
!-- ids         start index for i in domain
!-- ide         end index for i in domain
!-- jds         start index for j in domain
!-- jde         end index for j in domain
!-- kds         start index for k in domain
!-- kde         end index for k in domain
!-- ims         start index for i in memory
!-- ime         end index for i in memory
!-- jms         start index for j in memory
!-- jme         end index for j in memory
!-- kms         start index for k in memory
!-- kme         end index for k in memory
!-- its         start index for i in tile
!-- ite         end index for i in tile
!-- jts         start index for j in tile
!-- jte         end index for j in tile
!-- kts         start index for k in tile
!-- kte         end index for k in tile
!-------------------------------------------------------------------

      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) ::            & 
                                        BR,                             &
                                        CHS,                            &
                                        CHS2,                           &
                                        CPM,                            &
                                        CQS2,                           &
                                        FLHC,                           &
                                        FLQC,                           &
                                        GZ1OZ0,                         &
                                        HFX,                            &
                                        LH,                             &
                                        PSIM,                           &
                                        PSIH,                           &
                                        QFX,                            &
                                        QGH,                            &
                                        QSFC,                           &
                                        UST,                            &
                                        ZNT,                            &
                                        WSPD

      REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) ::                 &
                                        U10,                            &
                                        V10 


!--------------------------- LOCAL VARS ------------------------------

      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)=QSFC(I,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,                                   &
!WRF              SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX,      &
!WRF              SLRAD,SNOWMT,DELT,                                    &
                  Z0RL,                                                 &
!WRF              TG3,GFLUX,F10M,                                       &
                  U10M,V10M,T2M,Q2M,                                    &
!WRF              ZSOIL,                                                &
                  CM,CH,RB,                                             &
!WRF              RHSCNPY,RHSMC,AIM,BIM,CIM,                            &
                  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,                             &
!WRF     &                  SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,   &
!WRF     &                  DLWFLX,SLRAD,SNOWMT,DELT,                   &
     &                  Z0RL,                                           &
!WRF     &                  TG3,GFLUX,F10M,                             &
     &                  U10M,V10M,T2M,Q2M,                              &
!WRF     &                  ZSOIL,                                      &
     &                  CM, CH, RB,                                     &
!WRF     &                  RHSCNPY,RHSMC,AIM,BIM,CIM,                  &
     &                  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
!
!     include 'constant.h'
!
      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)
!
!     Locals
!
      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,                   &
!WRF     &                     fhs,   funcdf,   funckt,g,      hl0, hl0inf, &
     &                     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)!C CA IS THE VON KARMAN CONSTANT
      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)
!WRF      real(kind=kind_phys) KT1(IM),       KT2(IM),      KTSOIL,         &
      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)

!
!  the 13 vegetation types are:
!
!  1  ...  broadleave-evergreen trees (tropical forest)
!  2  ...  broadleave-deciduous trees
!  3  ...  broadleave and needle leave trees (mixed forest)
!  4  ...  needleleave-evergreen trees
!  5  ...  needleleave-deciduous trees (larch)
!  6  ...  broadleave trees with groundcover (savanna)
!  7  ...  groundcover only (perenial)
!  8  ...  broadleave shrubs with perenial groundcover
!  9  ...  broadleave shrubs with bare soil
! 10  ...  dwarf trees and shrubs with ground cover (trunda)
! 11  ...  bare soil
! 12  ...  cultivations (use parameters from type 7)
! 13  ...  glacial
!
      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/
!
!!!   save rsmax, rsmin, rgl, hs, smmax, smdry, smref, smwlt
!

!WRF      DELT2 = DELT * 2.
!
!     ESTIMATE SIGMA ** K AT 2 M
!
      SIG2K = 1. - 4. * G * 2. / (CP * 280.)
!
!  INITIALIZE VARIABLES. ALL UNITS ARE SUPPOSEDLY M.K.S. UNLESS SPECIFIE
!  PSURF IS IN PASCALS
!  WIND IS WIND SPEED, THETA1 IS ADIABATIC SURFACE TEMP FROM LEVEL 1
!  RHO IS DENSITY, QS1 IS SAT. HUM. AT LEVEL1 AND QSS IS SAT. HUM. AT
!  SURFACE
!  CONVERT SLRAD TO THE CIVILIZED UNIT FROM LANGLEY MINUTE-1 K-4
!  SURFACE ROUGHNESS LENGTH IS CONVERTED TO M FROM CM
!
!!
!     qs1 = fpvs(t1)
!     qss = fpvs(tskin)
      DO I=1,IM
        XRCL(I)  = SQRT(RCL(I))
        PSURF(I) = 1000. * PS(I)
        PS1(I)   = 1000. * PRSL1(I)
!       SLWD(I)  = SLRAD(I) * CONVRAD
!WRF        SLWD(I)  = SLRAD(I)
!
!  DLWFLX has been given a negative sign for downward longwave
!  snet is the net shortwave flux
!
!WRF        SNET(I) = -SLWD(I) - DLWFLX(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))
!jfe    QS1(I) = 1000. * FPVS(T1(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))
!jfe    QSS(I) = 1000. * FPVS(TSURF(I))
        qss(i) = fpvs(tskin(i))
        QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
!       RS = PLANTR
        RS(I) = 0.
!WRF        if(VEGTYPE(I).gt.0.) RS(I) = rsmin(VEGTYPE(I))
        Z0(I) = .01 * Z0RL(i)
!WRF        CANOPY(I)= MAX(CANOPY(I),0._kind_phys)
        DM(I) = 1.
!WRF
      GOTO 1111
!WRF
        FACTSNW(I) = 10.
        IF(SLIMSK(I).EQ.2.) FACTSNW(I) = 3.
!
!  SNOW DEPTH IN WATER EQUIVALENT IS CONVERTED FROM MM TO M UNIT
!
        SNOWD(I) = SHELEG(I) / 1000.
        FLAGSNW(I) = .FALSE.
!
!  WHEN SNOW DEPTH IS LESS THAN 1 MM, A PATCHY SNOW IS ASSUMED AND
!  SOIL IS ALLOWED TO INTERACT WITH THE ATMOSPHERE.
!  WE SHOULD EVENTUALLY MOVE TO A LINEAR COMBINATION OF SOIL AND
!  SNOW UNDER THE CONDITION OF PATCHY SNOW.
!
        IF(SNOWD(I).GT..001.OR.SLIMSK(I).EQ.2.) RS(I) = 0.
        IF(SNOWD(I).GT..001) FLAGSNW(I) = .TRUE.
!##DG  IF(LAT.EQ.LATD) THEN
!##DG    PRINT *, ' WIND,TV1,TVS,Q1,QS1,SNOW,SLIMSK=',
!##DG&   WIND,TV1,TVS,Q1,QS1,SNOWD,SLIMSK
!##DG    PRINT *, ' SNET, SLWD =', SNET, SLWD(I)
!##DG  ENDIF
        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
!WRF
1111  CONTINUE
!WRF
      ENDDO

!!
!WRF
      GOTO 2222
!WRF
      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
!WRF
2222  CONTINUE
!WRF
!!
      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
!
!  COMPUTE STABILITY DEPENDENT EXCHANGE COEFFICIENTS
!
!  THIS PORTION OF THE CODE IS PRESENTLY SUPPRESSED
!
      DO I=1,IM
        IF(SLIMSK(I).EQ.0.) THEN
          USTAR(I) = SQRT(G * Z0(I) / CHARNOCK)
        ENDIF
!
!  COMPUTE STABILITY INDICES (RB AND HLINF)
!

        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)
!         RESTAR = ALOG(RESTAR)
!         RESTAR = MIN(RESTAR,5.)
!         RESTAR = MAX(RESTAR,-5.)
!         RAT(I) = AA1 + BB1 * RESTAR + CC1 * RESTAR ** 2
!         RAT(I) = RAT(I) / (1. + BB2 * RESTAR
!    &                       + CC2 * RESTAR ** 2)
!  Rat taken from Zeng, Zhao and Dickinson 1997
          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

!##DG  IF(LAT.EQ.LATD) THEN
!##DG    PRINT *, ' z0max, ztmax, restar, RAT(I) =', 
!##DG &   z0max, ztmax, restar, RAT(I)
!##DG  ENDIF
      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((Z0MAX(I)+Z1(I)) / Z0MAX(I))
!        FH(I) = LOG((ZTMAX(I)+Z1(I)) / ZTMAX(I))
        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
!##DG  IF(LAT.EQ.LATD) THEN
!##DG    PRINT *, ' DTV, RB(I), FM(I), FH(I), HLINF =',
!##DG &   dtv, rb, FM(I), FH(I), hlinf
!##DG  ENDIF
!
!  STABLE CASE
!
      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
!
!  SECOND ITERATION
!
      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)
!         AA = SQRT(1. + 4. * ALPHA * HL12(I))
          BB = SQRT(1. + 4. * ALPHA * HL12(I))
          PH2(I) = BB0 - BB + LOG((BB + 1.) / (BB0 + 1.))
        ENDIF
      ENDDO
!!
!##DG  IF(LAT.EQ.LATD) THEN
!##DG    PRINT *, ' HL1(I), PM, PH =',
!##DG &   HL1(I),  pm, ph
!##DG  ENDIF
!
!  UNSTABLE CASE
!
!
!  CHECK FOR UNPHYSICAL OBUKHOV LENGTH
!
      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
!
!  GET PM AND PH
!
      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
!
!  FINISH THE EXCHANGE COEFFICIENT COMPUTATION TO PROVIDE FM AND FH
!
      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))
!       USTAR(I) = SQRT(CM(I) * WIND(I) * WIND(I))
      ENDDO
!##DG  IF(LAT.EQ.LATD) THEN
!##DG    PRINT *, ' FM, FH, CM, CH(I), USTAR =',
!##DG &   FM, FH, CM, ch, USTAR
!##DG  ENDIF
!
!  UPDATE Z0 OVER OCEAN
!
      DO I = 1, IM
        IF(SLIMSK(I).EQ.0.) THEN
          Z0(I) = (CHARNOCK / G) * USTAR(I) ** 2
!  NEW IMPLEMENTATION OF Z0
!         CC = USTAR(I) * Z0 / RNU
!         PP = CC / (1. + CC)
!         FF = G * ARNU / (CHARNOCK * USTAR(I) ** 3)
!         Z0 = ARNU / (USTAR(I) * FF ** PP)
          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
!
!  RCP = RHO CP CH V
!
      DO I = 1, IM
        RCH(I) = RHO(I) * CP * CH(I) * WIND(I)
      ENDDO


!
!  SENSIBLE AND LATENT HEAT FLUX OVER OPEN WATER
!
      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

        !
        !  COMPUTE SOIL/SNOW/ICE HEAT FLUX IN PREPARATION FOR SURFACE ENERGY
        !  BALANCE CALCULATION
        !
	    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
        !  DF FOR ICE IS TAKEN FROM MAYKUT AND UNTERSTEINER
        !  DF IS IN SI UNIT OF W K-1 M-1
	        DFT0(I) = 2.2
	      ENDIF
	    ENDDO
        !!
	    DO I=1,IM
	      IF(SLIMSK(I).NE.0.) THEN
        !         IF(SNOWD(I).GT..001) THEN
	        IF(FLAGSNW(I)) THEN
        !
        !  WHEN SNOW COVERED, GROUND HEAT FLUX COMES FROM SNOW
        !
		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
        !
        !  COMPUTE POTENTIAL EVAPORATION FOR LAND AND SEA ICE
        !
	    DO I = 1, IM
	      IF(FLAG(I)) THEN
	        T12 = T1(I) * T1(I)
	        T14 = T12 * T12
        !
        !  RCAP = FNET - SIGMA T**4 + GFLX - RHO CP CH V (T1-THETA1)
        !
	        RCAP(I) = -SLWD(I) - SIGMA * T14 + GFLUX(I)                   &
	   &              - RCH(I) * (T1(I) - THETA1(I))
        !
        !  RSMALL = 4 SIGMA T**3 / RCH(I) + 1
        !
	        RSMALL(I) = 4. * SIGMA * T1(I) * T12 / RCH(I) + 1.
        !
        !  DELTA = L / CP * DQS/DT
        !
	        DELTA(I) = ELOCP * EPS * HVAP * QS1(I) / (RD * T12)
        !
        !  POTENTIAL EVAPOTRANSPIRATION ( WATTS / M**2 ) AND
        !  POTENTIAL EVAPORATION
        !
	        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
        !
        !  ACTUAL EVAPORATION OVER LAND IN THREE PARTS : EDIR, ET, AND EC
        !
        !  DIRECT EVAPORATION FROM SOIL, THE UNIT GOES FROM M S-1 TO KG M-2 S-1
        !
	    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 = .75 * THSAT(SOILTYP(I))
	        TREF(I) = smref(SOILTYP(I))
        !         TWILT = TWLT(SOILTYP(I))
	        TWILT(I) = smwlt(SOILTYP(I))
	        smcdry = smdry(SOILTYP(I))
        !         FX(I) = -2. * DF1(I) * (SMC(I,1) - .23) / ZSOIL(I,1)
        !    &            - KT1(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)
        !
        !  SIGMAF IS THE FRACTION OF AREA COVERED BY VEGETATION
        !
	        EDIR(I) = FX(I) * (1. - SIGMAF(I)) * PARTLND(I)
	      ENDIF
	    ENDDO
        !
        !  calculate stomatal resistance
        !
	    DO I = 1, IM
	      if(FLAG(I)) then
        !
        !  resistance due to PAR. We use net solar flux as proxy at the present time
        !
	        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.
        !
        !  resistance due to thermal effect
        !
        !         rct = 1. - .0016 * (topt - theta1) ** 2
        !         rct = max(rct,.0001)
        !
        !  resistance due to humidity
        !
        !         rcq = 1. / (1. + hs(VEGTYPE(I)) * (QS1(I) - Q0(I)))
        !         rcq = max(rcq,.0001)
        !
        !  compute resistance without the effect of soil moisture
        !
	        RS(I) = RS(I) / (rcs * rct * rcq)
	      endif
	    ENDDO
        !
        !  TRANSPIRATION FROM ALL LEVELS OF THE SOIL
        !
	    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)
        !
        !  resistance due to soil moisture deficit
        !
	        rss = GX(I) * (ZSOIL(I,1) / ZSOIL(I,km))
	        rss = max(rss,.0001_kind_phys)
	        RSI = RS(I) / rss
        !
        !  transpiration a la Monteith
        !
	        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)
        !
        !  resistance due to soil moisture deficit
        !
	        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
        !
        !  transpiration a la Monteith
        !
	        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
        !
        !  CANOPY RE-EVAPORATION
        !
	    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
        !
        !  SUM UP TOTAL EVAPORATION
        !
	    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
        !!
        !
        !  RETURN EVAP UNIT FROM KG M-2 S-1 TO WATTS M-2
        !
	    DO I=1,IM
	      IF(FLAG(I)) THEN
	        EVAP(I) = MIN(EVAP(I)*HVAP,EP(I))
	      ENDIF
	    ENDDO
        !##DG  IF(LAT.EQ.LATD) THEN
        !##DG    PRINT *, 'FX(I), SIGMAF, EDIR(I), ETPFAC=', FX(I)*HVAP,SIGMAF,
        !##DG&          EDIR(I)*HVAP,ETPFAC*HVAP
        !##DG    PRINT *, ' ET =', (ET(K)*HVAP,K=1,KM)
        !##DG    PRINT *, ' CANFAC(I), EC(I), EVAP', CANFAC(I),EC(I)*HVAP,EVAP
        !##DG  ENDIF
        !
        !  EVAPORATION OVER BARE SEA ICE
        !
	    DO I = 1, IM
        !       IF(SLIMSK(I).EQ.2.AND.SNOWD(I).LE..001) THEN
	      IF(SLIMSK(I).EQ.2.) THEN
	        EVAP(I) = PARTLND(I) * EP(I)
	      ENDIF
	    ENDDO
        !
        !  TREAT DOWNWARD MOISTURE FLUX SITUATION
        !  (EVAP WAS PRESET TO ZERO SO NO UPDATE NEEDED)
        !  DEW IS CONVERTED FROM KG M-2 TO M TO CONFORM TO PRECIP UNIT
        !
	    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
        !
        !  SNOW COVERED LAND AND SEA ICE
        !
	    DO I = 1, IM
	      FLAG(I) = SLIMSK(I).NE.0..AND.SNOWD(I).GT.0.
	    ENDDO
        !
        !  CHANGE OF SNOW DEPTH DUE TO EVAPORATION OR SUBLIMATION
        !
        !  CONVERT EVAP FROM KG M-2 S-1 TO M S-1 TO DETERMINE THE REDUCTION OF S
        !
	    DO I = 1, IM
	      IF(FLAG(I)) THEN
	        BFACT = SNOWD(I) / (DELT * EP(I) / (HVAP * RHOH2O))
	        BFACT = MIN(BFACT,1._kind_phys)
        !
        !  THE EVAPORATION OF SNOW
        !
	        IF(EP(I).LE.0.) BFACT = 1.
	        IF(SNOWD(I).LE..001) THEN
        !           EVAP = (SNOWD(I)/.001)*BFACT*EP(I) + EVAP
        !           SNOEVP(I) = bfact * EP(I) * (1. - PARTLND(I))
        !           EVAP = EVAP + SNOEVP(I)
		SNOEVP(I) = bfact * EP(I)
        !           EVAP   = EVAP + SNOEVP(I) * (1. - PARTLND(I))
		EVAP(I)=EVAP(I)+SNOEVP(I)*(1.-PARTLND(I))
	        ELSE
        !           EVAP(I) = BFACT * EP(I)
		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))           &
        !    &           + THETA1 - T1                                          &
        !    &           - BFACT * EP(I)) / (RSMALL(I) * RCH(I)                 &
	   &           - SNOEVP(I)) / (RSMALL(I) * RCH(I)                     &
	   &           + DFSNOW / (FACTSNW(I)* MAX(SNOWD(I),.001_kind_phys)))
        !         SNOWD(I) = SNOWD(I) - BFACT * EP(I) * DELT / (RHOH2O * HVAP)
	        SNOWD(I) = SNOWD(I) - SNOEVP(I) * delt / (rhoh2o * hvap)
	        SNOWD(I) = MAX(SNOWD(I),0._kind_phys)
	      ENDIF
	    ENDDO
        !
        !  SNOW MELT (M)
        !
         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
        !
        !  We need to re-evaluate evaporation because of snow melt
        !    the skin temperature is now bounded to 0 deg C
        !
        !     qss = fpvs(tsurf)
	    DO I = 1, IM
        !       IF (SNOWD(I) .GT. 0.0) THEN
	      IF (SNOWD(I) .GT. snomin) THEN
        !jfe      QSS(I) = 1000. * FPVS(TSURF(I))
	        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
        !
        !  PREPARE TENDENCY TERMS FOR THE SOIL MOISTURE FIELD WITHOUT PRECIPITAT
        !  THE UNIT OF MOISTURE FLUX NEEDS TO BECOME M S-1 FOR SOIL MOISTURE
        !   HENCE THE FACTOR OF RHOH2O
        !
	    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, BIM, AND CIM ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
        !  IMPLICIT UPDATE OF THE SOIL MOISTURE
        !
	        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
        !
        !  UPDATE SOIL TEMPERATURE AND SEA ICE TEMPERATURE
        !
	    DO I=1,IM
	      FLAG(I) = SLIMSK(I).NE.0.
	    ENDDO
        !
        !  SURFACE TEMPERATURE IS PART OF THE UPDATE WHEN SNOW IS ABSENT
        !
	    DO I=1,IM
        !       IF(FLAG(I).AND.SNOWD(I).LE..001) THEN
	      IF(FLAG(I).AND..NOT.FLAGSNW(I)) THEN
	        YY(I) = T1(I) +                                               &
        !    &          (RCAP(I)-GFLUX(I) + THETA1 - 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.SNOWD(I).GT..001) THEN
	      IF(FLAG(I).AND.FLAGSNW(I)) THEN
	        YY(I) = STSOIL(I,1)
        !
        !  HEAT FLUX FROM SNOW IS EXPLICIT IN TIME
        !
	        ZZ(I) = 1.
	        XX(I) = DFSNOW * (STSOIL(I,1) - TSURF(I))                     &
	   &            / (-FACTSNW(I) * MAX(SNOWD(I),.001_kind_phys))
	      ENDIF
	    ENDDO
        !
        !  COMPUTE THE FORCING AND THE IMPLICIT MATRIX ELEMENTS FOR UPDATE
        !
        !  CH2O IS THE HEAT CAPACITY OF WATER AND CSOIL IS THE HEAT CAPACITY OF
        !
	    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, BI, AND CI ARE THE ELEMENTS OF THE TRIDIAGONAL MATRIX FOR THE
        !  IMPLICIT UPDATE OF THE SOIL TEMPERATURE
        !
	        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))  
        !         SS = DFT0(I) * (STSOIL(I,1) - YY(I))                          &
        !    &         / (.5 * ZSOIL(I,1) * ZZ(I))
        !         RHSTC(1) = (DFT1(I) * DTDZ1(I) - SS)
	        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
        !
        !  AT THE BOTTOM, CLIMATOLOGY IS ASSUMED AT 2M DEPTH FOR LAND AND
        !  FREEZING TEMPERATURE IS ASSUMED FOR SEA ICE AT Z(KM)
	        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
        !
        !  SOLVE THE TRI-DIAGONAL MATRIX
        !
	    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
        !  FORWARD ELIMINATION
	    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
        !!
        !  BACKWARD SUBSTITUTTION
	    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
        !
        !  UPDATE SOIL AND ICE TEMPERATURE
        !
	    DO K = 1, KM
	      DO I=1,IM
	        IF(FLAG(I)) THEN
		STSOIL(I,K) = STSOIL(I,K) + CI(I,K)
	        ENDIF
	      ENDDO
	    ENDDO
        !
        !  UPDATE SURFACE TEMPERATURE FOR SNOW FREE SURFACES
        !
	    DO I=1,IM
        !       IF(SLIMSK(I).NE.0..AND.SNOWD(I).LE..001) THEN
	      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.SNOWD(I).LE..001) THEN
	      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
        !
        !  TIME FILTER FOR SOIL AND SKIN TEMPERATURE
        !
	    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
        !
        !  GFLUX CALCULATION
        !
	    DO I=1,IM
	      FLAG(I) = SLIMSK(I).NE.0.                                       &
        !    &            .AND.SNOWD(I).GT..001                                 &
	   &            .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.SNOWD(I).LE..001) THEN
	      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

        !
        !  CALCULATE SENSIBLE HEAT FLUX
        !
!WRF	    DO I = 1, IM
!WRF      HFLX(I) = RCH(I) * (TSKIN(I) - THETA1(I))
!WRF	    ENDDO
        !
        !  THE REST OF THE OUTPUT
        !
!WRF	    DO I = 1, IM
!WRF	      QSURF(I) = Q1(I) + EVAP(I) / (ELOCP * RCH(I))
!WRF	      DM(I) = 1.
        !
        !  CONVERT SNOW DEPTH BACK TO MM OF WATER EQUIVALENT
        !
!WRF	      SHELEG(I) = SNOWD(I) * 1000.
!WRF	    ENDDO
        !

      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)
!WRF         T2M(I) = TSKIN(I) * (1. - FH2(I) / FH(I))                      &
!WRF     &          + THETA1(I) * FH2(I) / FH(I)
!WRF         T2M(I) = T2M(I) * SIG2K
!        Q2M(I) = QSURF(I) * (1. - FH2(I) / FH(I))                      &
!    &         + Q1(I) * FH2(I) / FH(I)
!       T2M(I) = T1
!       Q2M(I) = Q1
!WRF        IF(EVAP(I).GE.0.) THEN
!
!  IN CASE OF EVAPORATION, USE THE INFERRED QSURF TO DEDUCE Q2M
!
!WRF          Q2M(I) = QSURF(I) * (1. - FH2(I) / FH(I))                     &
!WRF     &         + Q1(I) * FH2(I) / FH(I)
!WRF        ELSE
!
!  FOR DEW FORMATION SITUATION, USE SATURATED Q AT TSKIN
!
!jfe      QSS(I) = 1000. * FPVS(TSKIN(I))
!WRF          qss(I) = fpvs(tskin(I))
!WRF          QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
!WRF          Q2M(I) = QSS(I) * (1. - FH2(I) / FH(I))                       &
!WRF     &         + Q1(I) * FH2(I) / FH(I)
!WRF        ENDIF
!jfe    QSS(I) = 1000. * FPVS(T2M(I))
!WRF        QSS(I) = fpvs(t2m(I))
!       QSS(I) = 1000. * T2MO(I)
!WRF        QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I))
!WRF        Q2M(I) = MIN(Q2M(I),QSS(I))
      ENDDO
!!
!     DO I = 1, IM
!       RNET(I) = -SLWD(I) - SIGMA * TSKIN(I) **4
!     ENDDO
!!
!
!WRF      do i=1,im
!WRF        tem     = 1.0 / rho(i)
!WRF        hflx(i) = hflx(i) * tem * cpinv
!WRF        evap(i) = evap(i) * tem * hvapi
!WRF      enddo


!
!##DG  IF(LAT.EQ.LATD) THEN
!C       RBAL = -SLWD-SIGMA*TSKIN**4+GFLUX
!C    &         -EVAP - HFLX
!##DG    PRINT 6000,HFLX,EVAP,GFLUX,
!##DG&             STC(1), STC(2),TSKIN,RNET,SLWD
!##DG    PRINT *, ' T1 =', T1
 6000 FORMAT(8(F8.2,','))
!C     PRINT *, ' EP, ETP,T2M(I) =', EP, ETP,T2M(I)
!C     PRINT *, ' FH, FH2 =', FH, FH2
!C     PRINT *, ' PH, PH2 =', PH, PH2
!C     PRINT *, ' CH, RCH =', CH, RCH
!C     PRINT *, ' TERM1, TERM2 =', TERM1, TERM2
!C     PRINT *, ' RS(I), PLANTR =', RS(I), PLANTR
!##DG  ENDIF

      RETURN
      END SUBROUTINE PROGTM
!
!  PROGT2 IS THE SECOND PART OF THE SOIL MODEL THAT IS EXECUTED
!  AFTER PRECIPITATION FOR THE TIME STEP HAS BEEN CALCULATED
!
!FPP$ NOCONCUR R
!FPP$ EXPAND(FUNCDF,FUNCKT,THSAT)
      SUBROUTINE PROGT2(IM,KM,RHSCNPY,                                  &
     &                  RHSMC,AI,BI,CI,SMC,SLIMSK,                      &
     &                  CANOPY,PRECIP,RUNOFF,SNOWMT,                    &
     &                  ZSOIL,SOILTYP,SIGMAF,DELT,me)
!c
      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,                       &
!WRF     &                     rzero, scanop, tdif, thsat, KSAT
     &                     rzero, scanop, tdif, KSAT
!
      LOGICAL FLAG(IM)
!c
      PARAMETER (SCANOP=.5, RHOH2O=1000.)
      PARAMETER (CTFIL1=.5, CTFIL2=1.-CTFIL1)
!     PARAMETER (CTFIL1=1., CTFIL2=1.-CTFIL1)
      PARAMETER (RFFACT=.15)
!
!##DG LATD = 44
      LOND = 353
      DELT2 = DELT * 2.
!
!  PRECIPITATION RATE IS NEEDED IN UNIT OF KG M-2 S-1
!
      DO I=1,IM
        PRCP(I) = RHOH2O * (PRECIP(I)+SNOWMT(I)) / DELT
        RUNOFF(I) = 0.
        CNPY(I) = CANOPY(I)
      ENDDO
!##DG  IF(LAT.EQ.LATD) THEN
!##DG    PRINT *, ' BEFORE RUNOFF CAL, RHSMC =', RHSMC(1)
!##DG  ENDIF
!
!  UPDATE CANOPY WATER CONTENT
!
      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
!
!  CALCULATE INFILTRATION RATE
!
          INF(I) = PRCP(I)
          TSAT(I) = THSAT(SOILTYP(I))
!         DSAT = FUNCDF(TSAT(I),SOILTYP(I))
!         KSAT = FUNCKT(TSAT(I),SOILTYP(I))
!         INFMAX(I) = -DSAT * (TSAT(I) - SMC(I,1))
!    &                / (.5 * ZSOIL(I,1))                               &
!    &                + KSAT 
          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(SMC(I,1).GE.TSAT(I)) INFMAX(I) = KSAT
!         IF(SMC(I,1).GE.TSAT(I)) INFMAX(I) = ZSOIL(I,1) * RHSMC(I,1)
          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
!!
!##DG  IF(LAT.EQ.LATD) THEN
!##DG    PRINT *, ' PRCP(I), INFMAX(I), RUNOFF =', PRCP(I),INFMAX(I),RUNOFF
!##DG    PRINT *, ' SMSOIL =', SMC(1), SMC(2)
!##DG    PRINT *, ' RHSMC =', RHSMC(1)
!##DG  ENDIF
!
!  WE CURRENTLY IGNORE THE EFFECT OF RAIN ON SEA ICE
!
      DO I=1,IM
        FLAG(I) = SLIMSK(I).EQ.1.
      ENDDO
!!
!
!  SOLVE THE TRI-DIAGONAL MATRIX
!
      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
!  FORWARD ELIMINATION
      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
!  BACKWARD SUBSTITUTTION
      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
!
!  UPDATE SOIL MOISTURE
!
      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
!       IF(FLAG(I)) THEN
!         CANOPY(I) = CTFIL1 * CANOPY(I) + CTFIL2 * CNPY(I)
!       ENDIF
!     I = 1
!     PRINT *, ' SMC'
!     PRINT 6000, SMC(1), SMC(2)
!6000 FORMAT(2(F8.5,','))
      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
!     USE module_progtm
      implicit none
      integer              ktype
      real(kind=kind_phys) twlt
!
      TWLT = .1
      RETURN
      END FUNCTION TWLT

 END MODULE module_sf_gfs
