/wrfv2_fire/phys/module_sf_noahlsm.F

Large files files are truncated, but you can click here to view the full file

MODULE module_sf_noahlsm
  USE module_model_constants, only : CP, R_D, XLF, XLV, RHOWATER, STBOLT, KARMAN

!   REAL, PARAMETER    :: CP = 1004.5
  REAL, PARAMETER      :: RD = 287.04, SIGMA = 5.67E-8,                 &
                          CPH2O = 4.218E+3,CPICE = 2.106E+3,            &
                          LSUBF = 3.335E+5,                             &
                          EMISSI_S = 0.95

! VEGETATION PARAMETERS
        INTEGER :: LUCATS , BARE
        INTEGER :: NATURAL
        integer, PARAMETER :: NLUS=50
        CHARACTER(LEN=256) LUTYPE
        INTEGER, DIMENSION(1:NLUS) :: NROTBL
        real, dimension(1:NLUS) ::  SNUPTBL, RSTBL, RGLTBL, HSTBL,                &
                                    SHDTBL, MAXALB,                               &
                                    EMISSMINTBL, EMISSMAXTBL,                     &
                                    LAIMINTBL, LAIMAXTBL,                         &
                                    Z0MINTBL, Z0MAXTBL,                           &
                                    ALBEDOMINTBL, ALBEDOMAXTBL
        REAL ::   TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA

! SOIL PARAMETERS
        INTEGER :: SLCATS
        INTEGER, PARAMETER :: NSLTYPE=30
        CHARACTER(LEN=256) SLTYPE
        REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,F11,                           &
        MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ

! LSM GENERAL PARAMETERS
        INTEGER :: SLPCATS
        INTEGER, PARAMETER :: NSLOPE=30
        REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
        REAL ::  SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA,           &
                 REFKDT_DATA,FRZK_DATA,ZBOT_DATA,  SMLOW_DATA,SMHIGH_DATA,        &
                        CZIL_DATA
        REAL ::  LVCOEF_DATA

        CHARACTER*256  :: err_message

        integer, private :: iloc, jloc
!
CONTAINS
!

      SUBROUTINE SFLX (IILOC,JJLOC,ICE,FFROZP,ISURBAN,DT,ZLVL,NSOIL,SLDPTH,         &    !C
                       LOCAL,                                           &    !L
                       LLANDUSE, LSOIL,                                 &    !CL
                       LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2,SFCSPD,  &    !F
                       COSZ,PRCPRAIN, SOLARDIRECT,                      &    !F
                       TH2,Q2SAT,DQSDT2,                                &    !I
                       VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHDMIN,SHDMAX,    &    !I
                       ALB, SNOALB,TBOT, Z0BRD, Z0, EMISSI, EMBRD,      &    !S
                       CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM,    &    !H
! ----------------------------------------------------------------------
! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
! MODEL).  OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
! ----------------------------------------------------------------------
                       ETA,SHEAT, ETA_KINEMATIC,FDOWN,                  &    !O
                       EC,EDIR,ET,ETT,ESNOW,DRIP,DEW,                   &    !O
                       BETA,ETP,SSOIL,                                  &    !O
                       FLX1,FLX2,FLX3,                                  &    !O
                       SNOMLT,SNCOVR,                                   &    !O
                       RUNOFF1,RUNOFF2,RUNOFF3,                         &    !O
                       RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL,             &    !O
                       SOILW,SOILM,Q1,SMAV,                             &    !D
                       RDLAI2D,USEMONALB,                               &
                       SNOTIME1,                                        &
                       RIBB,                                            &
                       SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT)                    !P
! ----------------------------------------------------------------------
! SUBROUTINE SFLX - UNIFIED NOAHLSM VERSION 1.0 JULY 2007
! ----------------------------------------------------------------------
! SUB-DRIVER FOR "Noah 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)
! ----------------------------------------------------------------------
! SFLX ARGUMENT LIST KEY:
! ----------------------------------------------------------------------
!  C  CONFIGURATION INFORMATION
!  L  LOGICAL
! CL  4-string character bearing logical meaning
!  F  FORCING DATA
!  I  OTHER (INPUT) FORCING DATA
!  S  SURFACE CHARACTERISTICS
!  H  HISTORY (STATE) VARIABLES
!  O  OUTPUT VARIABLES
!  D  DIAGNOSTIC OUTPUT
!  P  Parameters
!  Msic Miscellaneous terms passed from gridded driver
! ----------------------------------------------------------------------
! 1. CONFIGURATION INFORMATION (C):
! ----------------------------------------------------------------------
!   ICE        SEA-ICE FLAG  (=1: SEA-ICE, =0: LAND (NO ICE), =-1 LAND-ICE).
!   DT         TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
!                1800 SECS OR LESS)
!   ZLVL       HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
!   NSOIL      NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
!                PARAMETER NSOLD SET BELOW)
!   SLDPTH     THE THICKNESS OF EACH SOIL LAYER (M)
! ----------------------------------------------------------------------
! 2. LOGICAL:
! ----------------------------------------------------------------------
!   LCH       Exchange coefficient (Ch) calculation flag (false: using
!                ch-routine SFCDIF; true: Ch is brought in)
!   LOCAL      Flag for local-site simulation (where there is no
!              maps for albedo, veg fraction, and roughness
!             true:  all LSM parameters (inluding albedo, veg fraction and
!                    roughness length) will be defined by three tables
!   LLANDUSE  (=USGS, using USGS landuse classification)
!   LSOIL     (=STAS, using FAO/STATSGO soil texture classification)
! ----------------------------------------------------------------------
! 3. FORCING DATA (F):
! ----------------------------------------------------------------------
!   LWDN       LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
!   SOLDN      SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
!   SOLNET     NET DOWNWARD SOLAR RADIATION ((W M-2; POSITIVE)
!   SFCPRS     PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
!   PRCP       PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
!   SFCTMP     AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
!   TH2        AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
!   Q2         MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
!   COSZ       Solar zenith angle (not used for now)
!   PRCPRAIN   Liquid-precipitation rate (KG M-2 S-1) (not used)
! SOLARDIRECT  Direct component of downward solar radiation (W M-2) (not used)
!   FFROZP     FRACTION OF FROZEN PRECIPITATION
! ----------------------------------------------------------------------
! 4. OTHER FORCING (INPUT) DATA (I):
! ----------------------------------------------------------------------
!   SFCSPD     WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
!   Q2SAT      SAT SPECIFIC HUMIDITY AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
!   DQSDT2     SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
!                (KG KG-1 K-1)
! ----------------------------------------------------------------------
! 5. CANOPY/SOIL CHARACTERISTICS (S):
! ----------------------------------------------------------------------
!   VEGTYP     VEGETATION TYPE (INTEGER INDEX)
!   SOILTYP    SOIL TYPE (INTEGER INDEX)
!   SLOPETYP   CLASS OF SFC SLOPE (INTEGER INDEX)
!   SHDFAC     AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
!                (FRACTION= 0.0-1.0)
!   SHDMIN     MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
!                (FRACTION= 0.0-1.0) <= SHDFAC
!   PTU        PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
!                (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
!                VEG PARMS)
!   ALB        BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
!                DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
!                MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
!                INCLUDE DIURNAL SUN ANGLE EFFECT)
!   SNOALB     UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
!                ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
!   TBOT       BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
!                TEMPERATURE)
!   Z0BRD      Background fixed roughness length (M)
!   Z0         Time varying roughness length (M) as function of snow depth
!
!   EMBRD      Background surface emissivity (between 0 and 1)
!   EMISSI     Surface emissivity (between 0 and 1)
! ----------------------------------------------------------------------
! 6. HISTORY (STATE) VARIABLES (H):
! ----------------------------------------------------------------------
!  CMC         CANOPY MOISTURE CONTENT (M)
!  T1          GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
!  STC(NSOIL)  SOIL TEMP (K)
!  SMC(NSOIL)  TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
!  SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
!                NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
!  SNOWH       ACTUAL SNOW DEPTH (M)
!  SNEQV       LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
!                NOTE: SNOW DENSITY = SNEQV/SNOWH
!  ALBEDO      SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
!                =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
!                =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
!  CH          SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
!                (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
!                IT HAS BEEN MULTIPLIED BY WIND SPEED.
!  CM          SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
!                CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
!                MULTIPLIED BY WIND SPEED.
! ----------------------------------------------------------------------
! 7. OUTPUT (O):
! ----------------------------------------------------------------------
! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL.  FOR THIS APPLICATION,
! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
! NECESSARY.  OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
!   ETA        ACTUAL LATENT HEAT FLUX (W m-2: NEGATIVE, IF UP FROM
!              SURFACE)
!  ETA_KINEMATIC atctual latent heat flux in Kg m-2 s-1
!   SHEAT      SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
!              SURFACE)
!   FDOWN      Radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
! ----------------------------------------------------------------------
!   EC         CANOPY WATER EVAPORATION (W m-2)
!   EDIR       DIRECT SOIL EVAPORATION (W m-2)
!   ET(NSOIL)  PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
!                 (W m-2)
!   ETT        TOTAL PLANT TRANSPIRATION (W m-2)
!   ESNOW      SUBLIMATION FROM (OR DEPOSITION TO IF <0) SNOWPACK
!                (W m-2)
!   DRIP       THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
!                WATER-HOLDING CAPACITY (M)
!   DEW        DEWFALL (OR FROSTFALL FOR T<273.15) (M)
! ----------------------------------------------------------------------
!   BETA       RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
!   ETP        POTENTIAL EVAPORATION (W m-2)
!   SSOIL      SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
! ----------------------------------------------------------------------
!   FLX1       PRECIP-SNOW SFC (W M-2)
!   FLX2       FREEZING RAIN LATENT HEAT FLUX (W M-2)
!   FLX3       PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
! ----------------------------------------------------------------------
!   SNOMLT     SNOW MELT (M) (WATER EQUIVALENT)
!   SNCOVR     FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
! ----------------------------------------------------------------------
!   RUNOFF1    SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
!   RUNOFF2    SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
!                SOIL LAYER (BASEFLOW)
!   RUNOFF3    NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
!                FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP (M S-1).
! Note: the above RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
! ----------------------------------------------------------------------
!   RC         CANOPY RESISTANCE (S M-1)
!   PC         PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
!                = ACTUAL TRANSP
!   XLAI       LEAF AREA INDEX (DIMENSIONLESS)
!   RSMIN      MINIMUM CANOPY RESISTANCE (S M-1)
!   RCS        INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
!   RCT        AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
!   RCQ        ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
!   RCSOIL     SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
! ----------------------------------------------------------------------
! 8. DIAGNOSTIC OUTPUT (D):
! ----------------------------------------------------------------------
!   SOILW      AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
!              BETWEEN SMCWLT AND SMCMAX)
!   SOILM      TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
!   Q1         Effective mixing ratio at surface (kg kg-1), used for
!              diagnosing the mixing ratio at 2 meter for coupled model
!   SMAV       Soil Moisture Availability for each layer, as a fraction 
!              between SMCWLT and SMCMAX.
!  Documentation for SNOTIME1 and SNOABL2 ?????
!  What categories of arguments do these variables fall into ????
!  Documentation for RIBB ?????
!  What category of argument does RIBB fall into ?????
! ----------------------------------------------------------------------
! 9. PARAMETERS (P):
! ----------------------------------------------------------------------
!   SMCWLT     WILTING POINT (VOLUMETRIC)
!   SMCDRY     DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
!                LAYER ENDS (VOLUMETRIC)
!   SMCREF     SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
!                STRESS (VOLUMETRIC)
!   SMCMAX     POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
!                (VOLUMETRIC)
!   NROOT      NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
!              IN SUBROUTINE REDPRM.
! ----------------------------------------------------------------------


      IMPLICIT NONE
! ----------------------------------------------------------------------

! DECLARATIONS - LOGICAL AND CHARACTERS
! ----------------------------------------------------------------------

      INTEGER, INTENT(IN) :: IILOC, JJLOC
      LOGICAL, INTENT(IN)::  LOCAL
      LOGICAL            ::  FRZGRA, SNOWNG
      CHARACTER (LEN=256), INTENT(IN)::  LLANDUSE, LSOIL

! ----------------------------------------------------------------------
! 1. CONFIGURATION INFORMATION (C):
! ----------------------------------------------------------------------
      INTEGER,INTENT(IN) ::  ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
      INTEGER, INTENT(IN) :: ISURBAN
      INTEGER,INTENT(OUT)::  NROOT
      INTEGER  KZ, K, iout

! ----------------------------------------------------------------------
! 2. LOGICAL:
! ----------------------------------------------------------------------
      LOGICAL, INTENT(IN) :: RDLAI2D
      LOGICAL, INTENT(IN) :: USEMONALB

      REAL, INTENT(IN)   :: SHDMIN,SHDMAX,DT,DQSDT2,LWDN,PRCP,PRCPRAIN,     &
                            Q2,Q2SAT,SFCPRS,SFCSPD,SFCTMP, SNOALB,          &
                            SOLDN,SOLNET,TBOT,TH2,ZLVL,                            &
                            FFROZP
      REAL, INTENT(OUT)  :: EMBRD
      REAL, INTENT(OUT)  :: ALBEDO
      REAL, INTENT(INOUT):: COSZ, SOLARDIRECT,CH,CM,                        &
                            CMC,SNEQV,SNCOVR,SNOWH,T1,XLAI,SHDFAC,Z0BRD,    &
                            EMISSI, ALB
      REAL, INTENT(INOUT):: SNOTIME1
      REAL, INTENT(INOUT):: RIBB
      REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SLDPTH
      REAL, DIMENSION(1:NSOIL), INTENT(OUT):: ET
      REAL, DIMENSION(1:NSOIL), INTENT(OUT):: SMAV
      REAL, DIMENSION(1:NSOIL), INTENT(INOUT) ::  SH2O, SMC, STC
      REAL,DIMENSION(1:NSOIL)::   RTDIS, ZSOIL

      REAL,INTENT(OUT)   :: ETA_KINEMATIC,BETA,DEW,DRIP,EC,EDIR,ESNOW,ETA,  &
                            ETP,FLX1,FLX2,FLX3,SHEAT,PC,RUNOFF1,RUNOFF2,    &
                            RUNOFF3,RC,RSMIN,RCQ,RCS,RCSOIL,RCT,SSOIL,      &
                            SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT, SOILM,      &
                            SOILW,FDOWN,Q1
      REAL :: BEXP,CFACTR,CMCMAX,CSOIL,CZIL,DF1,DF1H,DF1A,DKSAT,DWSAT,      &
              DSOIL,DTOT,ETT,FRCSNO,FRCSOI,EPSCA,F1,FXEXP,FRZX,HS,          &
              KDT,LVH2O,PRCP1,PSISAT,QUARTZ,R,RCH,REFKDT,RR,RGL,            &
              RSMAX,                                                        &
              RSNOW,SNDENS,SNCOND,SBETA,SN_NEW,SLOPE,SNUP,SALP,SOILWM,      &
              SOILWW,T1V,T24,T2V,TH2V,TOPT,TFREEZ,TSNOW,ZBOT,Z0,PRCPF,      &
              ETNS,PTU,LSUBS
        REAL ::  LVCOEF
      REAL :: INTERP_FRACTION
      REAL :: LAIMIN,    LAIMAX
      REAL :: ALBEDOMIN, ALBEDOMAX
      REAL :: EMISSMIN,  EMISSMAX
      REAL :: Z0MIN,     Z0MAX

! ----------------------------------------------------------------------
! DECLARATIONS - PARAMETERS
! ----------------------------------------------------------------------
      PARAMETER (TFREEZ = 273.15)
      PARAMETER (LVH2O = 2.501E+6)
      PARAMETER (LSUBS = 2.83E+6)
      PARAMETER (R = 287.04)

! ----------------------------------------------------------------------
!   INITIALIZATION
! ----------------------------------------------------------------------
      ILOC = IILOC
      JLOC = JJLOC

      RUNOFF1 = 0.0
      RUNOFF2 = 0.0
      RUNOFF3 = 0.0
      SNOMLT = 0.0

! ----------------------------------------------------------------------
!  THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE / LAND-ICE / ICE-FREE LAND
!     SEA-ICE CASE,          ICE =  1
!     NON-GLACIAL LAND,      ICE =  0
!     GLACIAL ICE,           ICE = -1
      IF (ICE == 1) call wrf_error_fatal("Sea-ice point in Noah-LSM")
      IF (ICE == -1) SHDFAC = 0.0

! ----------------------------------------------------------------------
! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
!   EACH SOIL LAYER.  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

! ----------------------------------------------------------------------
! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
! ----------------------------------------------------------------------
         CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX,TOPT,   &
                       REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,    &
                         PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT,          &
                         SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP,      &
                         RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL,CZIL,              &
                         LAIMIN, LAIMAX, EMISSMIN, EMISSMAX, ALBEDOMIN,    &
                         ALBEDOMAX, Z0MIN, Z0MAX, CSOIL, PTU, LLANDUSE,    &
                         LSOIL,LOCAL,LVCOEF)

!urban
         IF(VEGTYP==ISURBAN)THEN
              SHDFAC=0.05
              RSMIN=400.0
              SMCMAX = 0.45
              SMCREF = 0.42
              SMCWLT = 0.40
              SMCDRY = 0.40
         ENDIF

         IF ( SHDFAC >= SHDMAX ) THEN
            EMBRD = EMISSMAX
            IF (.NOT. RDLAI2D) THEN
               XLAI  = LAIMAX
            ENDIF
            IF (.NOT. USEMONALB) THEN
               ALB   = ALBEDOMIN
            ENDIF
            Z0BRD = Z0MAX
         ELSE IF ( SHDFAC <= SHDMIN ) THEN
            EMBRD = EMISSMIN
            IF(.NOT. RDLAI2D) THEN
               XLAI  = LAIMIN
            ENDIF
            IF(.NOT. USEMONALB) then
               ALB   = ALBEDOMAX
            ENDIF
            Z0BRD = Z0MIN
         ELSE

            IF ( SHDMAX > SHDMIN ) THEN

               INTERP_FRACTION = ( SHDFAC - SHDMIN ) / ( SHDMAX - SHDMIN )
               ! Bound INTERP_FRACTION between 0 and 1
               INTERP_FRACTION = MIN ( INTERP_FRACTION, 1.0 )
               INTERP_FRACTION = MAX ( INTERP_FRACTION, 0.0 )
               ! Scale Emissivity and LAI between EMISSMIN and EMISSMAX by INTERP_FRACTION
               EMBRD = ( ( 1.0 - INTERP_FRACTION ) * EMISSMIN  ) + ( INTERP_FRACTION * EMISSMAX  )
               IF (.NOT. RDLAI2D) THEN
                  XLAI  = ( ( 1.0 - INTERP_FRACTION ) * LAIMIN    ) + ( INTERP_FRACTION * LAIMAX    )
               ENDIF
               if (.not. USEMONALB) then
                  ALB   = ( ( 1.0 - INTERP_FRACTION ) * ALBEDOMAX ) + ( INTERP_FRACTION * ALBEDOMIN )
               endif
               Z0BRD = ( ( 1.0 - INTERP_FRACTION ) * Z0MIN     ) + ( INTERP_FRACTION * Z0MAX     )

            ELSE

               EMBRD = 0.5 * EMISSMIN  + 0.5 * EMISSMAX
               IF (.NOT. RDLAI2D) THEN
                  XLAI  = 0.5 * LAIMIN    + 0.5 * LAIMAX
               ENDIF
               if (.not. USEMONALB) then
                  ALB   = 0.5 * ALBEDOMIN + 0.5 * ALBEDOMAX
               endif
               Z0BRD = 0.5 * Z0MIN     + 0.5 * Z0MAX

            ENDIF

         ENDIF
! ----------------------------------------------------------------------
!  INITIALIZE PRECIPITATION LOGICALS.
! ----------------------------------------------------------------------
         SNOWNG = .FALSE.
         FRZGRA = .FALSE.

! ----------------------------------------------------------------------
         IF ( ICE == -1 ) THEN
            !
            ! FOR GLACIAL ICE, IF S.W.E. (SNEQV) BELOW THRESHOLD LOWER
            ! BOUND (0.10 M FOR GLACIAL ICE), THEN SET AT LOWER BOUND.
            !
            IF ( SNEQV < 0.10 ) THEN
               SNEQV = 0.10
               SNOWH = 0.50
            ENDIF
            !
            ! FOR GLACIAL ICE, SET SMC AND SH20 VALUES = 1.0
            !
            DO KZ = 1,NSOIL
               SMC(KZ) = 1.0
               SH2O(KZ) = 1.0
            END DO
         ENDIF
! ----------------------------------------------------------------------
! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
!   SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
!   SUBROUTINE)
! ----------------------------------------------------------------------
         IF ( SNEQV <= 1.E-7 ) THEN ! safer IF	kmh (2008/03/25)
            SNEQV = 0.0
            SNDENS = 0.0
            SNOWH = 0.0
            SNCOND = 1.0
         ELSE
            SNDENS = SNEQV / SNOWH
            IF(SNDENS > 1.0) THEN
             CALL wrf_error_fatal ( 'Physical snow depth is less than snow water equiv.' )
            ENDIF
            CALL CSNOW (SNCOND,SNDENS)
         END IF
! ----------------------------------------------------------------------
! DETERMINE IF IT'S PRECIPITATING 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
! snow defined when fraction of frozen precip (FFROZP) > 0.5,
! passed in from model microphysics.
            IF (FFROZP .GT. 0.5) THEN
               SNOWNG = .TRUE.
            ELSE
               IF (T1 <= TFREEZ) FRZGRA = .TRUE.
            END IF
         END IF
! ----------------------------------------------------------------------
! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
! IT TO THE EXISTING SNOWPACK.
! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
! ----------------------------------------------------------------------
         IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
            SN_NEW = PRCP * DT * 0.001
            SNEQV = SNEQV + SN_NEW
            PRCPF = 0.0

! ----------------------------------------------------------------------
! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
! UPDATE SNOW THERMAL CONDUCTIVITY
! ----------------------------------------------------------------------
            CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
!
! kmh 09/04/2006 set Snow Density at 0.2 g/cm**3
! for "cold permanent ice" or new "dry" snow
!
           IF ( (ICE == -1) .and. (SNCOVR .GT. 0.99) ) THEN
!  if soil temperature less than 268.15 K, treat as typical Antarctic/Greenland snow firn
              IF ( STC(1) .LT. (TFREEZ - 5.) ) SNDENS = 0.2
              IF ( SNOWNG .AND. (T1.LT.273.) .AND. (SFCTMP.LT.273.) ) SNDENS=0.2
           ENDIF
!
            CALL CSNOW (SNCOND,SNDENS)

! ----------------------------------------------------------------------
! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
! ANY CANOPY "DRIP" ADDED TO THIS LATER)
! ----------------------------------------------------------------------
         ELSE
            PRCPF = PRCP
         ENDIF
! ----------------------------------------------------------------------
! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
! ----------------------------------------------------------------------
         IF (SNEQV  == 0.0) THEN
            SNCOVR = 0.0
            ALBEDO = ALB
            EMISSI = EMBRD
         ELSE
! ----------------------------------------------------------------------
! DETERMINE SNOW FRACTIONAL COVERAGE.
! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
! ----------------------------------------------------------------------
            CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)

            ! kmh 2008/03/25
            ! Don't limit snow cover fraction over permanent ice
            IF ( ICE == 0 ) SNCOVR = MIN(SNCOVR,0.98)

            CALL ALCALC (ALB,SNOALB,EMBRD,SHDFAC,SHDMIN,SNCOVR,T1, &
                 ALBEDO,EMISSI,DT,SNOWNG,SNOTIME1,LVCOEF)
         ENDIF
! ----------------------------------------------------------------------
! THERMAL CONDUCTIVITY FOR GLACIAL ICE CASE
! ----------------------------------------------------------------------
         IF ( ICE == -1 ) THEN
            DF1 = 2.2
            !
            ! kmh 09/03/2006
            ! kmh 03/25/2008  change SNCOVR threshold to 0.97
            !
            ! only apply (small) DF1 conductivity for permanent land ice
            !
            IF (  SNCOVR .GT. 0.97 ) THEN
               DF1 = SNCOND
            ENDIF
         ELSEIF ( ICE == 0 ) THEN
! ----------------------------------------------------------------------
! 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)
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
! ----------------------------------------------------------------------
            CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))

!urban
            IF ( VEGTYP == ISURBAN ) DF1=3.24

            DF1 = DF1 * EXP (SBETA * SHDFAC)
!
! kmh 09/03/2006
! kmh 03/25/2008  change SNCOVR threshold to 0.97
!
            IF (  SNCOVR .GT. 0.97 ) THEN
               DF1 = SNCOND
            ENDIF
!
! ----------------------------------------------------------------------
! 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
! ----------------------------------------------------------------------
         END IF

         DSOIL = - (0.5 * ZSOIL (1))
         IF (SNEQV == 0.) THEN
            SSOIL = DF1 * (T1- STC (1) ) / DSOIL
         ELSE
            DTOT = SNOWH + DSOIL
            FRCSNO = SNOWH / DTOT

! 1. HARMONIC MEAN (SERIES FLOW)
!        DF1 = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
            FRCSOI = DSOIL / DTOT
! 2. ARITHMETIC MEAN (PARALLEL FLOW)
!        DF1 = FRCSNO*SNCOND + FRCSOI*DF1
            DF1H = (SNCOND * DF1)/ (FRCSOI * SNCOND+ FRCSNO * DF1)

! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
!        DF1 = (SNCOND**FRCSNO)*(DF1**FRCSOI)
! weigh DF by snow fraction
!        DF1 = DF1H*SNCOVR + DF1A*(1.0-SNCOVR)
!        DF1 = DF1H*SNCOVR + DF1*(1.0-SNCOVR)
            DF1A = FRCSNO * SNCOND+ FRCSOI * DF1

! ----------------------------------------------------------------------
! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
! MID-LAYER SOIL TEMPERATURE
! ----------------------------------------------------------------------
            DF1 = DF1A * SNCOVR + DF1* (1.0- SNCOVR)
            IF ( ICE == -1 ) then
               !  kmh  12/15/2005  correct for too deep snow layer
               !  kmh  09/03/2006  adjust DTOT
               IF ( DTOT .GT. 2.*DSOIL ) then
                  DTOT = 2.*DSOIL
               ENDIF
            ENDIF
            SSOIL = DF1 * (T1- STC (1) ) / DTOT
         END IF
! ----------------------------------------------------------------------
! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
! THE PREVIOUS TIMESTEP.
! ----------------------------------------------------------------------
         IF (SNCOVR  > 0. ) THEN
            CALL SNOWZ0 (SNCOVR,Z0,Z0BRD,SNOWH)
         ELSE
            Z0=Z0BRD
         END IF
! ----------------------------------------------------------------------
! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
! HEAT AND MOISTURE.

! NOTE !!!
! 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.

! NOTE !!!
! 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.

! NOTE !!!
! ----------------------------------------------------------------------
! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT. Needed as a state variable
! for iterative/implicit solution of CH in SFCDIF
! ----------------------------------------------------------------------
!        IF(.NOT.LCH) THEN
!          T1V = T1 * (1.0+ 0.61 * Q2)
!          TH2V = TH2 * (1.0+ 0.61 * Q2)
!          CALL SFCDIF_off (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
!        ENDIF

! ----------------------------------------------------------------------
! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
! CALCULATIONS.
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
! PENMAN EP SUBROUTINE THAT FOLLOWS
! ----------------------------------------------------------------------
!         FDOWN = SOLDN * (1.0- ALBEDO) + LWDN
         FDOWN =  SOLNET + LWDN
! ----------------------------------------------------------------------
! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
! PENMAN.
         T2V = SFCTMP * (1.0+ 0.61 * Q2 )

         iout=0
         if(iout.eq.1) then
         print*,'before penman'
         print*,' SFCTMP',SFCTMP,'SFCPRS',SFCPRS,'CH',CH,'T2V',T2V,      &
       'TH2',TH2,'PRCP',PRCP,'FDOWN',FDOWN,'T24',T24,'SSOIL',SSOIL,      &
        'Q2',Q2,'Q2SAT',Q2SAT,'ETP',ETP,'RCH',RCH,                       &
        'EPSCA',EPSCA,'RR',RR  ,'SNOWNG',SNOWNG,'FRZGRA',FRZGRA,           &
        'DQSDT2',DQSDT2,'FLX2',FLX2,'SNOWH',SNOWH,'SNEQV',SNEQV,         &
        ' DSOIL',DSOIL,' FRCSNO',FRCSNO,' SNCOVR',SNCOVR,' DTOT',DTOT,   &
       ' ZSOIL (1)',ZSOIL(1),' DF1',DF1,'T1',T1,' STC1',STC(1),          &
        'ALBEDO',ALBEDO,'SMC',SMC,'STC',STC,'SH2O',SH2O
         endif

         CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL,     &
                       Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA,          &
!
! kmh 01/09/2007 add T1,ICE,SNCOVR to call
!
                         DQSDT2,FLX2,EMISSI,SNEQV,T1,ICE,SNCOVR)
!
! ----------------------------------------------------------------------
! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
! IF NONZERO GREENNESS FRACTION
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
!  FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
!  BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
! ----------------------------------------------------------------------
         IF (SHDFAC > 0.) THEN
            CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL,     &
                          SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2,  &
                          TOPT,RSMAX,RGL,HS,XLAI,                        &
                          RCS,RCT,RCQ,RCSOIL,EMISSI)
         ELSE
            RC = 0.0
         END IF
! ----------------------------------------------------------------------
! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
! EXISTS OR NOT:
! ----------------------------------------------------------------------
         ESNOW = 0.0
         IF (SNEQV  == 0.0) THEN
            CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT,                  &
                            SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,           &
                            SHDFAC,                                      &
                            SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI,  &
                            SSOIL,                                       &
                            STC,EPSCA,BEXP,PC,RCH,RR,CFACTR,             &
                            SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL,            &
                            DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2,       &
                            RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS,      &
                            QUARTZ,FXEXP,CSOIL,                          &
                            BETA,DRIP,DEW,FLX1,FLX3,VEGTYP,ISURBAN)
            ETA_KINEMATIC = ETA
         ELSE
            CALL SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT,    &
                         SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,              &
                         SBETA,DF1,                                      &
                         Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA,  &
                         SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS,&
                         SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT,               &
                         ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1,     &
                         RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT,    &
                         ICE,RTDIS,QUARTZ,FXEXP,CSOIL,                   &
                         BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI, &
                         RIBB,SOLDN,                                     &
                         ISURBAN,                                        &
                         VEGTYP)
            ETA_KINEMATIC =  ESNOW + ETNS
         END IF

!     Calculate effective mixing ratio at grnd level (skin)
!
!     Q1=Q2+ETA*CP/RCH
     Q1=Q2+ETA_KINEMATIC*CP/RCH
!
! ----------------------------------------------------------------------
! DETERMINE SENSIBLE HEAT (H) IN ENERGY UNITS (W M-2)
! ----------------------------------------------------------------------
         SHEAT = - (CH * CP * SFCPRS)/ (R * T2V) * ( TH2- T1 )

! ----------------------------------------------------------------------
! CONVERT EVAP TERMS FROM KINEMATIC (KG M-2 S-1) TO ENERGY UNITS (W M-2)
! ----------------------------------------------------------------------
      EDIR = EDIR * LVH2O
      EC = EC * LVH2O
      DO K=1,4
      ET(K) = ET(K) * LVH2O
      ENDDO
      ETT = ETT * LVH2O
      ESNOW = ESNOW * LSUBS
      ETP = ETP*((1.-SNCOVR)*LVH2O + SNCOVR*LSUBS)
      IF (ETP .GT. 0.) THEN
         ETA = EDIR + EC + ETT + ESNOW
      ELSE
        ETA = ETP
      ENDIF
! ----------------------------------------------------------------------
! DETERMINE BETA (RATIO OF ACTUAL TO POTENTIAL EVAP)
! ----------------------------------------------------------------------
      IF (ETP == 0.0) THEN
        BETA = 0.0
      ELSE
        BETA = ETA/ETP
      ENDIF

! ----------------------------------------------------------------------
! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
!   SSOIL>0: WARM THE SURFACE  (NIGHT TIME)
!   SSOIL<0: COOL THE SURFACE  (DAY TIME)
! ----------------------------------------------------------------------
         SSOIL = -1.0* SSOIL

! ----------------------------------------------------------------------
!  FOR THE CASE OF LAND (BUT NOT GLACIAL ICE):
!  CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
!  AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW.  RUNOFF2 IS ALREADY
!  A RATE AT THIS POINT
! ----------------------------------------------------------------------
         IF (ICE == 0) THEN
            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
            SOILWM = -1.0* (SMCMAX - SMCWLT)* ZSOIL (1)
            SOILWW = -1.0* (SMC (1) - SMCWLT)* ZSOIL (1)
!
            DO K = 1,NSOIL
               SMAV(K)=(SMC(K) - SMCWLT)/(SMCMAX - SMCWLT)
            END DO

            IF (NROOT >= 2) THEN
              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
            END IF
            IF (SOILWM .LT. 1.E-6) THEN
              SOILWM = 0.0
              SOILW  = 0.0
              SOILM  = 0.0
            ELSE
              SOILW = SOILWW / SOILWM
            END IF
         ELSEIF ( ICE == -1 ) THEN
! ----------------------------------------------------------------------
! FOR THE CASE OF GLACIAL ICE (ICE == -1), ADD ANY SNOWMELT DIRECTLY TO
! SURFACE RUNOFF (RUNOFF1) SINCE THERE IS NO SOIL MEDIUM, AND THUS NO
! CALL TO SUBROUTINE SMFLX (FOR SOIL MOISTURE TENDENCY).
! ----------------------------------------------------------------------
            RUNOFF1 = SNOMLT/DT
           SOILWM = 0.0
           SOILW  = 0.0
           SOILM  = 0.0
           DO K = 1,NSOIL
             SMAV(K)= 1.0
           END DO
         END IF

! ----------------------------------------------------------------------
  END SUBROUTINE SFLX
! ----------------------------------------------------------------------

      SUBROUTINE ALCALC (ALB,SNOALB,EMBRD,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO,EMISSI,   &
                         DT,SNOWNG,SNOTIME1,LVCOEF)

! ----------------------------------------------------------------------
! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
!   ALB     SNOWFREE ALBEDO
!   SNOALB  MAXIMUM (DEEP) SNOW ALBEDO
!   SHDFAC    AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
!   SHDMIN    MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
!   SNCOVR  FRACTIONAL SNOW COVER
!   ALBEDO  SURFACE ALBEDO INCLUDING SNOW EFFECT
!   TSNOW   SNOW SURFACE TEMPERATURE (K)
! ----------------------------------------------------------------------
      IMPLICIT NONE

! ----------------------------------------------------------------------
! 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)
! ----------------------------------------------------------------------
      REAL, INTENT(IN)  ::  ALB, SNOALB, EMBRD, SHDFAC, SHDMIN, SNCOVR, TSNOW
      REAL, INTENT(IN)  :: DT
      LOGICAL, INTENT(IN) :: SNOWNG
      REAL, INTENT(INOUT):: SNOTIME1
      REAL, INTENT(OUT) ::  ALBEDO, EMISSI
      REAL              :: SNOALB2
      REAL              :: TM,SNOALB1
      REAL, INTENT(IN)  :: LVCOEF
      REAL, PARAMETER   :: SNACCA=0.94,SNACCB=0.58,SNTHWA=0.82,SNTHWB=0.46
! turn of vegetation effect
!      ALBEDO = ALB + (1.0- (SHDFAC - SHDMIN))* SNCOVR * (SNOALB - ALB)
!      ALBEDO = (1.0-SNCOVR)*ALB + SNCOVR*SNOALB !this is equivalent to below
      ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
      EMISSI = EMBRD + SNCOVR*(EMISSI_S - EMBRD)

!     BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
!          IF (TSNOW.LE.263.16) THEN
!            ALBEDO=SNOALB
!          ELSE
!            IF (TSNOW.LT.273.16) THEN
!              TM=0.1*(TSNOW-263.16)
!              SNOALB1=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
!            ELSE
!              SNOALB1=0.67
!             IF(SNCOVR.GT.0.95) SNOALB1= 0.6
!             SNOALB1 = ALB + SNCOVR*(SNOALB-ALB)
!            ENDIF
!          ENDIF
!            ALBEDO = ALB + SNCOVR*(SNOALB1-ALB)

!     ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
!          SNOALB1 = SNOALB+COEF*(0.85-SNOALB)
!          SNOALB2=SNOALB1
!!m          LSTSNW=LSTSNW+1
!          SNOTIME1 = SNOTIME1 + DT
!          IF (SNOWNG) THEN
!             SNOALB2=SNOALB
!!m             LSTSNW=0
!             SNOTIME1 = 0.0
!          ELSE
!            IF (TSNOW.LT.273.16) THEN
!!              SNOALB2=SNOALB-0.008*LSTSNW*DT/86400
!!m              SNOALB2=SNOALB-0.008*SNOTIME1/86400
!              SNOALB2=(SNOALB2-0.65)*EXP(-0.05*DT/3600)+0.65
!!              SNOALB2=(ALBEDO-0.65)*EXP(-0.01*DT/3600)+0.65
!            ELSE
!              SNOALB2=(SNOALB2-0.5)*EXP(-0.0005*DT/3600)+0.5
!!              SNOALB2=(SNOALB-0.5)*EXP(-0.24*LSTSNW*DT/86400)+0.5
!!m              SNOALB2=(SNOALB-0.5)*EXP(-0.24*SNOTIME1/86400)+0.5
!            ENDIF
!          ENDIF
!
!!               print*,'SNOALB2',SNOALB2,'ALBEDO',ALBEDO,'DT',DT
!          ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
!          IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
!!m          LSTSNW1=LSTSNW
!!          SNOTIME = SNOTIME1

! formulation by Livneh
! ----------------------------------------------------------------------
! SNOALB IS CONSIDERED AS THE MAXIMUM SNOW ALBEDO FOR NEW SNOW, AT
! A VALUE OF 85%. SNOW ALBEDO CURVE DEFAULTS ARE FROM BRAS P.263. SHOULD
! NOT BE CHANGED EXCEPT FOR SERIOUS PROBLEMS WITH SNOW MELT.
! TO IMPLEMENT ACCUMULATIN PARAMETERS, SNACCA AND SNACCB, ASSERT THAT IT
! IS INDEED ACCUMULATION SEASON. I.E. THAT SNOW SURFACE TEMP IS BELOW
! ZERO AND THE DATE FALLS BETWEEN OCTOBER AND FEBRUARY
! ----------------------------------------------------------------------
         SNOALB1 = SNOALB+LVCOEF*(0.85-SNOALB)
         SNOALB2=SNOALB1
! ---------------- Initial LSTSNW --------------------------------------
          IF (SNOWNG) THEN
             SNOTIME1 = 0.
          ELSE
           SNOTIME1=SNOTIME1+DT
!               IF (TSNOW.LT.273.16) THEN
                   SNOALB2=SNOALB1*(SNACCA**((SNOTIME1/86400.0)**SNACCB))
!               ELSE
!                  SNOALB2 =SNOALB1*(SNTHWA**((SNOTIME1/86400.0)**SNTHWB))
!               ENDIF
          ENDIF
!
           SNOALB2 = MAX ( SNOALB2, ALB )
           ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
           IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2

!          IF (TSNOW.LT.273.16) THEN
!            ALBEDO=SNOALB-0.008*DT/86400
!          ELSE
!            ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
!          ENDIF

!      IF (ALBEDO > SNOALB) ALBEDO = SNOALB

! ----------------------------------------------------------------------
  END SUBROUTINE ALCALC
! ----------------------------------------------------------------------

      SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL,       &
                         SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2,    &
                         TOPT,RSMAX,RGL,HS,XLAI,                          &
                         RCS,RCT,RCQ,RCSOIL,EMISSI)

! ----------------------------------------------------------------------
! SUBROUTINE CANRES
! ----------------------------------------------------------------------
! CALCULATE 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
! ----------------------------------------------------------------------
! INPUT:
!   SOLAR   INCOMING SOLAR RADIATION
!   CH      SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
!   SFCTMP  AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
!   Q2      AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
!   Q2SAT   SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
!   DQSDT2  SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
!   SFCPRS  SURFACE PRESSURE
!   SMC     VOLUMETRIC SOIL MOISTURE
!   ZSOIL   SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
!   NSOIL   NO. OF SOIL LAYERS
!   NROOT   NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
!   XLAI    LEAF AREA INDEX
!   SMCWLT  WILTING POINT
!   SMCREF  REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
!             SETS IN)
! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
!   SURBOUTINE REDPRM
! OUTPUT:
!   PC  PLANT COEFFICIENT
!   RC  CANOPY RESISTANCE
! ----------------------------------------------------------------------

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NROOT,NSOIL
      INTEGER  K
      REAL,    INTENT(IN) :: CH,DQSDT2,HS,Q2,Q2SAT,RSMIN,RGL,RSMAX,        &
                             SFCPRS,SFCTMP,SMCREF,SMCWLT, SOLAR,TOPT,XLAI, &
                             EMISSI
      REAL,DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
      REAL,    INTENT(OUT):: PC,RC,RCQ,RCS,RCSOIL,RCT
      REAL                :: DELTA,FF,GX,P,RR
      REAL, DIMENSION(1:NSOIL) ::  PART
      REAL, PARAMETER     :: SLV = 2.501000E6


! ----------------------------------------------------------------------
! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
! ----------------------------------------------------------------------
      RCS = 0.0
      RCT = 0.0
      RCQ = 0.0
      RCSOIL = 0.0

! ----------------------------------------------------------------------
! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
! ----------------------------------------------------------------------
      RC = 0.0
      FF = 0.55*2.0* SOLAR / (RGL * XLAI)
      RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)

! ----------------------------------------------------------------------
! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
! ----------------------------------------------------------------------
      RCS = MAX (RCS,0.0001)
      RCT = 1.0- 0.0016* ( (TOPT - SFCTMP)**2.0)

! ----------------------------------------------------------------------
! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
! RCQ EXPRESSION FROM SSIB
! ----------------------------------------------------------------------
      RCT = MAX (RCT,0.0001)
      RCQ = 1.0/ (1.0+ HS * (Q2SAT - Q2))

! ----------------------------------------------------------------------…