sflx.f

 
 
      module sflx
      contains
      subroutine gfssflx                                                &!  ---  inputs:
     &     ( nsoil, couple, icein, ffrozp, dt, zlvl, sldpth,            &
     &       swdn, swnet, lwdn, sfcems, sfcprs, sfctmp,                 &
     &       sfcspd, prcp, q2, q2sat, dqsdt2, th2, ivegsrc,             &
     &       vegtyp, soiltyp, slopetyp, shdmin, alb, snoalb,            &
     &       rhonewsn, exticeden,                                       &
     &       bexpp, xlaip,                                              & !  sfc-perts, mgehne
     &       lheatstrg,                                                 &!  ---  input/outputs:
     &       tbot, cmc, t1, stc, smc, sh2o, sneqv, ch, cm,z0,           &!  ---  outputs:
     &       nroot, shdfac, snowh, albedo, eta, sheat, ec,              &
     &       edir, et, ett, esnow, drip, dew, beta, etp, ssoil,         &
     &       flx1, flx2, flx3, runoff1, runoff2, runoff3,               &
     &       snomlt, sncovr, rc, pc, rsmin, xlai, rcs, rct, rcq,        &
     &       rcsoil, soilw, soilm, smcwlt, smcdry, smcref, smcmax,      &
     &       errmsg, errflg )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine sflx - version 2.7:                                       !
!  sub-driver for "noah/osu lsm" family of physics subroutines for a    !
!  soil/veg/snowpack land-surface model to update soil moisture, soil   !
!  ice, soil temperature, skin temperature, snowpack water content,     !
!  snowdepth, and all terms of the surface energy balance and surface   !
!  water balance (excluding input atmospheric forcings of downward      !
!  radiation and precip)                                                !
!                                                                       !
!  usage:                                                               !
!                                                                       !
!      call sflx                                                        !
!  ---  inputs:                                                         !
!          ( nsoil, couple, icein, ffrozp, dt, zlvl, sldpth,            !
!            swdn, swnet, lwdn, sfcems, sfcprs, sfctmp,                 !
!            sfcspd, prcp, q2, q2sat, dqsdt2, th2,ivegsrc,              !
!            vegtyp, soiltyp, slopetyp, shdmin, alb, snoalb,            !
!  ---  input/outputs:                                                  !
!            tbot, cmc, t1, stc, smc, sh2o, sneqv, ch, cm,              !
!  ---  outputs:                                                        !
!            nroot, shdfac, snowh, albedo, eta, sheat, ec,              !
!            edir, et, ett, esnow, drip, dew, beta, etp, ssoil,         !
!            flx1, flx2, flx3, runoff1, runoff2, runoff3,               !
!            snomlt, sncovr, rc, pc, rsmin, xlai, rcs, rct, rcq,        !
!            rcsoil, soilw, soilm, smcwlt, smcdry, smcref, smcmax )     !
!                                                                       !
!                                                                       !
!  subprograms called:  redprm, snow_new, csnow, snfrac, alcalc,        !
!            tdfcnd, snowz0, sfcdif, penman, canres, nopac, snopac.     !
!                                                                       !
!                                                                       !
!  program history log:                                                 !
!    jun  2003  -- k. mitchell et. al -- created version 2.7            !
!         200x  -- sarah lu    modified the code including:             !
!                       added passing argument, couple; replaced soldn  !
!                       and solnet by radflx; call sfcdif if couple=0;  !
!                       apply time filter to stc and tskin; and the     !
!                       way of namelist inport.                         !
!    feb  2004 -- m. ek noah v2.7.1 non-linear weighting of snow vs     !
!                       non-snow covered portions of gridbox            !
!    apr  2009  -- y.-t. hou   added lw surface emissivity effect,      !
!                       streamlined and reformatted the code, and       !
!                       consolidated constents/parameters by using      !
!                       module physcons, and added program documentation!
!    sep  2009 -- s. moorthi minor fixes                                !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers (>=2 but <=nsold)  1    !
!     couple   - integer, =0:uncoupled (land model only)           1    !
!                         =1:coupled with parent atmos model            !
!     icein    - integer, sea-ice flag (=1: sea-ice, =0: land)     1    !
!     ffrozp   - real, fractional snow/rain                        1    !
!     dt       - real, time step (<3600 sec)                       1    !
!     zlvl     - real, height abv atmos ground forcing vars (m)    1    !
!     sldpth   - real, thickness of each soil layer (m)          nsoil  !
!     swdn     - real, downward sw radiation flux (w/m**2)         1    !
!     swnet    - real, downward sw net (dn-up) flux (w/m**2)       1    !
!     lwdn     - real, downward lw radiation flux (w/m**2)         1    !
!     sfcems   - real, sfc lw emissivity (fractional)              1    !
!     sfcprs   - real, pressure at height zlvl abv ground(pascals) 1    !
!     sfctmp   - real, air temp at height zlvl abv ground (k)      1    !
!     sfcspd   - real, wind speed at height zlvl abv ground (m/s)  1    !
!     prcp     - real, precip rate (kg m-2 s-1)                    1    !
!     q2       - real, mixing ratio at hght zlvl abv grnd (kg/kg)  1    !
!     q2sat    - real, sat mixing ratio at zlvl abv grnd (kg/kg)   1    !
!     dqsdt2   - real, slope of sat specific humidity curve at     1    !
!                      t=sfctmp (kg kg-1 k-1)                           !
!     th2      - real, air potential temp at zlvl abv grnd (k)     1    !
!     ivegsrc  - integer, sfc veg type data source umd or igbp          !
!     vegtyp   - integer, vegetation type (integer index)          1    !
!     soiltyp  - integer, soil type (integer index)                1    !
!     slopetyp - integer, class of sfc slope (integer index)       1    !
!     shdmin   - real, min areal coverage of green veg (fraction)  1    !
!     alb      - real, bkground snow-free sfc albedo (fraction)    1    !
!     snoalb   - real, max albedo over deep snow     (fraction)    1    !
!     lheatstrg- logical, flag for canopy heat storage             1    !
!                         parameterization                              !
!                                                                       !
!  input/outputs:                                                       !
!     tbot     - real, bottom soil temp (k)                        1    !
!                      (local yearly-mean sfc air temp)                 !
!     cmc      - real, canopy moisture content (m)                 1    !
!     t1       - real, ground/canopy/snowpack eff skin temp (k)    1    !
!     stc      - real, soil temp (k)                             nsoil  !
!     smc      - real, total soil moisture (vol fraction)        nsoil  !
!     sh2o     - real, unfrozen soil moisture (vol fraction)     nsoil  !
!                      note: frozen part = smc-sh2o                     !
!     sneqv    - real, water-equivalent snow depth (m)             1    !
!                      note: snow density = snwqv/snowh                 !
!     ch       - real, sfc exchange coeff for heat & moisture (m/s)1    !
!                      note: conductance since it's been mult by wind   !
!     cm       - real, sfc exchange coeff for momentum (m/s)       1    !
!                      note: conductance since it's been mult by wind   !
!                                                                       !
!  outputs:                                                             !
!     nroot    - integer, number of root layers                    1    !
!     shdfac   - real, aeral coverage of green veg (fraction)      1    !
!     snowh    - real, snow depth (m)                              1    !
!     albedo   - real, sfc albedo incl snow effect (fraction)      1    !
!     eta      - real, downward latent heat flux (w/m2)            1    !
!     sheat    - real, downward sensible heat flux (w/m2)          1    !
!     ec       - real, canopy water evaporation (w/m2)             1    !
!     edir     - real, direct soil evaporation (w/m2)              1    !
!     et       - real, plant transpiration     (w/m2)            nsoil  !
!     ett      - real, total plant transpiration (w/m2)            1    !
!     esnow    - real, sublimation from snowpack (w/m2)            1    !
!     drip     - real, through-fall of precip and/or dew in excess 1    !
!                      of canopy water-holding capacity (m)             !
!     dew      - real, dewfall (or frostfall for t<273.15) (m)     1    !
!     beta     - real, ratio of actual/potential evap              1    !
!     etp      - real, potential evaporation (w/m2)                1    !
!     ssoil    - real, upward soil heat flux (w/m2)                1    !
!     flx1     - real, precip-snow sfc flux  (w/m2)                1    !
!     flx2     - real, freezing rain latent heat flux (w/m2)       1    !
!     flx3     - real, phase-change heat flux from snowmelt (w/m2) 1    !
!     snomlt   - real, snow melt (m) (water equivalent)            1    !
!     sncovr   - real, fractional snow cover                       1    !
!     runoff1  - real, surface runoff (m/s) not infiltrating sfc   1    !
!     runoff2  - real, sub sfc runoff (m/s) (baseflow)             1    !
!     runoff3  - real, excess of porosity for a given soil layer   1    !
!     rc       - real, canopy resistance (s/m)                     1    !
!     pc       - real, plant coeff (fraction) where pc*etp=transpi 1    !
!     rsmin    - real, minimum canopy resistance (s/m)             1    !
!     xlai     - real, leaf area index  (dimensionless)            1    !
!     rcs      - real, incoming solar rc factor  (dimensionless)   1    !
!     rct      - real, air temp rc factor        (dimensionless)   1    !
!     rcq      - real, atoms vapor press deficit rc factor         1    !
!     rcsoil   - real, soil moisture rc factor   (dimensionless)   1    !
!     soilw    - real, available soil mois in root zone            1    !
!     soilm    - real, total soil column mois (frozen+unfrozen) (m)1    !
!     smcwlt   - real, wilting point (volumetric)                  1    !
!     smcdry   - real, dry soil mois threshold (volumetric)        1    !
!     smcref   - real, soil mois threshold     (volumetric)        1    !
!     smcmax   - real, porosity (sat val of soil mois)             1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
      use machine ,   only : kind_phys
!
      use physcons,   only : con_cp, con_rd, con_t0c, con_g, con_pi,    &
     &                       con_cliq, con_csol, con_hvap, con_hfus,    &
     &                       con_sbc
!
      implicit none
 
!  ---  constant parameters:
!      *** note: some of the constants are different in subprograms and need to
!          be consolidated with the standard def in module physcons at sometime
!          at the present time, those diverse values are kept temperately to
!          provide the same result as the original codes.  -- y.t.h.  may09
 
      integer,               parameter :: nsold   = 4
 
!     real (kind=kind_phys), parameter :: gs      = con_g       !< con_g   =9.80665
      real (kind=kind_phys), parameter :: gs1     = 9.8         
      real (kind=kind_phys), parameter :: gs2     = 9.81        
      real (kind=kind_phys), parameter :: tfreez  = con_t0c     
      real (kind=kind_phys), parameter :: lsubc   = 2.501e+6    
      real (kind=kind_phys), parameter :: lsubf   = 3.335e5     
      real (kind=kind_phys), parameter :: lsubs   = 2.83e+6     ! ? in sflx, snopac
      real (kind=kind_phys), parameter :: elcp    = 2.4888e+3   ! ? in penman
!     real (kind=kind_phys), parameter :: rd      = con_rd      ! con_rd  =287.05
      real (kind=kind_phys), parameter :: rd1     = 287.04      ! con_rd in sflx, penman, canres
      real (kind=kind_phys), parameter :: cp      = con_cp      ! con_cp  =1004.6
      real (kind=kind_phys), parameter :: cp1     = 1004.5      ! con_cp in sflx, canres
      real (kind=kind_phys), parameter :: cp2     = 1004.0      ! con_cp in htr
!     real (kind=kind_phys), parameter :: cph2o   = con_cliq    ! con_cliq=4.1855e+3
      real (kind=kind_phys), parameter :: cph2o1  = 4.218e+3    ! con_cliq in penman, snopac
      real (kind=kind_phys), parameter :: cph2o2  = 4.2e6       ! con_cliq in hrt *unit diff!
      real (kind=kind_phys), parameter :: cpice   = con_csol    ! con_csol=2.106e+3
      real (kind=kind_phys), parameter :: cpice1  = 2.106e6     ! con_csol in hrt *unit diff!
!     real (kind=kind_phys), parameter :: sigma   = con_sbc     ! con_sbc=5.6704e-8
      real (kind=kind_phys), parameter :: sigma1  = 5.67e-8     ! con_sbc in penman, nopac, snopac
 
!  ---  inputs:
      integer, intent(in) :: nsoil, couple, icein, vegtyp, soiltyp,     &
     &       slopetyp, ivegsrc
 
      real (kind=kind_phys), intent(in) :: ffrozp, dt, zlvl, lwdn,      &
     &       sldpth(nsoil), swdn, swnet, sfcems, sfcprs, sfctmp,        &
     &       sfcspd, prcp, q2, q2sat, dqsdt2, th2, shdmin, alb, snoalb, &
     &       bexpp, xlaip, rhonewsn                                     & !sfc-perts, mgehne
 
      logical, intent(in) :: lheatstrg, exticeden
 
!  ---  input/outputs:
      real (kind=kind_phys), intent(inout) :: tbot, cmc, t1, sneqv,     &
     &       stc(nsoil), smc(nsoil), sh2o(nsoil), ch, cm
 
!  ---  outputs:
      integer, intent(out) :: nroot
 
      real (kind=kind_phys), intent(out) :: shdfac, snowh, albedo,      &
     &       eta, sheat, ec, edir, et(nsoil), ett, esnow, drip, dew,    &
     &       beta, etp, ssoil, flx1, flx2, flx3, snomlt, sncovr,        &
     &       runoff1, runoff2, runoff3, rc, pc, rsmin, xlai, rcs,       &
     &       rct, rcq, rcsoil, soilw, soilm, smcwlt, smcdry, smcref,    &
     &       smcmax
      character(len=*),     intent(out) :: errmsg
      integer,              intent(out) :: errflg
 
!  ---  locals:
!     real (kind=kind_phys) ::  df1h,
      real (kind=kind_phys) ::  bexp, cfactr, cmcmax, csoil, czil,      &
     &       df1, df1a, dksat, dwsat, dsoil, dtot, frcsno,              &
     &       frcsoi, epsca, fdown, f1, fxexp, frzx, hs, kdt, prcp1,     &
     &       psisat, quartz, rch, refkdt, rr, rgl, rsmax, sndens,       &
     &       sncond, sbeta, sn_new, slope, snup, salp, soilwm, soilww,  &
     &       t1v, t24, t2v, th2v, topt, tsnow, zbot, z0
      
      real (kind=kind_phys) ::  shdfac0
      real (kind=kind_phys), dimension(nsold) :: rtdis, zsoil
 
      logical :: frzgra, snowng
 
      integer :: ice, k, kz
!
!===> ...  begin here
!
! Initialize CCPP error-handling
      errflg = 0
      errmsg = ''
 
!  --- ...  initialization
 
      runoff1 = 0.0
      runoff2 = 0.0
      runoff3 = 0.0
      snomlt  = 0.0
      rc      = 0.0
 
!  --- ...  define local variable ice to achieve:
!             sea-ice case,          ice =  1
!             non-glacial land,      ice =  0
!             glacial-ice land,      ice = -1
!             if vegtype=15 (glacial-ice), re-set ice flag = -1 (glacial-ice)
!    note - for open-sea, sflx should *not* have been called. set green
!           vegetation fraction (shdfac) = 0.
 
      shdfac0 = shdfac
      ice = icein
 
      if(ivegsrc == 2) then
       if (vegtyp == 13) then
        ice = -1
        shdfac = 0.0
       endif
      endif
 
      if(ivegsrc == 1) then
       if (vegtyp == 15) then
        ice = -1
        shdfac = 0.0
       endif
      endif
 
      if (ice == 1) then
 
        shdfac = 0.0
 
 
        do kz = 1, nsoil
          zsoil(kz) = -3.0 * float(kz) / float(nsoil)
        enddo
 
      else
 
!    note - sign of zsoil is negative (denoting below ground)
 
        zsoil(1) = -sldpth(1)
        do kz = 2, nsoil
          zsoil(kz) = -sldpth(kz) + zsoil(kz-1)
        end do
 
      endif   ! end if_ice_block
 
!  --- ...  next is crucial call to set the land-surface parameters,
!           including soil-type and veg-type dependent parameters.
!           set shdfac=0.0 for bare soil surfaces
 
      call redprm(errmsg, errflg)
      if(errflg/=0) return
 
        if(ivegsrc == 1) then
!only igbp type has urban
!urban
         if(vegtyp == 13)then
!             shdfac=0.05
!             rsmin=400.0
!             smcmax = 0.45
!             smcref = 0.42
!             smcwlt = 0.40
!             smcdry = 0.40
              rsmin=400.0*(1-shdfac0)+40.0*shdfac0   ! gvf
              shdfac=shdfac0                         ! gvf
              smcmax = 0.45*(1-shdfac0)+smcmax*shdfac0
              smcref = 0.42*(1-shdfac0)+smcref*shdfac0
              smcwlt = 0.40*(1-shdfac0)+smcwlt*shdfac0
              smcdry = 0.40*(1-shdfac0)+smcdry*shdfac0
         endif
        endif
 
!  ---  inputs:                                                            !
!          ( nsoil, vegtyp, soiltyp, slopetyp, sldpth, zsoil,              !
!  ---  outputs:                                                           !
!            cfactr, cmcmax, rsmin, rsmax, topt, refkdt, kdt,              !
!            sbeta, shdfac, rgl, hs, zbot, frzx, psisat, slope,            !
!            snup, salp, bexp, dksat, dwsat, smcmax, smcwlt,               !
!            smcref, smcdry, f1, quartz, fxexp, rtdis, nroot,              !
!            z0, czil, xlai, csoil )                                       !
 
!  --- ...  bexp sfc-perts, mgehne
 
      if( bexpp < 0.) then
         bexp = bexp * max(1.+bexpp, 0.)
      endif
      if( bexpp >= 0.) then
         bexp = bexp * min(1.+bexpp, 2.)
      endif
!  --- ...  lai sfc-perts, mgehne
      xlai = xlai * (1.+xlaip)
      xlai = max(xlai, .75)
 
 
      snowng = .false.
      frzgra = .false.
 
!    note - runoff2 is then a negative value (as a flag) over sea-ice or
!           glacial-ice, in order to achieve water balance.
 
      if (ice == 1) then
 
        if (sneqv < 0.01) then
          sneqv = 0.01
          snowh = 0.10
!         snowh = sneqv / sndens
        endif
 
      elseif (ice == -1) then
 
        if (sneqv < 0.10) then
!         sndens = sneqv / snowh
!         runoff2 = -(0.10 - sneqv) / dt
          sneqv = 0.10
          snowh = 1.00
!         snowh = sneqv / sndens
        endif
 
      endif   ! end if_ice_block
 
 
      if (ice /= 0) then
        do kz = 1, nsoil
          smc(kz) = 1.0
          sh2o(kz) = 1.0
        enddo
      endif
 
! (note that csnow is a function subroutine)
 
      if (sneqv .eq. 0.0) then
        sndens = 0.0
        snowh = 0.0
        sncond = 1.0
      else
        sndens = sneqv / snowh
        sndens = max( 0.0, min( 1.0, sndens ))   ! added by moorthi
 
        call csnow
!  ---  inputs:                                                         !
!          ( sndens,                                                    !
!  ---  outputs:                                                        !
!            sncond )                                                   !
 
      endif
 
 
      if (prcp > 0.0) then
        if (ffrozp > 0.) then
          snowng = .true.
        else
          if (t1 <= tfreez) frzgra = .true.
        endif
      endif
 
! determine new snowfall (converting precipitation rate from
! \f$kg m^{-2} s^{-1}\f$ to a liquid equiv snow depth in meters)
!  and add it to the existing snowpack.
 
      if (snowng .or. frzgra) then
 
!   snowfall
       if (snowng) then
        sn_new = ffrozp*prcp * dt * 0.001
        sneqv = sneqv + sn_new
        prcp1 = (1.-ffrozp)*prcp
       endif
!    freezing rain
       if (frzgra) then
        sn_new = prcp * dt * 0.001
        sneqv = sneqv + sn_new
        prcp1 = 0.0
       endif
 
        call snow_new
!  ---  inputs:                                                         !
!          ( sfctmp, sn_new, rhonewsn, exticeden,                       !
!  ---  input/outputs:                                                  !
!            snowh, sndens )                                            !
 
        call csnow
!  ---  inputs:                                                         !
!          ( sndens,                                                    !
!  ---  outputs:                                                        !
!            sncond )                                                   !
 
      else
 
 
        prcp1 = prcp
 
      endif   ! end if_snowng_block
 
 
      if (ice /= 0) then
 
!  --- ...  snow cover, albedo over sea-ice, glacial-ice
 
        sncovr = 1.0
        albedo = 0.65
 
      else
 
!  --- ...  non-glacial land
!           if snow depth=0, set snowcover fraction=0, albedo=snow free albedo.
 
        if (sneqv == 0.0) then
 
          sncovr = 0.0
          albedo = alb
 
        else
 
!  --- ...  determine snow fraction cover.
!           determine surface albedo modification due to snowdepth state.
          call snfrac
!  ---  inputs:                                                         !
!          ( sneqv, snup, salp, snowh,                                  !
!  ---  outputs:                                                        !
!            sncovr )                                                   !
 
          call alcalc
!  ---  inputs:                                                         !
!          ( alb, snoalb, shdfac, shdmin, sncovr, tsnow,                !
!  ---  outputs:                                                        !
!            albedo )                                                   !
 
        endif   ! end if_sneqv_block
 
      endif   ! end if_ice_block
 
!  --- ...  thermal conductivity for sea-ice case, glacial-ice case
 
      if (ice /= 0) then
 
        df1 = 2.2
 
      else
 
!  --- ...  next calculate the subsurface heat flux, which first requires
!           calculation of the thermal diffusivity.  treatment of the
!           latter follows that on pages 148-149 from "heat transfer in
!           cold climates", by v. j. lunardini (published in 1981
!           by van nostrand reinhold co.) i.e. treatment of two contiguous
!           "plane parallel" mediums (namely here the first soil layer
!           and the snowpack layer, if any). this diffusivity treatment
!           behaves well for both zero and nonzero snowpack, including the
!           limit of very thin snowpack.  this treatment also eliminates
!           the need to impose an arbitrary upper bound on subsurface
!           heat flux when the snowpack becomes extremely thin.
 
!  --- ...   first calculate thermal diffusivity of top soil layer, using
!            both the frozen and liquid soil moisture, following the
!            soil thermal diffusivity function of peters-lidard et al.
!            (1998,jas, vol 55, 1209-1224), which requires the specifying
!            the quartz content of the given soil class (see routine redprm)
 
        call tdfcnd                                                     &
!  ---  inputs:
     &     ( smc(1), quartz, smcmax, sh2o(1),                           &
!  ---  outputs:
     &       df1                                                        &
     &     )
!       if(ivegsrc == 1) then
!only igbp type has urban
!urban
!           if ( vegtyp == 13 ) df1=3.24
!       endif
 
!wz only urban for igbp type
!
!jhan urban canopy heat storage effect is included in pbl scheme
!
        if((.not.lheatstrg) .and.                                       &
     &      (ivegsrc == 1 .and. vegtyp == 13)) then
          df1 = 3.24*(1.-shdfac) + shdfac*df1*exp(sbeta*shdfac)
        else
          df1 = df1 * exp( sbeta*shdfac )
        endif
 
      endif   ! end if_ice_block
 
!  --- ...  finally "plane parallel" snowpack effect following
!           v.j. linardini reference cited above. note that dtot is
!           combined depth of snowdepth and thickness of first soil layer
 
      dsoil = -0.5 * zsoil(1)
 
      if (sneqv == 0.0) then
 
        ssoil = df1 * (t1 - stc(1)) / dsoil
 
      else
 
        dtot = snowh + dsoil
        frcsno = snowh / dtot
        frcsoi = dsoil / dtot
 
!  --- ...  1. harmonic mean (series flow)
 
!       df1  = (sncond*df1) / (frcsoi*sncond + frcsno*df1)
!       df1h = (sncond*df1) / (frcsoi*sncond + frcsno*df1)
 
!  --- ...  2. arithmetic mean (parallel flow)
 
!       df1  = frcsno*sncond + frcsoi*df1
        df1a = frcsno*sncond + frcsoi*df1
 
!  --- ...  3. geometric mean (intermediate between harmonic and arithmetic mean)
 
!       df1 = (sncond**frcsno) * (df1**frcsoi)
!       df1 = df1h*sncovr + df1a*(1.0-sncovr)
!       df1 = df1h*sncovr + df1 *(1.0-sncovr)
        df1 = df1a*sncovr + df1 *(1.0-sncovr)
 
 
        ssoil = df1 * (t1 - stc(1)) / dtot
 
      endif   ! end if_sneqv_block
 
 
!     if (couple == 0) then            ! uncoupled mode
        if (sncovr > 0.0) then
 
          call snowz0
!  ---  inputs:                                                         !
!          ( sncovr,                                                    !
!  ---  input/outputs:                                                  !
!            z0 )                                                       !
 
        endif
!     endif
 
 
      t2v = sfctmp * (1.0 + 0.61*q2)
 
!  --- ...  next call routine sfcdif to calculate the sfc exchange coef (ch)
!           for heat and moisture.
!    note - comment out call sfcdif, if sfcdif already called in calling
!           program (such as in coupled atmospheric model).
!         - do not call sfcdif until after above call to redprm, in case
!           alternative values of roughness length (z0) and zilintinkevich
!           coef (czil) are set there via namelist i/o.
!         - routine sfcdif returns a ch that represents the wind spd times
!           the "original" nondimensional "ch" typical in literature.  hence
!           the ch returned from sfcdif has units of m/s.  the important
!           companion coefficient of ch, carried here as "rch", is the ch
!           from sfcdif times air density and parameter "cp".  "rch" is
!           computed in "call penman". rch rather than ch is the coeff
!           usually invoked later in eqns.
!         - sfcdif also returns the surface exchange coefficient for momentum,
!           cm, also known as the surface drage coefficient, but cm is not
!           used here.
 
!  --- ...  key required radiation term is the total downward radiation
!           (fdown) = net solar (swnet) + downward longwave (lwdn),
!           for use in penman ep calculation (penman) and other surface
!           energy budget calcuations.  also need downward solar (swdn)
!           for canopy resistance routine (canres).
!    note - fdown, swdn are derived differently in the uncoupled and
!           coupled modes.
 
 
      if (couple == 0) then                      !......uncoupled mode
 
!  --- ...  uncoupled mode:
!           compute surface exchange coefficients
 
        t1v  = t1  * (1.0 + 0.61 * q2)
        th2v = th2 * (1.0 + 0.61 * q2)
 
        call sfcdif
!  ---  inputs:                                                         !
!          ( zlvl, z0, t1v, th2v, sfcspd, czil,                         !
!  ---  input/outputs:                                                  !
!            cm, ch )                                                   !
 
!     swnet = net solar radiation into the ground (w/m2; dn-up) from input
!     fdown  = net solar + downward lw flux at sfc (w/m2)
 
        fdown = swnet + lwdn
 
      else                                       !......coupled mode
 
!  --- ...  coupled mode (couple .ne. 0):
!           surface exchange coefficients computed externally and passed in,
!           hence subroutine sfcdif not called.
 
!     swnet = net solar radiation into the ground (w/m2; dn-up) from input
!     fdown  = net solar + downward lw flux at sfc (w/m2)
 
        fdown = swnet + lwdn
 
      endif   ! end if_couple_block
 
 
      call penman
!  ---  inputs:                                                         !
!          ( sfctmp, sfcprs, sfcems, ch, t2v, th2, prcp, fdown,         !
!            ssoil, q2, q2sat, dqsdt2, snowng, frzgra,                  !
!  ---  outputs:                                                        !
!            t24, etp, rch, epsca, rr, flx2 )                           !
 
 
      if (shdfac > 0.) then
 
!  --- ...  frozen ground extension: total soil water "smc" was replaced
!           by unfrozen soil water "sh2o" in call to canres below
 
        call canres
!  ---  inputs:                                                         !
!          ( nsoil, nroot, swdn, ch, q2, q2sat, dqsdt2, sfctmp,         !
!            sfcprs, sfcems, sh2o, smcwlt, smcref, zsoil, rsmin,        !
!            rsmax, topt, rgl, hs, xlai,                                !
!  ---  outputs:                                                        !
!            rc, pc, rcs, rct, rcq, rcsoil )                            !
 
      endif
 
 
      esnow = 0.0
 
      if (sneqv .eq. 0.0) then
        call nopac
!  ---  inputs:                                                         !
!          ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref,           !
!            smcdry, cmcmax, dt, shdfac, sbeta, sfctmp, sfcems,         !
!            t24, th2, fdown, epsca, bexp, pc, rch, rr, cfactr,         !
!            slope, kdt, frzx, psisat, zsoil, dksat, dwsat,             !
!            zbot, ice, rtdis, quartz, fxexp, csoil, lheatstrg,         !
!  ---  input/outputs:                                                  !
!            cmc, t1, stc, sh2o, tbot,                                  !
!  ---  outputs:                                                        !
!            eta, smc, ssoil, runoff1, runoff2, runoff3, edir,          !
!            ec, et, ett, beta, drip, dew, flx1, flx3 )                 !
 
      else
 
        call snopac
!  ---  inputs:                                                         !
!          ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref, smcdry,   !
!            cmcmax, dt, df1, sfcems, sfctmp, t24, th2, fdown, epsca,   !
!            bexp, pc, rch, rr, cfactr, slope, kdt, frzx, psisat,       !
!            zsoil, dwsat, dksat, zbot, shdfac, ice, rtdis, quartz,     !
!            fxexp, csoil, flx2, snowng, lheatstrg,                     !
!  ---  input/outputs:                                                  !
!            prcp1, cmc, t1, stc, sncovr, sneqv, sndens, snowh,         !
!            sh2o, tbot, beta,                                          !
!  ---  outputs:                                                        !
!            smc, ssoil, runoff1, runoff2, runoff3, edir, ec, et,       !
!            ett, snomlt, drip, dew, flx1, flx3, esnow )                !
 
!  run-total accumulated snow based on snowfall and snowmelt in [m]
 
      endif
 
 
 
      sheat = -(ch*cp1*sfcprs) / (rd1*t2v) * (th2 - t1)
 
!   convert eta from kg m-2 s-1 to w m-2
!     eta = eta * lsubc
!     etp = etp * lsubc
 
      edir = edir * lsubc
      ec = ec * lsubc
 
      do k = 1, 4
        et(k) = et(k) * lsubc
      enddo
 
      ett = ett * lsubc
      esnow = esnow * lsubs
      etp = etp * ((1.0 - sncovr)*lsubc + sncovr*lsubs)
 
      if (etp > 0.) then
        eta = edir + ec + ett + esnow
      else
        eta = etp
      endif
 
      if (etp == 0.0) then
        beta = 0.0
      else
        beta = eta/etp
      endif
 
 
      ssoil = -1.0 * ssoil
 
      if (ice == 0) then
 
 
        runoff3 = runoff3 / dt
        runoff2 = runoff2 + runoff3
 
      else
 
 
        runoff1 = snomlt / dt
      endif
 
 
      soilm = -1.0 * smc(1) * zsoil(1)
      do k = 2, nsoil
        soilm = soilm + smc(k)*(zsoil(k-1) - zsoil(k))
      enddo
 
      soilwm = -1.0 * (smcmax - smcwlt) * zsoil(1)
      soilww = -1.0 * (smc(1) - smcwlt) * zsoil(1)
      do k = 2, nroot
        soilwm = soilwm + (smcmax - smcwlt) * (zsoil(k-1) - zsoil(k))
        soilww = soilww + (smc(k) - smcwlt) * (zsoil(k-1) - zsoil(k))
      enddo
 
      soilw = soilww / soilwm
!
      return
 
 
! =================
      contains
! =================
 
!*************************************!
!  section-1  1st level subprograms   !
!*************************************!
 
!-----------------------------------
      subroutine alcalc
!...................................
!  ---  inputs:
!    &     ( alb, snoalb, shdfac, shdmin, sncovr, tsnow,                &
!  ---  outputs:
!    &       albedo                                                     &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine alcalc calculates albedo including snow effect (0 -> 1)   !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from calling program:                                  size   !
!     alb      - real, snowfree albedo                             1    !
!     snoalb   - real, maximum (deep) snow albedo                  1    !
!     shdfac   - real, areal fractional coverage of green veg.     1    !
!     shdmin   - real, minimum areal coverage of green veg.        1    !
!     sncovr   - real, fractional snow cover                       1    !
!     tsnow    - real, snow surface temperature (k)                1    !
!                                                                       !
!  outputs to calling program:                                          !
!     albedo   - real, surface albedo including snow effect        1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
!     real (kind=kind_phys), intent(in) :: alb, snoalb, shdfac,         &
!    &       shdmin, sncovr, tsnow
 
!  ---  outputs:
!     real (kind=kind_phys), intent(out) :: albedo
 
!  ---  locals: (none)
 
!
!===> ...  begin here
!
!  --- ...  snoalb is argument representing maximum albedo over deep snow,
!           as passed into sflx, and adapted from the satellite-based
!           maximum snow albedo fields provided by d. robinson and g. kukla
!           (1985, jcam, vol 24, 402-411)
 
!         albedo = alb + (1.0-(shdfac-shdmin))*sncovr*(snoalb-alb)
          albedo = alb + sncovr*(snoalb - alb)
 
          if (albedo > snoalb) albedo = snoalb
 
!  --- ...  base formulation (dickinson et al., 1986, cogley et al., 1990)
 
!     if (tsnow <= 263.16) then
!       albedo = snoalb
!     else
!       if (tsnow < 273.16) then
!         tm = 0.1 * (tsnow - 263.16)
!         albedo = 0.5 * ((0.9 - 0.2*(tm**3)) + (0.8 - 0.16*(tm**3)))
!       else
!         albedo = 0.67
!       endif
!     endif
 
!  --- ...  isba formulation (verseghy, 1991; baker et al., 1990)
 
!     if (tsnow < 273.16) then
!       albedo = snoalb - 0.008*dt/86400
!     else
!       albedo = (snoalb - 0.5) * exp( -0.24*dt/86400 ) + 0.5
!     endif
 
!
!...................................
      end subroutine alcalc
!-----------------------------------
 
 
!-----------------------------------
      subroutine canres
!  ---  inputs:
!    &     ( nsoil, nroot, swdn, ch, q2, q2sat, dqsdt2, sfctmp,         &
!    &       sfcprs, sfcems, sh2o, smcwlt, smcref, zsoil, rsmin,        &
!    &       rsmax, topt, rgl, hs, xlai,                                &
!  ---  outputs:
!    &       rc, pc, rcs, rct, rcq, rcsoil                              &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine canres calculates canopy resistance which depends on      !
!  incoming solar radiation, air temperature, atmospheric water vapor   !
!  pressure deficit at the lowest model level, and soil moisture        !
!  (preferably unfrozen soil moisture rather than total)                !
!                                                                       !
!  source:   jarvis (1976), noilhan and planton (1989, mwr), jacquemin  !
!            and noilhan (1990, blm)                                    !
!  see also: chen et al (1996, jgr, vol 101(d3), 7251-7268), eqns       !
!            12-14 and table 2 of sec. 3.1.2                            !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from calling program:                                  size   !
!     nsoil    - integer, no. of soil layers                       1    !
!     nroot    - integer, no. of soil layers in root zone (<nsoil) 1    !
!     swdn     - real, incoming solar radiation                    1    !
!     ch       - real, sfc exchange coeff for heat and moisture    1    !
!     q2       - real, air humidity at 1st level above ground      1    !
!     q2sat    - real, sat. air humidity at 1st level abv ground   1    !
!     dqsdt2   - real, slope of sat. humidity function wrt temp    1    !
!     sfctmp   - real, sfc temperature at 1st level above ground   1    !
!     sfcprs   - real, sfc pressure                                1    !
!     sfcems   - real, sfc emissivity for lw radiation             1    !
!     sh2o     - real, volumetric soil moisture                  nsoil  !
!     smcwlt   - real, wilting point                               1    !
!     smcref   - real, reference soil moisture                     1    !
!     zsoil    - real, soil depth (negative sign, as below grd)  nsoil  !
!     rsmin    - real, mimimum stomatal resistance                 1    !
!     rsmax    - real, maximum stomatal resistance                 1    !
!     topt     - real, optimum transpiration air temperature       1    !
!     rgl      - real, canopy resistance func (in solar rad term)  1    !
!     hs       - real, canopy resistance func (vapor deficit term) 1    !
!     xlai     - real, leaf area index                             1    !
!                                                                       !
!  outputs to calling program:                                          !
!     rc       - real, canopy resistance                           1    !
!     pc       - real, plant coefficient                           1    !
!     rcs      - real, incoming solar rc factor                    1    !
!     rct      - real, air temp rc factor                          1    !
!     rcq      - real, atoms vapor press deficit rc factor         1    !
!     rcsoil   - real, soil moisture rc factor                     1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
!     integer, intent(in) :: nsoil, nroot
 
!     real (kind=kind_phys), intent(in) :: swdn, ch, q2, q2sat,         &
!    &       dqsdt2, sfctmp, sfcprs, sfcems, smcwlt, smcref, rsmin,     &
!    &       rsmax, topt, rgl, hs, xlai, sh2o(nsoil), zsoil(nsoil)
 
!  ---  outputs:
!     real (kind=kind_phys), intent(out) :: rc, pc, rcs, rct, rcq,      &
!    &       rcsoil
 
!  ---  locals:
      real (kind=kind_phys) :: delta, ff, gx, rr, part(nsold)
 
      integer :: k
 
!
!===> ...  begin here
!
!  --- ...  initialize canopy resistance multiplier terms.
 
      rcs = 0.0
      rct = 0.0
      rcq = 0.0
      rcsoil = 0.0
      rc = 0.0
 
!  --- ...  contribution due to incoming solar radiation
 
      ff = 0.55 * 2.0 * swdn / (rgl*xlai)
      rcs = (ff + rsmin/rsmax) / (1.0 + ff)
      rcs = max( rcs, 0.0001 )
 
!  --- ...  contribution due to air temperature at first model level above ground
!           rct expression from noilhan and planton (1989, mwr).
 
      rct = 1.0 - 0.0016 * (topt - sfctmp)**2.0
      rct = max( rct, 0.0001 )
 
!  --- ...  contribution due to vapor pressure deficit at first model level.
!           rcq expression from ssib
 
      rcq = 1.0 / (1.0 + hs*(q2sat-q2))
      rcq = max( rcq, 0.01 )
 
!  --- ...  contribution due to soil moisture availability.
!           determine contribution from each soil layer, then add them up.
 
      gx = (sh2o(1) - smcwlt) / (smcref - smcwlt)
      gx = max( 0.0, min( 1.0, gx ) )
 
!  --- ...  use soil depth as weighting factor
      part(1) = (zsoil(1)/zsoil(nroot)) * gx
 
!  --- ...  use root distribution as weighting factor
!     part(1) = rtdis(1) * gx
 
      do k = 2, nroot
 
        gx = (sh2o(k) - smcwlt) / (smcref - smcwlt)
        gx = max( 0.0, min( 1.0, gx ) )
 
!  --- ...  use soil depth as weighting factor
        part(k) = ((zsoil(k) - zsoil(k-1)) / zsoil(nroot)) * gx
 
!  --- ...  use root distribution as weighting factor
!       part(k) = rtdis(k) * gx
 
      enddo
 
      do k = 1, nroot
        rcsoil = rcsoil + part(k)
      enddo
      rcsoil = max( rcsoil, 0.0001 )
 
!  --- ...  determine canopy resistance due to all factors.  convert canopy
!           resistance (rc) to plant coefficient (pc) to be used with
!           potential evap in determining actual evap.  pc is determined by:
!           pc * linerized penman potential evap = penman-monteith actual
!           evaporation (containing rc term).
 
      rc = rsmin / (xlai*rcs*rct*rcq*rcsoil)
      rr = (4.0*sfcems*sigma1*rd1/cp1) * (sfctmp**4.0)/(sfcprs*ch) + 1.0
      delta = (lsubc/cp1) * dqsdt2
 
      pc = (rr + delta) / (rr*(1.0 + rc*ch) + delta)
!
!...................................
      end subroutine canres
!-----------------------------------
 
 
!-----------------------------------
      subroutine csnow
!...................................
!  ---  inputs:
!    &     ( sndens,                                                    &
!  ---  outputs:
!    &       sncond                                                     &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine csnow calculates snow termal conductivity                 !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from the calling program:                              size   !
!     sndens   - real, snow density                                1    !
!                                                                       !
!  outputs to the calling program:                                      !
!     sncond   - real, snow termal conductivity                    1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  constant parameters:
      real (kind=kind_phys), parameter :: unit = 0.11631
 
!  ---  inputs:
!     real (kind=kind_phys), intent(in) :: sndens
 
!  ---  outputs:
!     real (kind=kind_phys), intent(out) :: sncond
 
!  ---  locals:
      real (kind=kind_phys) :: c
 
!
!===> ...  begin here
!
!  --- ...  sncond in units of cal/(cm*hr*c), returned in w/(m*c)
!           basic version is dyachkova equation (1960), for range 0.1-0.4
 
      c = 0.328 * 10**(2.25*sndens)
      sncond = unit * c
 
!  --- ...  de vaux equation (1933), in range 0.1-0.6
 
!      sncond = 0.0293 * (1.0 + 100.0*sndens**2)
 
!  --- ...  e. andersen from flerchinger
 
!      sncond = 0.021 + 2.51 * sndens**2
!
!...................................
      end subroutine csnow
!-----------------------------------
 
 
!-----------------------------------
      subroutine nopac
!...................................
!  ---  inputs:
!    &     ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref,           &
!    &       smcdry, cmcmax, dt, shdfac, sbeta, sfctmp, sfcems,         &
!    &       t24, th2, fdown, epsca, bexp, pc, rch, rr, cfactr,         &
!    &       slope, kdt, frzx, psisat, zsoil, dksat, dwsat,             &
!    &       zbot, ice, rtdis, quartz, fxexp, csoil, lheatstrg,         &
!  ---  input/outputs:
!    &       cmc, t1, stc, sh2o, tbot,                                  &
!  ---  outputs:
!    &       eta, smc, ssoil, runoff1, runoff2, runoff3, edir,          &
!    &       ec, et, ett, beta, drip, dew, flx1, flx3                   &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine nopac calculates soil moisture and heat flux values and   !
!  update soil moisture content and soil heat content values for the    !
!  case when no snow pack is present.                                   !
!                                                                       !
!                                                                       !
!  subprograms called:  evapo, smflx, tdfcnd, shflx                     !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from calling program:                                  size   !
!     nsoil    - integer, number of soil layers                    1    !
!     nroot    - integer, number of root layers                    1    !
!     etp      - real, potential evaporation                       1    !
!     prcp     - real, precip rate                                 1    !
!     smcmax   - real, porosity (sat val of soil mois)             1    !
!     smcwlt   - real, wilting point                               1    !
!     smcref   - real, soil mois threshold                         1    !
!     smcdry   - real, dry soil mois threshold                     1    !
!     cmcmax   - real, maximum canopy water parameters             1    !
!     dt       - real, time step                                   1    !
!     shdfac   - real, aeral coverage of green veg                 1    !
!     sbeta    - real, param to cal veg effect on soil heat flux   1    !
!     sfctmp   - real, air temp at height zlvl abv ground          1    !
!     sfcems   - real, sfc lw emissivity                           1    !
!     t24      - real, sfctmp**4                                   1    !
!     th2      - real, air potential temp at zlvl abv grnd         1    !
!     fdown    - real, net solar + downward lw flux at sfc         1    !
!     epsca    - real,                                             1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     pc       - real, plant coeff                                 1    !
!     rch      - real, companion coefficient of ch                 1    !
!     rr       - real,                                             1    !
!     cfactr   - real, canopy water parameters                     1    !
!     slope    - real, linear reservoir coefficient                1    !
!     kdt      - real,                                             1    !
!     frzx     - real, frozen ground parameter                     1    !
!     psisat   - real, saturated soil potential                    1    !
!     zsoil    - real, soil layer depth below ground (negative)  nsoil  !
!     dksat    - real, saturated soil hydraulic conductivity       1    !
!     dwsat    - real, saturated soil diffusivity                  1    !
!     zbot     - real, specify depth of lower bd soil              1    !
!     ice      - integer, sea-ice flag (=1: sea-ice, =0: land)     1    !
!     rtdis    - real, root distribution                         nsoil  !
!     quartz   - real, soil quartz content                         1    !
!     fxexp    - real, bare soil evaporation exponent              1    !
!     csoil    - real, soil heat capacity                          1    !
!     lheatstrg- logical, flag for canopy heat storage             1    !
!                         parameterization                              !
!                                                                       !
!  input/outputs from and to the calling program:                       !
!     cmc      - real, canopy moisture content                     1    !
!     t1       - real, ground/canopy/snowpack eff skin temp        1    !
!     stc      - real, soil temp                                 nsoil  !
!     sh2o     - real, unfrozen soil moisture                    nsoil  !
!     tbot     - real, bottom soil temp                            1    !
!                                                                       !
!  outputs to the calling program:                                      !
!     eta      - real, downward latent heat flux                   1    !
!     smc      - real, total soil moisture                       nsoil  !
!     ssoil    - real, upward soil heat flux                       1    !
!     runoff1  - real, surface runoff not infiltrating sfc         1    !
!     runoff2  - real, sub surface runoff (baseflow)               1    !
!     runoff3  - real, excess of porosity                          1    !
!     edir     - real, direct soil evaporation                     1    !
!     ec       - real, canopy water evaporation                    1    !
!     et       - real, plant transpiration                       nsoil  !
!     ett      - real, total plant transpiration                   1    !
!     beta     - real, ratio of actual/potential evap              1    !
!     drip     - real, through-fall of precip and/or dew           1    !
!     dew      - real, dewfall (or frostfall)                      1    !
!     flx1     - real, precip-snow sfc flux                        1    !
!     flx3     - real, phase-change heat flux from snowmelt        1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
!     integer, intent(in) :: nsoil, nroot, ice
 
!     real (kind=kind_phys), intent(in) :: etp, prcp, smcmax,           &
!    &       smcwlt, smcref, smcdry, cmcmax, dt, shdfac, sbeta,         &
!    &       sfctmp, sfcems, t24, th2, fdown, epsca, bexp, pc,          &
!    &       rch, rr, cfactr, slope, kdt, frzx, psisat,                 &
!    &       zsoil(nsoil), dksat, dwsat, zbot, rtdis(nsoil),            &
!    &       quartz, fxexp, csoil
 
!     logical, intent(in) :: lheatstrg
 
!  ---  input/outputs:
!     real (kind=kind_phys), intent(inout) :: cmc, t1, stc(nsoil),      &
!    &       sh2o(nsoil), tbot
 
!  ---  outputs:
!     real (kind=kind_phys), intent(out) :: eta, smc(nsoil), ssoil,     &
!    &       runoff1, runoff2, runoff3, edir, ec, et(nsoil), ett,       &
!    &       beta, drip, dew, flx1, flx3
 
!  ---  locals:
      real (kind=kind_phys) :: df1, eta1, etp1, prcp1, yy, yynum,       &
     &       zz1, ec1, edir1, et1(nsoil), ett1
 
      integer :: k
 
!
!===> ...  begin here
!
!  --- ...  convert etp from kg m-2 s-1 to ms-1 and initialize dew.
 
      prcp1= prcp * 0.001
      etp1 = etp  * 0.001
      dew  = 0.0
      edir = 0.0
      edir1= 0.0
      ec   = 0.0
      ec1  = 0.0
 
      do k = 1, nsoil
        et(k) = 0.0
        et1(k) = 0.0
      enddo
 
      ett  = 0.0
      ett1 = 0.0
 
      if (etp > 0.0) then
 
!  --- ...  convert prcp from 'kg m-2 s-1' to 'm s-1'.
 
        call evapo                                                      &
!  ---  inputs:
     &     ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil,                &
     &       sh2o, smcmax, smcwlt, smcref, smcdry, pc,                  &
     &       shdfac, cfactr, rtdis, fxexp,                              &
!  ---  outputs:
     &       eta1, edir1, ec1, et1, ett1                                &
     &     )
 
        call smflx                                                      &
!  ---  inputs:
     &     ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1,             &
     &       zsoil, slope, frzx, bexp, dksat, dwsat, shdfac,            &
     &       edir1, ec1, et1,                                           &
!  ---  input/outputs:
     &       cmc, sh2o, smc,                                            &
!  ---  outputs:
     &       runoff1, runoff2, runoff3, drip                            &
     &     )
 
      else
 
!  --- ...  if etp < 0, assume dew forms (transform etp1 into dew and
!           reinitialize etp1 to zero).
 
        eta1 = 0.0
        dew  = -etp1
 
!  --- ...  convert prcp from 'kg m-2 s-1' to 'm s-1' and add dew amount.
 
        prcp1 = prcp1 + dew
 
        call smflx                                                      &
!  ---  inputs:
     &     ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1,             &
     &       zsoil, slope, frzx, bexp, dksat, dwsat, shdfac,            &
     &       edir1, ec1, et1,                                           &
!  ---  input/outputs:
     &       cmc, sh2o, smc,                                            &
!  ---  outputs:
     &       runoff1, runoff2, runoff3, drip                            &
     &     )
 
      endif   ! end if_etp_block
 
!  --- ...  convert modeled evapotranspiration fm  m s-1  to  kg m-2 s-1
 
      eta  = eta1 * 1000.0
      edir = edir1 * 1000.0
      ec   = ec1 * 1000.0
 
      do k = 1, nsoil
        et(k) = et1(k) * 1000.0
      enddo
 
      ett = ett1 * 1000.0
 
!  --- ...  based on etp and e values, determine beta
 
      if ( etp <= 0.0 ) then
        beta = 0.0
        if ( etp < 0.0 ) then
          beta = 1.0
        endif
      else
        beta = eta / etp
      endif
 
!  --- ...  get soil thermal diffuxivity/conductivity for top soil lyr,
!           calc. adjusted top lyr soil temp and adjusted soil flux, then
!           call shflx to compute/update soil heat flux and soil temps.
 
      call tdfcnd                                                       &
!  ---  inputs:
     &     ( smc(1), quartz, smcmax, sh2o(1),                           &
!  ---  outputs:
     &       df1                                                        &
     &     )
!      if(ivegsrc == 1) then
!urban
!        if ( vegtyp == 13 ) df1=3.24
!      endif
 
!  --- ... vegetation greenness fraction reduction in subsurface heat
!          flux via reduction factor, which is convenient to apply here
!          to thermal diffusivity that is later used in hrt to compute
!          sub sfc heat flux (see additional comments on veg effect
!          sub-sfc heat flx in routine sflx)
!wz only urban for igbp type
!
!jhan urban canopy heat storage effect is included in pbl scheme
!
        if((.not.lheatstrg) .and.                                       &
     &      (ivegsrc == 1 .and. vegtyp == 13)) then
          df1 = 3.24*(1.-shdfac) + shdfac*df1*exp(sbeta*shdfac)
        else
          df1 = df1 * exp( sbeta*shdfac )
        endif
 
!  --- ...  compute intermediate terms passed to routine hrt (via routine
!           shflx below) for use in computing subsurface heat flux in hrt
 
      yynum = fdown - sfcems*sigma1*t24
      yy = sfctmp + (yynum/rch + th2 - sfctmp - beta*epsca)/rr
      zz1 = df1/(-0.5*zsoil(1)*rch*rr) + 1.0
 
      call shflx                                                        &
!  ---  inputs:
     &     ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot,              &
     &       psisat, bexp, df1, ice, quartz, csoil, vegtyp,             &
     &       shdfac, lheatstrg,                                         &
!  ---  input/outputs:
     &       stc, t1, tbot, sh2o,                                       &
!  ---  outputs:
     &       ssoil                                                      &
     &     )
 
!  --- ...  set flx1 and flx3 (snopack phase change heat fluxes) to zero since
!           they are not used here in snopac.  flx2 (freezing rain heat flux)
!           was similarly initialized in the penman routine.
 
      flx1 = 0.0
      flx3 = 0.0
!
!...................................
      end subroutine nopac
!-----------------------------------
 
 
!-----------------------------------
      subroutine penman
!...................................
!  ---  inputs:
!    &     ( sfctmp, sfcprs, sfcems, ch, t2v, th2, prcp, fdown,         &
!    &       ssoil, q2, q2sat, dqsdt2, snowng, frzgra,                  &
!  ---  outputs:
!    &       t24, etp, rch, epsca, rr, flx2                             &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine penman calculates potential evaporation for the current   !
!  point.  various partial sums/products are also calculated and passed !
!  back to the calling routine for later use.                           !
!                                                                       !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     sfctmp   - real, sfc temperature at 1st level above ground   1    !
!     sfcprs   - real, sfc pressure                                1    !
!     sfcems   - real, sfc emissivity for lw radiation             1    !
!     ch       - real, sfc exchange coeff for heat & moisture      1    !
!     t2v      - real, sfc virtual temperature                     1    !
!     th2      - real, air potential temp at zlvl abv grnd         1    !
!     prcp     - real, precip rate                                 1    !
!     fdown    - real, net solar + downward lw flux at sfc         1    !
!     ssoil    - real, upward soil heat flux                       1    !
!     q2       - real, mixing ratio at hght zlvl abv ground        1    !
!     q2sat    - real, sat mixing ratio at zlvl abv ground         1    !
!     dqsdt2   - real, slope of sat specific humidity curve        1    !
!     snowng   - logical, snow flag                                1    !
!     frzgra   - logical, freezing rain flag                       1    !
!                                                                       !
!  outputs:                                                             !
!     t24      - real, sfctmp**4                                   1    !
!     etp      - real, potential evaporation                       1    !
!     rch      - real, companion coefficient of ch                 1    !
!     epsca    - real,                                             1    !
!     rr       - real,                                             1    !
!     flx2     - real, freezing rain latent heat flux              1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
!     real (kind=kind_phys), intent(in) :: sfctmp, sfcprs, sfcems,      &
!    &       ch, t2v, th2, prcp, fdown, ssoil, q2, q2sat, dqsdt2
 
!     logical, intent(in) :: snowng, frzgra
 
!  ---  outputs:
!     real (kind=kind_phys), intent(out) :: t24, etp, rch, epsca,       &
!    &       rr, flx2
 
!  ---  locals:
      real (kind=kind_phys) :: a, delta, fnet, rad, rho
 
!
!===> ...  begin here
!
      flx2 = 0.0
 
!  --- ...  prepare partial quantities for penman equation.
 
      delta = elcp * dqsdt2
      t24 = sfctmp * sfctmp * sfctmp * sfctmp
      rr  = t24 * 6.48e-8 / (sfcprs*ch) + 1.0
      rho = sfcprs / (rd1*t2v)
      rch = rho * cp * ch
 
!  --- ...  adjust the partial sums / products with the latent heat
!           effects caused by falling precipitation.
 
      if (.not. snowng) then
        if (prcp > 0.0)  rr = rr + cph2o1*prcp/rch
      else
! ---- ...  fractional snowfall/rainfall
        rr = rr + (cpice*ffrozp+cph2o1*(1.-ffrozp))                      &
     &       *prcp/rch
      endif
 
      fnet = fdown - sfcems*sigma1*t24 - ssoil
 
!  --- ...  include the latent heat effects of frzng rain converting to ice
!           on impact in the calculation of flx2 and fnet.
 
      if (frzgra) then
        flx2 = -lsubf * prcp
        fnet = fnet - flx2
      endif
 
!  --- ...  finish penman equation calculations.
 
      rad = fnet/rch + th2 - sfctmp
      a = elcp * (q2sat - q2)
      epsca = (a*rr + rad*delta) / (delta + rr)
      etp = epsca * rch / lsubc
!
!...................................
      end subroutine penman
!-----------------------------------
 
 
!-----------------------------------
      subroutine redprm(errmsg, errflg)
!...................................
!  ---  inputs:
!    &     ( nsoil, vegtyp, soiltyp, slopetyp, sldpth, zsoil,              &
!  ---  outputs:
!    &       cfactr, cmcmax, rsmin, rsmax, topt, refkdt, kdt,              &
!    &       sbeta, shdfac, rgl, hs, zbot, frzx, psisat, slope,            &
!    &       snup, salp, bexp, dksat, dwsat, smcmax, smcwlt,               &
!    &       smcref, smcdry, f1, quartz, fxexp, rtdis, nroot,              &
!    &       z0, czil, xlai, csoil                                         &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine redprm internally sets(default valuess), or optionally    !
!  read-in via namelist i/o, all soil and vegetation parameters         !
!  required for the execusion of the noah lsm.                          !
!                                                                       !
!  optional non-default parameters can be read in, accommodating up to  !
!  30 soil, veg, or slope classes, if the default max number of soil,   !
!  veg, and/or slope types is reset.                                    !
!                                                                       !
!  future upgrades of routine redprm must expand to incorporate some    !
!  of the empirical parameters of the frozen soil and snowpack physics  !
!  (such as in routines frh2o, snowpack, and snow_new) not yet set in   !
!  this redprm routine, but rather set in lower level subroutines.      !
!                                                                       !
!  all soil, veg, slope, and universal parameters values are defined    !
!  externally (in subroutine "set_soilveg.f") and then accessed via     !
!  "use namelist_soilveg" (below) and then set here.                    !
!                                                                       !
!  soil types   zobler (1986)      cosby et al (1984) (quartz cont.(1)) !
!      1         coarse            loamy sand            (0.82)         !
!      2         medium            silty clay loam       (0.10)         !
!      3         fine              light clay            (0.25)         !
!      4         coarse-medium     sandy loam            (0.60)         !
!      5         coarse-fine       sandy clay            (0.52)         !
!      6         medium-fine       clay loam             (0.35)         !
!      7         coarse-med-fine   sandy clay loam       (0.60)         !
!      8         organic           loam                  (0.40)         !
!      9         glacial land ice  loamy sand            (na using 0.82)!
!     13: <old>- glacial land ice -<old>                                !
!     13:        glacial-ice (no longer use these parameters), now      !
!                treated as ice-only surface and sub-surface            !
!                (in subroutine hrtice)                                 !
!  upgraded to statsgo (19-type)
!     1: sand
!     2: loamy sand
!     3: sandy loam
!     4: silt loam
!     5: silt
!     6:loam
!     7:sandy clay loam
!     8:silty clay loam
!     9:clay loam
!     10:sandy clay
!     11: silty clay
!     12: clay
!     13: organic material
!     14: water
!     15: bedrock
!     16: other (land-ice)
!     17: playa
!     18: lava
!     19: white sand
!                                                                       !
!  ssib vegetation types (dorman and sellers, 1989; jam)                !
!      1:  broadleaf-evergreen trees  (tropical forest)                 !
!      2:  broadleaf-deciduous trees                                    !
!      3:  broadleaf and needleleaf trees (mixed forest)                !
!      4:  needleleaf-evergreen trees                                   !
!      5:  needleleaf-deciduous trees (larch)                           !
!      6:  broadleaf trees with groundcover (savanna)                   !
!      7:  groundcover only (perennial)                                 !
!      8:  broadleaf shrubs with perennial groundcover                  !
!      9:  broadleaf shrubs with bare soil                              !
!     10:  dwarf trees and shrubs with groundcover (tundra)             !
!     11:  bare soil                                                    !
!     12:  cultivations (the same parameters as for type 7)             !
!     13: <old>- glacial (the same parameters as for type 11) -<old>    !
!     13:  glacial-ice (no longer use these parameters), now treated as !
!          ice-only surface and sub-surface (in subroutine hrtice)      !
!  upgraded to IGBP (20-type)
!      1:Evergreen Needleleaf Forest
!      2:Evergreen Broadleaf Forest
!      3:Deciduous Needleleaf Forest
!      4:Deciduous Broadleaf Forest
!      5:Mixed Forests
!      6:Closed Shrublands
!      7:Open Shrublands
!      8:Woody Savannas
!      9:Savannas
!      10:Grasslands
!      11:Permanent wetlands
!      12:Croplands
!      13:Urban and Built-Up
!      14:Cropland/natural vegetation mosaic
!      15:Snow and Ice
!      16:Barren or Sparsely Vegetated
!      17:Water
!      18:Wooded Tundra
!      19:Mixed Tundra
!      20:Bare Ground Tundra
!                                                                       !
!  slopetyp is to estimate linear reservoir coefficient slope to the    !
!  baseflow runoff out of the bottom layer. lowest class (slopetyp=0)   !
!  means highest slope parameter = 1.                                   !
!                                                                       !
!  slope class       percent slope                                      !
!      1                0-8                                             !
!      2                8-30                                            !
!      3                > 30                                            !
!      4                0-30                                            !
!      5                0-8 & > 30                                      !
!      6                8-30 & > 30                                     !
!      7                0-8, 8-30, > 30                                 !
!      9                glacial ice                                     !
!    blank              ocean/sea                                       !
!                                                                       !
!  note: class 9 from zobler file should be replaced by 8 and 'blank' 9 !
!                                                                       !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from calling program:                                  size   !
!     nsoil    - integer, number of soil layers                    1    !
!     vegtyp   - integer, vegetation type (integer index)          1    !
!     soiltyp  - integer, soil type (integer index)                1    !
!     slopetyp - integer, class of sfc slope (integer index)       1    !
!     sldpth   - integer, thickness of each soil layer (m)       nsoil  !
!     zsoil    - integer, soil depth (negative sign) (m)         nsoil  !
!                                                                       !
!  outputs to the calling program:                                      !
!     cfactr   - real, canopy water parameters                     1    !
!     cmcmax   - real, maximum canopy water parameters             1    !
!     rsmin    - real, mimimum stomatal resistance                 1    !
!     rsmax    - real, maximum stomatal resistance                 1    !
!     topt     - real, optimum transpiration air temperature       1    !
!     refkdt   - real, =2.e-6 the sat. dk. val for soil type 2     1    !
!     kdt      - real,                                             1    !
!     sbeta    - real, param to cal veg effect on soil heat flux   1    !
!     shdfac   - real, vegetation greenness fraction               1    !
!     rgl      - real, canopy resistance func (in solar rad term)  1    !
!     hs       - real, canopy resistance func (vapor deficit term) 1    !
!     zbot     - real, specify depth of lower bd soil temp (m)     1    !
!     frzx     - real, frozen ground parameter, ice content        1    !
!                      threshold above which frozen soil is impermeable !
!     psisat   - real, saturated soil potential                    1    !
!     slope    - real, linear reservoir coefficient                1    !
!     snup     - real, threshold snow depth (water equi m)         1    !
!     salp     - real, snow cover shape parameter                  1    !
!                      from anderson's hydro-17 best fit salp = 2.6     !
!     bexp     - real, the 'b' parameter                           1    !
!     dksat    - real, saturated soil hydraulic conductivity       1    !
!     dwsat    - real, saturated soil diffusivity                  1    !
!     smcmax   - real, max soil moisture content (porosity)        1    !
!     smcwlt   - real, wilting pt soil moisture contents           1    !
!     smcref   - real, reference soil moisture (onset stress)      1    !
!     smcdry   - real, air dry soil moist content limits           1    !
!     f1       - real, used to comp soil diffusivity/conductivity  1    !
!     quartz   - real, soil quartz content                         1    !
!     fxexp    - real, bare soil evaporation exponent              1    !
!     rtdis    - real, root distribution                         nsoil  !
!     nroot    - integer, number of root layers                    1    !
!     z0       - real, roughness length (m)                        1    !
!     czil     - real, param to cal roughness length of heat       1    !
!     xlai     - real, leaf area index                             1    !
!     csoil    - real, soil heat capacity (j m-3 k-1)              1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
      use namelist_soilveg
 
!  ---  input:
!     integer, intent(in) :: nsoil, vegtyp, soiltyp, slopetyp
 
!     real (kind=kind_phys), intent(in) :: sldpth(nsoil), zsoil(nsoil)
 
!  ---  outputs:
!     real (kind=kind_phys), intent(out) :: cfactr, cmcmax, rsmin,      &
!    &       rsmax, topt, refkdt, kdt, sbeta, shdfac, rgl, hs, zbot,    &
!    &       frzx, psisat, slope, snup, salp, bexp, dksat, dwsat,       &
!    &       smcmax, smcwlt, smcref, smcdry, f1, quartz, fxexp, z0,     &
!    &       czil, xlai, csoil, rtdis(nsoil)
      character(len=*),     intent(out) :: errmsg
      integer,              intent(out) :: errflg
!     integer, intent(out) :: nroot
 
!  ---  locals:
      real (kind=kind_phys) :: frzfact, frzk, refdk
 
      integer :: i
 
!
!===> ...  begin here
!
! Initialize CCPP error-handling
      errflg = 0
      errmsg = ''
 
      if (soiltyp > defined_soil) then
        write(*,*) 'warning: too many soil types,soiltyp=',soiltyp,     &
     &   'defined_soil=',defined_soil
        errflg = 1
        errmsg = 'ERROR(sflx.f): too many soil types'
        return
      endif
 
      if (vegtyp > defined_veg) then
        write(*,*) 'warning: too many veg types'
        errflg = 1
        errmsg = 'ERROR(sflx.f): too many veg types'
        return
      endif
 
      if (slopetyp > defined_slope) then
        write(*,*) 'warning: too many slope types'
        errflg = 1
        errmsg = 'ERROR(sflx.f): too many slope types'
        return
      endif
 
!  --- ...  set-up universal parameters (not dependent on soiltyp, vegtyp
!           or slopetyp)
 
      zbot   = zbot_data
      salp   = salp_data
      cfactr = cfactr_data
      cmcmax = cmcmax_data
      sbeta  = sbeta_data
      rsmax  = rsmax_data
      topt   = topt_data
      refdk  = refdk_data
      frzk   = frzk_data
      fxexp  = fxexp_data
      refkdt = refkdt_data
      czil   = czil_data
      csoil  = csoil_data
 
!  --- ...  set-up soil parameters
 
      bexp  = bb(soiltyp)
      dksat = satdk(soiltyp)
      dwsat = satdw(soiltyp)
      f1    = f11(soiltyp)
      kdt   = refkdt * dksat / refdk
 
      psisat = satpsi(soiltyp)
      quartz = qtz(soiltyp)
      smcdry = drysmc(soiltyp)
      smcmax = maxsmc(soiltyp)
      smcref = refsmc(soiltyp)
      smcwlt = wltsmc(soiltyp)
 
      frzfact = (smcmax / smcref) * (0.412 / 0.468)
 
!  --- ...  to adjust frzk parameter to actual soil type: frzk * frzfact
 
      frzx = frzk * frzfact
 
!  --- ...  set-up vegetation parameters
 
      nroot = nroot_data(vegtyp)
      snup  = snupx(vegtyp)
      rsmin = rsmtbl(vegtyp)
 
      rgl = rgltbl(vegtyp)
      hs  = hstbl(vegtyp)
! roughness lengthe is defined in sfcsub
!     z0  = z0_data(vegtyp)
      xlai= lai_data(vegtyp)
 
      if (vegtyp == bare) shdfac = 0.0
 
      if (nroot > nsoil) then
        errflg = 1
        errmsg = 'ERROR(sflx.f): too many root layers'
        return
      endif
 
!  --- ...  calculate root distribution.  present version assumes uniform
!           distribution based on soil layer depths.
 
      do i = 1, nroot
        rtdis(i) = -sldpth(i) / zsoil(nroot)
      enddo
 
!  --- ...  set-up slope parameter
 
      slope = slope_data(slopetyp)
!
!...................................
      end subroutine redprm
!-----------------------------------
 
 
!-----------------------------------
      subroutine sfcdif
!...................................
!  ---  inputs:
!    &     ( zlvl, z0, t1v, th2v, sfcspd, czil,                         &
!  ---  input/outputs:
!    &       cm, ch                                                     &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine sfcdif calculates surface layer exchange coefficients     !
!  via iterative process. see chen et al (1997, blm)                    !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from the calling program:                              size   !
!     zlvl     - real, height abv atmos ground forcing vars (m)    1    !
!     z0       - real, roughness length (m)                        1    !
!     t1v      - real, surface exchange coefficient                1    !
!     th2v     - real, surface exchange coefficient                1    !
!     sfcspd   - real, wind speed at height zlvl abv ground (m/s)  1    !
!     czil     - real, param to cal roughness length of heat       1    !
!                                                                       !
!  input/outputs from and to the calling program:                       !
!     cm       - real, sfc exchange coeff for momentum (m/s)       1    !
!     ch       - real, sfc exchange coeff for heat & moisture (m/s)1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  --- constant parameters:
      integer,               parameter :: itrmx  = 5
      real (kind=kind_phys), parameter :: wwst   = 1.2
      real (kind=kind_phys), parameter :: wwst2  = wwst*wwst
      real (kind=kind_phys), parameter :: vkrm   = 0.40
      real (kind=kind_phys), parameter :: excm   = 0.001
      real (kind=kind_phys), parameter :: beta   = 1.0/270.0
      real (kind=kind_phys), parameter :: btg    = beta*gs1
      real (kind=kind_phys), parameter :: elfc   = vkrm*btg
      real (kind=kind_phys), parameter :: wold   = 0.15
      real (kind=kind_phys), parameter :: wnew   = 1.0-wold
      real (kind=kind_phys), parameter :: pihf   = 3.14159265/2.0  ! con_pi/2.0
 
      real (kind=kind_phys), parameter :: epsu2  = 1.e-4
      real (kind=kind_phys), parameter :: epsust = 0.07
      real (kind=kind_phys), parameter :: ztmin  = -5.0
      real (kind=kind_phys), parameter :: ztmax  = 1.0
      real (kind=kind_phys), parameter :: hpbl   = 1000.0
      real (kind=kind_phys), parameter :: sqvisc = 258.2
 
      real (kind=kind_phys), parameter :: ric    = 0.183
      real (kind=kind_phys), parameter :: rric   = 1.0/ric
      real (kind=kind_phys), parameter :: fhneu  = 0.8
      real (kind=kind_phys), parameter :: rfc    = 0.191
      real (kind=kind_phys), parameter :: rfac   = ric/(fhneu*rfc*rfc)
 
!  ---  inputs:
!     real (kind=kind_phys),  intent(in) :: zlvl, z0, t1v, th2v,        &
!    &       sfcspd, czil
 
!  ---  input/outputs:
!     real (kind=kind_phys),  intent(inout) :: cm, ch
 
!  ---  locals:
      real (kind=kind_phys) :: zilfc, zu, zt, rdz, cxch, dthv, du2,     &
     &       btgh, wstar2, ustar, zslu, zslt, rlogu, rlogt, rlmo,       &
     &       zetalt, zetalu, zetau, zetat, xlu4, xlt4, xu4, xt4,        &
     &       xlu, xlt, xu, xt, psmz, simm, pshz, simh, ustark,          &
     &       rlmn, rlma
 
      integer :: ilech, itr
 
!  ---  define local in-line functions:
 
      real (kind=kind_phys) :: pslmu, pslms, pslhu, pslhs, zz
      real (kind=kind_phys) :: pspmu, pspms, psphu, psphs, xx, yy
 
!  ...  1) lech's surface functions
 
      pslmu( zz ) = -0.96 * log( 1.0-4.5*zz )
      pslms( zz ) = zz*rric - 2.076*(1.0 - 1.0/(zz + 1.0))
      pslhu( zz ) = -0.96 * log( 1.0-4.5*zz )
      pslhs( zz ) = zz*rfac - 2.076*(1.0 - 1.0/(zz + 1.0))
 
!  ...  2) paulson's surface functions
 
      pspmu( xx ) = -2.0 * log( (xx + 1.0)*0.5 )                        &
     &            - log( (xx*xx + 1.0)*0.5 ) + 2.0*atan(xx) - pihf
      pspms( yy ) = 5.0 * yy
      psphu( xx ) = -2.0 * log( (xx*xx + 1.0)*0.5 )
      psphs( yy ) = 5.0 * yy
 
!
!===> ...  begin here
!
!  --- ...  this routine sfcdif can handle both over open water (sea, ocean) and
!           over solid surface (land, sea-ice).
 
      ilech = 0
 
!   --- ...  ztfc: ratio of zoh/zom  less or equal than 1
!            czil: constant c in zilitinkevich, s. s.1995,:note about zt
 
      zilfc = -czil * vkrm * sqvisc
 
      zu = z0
 
      rdz = 1.0 / zlvl
      cxch = excm * rdz
      dthv = th2v - t1v
      du2 = max( sfcspd*sfcspd, epsu2 )
 
!  --- ...  beljars correction of ustar
 
      btgh = btg * hpbl
 
!  --- ...  if statements to avoid tangent linear problems near zero
      if (btgh*ch*dthv /= 0.0) then
        wstar2 = wwst2 * abs( btgh*ch*dthv )**(2.0/3.0)
      else
        wstar2 = 0.0
      endif
 
      ustar = max( sqrt( cm*sqrt( du2+wstar2 ) ), epsust )
 
!  --- ...  zilitinkevitch approach for zt
 
      zt = exp( zilfc*sqrt( ustar*z0 ) ) * z0
 
      zslu = zlvl + zu
      zslt = zlvl + zt
 
!     print*,'zslt=',zslt
!     print*,'zlvl=',zvll
!     print*,'zt=',zt
 
      rlogu = log( zslu/zu )
      rlogt = log( zslt/zt )
 
      rlmo = elfc*ch*dthv / ustar**3
 
!     print*,'rlmo=',rlmo
!     print*,'elfc=',elfc
!     print*,'ch=',ch
!     print*,'dthv=',dthv
!     print*,'ustar=',ustar
 
      do itr = 1, itrmx
 
!  --- ...  1./ monin-obukkhov length-scale
 
        zetalt = max( zslt*rlmo, ztmin )
        rlmo   = zetalt / zslt
        zetalu = zslu * rlmo
        zetau  = zu * rlmo
        zetat  = zt * rlmo
 
        if (ilech == 0) then
 
          if (rlmo < 0.0) then
            xlu4 = 1.0 - 16.0 * zetalu
            xlt4 = 1.0 - 16.0 * zetalt
            xu4  = 1.0 - 16.0 * zetau
            xt4  = 1.0 - 16.0* zetat
 
            xlu = sqrt( sqrt( xlu4 ) )
            xlt = sqrt( sqrt( xlt4 ) )
            xu  = sqrt( sqrt( xu4  ) )
            xt  = sqrt( sqrt( xt4  ) )
 
            psmz = pspmu(xu)
 
!           print*,'-----------1------------'
!           print*,'psmz=',psmz
!           print*,'pspmu(zetau)=',pspmu( zetau )
!           print*,'xu=',xu
!           print*,'------------------------'
 
            simm = pspmu( xlu ) - psmz + rlogu
            pshz = psphu( xt  )
            simh = psphu( xlt ) - pshz + rlogt
          else
            zetalu = min( zetalu, ztmax )
            zetalt = min( zetalt, ztmax )
            psmz = pspms( zetau )
 
!           print*,'-----------2------------'
!           print*,'psmz=',psmz
!           print*,'pspms(zetau)=',pspms( zetau )
!           print*,'zetau=',zetau
!           print*,'------------------------'
 
            simm = pspms( zetalu ) - psmz + rlogu
            pshz = psphs( zetat  )
            simh = psphs( zetalt ) - pshz + rlogt
          endif   ! end if_rlmo_block
 
        else
 
!  --- ...  lech's functions
 
          if (rlmo < 0.0) then
            psmz = pslmu( zetau )
 
!           print*,'-----------3------------'
!           print*,'psmz=',psmz
!           print*,'pslmu(zetau)=',pslmu( zetau )
!           print*,'zetau=',zetau
!           print*,'------------------------'
 
            simm = pslmu( zetalu ) - psmz + rlogu
            pshz = pslhu( zetat  )
            simh = pslhu( zetalt ) - pshz + rlogt
          else
            zetalu = min( zetalu, ztmax )
            zetalt = min( zetalt, ztmax )
 
            psmz = pslms( zetau  )
 
!           print*,'-----------4------------'
!           print*,'psmz=',psmz
!           print*,'pslms(zetau)=',pslms( zetau )
!           print*,'zetau=',zetau
!           print*,'------------------------'
 
            simm = pslms( zetalu ) - psmz + rlogu
            pshz = pslhs( zetat  )
            simh = pslhs( zetalt ) - pshz + rlogt
          endif   ! end if_rlmo_block
 
        endif   ! end if_ilech_block
 
!  --- ...  beljaars correction for ustar
 
        ustar = max( sqrt( cm*sqrt( du2+wstar2 ) ), epsust )
 
!  --- ...  zilitinkevitch fix for zt
 
        zt = exp( zilfc*sqrt( ustar*z0 ) ) * z0
 
        zslt = zlvl + zt
        rlogt = log( zslt/zt )
 
        ustark = ustar * vkrm
        cm = max( ustark/simm, cxch )
        ch = max( ustark/simh, cxch )
 
!  --- ...  if statements to avoid tangent linear problems near zero
 
        if (btgh*ch*dthv /= 0.0) then
          wstar2 = wwst2 * abs(btgh*ch*dthv) ** (2.0/3.0)
        else
          wstar2 = 0.0
        endif
 
        rlmn = elfc*ch*dthv / ustar**3
        rlma = rlmo*wold + rlmn*wnew
 
        rlmo = rlma
 
      enddo   ! end do_itr_loop
 
!     print*,'----------------------------'
!     print*,'sfcdif output !  ! ! ! ! ! ! ! !  !   !    !'
!
!     print*,'zlvl=',zlvl
!     print*,'z0=',z0
!     print*,'t1v=',t1v
!     print*,'th2v=',th2v
!     print*,'sfcspd=',sfcspd
!     print*,'czil=',czil
!     print*,'cm=',cm
!     print*,'ch=',ch
!     print*,'----------------------------'
!
!...................................
      end subroutine sfcdif
!-----------------------------------
 
 
!-----------------------------------
      subroutine snfrac
!...................................
!  ---  inputs:
!    &     ( sneqv, snup, salp, snowh,                                  &
!  ---  outputs:
!    &       sncovr                                                     &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine snfrac calculatexsnow fraction (0 -> 1)                   !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from the calling program:                              size   !
!     sneqv    - real, snow water equivalent (m)                   1    !
!     snup     - real, threshold sneqv depth above which sncovr=1  1    !
!     salp     - real, tuning parameter                            1    !
!     snowh    - real, snow depth (m)                              1    !
!                                                                       !
!  outputs to the calling program:                                      !
!     sncovr   - real, fractional snow cover                       1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
!     real (kind=kind_phys),  intent(in) :: sneqv, snup, salp, snowh
 
!  ---  outputs:
!     real (kind=kind_phys),  intent(out) :: sncovr
 
!  ---  locals:
      real (kind=kind_phys) :: rsnow, z0n
 
!
!===> ...  begin here
!
!  --- ...  snup is veg-class dependent snowdepth threshhold (set in routine
!           redprm) above which snocvr=1.
 
          if (sneqv < snup) then
            rsnow = sneqv / snup
            sncovr = 1.0 - (exp(-salp*rsnow) - rsnow*exp(-salp))
          else
            sncovr = 1.0
          endif
 
          z0n = 0.035
 
!  --- ...  formulation of dickinson et al. 1986
 
!       sncovr = snowh / (snowh + 5.0*z0n)
 
!  --- ...  formulation of marshall et al. 1994
 
!       sncovr = sneqv / (sneqv + 2.0*z0n)
 
!
!...................................
      end subroutine snfrac
!-----------------------------------
 
 
!-----------------------------------
      subroutine snopac
!...................................
!  ---  inputs:
!    &     ( nsoil, nroot, etp, prcp, smcmax, smcwlt, smcref, smcdry,   &
!    &       cmcmax, dt, df1, sfcems, sfctmp, t24, th2, fdown, epsca,   &
!    &       bexp, pc, rch, rr, cfactr, slope, kdt, frzx, psisat,       &
!    &       zsoil, dwsat, dksat, zbot, shdfac, ice, rtdis, quartz,     &
!    &       fxexp, csoil, flx2, snowng, lheatstrg,                     &
!  ---  input/outputs:
!    &       prcp1, cmc, t1, stc, sncovr, sneqv, sndens, snowh,         &
!    &       sh2o, tbot, beta,                                          &
!  ---  outputs:
!    &       smc, ssoil, runoff1, runoff2, runoff3, edir, ec, et,       &
!    &       ett, snomlt, drip, dew, flx1, flx3, esnow                  &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine snopac calculates soil moisture and heat flux values and  !
!  update soil moisture content and soil heat content values for the    !
!  case when a snow pack is present.                                    !
!                                                                       !
!                                                                       !
!  subprograms called:  evapo, smflx, shflx, snowpack
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from the calling program:                              size   !
!     nsoil    - integer, number of soil layers                    1    !
!     nroot    - integer, number of root layers                    1    !
!     etp      - real, potential evaporation                       1    !
!     prcp     - real, precip rate                                 1    !
!     smcmax   - real, porosity                                    1    !
!     smcwlt   - real, wilting point                               1    !
!     smcref   - real, soil mois threshold                         1    !
!     smcdry   - real, dry soil mois threshold                     1    !
!     cmcmax   - real, maximum canopy water parameters             1    !
!     dt       - real, time step                                   1    !
!     df1      - real, thermal diffusivity                         m    !
!     sfcems   - real, lw surface emissivity                       1    !
!     sfctmp   - real, sfc temperature                             1    !
!     t24      - real, sfctmp**4                                   1    !
!     th2      - real, sfc air potential temperature               1    !
!     fdown    - real, net solar + downward lw flux at sfc         1    !
!     epsca    - real,                                             1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     pc       - real, plant coeff                                 1    !
!     rch      - real, companion coefficient of ch                 1    !
!     rr       - real,                                             1    !
!     cfactr   - real, canopy water parameters                     1    !
!     slope    - real, linear reservoir coefficient                1    !
!     kdt      - real,                                             1    !
!     frzx     - real, frozen ground parameter                     1    !
!     psisat   - real, saturated soil potential                    1    !
!     zsoil    - real, soil layer depth below ground (negative)  nsoil  !
!     dwsat    - real, saturated soil diffusivity                  1    !
!     dksat    - real, saturated soil hydraulic conductivity       1    !
!     zbot     - real, specify depth of lower bd soil              1    !
!     shdfac   - real, aeral coverage of green vegetation          1    !
!     ice      - integer, sea-ice flag (=1: sea-ice, =0: land)     1    !
!     rtdis    - real, root distribution                         nsoil  !
!     quartz   - real, soil quartz content                         1    !
!     fxexp    - real, bare soil evaporation exponent              1    !
!     csoil    - real, soil heat capacity                          1    !
!     flx2     - real, freezing rain latent heat flux              1    !
!     snowng   - logical, snow flag                                1    !
!     lheatstrg- logical, flag for canopy heat storage             1    !
!                         parameterization                              !
!                                                                       !
!  input/outputs from and to the calling program:                       !
!     prcp1    - real, effective precip                            1    !
!     cmc      - real, canopy moisture content                     1    !
!     t1       - real, ground/canopy/snowpack eff skin temp        1    !
!     stc      - real, soil temperature                          nsoil  !
!     sncovr   - real, snow cover                                  1    !
!     sneqv    - real, water-equivalent snow depth                 1    !
!     sndens   - real, snow density                                1    !
!     snowh    - real, snow depth                                  1    !
!     sh2o     - real, unfrozen soil moisture                    nsoil  !
!     tbot     - real, bottom soil temperature                     1    !
!     beta     - real, ratio of actual/potential evap              1    !
!                                                                       !
!  outputs to the calling program:                                      !
!     smc      - real, total soil moisture                       nsoil  !
!     ssoil    - real, upward soil heat flux                       1    !
!     runoff1  - real, surface runoff not infiltrating sfc         1    !
!     runoff2  - real, sub surface runoff                          1    !
!     runoff3  - real, excess of porosity for a given soil layer   1    !
!     edir     - real, direct soil evaporation                     1    !
!     ec       - real, canopy water evaporation                    1    !
!     et       - real, plant transpiration                       nsoil  !
!     ett      - real, total plant transpiration                   1    !
!     snomlt   - real, snow melt water equivalent                  1    !
!     drip     - real, through-fall of precip                      1    !
!     dew      - real, dewfall (or frostfall)                      1    !
!     flx1     - real, precip-snow sfc flux                        1    !
!     flx3     - real, phase-change heat flux from snowmelt        1    !
!     esnow    - real, sublimation from snowpack                   1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  constant parameters:
      real, parameter :: esdmin = 1.e-6
 
!  ---  inputs:
!     integer, intent(in) :: nsoil, nroot, ice
 
!     real (kind=kind_phys), intent(in) :: etp, prcp, smcmax, smcref,   &
!    &       smcwlt, smcdry, cmcmax, dt, df1, sfcems, sfctmp, t24,      &
!    &       th2, fdown, epsca, bexp, pc, rch, rr, cfactr, slope, kdt,  &
!    &       frzx, psisat, dwsat, dksat, zbot, shdfac, quartz,          &
!    &       csoil, fxexp, flx2, zsoil(nsoil), rtdis(nsoil)
 
!     logical, intent(in) :: snowng
!
!     logical, intent(in) :: lheatstrg
!
 
!  ---  input/outputs:
!     real (kind=kind_phys), intent(inout) :: prcp1, t1, sncovr, sneqv, &
!    &       sndens, snowh, cmc, tbot, beta, sh2o(nsoil), stc(nsoil)
 
!  ---  outputs:
!     real (kind=kind_phys), intent(out) :: ssoil, runoff1, runoff2,    &
!    &       runoff3, edir, ec, et(nsoil), ett, snomlt, drip, dew,      &
!    &       flx1, flx3, esnow, smc(nsoil)
 
!  ---  locals:
      real (kind=kind_phys):: denom, dsoil, dtot, etp1, ssoil1,         &
     &       snoexp, ex, t11, t12, t12a, t12b, yy, zz1, seh, t14,       &
     &       ec1, edir1, ett1, etns, etns1, esnow1, esnow2, etanrg,     &
     &       et1(nsoil)
 
      integer k
 
!     data snoexp /1.0/    !!! <----- for noah v2.7
      data snoexp /2.0/    !!! <----- for noah v2.7.1
 
!  --- ...  convert potential evap (etp) from kg m-2 s-1 to m s-1 and then to an
!           amount (m) given timestep (dt) and call it an effective snowpack
!           reduction amount, esnow2 (m) for a snowcover fraction = 1.0.  this is
!           the amount the snowpack would be reduced due to sublimation from the
!           snow sfc during the timestep.  sublimation will proceed at the
!           potential rate unless the snow depth is less than the expected
!           snowpack reduction.  for snowcover fraction = 1.0, 0=edir=et=ec, and
!           hence total evap = esnow = sublimation (potential evap rate)
 
!  --- ...  if sea-ice (ice=1) or glacial-ice (ice=-1), snowcover fraction = 1.0,
!           and sublimation is at the potential rate.
!           for non-glacial land (ice=0), if snowcover fraction < 1.0, total
!           evaporation < potential due to non-potential contribution from
!           non-snow covered fraction.
 
      prcp1 = prcp1 * 0.001
 
      edir  = 0.0
      edir1 = 0.0
 
      ec  = 0.0
      ec1 = 0.0
 
      do k = 1, nsoil
        et(k) = 0.0
        et1(k) = 0.0
      enddo
 
      ett   = 0.0
      ett1  = 0.0
      etns  = 0.0
      etns1 = 0.0
      esnow = 0.0
      esnow1= 0.0
      esnow2= 0.0
 
      dew = 0.0
      etp1 = etp * 0.001
 
      if (etp < 0.0) then
 
!  --- ...  if etp<0 (downward) then dewfall (=frostfall in this case).
 
        dew = -etp1
        esnow2 = etp1 * dt
        etanrg = etp * ((1.0-sncovr)*lsubc + sncovr*lsubs)
 
      else
 
!  --- ...  etp >= 0, upward moisture flux
 
        if (ice /= 0) then           ! for sea-ice and glacial-ice case
 
          esnow = etp
          esnow1 = esnow * 0.001
          esnow2 = esnow1 * dt
          etanrg = esnow * lsubs
 
        else                         ! for non-glacial land case
 
          if (sncovr < 1.0) then
 
            call evapo                                                  &
!  ---  inputs:
     &     ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil,                &
     &       sh2o, smcmax, smcwlt, smcref, smcdry, pc,                  &
     &       shdfac, cfactr, rtdis, fxexp,                              &
!  ---  outputs:
     &       etns1, edir1, ec1, et1, ett1                               &
     &     )
 
            edir1 = edir1 * (1.0 - sncovr)
            ec1 = ec1 * (1.0 - sncovr)
 
            do k = 1, nsoil
              et1(k) = et1(k) * (1.0 - sncovr)
            enddo
 
            ett1  = ett1  * (1.0 - sncovr)
            etns1 = etns1 * (1.0 - sncovr)
 
            edir = edir1 * 1000.0
            ec = ec1 * 1000.0
 
            do k = 1, nsoil
              et(k) = et1(k) * 1000.0
            enddo
 
            ett = ett1 * 1000.0
            etns = etns1 * 1000.0
 
          endif   ! end if_sncovr_block
 
          esnow = etp * sncovr
!         esnow1 = etp * 0.001
          esnow1 = esnow * 0.001
          esnow2 = esnow1 * dt
          etanrg = esnow*lsubs + etns*lsubc
 
        endif   ! end if_ice_block
 
      endif   ! end if_etp_block
 
!  --- ...  if precip is falling, calculate heat flux from snow sfc to newly
!           accumulating precip.  note that this reflects the flux appropriate for
!           the not-yet-updated skin temperature (t1).  assumes temperature of the
!           snowfall striking the gound is =sfctmp (lowest model level air temp).
 
      flx1 = 0.0
      if ( snowng ) then
!  --- ... fractional snowfall/rainfall
        flx1 = (cpice* ffrozp + cph2o1*(1.-ffrozp))                     &
     &         * prcp * (t1 - sfctmp)
      else
        if (prcp > 0.0) flx1 = cph2o1 * prcp * (t1 - sfctmp)
      endif
 
!  --- ...  calculate an 'effective snow-grnd sfc temp' (t12) based on heat
!           fluxes between the snow pack and the soil and on net radiation.
!           include flx1 (precip-snow sfc) and flx2 (freezing rain latent
!           heat) fluxes.
!           flx2 reflects freezing rain latent heat flux using t1 calculated
!           in penman.
 
      dsoil = -0.5 * zsoil(1)
      dtot = snowh + dsoil
      denom = 1.0 + df1 / (dtot * rr * rch)
 
!     t12a = ( (fdown - flx1 - flx2 - sigma1*t24) / rch                 &
!    &     + th2 - sfctmp - beta*epsca ) / rr
      t12a = ( (fdown - flx1 - flx2 - sfcems*sigma1*t24) / rch          &
     &     + th2 - sfctmp - etanrg/rch ) / rr
 
      t12b = df1 * stc(1) / (dtot * rr * rch)
      t12 = (sfctmp + t12a + t12b) / denom
 
!  --- ...  if the 'effective snow-grnd sfc temp' is at or below freezing, no snow
!           melt will occur.  set the skin temp to this effective temp.  reduce
!           (by sublimination ) or increase (by frost) the depth of the snowpack,
!           depending on sign of etp.
!           update soil heat flux (ssoil) using new skin temperature (t1)
!           since no snowmelt, set accumulated snowmelt to zero, set 'effective'
!           precip from snowmelt to zero, set phase-change heat flux from snowmelt
!           to zero.
 
      if (t12 <= tfreez) then
 
        t1 = t12
        ssoil = df1 * (t1 - stc(1)) / dtot
!wz     ssoil = (t1 - stc (1)) * max(7.0, df1/dtot)
        sneqv = max(0.0, sneqv-esnow2)
        flx3 = 0.0
        ex = 0.0
        snomlt = 0.0
 
      else
 
!  --- ...  if the 'effective snow-grnd sfc temp' is above freezing, snow melt
!           will occur.  call the snow melt rate,ex and amt, snomlt.  revise the
!           effective snow depth.  revise the skin temp because it would have chgd
!           due to the latent heat released by the melting. calc the latent heat
!           released, flx3. set the effective precip, prcp1 to the snow melt rate,
!           ex for use in smflx.  adjustment to t1 to account for snow patches.
!           calculate qsat valid at freezing point.  note that esat (saturation
!           vapor pressure) value of 6.11e+2 used here is that valid at frzzing
!           point.  note that etp from call penman in sflx is ignored here in
!           favor of bulk etp over 'open water' at freezing temp.
!           update soil heat flux (s) using new skin temperature (t1)
 
!  --- ...  noah v2.7.1   mek feb2004
!           non-linear weighting of snow vs non-snow covered portions of gridbox
!           so with snoexp = 2.0 (>1), surface skin temperature is higher than
!           for the linear case (snoexp = 1).
 
!       t1 = tfreez * sncovr**snoexp + t12 * (1.0 - sncovr**snoexp)
        t1 = tfreez * max(0.01,sncovr**snoexp) +                          &
     &          t12 * (1.0 - max(0.01,sncovr**snoexp))
 
        beta = 1.0
        ssoil = df1 * (t1 - stc(1)) / dtot
 
!  --- ...  if potential evap (sublimation) greater than depth of snowpack.
!           beta<1
!           snowpack has sublimated away, set depth to zero.
 
        if (sneqv-esnow2 <= esdmin) then
 
          sneqv = 0.0
          ex = 0.0
          snomlt = 0.0
          flx3 = 0.0
 
        else
 
!  --- ...  potential evap (sublimation) less than depth of snowpack, retain
!           beta=1.
 
          sneqv = sneqv - esnow2
          seh = rch * (t1 - th2)
 
          t14 = t1 * t1
          t14 = t14 * t14
 
          flx3 = fdown - flx1 - flx2 - sfcems*sigma1*t14                &
     &         - ssoil - seh - etanrg
          if (flx3 <= 0.0) flx3 = 0.0
 
          ex = flx3 * 0.001 / lsubf
 
!  --- ...  snowmelt reduction depending on snow cover
!           if snow cover less than 5% no snowmelt reduction
!     note: does 'if' below fail to match the melt water with the melt
!           energy?
 
!         if (sncovr > 0.05) ex = ex * sncovr
          snomlt = ex * dt
 
!  --- ...  esdmin represents a snowpack depth threshold value below which we
!           choose not to retain any snowpack, and instead include it in snowmelt.
 
          if (sneqv-snomlt >= esdmin) then
 
            sneqv = sneqv - snomlt
 
          else
 
!  --- ...  snowmelt exceeds snow depth
 
            ex = sneqv / dt
            flx3 = ex * 1000.0 * lsubf
            snomlt = sneqv
            sneqv = 0.0
 
          endif   ! end if_sneqv-snomlt_block
 
        endif   ! end if_sneqv-esnow2_block
 
!       prcp1 = prcp1 + ex
 
!  --- ...  if non-glacial land, add snowmelt rate (ex) to precip rate to be used
!           in subroutine smflx (soil moisture evolution) via infiltration.
 
!  --- ...  for sea-ice and glacial-ice, the snowmelt will be added to subsurface
!           runoff/baseflow later near the end of sflx (after return from call to
!           subroutine snopac)
 
        if (ice == 0) prcp1 = prcp1 + ex
 
      endif   ! end if_t12<=tfreez_block
 
!  --- ...  final beta now in hand, so compute evaporation.  evap equals etp
!           unless beta<1.
 
!      eta = beta * etp
 
!  --- ...  smflx returns updated soil moisture values for non-glacial land.
!           if sea-ice (ice=1) or glacial-ice (ice=-1), skip call to smflx, since
!           no soil medium for sea-ice or glacial-ice
 
      if (ice == 0) then
 
        call smflx                                                        &
!  ---  inputs:
     &     ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1,               &
     &       zsoil, slope, frzx, bexp, dksat, dwsat, shdfac,              &
     &       edir1, ec1, et1,                                             &
!  ---  input/outputs:
     &       cmc, sh2o, smc,                                              &
!  ---  outputs:
     &       runoff1, runoff2, runoff3, drip                              &
     &     )
 
      endif
 
!  --- ...  before call shflx in this snowpack case, set zz1 and yy arguments to
!           special values that ensure that ground heat flux calculated in shflx
!           matches that already computed for below the snowpack, thus the sfc
!           heat flux to be computed in shflx will effectively be the flux at the
!           snow top surface.  t11 is a dummy arguement so we will not use the
!           skin temp value as revised by shflx.
 
      zz1 = 1.0
      yy  = stc(1) - 0.5*ssoil*zsoil(1)*zz1 / df1
      t11 = t1
 
!  --- ...  shflx will calc/update the soil temps.  note:  the sub-sfc heat flux
!           (ssoil1) and the skin temp (t11) output from this shflx call are not
!           used  in any subsequent calculations. rather, they are dummy variables
!           here in the snopac case, since the skin temp and sub-sfc heat flux are
!           updated instead near the beginning of the call to snopac.
 
      call shflx                                                        &
!  ---  inputs:
     &     ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot,              &
     &       psisat, bexp, df1, ice, quartz, csoil, vegtyp,             &
     &       shdfac, lheatstrg,                                         &
!  ---  input/outputs:
     &       stc, t11, tbot, sh2o,                                      &
!  ---  outputs:
     &       ssoil1                                                     &
     &     )
 
!  --- ...  snow depth and density adjustment based on snow compaction.  yy is
!           assumed to be the soil temperture at the top of the soil column.
 
      if (ice == 0) then              ! for non-glacial land
 
        if (sneqv > 0.0) then
 
          call snowpack                                                 &
!  ---  inputs:
     &     ( sneqv, dt, t1, yy,                                         &
!  ---  input/outputs:
     &       snowh, sndens                                              &
     &     )
 
        else
 
          sneqv = 0.0
          snowh = 0.0
          sndens = 0.0
!         sncond = 1.0
          sncovr = 0.0
 
        endif   ! end if_sneqv_block
 
!  --- ...  over sea-ice or glacial-ice, if s.w.e. (sneqv) below threshold lower
!           bound (0.01 m for sea-ice, 0.10 m for glacial-ice), then set at
!           lower bound and store the source increment in subsurface runoff/
!           baseflow (runoff2).  note:  runoff2 is then a negative value (as
!           a flag) over sea-ice or glacial-ice, in order to achieve water balance.
 
      elseif (ice == 1) then          ! for sea-ice
 
        if (sneqv >= 0.01) then
 
          call snowpack                                                 &
!  ---  inputs:
     &     ( sneqv, dt, t1, yy,                                         &
!  ---  input/outputs:
     &       snowh, sndens                                              &
     &     )
 
        else
 
!         sndens = sneqv / snowh
!         runoff2 = -(0.01 - sneqv) / dt
          sneqv = 0.01
          snowh = 0.05
          sncovr = 1.0
!         snowh = sneqv / sndens
 
        endif   ! end if_sneqv_block
 
      else                            ! for glacial-ice
 
        if (sneqv >= 0.10) then
 
          call snowpack                                                 &
!  ---  inputs:
     &     ( sneqv, dt, t1, yy,                                         &
!  ---  input/outputs:
     &       snowh, sndens                                              &
     &     )
 
        else
 
!         sndens = sneqv / snowh
!         runoff2 = -(0.10 - sneqv) / dt
          sneqv = 0.10
          snowh = 0.50
          sncovr = 1.0
!         snowh = sneqv / sndens
 
        endif   ! end if_sneqv_block
 
      endif   ! end if_ice_block
 
!
!...................................
      end subroutine snopac
!-----------------------------------
 
 
!-----------------------------------
      subroutine snow_new
!...................................
!  ---  inputs:
!    &     ( sfctmp, sn_new, rhonewsn, exticeden,                       &
!  ---  input/outputs:
!    &       snowh, sndens                                              &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!    subroutine snow_new calculates snow depth and densitity to account !
!    for the new snowfall. new values of snow depth & density returned. !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from the calling program:                              size   !
!     sfctmp   - real, surface air temperature (k)                 1    !
!     sn_new   - real, new snowfall (m)                            1    !
!                                                                       !
!  input/outputs from and to the calling program:                       !
!     snowh    - real, snow depth (m)                              1    !
!     sndens   - real, snow density                                1    !
!                      (g/cm3=dimensionless fraction of h2o density)    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
!     real(kind=kind_phys), intent(in) :: sfctmp, sn_new
 
!  ---  input/outputs:
!     real(kind=kind_phys), intent(inout) :: snowh, sndens
 
!  ---  locals:
      real(kind=kind_phys) :: dsnew, snowhc, hnewc, newsnc, tempc
 
!
!===> ...  begin here
!
!  --- ...  conversion into simulation units
 
      snowhc = snowh * 100.0
      newsnc = sn_new * 100.0
      tempc  = sfctmp - tfreez
 
!  --- ...  calculating new snowfall density depending on temperature
!           equation from gottlib l. 'a general runoff model for
!           snowcovered and glacierized basin', 6th nordic hydrological
!           conference, vemadolen, sweden, 1980, 172-177pp.
 
      if(.not. exticeden) then
         if (tempc <= -15.0) then
            dsnew = 0.05
         else
            dsnew = 0.05 + 0.0017*(tempc + 15.0)**1.5
         endif
      else
         dsnew = rhonewsn*0.001
      endif
 
!  --- ...  adjustment of snow density depending on new snowfall
 
      hnewc  = newsnc / dsnew
      sndens = (snowhc*sndens + hnewc*dsnew) / (snowhc + hnewc)
      snowhc = snowhc + hnewc
      snowh  = snowhc * 0.01
!
!...................................
      end subroutine snow_new
!-----------------------------------
 
 
!-----------------------------------
      subroutine snowz0
!...................................
!  ---  inputs:
!    &     ( sncovr,                                                    &
!  ---  input/outputs:
!    &       z0                                                         &
!    &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!    subroutine snowz0 calculates total roughness length over snow      !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from the calling program:                              size   !
!     sncovr   - real, fractional snow cover                       1    !
!                                                                       !
!  input/outputs from and to the calling program:                       !
!     z0       - real, roughness length (m)                        1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
!     real(kind=kind_phys), intent(in) :: sncovr
 
!  ---  input/outputs:
!     real(kind=kind_phys), intent(inout) :: z0
 
!  ---  locals:
      real(kind=kind_phys) :: z0s
!
!===> ...  begin here
!
!     z0s = 0.001                     ! snow roughness length:=0.001 (m)
!  --- ...  current noah lsm condition - mbek, 09-oct-2001
      z0s = z0
 
      z0 = (1.0 - sncovr)*z0 + sncovr*z0s
 
!
!...................................
      end subroutine snowz0
!-----------------------------------
 
 
!-----------------------------------
      subroutine tdfcnd                                                 &
!  ---  inputs:
     &     ( smc, qz, smcmax, sh2o,                                     &
!  ---  outputs:
     &       df                                                         &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!    subroutine tdfcnd calculates thermal diffusivity and conductivity  !
!    of the soil for a given point and time.                            !
!                                                                       !
!    peters-lidard approach (peters-lidard et al., 1998)                !
!    june 2001 changes: frozen soil condition.                          !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!  use as in peters-lidard, 1998 (modif. from johansen, 1975).          !
!                                 pablo grunmann, 08/17/98              !
!  refs.:                                                               !
!    farouki, o.t.,1986: thermal properties of soils. series on rock    !
!             and soil mechanics, vol. 11, trans tech, 136 pp.          !
!    johansen, o., 1975: thermal conductivity of soils. ph.d. thesis,   !
!             university of trondheim,                                  !
!    peters-lidard, c. d., et al., 1998: the effect of soil thermal     !
!             conductivity parameterization on surface energy fluxes    !
!             and temperatures. journal of the atmospheric sciences,    !
!             vol. 55, pp. 1209-1224.                                   !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     smc      - real, top layer total soil moisture               1    !
!     qz       - real, quartz content (soil type dependent)        1    !
!     smcmax   - real, porosity                                    1    !
!     sh2o     - real, top layer unfrozen soil moisture            1    !
!                                                                       !
!  outputs:                                                             !
!     df       - real, soil thermal diffusivity and conductivity   1    !
!                                                                       !
!  locals:                                                              !
!     thkw     - water thermal conductivity                        1    !
!     thkqtz   - thermal conductivity for quartz                   1    !
!     thko     - thermal conductivity for other soil components    1    !
!     thkqtz   - thermal conductivity for the solids combined      1    !
!     thkice   - ice thermal conductivity                          1    !
!                                                                       !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  input:
      real (kind=kind_phys), intent(in) :: smc, qz, smcmax, sh2o
 
!  ---  output:
      real (kind=kind_phys), intent(out) :: df
 
!  ---  locals:
      real (kind=kind_phys) :: gammd, thkdry, ake, thkice, thko,        &
     &       thkqtz, thksat, thks, thkw, satratio, xu, xunfroz
!
!===> ...  begin here
!
!  --- ...  if the soil has any moisture content compute a partial sum/product
!           otherwise use a constant value which works well with most soils
 
!  --- ...  saturation ratio:
 
      satratio = smc / smcmax
 
!  --- ...  parameters  w/(m.k)
      thkice = 2.2
      thkw   = 0.57
      thko   = 2.0
!     if (qz <= 0.2) thko = 3.0
      thkqtz = 7.7
 
!  --- ...  solids' conductivity
 
      thks = (thkqtz**qz) * (thko**(1.0-qz))
 
!  --- ...  unfrozen fraction (from 1., i.e., 100%liquid, to 0. (100% frozen))
 
      xunfroz = (sh2o + 1.e-9) / (smc + 1.e-9)
 
!  --- ...  unfrozen volume for saturation (porosity*xunfroz)
 
      xu=xunfroz*smcmax
 
!  --- ...  saturated thermal conductivity
 
      thksat = thks**(1.-smcmax) * thkice**(smcmax-xu) * thkw**(xu)
 
!  --- ...  dry density in kg/m3
 
      gammd = (1.0 - smcmax) * 2700.0
 
!  --- ...  dry thermal conductivity in w.m-1.k-1
 
      thkdry = (0.135*gammd + 64.7) / (2700.0 - 0.947*gammd)
 
      if ( sh2o+0.0005 < smc ) then         ! frozen
 
        ake = satratio
 
      else                                  ! unfrozen
 
!  --- ...  range of validity for the kersten number (ake)
        if ( satratio > 0.1 ) then
 
!  --- ...  kersten number (using "fine" formula, valid for soils containing
!           at least 5% of particles with diameter less than 2.e-6 meters.)
!           (for "coarse" formula, see peters-lidard et al., 1998).
 
          ake = log10( satratio ) + 1.0
 
        else
 
!  --- ...  use k = kdry
          ake = 0.0
 
        endif   ! end if_satratio_block
 
      endif   ! end if_sh2o+0.0005_block
 
!  --- ...  thermal conductivity
 
      df = ake * (thksat - thkdry) + thkdry
!
!...................................
      end subroutine tdfcnd
!-----------------------------------
 
 
!*********************************************!
!  section-2  2nd level subprograms           !
!*********************************************!
 
 
!-----------------------------------
      subroutine evapo                                                  &
!  ---  inputs:
     &     ( nsoil, nroot, cmc, cmcmax, etp1, dt, zsoil,                &
     &       sh2o, smcmax, smcwlt, smcref, smcdry, pc,                  &
     &       shdfac, cfactr, rtdis, fxexp,                              &
!  ---  outputs:
     &       eta1, edir1, ec1, et1, ett1                                &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine evapo calculates soil moisture flux.  the soil moisture   !
!  content (smc - a per unit volume measurement) is a dependent variable!
!  that is updated with prognostic eqns. the canopy moisture content    !
!  (cmc) is also updated. frozen ground version:  new states added:     !
!  sh2o, and frozen ground correction factor, frzfact and parameter     !
!  slope.                                                               !
!                                                                       !
!                                                                       !
!  subprogram called:  devap, transp                                    !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs from calling program:                                  size   !
!     nsoil    - integer, number of soil layers                    1    !
!     nroot    - integer, number of root layers                    1    !
!     cmc      - real, canopy moisture content                     1    !
!     cmcmax   - real, maximum canopy water parameters             1    !
!     etp1     - real, potential evaporation                       1    !
!     dt       - real, time step                                   1    !
!     zsoil    - real, soil layer depth below ground             nsoil  !
!     sh2o     - real, unfrozen soil moisture                    nsoil  !
!     smcmax   - real, porosity                                    1    !
!     smcwlt   - real, wilting point                               1    !
!     smcref   - real, soil mois threshold                         1    !
!     smcdry   - real, dry soil mois threshold                     1    !
!     pc       - real, plant coeff                                 1    !
!     cfactr   - real, canopy water parameters                     1    !
!     rtdis    - real, root distribution                         nsoil  !
!     fxexp    - real, bare soil evaporation exponent              1    !
!                                                                       !
!  outputs to calling program:                                          !
!     eta1     - real, latent heat flux                            1    !
!     edir1    - real, direct soil evaporation                     1    !
!     ec1      - real, canopy water evaporation                    1    !
!     et1      - real, plant transpiration                       nsoil  !
!     ett1     - real, total plant transpiration                   1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
      integer, intent(in) :: nsoil, nroot
 
      real (kind=kind_phys),  intent(in) :: cmc, cmcmax, etp1, dt, pc,  &
     &       smcmax, smcwlt, smcref, smcdry, shdfac, cfactr, fxexp,     &
     &       zsoil(nsoil), sh2o(nsoil), rtdis(nsoil)
 
!  ---  outputs:
      real (kind=kind_phys),  intent(out) :: eta1, edir1, ec1, ett1,    &
     &       et1(nsoil)
 
!  ---  locals:
      real (kind=kind_phys) :: cmc2ms
 
      integer :: i, k
 
!
!===> ...  begin here
!
!  --- ...  executable code begins here if the potential evapotranspiration
!           is greater than zero.
 
      edir1 = 0.0
      ec1 = 0.0
 
      do k = 1, nsoil
        et1(k) = 0.0
      enddo
      ett1 = 0.0
 
      if (etp1 > 0.0) then
 
!  --- ...  retrieve direct evaporation from soil surface.  call this function
!           only if veg cover not complete.
!           frozen ground version:  sh2o states replace smc states.
 
        if (shdfac < 1.0) then
 
          call devap                                                    &
!  ---  inputs:
     &     ( etp1, sh2o(1), shdfac, smcmax, smcdry, fxexp,              &
!  ---  outputs:
     &       edir1                                                      &
     &     )
 
        endif
 
!  --- ...  initialize plant total transpiration, retrieve plant transpiration,
!           and accumulate it for all soil layers.
 
        if (shdfac > 0.0) then
 
          call transp                                                   &
!  ---  inputs:
     &     ( nsoil, nroot, etp1, sh2o, smcwlt, smcref,                   &
     &       cmc, cmcmax, zsoil, shdfac, pc, cfactr, rtdis,             &
!  ---  outputs:
     &       et1                                                        &
     &     )
 
          do k = 1, nsoil
            ett1 = ett1 + et1(k)
          enddo
 
!  --- ...  calculate canopy evaporation.
!           if statements to avoid tangent linear problems near cmc=0.0.
 
          if (cmc > 0.0) then
            ec1 = shdfac * ( (cmc/cmcmax)**cfactr ) * etp1
          else
            ec1 = 0.0
          endif
 
!  --- ...  ec should be limited by the total amount of available water
!           on the canopy.  -f.chen, 18-oct-1994
 
          cmc2ms = cmc / dt
          ec1 = min( cmc2ms, ec1 )
        endif
 
      endif   ! end if_etp1_block
 
!  --- ...  total up evap and transp types to obtain actual evapotransp
 
      eta1 = edir1 + ett1 + ec1
 
!
!...................................
      end subroutine evapo
!-----------------------------------
 
 
!-----------------------------------
      subroutine shflx                                                  &
!  ---  inputs:
     &     ( nsoil, smc, smcmax, dt, yy, zz1, zsoil, zbot,              &
     &       psisat, bexp, df1, ice, quartz, csoil, vegtyp,             &
     &       shdfac, lheatstrg,                                         &
!  ---  input/outputs:
     &       stc, t1, tbot, sh2o,                                       &
!  ---  outputs:
     &       ssoil                                                      &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine shflx updates the temperature state of the soil column    !
!  based on the thermal diffusion equation and update the frozen soil   !
!  moisture content based on the temperature.                           !
!                                                                       !
!  subprogram called:  hstep, hrtice, hrt                               !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     smc      - real, total soil moisture                       nsoil  !
!     smcmax   - real, porosity (sat val of soil mois)             1    !
!     dt       - real, time step                                   1    !
!     yy       - real, soil temperature at the top of column       1    !
!     zz1      - real,                                             1    !
!     zsoil    - real, soil layer depth below ground (negative)  nsoil  !
!     zbot     - real, specify depth of lower bd soil              1    !
!     psisat   - real, saturated soil potential                    1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     df1      - real, thermal diffusivity and conductivity        1    !
!     ice      - integer, sea-ice flag (=1: sea-ice, =0: land)     1    !
!     quartz   - real, soil quartz content                         1    !
!     csoil    - real, soil heat capacity                          1    !
!     vegtyp   - integer, vegtation type                           1    !
!     shdfac   - real, aeral coverage of green vegetation          1    !
!     lheatstrg- logical, flag for canopy heat storage             1    ! 
!                         parameterization                              !
!                                                                       !
!  input/outputs:                                                       !
!     stc      - real, soil temp                                 nsoil  !
!     t1       - real, ground/canopy/snowpack eff skin temp        1    !
!     tbot     - real, bottom soil temp                            1    !
!     sh2o     - real, unfrozen soil moisture                    nsoil  !
!                                                                       !
!  outputs:                                                             !
!     ssoil    - real, upward soil heat flux                       1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  parameter constants:
      real (kind=kind_phys), parameter :: ctfil1 = 0.5
      real (kind=kind_phys), parameter :: ctfil2 = 1.0 - ctfil1
 
!  ---  inputs:
      integer, intent(in) :: nsoil, ice, vegtyp
 
      real (kind=kind_phys), intent(in) :: smc(nsoil), smcmax, dt, yy,  &
     &       zz1, zsoil(nsoil), zbot, psisat, bexp, df1, quartz, csoil, &
     &       shdfac
 
      logical, intent(in) :: lheatstrg
 
!  ---  input/outputs:
      real (kind=kind_phys), intent(inout) :: stc(nsoil), t1, tbot,     &
     &       sh2o(nsoil)
 
!  ---  outputs:
      real (kind=kind_phys), intent(out) :: ssoil
 
!  ---  locals:
      real (kind=kind_phys) :: ai(nsold), bi(nsold), ci(nsold), oldt1,  &
     &       rhsts(nsold), stcf(nsold), stsoil(nsoil)
 
      integer :: i
 
!
!===> ...  begin here
!
      oldt1 = t1
      do i = 1, nsoil
         stsoil(i) = stc(i)
      enddo
 
!  --- ...  hrt routine calcs the right hand side of the soil temp dif eqn
 
      if (ice /= 0) then
 
!  --- ...  sea-ice case, glacial-ice case
 
        call hrtice                                                     &
!  ---  inputs:
     &     ( nsoil, stc, zsoil, yy, zz1, df1, ice,                      &
!  ---  input/outputs:
     &       tbot,                                                      &
!  ---  outputs:
     &       rhsts, ai, bi, ci                                          &
     &     )
 
        call hstep                                                      &
!  ---  inputs:
     &     ( nsoil, stc, dt,                                            &
!  ---  input/outputs:
     &       rhsts, ai, bi, ci,                                         &
!  ---  outputs:
     &       stcf                                                       &
     &     )
 
      else
 
!  --- ...  land-mass case
 
        call hrt                                                        &
!  ---  inputs:
     &     ( nsoil, stc, smc, smcmax, zsoil, yy, zz1, tbot,             &
     &       zbot, psisat, dt, bexp, df1, quartz, csoil,vegtyp,         &
     &       shdfac, lheatstrg,                                         &
!  ---  input/outputs:
     &       sh2o,                                                      &
!  ---  outputs:
     &       rhsts, ai, bi, ci                                          &
     &     )
 
        call hstep                                                      &
!  ---  inputs:
     &     ( nsoil, stc, dt,                                            &
!  ---  input/outputs:
     &       rhsts, ai, bi, ci,                                         &
!  ---  outputs:
     &       stcf                                                       &
     &     )
 
      endif
 
      do i = 1, nsoil
         stc(i) = stcf(i)
      enddo
 
!  --- ...  in the no snowpack case (via routine nopac branch,) update the grnd
!           (skin) temperature here in response to the updated soil temperature
!           profile above.  (note: inspection of routine snopac shows that t1
!           below is a dummy variable only, as skin temperature is updated
!           differently in routine snopac)
 
      t1 = (yy + (zz1 - 1.0)*stc(1)) / zz1
      t1 = ctfil1*t1 + ctfil2*oldt1
 
      do i = 1, nsoil
        stc(i) = ctfil1*stc(i) + ctfil2*stsoil(i)
      enddo
 
!  --- ...  calculate surface soil heat flux
 
      ssoil = df1*(stc(1) - t1) / (0.5*zsoil(1))
 
!
!...................................
      end subroutine shflx
!-----------------------------------
 
 
 
!-----------------------------------
      subroutine smflx                                                    &
!  ---  inputs:
     &     ( nsoil, dt, kdt, smcmax, smcwlt, cmcmax, prcp1,               &
     &       zsoil, slope, frzx, bexp, dksat, dwsat, shdfac,              &
     &       edir1, ec1, et1,                                             &
!  ---  input/outputs:
     &       cmc, sh2o, smc,                                              &
!  ---  outputs:
     &       runoff1, runoff2, runoff3, drip                              &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine smflx calculates soil moisture flux.  the soil moisture   !
!  content (smc - a per unit volume measurement) is a dependent variable!
!  that is updated with prognostic eqns. the canopy moisture content    !
!  (cmc) is also updated. frozen ground version:  new states added: sh2o!
!  and frozen ground correction factor, frzx and parameter slope.       !
!                                                                       !
!                                                                       !
!  subprogram called:  srt, sstep                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     dt       - real, time step                                   1    !
!     kdt      - real,                                             1    !
!     smcmax   - real, porosity                                    1    !
!     smcwlt   - real, wilting point                               1    !
!     cmcmax   - real, maximum canopy water parameters             1    !
!     prcp1    - real, effective precip                            1    !
!     zsoil    - real, soil layer depth below ground (negative)  nsoil  !
!     slope    - real, linear reservoir coefficient                1    !
!     frzx     - real, frozen ground parameter                     1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     dksat    - real, saturated soil hydraulic conductivity       1    !
!     dwsat    - real, saturated soil diffusivity                  1    !
!     shdfac   - real, aeral coverage of green veg                 1    !
!     edir1    - real, direct soil evaporation                     1    !
!     ec1      - real, canopy water evaporation                    1    !
!     et1      - real, plant transpiration                       nsoil  !
!                                                                       !
!  input/outputs:                                                       !
!     cmc      - real, canopy moisture content                     1    !
!     sh2o     - real, unfrozen soil moisture                    nsoil  !
!     smc      - real, total soil moisture                       nsoil  !
!                                                                       !
!  outputs:                                                             !
!     runoff1  - real, surface runoff not infiltrating sfc         1    !
!     runoff2  - real, sub surface runoff (baseflow)               1    !
!     runoff3  - real, excess of porosity                          1    !
!     drip     - real, through-fall of precip and/or dew           1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
      integer, intent(in) :: nsoil
 
      real (kind=kind_phys),  intent(in) :: dt, kdt, smcmax, smcwlt,    &
     &       cmcmax, prcp1, slope, frzx, bexp, dksat, dwsat, shdfac,    &
     &       edir1, ec1, et1(nsoil), zsoil(nsoil)
 
!  ---  input/outputs:
      real (kind=kind_phys),  intent(inout) :: cmc, sh2o(nsoil),        &
     &       smc(nsoil)
 
!  ---  outputs:
      real (kind=kind_phys),  intent(out) :: runoff1, runoff2,          &
     &       runoff3, drip
 
!  ---  locals:
      real (kind=kind_phys) :: dummy, excess, pcpdrp, rhsct, trhsct,    &
     &       rhstt(nsold), sice(nsold), sh2oa(nsold), sh2ofg(nsold),    &
     &       ai(nsold), bi(nsold), ci(nsold)
 
      integer :: i, k
!
!===> ...  begin here
!
!  --- ...  executable code begins here.
 
      dummy = 0.0
 
!  --- ...  compute the right hand side of the canopy eqn term ( rhsct )
 
      rhsct = shdfac*prcp1 - ec1
 
!  --- ...  convert rhsct (a rate) to trhsct (an amount) and add it to
!           existing cmc.  if resulting amt exceeds max capacity, it becomes
!           drip and will fall to the grnd.
 
      drip = 0.0
      trhsct = dt * rhsct
      excess = cmc + trhsct
 
      if (excess > cmcmax) drip = excess - cmcmax
 
!  --- ...  pcpdrp is the combined prcp1 and drip (from cmc) that goes into
!           the soil
 
      pcpdrp = (1.0 - shdfac)*prcp1 + drip/dt
 
!  --- ...  store ice content at each soil layer before calling srt & sstep
 
      do i = 1, nsoil
        sice(i) = smc(i) - sh2o(i)
      enddo
 
!  --- ...  call subroutines srt and sstep to solve the soil moisture
!           tendency equations.
 
!  ---  if the infiltrating precip rate is nontrivial,
!         (we consider nontrivial to be a precip total over the time step
!         exceeding one one-thousandth of the water holding capacity of
!         the first soil layer)
!       then call the srt/sstep subroutine pair twice in the manner of
!         time scheme "f" (implicit state, averaged coefficient)
!         of section 2 of kalnay and kanamitsu (1988, mwr, vol 116,
!         pages 1945-1958)to minimize 2-delta-t oscillations in the
!         soil moisture value of the top soil layer that can arise because
!         of the extreme nonlinear dependence of the soil hydraulic
!         diffusivity coefficient and the hydraulic conductivity on the
!         soil moisture state
!       otherwise call the srt/sstep subroutine pair once in the manner of
!         time scheme "d" (implicit state, explicit coefficient)
!         of section 2 of kalnay and kanamitsu
!       pcpdrp is units of kg/m**2/s or mm/s, zsoil is negative depth in m
 
!     if ( pcpdrp .gt. 0.0 ) then
      if ( (pcpdrp*dt) > (0.001*1000.0*(-zsoil(1))*smcmax) ) then
 
!  --- ...  frozen ground version:
!           smc states replaced by sh2o states in srt subr.  sh2o & sice states
!           included in sstep subr.  frozen ground correction factor, frzx
!           added.  all water balance calculations using unfrozen water
 
        call srt                                                        &
!  ---  inputs:
     &     ( nsoil, edir1, et1, sh2o, sh2o, pcpdrp, zsoil, dwsat,       &
     &       dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice,   &
!  ---  outputs:
     &       rhstt, runoff1, runoff2, ai, bi, ci                        &
     &     )
 
        call sstep                                                      &
!  ---  inputs:
     &     ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice,       &
!  ---  input/outputs:
     &       dummy, rhstt, ai, bi, ci,                                  &
!  ---  outputs:
     &       sh2ofg, runoff3, smc                                       &
     &     )
 
        do k = 1, nsoil
          sh2oa(k) = (sh2o(k) + sh2ofg(k)) * 0.5
        enddo
 
        call srt                                                        &
!  ---  inputs:
     &     ( nsoil, edir1, et1, sh2o, sh2oa, pcpdrp, zsoil, dwsat,      &
     &       dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice,   &
!  ---  outputs:
     &       rhstt, runoff1, runoff2, ai, bi, ci                        &
     &     )
 
        call sstep                                                      &
!  ---  inputs:
     &     ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice,       &
!  ---  input/outputs:
     &       cmc, rhstt, ai, bi, ci,                                    &
!  ---  outputs:
     &       sh2o, runoff3, smc                                         &
     &     )
 
      else
 
        call srt                                                        &
!  ---  inputs:
     &     ( nsoil, edir1, et1, sh2o, sh2o, pcpdrp, zsoil, dwsat,       &
     &       dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice,   &
!  ---  outputs:
     &       rhstt, runoff1, runoff2, ai, bi, ci                        &
     &     )
 
        call sstep                                                      &
!  ---  inputs:
     &     ( nsoil, sh2o, rhsct, dt, smcmax, cmcmax, zsoil, sice,       &
!  ---  input/outputs:
     &       cmc, rhstt, ai, bi, ci,                                    &
!  ---  outputs:
     &       sh2o, runoff3, smc                                         &
     &     )
 
      endif
 
!     runof = runoff
!
!...................................
      end subroutine smflx
!-----------------------------------
 
 
!-----------------------------------
      subroutine snowpack                                               &
!  ---  inputs:
     &     ( esd, dtsec, tsnow, tsoil,                                  &
!  ---  input/outputs:
     &       snowh, sndens                                              &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!    subroutine snowpack calculates compaction of snowpack under        !
!    conditions of increasing snow density, as obtained from an         !
!    approximate solution of e. anderson's differential equation (3.29),!
!    noaa technical report nws 19, by victor koren, 03/25/95.           !
!    subroutine will return new values of snowh and sndens              !
!                                                                       !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     esd      - real, water equivalent of snow (m)                1    !
!     dtsec    - real, time step (sec)                             1    !
!     tsnow    - real, snow surface temperature (k)                1    !
!     tsoil    - real, soil surface temperature (k)                1    !
!                                                                       !
!  input/outputs:                                                       !
!     snowh    - real, snow depth (m)                              1    !
!     sndens   - real, snow density                                1    !
!                      (g/cm3=dimensionless fraction of h2o density)    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  parameter constants:
      real (kind=kind_phys), parameter :: c1 = 0.01
      real (kind=kind_phys), parameter :: c2 = 21.0
 
!  ---  inputs:
      real (kind=kind_phys), intent(in) :: esd, dtsec, tsnow, tsoil
 
!  ---  input/outputs:
      real (kind=kind_phys), intent(inout) :: snowh, sndens
 
!  ---  locals:
      real (kind=kind_phys) :: bfac, dsx, dthr, dw, snowhc, pexp,       &
     &       tavgc, tsnowc, tsoilc, esdc, esdcx
 
      integer :: ipol, j
!
!===> ...  begin here
!
!  --- ...  conversion into simulation units
 
      snowhc = snowh * 100.0
      esdc   = esd * 100.0
      dthr   = dtsec / 3600.0
      tsnowc = tsnow - tfreez
      tsoilc = tsoil - tfreez
 
!  --- ...  calculating of average temperature of snow pack
 
      tavgc = 0.5 * (tsnowc + tsoilc)
 
!  --- ...  calculating of snow depth and density as a result of compaction
!           sndens=ds0*(exp(bfac*esd)-1.)/(bfac*esd)
!           bfac=dthr*c1*exp(0.08*tavgc-c2*ds0)
!     note: bfac*esd in sndens eqn above has to be carefully treated
!           numerically below:
!        c1 is the fractional increase in density (1/(cm*hr))
!        c2 is a constant (cm3/g) kojima estimated as 21 cms/g
 
      if (esdc > 1.e-2) then
        esdcx = esdc
      else
        esdcx = 1.e-2
      endif
 
      bfac = dthr*c1 * exp(0.08*tavgc - c2*sndens)
 
!     dsx = sndens * ((dexp(bfac*esdc)-1.0) / (bfac*esdc))
 
!  --- ...  the function of the form (e**x-1)/x imbedded in above expression
!           for dsx was causing numerical difficulties when the denominator "x"
!           (i.e. bfac*esdc) became zero or approached zero (despite the fact
!           that the analytical function (e**x-1)/x has a well defined limit
!           as "x" approaches zero), hence below we replace the (e**x-1)/x
!           expression with an equivalent, numerically well-behaved
!           polynomial expansion.
 
!  --- ...  number of terms of polynomial expansion, and hence its accuracy,
!           is governed by iteration limit "ipol".
!           ipol greater than 9 only makes a difference on double
!           precision (relative errors given in percent %).
!       ipol=9, for rel.error <~ 1.6 e-6 % (8 significant digits)
!       ipol=8, for rel.error <~ 1.8 e-5 % (7 significant digits)
!       ipol=7, for rel.error <~ 1.8 e-4 % ...
 
      ipol = 4
      pexp = 0.0
 
      do j = ipol, 1, -1
!       pexp = (1.0 + pexp)*bfac*esdc /real(j+1)
        pexp = (1.0 + pexp)*bfac*esdcx/real(j+1)
      enddo
      pexp = pexp + 1.
 
      dsx = sndens * pexp
 
!  --- ...  above line ends polynomial substitution
!           end of koren formulation
 
!! --- ...  base formulation (cogley et al., 1990)
!           convert density from g/cm3 to kg/m3
 
!!      dsm = sndens * 1000.0
 
!!      dsx = dsm + dtsec*0.5*dsm*gs2*esd /                             &
!!   &        (1.e7*exp(-0.02*dsm + kn/(tavgc+273.16)-14.643))
 
!! --- ...  convert density from kg/m3 to g/cm3
 
!!      dsx = dsx / 1000.0
 
!! --- ...  end of cogley et al. formulation
 
!  --- ...  set upper/lower limit on snow density
 
      dsx = max( min( dsx, 0.40 ), 0.05 )
      sndens = dsx
 
!  --- ...  update of snow depth and density depending on liquid water
!           during snowmelt.  assumed that 13% of liquid water can be
!           stored in snow per day during snowmelt till snow density 0.40.
 
      if (tsnowc >= 0.0) then
        dw = 0.13 * dthr / 24.0
        sndens = sndens*(1.0 - dw) + dw
        if (sndens > 0.40) sndens = 0.40
      endif
 
!  --- ...  calculate snow depth (cm) from snow water equivalent and snow
!           density. change snow depth units to meters
 
      snowhc = esdc / sndens
      snowh  = snowhc * 0.01
 
!
!...................................
      end subroutine snowpack
!-----------------------------------
 
 
!*********************************************!
!  section-3  3rd or lower level subprograms  !
!*********************************************!
 
 
!-----------------------------------
      subroutine devap                                                  &
!  ---  inputs:
     &     ( etp1, smc, shdfac, smcmax, smcdry, fxexp,                  &
!  ---  outputs:
     &       edir1                                                      &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine devap calculates direct soil evaporation                  !
!                                                                       !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     etp1     - real, potential evaporation                       1    !
!     smc      - real, unfrozen soil moisture                      1    !
!     shdfac   - real, aeral coverage of green vegetation          1    !
!     smcmax   - real, porosity (sat val of soil mois)             1    !
!     smcdry   - real, dry soil mois threshold                     1    !
!     fxexp    - real, bare soil evaporation exponent              1    !
!                                                                       !
!  outputs:                                                             !
!     edir1    - real, direct soil evaporation                     1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
      real (kind=kind_phys), intent(in) :: etp1, smc, shdfac, smcmax,   &
     &       smcdry, fxexp
 
!  ---  outputs:
      real (kind=kind_phys), intent(out) :: edir1
 
!  ---  locals:
      real (kind=kind_phys) :: fx, sratio
!
!===> ...  begin here
!
!  --- ...  direct evap a function of relative soil moisture availability,
!           linear when fxexp=1.
!           fx > 1 represents demand control
!           fx < 1 represents flux control
 
      sratio = (smc - smcdry) / (smcmax - smcdry)
 
      if (sratio > 0.0) then
        fx = sratio**fxexp
        fx = max( min( fx, 1.0 ), 0.0 )
      else
        fx = 0.0
      endif
 
!  --- ...  allow for the direct-evap-reducing effect of shade
 
      edir1 = fx * ( 1.0 - shdfac ) * etp1
!
!...................................
      end subroutine devap
!-----------------------------------
 
 
!-----------------------------------
      subroutine frh2o                                                  &
!  ---  inputs:
     &     ( tkelv, smc, sh2o, smcmax, bexp, psis,                      &
!  ---  outputs:
     &       liqwat                                                     &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine frh2o calculates amount of supercooled liquid soil water  !
!  content if temperature is below 273.15k (t0).  requires newton-type  !
!  iteration to solve the nonlinear implicit equation given in eqn 17   !
!  of koren et al (1999, jgr, vol 104(d16), 19569-19585).               !
!                                                                       !
!  new version (june 2001): much faster and more accurate newton        !
!  iteration achieved by first taking log of eqn cited above -- less    !
!  than 4 (typically 1 or 2) iterations achieves convergence.  also,    !
!  explicit 1-step solution option for special case of parameter ck=0,  !
!  which reduces the original implicit equation to a simpler explicit   !
!  form, known as the "flerchinger eqn". improved handling of solution  !
!  in the limit of freezing point temperature t0.                       !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     tkelv    - real, temperature (k)                             1    !
!     smc      - real, total soil moisture content (volumetric)    1    !
!     sh2o     - real, liquid soil moisture content (volumetric)   1    !
!     smcmax   - real, saturation soil moisture content            1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     psis     - real, saturated soil matric potential             1    !
!                                                                       !
!  outputs:                                                             !
!     liqwat   - real, supercooled liquid water content            1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  constant parameters:
      real (kind=kind_phys), parameter :: ck    = 8.0
!     real (kind=kind_phys), parameter :: ck    = 0.0
      real (kind=kind_phys), parameter :: blim  = 5.5
      real (kind=kind_phys), parameter :: error = 0.005
 
!  ---  inputs:
      real (kind=kind_phys), intent(in) :: tkelv, smc, sh2o, smcmax,    &
     &       bexp, psis
 
!  ---  outputs:
      real (kind=kind_phys), intent(out) :: liqwat
 
!  ---  locals:
      real (kind=kind_phys) :: bx, denom, df, dswl, fk, swl, swlk
 
      integer :: nlog, kcount
!
!===> ...  begin here
!
!  --- ...  limits on parameter b: b < 5.5  (use parameter blim)
!           simulations showed if b > 5.5 unfrozen water content is
!           non-realistically high at very low temperatures.
 
      bx = bexp
      if (bexp > blim)  bx = blim
 
!  --- ...  initializing iterations counter and iterative solution flag.
 
      nlog  = 0
      kcount= 0
 
!  --- ...  if temperature not significantly below freezing (t0), sh2o = smc
 
      if (tkelv > (tfreez-1.e-3)) then
 
        liqwat = smc
 
      else
 
        if (ck /= 0.0) then
 
!  --- ...  option 1: iterated solution for nonzero ck
!                     in koren et al, jgr, 1999, eqn 17
 
!  --- ...  initial guess for swl (frozen content)
 
          swl = smc - sh2o
 
!  --- ...  keep within bounds.
 
          swl = max( min( swl, smc-0.02 ), 0.0 )
 
!  --- ...  start of iterations
 
          do while ( (nlog < 10) .and. (kcount == 0) )
            nlog = nlog + 1
 
            df = log( (psis*gs2/lsubf) * ( (1.0 + ck*swl)**2.0 )        &
     &         * (smcmax/(smc-swl))**bx ) - log(-(tkelv-tfreez)/tkelv)
 
            denom = 2.0*ck/(1.0 + ck*swl) + bx/(smc - swl)
            swlk  = swl - df/denom
 
!  --- ...  bounds useful for mathematical solution.
 
            swlk = max( min( swlk, smc-0.02 ), 0.0 )
 
!  --- ...  mathematical solution bounds applied.
 
            dswl = abs(swlk - swl)
            swl = swlk
 
!  --- ...  if more than 10 iterations, use explicit method (ck=0 approx.)
!           when dswl less or eq. error, no more iterations required.
 
            if ( dswl <= error )  then
              kcount = kcount + 1
            endif
          enddo   !  end do_while_loop
 
!  --- ...  bounds applied within do-block are valid for physical solution.
 
          liqwat = smc - swl
 
        endif   ! end if_ck_block
 
!  --- ...  option 2: explicit solution for flerchinger eq. i.e. ck=0
!                     in koren et al., jgr, 1999, eqn 17
!           apply physical bounds to flerchinger solution
 
        if (kcount == 0) then
          fk = ( ( (lsubf/(gs2*(-psis)))                                &
     &       * ((tkelv-tfreez)/tkelv) )**(-1/bx) ) * smcmax
 
          fk = max( fk, 0.02 )
 
          liqwat = min( fk, smc )
        endif
 
      endif   ! end if_tkelv_block
!
!...................................
      end subroutine frh2o
!-----------------------------------
 
 
!-----------------------------------
      subroutine hrt                                                    &
!  ---  inputs:
     &     ( nsoil, stc, smc, smcmax, zsoil, yy, zz1, tbot,             &
     &       zbot, psisat, dt, bexp, df1, quartz, csoil, vegtyp,        &
     &       shdfac, lheatstrg,                                         &
!  ---  input/outputs:
     &       sh2o,                                                      &
!  ---  outputs:
     &       rhsts, ai, bi, ci                                          &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine hrt calculates the right hand side of the time tendency   !
!  term of the soil thermal diffusion equation.  also to compute        !
!  (prepare) the matrix coefficients for the tri-diagonal matrix of     !
!  the implicit time scheme.                                            !
!                                                                       !
!  subprogram called:  tbnd, snksrc, tmpavg                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     stc      - real, soil temperature                          nsoil  !
!     smc      - real, total soil moisture                       nsoil  !
!     smcmax   - real, porosity                                    1    !
!     zsoil    - real, soil layer depth below ground (negative)  nsoil  !
!     yy       - real,                                             1    !
!     zz1      - real, soil temperture at the top soil column      1    !
!     tbot     - real, bottom soil temp                            1    !
!     zbot     - real, specify depth of lower bd soil              1    !
!     psisat   - real, saturated soil potential                    1    !
!     dt       - real, time step                                   1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     df1      - real, thermal diffusivity                         1    !
!     quartz   - real, soil quartz content                         1    !
!     csoil    - real, soil heat capacity                          1    !
!     vegtyp   - integer, vegetation type                          1    !
!     shdfac   - real, aeral coverage of green vegetation          1    !
!     lheatstrg- logical, flag for canopy heat storage             1    !
!                         parameterization                              !
!                                                                       !
!  input/outputs:                                                       !
!     sh2o     - real, unfrozen soil moisture                    nsoil  !
!                                                                       !
!  outputs:                                                             !
!     rhsts    - real, time tendency of soil thermal diffusion   nsoil  !
!     ai       - real, matrix coefficients                       nsold  !
!     bi       - real, matrix coefficients                       nsold  !
!     ci       - real, matrix coefficients                       nsold  !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
      integer, intent(in) :: nsoil, vegtyp
 
      real (kind=kind_phys),  intent(in) :: stc(nsoil), smc(nsoil),     &
     &       smcmax, zsoil(nsoil), yy, zz1, tbot, zbot, psisat, dt,     &
     &       bexp, df1, quartz, csoil, shdfac
 
      logical, intent(in) :: lheatstrg
 
!  ---  input/outputs:
      real (kind=kind_phys),  intent(inout) :: sh2o(nsoil)
 
!  ---  outputs:
      real (kind=kind_phys),  intent(out) :: rhsts(nsoil), ai(nsold),   &
     &       bi(nsold), ci(nsold)
 
!  ---  locals:
      real (kind=kind_phys) :: ddz, ddz2, denom, df1n, df1k, dtsdz,     &
     &       dtsdz2, hcpct, qtot, ssoil, sice, tavg, tbk, tbk1,         &
     &       tsnsr, tsurf, csoil_loc
 
      integer :: i, k
 
      logical :: itavg
 
!
!===> ...  begin here
!
        csoil_loc=csoil
 
       if (.not.lheatstrg .and. ivegsrc == 1)then
!urban
!
!jhan urban canopy heat storage effect is included in pbl scheme
!
        if( vegtyp == 13 ) then
!           csoil_loc=3.0e6
            csoil_loc=3.0e6*(1.-shdfac)+csoil*shdfac  ! gvf
        endif
       endif
 
!  --- ...  initialize logical for soil layer temperature averaging.
 
      itavg = .true.
!     itavg = .false.
 
!  ===  begin section for top soil layer
 
!  --- ...  calc the heat capacity of the top soil layer
 
      hcpct = sh2o(1)*cph2o2 + (1.0 - smcmax)*csoil_loc                 &
     &      + (smcmax - smc(1))*cp2 + (smc(1) - sh2o(1))*cpice1
 
!  --- ...  calc the matrix coefficients ai, bi, and ci for the top layer
 
      ddz = 1.0 / ( -0.5*zsoil(2) )
      ai(1) = 0.0
      ci(1) = (df1*ddz) / ( zsoil(1)*hcpct )
      bi(1) = -ci(1) + df1 / ( 0.5*zsoil(1)*zsoil(1)*hcpct*zz1 )
 
!  --- ...  calculate the vertical soil temp gradient btwn the 1st and 2nd soil
!           layers.  then calculate the subsurface heat flux. use the temp
!           gradient and subsfc heat flux to calc "right-hand side tendency
!           terms", or "rhsts", for top soil layer.
 
      dtsdz = (stc(1) - stc(2)) / (-0.5*zsoil(2))
      ssoil = df1 * (stc(1) - yy) / (0.5*zsoil(1)*zz1)
      rhsts(1) = (df1*dtsdz - ssoil) / (zsoil(1)*hcpct)
 
!  --- ...  next capture the vertical difference of the heat flux at top and
!           bottom of first soil layer for use in heat flux constraint applied to
!           potential soil freezing/thawing in routine snksrc.
 
      qtot = ssoil - df1*dtsdz
 
!  --- ...  if temperature averaging invoked (itavg=true; else skip):
!           set temp "tsurf" at top of soil column (for use in freezing soil
!           physics later in subroutine snksrc).  if snowpack content is
!           zero, then tsurf expression below gives tsurf = skin temp.  if
!           snowpack is nonzero (hence argument zz1=1), then tsurf expression
!           below yields soil column top temperature under snowpack.  then
!           calculate temperature at bottom interface of 1st soil layer for use
!           later in subroutine snksrc
 
      if (itavg) then
 
        tsurf = (yy + (zz1-1)*stc(1)) / zz1
 
        call tbnd                                                       &
!  ---  inputs:
     &     ( stc(1), stc(2), zsoil, zbot, 1, nsoil,                     &
!  ---  outputs:
     &       tbk                                                        &
     &     )
 
      endif
 
!  --- ...  calculate frozen water content in 1st soil layer.
 
      sice = smc(1) - sh2o(1)
 
!  --- ...  if frozen water present or any of layer-1 mid-point or bounding
!           interface temperatures below freezing, then call snksrc to
!           compute heat source/sink (and change in frozen water content)
!           due to possible soil water phase change
 
      if ( (sice > 0.0) .or. (tsurf < tfreez) .or.                      &
     &     (stc(1) < tfreez) .or. (tbk < tfreez) ) then
 
        if (itavg) then
 
          call tmpavg                                                   &
!  ---  inputs:
     &     ( tsurf, stc(1), tbk, zsoil, nsoil, 1,                       &
!  ---  outputs:
     &       tavg                                                       &
     &     )
 
        else
 
          tavg = stc(1)
 
        endif   ! end if_itavg_block
 
        call snksrc                                                     &
!  ---  inputs:
     &     ( nsoil, 1, tavg, smc(1), smcmax, psisat, bexp, dt,          &
     &       qtot, zsoil,                                               &
!  ---  input/outputs:
     &       sh2o(1),                                                   &
!  ---  outputs:
     &       tsnsr                                                      &
     &     )
 
 
        rhsts(1) = rhsts(1) - tsnsr / ( zsoil(1)*hcpct )
 
      endif   ! end if_sice_block
 
!  ===  this ends section for top soil layer.
 
!  --- ...  initialize ddz2
 
      ddz2 = 0.0
 
!  --- ...  loop thru the remaining soil layers, repeating the above process
!           (except subsfc or "ground" heat flux not repeated in lower layers)
 
      df1k = df1
 
      do k = 2, nsoil
 
!  --- ...  calculate heat capacity for this soil layer.
 
        hcpct = sh2o(k)*cph2o2 + (1.0 - smcmax)*csoil_loc               &
     &        + (smcmax - smc(k))*cp2 + (smc(k) - sh2o(k))*cpice1
 
        if (k /= nsoil) then
 
!  --- ...  this section for layer 2 or greater, but not last layer.
!           calculate thermal diffusivity for this layer.
 
          call tdfcnd                                                   &
!  ---  inputs:
     &     ( smc(k), quartz, smcmax, sh2o(k),                           &
!  ---  outputs:
     &       df1n                                                       &
     &     )
!urban
!     if (ivegsrc == 1)then
!      if ( vegtyp == 13 ) df1n = 3.24
!     endif
!wz only urban for igbp type
!
!jhan urban canopy heat storage effect is included in pbl scheme
!
        if((.not.lheatstrg) .and.
     &      (ivegsrc == 1 .and. vegtyp == 13)) then
          df1n = 3.24*(1.-shdfac) + shdfac*df1n
        endif
 
!  --- ...  calc the vertical soil temp gradient thru this layer
 
          denom = 0.5 * (zsoil(k-1) - zsoil(k+1))
          dtsdz2 = (stc(k) - stc(k+1)) / denom
 
!  --- ...  calc the matrix coef, ci, after calc'ng its partial product
 
          ddz2 = 2.0 / (zsoil(k-1) - zsoil(k+1))
          ci(k) = -df1n*ddz2 / ((zsoil(k-1) - zsoil(k)) * hcpct)
 
!  --- ...  if temperature averaging invoked (itavg=true; else skip):
!           calculate temp at bottom of layer.
 
          if (itavg) then
 
            call tbnd                                                   &
!  ---  inputs:
     &     ( stc(k), stc(k+1), zsoil, zbot, k, nsoil,                   &
!  ---  outputs:
     &       tbk1                                                       &
     &     )
 
          endif
 
        else
 
!  --- ...  special case of bottom soil layer:  calculate thermal diffusivity
!           for bottom layer.
 
          call tdfcnd                                                   &
!  ---  inputs:
     &     ( smc(k), quartz, smcmax, sh2o(k),                           &
!  ---  outputs:
     &       df1n                                                       &
     &     )
!urban
!     if (ivegsrc == 1)then
!      if ( vegtyp == 13 ) df1n = 3.24
!     endif
!wz only urban for igbp type
!
!jhan urban canopy heat storage effect is included in pbl scheme
!
        if((.not.lheatstrg) .and.
     &      (ivegsrc == 1 .and. vegtyp == 13)) then
          df1n = 3.24*(1.-shdfac) + shdfac*df1n
        endif
 
!  --- ...  calc the vertical soil temp gradient thru bottom layer.
 
          denom = 0.5 * (zsoil(k-1) + zsoil(k)) - zbot
          dtsdz2 = (stc(k) - tbot) / denom
 
!  --- ...  set matrix coef, ci to zero if bottom layer.
 
          ci(k) = 0.0
 
!  --- ...  if temperature averaging invoked (itavg=true; else skip):
!           calculate temp at bottom of last layer.
 
          if (itavg) then
 
            call tbnd                                                   &
!  ---  inputs:
     &     ( stc(k), tbot, zsoil, zbot, k, nsoil,                       &
!  ---  outputs:
     &       tbk1                                                       &
     &     )
 
          endif
 
        endif   ! end if_k_block
 
!  --- ...  calculate rhsts for this layer after calc'ng a partial product.
 
        denom = (zsoil(k) - zsoil(k-1)) * hcpct
        rhsts(k) = ( df1n*dtsdz2 - df1k*dtsdz ) / denom
 
        qtot = -1.0 * denom * rhsts(k)
        sice = smc(k) - sh2o(k)
 
        if ( (sice > 0.0) .or. (tbk < tfreez) .or.                      &
     &       (stc(k) < tfreez) .or. (tbk1 < tfreez) ) then
 
          if (itavg) then
 
            call tmpavg                                                 &
!  ---  inputs:
     &     ( tbk, stc(k), tbk1, zsoil, nsoil, k,                        &
!  ---  outputs:
     &       tavg                                                       &
     &     )
 
          else
            tavg = stc(k)
          endif
 
          call snksrc                                                   &
!  ---  inputs:
     &     ( nsoil, k, tavg, smc(k), smcmax, psisat, bexp, dt,          &
     &       qtot, zsoil,                                               &
!  ---  input/outputs:
     &       sh2o(k),                                                   &
!  ---  outputs:
     &       tsnsr                                                      &
     &     )
 
          rhsts(k) = rhsts(k) - tsnsr/denom
        endif
 
!  --- ...  calc matrix coefs, ai, and bi for this layer.
 
        ai(k) = - df1 * ddz / ((zsoil(k-1) - zsoil(k)) * hcpct)
        bi(k) = -(ai(k) + ci(k))
 
!  --- ...  reset values of df1, dtsdz, ddz, and tbk for loop to next soil layer.
 
        tbk   = tbk1
        df1k  = df1n
        dtsdz = dtsdz2
        ddz   = ddz2
 
      enddo   ! end do_k_loop
 
!
!...................................
      end subroutine hrt
!-----------------------------------
 
 
!-----------------------------------
      subroutine hrtice                                                 &
!  ---  inputs:
     &     ( nsoil, stc, zsoil, yy, zz1, df1, ice,                      &
!  ---  input/outputs:
     &       tbot,                                                      &
!  ---  outputs:
     &       rhsts, ai, bi, ci                                          &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine hrtice calculates the right hand side of the time tendency!
!  term of the soil thermal diffusion equation for sea-ice (ice = 1) or !
!  glacial-ice (ice). compute (prepare) the matrix coefficients for the !
!  tri-diagonal matrix of the implicit time scheme.                     !
!  (note:  this subroutine only called for sea-ice or glacial ice, but  !
!  not for non-glacial land (ice = 0).                                  !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     stc      - real, soil temperature                          nsoil  !
!     zsoil    - real, soil depth (negative sign, as below grd)  nsoil  !
!     yy       - real, soil temperature at the top of column       1    !
!     zz1      - real,                                             1    !
!     df1      - real, thermal diffusivity and conductivity        1    !
!     ice      - integer, sea-ice flag (=1: sea-ice, =0: land)     1    !
!                                                                       !
!  input/outputs:                                                       !
!     tbot     - real, bottom soil temperature                     1    !
!                                                                       !
!  outputs:                                                             !
!     rhsts    - real, time tendency of soil thermal diffusion   nsoil  !
!     ai       - real, matrix coefficients                       nsold  !
!     bi       - real, matrix coefficients                       nsold  !
!     ci       - real, matrix coefficients                       nsold  !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
      integer, intent(in) :: nsoil, ice
 
      real (kind=kind_phys), intent(in) :: stc(nsoil), zsoil(nsoil),    &
     &       yy, zz1, df1
 
!  ---  input/outputs:
      real (kind=kind_phys), intent(inout) :: tbot
 
!  ---  outputs:
      real (kind=kind_phys), intent(out) :: rhsts(nsoil), ai(nsold),    &
     &       bi(nsold), ci(nsold)
 
!  ---  locals:
      real (kind=kind_phys) :: ddz, ddz2, denom, dtsdz, dtsdz2,         &
     &       hcpct, ssoil, zbot
 
      integer :: k
 
!
!===> ...  begin here
!
!  --- ...  set a nominal universal value of the sea-ice specific heat capacity,
!           hcpct = 1880.0*917.0 = 1.72396e+6 (source:  fei chen, 1995)
!           set bottom of sea-ice pack temperature: tbot = 271.16
!           set a nominal universal value of glacial-ice specific heat capacity,
!           hcpct = 2100.0*900.0 = 1.89000e+6 (source:  bob grumbine, 2005)
!           tbot passed in as argument, value from global data set
 
      if (ice == 1) then
!  --- ...  sea-ice
        hcpct = 1.72396e+6
        tbot = 271.16
      else
!  --- ...  glacial-ice
        hcpct = 1.89000e+6
      endif
 
!  --- ...  the input argument df1 is a universally constant value of sea-ice
!           and glacial-ice thermal diffusivity, set in sflx as df1 = 2.2.
 
!  --- ...  set ice pack depth.  use tbot as ice pack lower boundary temperature
!           (that of unfrozen sea water at bottom of sea ice pack).  assume ice
!           pack is of n=nsoil layers spanning a uniform constant ice pack
!           thickness as defined by zsoil(nsoil) in routine sflx.
!           if glacial-ice, set zbot = -25 meters
 
      if (ice == 1) then
!  --- ...  sea-ice
        zbot = zsoil(nsoil)
      else
!  --- ...  glacial-ice
        zbot = -25.0
      endif
 
!  --- ...  calc the matrix coefficients ai, bi, and ci for the top layer
 
      ddz = 1.0 / (-0.5*zsoil(2))
      ai(1) = 0.0
      ci(1) = (df1*ddz) / (zsoil(1)*hcpct)
      bi(1) = -ci(1) + df1 / (0.5*zsoil(1)*zsoil(1)*hcpct*zz1)
 
!  --- ...  calc the vertical soil temp gradient btwn the top and 2nd soil
!           layers. recalc/adjust the soil heat flux.  use the gradient and
!           flux to calc rhsts for the top soil layer.
 
      dtsdz = (stc(1) - stc(2)) / (-0.5*zsoil(2))
      ssoil = df1 * (stc(1) - yy) / (0.5*zsoil(1)*zz1)
      rhsts(1) = (df1*dtsdz - ssoil) / (zsoil(1)*hcpct)
 
!  --- ...  initialize ddz2
 
      ddz2 = 0.0
 
!  --- ...  loop thru the remaining soil layers, repeating the above process
 
      do k = 2, nsoil
 
        if (k /= nsoil) then
 
!  --- ...  calc the vertical soil temp gradient thru this layer.
 
          denom = 0.5 * (zsoil(k-1) - zsoil(k+1))
          dtsdz2 = (stc(k) - stc(k+1)) / denom
 
!  --- ...  calc the matrix coef, ci, after calc'ng its partial product.
 
          ddz2 = 2.0 / (zsoil(k-1) - zsoil(k+1))
          ci(k) = -df1*ddz2 / ((zsoil(k-1) - zsoil(k))*hcpct)
 
        else
 
!  --- ...  calc the vertical soil temp gradient thru the lowest layer.
 
          dtsdz2 = (stc(k) - tbot)                                      &
     &           / (0.5*(zsoil(k-1) + zsoil(k)) - zbot)
 
!  --- ...  set matrix coef, ci to zero.
 
          ci(k) = 0.0
 
        endif   ! end if_k_block
 
!  --- ...  calc rhsts for this layer after calc'ng a partial product.
 
        denom = (zsoil(k) - zsoil(k-1)) * hcpct
        rhsts(k) = (df1*dtsdz2 - df1*dtsdz) / denom
 
!  --- ...  calc matrix coefs, ai, and bi for this layer.
 
        ai(k) = - df1*ddz / ((zsoil(k-1) - zsoil(k)) * hcpct)
        bi(k) = -(ai(k) + ci(k))
 
!  --- ...  reset values of dtsdz and ddz for loop to next soil lyr.
 
        dtsdz = dtsdz2
        ddz   = ddz2
 
      enddo   ! end do_k_loop
!
!...................................
      end subroutine hrtice
!-----------------------------------
 
 
!-----------------------------------
      subroutine hstep                                                  &
!  ---  inputs:
     &     ( nsoil, stcin, dt,                                          &
!  ---  input/outputs:
     &       rhsts, ai, bi, ci,                                         &
!  ---  outputs:
     &       stcout                                                     &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine hstep calculates/updates the soil temperature field.      !
!                                                                       !
!  subprogram called:  rosr12                                           !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     stcin    - real, soil temperature                          nsoil  !
!     dt       - real, time step                                   1    !
!                                                                       !
!  input/outputs:                                                       !
!     rhsts    - real, time tendency of soil thermal diffusion   nsoil  !
!     ai       - real, matrix coefficients                       nsold  !
!     bi       - real, matrix coefficients                       nsold  !
!     ci       - real, matrix coefficients                       nsold  !
!                                                                       !
!  outputs:                                                             !
!     stcout   - real, updated soil temperature                  nsoil  !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
      integer, intent(in) :: nsoil
 
      real (kind=kind_phys),  intent(in) :: stcin(nsoil), dt
 
!  ---  input/outputs:
      real (kind=kind_phys),  intent(inout) :: rhsts(nsoil),            &
     &      ai(nsold), bi(nsold), ci(nsold)
 
!  ---  outputs:
      real (kind=kind_phys),  intent(out) :: stcout(nsoil)
 
!  ---  locals:
      integer :: k
 
      real (kind=kind_phys) :: ciin(nsold), rhstsin(nsoil)
 
!
!===> ...  begin here
!
!  --- ...  create finite difference values for use in rosr12 routine
 
      do k = 1, nsoil
        rhsts(k) = rhsts(k) * dt
        ai(k) = ai(k) * dt
        bi(k) = 1.0 + bi(k)*dt
        ci(k) = ci(k) * dt
      enddo
 
!  --- ...  copy values for input variables before call to rosr12
 
      do k = 1, nsoil
         rhstsin(k) = rhsts(k)
      enddo
 
      do k = 1, nsold
        ciin(k) = ci(k)
      enddo
 
!  --- ...  solve the tri-diagonal matrix equation
 
      call rosr12                                                       &
!  ---  inputs:
     &     ( nsoil, ai, bi, rhstsin,                                    &
!  ---  input/outputs:
     &       ciin,                                                      &
!  ---  outputs:
     &       ci, rhsts                                                  &
     &     )
 
!  --- ...  calc/update the soil temps using matrix solution
 
      do k = 1, nsoil
        stcout(k) = stcin(k) + ci(k)
      enddo
!
!...................................
      end subroutine hstep
!-----------------------------------
 
 
!-----------------------------------
      subroutine rosr12                                                 &
!  ---  inputs:
     &     ( nsoil, a, b, d,                                            &
!  ---  input/outputs:
     &       c,                                                         &
!  ---  outputs:
     &       p, delta                                                   &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine rosr12 inverts (solve) the tri-diagonal matrix problem    !
!  shown below:                                                         !
!                                                                       !
! ###                                            ### ###  ###   ###  ###!
! #b(1), c(1),  0  ,  0  ,  0  ,   . . .  ,    0   # #      #   #      #!
! #a(2), b(2), c(2),  0  ,  0  ,   . . .  ,    0   # #      #   #      #!
! # 0  , a(3), b(3), c(3),  0  ,   . . .  ,    0   # #      #   # d(3) #!
! # 0  ,  0  , a(4), b(4), c(4),   . . .  ,    0   # # p(4) #   # d(4) #!
! # 0  ,  0  ,  0  , a(5), b(5),   . . .  ,    0   # # p(5) #   # d(5) #!
! # .                                          .   # #  .   # = #   .  #!
! # .                                          .   # #  .   #   #   .  #!
! # .                                          .   # #  .   #   #   .  #!
! # 0  , . . . , 0 , a(m-2), b(m-2), c(m-2),   0   # #p(m-2)#   #d(m-2)#!
! # 0  , . . . , 0 ,   0   , a(m-1), b(m-1), c(m-1)# #p(m-1)#   #d(m-1)#!
! # 0  , . . . , 0 ,   0   ,   0   ,  a(m) ,  b(m) # # p(m) #   # d(m) #!
! ###                                            ### ###  ###   ###  ###!
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     a        - real, matrix coefficients                       nsoil  !
!     b        - real, matrix coefficients                       nsoil  !
!     d        - real, soil water time tendency                  nsoil  !
!                                                                       !
!  input/outputs:                                                       !
!     c        - real, matrix coefficients                       nsoil  !
!                                                                       !
!  outputs:                                                             !
!     p        - real,                                           nsoil  !
!     delta    - real,                                           nsoil  !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
      integer, intent(in) :: nsoil
 
      real (kind=kind_phys),  dimension(nsoil), intent(in) :: a, b, d
 
!  ---  input/outputs:
      real (kind=kind_phys),  dimension(nsoil), intent(inout) :: c
 
!  ---  outputs:
      real (kind=kind_phys),  dimension(nsoil), intent(out) :: p, delta
 
!  ---  locals:
      integer :: k, kk
 
!
!===> ...  begin here
!
!  --- ...  initialize eqn coef c for the lowest soil layer
 
      c(nsoil) = 0.0
 
!  --- ...  solve the coefs for the 1st soil layer
 
      p(1) = -c(1) / b(1)
      delta(1) = d(1) / b(1)
 
!  --- ...  solve the coefs for soil layers 2 thru nsoil
 
      do k = 2, nsoil
        p(k) = -c(k) * ( 1.0 / (b(k) + a(k)*p(k-1)) )
        delta(k) = (d(k) - a(k)*delta(k-1))                              &
     &           * ( 1.0 / (b(k) + a(k)*p(k-1)) )
      enddo
 
!  --- ...  set p to delta for lowest soil layer
 
      p(nsoil) = delta(nsoil)
 
!  --- ...  adjust p for soil layers 2 thru nsoil
 
      do k = 2, nsoil
         kk = nsoil - k + 1
         p(kk) = p(kk)*p(kk+1) + delta(kk)
      enddo
!
!...................................
      end subroutine rosr12
!-----------------------------------
 
 
!-----------------------------------
      subroutine snksrc                                                 &
!  ---  inputs:
     &     ( nsoil, k, tavg, smc, smcmax, psisat, bexp, dt,             &
     &       qtot, zsoil,                                               &
!  ---  input/outputs:
     &       sh2o,                                                      &
!  ---  outputs:
     &       tsrc                                                       &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!  subroutine snksrc calculates sink/source term of the termal          !
!  diffusion equation.                                                  !
!                                                                       !
!  subprograms called: frh2o                                            !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     k        - integer, index of soil layers                     1    !
!     tavg     - real, soil layer average temperature              1    !
!     smc      - real, total soil moisture                         1    !
!     smcmax   - real, porosity                                    1    !
!     psisat   - real, saturated soil potential                    1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     dt       - real, time step                                   1    !
!     qtot     - real, tot vertical diff of heat flux              1    !
!     zsoil    - real, soil layer depth below ground (negative)  nsoil  !
!                                                                       !
!  input/outputs:                                                       !
!     sh2o     - real, available liqued water                      1    !
!                                                                       !
!  outputs:                                                             !
!     tsrc     - real, heat source/sink                            1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  parameter constants:
      real (kind=kind_phys), parameter :: dh2o = 1.0000e3
 
!  ---  inputs:
      integer, intent(in) :: nsoil, k
 
      real (kind=kind_phys), intent(in) :: tavg, smc, smcmax, psisat,   &
     &       bexp, dt, qtot, zsoil(nsoil)
 
!  ---  input/outputs:
      real (kind=kind_phys), intent(inout) :: sh2o
 
!  ---  outputs:
      real (kind=kind_phys), intent(out) :: tsrc
 
!  ---  locals:
      real (kind=kind_phys) :: dz, free, xh2o
 
!  ---  external functions:
!     real (kind=kind_phys) :: frh2o
!
!===> ...  begin here
!
      if (k == 1) then
        dz = -zsoil(1)
      else
        dz = zsoil(k-1) - zsoil(k)
      endif
 
!  --- ...  via function frh2o, compute potential or 'equilibrium' unfrozen
!           supercooled free water for given soil type and soil layer temperature.
!           function frh20 invokes eqn (17) from v. koren et al (1999, jgr, vol.
!           104, pg 19573).  (aside:  latter eqn in journal in centigrade units.
!           routine frh2o use form of eqn in kelvin units.)
 
!     free = frh2o( tavg,smc,sh2o,smcmax,bexp,psisat )
 
      call frh2o                                                        &
!  ---  inputs:
     &     ( tavg, smc, sh2o, smcmax, bexp, psisat,                     &
!  ---  outputs:
     &       free                                                       &
     &     )
 
 
!  --- ...  in next block of code, invoke eqn 18 of v. koren et al (1999, jgr,
!           vol. 104, pg 19573.)  that is, first estimate the new amountof liquid
!           water, 'xh2o', implied by the sum of (1) the liquid water at the begin
!           of current time step, and (2) the freeze of thaw change in liquid
!           water implied by the heat flux 'qtot' passed in from routine hrt.
!           second, determine if xh2o needs to be bounded by 'free' (equil amt) or
!           if 'free' needs to be bounded by xh2o.
 
      xh2o = sh2o + qtot*dt / (dh2o*lsubf*dz)
 
!  --- ...  first, if freezing and remaining liquid less than lower bound, then
!           reduce extent of freezing, thereby letting some or all of heat flux
!           qtot cool the soil temp later in routine hrt.
 
      if ( xh2o < sh2o .and. xh2o < free) then
        if ( free > sh2o ) then
          xh2o = sh2o
        else
          xh2o = free
        endif
      endif
 
!  --- ...  second, if thawing and the increase in liquid water greater than
!           upper bound, then reduce extent of thaw, thereby letting some or
!           all of heat flux qtot warm the soil temp later in routine hrt.
 
      if ( xh2o > sh2o .and. xh2o > free )  then
        if ( free < sh2o ) then
          xh2o = sh2o
        else
          xh2o = free
        endif
      endif
 
      xh2o = max( min( xh2o, smc ), 0.0 )
 
!  --- ...  calculate phase-change heat source/sink term for use in routine hrt
!           and update liquid water to reflcet final freeze/thaw increment.
 
      tsrc = -dh2o * lsubf * dz * (xh2o - sh2o) / dt
      sh2o = xh2o
!
!...................................
      end subroutine snksrc
!-----------------------------------
 
 
!-----------------------------------
      subroutine srt                                                    &
!  ---  inputs:
     &     ( nsoil, edir, et, sh2o, sh2oa, pcpdrp, zsoil, dwsat,        &
     &       dksat, smcmax, bexp, dt, smcwlt, slope, kdt, frzx, sice,   &
!  ---  outputs:
     &       rhstt, runoff1, runoff2, ai, bi, ci                        &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!    subroutine srt calculates the right hand side of the time tendency !
!    term of the soil water diffusion equation.  also to compute        !
!    ( prepare ) the matrix coefficients for the tri-diagonal matrix    !
!    of the implicit time scheme.                                       !
!                                                                       !
!  subprogram called:  wdfcnd                                           !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     edir     - real, direct soil evaporation                     1    !
!     et       - real, plant transpiration                       nsoil  !
!     sh2o     - real, unfrozen soil moisture                    nsoil  !
!     sh2oa    - real,                                           nsoil  !
!     pcpdrp   - real, combined prcp and drip                      1    !
!     zsoil    - real, soil layer depth below ground             nsoil  !
!     dwsat    - real, saturated soil diffusivity                  1    !
!     dksat    - real, saturated soil hydraulic conductivity       1    !
!     smcmax   - real, porosity                                    1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     dt       - real, time step                                   1    !
!     smcwlt   - real, wilting point                               1    !
!     slope    - real, linear reservoir coefficient                1    !
!     kdt      - real,                                             1    !
!     frzx     - real, frozen ground parameter                     1    !
!     sice     - real, ice content at each soil layer            nsoil  !
!                                                                       !
!  outputs:                                                             !
!     rhstt    - real, soil water time tendency                  nsoil  !
!     runoff1  - real, surface runoff not infiltrating sfc         1    !
!     runoff2  - real, sub surface runoff (baseflow)               1    !
!     ai       - real, matrix coefficients                       nsold  !
!     bi       - real, matrix coefficients                       nsold  !
!     ci       - real, matrix coefficients                       nsold  !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  inputs:
      integer, intent(in) :: nsoil
 
      real (kind=kind_phys), dimension(nsoil), intent(in) :: et,        &
     &       sh2o, sh2oa, zsoil, sice
 
      real (kind=kind_phys), intent(in) :: edir, pcpdrp, dwsat, dksat,  &
     &       smcmax, smcwlt, bexp, dt, slope, kdt, frzx
 
!  --- outputs:
      real (kind=kind_phys), intent(out) :: runoff1, runoff2,           &
     &      rhstt(nsoil), ai(nsold), bi(nsold), ci(nsold)
 
 
!  ---  locals:
      real (kind=kind_phys) :: acrt, dd, ddt, ddz, ddz2, denom, denom2, &
     &       dice, dsmdz, dsmdz2, dt1, fcr, infmax, mxsmc, mxsmc2, px,  &
     &       numer, pddum, sicemax, slopx, smcav, sstt, sum, val, wcnd, &
     &       wcnd2, wdf, wdf2, dmax(nsold)
 
      integer :: cvfrz, ialp1, iohinf, j, jj, k, ks
!
!===> ...  begin here
!
!  --- ...  frozen ground version:
!           reference frozen ground parameter, cvfrz, is a shape parameter
!           of areal distribution function of soil ice content which equals
!           1/cv. cv is a coefficient of spatial variation of soil ice content.
!           based on field data cv depends on areal mean of frozen depth, and
!           it close to constant = 0.6 if areal mean frozen depth is above 20 cm.
!           that is why parameter cvfrz = 3 (int{1/0.6*0.6}). current logic
!           doesn't allow cvfrz be bigger than 3
 
        parameter(cvfrz = 3)
 
c ----------------------------------------------------------------------
!  --- ...  determine rainfall infiltration rate and runoff.  include
!           the infiltration formule from schaake and koren model.
!           modified by q duan
 
      iohinf = 1
 
!  --- ... let sicemax be the greatest, if any, frozen water content within
!          soil layers.
 
      sicemax = 0.0
      do ks = 1, nsoil
       if (sice(ks) > sicemax) sicemax = sice(ks)
      enddo
 
!  --- ...  determine rainfall infiltration rate and runoff
 
      pddum = pcpdrp
      runoff1 = 0.0
 
      if (pcpdrp /= 0.0) then
 
!  --- ...  modified by q. duan, 5/16/94
 
        dt1 = dt/86400.
        smcav = smcmax - smcwlt
        dmax(1) = -zsoil(1) * smcav
 
!  --- ...  frozen ground version:
 
        dice = -zsoil(1) * sice(1)
 
        dmax(1) = dmax(1)*(1.0 - (sh2oa(1)+sice(1)-smcwlt)/smcav)
        dd = dmax(1)
 
        do ks = 2, nsoil
 
!  --- ...  frozen ground version:
 
          dice = dice + ( zsoil(ks-1) - zsoil(ks) ) * sice(ks)
 
          dmax(ks) = (zsoil(ks-1)-zsoil(ks))*smcav
          dmax(ks) = dmax(ks)*(1.0 - (sh2oa(ks)+sice(ks)-smcwlt)/smcav)
          dd = dd + dmax(ks)
        enddo
 
!  --- ...  val = (1.-exp(-kdt*sqrt(dt1)))
!           in below, remove the sqrt in above
 
        val = 1.0 - exp(-kdt*dt1)
        ddt = dd * val
 
        px = pcpdrp * dt
        if (px < 0.0) px = 0.0
 
        infmax = (px*(ddt/(px+ddt)))/dt
 
!  --- ...  frozen ground version:
!           reduction of infiltration based on frozen ground parameters
 
        fcr = 1.
 
        if (dice > 1.e-2) then
          acrt = cvfrz * frzx / dice
          sum = 1.
 
          ialp1 = cvfrz - 1
          do j = 1, ialp1
            k = 1
 
            do jj = j+1,ialp1
              k = k * jj
            enddo
 
            sum = sum + (acrt**( cvfrz-j)) / float(k)
          enddo
 
          fcr = 1.0 - exp(-acrt) * sum
        endif
 
        infmax = infmax * fcr
 
!  --- ...  correction of infiltration limitation:
!           if infmax .le. hydrolic conductivity assign infmax the value
!           of hydrolic conductivity
 
!       mxsmc = max ( sh2oa(1), sh2oa(2) )
        mxsmc = sh2oa(1)
 
        call wdfcnd                                                     &
!  ---  inputs:
     &     ( mxsmc, smcmax, bexp, dksat, dwsat, sicemax,                &
!  ---  outputs:
     &       wdf, wcnd                                                  &
     &     )
 
        infmax = max( infmax, wcnd )
        infmax = min( infmax, px )
 
        if (pcpdrp > infmax) then
          runoff1 = pcpdrp - infmax
          pddum   = infmax
        endif
 
      endif   ! end if_pcpdrp_block
 
!  --- ... to avoid spurious drainage behavior, 'upstream differencing'
!          in line below replaced with new approach in 2nd line:
!          'mxsmc = max(sh2oa(1), sh2oa(2))'
 
      mxsmc = sh2oa(1)
 
      call wdfcnd                                                       &
!  ---  inputs:
     &     ( mxsmc, smcmax, bexp, dksat, dwsat, sicemax,                &
!  ---  outputs:
     &       wdf, wcnd                                                  &
     &     )
 
!  --- ...  calc the matrix coefficients ai, bi, and ci for the top layer
 
      ddz = 1.0 / ( -.5*zsoil(2) )
      ai(1) = 0.0
      bi(1) = wdf * ddz / ( -zsoil(1) )
      ci(1) = -bi(1)
 
!  --- ...  calc rhstt for the top layer after calc'ng the vertical soil
!           moisture gradient btwn the top and next to top layers.
 
      dsmdz = ( sh2o(1) - sh2o(2) ) / ( -.5*zsoil(2) )
      rhstt(1) = (wdf*dsmdz + wcnd - pddum + edir + et(1)) / zsoil(1)
      sstt = wdf * dsmdz + wcnd + edir + et(1)
 
!  --- ...  initialize ddz2
 
      ddz2 = 0.0
 
!  --- ...  loop thru the remaining soil layers, repeating the abv process
 
      do k = 2, nsoil
        denom2 = (zsoil(k-1) - zsoil(k))
 
        if (k /= nsoil) then
          slopx = 1.0
 
!  --- ...  again, to avoid spurious drainage behavior, 'upstream differencing'
!           in line below replaced with new approach in 2nd line:
!           'mxsmc2 = max (sh2oa(k), sh2oa(k+1))'
 
          mxsmc2 = sh2oa(k)
 
          call wdfcnd                                                   &
!  ---  inputs:
     &     ( mxsmc2, smcmax, bexp, dksat, dwsat, sicemax,               &
!  ---  outputs:
     &       wdf2, wcnd2                                                &
     &     )
 
!  --- ...  calc some partial products for later use in calc'ng rhstt
 
          denom = (zsoil(k-1) - zsoil(k+1))
          dsmdz2 = (sh2o(k) - sh2o(k+1)) / (denom * 0.5)
 
!  --- ...  calc the matrix coef, ci, after calc'ng its partial product
 
          ddz2 = 2.0 / denom
          ci(k) = -wdf2 * ddz2 / denom2
 
        else   ! if_k_block
 
!  --- ...  slope of bottom layer is introduced
 
          slopx = slope
 
!  --- ...  retrieve the soil water diffusivity and hydraulic conductivity
!           for this layer
 
          call wdfcnd                                                   &
!  ---  inputs:
     &     ( sh2oa(nsoil), smcmax, bexp, dksat, dwsat, sicemax,         &
!  ---  outputs:
     &       wdf2, wcnd2                                                &
     &     )
 
!  --- ...  calc a partial product for later use in calc'ng rhstt
          dsmdz2 = 0.0
 
!  --- ...  set matrix coef ci to zero
 
          ci(k) = 0.0
 
        endif   ! end if_k_block
 
!  --- ...  calc rhstt for this layer after calc'ng its numerator
 
        numer = wdf2*dsmdz2 + slopx*wcnd2 - wdf*dsmdz - wcnd + et(k)
        rhstt(k) = numer / (-denom2)
 
!  --- ...  calc matrix coefs, ai, and bi for this layer
 
        ai(k) = -wdf * ddz / denom2
        bi(k) = -( ai(k) + ci(k) )
 
!  --- ...  reset values of wdf, wcnd, dsmdz, and ddz for loop to next lyr
!      runoff2:  sub-surface or baseflow runoff
 
        if (k == nsoil) then
          runoff2 = slopx * wcnd2
        endif
 
        if (k /= nsoil) then
          wdf  = wdf2
          wcnd = wcnd2
          dsmdz= dsmdz2
          ddz  = ddz2
        endif
      enddo   ! end do_k_loop
!
!...................................
      end subroutine srt
!-----------------------------------
 
 
!-----------------------------------
      subroutine sstep                                                  &
!  ---  inputs:
     &     ( nsoil, sh2oin, rhsct, dt, smcmax, cmcmax, zsoil, sice,     &
!  ---  input/outputs:
     &       cmc, rhstt, ai, bi, ci,                                    &
!  ---  outputs:
     &       sh2oout, runoff3, smc                                      &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!    subroutine sstep calculates/updates soil moisture content values   !
!    and canopy moisture content values.                                !
!                                                                       !
!  subprogram called:  rosr12                                           !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     sh2oin   - real, unfrozen soil moisture                    nsoil  !
!     rhsct    - real,                                             1    !
!     dt       - real, time step                                   1    !
!     smcmax   - real, porosity                                    1    !
!     cmcmax   - real, maximum canopy water parameters             1    !
!     zsoil    - real, soil layer depth below ground             nsoil  !
!     sice     - real, ice content at each soil layer            nsoil  !
!                                                                       !
!  input/outputs:                                                       !
!     cmc      - real, canopy moisture content                     1    !
!     rhstt    - real, soil water time tendency                  nsoil  !
!     ai       - real, matrix coefficients                       nsold  !
!     bi       - real, matrix coefficients                       nsold  !
!     ci       - real, matrix coefficients                       nsold  !
!                                                                       !
!  outputs:                                                             !
!     sh2oout  - real, updated soil moisture content             nsoil  !
!     runoff3  - real, excess of porosity                          1    !
!     smc      - real, total soil moisture                       nsoil  !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  input:
      integer, intent(in) :: nsoil
 
      real (kind=kind_phys), dimension(nsoil), intent(in) :: sh2oin,    &
     &       zsoil, sice
 
      real (kind=kind_phys), intent(in) :: rhsct, dt, smcmax, cmcmax
 
!  ---  inout/outputs:
      real (kind=kind_phys), intent(inout) :: cmc, rhstt(nsoil),        &
     &       ai(nsold), bi(nsold), ci(nsold)
 
!  ---  outputs:
      real (kind=kind_phys), intent(out) :: sh2oout(nsoil), runoff3,    &
     &       smc(nsoil)
 
!  ---  locals:
      real (kind=kind_phys) :: ciin(nsold), rhsttin(nsoil), ddz, stot,  &
     &       wplus
 
      integer :: i, k, kk11
!
!===> ...  begin here
!
!  --- ...  create 'amount' values of variables to be input to the
!           tri-diagonal matrix routine.
 
      do k = 1, nsoil
        rhstt(k) = rhstt(k) * dt
        ai(k) = ai(k) * dt
        bi(k) = 1. + bi(k) * dt
        ci(k) = ci(k) * dt
      enddo
 
!  --- ...  copy values for input variables before call to rosr12
 
      do k = 1, nsoil
        rhsttin(k) = rhstt(k)
      enddo
 
      do k = 1, nsold
        ciin(k) = ci(k)
      enddo
 
!  --- ...  call rosr12 to solve the tri-diagonal matrix
 
      call rosr12                                                       &
!  ---  inputs:
     &     ( nsoil, ai, bi, rhsttin,                                    &
!  ---  input/outputs:
     &       ciin,                                                      &
!  ---  outputs:
     &       ci, rhstt                                                  &
     &     )
 
!  --- ...  sum the previous smc value and the matrix solution to get
!           a new value.  min allowable value of smc will be 0.02.
!      runoff3: runoff within soil layers
 
      wplus   = 0.0
      runoff3 = 0.0
      ddz     = -zsoil(1)
 
      do k = 1, nsoil
        if (k /= 1) ddz = zsoil(k - 1) - zsoil(k)
 
        sh2oout(k) = sh2oin(k) + ci(k) + wplus/ddz
 
        stot = sh2oout(k) + sice(k)
        if (stot > smcmax) then
          if (k == 1) then
            ddz = -zsoil(1)
          else
            kk11 = k - 1
            ddz = -zsoil(k) + zsoil(kk11)
          endif
 
          wplus = (stot - smcmax) * ddz
        else
          wplus = 0.0
        endif
 
        smc(k) = max( min( stot, smcmax ), 0.02 )
        sh2oout(k) = max( smc(k)-sice(k), 0.0 )
      enddo
 
      runoff3 = wplus
 
!  --- ...  update canopy water content/interception (cmc).  convert rhsct to
!           an 'amount' value and add to previous cmc value to get new cmc.
 
      cmc = cmc + dt*rhsct
      if (cmc < 1.e-20) cmc = 0.0
      cmc = min( cmc, cmcmax )
!
!...................................
      end subroutine sstep
!-----------------------------------
 
 
!-----------------------------------
      subroutine tbnd                                                   &
!  ---  inputs:
     &     ( tu, tb, zsoil, zbot, k, nsoil,                             &
!  ---  outputs:
     &       tbnd1                                                      &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!                                                                       !
!    subroutine tbnd calculates temperature on the boundary of the      !
!    layer by interpolation of the middle layer temperatures            !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     tu       - real, soil temperature                            1    !
!     tb       - real, bottom soil temp                            1    !
!     zsoil    - real, soil layer depth                          nsoil  !
!     zbot     - real, specify depth of lower bd soil              1    !
!     k        - integer, soil layer index                         1    !
!     nsoil    - integer, number of soil layers                    1    !
!                                                                       !
!  outputs:                                                             !
!     tbnd1    - real, temperature at bottom interface of the lyr  1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  input:
      integer, intent(in) :: k, nsoil
 
      real (kind=kind_phys), intent(in) :: tu, tb, zbot, zsoil(nsoil)
 
!  ---  output:
      real (kind=kind_phys), intent(out) :: tbnd1
 
!  ---  locals:
      real (kind=kind_phys) :: zb, zup
 
!  --- ...  use surface temperature on the top of the first layer
 
      if (k == 1) then
        zup = 0.0
      else
        zup = zsoil(k-1)
      endif
 
!  --- ...  use depth of the constant bottom temperature when interpolate
!           temperature into the last layer boundary
 
      if (k == nsoil) then
        zb = 2.0*zbot - zsoil(k)
      else
        zb = zsoil(k+1)
      endif
 
!  --- ...  linear interpolation between the average layer temperatures
 
      tbnd1 = tu + (tb-tu)*(zup-zsoil(k))/(zup-zb)
!
!...................................
      end subroutine tbnd
!-----------------------------------
 
 
!-----------------------------------
      subroutine tmpavg                                                 &
!  ---  inputs:
     &     ( tup, tm, tdn, zsoil, nsoil, k,                             &
!  ---  outputs:
     &       tavg                                                       &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!    subroutine tmpavg calculates soil layer average temperature (tavg) !
!    in freezing/thawing layer using up, down, and middle layer         !
!    temperatures (tup, tdn, tm), where tup is at top boundary of       !
!    layer, tdn is at bottom boundary of layer.  tm is layer prognostic !
!    state temperature.                                                 !
!                                                                       !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     tup      - real, temperature ar top boundary of layer        1    !
!     tm       - real, layer prognostic state temperature          1    !
!     tdn      - real, temperature ar bottom boundary of layer     1    !
!     zsoil    - real, soil layer depth                          nsoil  !
!     nsoil    - integer, number of soil layers                    1    !
!     k        - integer, layer index                              1    !
!  outputs:                                                             !
!     tavg     - real, soil layer average temperature              1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  input:
      integer, intent(in) :: nsoil, k
 
      real (kind=kind_phys), intent(in) :: tup, tm, tdn, zsoil(nsoil)
 
!  ---  output:
      real (kind=kind_phys), intent(out) :: tavg
 
!  ---  locals:
      real (kind=kind_phys) :: dz, dzh, x0, xdn, xup
!
!===> ...  begin here
!
      if (k == 1) then
        dz = -zsoil(1)
      else
        dz = zsoil(k-1) - zsoil(k)
      endif
 
      dzh = dz * 0.5
 
      if (tup < tfreez) then
 
        if (tm < tfreez) then
          if (tdn < tfreez) then               ! tup, tm, tdn < t0
            tavg = (tup + 2.0*tm + tdn) / 4.0
          else                                 ! tup & tm < t0,  tdn >= t0
            x0 = (tfreez - tm) * dzh / (tdn - tm)
            tavg = 0.5*(tup*dzh + tm*(dzh+x0)+tfreez*(2.*dzh-x0)) / dz
          endif
        else
          if (tdn < tfreez) then               ! tup < t0, tm >= t0, tdn < t0
            xup  = (tfreez-tup) * dzh / (tm-tup)
            xdn  = dzh - (tfreez-tm) * dzh / (tdn-tm)
            tavg = 0.5*(tup*xup + tfreez*(2.*dz-xup-xdn)+tdn*xdn) / dz
          else                                 ! tup < t0, tm >= t0, tdn >= t0
            xup  = (tfreez-tup) * dzh / (tm-tup)
            tavg = 0.5*(tup*xup + tfreez*(2.*dz-xup)) / dz
          endif
        endif
 
      else    ! if_tup_block
 
        if (tm < tfreez) then
          if (tdn < tfreez) then               ! tup >= t0, tm < t0, tdn < t0
            xup  = dzh - (tfreez-tup) * dzh / (tm-tup)
            tavg = 0.5*(tfreez*(dz-xup) + tm*(dzh+xup)+tdn*dzh) / dz
          else                                 ! tup >= t0, tm < t0, tdn >= t0
            xup  = dzh - (tfreez-tup) * dzh / (tm-tup)
            xdn  = (tfreez-tm) * dzh / (tdn-tm)
            tavg = 0.5 * (tfreez*(2.*dz-xup-xdn) + tm*(xup+xdn)) / dz
          endif
        else
          if (tdn < tfreez) then               ! tup >= t0, tm >= t0, tdn < t0
            xdn  = dzh - (tfreez-tm) * dzh / (tdn-tm)
            tavg = (tfreez*(dz-xdn) + 0.5*(tfreez+tdn)*xdn) / dz
          else                                 ! tup >= t0, tm >= t0, tdn >= t0
            tavg = (tup + 2.0*tm + tdn) / 4.0
          endif
        endif
 
      endif   ! end if_tup_block
!
!...................................
      end subroutine tmpavg
!-----------------------------------
 
 
!-----------------------------------
      subroutine transp                                                 &
!  ---  inputs:
     &     ( nsoil, nroot, etp1, smc, smcwlt, smcref,                   &
     &       cmc, cmcmax, zsoil, shdfac, pc, cfactr, rtdis,             &
!  ---  outputs:
     &       et1                                                        &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!     subroutine transp calculates transpiration for the veg class.     !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     nsoil    - integer, number of soil layers                    1    !
!     nroot    - integer, number of root layers                    1    !
!     etp1     - real, potential evaporation                       1    !
!     smc      - real, unfrozen soil moisture                    nsoil  !
!     smcwlt   - real, wilting point                               1    !
!     smcref   - real, soil mois threshold                         1    !
!     cmc      - real, canopy moisture content                     1    !
!     cmcmax   - real, maximum canopy water parameters             1    !
!     zsoil    - real, soil layer depth below ground             nsoil  !
!     shdfac   - real, aeral coverage of green vegetation          1    !
!     pc       - real, plant coeff                                 1    !
!     cfactr   - real, canopy water parameters                     1    !
!     rtdis    - real, root distribution                         nsoil  !
!                                                                       !
!  outputs:                                                             !
!     et1      - real, plant transpiration                       nsoil  !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  input:
      integer, intent(in) :: nsoil, nroot
 
      real (kind=kind_phys), intent(in) :: etp1, smcwlt, smcref,        &
     &       cmc, cmcmax, shdfac, pc, cfactr
 
      real (kind=kind_phys), dimension(nsoil), intent(in) :: smc,       &
     &       zsoil, rtdis
 
!  ---  output:
      real (kind=kind_phys), dimension(nsoil), intent(out) :: et1
 
!  ---  locals:
      real (kind=kind_phys) :: denom, etp1a, rtx, sgx, gx(7)
 
      integer :: i, k
!
!===> ...  begin here
!
!  --- ...  initialize plant transp to zero for all soil layers.
 
      do k = 1, nsoil
        et1(k) = 0.0
      enddo
 
!  --- ...  calculate an 'adjusted' potential transpiration
!           if statement below to avoid tangent linear problems near zero
!           note: gx and other terms below redistribute transpiration by layer,
!           et(k), as a function of soil moisture availability, while preserving
!           total etp1a.
 
      if (cmc /= 0.0) then
        etp1a = shdfac * pc * etp1 * (1.0 - (cmc /cmcmax) ** cfactr)
      else
        etp1a = shdfac * pc * etp1
      endif
 
      sgx = 0.0
      do i = 1, nroot
        gx(i) = ( smc(i) - smcwlt ) / ( smcref - smcwlt )
        gx(i) = max( min( gx(i), 1.0 ), 0.0 )
        sgx = sgx + gx(i)
      enddo
      sgx = sgx / nroot
 
      denom = 0.0
      do i = 1, nroot
        rtx = rtdis(i) + gx(i) - sgx
        gx(i) = gx(i) * max( rtx, 0.0 )
        denom = denom + gx(i)
      enddo
      if (denom <= 0.0) denom = 1.0
 
      do i = 1, nroot
        et1(i) = etp1a * gx(i) / denom
      enddo
 
!  --- ...  above code assumes a vertically uniform root distribution
!           code below tests a variable root distribution
 
!     et(1) = ( zsoil(1) / zsoil(nroot) ) * gx * etp1a
!     et(1) = ( zsoil(1) / zsoil(nroot) ) * etp1a
 
!  --- ...  using root distribution as weighting factor
 
!     et(1) = rtdis(1) * etp1a
!     et(1) = etp1a * part(1)
 
!  --- ...  loop down thru the soil layers repeating the operation above,
!           but using the thickness of the soil layer (rather than the
!           absolute depth of each layer) in the final calculation.
 
!     do k = 2, nroot
!       gx = ( smc(k) - smcwlt ) / ( smcref - smcwlt )
!       gx = max ( min ( gx, 1.0 ), 0.0 )
!  --- ...  test canopy resistance
!       gx = 1.0
!       et(k) = ((zsoil(k)-zsoil(k-1))/zsoil(nroot))*gx*etp1a
!       et(k) = ((zsoil(k)-zsoil(k-1))/zsoil(nroot))*etp1a
 
!  --- ...  using root distribution as weighting factor
 
!       et(k) = rtdis(k) * etp1a
!       et(k) = etp1a*part(k)
!     enddo
 
!
!...................................
      end subroutine transp
!-----------------------------------
 
 
!-----------------------------------
      subroutine wdfcnd                                                 &
!  ---  inputs:
     &     ( smc, smcmax, bexp, dksat, dwsat, sicemax,                  &
!  ---  outputs:
     &       wdf, wcnd                                                  &
     &     )
 
! ===================================================================== !
!  description:                                                         !
!     subroutine wdfcnd calculates soil water diffusivity and soil      !
!     hydraulic conductivity.                                           !
!                                                                       !
!  subprogram called:  none                                             !
!                                                                       !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
!  inputs:                                                       size   !
!     smc      - real, layer total soil moisture                   1    !
!     smcmax   - real, porosity                                    1    !
!     bexp     - real, soil type "b" parameter                     1    !
!     dksat    - real, saturated soil hydraulic conductivity       1    !
!     dwsat    - real, saturated soil diffusivity                  1    !
!     sicemax  - real, max frozen water content in soil layer      1    !
!                                                                       !
!  outputs:                                                             !
!     wdf      - real, soil water diffusivity                      1    !
!     wcnd     - real, soil hydraulic conductivity                 1    !
!                                                                       !
!  ====================    end of description    =====================  !
!
!  ---  input:
      real (kind=kind_phys), intent(in)  :: smc, smcmax, bexp, dksat,   &
     &       dwsat, sicemax
 
!  ---  output:
      real (kind=kind_phys), intent(out) :: wdf, wcnd
 
!  ---  locals:
      real (kind=kind_phys) :: expon, factr1, factr2, vkwgt
!
!===> ...  begin here
!
!  --- ...  calc the ratio of the actual to the max psbl soil h2o content
 
      factr1 = min(1.0, max(0.0, 0.2/smcmax))
      factr2 = min(1.0, max(0.0, smc/smcmax))
 
!  --- ...  prep an expntl coef and calc the soil water diffusivity
 
      expon = bexp + 2.0
      wdf = dwsat * factr2 ** expon
 
!  --- ...  frozen soil hydraulic diffusivity.  very sensitive to the vertical
!           gradient of unfrozen water. the latter gradient can become very
!           extreme in freezing/thawing situations, and given the relatively
!           few and thick soil layers, this gradient sufferes serious
!           trunction errors yielding erroneously high vertical transports of
!           unfrozen water in both directions from huge hydraulic diffusivity.
!           therefore, we found we had to arbitrarily constrain wdf
!
!           version d_10cm:  .......  factr1 = 0.2/smcmax
!           weighted approach.......  pablo grunmann, 28_sep_1999.
 
      if (sicemax > 0.0) then
        vkwgt = 1.0 / (1.0 + (500.0*sicemax)**3.0)
        wdf = vkwgt*wdf + (1.0- vkwgt)*dwsat*factr1**expon
      endif
 
!  --- ...  reset the expntl coef and calc the hydraulic conductivity
 
      expon = (2.0 * bexp) + 3.0
      wcnd = dksat * factr2 ** expon
!
!...................................
      end subroutine wdfcnd
!-----------------------------------
 
! =========================== !
!     end contain programs    !
! =========================== !
 
!...................................
      end subroutine gfssflx
!-----------------------------------
      end module sflx