SUBROUTINE sflx(ice,satflg,dt,z,  &
     f,soldn,sfcprs,prcp,sfctmp,th2,q2,q2sat,dqsdt2,tbot,ch,chkff,  &
     vegtyp,soityp,shdfac,  &
     etp,eta,h,s,runoff1,runoff2,q1,snmax,  &
     t1,cmc,smc,stc,snodep,soilw,soilm)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !    purpose:  to calculate and accumulate evapotranspiration values
  !    =======   based upon a soil hydrology model from oregon state
  !              university. the routine also calculates soil temps
  !              and moisture values at two soil depths.
  !
  !some universal veg/soil parameters
  IMPLICIT NONE
  REAL, PARAMETER       :: cfactr=0.5
  REAL, PARAMETER       :: cmcmax=0.5e-3
  REAL, PARAMETER       :: topt=298.0
  REAL, PARAMETER       :: rsmax=5000.0
  INTEGER, PARAMETER    :: nroot=2
  INTEGER, PARAMETER    :: nsoil=2
  INTEGER, PARAMETER    :: nsold=2

  INTEGER, INTENT(in out)  :: ice, vegtyp,soityp
  REAL, INTENT(in out)     :: satflg,z,f,soldn,q2sat,dqsdt2
  REAL, INTENT(in out)     :: tbot,shdfac,cmc, snodep
  REAL, INTENT(in)         :: dt, sfcprs,prcp, sfctmp, th2 
  REAL, INTENT(in)         :: q2, ch,chkff , t1, smc(nsold), stc(nsold)
  REAL, INTENT(out)        :: etp,eta,h,s,runoff1,runoff2,q1
  REAL, INTENT(out)        :: snmax, soilw, soilm

  LOGICAL :: snowng
  LOGICAL :: frzgra
  REAL :: b, df1, dksat, dwsat, epsca, esd, f1, pc, rch, rr
  REAL :: smcdry, smcmax,smcref,smcwlt,t1min,t1v,t24,t2v,th2v,z0
  REAL :: zsoil  ( nsold )
  REAL :: beta,drip,ec,edir,ett,flx1,flx2,flx3,rho,runof, dew,rib
  REAL :: runoff3,albedo,rcmin,rgl,hs,prcp1,rc,edir1,ec1
  REAL :: ett1,heat,t14,soidif,soilwm,soilww
  INTEGER :: kz,k
  COMMON/rite/ beta,drip,ec,edir,ett,flx1,flx2,flx3,rho,runof, dew,rib

  !   initialization
  runoff1 = 0.0
  runoff2 = 0.0
  runoff3 = 0.0
  snmax = 0.0

  !   specify the veg/soil parameters
  CALL redprm(vegtyp,soityp,  &
       cmcmax,topt,rsmax,albedo,z0,shdfac,nroot,rcmin,rgl,hs,  &
       b,dksat,dwsat,smcmax,smcwlt,smcref,smcdry,f1)

  IF(ice == 1) THEN
     DO kz = 1, nsoil
        zsoil(kz)=-3.*float(kz)/float(nsoil)
     END DO
  ELSE
     zsoil(1)=-0.1
     zsoil(2)=-2.0
  END IF

  !executable code begins here...initialize misc variables.
  s = 0.0
  snowng = .FALSE.
  frzgra = .FALSE.

  !convert the input snow depth (meters) from a snowboard type
  !measurement to a liquid equivalent snow depth (esd, in meters)
  !using a 10 to 1 ratio.
  !      esd = snodep * 0.1
  esd = snodep
  IF(ice == 1) THEN
     esd = 0.01
  END IF

  !calc virtual temps and potential temps at grnd (sub 1) and at
  !the 1st mdl lvl abv the grnd (sub 2). expon is cp divd by r.
  t1v = t1 * (1.0 + 0.61 * q2 )
  t2v = sfctmp * ( 1.0 + 0.61 * q2 )
  !      th2 = sfctmp + ( 0.0097545 * z )
  th2v = th2 * (1.0 + 0.61 * q2 )

  !determine if it's preciping and what kind of precip it is.
  !if it's prcping and the air temp is colder than 0 c, it's snowing!
  !if it's prcping and the air temp is warmer than 0 c, but the grnd
  !temp is colder than 0 c, freezing rain is presumed to be falling.
  IF ( prcp > 0.0 ) THEN
     IF ( sfctmp <= 273.16 ) THEN
        snowng = .TRUE.
     ELSE
        IF ( sfctmp <= 273.16 ) frzgra = .TRUE.
     END IF
  END IF

  !if either prcp flag is set, convert the prcp rate from
  !kg m-2 s-1 to a liquid equiv snow depth (meters) and accum it.
  IF ( ( snowng ) .OR. ( frzgra ) ) THEN
     esd =esd + prcp * dt * 0.001
     prcp1 = 0.0
  ELSE
     prcp1 = prcp
  END IF

  !check: if esd < 1.e-6 m. the snow cover will be ignored.
  IF(esd < 1.e-6) THEN
     snmax=esd
     esd = 0.0
     snowng = .FALSE.
  END IF

  !calc grnd heat flux for snow pack case (will effectively
  !calc the heat flux from the snow pack).  remember, if
  !t1 > 273.16 k, there will be snow melt and the skin temp will
  !be 273.16 k. const (0.35) is thermal diffusivity for snow.
  IF ( esd > 0.0 ) THEN
     t1min = amin1 ( t1, 273.16 )
     s = 0.35 * ( t1min - stc(1) ) / ( 10.0 * esd )
     !set the minimum flux within snow to -200 w m-2
     !for numerical reasons
     s = MIN ( MAX ( s,-100.0 ), 100.0 )
  ELSE
     !calculate grnd heat flux for case with no snow pack
     CALL tdfcnd ( df1, smc(1), b, f1 ) 
     s = df1 * ( stc(1) - t1 ) / ( 0.5 * zsoil(1) )
  END IF

  !call penman subroutine to calculate potential evaporation
  !and other partial products and sums for later calculations.
  CALL penman ( sfctmp,sfcprs,ch,t2v,th2,prcp,f,t24,s,q2,  &
       q2sat,etp,rch,epsca,rr,snowng,frzgra,dqsdt2)
  IF(satflg <= 0.5) etp = 0.0

  !call canres to calculate the canopy resistance and convert it
  !into pc
  IF(shdfac > 1.e-6) THEN
     CALL canres(soldn,ch,sfctmp,q2,sfcprs,smc,zsoil,nsoil,  &
          smcwlt,smcref,rcmin,rc,pc,nroot,q2sat,dqsdt2, topt,rsmax,rgl,hs)
  END IF

  !if point is to be fully processed....
  !based on whether or not there is a snow pack, call snopac or
  !nopac to calculate an estimate of actual evapotranspiration,
  !update soil moisture, and update soil temperature.
  IF ( esd == 0.0 ) THEN
     CALL nopac ( etp, eta, prcp, smc, smcmax, smcwlt,  &
          smcref,smcdry, cmc, cmcmax, nsoil, dt, shdfac,  &
          q1,q2,t1,sfctmp,t24,th2,f,f1,s,stc, epsca, b, pc, rch, rr,  cfactr,  &
          zsoil, dksat, dwsat, tbot,runoff1,runoff2,  &
          runoff3, edir1, ec1, ett1,nroot,ice)
  ELSE
     CALL snopac ( etp,eta,prcp,prcp1,snowng,smc,smcmax,smcwlt,  &
          smcref, smcdry, cmc, cmcmax, nsoil, dt, q1,  &
          q2,t1,sfctmp,t24,th2,f,f1,s,stc,epsca,sfcprs,  &
          b, pc, rch, rr, cfactr, esd, zsoil, dwsat, dksat, tbot, shdfac,runoff1,  &
          runoff2,runoff3,edir1,ec1,ett1,nroot,snmax,ice)
  END IF

  !reconvert the liquid equivalent snow depth (meters) to an
  !output snow depth (meters) using a 10 to 1 ratio.
  !snodep = esd * 10.0
  snodep = esd

  !calc sub-products then use to calc sensible heat flux.
  heat = -chkff * ( th2 - t1 )
  t14 = t1 * t1 * t1 * t1
  rr = t14 * 6.48e-8 / ( sfcprs * ch ) + 1.0
  h = heat

  !  convert eta from kg m-2 s-1 to w m-2
  eta = eta*2.5e+6
  etp = etp*2.5e+6

  ! convert the sign of soil heat flux so that:
  !         s>0: warm the surface  (night time)
  !         s<0: cool the surface  (day time)
  s=-1.0*s

  !  convert runoff3 from m to m s-1
  runoff3 = runoff3/dt
  runoff2 = runoff2+runoff3
  soilm=-1.0*smc(1)*zsoil(1)
  DO k = 2, nsoil
     soilm=soilm+smc(k)*(zsoil(k-1)-zsoil(k))
  END DO
  soidif = MIN(smcdry, smcwlt)
  soilwm=-1.0*(smcmax-soidif)*zsoil(1)
  soilww=-1.0*(smc(1)-soidif)*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))
  END DO
  soilw=soilww/soilwm
END SUBROUTINE sflx


SUBROUTINE canres(solar,ch,sfctmp,q2,sfcprs,smc,zsoil,nsoil,  &
     smcwlt,smcref,rcmin,rc,pc,nroot,q2sat,dqsdt2, topt,rsmax,rgl,hs)
  !this routine calculates the canopy resistance which depends
  !on incoming solar radiation, air temperature, atmospheric
  !water vapor pressure deficit at the lowest model level, ans soil
  !moisture.
  ! ----------------------------------------------------------------------
  ! source:  jarvis (1976), jacquemin and noilhan (1990 blm)
  ! ----------------------------------------------------------------------
  ! ----------------------------------------------------------------------
  ! input:   solar: incoming solar radiation
  !                ch: surface exchange coefficient for heat and moisture
  !                sfctmp: air temperature at 1st level
  !                q2:  air humidity at 1st level
  !                sfcprs: surface pressure
  !                smc: volumetric soil moisture
  !                zsoil: soil depth
  !                nsoil: number of soil layers
  !                smcwlt: wilting point
  !                smcref: reference soil moisture (field capacity)

  ! output:  pc: plant coefficient
  !                rc: canopy resistance
  ! ----------------------------------------------------------------------
  ! **********************************************************************
  IMPLICIT NONE
  REAL, INTENT(in)         :: solar
  REAL, INTENT(in out)     :: ch
  REAL, INTENT(in)         :: sfctmp
  REAL, INTENT(in out)     :: q2
  REAL, INTENT(in)         :: sfcprs
  REAL, INTENT(in out)     :: smc(nsoil)
  REAL, INTENT(in)         :: zsoil(nsoil)
  INTEGER, INTENT(in out)  :: nsoil
  REAL, INTENT(in out)     :: smcwlt
  REAL, INTENT(in out)     :: smcref
  REAL, INTENT(in)         :: rcmin
  REAL, INTENT(out)        :: rc
  REAL, INTENT(out)        :: pc
  INTEGER, INTENT(in)      :: nroot
  REAL, INTENT(in)         :: q2sat
  REAL, INTENT(in)         :: dqsdt2
  REAL, INTENT(in out)     :: topt
  REAL, INTENT(in out)     :: rsmax
  REAL, INTENT(in)         :: rgl
  REAL, INTENT(in)         :: hs
  INTEGER, PARAMETER :: nsold = 2
  REAL, DIMENSION(nsold) :: part
  REAL, PARAMETER :: sigma=5.67e-8
  REAL, PARAMETER :: rd=287.0
  REAL, PARAMETER :: cp=1004.5
  REAL, PARAMETER :: slv=2.5e6
  REAL :: rcs,rct,rcq,rcsoil,ff,p,qs,gx,tair4,st1,slvcp
  REAL :: rr,delta
  INTEGER :: k

  ! for low crops use these values from jacquemin and noilhan (90 blm)
  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*solar/rgl
  rcs = (ff + rcmin/rsmax) / (1.0 + ff)
  rcs = amax1(rcs,0.0001)

  ! contribution due to the temperature at the first model level.
  !
  !if (topt.lt.sfctmp) print*,' > in sflx,topt < sfctmp ',topt,sfctmp
  !!(r1=topt-sfctmp,  r1 nao pode ser negativo para sx4
  !!fisicamente e' permitido topt < sfctmp?  nao esquecer de checar isto
  !!rct = 1.0 - 0.0016*((topt-sfctmp)**2.0)
  !
  rct = 1.0 - 0.0016*((topt-sfctmp)*(topt-sfctmp))
  rct = amax1(rct,0.0001)

  ! contribution due to vapor pressure deficit at the first model level.
  p = sfcprs
  qs = q2sat
  ! rcq expression from ssib
  rcq = 1.0/(1.0+hs*(qs-q2))
  rcq = amax1(rcq,0.01)

  ! contribution due to soil moisture availability.
  ! determine contribution from each soil layer, then add them up.
  gx = (smc(1) - smcwlt) / (smcref - smcwlt)
  IF (gx > 1.) gx = 1.
  IF (gx < 0.) gx = 0.
  part(1) = (zsoil(1)/zsoil(nroot)) * gx

  DO  k = 2, nroot
     gx = (smc(k) - smcwlt) / (smcref - smcwlt)
     IF (gx > 1.) gx = 1.
     IF (gx < 0.) gx = 0.
     part(k) = ((zsoil(k)-zsoil(k-1))/zsoil(nroot)) * gx
  END DO
  DO  k = 1, nroot
     rcsoil = rcsoil+part(k)
  END DO
  rcsoil = amax1(rcsoil,0.0001)

  ! determine canopy resistance due to all factors.
  ! convert canopy resistance (rc) to plant coefficient (pc).
  rc = rcmin/(rcs*rct*rcq*rcsoil)
  tair4 = sfctmp**4
  st1 = (4.*sigma*rd)/cp
  slvcp = slv/cp
  rr = st1*tair4/(sfcprs*ch) + 1.0
  delta = slvcp*dqsdt2

  pc = (rr+delta)/(rr*(1.+rc*ch)+delta)
END SUBROUTINE canres


FUNCTION devap ( etp1, smc, zsoil, shdfac, smcmax, b,  &
     dksat, dwsat, smcdry, smcref, smcwlt)
  !name:  direct evaporation (devap) function  version: n/a


  IMPLICIT NONE
  REAL, INTENT(in)                         :: etp1
  REAL, INTENT(in out)                     :: smc
  REAL, INTENT(in out)                     :: zsoil
  REAL, INTENT(in)                         :: shdfac
  REAL, INTENT(in out)                     :: smcmax
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(in out)                     :: dksat
  REAL, INTENT(in out)                     :: dwsat
  REAL, INTENT(in out)                     :: smcdry
  REAL, INTENT(in out)                     :: smcref
  REAL, INTENT(in out)                     :: smcwlt

  REAL :: devap
  REAL :: fx
  REAL :: wcnd
  REAL :: wdf

  !     get soil h2o diffusivity and hydraulic conductivity given smc
  CALL wdfcnd ( wdf, wcnd, smc, smcmax, b, dksat, dwsat )

  !fx represents soil moisture flux / atmospheric demand
  !fx = ( -2. * wdf * ( smc - smcdry ) / zsoil - wcnd ) / etp1
  !*** test of beta methode, f. chen 7/11/96 ***
  fx = (smc - smcwlt) / (smcref - smcwlt)

  !fx > 1 represents demand control
  !fx < 1 represents flux control
  fx = MAX ( MIN ( fx, 1. ) ,0. )

  !allow for the direct-evap-reducing effect of shade
  devap = fx * ( 1.0 - shdfac ) * etp1
END FUNCTION devap

SUBROUTINE hrt ( rhsts, stc, smc, smcmax, nsoil, zsoil, yy, zz1,  &
     tbot, b, f1, df1 )
  !c    purpose:  to calculate the right hand side of the time tendency
  !c    =======   term of the soil thermal diffusion equation.  also to
  !c              compute ( prepare ) the matrix coefficients for the
  !c              tri-diagonal matrix of the implicit time scheme.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: rhsts( nsoil )
  REAL, INTENT(in)                         :: stc( nsoil )
  REAL, INTENT(in)                         :: smc( nsoil )
  REAL, INTENT(in out)                     :: smcmax
  INTEGER, INTENT(in)                      :: nsoil
  REAL, INTENT(in)                         :: zsoil( nsoil )
  REAL, INTENT(in)                         :: yy
  REAL, INTENT(in)                         :: zz1
  REAL, INTENT(in out)                     :: tbot
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(in out)                     :: f1
  REAL, INTENT(in)                         :: df1
  INTEGER, PARAMETER :: nsold = 2

  REAL :: ai    ( nsold )
  REAL :: bi    ( nsold )
  REAL :: ci    ( nsold )

  REAL :: ddz
  REAL :: ddz2
  REAL :: denom
  REAL :: dfmax1
  REAL :: dfmaxn
  REAL :: dtsdz
  REAL :: dtsdz2
  REAL :: hcpct
  INTEGER :: k
  REAL :: mxsmc
  REAL :: s

  COMMON /abci/ ai, bi, ci
  REAL, PARAMETER :: cair=1004.0
  REAL, PARAMETER :: ch2o=4.2e6
  REAL, PARAMETER :: csoil=1.26e6
  REAL, PARAMETER :: zbot=-3.0

  ! retrieve the thermal diffusivity for the wetter of the top soil
  ! layer or the next lower soil layer
  mxsmc = MAX ( smc(1), smc(2) )
  CALL tdfcnd ( dfmax1, mxsmc, b, f1 )

  ! calc the heat capacity of the top soil layer
  hcpct = smc(1)*ch2o + (1.0-smcmax)*csoil + (smcmax-smc(1))*cair

  !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) =  ( dfmax1 * 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) )
  s = df1 * ( stc(1) - yy ) / ( 0.5 * zsoil(1) * zz1 )
  rhsts(1) = ( dfmax1 * dtsdz - s ) / ( zsoil(1) * hcpct )

  ! initialize ddz2
  ddz2 = 0.0

  ! loop thru the remaining soil layers, repeating the above process
  DO  k = 2, nsoil

     !calc this soil layer's heat capacity
     hcpct = smc(k)*ch2o +(1.0-smcmax)*csoil +(smcmax-smc(k))*cair
     IF ( k /= nsoil ) THEN

        !retrieve the thermal diffusivity for the wetter of the
        !current or the next lower soil layer
        mxsmc = MAX ( smc(k), smc(k+1) )
        CALL tdfcnd ( dfmaxn, mxsmc, b, f1 )
        !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. / (zsoil(k-1) - zsoil(k+1))
        ci(k) = -dfmaxn * ddz2 / ((zsoil(k-1) - zsoil(k)) * hcpct)

     ELSE

        !retrieve the thermal diffusivity for the lowest soil layer
        CALL tdfcnd ( dfmaxn, smc(k), b, f1 )
        !calc the vertical soil temp gradient thru the lowest layer
        dtsdz2 = (stc(k)-tbot) / (.5 * (zsoil(k-1) + zsoil(k))-zbot)
        !set matrix coef, ci to zero
        ci(k) = 0.

     END IF

     !calc rhsts for this layer after calc'ng a partial product
     denom = ( zsoil(k) - zsoil(k-1) ) * hcpct
     rhsts(k) = ( dfmaxn * dtsdz2 - dfmax1 * dtsdz ) / denom
     !calc matrix coefs, ai, and bi for this layer.
     ai(k) = - dfmax1 * ddz / ((zsoil(k-1) - zsoil(k)) * hcpct)
     bi(k) = -(ai(k) + ci(k))
     !reset values of dfmax1, dtsdz, and ddz for loop to next soil lyr
     dfmax1 = dfmaxn
     dtsdz = dtsdz2
     ddz = ddz2
  END DO
END SUBROUTINE hrt

SUBROUTINE hrtice (rhsts,stc,smc,smcmax,nsoil,zsoil,yy,zz1, tbot, b, f1, df1 )
  !    purpose:  to calculate 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.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: rhsts( nsoil )
  REAL, INTENT(in)                         :: stc( nsoil )
  REAL, INTENT(in out)                     :: smc( nsoil )
  REAL, INTENT(in out)                     :: smcmax
  INTEGER, INTENT(in)                      :: nsoil
  REAL, INTENT(in)                         :: zsoil( nsoil )
  REAL, INTENT(in)                         :: yy
  REAL, INTENT(in)                         :: zz1
  REAL, INTENT(in out)                     :: tbot
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(in out)                     :: f1
  REAL, INTENT(in)                         :: df1
  INTEGER, PARAMETER :: nsold = 2

  REAL :: ai    ( nsold )
  REAL :: bi    ( nsold )
  REAL :: ci    ( nsold )
  REAL :: ddz
  REAL :: ddz2
  REAL :: denom
  REAL :: dfmax1
  REAL :: dfmaxn
  REAL :: dtsdz
  REAL :: dtsdz2
  REAL :: hcpct
  INTEGER :: k
  REAL :: s
  REAL :: zbot
  COMMON /abci/ ai, bi, ci
  REAL, PARAMETER :: cair=1004.0
  REAL, PARAMETER :: ch2o=4.2e6
  REAL, PARAMETER :: csoil=1.26e6

  zbot=zsoil(nsoil)
  dfmax1=2.2
  hcpct=1880.0*917.0

  !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) =  ( dfmax1 * 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) )
  s = df1 * ( stc(1) - yy ) / ( 0.5 * zsoil(1) * zz1 )
  rhsts(1) = ( dfmax1 * dtsdz - s ) / ( 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
        dfmaxn=2.2
        !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. / (zsoil(k-1) - zsoil(k+1))
        ci(k) = -dfmaxn * ddz2 / ((zsoil(k-1) - zsoil(k)) * hcpct)
     ELSE
        dfmaxn=2.2
        !calc the vertical soil temp gradient thru the lowest layer
        dtsdz2 = (stc(k)-tbot)/(.5 * (zsoil(k-1) + zsoil(k))-zbot)
        !set matrix coef, ci to zero
        ci(k) = 0.
     END IF
     !calc rhsts for this layer after calc'ng a partial product
     denom = ( zsoil(k) - zsoil(k-1) ) * hcpct
     rhsts(k) = ( dfmaxn * dtsdz2 - dfmax1 * dtsdz ) / denom
     !calc matrix coefs, ai, and bi for this layer.
     ai(k) = - dfmax1 * ddz / ((zsoil(k-1) - zsoil(k)) * hcpct)
     bi(k) = -(ai(k) + ci(k))
     !reset values of dfmax1, dtsdz, and ddz for loop to next soil lyr
     dfmax1 = dfmaxn
     dtsdz = dtsdz2
     ddz = ddz2
  END DO
END SUBROUTINE hrtice

SUBROUTINE hstep ( stc, rhsts, dt, nsoil )
  ! purpose:  to calculate/update the soil temperature field.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: stc( nsoil )
  REAL, INTENT(out)                        :: rhsts( nsoil )
  REAL, INTENT(in)                         :: dt
  INTEGER, INTENT(in)                      :: nsoil
  INTEGER, PARAMETER :: nsold=2
  INTEGER :: k
  REAL :: ai    ( nsold )
  REAL :: bi    ( nsold )
  REAL :: ci    ( nsold )
  COMMON /abci/ ai, bi, ci

  !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. + bi(k) * dt
     ci(k) = ci(k) * dt
  END DO

  !solve the tri-diagonal matrix equation
  CALL rosr12 ( ci,ai,bi,ci,rhsts,rhsts,nsoil )

  !calc/update the soil temps using matrix solution
  DO  k = 1 , nsoil
     stc(k) = stc(k) + ci(k)
  END DO
END SUBROUTINE hstep


SUBROUTINE nopac ( etp, eta, prcp, smc, smcmax, smcwlt,  &
     smcref,smcdry,cmc,cmcmax, nsoil, dt, shdfac,  &
     q1, q2, t1, sfctmp, t24, th2, f, f1, s, stc,  &
     epsca, b, pc, rch, rr,  cfactr, zsoil, dksat,  &
     dwsat, tbot, runoff1, runoff2,runoff3, edir1, ec1, ett1, nroot, ice)
  !    purpose:  to calculate 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.

  IMPLICIT NONE
  REAL, INTENT(in)                         :: etp
  REAL, INTENT(out)                        :: eta
  REAL, INTENT(in)                         :: prcp
  REAL, INTENT(in out)                     :: smc( nsoil )
  REAL, INTENT(in out)                     :: smcmax
  REAL, INTENT(in out)                     :: smcwlt
  REAL, INTENT(in out)                     :: smcref
  REAL, INTENT(in out)                     :: smcdry
  REAL, INTENT(in out)                     :: cmc
  REAL, INTENT(in out)                     :: cmcmax
  INTEGER, INTENT(in out)                  :: nsoil
  REAL, INTENT(in out)                     :: dt
  REAL, INTENT(in out)                     :: shdfac
  REAL, INTENT(in out)                     :: q1
  REAL, INTENT(in out)                     :: q2
  REAL, INTENT(in out)                     :: t1
  REAL, INTENT(in)                         :: sfctmp
  REAL, INTENT(in)                         :: t24
  REAL, INTENT(in)                         :: th2
  REAL, INTENT(in)                         :: f
  REAL, INTENT(in out)                     :: f1
  REAL, INTENT(in out)                     :: s
  REAL, INTENT(in out)                     :: stc( nsoil )
  REAL, INTENT(in out)                     :: epsca
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(in out)                     :: pc
  REAL, INTENT(in)                         :: rch
  REAL, INTENT(in)                         :: rr
  REAL, INTENT(in out)                     :: cfactr
  REAL, INTENT(in)                         :: zsoil( nsoil )
  REAL, INTENT(in out)                     :: dksat
  REAL, INTENT(in out)                     :: dwsat
  REAL, INTENT(in out)                     :: tbot
  REAL, INTENT(in out)                     :: runoff1
  REAL, INTENT(in out)                     :: runoff2
  REAL, INTENT(in out)                     :: runoff3
  REAL, INTENT(in out)                     :: edir1
  REAL, INTENT(in out)                     :: ec1
  REAL, INTENT(in out)                     :: ett1
  INTEGER, INTENT(in out)                  :: nroot
  INTEGER, INTENT(in out)                  :: ice

  REAL :: beta
  REAL :: dew
  REAL :: df1
  REAL :: drip
  REAL :: ec
  REAL :: edir
  REAL :: eta1
  REAL :: etp1
  REAL :: ett
  REAL :: flx1
  REAL :: flx2
  REAL :: flx3
  REAL :: prcp1
  REAL :: rho
  REAL :: rib
  REAL :: runoff
  REAL :: yy
  REAL :: yynum
  REAL :: zz1
  COMMON/rite/ beta,drip,ec,edir,ett,flx1,flx2,flx3,rho,runoff, dew,rib
  REAL, PARAMETER :: cp=1004.5
  REAL, PARAMETER :: sigma=5.67e-8

  !     executable code begins 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
  IF ( etp > 0.0 ) THEN
     !convert prcp from  kg m-2 s-1  to  m s-1
     CALL smflx ( eta1,smc,nsoil,cmc,etp1,dt,prcp1,zsoil,  &
          smcmax,b,pc,smcwlt,dksat,dwsat,smcref,shdfac,  &
          cmcmax,smcdry,cfactr, runoff1,runoff2, runoff3,  &
          edir1, ec1, ett1, sfctmp,q2,nroot)
     !convert modeled evapotranspiration fm  m s-1  to  kg m-2 s-1
     eta = eta1 * 1000.0
  ELSE
     !if etp < 0, assume dew forms (transform etp1 into dew
     !and reinitialize etp1 to zero)
     dew = -etp1
     etp1 = 0.0
     !convert prcp from  kg m-2 s-1  to  m s-1  and add dew amt
     prcp1 = prcp1 + dew
     CALL smflx ( eta1,smc,nsoil,cmc,etp1,dt,prcp1,zsoil,  &
          smcmax,b,pc,smcwlt,dksat,dwsat,smcref,shdfac,  &
          cmcmax,smcdry,cfactr, runoff1,runoff2, runoff3,  &
          edir1, ec1, ett1, sfctmp, q2, nroot)
     !convert modeled evapotranspiration fm  m s-1  to  kg m-2 s-1
     eta = eta1 * 1000.0
  END IF

  !based on etp and e values, determine beta
  IF ( etp <= 0.0 ) THEN
     beta = 0.0
     IF ( etp < 0.0 ) THEN
        beta = 1.0
        eta = etp
     END IF
  ELSE
     beta = eta / etp
  END IF

  ! calc mixing ratio at grnd level (skin)
  !q1 = q2 + eta * cp / rch

  !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 ( df1, smc(1), b, f1 )
  yynum = f - sigma * t24
  yy = sfctmp + (yynum/rch+th2-sfctmp-beta*epsca) / rr
  zz1 = df1 / ( -0.5 * zsoil(1) * rch * rr ) + 1
  CALL shflx ( s,stc,smc,smcmax,nsoil,t1,dt,yy,zz1,zsoil,tbot, b,f1,df1, ice )
  !set flx1, and flx3 to zero since they are not used.  flx2
  !was similarly initialized in the penman routine.
  flx1 = 0.0
  flx3 = 0.0
END SUBROUTINE nopac

SUBROUTINE penman(sfctmp,sfcprs,ch,t2v,th2,prcp,f,t24,s,q2,  &
     q2sat,etp,rch,epsca,rr,snowng,frzgra,dqsdt2)
  !    purpose:  to calculate potential evaporation for the current point.
  !    =======   various partial sums/products are also calculated and
  !              passed back to the calling routine for later use.

  IMPLICIT NONE
  REAL, INTENT(in)                         :: sfctmp
  REAL, INTENT(in)                         :: sfcprs
  REAL, INTENT(in)                         :: ch
  REAL, INTENT(in)                         :: t2v
  REAL, INTENT(in)                         :: th2
  REAL, INTENT(in)                         :: prcp
  REAL, INTENT(in)                         :: f
  REAL, INTENT(out)                        :: t24
  REAL, INTENT(in)                         :: s
  REAL, INTENT(in)                         :: q2
  REAL, INTENT(in)                         :: q2sat
  REAL, INTENT(out)                        :: etp
  REAL, INTENT(out)                        :: rch
  REAL, INTENT(out)                        :: epsca
  REAL, INTENT(out)                        :: rr
  LOGICAL, INTENT(in)                      :: snowng
  LOGICAL, INTENT(in)                      :: frzgra
  REAL, INTENT(in)                         :: dqsdt2
  REAL :: t,tt
  REAL :: a
  REAL :: beta
  REAL :: delta
  REAL :: dew
  REAL :: drip
  REAL :: ec
  REAL :: edir
  REAL :: ett
  REAL :: flx1
  REAL :: flx2
  REAL :: flx3
  REAL :: fnet
  REAL :: rad
  REAL :: rho
  REAL :: rib
  REAL :: runoff
  COMMON/rite/ beta,drip,ec,edir,ett,flx1,flx2,flx3,rho,runoff, dew,rib
  REAL, PARAMETER :: cp=1004.6
  REAL, PARAMETER :: cph2o=4.218e+3
  REAL, PARAMETER :: cpice=2.106e+3
  REAL, PARAMETER :: r=287.04
  REAL, PARAMETER :: elcp=2.4888e+3
  REAL, PARAMETER :: lsubf=3.335e+5
  REAL, PARAMETER :: lsubc=2.5e+6
  REAL, PARAMETER :: sigma=5.67e-8

  !executable code begins here...define function.
  tt(t) = t * t
  flx2 = 0.0

  !prepare partial quantities for penman equation.
  delta = elcp * dqsdt2
  t24 = tt( sfctmp ) * tt( sfctmp )
  rr = t24 * 6.48e-8 / ( sfcprs * ch ) + 1.0
  rho = sfcprs / ( r * 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 + cph2o * prcp / rch
  ELSE
     rr = rr + cpice * prcp / rch
  END IF
  fnet = f - sigma * t24 - s

  !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
  END IF

  !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 rosr12 ( p, a, b, c, d, delta, nsoil )
  !    purpose:  to invert (solve) the tri-diagonal matrix problem shown
  !    =======   below:

  IMPLICIT NONE
  REAL, INTENT(out)                        :: p(*)
  REAL, INTENT(in)                         :: a(*)
  REAL, INTENT(in)                         :: b(*)
  REAL, INTENT(out)                        :: c(*)
  REAL, INTENT(in)                         :: d(*)
  REAL, INTENT(out)                        :: delta(*)
  INTEGER, INTENT(in)                      :: nsoil
  INTEGER :: k
  INTEGER :: kk

  !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)))
  END DO

  !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)
  END DO
END SUBROUTINE rosr12

SUBROUTINE shflx ( s,stc,smc,smcmax,nsoil,t1,dt,yy,zz1,zsoil,tbot,  &
     b,f1,df1,ice )
  !    purpose:  to calculate soil heat flux.  soil temperature content
  !    =======   (an array containing soil temperatures) is also updated.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: s
  REAL, INTENT(in)                         :: stc( nsoil )
  REAL, INTENT(in out)                     :: smc( nsoil )
  REAL, INTENT(in out)                     :: smcmax
  INTEGER, INTENT(in out)                  :: nsoil
  REAL, INTENT(out)                        :: t1
  REAL, INTENT(in out)                     :: dt
  REAL, INTENT(in out)                     :: yy
  REAL, INTENT(in)                         :: zz1
  REAL, INTENT(in)                         :: zsoil( nsoil )
  REAL, INTENT(in out)                     :: tbot
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(in out)                     :: f1
  REAL, INTENT(in)                         :: df1
  INTEGER, INTENT(in out)                  :: ice
  INTEGER, PARAMETER :: nsold = 2

  REAL :: rhsts ( nsold )

  !hrt routine calcs the right hand side of the soil temp dif eqn
  IF(ice == 1) THEN
     CALL hrtice(rhsts,stc,smc,smcmax,nsoil,zsoil,yy,zz1,tbot,b, f1,df1)
  ELSE
     CALL hrt(rhsts,stc,smc,smcmax,nsoil,zsoil,yy,zz1,tbot,b,f1,df1)
  END IF

  !hstep routine calcs/updates soil temps based on rhsts
  CALL hstep ( stc,rhsts,dt,nsoil )

  !calc the grnd (skin) temperature
  t1 = (yy + (zz1 - 1.0) * stc(1)) / zz1

  !calc the sfc soil heat flux
  s = df1 * (stc(1) - t1) / (0.5 * zsoil(1))
END SUBROUTINE shflx

SUBROUTINE smflx ( eta1,smc,nsoil,cmc,etp1,dt,prcp1,zsoil,  &
     smcmax,b,pc,smcwlt,dksat,dwsat,smcref,shdfac,cmcmax,  &
     smcdry,cfactr, runoff1,runoff2, runoff3, edir1, ec1, ett1, sfctmp,q2,nroot)
  !    purpose:  to calculate 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.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: eta1
  REAL, INTENT(in)                         :: smc( nsoil )
  INTEGER, INTENT(in)                      :: nsoil
  REAL, INTENT(in)                         :: cmc
  REAL, INTENT(in)                         :: etp1
  REAL, INTENT(in)                         :: dt
  REAL, INTENT(in)                         :: prcp1
  REAL, INTENT(in)                         :: zsoil( nsoil )
  REAL, INTENT(in)                         :: smcmax
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(in out)                     :: pc
  REAL, INTENT(in out)                     :: smcwlt
  REAL, INTENT(in out)                     :: dksat
  REAL, INTENT(in out)                     :: dwsat
  REAL, INTENT(in out)                     :: smcref
  REAL, INTENT(in)                         :: shdfac
  REAL, INTENT(in)                         :: cmcmax
  REAL, INTENT(in out)                     :: smcdry
  REAL, INTENT(in)                         :: cfactr
  REAL, INTENT(in out)                     :: runoff1
  REAL, INTENT(in out)                     :: runoff2
  REAL, INTENT(in out)                     :: runoff3
  REAL, INTENT(out)                        :: edir1
  REAL, INTENT(out)                        :: ec1
  REAL, INTENT(out)                        :: ett1
  REAL, INTENT(in out)                     :: sfctmp
  REAL, INTENT(in out)                     :: q2
  INTEGER, INTENT(in out)                  :: nroot
  INTEGER, PARAMETER :: nsold = 2
  INTEGER :: k
  REAL :: dummy,devap,cmc2ms
  REAL :: beta
  REAL :: dew
  REAL :: drip
  REAL :: ec
  REAL :: edir
  REAL :: et     ( nsold )
  REAL :: ett
  REAL :: excess
  REAL :: flx1
  REAL :: flx2
  REAL :: flx3
  REAL :: pcpdrp
  REAL :: rho
  REAL :: rhsct
  REAL :: rhstt  ( nsold )
  REAL :: rib
  REAL :: runof
  REAL :: runoff
  REAL :: smca   ( nsold )
  REAL :: smcfg  ( nsold )
  REAL :: trhsct
  COMMON/rite/ beta,drip,ec,edir,ett,flx1,flx2,flx3,rho,runof, dew,rib

  !executable code begins here....if the potential evapotrans-
  !piration is greater than zero...
  dummy=0.
  edir = 0.
  ec = 0.
  ett = 0.
  DO  k = 1, nsoil
     et ( k ) = 0.
  END DO
  IF ( etp1 > 0.0 ) THEN
     !retrieve direct evaporation from soil surface
     edir = devap ( etp1, smc(1), zsoil(1), shdfac, smcmax,  &
          b, dksat, dwsat, smcdry,smcref, smcwlt)
     !initialize plant total transpiration, retrieve plant
     !transpiration, and accumulate it for all soil layers.
     ett = 0.
     IF(shdfac > 0.0) THEN
        CALL transp ( et, nsoil, etp1, smc, cmc, zsoil, shdfac,  &
             smcwlt, cmcmax, pc, cfactr, smcref, sfctmp,q2,nroot)

        DO  k = 1 , nsoil
           ett = ett + et ( k )
        END DO
     END IF
     !calculate canopy evaporation
     ec = shdfac * ( ( cmc / cmcmax ) ** cfactr ) * etp1
     cmc2ms=cmc/dt
  END IF

  !total up evap and transp types to obtain actual evapotransp
  edir1=edir
  ec1=ec
  ett1=ett
  eta1 = edir + ett + ec

  !compute the right hand side of the canopy eqn term ( rhsct )
  rhsct = shdfac * prcp1 - ec

  !convert rhsct (a rate) to trhsct (an amt) and add it to existing
  !cmc. if resulting amt exceeds max capacity, it becomes drip
  !and will fall to the grnd.
  drip = 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. - shdfac) * prcp1 + drip / dt

  !call subroutines srt and sstep to solve the soil moisture
  !tendency equations. iterate twice thru the calls. this eliminates
  !2-delta-t oscillations in the smc values even with timesteps of
  !up to 1 hour. the 1st set of calls yields 1st guess smc values,
  !the 2nd set of calls yields final smc values.
  IF ( pcpdrp > 0.0 ) THEN
     CALL srt ( rhstt,runoff,edir,et,smc,smc,nsoil,pcpdrp,zsoil,  &
          dwsat, dksat, smcmax, b, runoff1, runoff2, dt, smcwlt)
     CALL sstep ( smcfg,smc,dummy,rhstt,rhsct,dt,nsoil,smcmax,  &
          cmcmax, runoff3, zsoil)
     DO  k = 1, nsoil
        smca(k) = ( smc(k) + smcfg(k) ) * 0.5
     END DO
     CALL srt ( rhstt,runoff,edir,et,smc,smca,nsoil,pcpdrp,zsoil,  &
          dwsat, dksat, smcmax, b, runoff1, runoff2, dt, smcwlt)
     CALL sstep ( smc,smc,cmc,rhstt,rhsct,dt,nsoil,smcmax,  &
          cmcmax, runoff3, zsoil)
  ELSE
     CALL srt ( rhstt,runoff,edir,et,smc,smc,nsoil,pcpdrp,zsoil,  &
          dwsat, dksat, smcmax, b, runoff1, runoff2, dt, smcwlt)
     CALL sstep ( smc,smc,cmc,rhstt,rhsct,dt,nsoil,smcmax,  &
          cmcmax, runoff3, zsoil)
  END IF
  runof = runoff
END SUBROUTINE smflx

SUBROUTINE snopac ( etp,eta,prcp,prcp1,snowng,smc,smcmax,smcwlt,  &
     smcref, smcdry, cmc, cmcmax, nsoil, dt, q1,  &
     q2,t1,sfctmp,t24,th2,f,f1,s,stc,epsca,sfcprs, b, pc, rch, rr, cfactr, esd,  &
     zsoil, dwsat, dksat, tbot, shdfac, runoff1,  &
     runoff2,runoff3,edir1,ec1,ett1,nroot,snmax,ice)
  !    purpose:  to calculate 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.

  IMPLICIT NONE
  REAL, INTENT(in out)                     :: etp
  REAL, INTENT(out)                        :: eta
  REAL, INTENT(in)                         :: prcp
  REAL, INTENT(out)                        :: prcp1
  LOGICAL, INTENT(in)                      :: snowng
  REAL, INTENT(in out)                     :: smc( nsoil )
  REAL, INTENT(in out)                     :: smcmax
  REAL, INTENT(in out)                     :: smcwlt
  REAL, INTENT(in out)                     :: smcref
  REAL, INTENT(in out)                     :: smcdry
  REAL, INTENT(in out)                     :: cmc
  REAL, INTENT(in out)                     :: cmcmax
  INTEGER, INTENT(in out)                  :: nsoil
  REAL, INTENT(in)                         :: dt
  REAL, INTENT(in out)                     :: q1
  REAL, INTENT(in out)                     :: q2
  REAL, INTENT(out)                        :: t1
  REAL, INTENT(in)                         :: sfctmp
  REAL, INTENT(in out)                     :: t24
  REAL, INTENT(in)                         :: th2
  REAL, INTENT(in)                         :: f
  REAL, INTENT(in out)                     :: f1
  REAL, INTENT(out)                        :: s
  REAL, INTENT(in)                         :: stc( nsoil )
  REAL, INTENT(in out)                     :: epsca
  REAL, INTENT(in out)                     :: sfcprs
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(in out)                     :: pc
  REAL, INTENT(in)                         :: rch
  REAL, INTENT(in)                         :: rr
  REAL, INTENT(in out)                     :: cfactr
  REAL, INTENT(in out)                     :: esd
  REAL, INTENT(in)                         :: zsoil( nsoil )
  REAL, INTENT(in out)                     :: dwsat
  REAL, INTENT(in out)                     :: dksat
  REAL, INTENT(in out)                     :: tbot
  REAL, INTENT(in out)                     :: shdfac
  REAL, INTENT(in out)                     :: runoff1
  REAL, INTENT(in out)                     :: runoff2
  REAL, INTENT(in out)                     :: runoff3
  REAL, INTENT(in out)                     :: edir1
  REAL, INTENT(in out)                     :: ec1
  REAL, INTENT(in out)                     :: ett1
  INTEGER, INTENT(in out)                  :: nroot
  REAL, INTENT(out)                        :: snmax
  INTEGER, INTENT(in out)                  :: ice
  REAL :: qsat,etp3,seh,t14
  REAL :: beta
  REAL :: denom
  REAL :: dew
  REAL :: df1
  REAL :: drip
  REAL :: ec
  REAL :: edir
  REAL :: eta1
  REAL :: etp1
  REAL :: etp2
  REAL :: ett
  REAL :: ex
  REAL :: flx1
  REAL :: flx2
  REAL :: flx3
  REAL :: rho
  REAL :: rib
  REAL :: runoff
  REAL :: s1
  REAL :: t11
  REAL :: t12
  REAL :: t12a
  REAL :: t12b
  REAL :: yy
  REAL :: zz1
  COMMON/rite/ beta,drip,ec,edir,ett,flx1,flx2,flx3,rho,runoff, dew,rib
  REAL, PARAMETER :: cp=1004.5
  REAL, PARAMETER :: cph2o=4.218e+3
  REAL, PARAMETER :: cpice=2.106e+3
  REAL, PARAMETER :: lsubf=3.335e+5
  REAL, PARAMETER :: lsubc=2.5e+6
  REAL, PARAMETER :: lsubs=2.83e+6
  REAL, PARAMETER :: sigma=5.67e-8

  !executable code begins here...
  !convert etp fm  kg m-2 s-1  to  m s-1  and then to an amount
  !(m) and call it an effective snowpack reduction amt, etp2 (m).
  !this is the amt the snowpack would be reduced due to evaporation
  !from the snow sfc during the timestep.  evaporation will
  !proceed at the potntl rate unless the snow depth is less than
  !the expected snowpack reduction.
  prcp1=prcp1*0.001
  etp2 = etp * 0.001 * dt
  beta = 1.0
  IF(ice /= 1) THEN
     IF ( esd < etp2 ) THEN
        beta = esd / etp2
     END IF
  END IF
  dew = 0.0
  IF ( etp < 0.0 ) THEN
     dew = -etp * 0.001
  END IF

  !calc an 'effective snow-grnd sfc temp' based on heat fluxes
  !btwn the snow pack and the soil and on net radiation. the
  !constant 0.35 (w m-1 k-1) is the thermal diffusivity of snow,
  t12a = ((f - sigma * t24) / rch+th2-sfctmp-beta*epsca) / rr
  t12b = 0.35 * stc(1) / ( 10.0 * esd * rr * rch )
  denom = 1.0 + 0.35 / ( 10.0 * esd * rr * rch )
  t12 = (sfctmp + t12a + t12b ) / denom

  !recalculate the heat flux into or out of the soil,
  !and set the effective potential evapotransp to zero.
  !the constant 0.35 (w m-1 k-1) is the thermal diffusivity of snow.
  s = 0.35 * ( t12 - stc(1) ) / ( 10.0 * esd )
  s = MIN ( MAX ( s,-100.0 ), 100.0 )
  IF ( t12 <= 273.16 ) THEN
     !if the 'effective snow-grnd sfc temp' is at or blo freezing,
     !no snow melt will occur.  set the skin temp to this
     !effective temp and set the effective precip to zero.
     esd = amax1 ( 0.0 , esd - etp2 )
     t1 = t12
     flx3 = 0.0
     ex = 0.0
     snmax = 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, snmax.
     !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.
     !this is the old version of snow melt
     !  ex = 0.001 * rch * rr * ( t12 - 273.16 ) / lsubf
     !  snmax = ex * dt
     !  ex = amin1 ( esd / dt, ex )
     !  esd = amax1 ( esd - snmax, 0.0 )
     !  t1 = amax1 ( 273.16, (t12-(lsubf * ex * 1000.0 / (rch*rr))))
     !  flx3 = lsubf * ex * 1000.0
     !  prcp1 = prcp1 + ex
     !below is the correction of evaporation considering snow melting
     !                                              f. chen (96/3/18)
     t1 = 273.16
     qsat = (0.622*6.11e2)/(sfcprs-0.378*6.11e2)
     etp = rch*(qsat-q2)/cp
     etp2 = etp*0.001*dt
     beta = 1.0
     IF ( esd <= etp2 ) THEN
        beta = esd / etp2
        esd = 0.0
        snmax = 0.0
        ex = 0.0
     ELSE
        esd = MAX ( 0.0 , esd - etp2 )
        etp3 = etp*lsubc
        s = 0.35 * ( t1 - stc(1) ) / ( 10.0 * esd )
        s = MIN ( MAX ( s,-100.0 ), 100.0 )
        seh = rch*(t1-th2)
        t14 = t1*t1
        t14 = t14*t14
        flx3 = f - sigma*t14 - s - seh - etp3
        IF(flx3 <= 0.0) flx3=0.0
        ex = flx3*0.001/lsubf
        snmax = ex * dt
     END IF
     IF(snmax < esd) THEN
        esd = esd - snmax
     ELSE
        ex = esd/dt
        snmax = esd
        esd = 0.0
        flx3 = ex*1000.0*lsubf
        !call tdfcnd ( df1, smc(1), b, f1 )
        !yynum = f - sigma * t24
        !yy = sfctmp + (yynum/rch+th2-sfctmp-(flx3+etp3)) / rr
        !zz1 = df1 / ( -0.5 * zsoil(1) * rch * rr ) + 1.0
        !t1 = (yy + (zz1 - 1.0) * stc(1)) / zz1
        !s = df1 * ( stc(1) - t1 ) / ( 0.5 * zsoil(1) )
     END IF
     prcp1 = prcp1 + ex
  END IF

  ! set the effective potential evapotransp to zero
  etp1 = 0.0

  !smflx returns soil moisture values and preliminary values of
  !evapotranspiration.  in this, the snow pack case, the prelim-
  !values (eta1) are not used in subsequent calculation of e.
  IF(ice /= 1) THEN
     CALL smflx ( eta1,smc,nsoil,cmc,etp1,dt,prcp1,zsoil,  &
          smcmax,b,pc,smcwlt,dksat,dwsat,  &
          smcref,shdfac,cmcmax,smcdry,cfactr,runoff1,runoff2,  &
          runoff3, edir1, ec1, ett1,sfctmp,q2,nroot)
  END IF
  eta = beta * etp

  !calc/update the grnd sfc (skin) mixing ratio.
  !q1 = q2 + eta * cp / rch

  !retrieve the soil temp diffusivity/conductivity.
  CALL tdfcnd ( df1, smc(1), b, f1 )
  IF(ice == 1) df1=2.2

  !the 'adjusted top soil lyr temp' ( yy ) and the 'adjusted soil
  !heat flux' ( zz1 ) are set to the top soil lyr temp, and 1,
  !respectively.  these are close-enough approximations because
  !the sfc heat flux to be computed in shflx will effectively be
  !the flux at the snow top surface.  t11 is a dummy arguement
  !since we will not use its value as revised by shflx.
  zz1 = 1.0
  yy = stc(1)-0.5*s*zsoil(1)*zz1/df1
  t11 = t1

  !shflx will calc/update the soil temps. note the soil heat flux
  !( s1 ) and the skin temp ( t11 ) output from this call are not
  !used in any subsequent calculations.
  CALL shflx ( s1,stc,smc,smcmax,nsoil,t11,dt,yy,zz1,zsoil,tbot, b,f1,df1,ice )

  !if precip is falling, calc heat flux from snow sfc to newly
  !accumulating precip.
  flx1 = 0.0
  IF ( snowng ) THEN
     flx1 = cpice * prcp * ( t1 - sfctmp )
  ELSE
     IF ( prcp > 0.0 ) flx1 = cph2o * prcp * ( t1 - sfctmp )
  END IF
END SUBROUTINE snopac


SUBROUTINE srt (rhstt,runoff,edir,et,smc,smca,nsoil,pcpdrp,  &
     zsoil,dwsat,dksat,smcmax,b, runoff1, runoff2,dt, smcwlt)
  !    purpose:  to calculate 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.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: rhstt( nsoil )
  REAL, INTENT(in out)                     :: runoff
  REAL, INTENT(in)                         :: edir
  REAL, INTENT(in)                         :: et( nsoil )
  REAL, INTENT(in)                         :: smc( nsoil )
  REAL, INTENT(in)                         :: smca( nsoil )
  INTEGER, INTENT(in)                     :: nsoil
  REAL, INTENT(in)                         :: pcpdrp
  REAL, INTENT(in)                         :: zsoil( nsoil )
  REAL, INTENT(in out)                     :: dwsat
  REAL, INTENT(in)                         :: dksat
  REAL, INTENT(in)                         :: smcmax
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(out)                        :: runoff1
  REAL, INTENT(out)                        :: runoff2
  REAL, INTENT(in)                         :: dt
  REAL, INTENT(in)                         :: smcwlt
  INTEGER, PARAMETER :: nsold = 2
  INTEGER :: k,iohinf,kdt
  REAL :: refkdt,refdk,dt1,smcav,d1max,d2max
  REAL :: d1,d2,dd,val,ddt,px,sstt
  REAL :: ai     ( nsold )
  REAL :: bi     ( nsold )
  REAL :: ci     ( nsold )
  REAL :: ddz
  REAL :: ddz2
  REAL :: denom
  REAL :: denom2
  REAL :: dsmdz
  REAL :: dsmdz2
  REAL :: infmax
  REAL :: mxsmc
  REAL :: mxsmc2
  REAL :: numer
  REAL :: pddum
  REAL :: wcnd
  REAL :: wcnd2
  REAL :: wdf
  REAL :: wdf2
  COMMON /abci/ ai, bi, ci

  !determine rainfall infiltration rate and runoff
  ! *************  include the infiltration formule from schaake and koren
  !           model
  !    modified by q duan
  !    replace the maximum infiltration equation in osu model
  !    by schaake/koren expression

  !    initialize parameters and variables
  iohinf=1

  !determine rainfall infiltration rate and runoff
  pddum = pcpdrp
  runoff1 = 0.0
  IF ( pcpdrp /= 0.0 ) THEN
     !c++  modified by q. duan, 5/16/94
     !        if (iohinf .eq. 1) then
     refkdt = 3.0
     refdk = 3.4341e-6
     dt1 = dt/86400.
     smcav = smcmax - smcwlt
     d1max = -zsoil(1)*smcav
     d2max = (zsoil(1)-zsoil(2))*smcav
     d1 = d1max * (1.0 - (smca(1)-smcwlt)/smcav)
     d2 = d2max * (1.0 - (smca(2)-smcwlt)/smcav)
     dd = d1 + d2
     kdt = refkdt * dksat/refdk
     !val = (1.-exp(-kdt*sqrt(dt1)))
     !remove the sqrt by f. chen, 07/08/96
     val = (1.-EXP(-kdt*dt1))
     ddt = dd*val
     px = pcpdrp*dt
     IF(px < 0.0) px = 0.0
     infmax = (px*(ddt/(px+ddt)))/dt
     infmax = MIN(infmax, dksat)
     !else
     !infmax = -dwsat*(smcmax-smca(1))/(0.5*zsoil(1))+dksat
     !if (smca(1).ge.smcmax) infmax = dksat
     !end if

     !runoff1: surface runoff
     IF ( pcpdrp > infmax ) THEN
        runoff1 = pcpdrp - infmax
        pddum = infmax
     END IF

  END IF


  !below is the original osu model
  !      pddum = pcpdrp
  !      print*,' ************** in srt *******************'
  !      print*,' pcpdrp=', pcpdrp
  !      runoff1 = 0.0
  !      if ( pcpdrp .ne. 0.0 ) then
  !        infmax = - dwsat * (smcmax - smca(1))/(0.5 * zsoil(1)) + dksat
  !        if ( smca(1) .ge. smcmax ) infmax = dksat
  !        if ( pcpdrp .gt. infmax ) then
  !          runoff1 = pcpdrp - infmax
  ! runoff1: surface runoff *******
  !          pddum = infmax
  !        end if
  !      end if

  !determine soil moisture diffusivity (wdf) & conductivity (wcnd)
  !coefs for top soil lyr using the higher of smc for the top
  !soil layer or smc for the next lower soil layer.
  mxsmc = MAX( smca(1), smca(2) )
  CALL wdfcnd ( wdf,wcnd,mxsmc,smcmax,b,dksat,dwsat )

  !calc the matrix coefficients ai, bi, and ci for the top layer
  ddz = 1. / ( -.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 = ( smc(1) - smc(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
        !retrieve the soil water diffusivity for the wetter of the
        !current or the next lower soil layer
        !mxsmc2 = max ( smca(k), smca(k+1) )

        !pilps test
        mxsmc2 =  smca(k)
        CALL wdfcnd ( wdf2,wcnd2,mxsmc2,smcmax,b,dksat,dwsat )

        !calc some partial products for later use in calc'ng rhstt
        denom = ( zsoil(k-1) - zsoil(k+1) )
        dsmdz2 = ( smc(k) - smc(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

        !retrieve the soil water diffusivity and hydraulic
        !conductivity for this layer
        CALL wdfcnd( wdf2,wcnd2,smca(nsoil),smcmax,b,dksat,dwsat )

        !calc a partial product for later use in calc'ng rhstt
        dsmdz2 = 0.0

        !set matrix coef ci to zero
        ci(k) = 0.0
     END IF
     !calc rhstt for this layer after calc'ng its numerator
     numer = (wdf2 * dsmdz2) + 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
     IF(k == nsoil) THEN
        !runoff2: ground water runoff 
        runoff2 = wcnd2
     END IF
     IF ( k /= nsoil ) THEN
        wdf = wdf2
        wcnd = wcnd2
        dsmdz = dsmdz2
        ddz = ddz2
     END IF
  END DO
END SUBROUTINE srt

SUBROUTINE sstep ( smcout, smcin, cmc, rhstt, rhsct, dt,  &
     nsoil, smcmax, cmcmax, runoff3, zsoil )
  !    purpose:  to calculate/update the soil moisture content values
  !    =======   and the canopy moisture content values.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: smcout( nsoil )
  REAL, INTENT(in)                         :: smcin( nsoil )
  REAL, INTENT(out)                        :: cmc
  REAL, INTENT(out)                        :: rhstt( nsoil )
  REAL, INTENT(in)                         :: rhsct
  REAL, INTENT(in)                         :: dt
  INTEGER, INTENT(in)                      :: nsoil
  REAL, INTENT(in)                         :: smcmax
  REAL, INTENT(in)                         :: cmcmax
  REAL, INTENT(out)                        :: runoff3
  REAL, INTENT(in)                         :: zsoil(nsoil)
  INTEGER, PARAMETER :: nsold = 2
  INTEGER :: k
  REAL :: ddz,kk11
  REAL :: ai     ( nsold )
  REAL :: bi     ( nsold )
  REAL :: ci     ( nsold )
  COMMON /abci/ ai, bi, ci

  !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
  END DO

  !call rosr12 to solve the tri-diagonal matrix
  CALL rosr12 ( ci, ai, bi, ci, rhstt, rhstt, nsoil )

  !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
  runoff3 = 0.
  DO  k = 1 , nsoil
     smcout(k) = smcin(k) + ci(k)
     IF(smcout(k) > smcmax) THEN
        IF(k == 1) THEN
           ddz = -zsoil(1)
        ELSE
           kk11 = k -1
           ddz = -zsoil(k) + zsoil(kk11)
        END IF
        runoff3 = runoff3 - ((smcmax-smcout(k))*ddz)
     END IF
     smcout(k) = MAX ( MIN ( smcout(k), smcmax ), 0.02 )
  END DO

  !convert rhsct to an 'amount' value and add to previous
  !cmc value to get new cmc value.
  cmc = cmc + dt * rhsct
  cmc = MIN ( MAX ( cmc, 0. ), cmcmax )
END SUBROUTINE sstep

SUBROUTINE tdfcnd ( df, smc, b, f1 )
  !    purpose:  to calculate thermal diffusivity and conductivity of
  !    =======   the soil for a given point and time.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: df
  REAL, INTENT(in)                         :: smc
  REAL, INTENT(in)                         :: b
  REAL, INTENT(in)                         :: f1

  REAL :: pf

  !if the soil has any moisture content compute a partial sum/product
  !otherwise use a constant value which works well with most soils
  IF ( smc > 0. ) THEN
     pf = f1 - b * alog10 ( smc )
  ELSE
     pf = 5.2
  END IF

  !if the partial sum/product is lower than or equal to 5.1, compute
  !a thermal conductivity/diffusivity.  cap its value at 1.9.
  !otherwise set the thermal conductivity/diffusivity to a constant.
  IF ( pf <= 5.1 ) THEN
     df = EXP ( -( 2.7 + pf ) ) * 420.0
     df = amin1( df, 1.9 )
  ELSE
     df = 0.1744
  END IF

  !snesitivity test
  !df = 1.2*df
END SUBROUTINE tdfcnd

SUBROUTINE transp ( et, nsoil, etp1, smc, cmc, zsoil,  &
     shdfac,smcwlt,cmcmax,pc,cfactr,smcref,sfctmp,q2,nroot)
  !    purpose:  to calculate transpiration from the vegtyp for this pt.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: et( nsoil )
  INTEGER, INTENT(in)                      :: nsoil
  REAL, INTENT(in)                         :: etp1
  REAL, INTENT(in)                         :: smc( nsoil )
  REAL, INTENT(in out)                     :: cmc
  REAL, INTENT(in)                         :: zsoil( nsoil )
  REAL, INTENT(in)                         :: shdfac
  REAL, INTENT(in)                         :: smcwlt
  REAL, INTENT(in out)                     :: cmcmax
  REAL, INTENT(in)                         :: pc
  REAL, INTENT(in out)                     :: cfactr
  REAL, INTENT(in)                         :: smcref
  REAL, INTENT(in out)                     :: sfctmp
  REAL, INTENT(in out)                     :: q2
  INTEGER, INTENT(in)                      :: nroot
  INTEGER :: k

  REAL :: etp1a
  REAL :: gx

  !initialize  plant transp to zero for all soil layers.
  DO  k = 1, nsoil
     et(k) = 0.
  END DO

  !calc an 'adjusted' potntl transpiration
  etp1a = shdfac * pc * etp1 * (1.0 - (cmc /cmcmax) ** cfactr)

  !calc the plant transp.  first, calc gx and et for the 1st
  !(top) soil layer.  gx is the simple transpiration function, a
  !ratio of the actual smc surplus to the reference smc surplus.
  !use the absolute depth of the 1st layer (zsoil(1)) to get
  !the proportion of the potential in the final calculation.
  gx = ( smc(1) - smcwlt ) / ( smcref - smcwlt )
  gx = MAX ( MIN ( gx, 1. ), 0. )
  !test canopy resistance
  gx = 1.0
  et(1) = ( zsoil(1) / zsoil(nroot) ) * gx * etp1a

  !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. )
     !test canopy resistance
     gx = 1.0
     et(k) = ((zsoil(k)-zsoil(k-1))/zsoil(nroot))*gx*etp1a
  END DO
END SUBROUTINE transp

SUBROUTINE wdfcnd ( wdf,wcnd,smc,smcmax,b,dksat,dwsat )
  !    purpose:  to calculate soil water diffusivity and soil
  !    =======   hydraulic conductivity.

  IMPLICIT NONE
  REAL, INTENT(out)                        :: wdf
  REAL, INTENT(out)                        :: wcnd
  REAL, INTENT(out)                        :: smc
  REAL, INTENT(out)                        :: smcmax
  REAL, INTENT(in)                         :: b
  REAL, INTENT(in)                         :: dksat
  REAL, INTENT(in)                         :: dwsat

  REAL :: expon
  REAL :: factr2

  !calc the ratio of the actual to the max psbl soil h2o content
  smc = smc
  smcmax = smcmax
  factr2 = smc / smcmax

  !prep an expntl coef and calc the soil water diffusivity
  expon = b + 2.0
  wdf = dwsat * factr2 ** expon

  !reset the expntl coef and calc the hydraulic conductivity
  expon = ( 2.0 * b ) + 3.0
  wcnd = dksat * factr2 ** expon
END SUBROUTINE wdfcnd

SUBROUTINE redprm(vegtyp,soityp,  &
     cmcmax,topt,rsmax,albedo,z0,shdfac,nroot,rcmin,rgl,hs,  &
     b,dksat,dwsat,smcmax,smcwlt,smcref,smcdry,f1)
  !this subroutine reads the soil and vegetation parameters
  !required for the execusion of the osu land-surface scheme.
  !                                        by f. chen 3/15/95

  IMPLICIT NONE
  INTEGER, INTENT(in out)                  :: vegtyp
  INTEGER, INTENT(in out)                  :: soityp
  REAL, INTENT(in out)                     :: cmcmax
  REAL, INTENT(in out)                     :: topt
  REAL, INTENT(in out)                     :: rsmax
  REAL, INTENT(in out)                     :: albedo
  REAL, INTENT(in out)                     :: z0
  REAL, INTENT(in out)                     :: shdfac
  INTEGER, INTENT(in out)                  :: nroot
  REAL, INTENT(in out)                     :: rcmin
  REAL, INTENT(in out)                     :: rgl
  REAL, INTENT(in out)                     :: hs
  REAL, INTENT(in out)                     :: b
  REAL, INTENT(in out)                     :: dksat
  REAL, INTENT(in out)                     :: dwsat
  REAL, INTENT(in out)                     :: smcmax
  REAL, INTENT(in out)                     :: smcwlt
  REAL, INTENT(in out)                     :: smcref
  REAL, INTENT(in out)                     :: smcdry
  REAL, INTENT(in out)                     :: f1
  REAL :: psisat
  !set-up soil parameters
  CALL prmsoi(soityp,b,smcdry,f1,smcmax,smcref, psisat,dksat,dwsat,smcwlt)

  !set-up vegetation parameters
  CALL  prmveg(vegtyp,cmcmax,topt,rsmax,albedo,z0, shdfac,nroot,rcmin,rgl,hs)
END SUBROUTINE redprm


SUBROUTINE  prmveg(vegtyp,cmcmax,topt,rsmax,albedo,z0,  &
     shdfac,nroot,rcmin,rgl,hs)
  !set-up vegetation parameters for a given vegetaion type
  !  input:
  !        vegtyp: vegetation type
  ! ouput:
  !       vegetation parameters:
  !   albedo: sfc albedo
  !  cmxtbl: max cnpy capacity
  !      z0: roughness length
  !  shdfac: plant shade factor
  !                nroot: rooting depth
  !                rcmin: mimimum stomatal resistance
  !        rsmax and rgl: parameters used in radiation stress function
  !                   hs: parameter used in vapor pressure deficit function
  !                 topt: parameter used in temperature stress function
  !               cmcmax: maximum canopy water capacity
  ! *************************************************************************
  !     ssib vegetation types (dorman and sellers, 1989; jam)
  !    1:   broadleaf-evergreen trees  (tropical forest)
  !    2:   broadleaf-deciduous tress
  !    3:   broadleaf and needleleaf tress (mixed forest)
  !    4:   needleleaf-evergreen trees
  !    5:   needleleaf-deciduous tress (larch)
  !    6:   broadleaf tress 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 for the type 7)
  !   13:   glacial
  !***************************************************************************

  IMPLICIT NONE
  INTEGER, INTENT(in out)                  :: vegtyp
  REAL, INTENT(in out)                     :: cmcmax
  REAL, INTENT(in out)                     :: topt
  REAL, INTENT(in out)                     :: rsmax
  REAL, INTENT(in out)                     :: albedo
  REAL, INTENT(in out)                     :: z0
  REAL, INTENT(out)                        :: shdfac
  INTEGER, INTENT(in out)                  :: nroot
  REAL, INTENT(out)                        :: rcmin
  REAL, INTENT(out)                        :: rgl
  REAL, INTENT(out)                        :: hs
  INTEGER :: nrotbl(13)
  REAL :: albtbl(13),z0tbl(13),shdtbl(13),rsmtbl(13), rgltbl(13),hstbl(13)

  DATA albtbl/0.11, 0.19,  0.16,  0.13,  0.19,  0.19,  0.19,  &
       0.29, 0.29,  0.14,  0.15,  0.19,  0.15/
  DATA  z0tbl/2.653,0.826, 0.563, 1.089, 0.854, 0.856, 0.075,  &
       0.238,0.065, 0.076, 0.011, 0.075, 0.011/
  DATA shdtbl/0.90, 0.80,  0.80,  0.85,  0.40,  0.50,  0.80,  &
       0.20, 0.15,  0.30,  0.0,   0.80,  0.0/
  DATA nrotbl/13*2/
  DATA rsmtbl/150.0,100.0, 125.0, 150.0, 100.0, 70.0,  40.0,  &
       300.0,400.0, 150.0, 999.0, 40.0,  999.0/
  DATA rgltbl/30.0, 30.0,  30.0,  30.0,  30.0,  65.0,  100.0,  &
       100.0,100.0, 100.0, 999.0, 100.0, 999.0/
  DATA hstbl/41.69, 54.53, 51.93, 47.35, 47.35, 54.53, 36.35,  &
       42.00, 42.00, 42.00, 999.0, 36.35, 999.0/
  SAVE

  ! albedo = albtbl(vegtyp)
  !        z0 = z0tbl(vegtyp)
  !        nroot = nrotbl(vegtyp)
  !        shdfac = shdtbl(vegtyp)
  rcmin = rsmtbl(vegtyp)
  rgl = rgltbl(vegtyp)
  hs = hstbl(vegtyp)
  IF(vegtyp == 11) shdfac=0.0
END SUBROUTINE  prmveg


SUBROUTINE  prmsoi(soltyp,sbb,sdrysmc,sf11,smaxsmc,srefsmc,  &
     ssatpsi,ssatdk,ssatdw,swltsmc)
  !set-up soil parameters for given soil type
  !  input:
  !        soltyp: soil type
  ! ouput:
  !       soil parameters:
  !  smaxsmc: max soil moisture content
  !  srefsmc: reference soil moisture
  !  swltsmc: wilting pt soil moisture contents
  !  sdrysmc: air dry soil moist content limits
  !  ssatpsi: sat soil potential coefs.
  !   ssatdk:  sat soil diffusivity/conductivity coefs.
  !   sbb:     soil diffusivity/conductivity coef.
  !  ssatdw:  sat soil diffusivity/conductivity coefs.
  !  sf11:    soil diffusivity/conductivity coef.
  ! *************************************************************************
  !  soil types     zobler (1986)         cosby et al (1984)
  !             1          coarse                loamy sand
  !             2          medium                silty clay loam
  !             3          fine                  light clay
  !             4          coarse-medium         sandy loam
  !             5          coarse-fine           sandy clay
  !             6          medium-fine           clay loam
  !             7          coarse-med-fine       sandy clay loam
  !             8          organic               loam
  !             9          land ice              loamy sand
  !***************************************************************************

  IMPLICIT NONE
  INTEGER, INTENT(in out)                  :: soltyp
  REAL, INTENT(out)                        :: sbb
  REAL, INTENT(out)                        :: sdrysmc
  REAL, INTENT(out)                        :: sf11
  REAL, INTENT(out)                        :: smaxsmc
  REAL, INTENT(out)                        :: srefsmc
  REAL, INTENT(out)                        :: ssatpsi
  REAL, INTENT(out)                        :: ssatdk
  REAL, INTENT(out)                        :: ssatdw
  REAL, INTENT(out)                        :: swltsmc

  REAL :: bb(9),drysmc(9),f11(9),maxsmc(9),refsmc(9),satpsi(9),  &
       satdk(9),satdw(9),wltsmc(9)

  DATA maxsmc/0.421, 0.464, 0.468, 0.434, 0.406, 0.465, 0.404, 0.439, 0.421/
  DATA drysmc/0.07,  0.14,  0.22,  0.08,  0.18,  0.16,  0.120, 0.10,  0.07/
  DATA satpsi/0.04,  0.62, 0.47, 0.14, 0.10, 0.26, 0.14, 0.36,  0.04/
  DATA satdk /1.41e-5, 0.20e-5, 0.10e-5, 0.52e-5, 0.72e-5,  &
       0.25e-5, 0.45e-5, 0.34e-5, 1.41e-5/
  DATA bb    /4.26,  8.72, 11.55, 4.74, 10.73, 8.17, 6.77, 5.25,  4.26/
  DATA refsmc/0.283, 0.387, 0.412, 0.312, 0.338, 0.382, 0.315, 0.329, 0.283/
  DATA wltsmc/0.029, 0.119, 0.139, 0.047, 0.010, 0.103, 0.069, 0.066, 0.029/
  DATA satdw  /5.71e-6, 2.33e-5, 1.16e-5, 7.95e-6, 1.90e-5,  &
       1.14e-5, 1.06e-5, 1.46e-5, 5.71e-6/
  DATA f11    /-0.999, -1.116, -2.137, -0.572, -3.201, -1.302,  &
       -1.519, -0.329, -0.999/

  ! note: satdw = bb*satdk*(satpsi/maxsmc)
  !         f11 = alog10(satpsi) + bb*alog10(maxsmc) + 2.0
  !       refsmc1=maxsmc*(5.79e-9/satdk)**(1/(2*bb+3)) 5.79e-9 m/s= 0.5 mm/day
  !       refsmc=refsmc1+1./3.(maxsmc-refsmc1)
  !       wltsmc1=maxsmc*(200./satpsi)**(-1./bb)    (wetzel and chang, 1987)
  !       wltsmc=wltsmc1-0.5*wltsmc1
  SAVE

  !soil parameters
  sbb = bb(soltyp)
  sdrysmc = drysmc(soltyp)
  sf11 = f11(soltyp)
  smaxsmc = maxsmc(soltyp)
  srefsmc = refsmc(soltyp)
  ssatpsi = satpsi(soltyp)
  ssatdk = satdk(soltyp)
  ssatdw = satdw(soltyp)
  swltsmc = wltsmc(soltyp)
END SUBROUTINE  prmsoi

SUBROUTINE dcoef ( z, z0, t1v, th2v, sfcspd, cm, ch )
  !    name:  determine coefficients (dcoef)       version: n/a

  IMPLICIT NONE
  REAL, INTENT(in)   :: z
  REAL, INTENT(in)   :: z0
  REAL, INTENT(in)   :: t1v
  REAL, INTENT(in)   :: th2v
  REAL, INTENT(out)  :: sfcspd
  REAL, INTENT(out)  :: cm
  REAL, INTENT(out)  :: ch
  REAL :: a,ah,am,beta,b1,b2,ctp,cup,cus,dew,drip,ec
  REAL :: edir,ett,exmch,flx1,flx2,flx3,g,pr,rho,rib
  REAL :: runoff,vk,z0h
  COMMON/rite/ beta,drip,ec,edir,ett,flx1,flx2,flx3,rho,runoff, dew,rib

  DATA b1     / 9.4 /
  DATA b2     / 15.0 /
  DATA cus    / 7.4 /
  DATA exmch  / -1. /
  DATA g      / 9.806 /

  !data pr     / .74 /
  !set pr=1 the same as in ek's version for pilps 1994
  DATA pr     / 1.0 /
  DATA vk     / .4 /

  sfcspd=MAX(sfcspd, 0.01)

  !set z0h as function of z0 as in beljaars and betts (1992?)
  z0h = z0/10.

  !calc a friction velocity for use in calculating the drag
  !coef for momentum and one for use in calculating the drag
  !coef for heat.
  a = vk / alog( z / z0 )
  am = a * a
  ah = a * vk / alog( z / z0h )

  !calc a bulk richardson number.  constrain its value in the
  !stable case to a maximum of 1.0 to avoid creating exchange
  !coefficients that approach zero.
  rib = g * z * ( th2v - t1v ) / ( th2v * sfcspd * sfcspd )
  rib = MIN( rib, 1.0 )

  IF ( rib >= 0. ) THEN
     !if the richardson number is .ge. zero, the air is stably
     !stratified.  calc the drag coefficients using a method
     !developed by mahrt (monthly weather review, 1987).  the
     !smallest allowable cm value is 1.0e-6.
     cm = am * sfcspd * EXP( exmch * rib )
     IF ( cm < 1.0e-6) THEN
        cm = 1.0e-6
        ch = 1.0e-6
     ELSE
        ch = ( cm * ah / am ) / pr
     END IF
  ELSE

     !otherwise, the air is unstably stratified and the drag
     !coefficients will be calculated as follows.
     cup = cus * am * b1 * SQRT( -rib * z / z0 )
     ctp = cus * ah * b1 * SQRT( -rib * z / z0h )
     cm = ( 1.0 - (b1 * rib)/(1.0 + cup) ) * am * sfcspd
     ch = ( 1.0 - (b2 * rib)/(1.0 + ctp) ) * ah * sfcspd / pr
  END IF
END SUBROUTINE dcoef
