wrfv2_fire/phys/module_sf_noahlsm.F

   1 MODULE module_sf_noahlsm
   2   USE module_model_constants
   3
   4 !   REAL, PARAMETER    :: CP = 1004.5
   5   REAL, PARAMETER      :: RD = 287.04, SIGMA = 5.67E-8,                 &
   6                           CPH2O = 4.218E+3,CPICE = 2.106E+3,            &
   7                           LSUBF = 3.335E+5,                             &
   8                           EMISSI_S = 0.95
   9
  10 ! VEGETATION PARAMETERS
  11         INTEGER :: LUCATS , BARE
  12         INTEGER :: NATURAL
  13         integer, PARAMETER :: NLUS=50
  14         CHARACTER(LEN=256) LUTYPE
  15         INTEGER, DIMENSION(1:NLUS) :: NROTBL
  16         real, dimension(1:NLUS) ::  SNUPTBL, RSTBL, RGLTBL, HSTBL,                &
  17                                     SHDTBL, MAXALB,                               &
  18                                     EMISSMINTBL, EMISSMAXTBL,                     &
  19                                     LAIMINTBL, LAIMAXTBL,                         &
  20                                     Z0MINTBL, Z0MAXTBL,                           &
  21                                     ALBEDOMINTBL, ALBEDOMAXTBL
  22         REAL ::   TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA
  23
  24 ! SOIL PARAMETERS
  25         INTEGER :: SLCATS
  26         INTEGER, PARAMETER :: NSLTYPE=30
  27         CHARACTER(LEN=256) SLTYPE
  28         REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,F11,                           &
  29         MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ
  30
  31 ! LSM GENERAL PARAMETERS
  32         INTEGER :: SLPCATS
  33         INTEGER, PARAMETER :: NSLOPE=30
  34         REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
  35         REAL ::  SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA,           &
  36                  REFKDT_DATA,FRZK_DATA,ZBOT_DATA,  SMLOW_DATA,SMHIGH_DATA,        &
  37                         CZIL_DATA
  38         REAL ::  LVCOEF_DATA
  39
  40         CHARACTER*256  :: err_message
  41
  42 !
  43 CONTAINS
  44 !
  45
  46       SUBROUTINE SFLX (FFROZP,ICE,ISURBAN,DT,ZLVL,NSOIL,SLDPTH,         &    !C
  47                        LOCAL,                                           &    !L
  48                        LLANDUSE, LSOIL,                                 &    !CL
  49                        LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2,SFCSPD,  &    !F
  50                        COSZ,PRCPRAIN, SOLARDIRECT,                      &    !F
  51                        TH2,Q2SAT,DQSDT2,                                &    !I
  52                        VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHDMIN,SHDMAX,    &    !I
  53                        ALB, SNOALB,TBOT, Z0BRD, Z0, EMISSI, EMBRD,      &    !S
  54                        CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM,    &    !H
  55 ! ----------------------------------------------------------------------
  56 ! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
  57 ! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
  58 ! MODEL).  OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
  59 ! ----------------------------------------------------------------------
  60                        ETA,SHEAT, ETA_KINEMATIC,FDOWN,                  &    !O
  61                        EC,EDIR,ET,ETT,ESNOW,DRIP,DEW,                   &    !O
  62                        BETA,ETP,SSOIL,                                  &    !O
  63                        FLX1,FLX2,FLX3,                                  &    !O
  64                        SNOMLT,SNCOVR,                                   &    !O
  65                        RUNOFF1,RUNOFF2,RUNOFF3,                         &    !O
  66                        RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL,             &    !O
  67                        SOILW,SOILM,Q1,                                  &    !D
  68                        RDLAI2D,USEMONALB,                               &
  69                        SNOTIME1,                                        &
  70                        RIBB,                                            &
  71                        SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT)                    !P
  72 ! ----------------------------------------------------------------------
  73 ! SUBROUTINE SFLX - UNIFIED NOAHLSM VERSION 1.0 JULY 2007
  74 ! ----------------------------------------------------------------------
  75 ! SUB-DRIVER FOR "Noah LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
  76 ! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL MOISTURE, SOIL
  77 ! ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, SNOWPACK WATER CONTENT,
  78 ! SNOWDEPTH, AND ALL TERMS OF THE SURFACE ENERGY BALANCE AND SURFACE
  79 ! WATER BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF DOWNWARD
  80 ! RADIATION AND PRECIP)
  81 ! ----------------------------------------------------------------------
  82 ! SFLX ARGUMENT LIST KEY:
  83 ! ----------------------------------------------------------------------
  84 !  C  CONFIGURATION INFORMATION
  85 !  L  LOGICAL
  86 ! CL  4-string character bearing logical meaning
  87 !  F  FORCING DATA
  88 !  I  OTHER (INPUT) FORCING DATA
  89 !  S  SURFACE CHARACTERISTICS
  90 !  H  HISTORY (STATE) VARIABLES
  91 !  O  OUTPUT VARIABLES
  92 !  D  DIAGNOSTIC OUTPUT
  93 !  P  Parameters
  94 !  Msic Miscellaneous terms passed from gridded driver
  95 ! ----------------------------------------------------------------------
  96 ! 1. CONFIGURATION INFORMATION (C):
  97 ! ----------------------------------------------------------------------
  98 !   ICE        SEA-ICE FLAG  (=1: SEA-ICE, =0: LAND (NO ICE), --1 LAND-ICE).
  99 !   DT         TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
 100 !                1800 SECS OR LESS)
 101 !   ZLVL       HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
 102 !   NSOIL      NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
 103 !                PARAMETER NSOLD SET BELOW)
 104 !   SLDPTH     THE THICKNESS OF EACH SOIL LAYER (M)
 105 ! ----------------------------------------------------------------------
 106 ! 2. LOGICAL:
 107 ! ----------------------------------------------------------------------
 108 !   LCH       Exchange coefficient (Ch) calculation flag (false: using
 109 !                ch-routine SFCDIF; true: Ch is brought in)
 110 !   LOCAL      Flag for local-site simulation (where there is no
 111 !              maps for albedo, veg fraction, and roughness
 112 !             true:  all LSM parameters (inluding albedo, veg fraction and
 113 !                    roughness length) will be defined by three tables
 114 !   LLANDUSE  (=USGS, using USGS landuse classification)
 115 !   LSOIL     (=STAS, using FAO/STATSGO soil texture classification)
 116 ! ----------------------------------------------------------------------
 117 ! 3. FORCING DATA (F):
 118 ! ----------------------------------------------------------------------
 119 !   LWDN       LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
 120 !   SOLDN      SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
 121 !   SOLNET     NET DOWNWARD SOLAR RADIATION ((W M-2; POSITIVE)
 122 !   SFCPRS     PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
 123 !   PRCP       PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
 124 !   SFCTMP     AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
 125 !   TH2        AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
 126 !   Q2         MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
 127 !   COSZ       Solar zenith angle (not used for now)
 128 !   PRCPRAIN   Liquid-precipitation rate (KG M-2 S-1) (not used)
 129 ! SOLARDIRECT  Direct component of downward solar radiation (W M-2) (not used)
 130 !   FFROZP     FRACTION OF FROZEN PRECIPITATION
 131 ! ----------------------------------------------------------------------
 132 ! 4. OTHER FORCING (INPUT) DATA (I):
 133 ! ----------------------------------------------------------------------
 134 !   SFCSPD     WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
 135 !   Q2SAT      SAT SPECIFIC HUMIDITY AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
 136 !   DQSDT2     SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
 137 !                (KG KG-1 K-1)
 138 ! ----------------------------------------------------------------------
 139 ! 5. CANOPY/SOIL CHARACTERISTICS (S):
 140 ! ----------------------------------------------------------------------
 141 !   VEGTYP     VEGETATION TYPE (INTEGER INDEX)
 142 !   SOILTYP    SOIL TYPE (INTEGER INDEX)
 143 !   SLOPETYP   CLASS OF SFC SLOPE (INTEGER INDEX)
 144 !   SHDFAC     AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
 145 !                (FRACTION= 0.0-1.0)
 146 !   SHDMIN     MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
 147 !                (FRACTION= 0.0-1.0) <= SHDFAC
 148 !   PTU        PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
 149 !                (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
 150 !                VEG PARMS)
 151 !   ALB        BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
 152 !                DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
 153 !                MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
 154 !                INCLUDE DIURNAL SUN ANGLE EFFECT)
 155 !   SNOALB     UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
 156 !                ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
 157 !   TBOT       BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
 158 !                TEMPERATURE)
 159 !   Z0BRD      Background fixed roughness length (M)
 160 !   Z0         Time varying roughness length (M) as function of snow depth
 161 !
 162 !   EMBRD      Background surface emissivity (between 0 and 1)
 163 !   EMISSI     Surface emissivity (between 0 and 1)
 164 ! ----------------------------------------------------------------------
 165 ! 6. HISTORY (STATE) VARIABLES (H):
 166 ! ----------------------------------------------------------------------
 167 !  CMC         CANOPY MOISTURE CONTENT (M)
 168 !  T1          GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
 169 !  STC(NSOIL)  SOIL TEMP (K)
 170 !  SMC(NSOIL)  TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
 171 !  SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
 172 !                NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
 173 !  SNOWH       ACTUAL SNOW DEPTH (M)
 174 !  SNEQV       LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
 175 !                NOTE: SNOW DENSITY = SNEQV/SNOWH
 176 !  ALBEDO      SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
 177 !                =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
 178 !                =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
 179 !  CH          SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
 180 !                (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
 181 !                IT HAS BEEN MULTIPLIED BY WIND SPEED.
 182 !  CM          SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
 183 !                CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
 184 !                MULTIPLIED BY WIND SPEED.
 185 ! ----------------------------------------------------------------------
 186 ! 7. OUTPUT (O):
 187 ! ----------------------------------------------------------------------
 188 ! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
 189 ! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL.  FOR THIS APPLICATION,
 190 ! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
 191 ! NECESSARY.  OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
 192 !   ETA        ACTUAL LATENT HEAT FLUX (W m-2: NEGATIVE, IF UP FROM
 193 !              SURFACE)
 194 !  ETA_KINEMATIC atctual latent heat flux in Kg m-2 s-1
 195 !   SHEAT      SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
 196 !              SURFACE)
 197 !   FDOWN      Radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
 198 ! ----------------------------------------------------------------------
 199 !   EC         CANOPY WATER EVAPORATION (W m-2)
 200 !   EDIR       DIRECT SOIL EVAPORATION (W m-2)
 201 !   ET(NSOIL)  PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
 202 !                 (W m-2)
 203 !   ETT        TOTAL PLANT TRANSPIRATION (W m-2)
 204 !   ESNOW      SUBLIMATION FROM (OR DEPOSITION TO IF <0) SNOWPACK
 205 !                (W m-2)
 206 !   DRIP       THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
 207 !                WATER-HOLDING CAPACITY (M)
 208 !   DEW        DEWFALL (OR FROSTFALL FOR T<273.15) (M)
 209 ! ----------------------------------------------------------------------
 210 !   BETA       RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
 211 !   ETP        POTENTIAL EVAPORATION (W m-2)
 212 !   SSOIL      SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
 213 ! ----------------------------------------------------------------------
 214 !   FLX1       PRECIP-SNOW SFC (W M-2)
 215 !   FLX2       FREEZING RAIN LATENT HEAT FLUX (W M-2)
 216 !   FLX3       PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
 217 ! ----------------------------------------------------------------------
 218 !   SNOMLT     SNOW MELT (M) (WATER EQUIVALENT)
 219 !   SNCOVR     FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
 220 ! ----------------------------------------------------------------------
 221 !   RUNOFF1    SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
 222 !   RUNOFF2    SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
 223 !                SOIL LAYER (BASEFLOW)
 224 !   RUNOFF3    NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
 225 !                FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP (M S-1).
 226 ! Note: the above RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
 227 ! ----------------------------------------------------------------------
 228 !   RC         CANOPY RESISTANCE (S M-1)
 229 !   PC         PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
 230 !                = ACTUAL TRANSP
 231 !   XLAI       LEAF AREA INDEX (DIMENSIONLESS)
 232 !   RSMIN      MINIMUM CANOPY RESISTANCE (S M-1)
 233 !   RCS        INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
 234 !   RCT        AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
 235 !   RCQ        ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
 236 !   RCSOIL     SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
 237 ! ----------------------------------------------------------------------
 238 ! 8. DIAGNOSTIC OUTPUT (D):
 239 ! ----------------------------------------------------------------------
 240 !   SOILW      AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
 241 !              BETWEEN SMCWLT AND SMCMAX)
 242 !   SOILM      TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
 243 !   Q1         Effective mixing ratio at surface (kg kg-1), used for
 244 !              diagnosing the mixing ratio at 2 meter for coupled model
 245 !  Documentation for SNOTIME1 and SNOABL2 ?????
 246 !  What categories of arguments do these variables fall into ????
 247 !  Documentation for RIBB ?????
 248 !  What category of argument does RIBB fall into ?????
 249 ! ----------------------------------------------------------------------
 250 ! 9. PARAMETERS (P):
 251 ! ----------------------------------------------------------------------
 252 !   SMCWLT     WILTING POINT (VOLUMETRIC)
 253 !   SMCDRY     DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
 254 !                LAYER ENDS (VOLUMETRIC)
 255 !   SMCREF     SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
 256 !                STRESS (VOLUMETRIC)
 257 !   SMCMAX     POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
 258 !                (VOLUMETRIC)
 259 !   NROOT      NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
 260 !              IN SUBROUTINE REDPRM.
 261 ! ----------------------------------------------------------------------
 262
 263
 264       IMPLICIT NONE
 265 ! ----------------------------------------------------------------------
 266
 267 ! DECLARATIONS - LOGICAL AND CHARACTERS
 268 ! ----------------------------------------------------------------------
 269       LOGICAL, INTENT(IN)::  LOCAL
 270       LOGICAL            ::  FRZGRA, SNOWNG
 271       CHARACTER (LEN=256), INTENT(IN)::  LLANDUSE, LSOIL
 272
 273 ! ----------------------------------------------------------------------
 274 ! 1. CONFIGURATION INFORMATION (C):
 275 ! ----------------------------------------------------------------------
 276       INTEGER,INTENT(IN) ::  ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
 277       INTEGER, INTENT(IN) :: ISURBAN
 278       INTEGER,INTENT(OUT)::  NROOT
 279       INTEGER  KZ, K, iout
 280
 281 ! ----------------------------------------------------------------------
 282 ! 2. LOGICAL:
 283 ! ----------------------------------------------------------------------
 284       LOGICAL, INTENT(IN) :: RDLAI2D
 285       LOGICAL, INTENT(IN) :: USEMONALB
 286
 287       REAL, INTENT(IN)   :: SHDMIN,SHDMAX,DT,DQSDT2,LWDN,PRCP,PRCPRAIN,     &
 288                             Q2,Q2SAT,SFCPRS,SFCSPD,SFCTMP, SNOALB,          &
 289                             SOLDN,SOLNET,TBOT,TH2,ZLVL,                            &
 290                             FFROZP
 291       REAL, INTENT(OUT)  :: EMBRD
 292       REAL, INTENT(OUT)  :: ALBEDO
 293       REAL, INTENT(INOUT):: COSZ, SOLARDIRECT,CH,CM,                        &
 294                             CMC,SNEQV,SNCOVR,SNOWH,T1,XLAI,SHDFAC,Z0BRD,    &
 295                             EMISSI, ALB
 296       REAL, INTENT(INOUT):: SNOTIME1
 297       REAL, INTENT(INOUT):: RIBB
 298       REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SLDPTH
 299       REAL, DIMENSION(1:NSOIL), INTENT(OUT):: ET
 300       REAL, DIMENSION(1:NSOIL), INTENT(INOUT) ::  SH2O, SMC, STC
 301       REAL,DIMENSION(1:NSOIL)::   RTDIS, ZSOIL
 302
 303       REAL,INTENT(OUT)   :: ETA_KINEMATIC,BETA,DEW,DRIP,EC,EDIR,ESNOW,ETA,  &
 304                             ETP,FLX1,FLX2,FLX3,SHEAT,PC,RUNOFF1,RUNOFF2,    &
 305                             RUNOFF3,RC,RSMIN,RCQ,RCS,RCSOIL,RCT,SSOIL,      &
 306                             SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT, SOILM,      &
 307                             SOILW,FDOWN,Q1
 308       REAL :: BEXP,CFACTR,CMCMAX,CSOIL,CZIL,DF1,DF1H,DF1A,DKSAT,DWSAT,      &
 309               DSOIL,DTOT,ETT,FRCSNO,FRCSOI,EPSCA,F1,FXEXP,FRZX,HS,          &
 310               KDT,LVH2O,PRCP1,PSISAT,QUARTZ,R,RCH,REFKDT,RR,RGL,            &
 311               RSMAX,                                                        &
 312               RSNOW,SNDENS,SNCOND,SBETA,SN_NEW,SLOPE,SNUP,SALP,SOILWM,      &
 313               SOILWW,T1V,T24,T2V,TH2V,TOPT,TFREEZ,TSNOW,ZBOT,Z0,PRCPF,      &
 314               ETNS,PTU,LSUBS
 315         REAL ::  LVCOEF
 316       REAL :: INTERP_FRACTION
 317       REAL :: LAIMIN,    LAIMAX
 318       REAL :: ALBEDOMIN, ALBEDOMAX
 319       REAL :: EMISSMIN,  EMISSMAX
 320       REAL :: Z0MIN,     Z0MAX
 321
 322 ! ----------------------------------------------------------------------
 323 ! DECLARATIONS - PARAMETERS
 324 ! ----------------------------------------------------------------------
 325       PARAMETER (TFREEZ = 273.15)
 326       PARAMETER (LVH2O = 2.501E+6)
 327       PARAMETER (LSUBS = 2.83E+6)
 328       PARAMETER (R = 287.04)
 329 ! ----------------------------------------------------------------------
 330 !   INITIALIZATION
 331 ! ----------------------------------------------------------------------
 332          RUNOFF1 = 0.0
 333          RUNOFF2 = 0.0
 334          RUNOFF3 = 0.0
 335          SNOMLT = 0.0
 336
 337 ! ----------------------------------------------------------------------
 338 !  THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE / LAND-ICE / ICE-FREE LAND
 339 !     SEA-ICE CASE,          ICE =  1
 340 !     NON-GLACIAL LAND,      ICE =  0
 341 !     GLACIAL-ICE LAND,      ICE = -1
 342       IF (ICE /= 0) SHDFAC = 0.0
 343 ! ----------------------------------------------------------------------
 344 ! SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS
 345 ! ----------------------------------------------------------------------
 346          IF (ICE == 1) THEN
 347             DO KZ = 1,NSOIL
 348                ZSOIL (KZ) = -3.* FLOAT (KZ)/ FLOAT (NSOIL)
 349             END DO
 350
 351 ! ----------------------------------------------------------------------
 352 ! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
 353 !   EACH SOIL LAYER.  NOTE:  SIGN OF ZSOIL IS NEGATIVE (DENOTING BELOW
 354 !   GROUND)
 355 ! ----------------------------------------------------------------------
 356          ELSE
 357             ZSOIL (1) = - SLDPTH (1)
 358             DO KZ = 2,NSOIL
 359                ZSOIL (KZ) = - SLDPTH (KZ) + ZSOIL (KZ -1)
 360             END DO
 361          END IF
 362 ! ----------------------------------------------------------------------
 363 ! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
 364 ! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
 365 ! ----------------------------------------------------------------------
 366          CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX,TOPT,   &
 367                        REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,    &
 368                          PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT,          &
 369                          SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP,      &
 370                          RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL,CZIL,              &
 371                          LAIMIN, LAIMAX, EMISSMIN, EMISSMAX, ALBEDOMIN,    &
 372                          ALBEDOMAX, Z0MIN, Z0MAX, CSOIL, PTU, LLANDUSE,    &
 373                          LSOIL,LOCAL,LVCOEF)
 374
 375 !urban
 376          IF(VEGTYP==ISURBAN)THEN
 377               SHDFAC=0.05
 378               RSMIN=400.0
 379               SMCMAX = 0.45
 380               SMCREF = 0.42
 381               SMCWLT = 0.40
 382               SMCDRY = 0.40
 383          ENDIF
 384
 385          IF ( SHDFAC >= SHDMAX ) THEN
 386             EMBRD = EMISSMAX
 387             IF (.NOT. RDLAI2D) THEN
 388                XLAI  = LAIMAX
 389             ENDIF
 390             IF (.NOT. USEMONALB) THEN
 391                ALB   = ALBEDOMIN
 392             ENDIF
 393             Z0BRD = Z0MAX
 394          ELSE IF ( SHDFAC <= SHDMIN ) THEN
 395             EMBRD = EMISSMIN
 396             IF(.NOT. RDLAI2D) THEN
 397                XLAI  = LAIMIN
 398             ENDIF
 399             IF(.NOT. USEMONALB) then
 400                ALB   = ALBEDOMAX
 401             ENDIF
 402             Z0BRD = Z0MIN
 403          ELSE
 404
 405             IF ( SHDMAX > SHDMIN ) THEN
 406
 407                INTERP_FRACTION = ( SHDFAC - SHDMIN ) / ( SHDMAX - SHDMIN )
 408                ! Bound INTERP_FRACTION between 0 and 1
 409                INTERP_FRACTION = MIN ( INTERP_FRACTION, 1.0 )
 410                INTERP_FRACTION = MAX ( INTERP_FRACTION, 0.0 )
 411                ! Scale Emissivity and LAI between EMISSMIN and EMISSMAX by INTERP_FRACTION
 412                EMBRD = ( ( 1.0 - INTERP_FRACTION ) * EMISSMIN  ) + ( INTERP_FRACTION * EMISSMAX  )
 413                IF (.NOT. RDLAI2D) THEN
 414                   XLAI  = ( ( 1.0 - INTERP_FRACTION ) * LAIMIN    ) + ( INTERP_FRACTION * LAIMAX    )
 415                ENDIF
 416                if (.not. USEMONALB) then
 417                   ALB   = ( ( 1.0 - INTERP_FRACTION ) * ALBEDOMAX ) + ( INTERP_FRACTION * ALBEDOMIN )
 418                endif
 419                Z0BRD = ( ( 1.0 - INTERP_FRACTION ) * Z0MIN     ) + ( INTERP_FRACTION * Z0MAX     )
 420
 421             ELSE
 422
 423                EMBRD = 0.5 * EMISSMIN  + 0.5 * EMISSMAX
 424                IF (.NOT. RDLAI2D) THEN
 425                   XLAI  = 0.5 * LAIMIN    + 0.5 * LAIMAX
 426                ENDIF
 427                if (.not. USEMONALB) then
 428                   ALB   = 0.5 * ALBEDOMIN + 0.5 * ALBEDOMAX
 429                endif
 430                Z0BRD = 0.5 * Z0MIN     + 0.5 * Z0MAX
 431
 432             ENDIF
 433
 434          ENDIF
 435 ! ----------------------------------------------------------------------
 436 !  INITIALIZE PRECIPITATION LOGICALS.
 437 ! ----------------------------------------------------------------------
 438          SNOWNG = .FALSE.
 439          FRZGRA = .FALSE.
 440
 441 ! ----------------------------------------------------------------------
 442 ! OVER SEA-ICE OR GLACIAL-ICE, IF S.W.E. (SNEQV) BELOW THRESHOLD LOWER
 443 ! BOUND (0.01 M FOR SEA-ICE, 0.10 M FOR GLACIAL-ICE), THEN SET AT LOWER
 444 ! BOUND
 445 ! ----------------------------------------------------------------------
 446 ! IF SEA-ICE CASE, ASSIGN DEFAULT WATER-EQUIV SNOW ON TOP
 447 ! ----------------------------------------------------------------------
 448          IF (ICE == 1) THEN
 449             ! Sea-ice case
 450             IF ( SNEQV < 0.01 ) THEN
 451                SNEQV = 0.01
 452                SNOWH = 0.05
 453             ENDIF
 454          ELSE IF ( ICE == -1 ) THEN
 455             ! Land-ice case
 456             IF ( SNEQV < 0.10 ) THEN
 457                SNEQV = 0.10
 458                SNOWH = 0.50
 459             ENDIF
 460          END IF
 461 ! ----------------------------------------------------------------------
 462 ! FOR SEA-ICE AND GLACIAL-ICE CASES, SET SMC AND SH20 VALUES = 1.0
 463 ! ----------------------------------------------------------------------
 464          IF ( ICE /= 0 ) THEN
 465             DO KZ = 1,NSOIL
 466                SMC(KZ) = 1.0
 467                SH2O(KZ) = 1.0
 468             END DO
 469          ENDIF
 470 ! ----------------------------------------------------------------------
 471 ! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
 472 !   SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
 473 !   SUBROUTINE)
 474 ! ----------------------------------------------------------------------
 475          IF ( SNEQV <= 1.E-7 ) THEN ! safer IF  kmh (2008/03/25)
 476             SNEQV = 0.0
 477             SNDENS = 0.0
 478             SNOWH = 0.0
 479             SNCOND = 1.0
 480          ELSE
 481             SNDENS = SNEQV / SNOWH
 482             IF(SNDENS > 1.0) THEN
 483              CALL wrf_error_fatal ( 'Physical snow depth is less than snow water equiv.' )
 484             ENDIF
 485             CALL CSNOW (SNCOND,SNDENS)
 486          END IF
 487 ! ----------------------------------------------------------------------
 488 ! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
 489 ! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
 490 ! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
 491 ! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
 492 ! ----------------------------------------------------------------------
 493          IF (PRCP > 0.0) THEN
 494 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
 495 ! passed in from model microphysics.
 496             IF (FFROZP .GT. 0.5) THEN
 497                SNOWNG = .TRUE.
 498             ELSE
 499                IF (T1 <= TFREEZ) FRZGRA = .TRUE.
 500             END IF
 501          END IF
 502 ! ----------------------------------------------------------------------
 503 ! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
 504 ! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
 505 ! IT TO THE EXISTING SNOWPACK.
 506 ! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
 507 ! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
 508 ! ----------------------------------------------------------------------
 509          IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
 510             SN_NEW = PRCP * DT * 0.001
 511             SNEQV = SNEQV + SN_NEW
 512             PRCPF = 0.0
 513
 514 ! ----------------------------------------------------------------------
 515 ! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
 516 ! UPDATE SNOW THERMAL CONDUCTIVITY
 517 ! ----------------------------------------------------------------------
 518             CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
 519 !
 520 ! kmh 09/04/2006 set Snow Density at 0.2 g/cm**3
 521 ! for "cold permanent ice" or new "dry" snow
 522 !
 523            IF ( (ICE /= 0) .and. SNCOVR .GT. 0.99 ) THEN
 524 !  if soil temperature less than 268.15 K, treat as typical Antarctic/Greenland snow firn
 525               IF ( STC(1) .LT. (TFREEZ - 5.) ) SNDENS = 0.2
 526               IF ( SNOWNG .AND. (T1.LT.273.) .AND. (SFCTMP.LT.273.) ) SNDENS=0.2
 527            ENDIF
 528 !
 529             CALL CSNOW (SNCOND,SNDENS)
 530
 531 ! ----------------------------------------------------------------------
 532 ! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
 533 ! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
 534 ! ANY CANOPY "DRIP" ADDED TO THIS LATER)
 535 ! ----------------------------------------------------------------------
 536          ELSE
 537             PRCPF = PRCP
 538          END IF
 539 ! ----------------------------------------------------------------------
 540 ! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
 541 ! ----------------------------------------------------------------------
 542 ! ----------------------------------------------------------------------
 543 ! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
 544 ! ----------------------------------------------------------------------
 545          IF (ICE == 0 .OR. ICE == -1) THEN
 546             IF (SNEQV  == 0.0) THEN
 547                SNCOVR = 0.0
 548                ALBEDO = ALB
 549                EMISSI = EMBRD
 550             ELSE
 551 ! ----------------------------------------------------------------------
 552 ! DETERMINE SNOW FRACTIONAL COVERAGE.
 553 ! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
 554 ! ----------------------------------------------------------------------
 555               CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
 556 !   Don't limit snow cover fraction over permanent ice kmh 2008/03/25
 557             if ( ICE == 0 ) then
 558               SNCOVR = MIN(SNCOVR,0.98)
 559             endif
 560               CALL ALCALC (ALB,SNOALB,EMBRD,SHDFAC,SHDMIN,SNCOVR,T1,ALBEDO,EMISSI,   &
 561                          DT,SNOWNG,SNOTIME1,LVCOEF)
 562             END IF
 563 ! ----------------------------------------------------------------------
 564 ! SNOW COVER, ALBEDO OVER SEA-ICE, GLACIAL ICE
 565 ! ----------------------------------------------------------------------
 566          ELSE
 567             SNCOVR = 1.0
 568 !
 569 ! Albedo of sea ice
 570 !
 571 ! This value should vary seasonally. 0.65 may be good for Arctic Ocean summer bare ice
 572 !   value could be as low as 0.4 for Arctic bare ice and melt pond combo (Perovich data)
 573 !   0.82 may be good for Arctic spring/fall sea ice (Perovich data)
 574 !   0.81 may be good for Antarctic sea ice (Wendler et al. December cruise data)
 575 !
 576             ALBEDO = 0.80
 577 !
 578             EMISSI = 0.98
 579          END IF
 580 ! ----------------------------------------------------------------------
 581 ! THERMAL CONDUCTIVITY FOR SEA-ICE CASE, GLACIAL-ICE CASE
 582 ! ----------------------------------------------------------------------
 583          IF ( (ICE == 1) .or. (ICE == -1) ) THEN
 584             DF1 = 2.2
 585 !
 586 ! kmh 09/03/2006
 587 ! kmh 03/25/2008  change SNCOVR threshold to 0.97
 588 !
 589             IF (  SNCOVR .GT. 0.97 ) THEN
 590                DF1 = SNCOND
 591             ENDIF
 592 !
 593          ELSE
 594 ! ----------------------------------------------------------------------
 595 ! NEXT CALCULATE THE SUBSURFACE HEAT FLUX, WHICH FIRST REQUIRES
 596 ! CALCULATION OF THE THERMAL DIFFUSIVITY.  TREATMENT OF THE
 597 ! LATTER FOLLOWS THAT ON PAGES 148-149 FROM "HEAT TRANSFER IN
 598 ! COLD CLIMATES", BY V. J. LUNARDINI (PUBLISHED IN 1981
 599 ! BY VAN NOSTRAND REINHOLD CO.) I.E. TREATMENT OF TWO CONTIGUOUS
 600 ! "PLANE PARALLEL" MEDIUMS (NAMELY HERE THE FIRST SOIL LAYER
 601 ! AND THE SNOWPACK LAYER, IF ANY). THIS DIFFUSIVITY TREATMENT
 602 ! BEHAVES WELL FOR BOTH ZERO AND NONZERO SNOWPACK, INCLUDING THE
 603 ! LIMIT OF VERY THIN SNOWPACK.  THIS TREATMENT ALSO ELIMINATES
 604 ! THE NEED TO IMPOSE AN ARBITRARY UPPER BOUND ON SUBSURFACE
 605 ! HEAT FLUX WHEN THE SNOWPACK BECOMES EXTREMELY THIN.
 606 ! ----------------------------------------------------------------------
 607 ! FIRST CALCULATE THERMAL DIFFUSIVITY OF TOP SOIL LAYER, USING
 608 ! BOTH THE FROZEN AND LIQUID SOIL MOISTURE, FOLLOWING THE
 609 ! SOIL THERMAL DIFFUSIVITY FUNCTION OF PETERS-LIDARD ET AL.
 610 ! (1998,JAS, VOL 55, 1209-1224), WHICH REQUIRES THE SPECIFYING
 611 ! THE QUARTZ CONTENT OF THE GIVEN SOIL CLASS (SEE ROUTINE REDPRM)
 612 ! ----------------------------------------------------------------------
 613 ! ----------------------------------------------------------------------
 614 ! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
 615 ! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
 616 ! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
 617 ! ----------------------------------------------------------------------
 618             CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
 619
 620 !urban
 621             IF ( VEGTYP == ISURBAN ) DF1=3.24
 622
 623             DF1 = DF1 * EXP (SBETA * SHDFAC)
 624 !
 625 ! kmh 09/03/2006
 626 ! kmh 03/25/2008  change SNCOVR threshold to 0.97
 627 !
 628             IF (  SNCOVR .GT. 0.97 ) THEN
 629                DF1 = SNCOND
 630             ENDIF
 631 !
 632 ! ----------------------------------------------------------------------
 633 ! FINALLY "PLANE PARALLEL" SNOWPACK EFFECT FOLLOWING
 634 ! V.J. LINARDINI REFERENCE CITED ABOVE. NOTE THAT DTOT IS
 635 ! COMBINED DEPTH OF SNOWDEPTH AND THICKNESS OF FIRST SOIL LAYER
 636 ! ----------------------------------------------------------------------
 637          END IF
 638
 639          DSOIL = - (0.5 * ZSOIL (1))
 640          IF (SNEQV == 0.) THEN
 641             SSOIL = DF1 * (T1- STC (1) ) / DSOIL
 642          ELSE
 643             DTOT = SNOWH + DSOIL
 644             FRCSNO = SNOWH / DTOT
 645
 646 ! 1. HARMONIC MEAN (SERIES FLOW)
 647 !        DF1 = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
 648             FRCSOI = DSOIL / DTOT
 649 ! 2. ARITHMETIC MEAN (PARALLEL FLOW)
 650 !        DF1 = FRCSNO*SNCOND + FRCSOI*DF1
 651             DF1H = (SNCOND * DF1)/ (FRCSOI * SNCOND+ FRCSNO * DF1)
 652
 653 ! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
 654 !        DF1 = (SNCOND**FRCSNO)*(DF1**FRCSOI)
 655 ! weigh DF by snow fraction
 656 !        DF1 = DF1H*SNCOVR + DF1A*(1.0-SNCOVR)
 657 !        DF1 = DF1H*SNCOVR + DF1*(1.0-SNCOVR)
 658             DF1A = FRCSNO * SNCOND+ FRCSOI * DF1
 659
 660 ! ----------------------------------------------------------------------
 661 ! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
 662 ! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
 663 ! MID-LAYER SOIL TEMPERATURE
 664 ! ----------------------------------------------------------------------
 665             DF1 = DF1A * SNCOVR + DF1* (1.0- SNCOVR)
 666             IF ( ICE /= 0 ) then
 667                !  kmh  12/15/2005  correct for too deep snow layer
 668                !  kmh  09/03/2006  adjust DTOT
 669                IF ( DTOT .GT. 2.*DSOIL ) then
 670                   DTOT = 2.*DSOIL
 671                ENDIF
 672             ENDIF
 673             SSOIL = DF1 * (T1- STC (1) ) / DTOT
 674          END IF
 675 ! ----------------------------------------------------------------------
 676 ! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
 677 ! THE PREVIOUS TIMESTEP.
 678 ! ----------------------------------------------------------------------
 679          IF (SNCOVR  > 0. ) THEN
 680             CALL SNOWZ0 (SNCOVR,Z0,Z0BRD,SNOWH)
 681          ELSE
 682             Z0=Z0BRD
 683          END IF
 684 ! ----------------------------------------------------------------------
 685 ! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
 686 ! HEAT AND MOISTURE.
 687
 688 ! NOTE !!!
 689 ! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, IN CASE
 690 ! ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND ZILINTINKEVICH COEF
 691 ! (CZIL) ARE SET THERE VIA NAMELIST I/O.
 692
 693 ! NOTE !!!
 694 ! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD TIMES THE
 695 ! "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE.  HENCE THE CH
 696 ! RETURNED FROM SFCDIF HAS UNITS OF M/S.  THE IMPORTANT COMPANION
 697 ! COEFFICIENT OF CH, CARRIED HERE AS "RCH", IS THE CH FROM SFCDIF TIMES
 698 ! AIR DENSITY AND PARAMETER "CP".  "RCH" IS COMPUTED IN "CALL PENMAN".
 699 ! RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATER IN EQNS.
 700
 701 ! NOTE !!!
 702 ! ----------------------------------------------------------------------
 703 ! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
 704 ! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT. Needed as a state variable
 705 ! for iterative/implicit solution of CH in SFCDIF
 706 ! ----------------------------------------------------------------------
 707 !        IF(.NOT.LCH) THEN
 708 !          T1V = T1 * (1.0+ 0.61 * Q2)
 709 !          TH2V = TH2 * (1.0+ 0.61 * Q2)
 710 !          CALL SFCDIF_off (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
 711 !        ENDIF
 712
 713 ! ----------------------------------------------------------------------
 714 ! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
 715 ! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
 716 ! CALCULATIONS.
 717 ! ----------------------------------------------------------------------
 718 ! ----------------------------------------------------------------------
 719 ! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
 720 ! PENMAN EP SUBROUTINE THAT FOLLOWS
 721 ! ----------------------------------------------------------------------
 722 !         FDOWN = SOLDN * (1.0- ALBEDO) + LWDN
 723          FDOWN =  SOLNET + LWDN
 724 ! ----------------------------------------------------------------------
 725 ! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
 726 ! PENMAN.
 727          T2V = SFCTMP * (1.0+ 0.61 * Q2 )
 728
 729          iout=0
 730          if(iout.eq.1) then
 731          print*,'before penman'
 732          print*,' SFCTMP',SFCTMP,'SFCPRS',SFCPRS,'CH',CH,'T2V',T2V,      &
 733        'TH2',TH2,'PRCP',PRCP,'FDOWN',FDOWN,'T24',T24,'SSOIL',SSOIL,      &
 734         'Q2',Q2,'Q2SAT',Q2SAT,'ETP',ETP,'RCH',RCH,                       &
 735         'EPSCA',EPSCA,'RR',RR  ,'SNOWNG',SNOWNG,'FRZGRA',FRZGRA,           &
 736         'DQSDT2',DQSDT2,'FLX2',FLX2,'SNOWH',SNOWH,'SNEQV',SNEQV,         &
 737         ' DSOIL',DSOIL,' FRCSNO',FRCSNO,' SNCOVR',SNCOVR,' DTOT',DTOT,   &
 738        ' ZSOIL (1)',ZSOIL(1),' DF1',DF1,'T1',T1,' STC1',STC(1),          &
 739         'ALBEDO',ALBEDO,'SMC',SMC,'STC',STC,'SH2O',SH2O
 740          endif
 741
 742          CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL,     &
 743                        Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA,          &
 744 !
 745 ! kmh 01/09/2007 add T1,ICE,SNCOVR to call
 746 !
 747                          DQSDT2,FLX2,EMISSI,SNEQV,T1,ICE,SNCOVR)
 748 !
 749 ! ----------------------------------------------------------------------
 750 ! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
 751 ! IF NONZERO GREENNESS FRACTION
 752 ! ----------------------------------------------------------------------
 753
 754 ! ----------------------------------------------------------------------
 755 !  FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
 756 !  BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
 757 ! ----------------------------------------------------------------------
 758          IF (SHDFAC > 0.) THEN
 759             CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL,     &
 760                           SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2,  &
 761                           TOPT,RSMAX,RGL,HS,XLAI,                        &
 762                           RCS,RCT,RCQ,RCSOIL,EMISSI)
 763          END IF
 764 ! ----------------------------------------------------------------------
 765 ! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
 766 ! EXISTS OR NOT:
 767 ! ----------------------------------------------------------------------
 768          ESNOW = 0.0
 769          IF (SNEQV  == 0.0) THEN
 770             CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT,                  &
 771                             SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,           &
 772                             SHDFAC,                                      &
 773                             SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI,  &
 774                             SSOIL,                                       &
 775                             STC,EPSCA,BEXP,PC,RCH,RR,CFACTR,             &
 776                             SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL,            &
 777                             DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2,       &
 778                             RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS,      &
 779                             QUARTZ,FXEXP,CSOIL,                          &
 780                             BETA,DRIP,DEW,FLX1,FLX3,VEGTYP,ISURBAN)
 781             ETA_KINEMATIC = ETA
 782          ELSE
 783             CALL SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT,    &
 784                          SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,              &
 785                          SBETA,DF1,                                      &
 786                          Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA,  &
 787                          SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS,&
 788                          SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT,               &
 789                          ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1,     &
 790                          RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT,    &
 791                          ICE,RTDIS,QUARTZ,FXEXP,CSOIL,                   &
 792                          BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI, &
 793                          RIBB,SOLDN,                                     &
 794                          ISURBAN,                                        &
 795                          VEGTYP)
 796             ETA_KINEMATIC =  ESNOW + ETNS
 797          END IF
 798
 799 !     Calculate effective mixing ratio at grnd level (skin)
 800 !
 801 !     Q1=Q2+ETA*CP/RCH
 802      Q1=Q2+ETA_KINEMATIC*CP/RCH
 803 !
 804 ! ----------------------------------------------------------------------
 805 ! DETERMINE SENSIBLE HEAT (H) IN ENERGY UNITS (W M-2)
 806 ! ----------------------------------------------------------------------
 807          SHEAT = - (CH * CP * SFCPRS)/ (R * T2V) * ( TH2- T1 )
 808
 809 ! ----------------------------------------------------------------------
 810 ! CONVERT EVAP TERMS FROM KINEMATIC (KG M-2 S-1) TO ENERGY UNITS (W M-2)
 811 ! ----------------------------------------------------------------------
 812       EDIR = EDIR * LVH2O
 813       EC = EC * LVH2O
 814       DO K=1,4
 815       ET(K) = ET(K) * LVH2O
 816       ENDDO
 817       ETT = ETT * LVH2O
 818       ESNOW = ESNOW * LSUBS
 819       ETP = ETP*((1.-SNCOVR)*LVH2O + SNCOVR*LSUBS)
 820       IF (ETP .GT. 0.) THEN
 821          ETA = EDIR + EC + ETT + ESNOW
 822       ELSE
 823         ETA = ETP
 824       ENDIF
 825 ! ----------------------------------------------------------------------
 826 ! DETERMINE BETA (RATIO OF ACTUAL TO POTENTIAL EVAP)
 827 ! ----------------------------------------------------------------------
 828       IF (ETP == 0.0) THEN
 829         BETA = 0.0
 830       ELSE
 831         BETA = ETA/ETP
 832       ENDIF
 833
 834 ! ----------------------------------------------------------------------
 835 ! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
 836 !   SSOIL>0: WARM THE SURFACE  (NIGHT TIME)
 837 !   SSOIL<0: COOL THE SURFACE  (DAY TIME)
 838 ! ----------------------------------------------------------------------
 839          SSOIL = -1.0* SSOIL
 840
 841 ! ----------------------------------------------------------------------
 842 !  FOR THE CASE OF LAND (BUT NOT GLACIAL-ICE):
 843 !  CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
 844 !  AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW.  RUNOFF2 IS ALREADY
 845 !  A RATE AT THIS POINT
 846 ! ----------------------------------------------------------------------
 847          IF (ICE == 0) THEN
 848             RUNOFF3 = RUNOFF3/ DT
 849             RUNOFF2 = RUNOFF2+ RUNOFF3
 850             SOILM = -1.0* SMC (1)* ZSOIL (1)
 851             DO K = 2,NSOIL
 852               SOILM = SOILM + SMC (K)* (ZSOIL (K -1) - ZSOIL (K))
 853             END DO
 854             SOILWM = -1.0* (SMCMAX - SMCWLT)* ZSOIL (1)
 855             SOILWW = -1.0* (SMC (1) - SMCWLT)* ZSOIL (1)
 856
 857             IF (NROOT >= 2) THEN
 858               DO K = 2,NROOT
 859                SOILWM = SOILWM + (SMCMAX - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
 860                SOILWW = SOILWW + (SMC(K) - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
 861               END DO
 862             END IF
 863             IF (SOILWM .LT. 1.E-6) THEN
 864               SOILWM = 0.0
 865               SOILW  = 0.0
 866               SOILM  = 0.0
 867             ELSE
 868               SOILW = SOILWW / SOILWM
 869             END IF
 870          ELSE
 871 ! ----------------------------------------------------------------------
 872 ! FOR THE CASE OF SEA-ICE (ICE=1) OR GLACIAL-ICE (ICE=-1), ADD ANY
 873 ! SNOWMELT DIRECTLY TO SURFACE RUNOFF (RUNOFF1) SINCE THERE IS NO
 874 ! SOIL MEDIUM, AND THUS NO CALL TO SUBROUTINE SMFLX (FOR SOIL MOISTURE
 875 ! TENDENCY).
 876 ! ----------------------------------------------------------------------
 877             RUNOFF1 = SNOMLT/DT
 878            SOILWM = 0.0
 879            SOILW  = 0.0
 880            SOILM  = 0.0
 881          END IF
 882
 883 ! ----------------------------------------------------------------------
 884   END SUBROUTINE SFLX
 885 ! ----------------------------------------------------------------------
 886
 887       SUBROUTINE ALCALC (ALB,SNOALB,EMBRD,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO,EMISSI,   &
 888                          DT,SNOWNG,SNOTIME1,LVCOEF)
 889
 890 ! ----------------------------------------------------------------------
 891 ! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
 892 !   ALB     SNOWFREE ALBEDO
 893 !   SNOALB  MAXIMUM (DEEP) SNOW ALBEDO
 894 !   SHDFAC    AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
 895 !   SHDMIN    MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
 896 !   SNCOVR  FRACTIONAL SNOW COVER
 897 !   ALBEDO  SURFACE ALBEDO INCLUDING SNOW EFFECT
 898 !   TSNOW   SNOW SURFACE TEMPERATURE (K)
 899 ! ----------------------------------------------------------------------
 900       IMPLICIT NONE
 901
 902 ! ----------------------------------------------------------------------
 903 ! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
 904 ! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
 905 ! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
 906 ! (1985, JCAM, VOL 24, 402-411)
 907 ! ----------------------------------------------------------------------
 908       REAL, INTENT(IN)  ::  ALB, SNOALB, EMBRD, SHDFAC, SHDMIN, SNCOVR, TSNOW
 909       REAL, INTENT(IN)  :: DT
 910       LOGICAL, INTENT(IN) :: SNOWNG
 911       REAL, INTENT(INOUT):: SNOTIME1
 912       REAL, INTENT(OUT) ::  ALBEDO, EMISSI
 913       REAL              :: SNOALB2
 914       REAL              :: TM,SNOALB1
 915       REAL, INTENT(IN)  :: LVCOEF
 916       REAL, PARAMETER   :: SNACCA=0.94,SNACCB=0.58,SNTHWA=0.82,SNTHWB=0.46
 917 ! turn of vegetation effect
 918 !      ALBEDO = ALB + (1.0- (SHDFAC - SHDMIN))* SNCOVR * (SNOALB - ALB)
 919 !      ALBEDO = (1.0-SNCOVR)*ALB + SNCOVR*SNOALB !this is equivalent to below
 920       ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
 921       EMISSI = EMBRD + SNCOVR*(EMISSI_S - EMBRD)
 922
 923 !     BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
 924 !          IF (TSNOW.LE.263.16) THEN
 925 !            ALBEDO=SNOALB
 926 !          ELSE
 927 !            IF (TSNOW.LT.273.16) THEN
 928 !              TM=0.1*(TSNOW-263.16)
 929 !              SNOALB1=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
 930 !            ELSE
 931 !              SNOALB1=0.67
 932 !             IF(SNCOVR.GT.0.95) SNOALB1= 0.6
 933 !             SNOALB1 = ALB + SNCOVR*(SNOALB-ALB)
 934 !            ENDIF
 935 !          ENDIF
 936 !            ALBEDO = ALB + SNCOVR*(SNOALB1-ALB)
 937
 938 !     ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
 939 !          SNOALB1 = SNOALB+COEF*(0.85-SNOALB)
 940 !          SNOALB2=SNOALB1
 941 !!m          LSTSNW=LSTSNW+1
 942 !          SNOTIME1 = SNOTIME1 + DT
 943 !          IF (SNOWNG) THEN
 944 !             SNOALB2=SNOALB
 945 !!m             LSTSNW=0
 946 !             SNOTIME1 = 0.0
 947 !          ELSE
 948 !            IF (TSNOW.LT.273.16) THEN
 949 !!              SNOALB2=SNOALB-0.008*LSTSNW*DT/86400
 950 !!m              SNOALB2=SNOALB-0.008*SNOTIME1/86400
 951 !              SNOALB2=(SNOALB2-0.65)*EXP(-0.05*DT/3600)+0.65
 952 !!              SNOALB2=(ALBEDO-0.65)*EXP(-0.01*DT/3600)+0.65
 953 !            ELSE
 954 !              SNOALB2=(SNOALB2-0.5)*EXP(-0.0005*DT/3600)+0.5
 955 !!              SNOALB2=(SNOALB-0.5)*EXP(-0.24*LSTSNW*DT/86400)+0.5
 956 !!m              SNOALB2=(SNOALB-0.5)*EXP(-0.24*SNOTIME1/86400)+0.5
 957 !            ENDIF
 958 !          ENDIF
 959 !
 960 !!               print*,'SNOALB2',SNOALB2,'ALBEDO',ALBEDO,'DT',DT
 961 !          ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
 962 !          IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
 963 !!m          LSTSNW1=LSTSNW
 964 !!          SNOTIME = SNOTIME1
 965
 966 ! formulation by Livneh
 967 ! ----------------------------------------------------------------------
 968 ! SNOALB IS CONSIDERED AS THE MAXIMUM SNOW ALBEDO FOR NEW SNOW, AT
 969 ! A VALUE OF 85%. SNOW ALBEDO CURVE DEFAULTS ARE FROM BRAS P.263. SHOULD
 970 ! NOT BE CHANGED EXCEPT FOR SERIOUS PROBLEMS WITH SNOW MELT.
 971 ! TO IMPLEMENT ACCUMULATIN PARAMETERS, SNACCA AND SNACCB, ASSERT THAT IT
 972 ! IS INDEED ACCUMULATION SEASON. I.E. THAT SNOW SURFACE TEMP IS BELOW
 973 ! ZERO AND THE DATE FALLS BETWEEN OCTOBER AND FEBRUARY
 974 ! ----------------------------------------------------------------------
 975          SNOALB1 = SNOALB+LVCOEF*(0.85-SNOALB)
 976          SNOALB2=SNOALB1
 977 ! ---------------- Initial LSTSNW --------------------------------------
 978           IF (SNOWNG) THEN
 979              SNOTIME1 = 0.
 980           ELSE
 981            SNOTIME1=SNOTIME1+DT
 982 !               IF (TSNOW.LT.273.16) THEN
 983                    SNOALB2=SNOALB1*(SNACCA**((SNOTIME1/86400.0)**SNACCB))
 984 !               ELSE
 985 !                  SNOALB2 =SNOALB1*(SNTHWA**((SNOTIME1/86400.0)**SNTHWB))
 986 !               ENDIF
 987           ENDIF
 988 !
 989            SNOALB2 = MAX ( SNOALB2, ALB )
 990            ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
 991            IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
 992
 993 !          IF (TSNOW.LT.273.16) THEN
 994 !            ALBEDO=SNOALB-0.008*DT/86400
 995 !          ELSE
 996 !            ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
 997 !          ENDIF
 998
 999 !      IF (ALBEDO > SNOALB) ALBEDO = SNOALB
1000
1001 ! ----------------------------------------------------------------------
1002   END SUBROUTINE ALCALC
1003 ! ----------------------------------------------------------------------
1004
1005       SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL,       &
1006                          SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2,    &
1007                          TOPT,RSMAX,RGL,HS,XLAI,                          &
1008                          RCS,RCT,RCQ,RCSOIL,EMISSI)
1009
1010 ! ----------------------------------------------------------------------
1011 ! SUBROUTINE CANRES
1012 ! ----------------------------------------------------------------------
1013 ! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
1014 ! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
1015 ! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
1016 ! MOISTURE RATHER THAN TOTAL)
1017 ! ----------------------------------------------------------------------
1018 ! SOURCE:  JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
1019 ! NOILHAN (1990, BLM)
1020 ! SEE ALSO:  CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
1021 ! AND TABLE 2 OF SEC. 3.1.2
1022 ! ----------------------------------------------------------------------
1023 ! INPUT:
1024 !   SOLAR   INCOMING SOLAR RADIATION
1025 !   CH      SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1026 !   SFCTMP  AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
1027 !   Q2      AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1028 !   Q2SAT   SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1029 !   DQSDT2  SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
1030 !   SFCPRS  SURFACE PRESSURE
1031 !   SMC     VOLUMETRIC SOIL MOISTURE
1032 !   ZSOIL   SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
1033 !   NSOIL   NO. OF SOIL LAYERS
1034 !   NROOT   NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
1035 !   XLAI    LEAF AREA INDEX
1036 !   SMCWLT  WILTING POINT
1037 !   SMCREF  REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
1038 !             SETS IN)
1039 ! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
1040 !   SURBOUTINE REDPRM
1041 ! OUTPUT:
1042 !   PC  PLANT COEFFICIENT
1043 !   RC  CANOPY RESISTANCE
1044 ! ----------------------------------------------------------------------
1045
1046       IMPLICIT NONE
1047       INTEGER, INTENT(IN) :: NROOT,NSOIL
1048       INTEGER  K
1049       REAL,    INTENT(IN) :: CH,DQSDT2,HS,Q2,Q2SAT,RSMIN,RGL,RSMAX,        &
1050                              SFCPRS,SFCTMP,SMCREF,SMCWLT, SOLAR,TOPT,XLAI, &
1051                              EMISSI
1052       REAL,DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
1053       REAL,    INTENT(OUT):: PC,RC,RCQ,RCS,RCSOIL,RCT
1054       REAL                :: DELTA,FF,GX,P,RR
1055       REAL, DIMENSION(1:NSOIL) ::  PART
1056       REAL, PARAMETER     :: SLV = 2.501000E6
1057
1058
1059 ! ----------------------------------------------------------------------
1060 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
1061 ! ----------------------------------------------------------------------
1062       RCS = 0.0
1063       RCT = 0.0
1064       RCQ = 0.0
1065       RCSOIL = 0.0
1066
1067 ! ----------------------------------------------------------------------
1068 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
1069 ! ----------------------------------------------------------------------
1070       RC = 0.0
1071       FF = 0.55*2.0* SOLAR / (RGL * XLAI)
1072       RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)
1073
1074 ! ----------------------------------------------------------------------
1075 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
1076 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
1077 ! ----------------------------------------------------------------------
1078       RCS = MAX (RCS,0.0001)
1079       RCT = 1.0- 0.0016* ( (TOPT - SFCTMP)**2.0)
1080
1081 ! ----------------------------------------------------------------------
1082 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
1083 ! RCQ EXPRESSION FROM SSIB
1084 ! ----------------------------------------------------------------------
1085       RCT = MAX (RCT,0.0001)
1086       RCQ = 1.0/ (1.0+ HS * (Q2SAT - Q2))
1087
1088 ! ----------------------------------------------------------------------
1089 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
1090 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
1091 ! ----------------------------------------------------------------------
1092       RCQ = MAX (RCQ,0.01)
1093       GX = (SMC (1) - SMCWLT) / (SMCREF - SMCWLT)
1094       IF (GX  >  1.) GX = 1.
1095       IF (GX  <  0.) GX = 0.
1096
1097 ! ----------------------------------------------------------------------
1098 ! USE SOIL DEPTH AS WEIGHTING FACTOR
1099 ! ----------------------------------------------------------------------
1100 ! ----------------------------------------------------------------------
1101 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1102 !      PART(1) = RTDIS(1) * GX
1103 ! ----------------------------------------------------------------------
1104       PART (1) = (ZSOIL (1)/ ZSOIL (NROOT)) * GX
1105       DO K = 2,NROOT
1106          GX = (SMC (K) - SMCWLT) / (SMCREF - SMCWLT)
1107          IF (GX >  1.) GX = 1.
1108          IF (GX <  0.) GX = 0.
1109 ! ----------------------------------------------------------------------
1110 ! USE SOIL DEPTH AS WEIGHTING FACTOR
1111 ! ----------------------------------------------------------------------
1112 ! ----------------------------------------------------------------------
1113 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1114 !        PART(K) = RTDIS(K) * GX
1115 ! ----------------------------------------------------------------------
1116          PART (K) = ( (ZSOIL (K) - ZSOIL (K -1))/ ZSOIL (NROOT)) * GX
1117       END DO
1118       DO K = 1,NROOT
1119          RCSOIL = RCSOIL + PART (K)
1120       END DO
1121
1122 ! ----------------------------------------------------------------------
1123 ! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS.  CONVERT CANOPY
1124 ! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
1125 ! EVAP IN DETERMINING ACTUAL EVAP.  PC IS DETERMINED BY:
1126 !   PC * LINERIZED PENMAN POTENTIAL EVAP =
1127 !   PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
1128 ! ----------------------------------------------------------------------
1129       RCSOIL = MAX (RCSOIL,0.0001)
1130
1131       RC = RSMIN / (XLAI * RCS * RCT * RCQ * RCSOIL)
1132 !      RR = (4.* SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) + 1.0
1133       RR = (4.* EMISSI *SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) &
1134              + 1.0
1135
1136       DELTA = (SLV / CP)* DQSDT2
1137
1138       PC = (RR + DELTA)/ (RR * (1. + RC * CH) + DELTA)
1139
1140 ! ----------------------------------------------------------------------
1141   END SUBROUTINE CANRES
1142 ! ----------------------------------------------------------------------
1143
1144       SUBROUTINE CSNOW (SNCOND,DSNOW)
1145
1146 ! ----------------------------------------------------------------------
1147 ! SUBROUTINE CSNOW
1148 ! FUNCTION CSNOW
1149 ! ----------------------------------------------------------------------
1150 ! CALCULATE SNOW TERMAL CONDUCTIVITY
1151 ! ----------------------------------------------------------------------
1152       IMPLICIT NONE
1153       REAL, INTENT(IN) :: DSNOW
1154       REAL, INTENT(OUT):: SNCOND
1155       REAL             :: C
1156       REAL, PARAMETER  :: UNIT = 0.11631
1157
1158 ! ----------------------------------------------------------------------
1159 ! SNCOND IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
1160 ! CSNOW IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
1161 ! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
1162 ! ----------------------------------------------------------------------
1163       C = 0.328*10** (2.25* DSNOW)
1164 !      CSNOW=UNIT*C
1165
1166 ! ----------------------------------------------------------------------
1167 ! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
1168 ! ----------------------------------------------------------------------
1169 !      SNCOND=0.0293*(1.+100.*DSNOW**2)
1170 !      CSNOW=0.0293*(1.+100.*DSNOW**2)
1171
1172 ! ----------------------------------------------------------------------
1173 ! E. ANDERSEN FROM FLERCHINGER
1174 ! ----------------------------------------------------------------------
1175 !      SNCOND=0.021+2.51*DSNOW**2
1176 !      CSNOW=0.021+2.51*DSNOW**2
1177
1178 !      SNCOND = UNIT * C
1179 ! double snow thermal conductivity
1180       SNCOND = 2.0 * UNIT * C
1181
1182 ! ----------------------------------------------------------------------
1183   END SUBROUTINE CSNOW
1184 ! ----------------------------------------------------------------------
1185
1186       SUBROUTINE DEVAP (EDIR,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP,         &
1187                         DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1188
1189 ! ----------------------------------------------------------------------
1190 ! SUBROUTINE DEVAP
1191 ! FUNCTION DEVAP
1192 ! ----------------------------------------------------------------------
1193 ! CALCULATE DIRECT SOIL EVAPORATION
1194 ! ----------------------------------------------------------------------
1195       IMPLICIT NONE
1196       REAL, INTENT(IN) :: ETP1,SMC,BEXP,DKSAT,DWSAT,FXEXP,              &
1197                           SHDFAC,SMCDRY,SMCMAX,ZSOIL,SMCREF,SMCWLT
1198       REAL, INTENT(OUT):: EDIR
1199       REAL             :: FX, SRATIO
1200
1201
1202 ! ----------------------------------------------------------------------
1203 ! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
1204 ! WHEN FXEXP=1.
1205 ! ----------------------------------------------------------------------
1206 ! ----------------------------------------------------------------------
1207 ! FX > 1 REPRESENTS DEMAND CONTROL
1208 ! FX < 1 REPRESENTS FLUX CONTROL
1209 ! ----------------------------------------------------------------------
1210
1211       SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
1212       IF (SRATIO > 0.) THEN
1213         FX = SRATIO**FXEXP
1214         FX = MAX ( MIN ( FX, 1. ) ,0. )
1215       ELSE
1216         FX = 0.
1217       ENDIF
1218
1219 ! ----------------------------------------------------------------------
1220 ! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
1221 ! ----------------------------------------------------------------------
1222       EDIR = FX * ( 1.0- SHDFAC ) * ETP1
1223
1224 ! ----------------------------------------------------------------------
1225   END SUBROUTINE DEVAP
1226 ! ----------------------------------------------------------------------
1227
1228       SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,               &
1229                          SH2O,                                          &
1230                          SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,             &
1231                          SMCREF,SHDFAC,CMCMAX,                          &
1232                          SMCDRY,CFACTR,                                 &
1233                          EDIR,EC,ET,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
1234
1235 ! ----------------------------------------------------------------------
1236 ! SUBROUTINE EVAPO
1237 ! ----------------------------------------------------------------------
1238 ! CALCULATE SOIL MOISTURE FLUX.  THE SOIL MOISTURE CONTENT (SMC - A PER
1239 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
1240 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
1241 ! FROZEN GROUND VERSION:  NEW STATES ADDED: SH2O, AND FROZEN GROUND
1242 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
1243 ! ----------------------------------------------------------------------
1244       IMPLICIT NONE
1245       INTEGER, INTENT(IN)   :: NSOIL, NROOT
1246       INTEGER               :: I,K
1247       REAL,    INTENT(IN)   :: BEXP, CFACTR,CMC,CMCMAX,DKSAT,           &
1248                                  DT,DWSAT,ETP1,FXEXP,PC,Q2,SFCTMP,      &
1249                                  SHDFAC,SMCDRY,SMCMAX,SMCREF,SMCWLT
1250       REAL,    INTENT(OUT)  :: EC,EDIR,ETA1,ETT
1251       REAL                  :: CMC2MS
1252       REAL,DIMENSION(1:NSOIL), INTENT(IN) :: RTDIS, SMC, SH2O, ZSOIL
1253       REAL,DIMENSION(1:NSOIL), INTENT(OUT) :: ET
1254
1255 ! ----------------------------------------------------------------------
1256 ! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
1257 ! GREATER THAN ZERO.
1258 ! ----------------------------------------------------------------------
1259       EDIR = 0.
1260       EC = 0.
1261       ETT = 0.
1262       DO K = 1,NSOIL
1263          ET (K) = 0.
1264       END DO
1265
1266 ! ----------------------------------------------------------------------
1267 ! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE.  CALL THIS FUNCTION
1268 ! ONLY IF VEG COVER NOT COMPLETE.
1269 ! FROZEN GROUND VERSION:  SH2O STATES REPLACE SMC STATES.
1270 ! ----------------------------------------------------------------------
1271       IF (ETP1 > 0.0) THEN
1272          IF (SHDFAC <  1.) THEN
1273              CALL DEVAP (EDIR,ETP1,SMC (1),ZSOIL (1),SHDFAC,SMCMAX,      &
1274                          BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1275          END IF
1276 ! ----------------------------------------------------------------------
1277 ! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
1278 ! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
1279 ! ----------------------------------------------------------------------
1280
1281          IF (SHDFAC > 0.0) THEN
1282             CALL TRANSP (ET,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT,     &
1283                           CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
1284             DO K = 1,NSOIL
1285                ETT = ETT + ET ( K )
1286             END DO
1287 ! ----------------------------------------------------------------------
1288 ! CALCULATE CANOPY EVAPORATION.
1289 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
1290 ! ----------------------------------------------------------------------
1291             IF (CMC > 0.0) THEN
1292                EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
1293             ELSE
1294                EC = 0.0
1295             END IF
1296 ! ----------------------------------------------------------------------
1297 ! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
1298 ! CANOPY.  -F.CHEN, 18-OCT-1994
1299 ! ----------------------------------------------------------------------
1300             CMC2MS = CMC / DT
1301             EC = MIN ( CMC2MS, EC )
1302          END IF
1303       END IF
1304 ! ----------------------------------------------------------------------
1305 ! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
1306 ! ----------------------------------------------------------------------
1307       ETA1 = EDIR + ETT + EC
1308
1309 ! ----------------------------------------------------------------------
1310   END SUBROUTINE EVAPO
1311 ! ----------------------------------------------------------------------
1312
1313   SUBROUTINE FAC2MIT(SMCMAX,FLIMIT)
1314     IMPLICIT NONE
1315     REAL, INTENT(IN)  :: SMCMAX
1316     REAL, INTENT(OUT) :: FLIMIT
1317
1318     FLIMIT = 0.90
1319
1320     IF ( SMCMAX == 0.395 ) THEN
1321        FLIMIT = 0.59
1322     ELSE IF ( ( SMCMAX == 0.434 ) .OR. ( SMCMAX == 0.404 ) ) THEN
1323        FLIMIT = 0.85
1324     ELSE IF ( ( SMCMAX == 0.465 ) .OR. ( SMCMAX == 0.406 ) ) THEN
1325        FLIMIT = 0.86
1326     ELSE IF ( ( SMCMAX == 0.476 ) .OR. ( SMCMAX == 0.439 ) ) THEN
1327        FLIMIT = 0.74
1328     ELSE IF ( ( SMCMAX == 0.200 ) .OR. ( SMCMAX == 0.464 ) ) THEN
1329        FLIMIT = 0.80
1330     ENDIF
1331
1332 ! ----------------------------------------------------------------------
1333   END SUBROUTINE FAC2MIT
1334 ! ----------------------------------------------------------------------
1335
1336       SUBROUTINE FRH2O (FREE,TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
1337
1338 ! ----------------------------------------------------------------------
1339 ! SUBROUTINE FRH2O
1340 ! ----------------------------------------------------------------------
1341 ! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
1342 ! TEMPERATURE IS BELOW 273.15K (T0).  REQUIRES NEWTON-TYPE ITERATION TO
1343 ! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
1344 ! (1999, JGR, VOL 104(D16), 19569-19585).
1345 ! ----------------------------------------------------------------------
1346 ! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
1347 ! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
1348 ! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE.  ALSO, EXPLICIT
1349 ! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
1350 ! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
1351 ! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
1352 ! LIMIT OF FREEZING POINT TEMPERATURE T0.
1353 ! ----------------------------------------------------------------------
1354 ! INPUT:
1355
1356 !   TKELV.........TEMPERATURE (Kelvin)
1357 !   SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
1358 !   SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
1359 !   SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
1360 !   B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
1361 !   PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
1362
1363 ! OUTPUT:
1364 !   FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
1365 !   FREE..........SUPERCOOLED LIQUID WATER CONTENT
1366 ! ----------------------------------------------------------------------
1367       IMPLICIT NONE
1368       REAL, INTENT(IN)     :: BEXP,PSIS,SH2O,SMC,SMCMAX,TKELV
1369       REAL, INTENT(OUT)    :: FREE
1370       REAL                 :: BX,DENOM,DF,DSWL,FK,SWL,SWLK
1371       INTEGER              :: NLOG,KCOUNT
1372 !      PARAMETER(CK = 0.0)
1373       REAL, PARAMETER      :: CK = 8.0, BLIM = 5.5, ERROR = 0.005,       &
1374                               HLICE = 3.335E5, GS = 9.81,DICE = 920.0,   &
1375                               DH2O = 1000.0, T0 = 273.15
1376
1377 ! ----------------------------------------------------------------------
1378 ! LIMITS ON PARAMETER B: B < 5.5  (use parameter BLIM)
1379 ! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
1380 ! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
1381 ! ----------------------------------------------------------------------
1382       BX = BEXP
1383
1384 ! ----------------------------------------------------------------------
1385 ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
1386 ! ----------------------------------------------------------------------
1387       IF (BEXP >  BLIM) BX = BLIM
1388       NLOG = 0
1389
1390 ! ----------------------------------------------------------------------
1391 !  IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
1392 ! ----------------------------------------------------------------------
1393       KCOUNT = 0
1394 !      FRH2O = SMC
1395       IF (TKELV > (T0- 1.E-3)) THEN
1396           FREE = SMC
1397       ELSE
1398
1399 ! ----------------------------------------------------------------------
1400 ! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
1401 ! IN KOREN ET AL, JGR, 1999, EQN 17
1402 ! ----------------------------------------------------------------------
1403 ! INITIAL GUESS FOR SWL (frozen content)
1404 ! ----------------------------------------------------------------------
1405          IF (CK /= 0.0) THEN
1406             SWL = SMC - SH2O
1407 ! ----------------------------------------------------------------------
1408 ! KEEP WITHIN BOUNDS.
1409 ! ----------------------------------------------------------------------
1410             IF (SWL > (SMC -0.02)) SWL = SMC -0.02
1411
1412 ! ----------------------------------------------------------------------
1413 !  START OF ITERATIONS
1414 ! ----------------------------------------------------------------------
1415             IF (SWL < 0.) SWL = 0.
1416  1001       Continue
1417               IF (.NOT.( (NLOG < 10) .AND. (KCOUNT == 0)))   goto 1002
1418               NLOG = NLOG +1
1419               DF = ALOG ( ( PSIS * GS / HLICE ) * ( ( 1. + CK * SWL )**2.) * &
1420                    ( SMCMAX / (SMC - SWL) )** BX) - ALOG ( - (               &
1421                    TKELV - T0)/ TKELV)
1422               DENOM = 2. * CK / ( 1. + CK * SWL ) + BX / ( SMC - SWL )
1423               SWLK = SWL - DF / DENOM
1424 ! ----------------------------------------------------------------------
1425 ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
1426 ! ----------------------------------------------------------------------
1427               IF (SWLK > (SMC -0.02)) SWLK = SMC - 0.02
1428               IF (SWLK < 0.) SWLK = 0.
1429
1430 ! ----------------------------------------------------------------------
1431 ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
1432 ! ----------------------------------------------------------------------
1433               DSWL = ABS (SWLK - SWL)
1434
1435 ! ----------------------------------------------------------------------
1436 ! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
1437 ! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
1438 ! ----------------------------------------------------------------------
1439               SWL = SWLK
1440               IF ( DSWL <= ERROR ) THEN
1441                     KCOUNT = KCOUNT +1
1442               END IF
1443 ! ----------------------------------------------------------------------
1444 !  END OF ITERATIONS
1445 ! ----------------------------------------------------------------------
1446 ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
1447 ! ----------------------------------------------------------------------
1448 !          FRH2O = SMC - SWL
1449            goto 1001
1450  1002   continue
1451            FREE = SMC - SWL
1452          END IF
1453 ! ----------------------------------------------------------------------
1454 ! END OPTION 1
1455 ! ----------------------------------------------------------------------
1456 ! ----------------------------------------------------------------------
1457 ! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
1458 ! IN KOREN ET AL., JGR, 1999, EQN 17
1459 ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
1460 ! ----------------------------------------------------------------------
1461          IF (KCOUNT == 0) THEN
1462              PRINT *,'Flerchinger USEd in NEW version. Iterations=',NLOG
1463                   FK = ( ( (HLICE / (GS * ( - PSIS)))*                    &
1464                        ( (TKELV - T0)/ TKELV))** ( -1/ BX))* SMCMAX
1465 !            FRH2O = MIN (FK, SMC)
1466              IF (FK < 0.02) FK = 0.02
1467              FREE = MIN (FK, SMC)
1468 ! ----------------------------------------------------------------------
1469 ! END OPTION 2
1470 ! ----------------------------------------------------------------------
1471          END IF
1472       END IF
1473 ! ----------------------------------------------------------------------
1474   END SUBROUTINE FRH2O
1475 ! ----------------------------------------------------------------------
1476
1477       SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,          &
1478                        TBOT,ZBOT,PSISAT,SH2O,DT,BEXP,                   &
1479                        F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP,ISURBAN)
1480
1481 ! ----------------------------------------------------------------------
1482 ! SUBROUTINE HRT
1483 ! ----------------------------------------------------------------------
1484 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1485 ! THERMAL DIFFUSION EQUATION.  ALSO TO COMPUTE ( PREPARE ) THE MATRIX
1486 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1487 ! ----------------------------------------------------------------------
1488       IMPLICIT NONE
1489       LOGICAL              :: ITAVG
1490       INTEGER, INTENT(IN)  :: NSOIL, VEGTYP
1491       INTEGER, INTENT(IN)  :: ISURBAN
1492       INTEGER              :: I, K
1493
1494       REAL, INTENT(IN)     :: BEXP, CSOIL, DF1, DT,F1,PSISAT,QUARTZ,     &
1495                               SMCMAX ,TBOT,YY,ZZ1, ZBOT
1496       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: SMC,STC,ZSOIL
1497       REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SH2O
1498       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: RHSTS
1499       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: AI, BI,CI
1500       REAL                 :: DDZ, DDZ2, DENOM, DF1N, DF1K, DTSDZ,       &
1501                               DTSDZ2,HCPCT,QTOT,SSOIL,SICE,TAVG,TBK,     &
1502                               TBK1,TSNSR,TSURF,CSOIL_LOC
1503       REAL, PARAMETER      :: T0 = 273.15, CAIR = 1004.0, CICE = 2.106E6,&
1504                               CH2O = 4.2E6
1505
1506
1507 !urban
1508         IF( VEGTYP == ISURBAN ) then
1509             CSOIL_LOC=3.0E6
1510         ELSE
1511             CSOIL_LOC=CSOIL
1512         ENDIF
1513
1514 ! ----------------------------------------------------------------------
1515 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
1516 ! ----------------------------------------------------------------------
1517        ITAVG = .TRUE.
1518 ! ----------------------------------------------------------------------
1519 ! BEGIN SECTION FOR TOP SOIL LAYER
1520 ! ----------------------------------------------------------------------
1521 ! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
1522 ! ----------------------------------------------------------------------
1523       HCPCT = SH2O (1)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC (1))&
1524        * CAIR                                                           &
1525               + ( SMC (1) - SH2O (1) )* CICE
1526
1527 ! ----------------------------------------------------------------------
1528 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1529 ! ----------------------------------------------------------------------
1530       DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
1531       AI (1) = 0.0
1532       CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
1533
1534 ! ----------------------------------------------------------------------
1535 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
1536 ! LAYERS.  THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
1537 ! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
1538 ! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
1539 ! ----------------------------------------------------------------------
1540       BI (1) = - CI (1) + DF1 / (0.5 * ZSOIL (1) * ZSOIL (1)* HCPCT *    &
1541        ZZ1)
1542       DTSDZ = (STC (1) - STC (2)) / ( -0.5 * ZSOIL (2))
1543       SSOIL = DF1 * (STC (1) - YY) / (0.5 * ZSOIL (1) * ZZ1)
1544 !      RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
1545       DENOM = (ZSOIL (1) * HCPCT)
1546
1547 ! ----------------------------------------------------------------------
1548 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
1549 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
1550 ! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
1551 ! ----------------------------------------------------------------------
1552 !      QTOT = SSOIL - DF1*DTSDZ
1553       RHSTS (1) = (DF1 * DTSDZ - SSOIL) / DENOM
1554
1555 ! ----------------------------------------------------------------------
1556 ! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
1557 ! ----------------------------------------------------------------------
1558       QTOT = -1.0* RHSTS (1)* DENOM
1559
1560 ! ----------------------------------------------------------------------
1561 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
1562 ! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
1563 ! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC).  IF SNOWPACK CONTENT IS
1564 ! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP.  IF
1565 ! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
1566 ! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK.  THEN
1567 ! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
1568 ! LATER IN FUNCTION SUBROUTINE SNKSRC
1569 ! ----------------------------------------------------------------------
1570       SICE = SMC (1) - SH2O (1)
1571       IF (ITAVG) THEN
1572          TSURF = (YY + (ZZ1-1) * STC (1)) / ZZ1
1573 ! ----------------------------------------------------------------------
1574 ! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
1575 ! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
1576 ! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
1577 ! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
1578 ! ----------------------------------------------------------------------
1579          CALL TBND (STC (1),STC (2),ZSOIL,ZBOT,1,NSOIL,TBK)
1580          IF ( (SICE > 0.) .OR. (STC (1) < T0) .OR.                         &
1581             (TSURF < T0) .OR. (TBK < T0) ) THEN
1582 !          TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1),
1583             CALL TMPAVG (TAVG,TSURF,STC (1),TBK,ZSOIL,NSOIL,1)
1584             CALL SNKSRC (TSNSR,TAVG,SMC (1),SH2O (1),                      &
1585                           ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1586 !          RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1587             RHSTS (1) = RHSTS (1) - TSNSR / DENOM
1588          END IF
1589       ELSE
1590 !          TSNSR = SNKSRC (STC(1),SMC(1),SH2O(1),
1591          IF ( (SICE > 0.) .OR. (STC (1) < T0) ) THEN
1592             CALL SNKSRC (TSNSR,STC (1),SMC (1),SH2O (1),                   &
1593                           ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1594 !          RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1595             RHSTS (1) = RHSTS (1) - TSNSR / DENOM
1596          END IF
1597 ! ----------------------------------------------------------------------
1598 ! THIS ENDS SECTION FOR TOP SOIL LAYER.
1599 ! ----------------------------------------------------------------------
1600       END IF
1601
1602 ! INITIALIZE DDZ2
1603 ! ----------------------------------------------------------------------
1604
1605       DDZ2 = 0.0
1606       DF1K = DF1
1607
1608 ! ----------------------------------------------------------------------
1609 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
1610 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
1611 ! ----------------------------------------------------------------------
1612 ! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
1613 ! ----------------------------------------------------------------------
1614       DO K = 2,NSOIL
1615          HCPCT = SH2O (K)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC (  &
1616                 K))* CAIR + ( SMC (K) - SH2O (K) )* CICE
1617 ! ----------------------------------------------------------------------
1618 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
1619 ! ----------------------------------------------------------------------
1620 ! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
1621 ! ----------------------------------------------------------------------
1622          IF (K /= NSOIL) THEN
1623
1624 ! ----------------------------------------------------------------------
1625 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
1626 ! ----------------------------------------------------------------------
1627             CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
1628
1629 !urban
1630        IF ( VEGTYP == ISURBAN ) DF1N = 3.24
1631
1632             DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
1633
1634 ! ----------------------------------------------------------------------
1635 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
1636 ! ----------------------------------------------------------------------
1637             DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
1638             DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
1639
1640 ! ----------------------------------------------------------------------
1641 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):  CALCULATE
1642 ! TEMP AT BOTTOM OF LAYER.
1643 ! ----------------------------------------------------------------------
1644             CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) *      &
1645        HCPCT)
1646             IF (ITAVG) THEN
1647                CALL TBND (STC (K),STC (K +1),ZSOIL,ZBOT,K,NSOIL,TBK1)
1648             END IF
1649
1650          ELSE
1651 ! ----------------------------------------------------------------------
1652 ! SPECIAL CASE OF BOTTOM SOIL LAYER:  CALCULATE THERMAL DIFFUSIVITY FOR
1653 ! BOTTOM LAYER.
1654 ! ----------------------------------------------------------------------
1655
1656 ! ----------------------------------------------------------------------
1657 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
1658 ! ----------------------------------------------------------------------
1659             CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
1660
1661
1662 !urban
1663        IF ( VEGTYP == ISURBAN ) DF1N = 3.24
1664
1665             DENOM = .5 * (ZSOIL (K -1) + ZSOIL (K)) - ZBOT
1666
1667 ! ----------------------------------------------------------------------
1668 ! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
1669 ! ----------------------------------------------------------------------
1670             DTSDZ2 = (STC (K) - TBOT) / DENOM
1671
1672 ! ----------------------------------------------------------------------
1673 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):  CALCULATE
1674 ! TEMP AT BOTTOM OF LAST LAYER.
1675 ! ----------------------------------------------------------------------
1676             CI (K) = 0.
1677             IF (ITAVG) THEN
1678                CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
1679             END IF
1680 ! ----------------------------------------------------------------------
1681 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
1682          END IF
1683 ! ----------------------------------------------------------------------
1684 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
1685 ! ----------------------------------------------------------------------
1686          DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
1687          RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
1688          QTOT = -1.0* DENOM * RHSTS (K)
1689
1690          SICE = SMC (K) - SH2O (K)
1691          IF (ITAVG) THEN
1692             CALL TMPAVG (TAVG,TBK,STC (K),TBK1,ZSOIL,NSOIL,K)
1693 !            TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL,
1694             IF ( (SICE > 0.) .OR. (STC (K) < T0) .OR.                   &
1695                (TBK .lt. T0) .OR. (TBK1 .lt. T0) ) THEN
1696                CALL SNKSRC (TSNSR,TAVG,SMC (K),SH2O (K),ZSOIL,NSOIL,    &
1697                              SMCMAX,PSISAT,BEXP,DT,K,QTOT)
1698                RHSTS (K) = RHSTS (K) - TSNSR / DENOM
1699             END IF
1700          ELSE
1701 !            TSNSR = SNKSRC(STC(K),SMC(K),SH2O(K),ZSOIL,NSOIL,
1702             IF ( (SICE > 0.) .OR. (STC (K) < T0) ) THEN
1703                CALL SNKSRC (TSNSR,STC (K),SMC (K),SH2O (K),ZSOIL,NSOIL, &
1704                              SMCMAX,PSISAT,BEXP,DT,K,QTOT)
1705                RHSTS (K) = RHSTS (K) - TSNSR / DENOM
1706             END IF
1707          END IF
1708
1709 ! ----------------------------------------------------------------------
1710 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
1711 ! ----------------------------------------------------------------------
1712          AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
1713
1714 ! ----------------------------------------------------------------------
1715 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
1716 ! ----------------------------------------------------------------------
1717          BI (K) = - (AI (K) + CI (K))
1718          TBK = TBK1
1719          DF1K = DF1N
1720          DTSDZ = DTSDZ2
1721          DDZ = DDZ2
1722       END DO
1723 ! ----------------------------------------------------------------------
1724   END SUBROUTINE HRT
1725 ! ----------------------------------------------------------------------
1726
1727       SUBROUTINE HRTICE (RHSTS,STC,TBOT,ICE,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
1728
1729 ! ----------------------------------------------------------------------
1730 ! SUBROUTINE HRTICE
1731 ! ----------------------------------------------------------------------
1732 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1733 ! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE (ICE=1) OR GLACIAL
1734 ! ICE (ICE=-1). COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE
1735 ! TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1736 !
1737 ! (NOTE:  THIS SUBROUTINE ONLY CALLED FOR SEA-ICE OR GLACIAL ICE, BUT
1738 ! NOT FOR NON-GLACIAL LAND (ICE = 0).
1739 ! ----------------------------------------------------------------------
1740       IMPLICIT NONE
1741
1742
1743       INTEGER, INTENT(IN)    :: NSOIL
1744       INTEGER                :: K
1745
1746       REAL,    INTENT(IN)    :: DF1,YY,ZZ1
1747       REAL, DIMENSION(1:NSOIL), INTENT(OUT):: AI, BI,CI
1748       REAL, DIMENSION(1:NSOIL), INTENT(IN) :: STC, ZSOIL
1749       REAL, DIMENSION(1:NSOIL), INTENT(OUT):: RHSTS
1750       REAL,                     INTENT(IN) :: TBOT
1751       INTEGER,                  INTENT(IN) :: ICE
1752       REAL                   :: DDZ,DDZ2,DENOM,DTSDZ,DTSDZ2,SSOIL,       &
1753                                 ZBOT
1754       REAL                   :: HCPCT
1755       REAL :: DF1K
1756       REAL :: DF1N
1757       REAL :: ZMD
1758
1759 ! ----------------------------------------------------------------------
1760 ! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
1761 ! HCPCT = 1880.0*917.0.
1762 ! ----------------------------------------------------------------------
1763       IF ( ICE == 1 ) THEN
1764          ! Sea-ice values
1765          HCPCT = 1.72396E+6
1766       ELSEIF (ICE == -1) THEN
1767 ! SET A NOMINAL UNIVERSAL VALUE OF GLACIAL-ICE SPECIFIC HEAT CAPACITY,
1768 !   HCPCT = 2100.0*900.0 = 1.89000E+6 (SOURCE:  BOB GRUMBINE, 2005)
1769 !   TBOT PASSED IN AS ARGUMENT, VALUE FROM GLOBAL DATA SET
1770          !
1771          ! A least-squares fit for the four points provided by
1772          ! Keith Hines for the Yen (1981) values for Antarctic
1773          ! snow firn.
1774          !
1775          HCPCT = 1.E6 * (0.8194 - 0.1309*0.5*ZSOIL(1))
1776          DF1K = DF1
1777       ENDIF
1778
1779 ! ----------------------------------------------------------------------
1780 ! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
1781 ! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
1782 ! ----------------------------------------------------------------------
1783 ! SET ICE PACK DEPTH.  USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
1784 ! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK).  ASSUME ICE
1785 ! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
1786 ! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
1787 ! ----------------------------------------------------------------------
1788 ! ----------------------------------------------------------------------
1789 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1790 ! ----------------------------------------------------------------------
1791       IF (ICE == 1) THEN
1792          ZBOT = ZSOIL (NSOIL)
1793       ELSE IF (ICE == -1) THEN
1794          ZBOT = -25.0
1795       ENDIF
1796       DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
1797       AI (1) = 0.0
1798       CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
1799
1800 ! ----------------------------------------------------------------------
1801 ! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
1802 ! RECALC/ADJUST THE SOIL HEAT FLUX.  USE THE GRADIENT AND FLUX TO CALC
1803 ! RHSTS FOR THE TOP SOIL LAYER.
1804 ! ----------------------------------------------------------------------
1805       BI (1) = - CI (1) + DF1/ (0.5 * ZSOIL (1) * ZSOIL (1) * HCPCT *    &
1806        ZZ1)
1807       DTSDZ = ( STC (1) - STC (2) ) / ( -0.5 * ZSOIL (2) )
1808       SSOIL = DF1 * ( STC (1) - YY ) / ( 0.5 * ZSOIL (1) * ZZ1 )
1809
1810 ! ----------------------------------------------------------------------
1811 ! INITIALIZE DDZ2
1812 ! ----------------------------------------------------------------------
1813       RHSTS (1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL (1) * HCPCT )
1814
1815 ! ----------------------------------------------------------------------
1816 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
1817 ! ----------------------------------------------------------------------
1818       DDZ2 = 0.0
1819       DF1K = DF1
1820       DF1N = DF1
1821       DO K = 2,NSOIL
1822
1823        IF ( ICE == -1 ) THEN
1824           ZMD = 0.5 * (ZSOIL(K)+ZSOIL(K-1))
1825           ! For the land-ice case
1826 ! kmh 09/03/2006 use Yen (1981)'s values for Antarctic snow firn
1827 !         IF ( K .eq. 2 ) HCPCT = 0.855108E6
1828 !         IF ( K .eq. 3 ) HCPCT = 0.922906E6
1829 !         IF ( K .eq. 4 ) HCPCT = 1.009986E6
1830
1831           ! Least squares fit to the four points supplied by Keith Hines
1832           ! from Yen (1981) for Antarctic snow firn.  Not optimal, but
1833           ! probably better than just a constant.
1834           HCPCT = 1.E6 * ( 0.8194 - 0.1309*ZMD )
1835
1836 !         IF ( K .eq. 2 ) DF1N = 0.345356
1837 !         IF ( K .eq. 3 ) DF1N = 0.398777
1838 !         IF ( K .eq. 4 ) DF1N = 0.472653
1839
1840           ! Least squares fit to the three points supplied by Keith Hines
1841           ! from Yen (1981) for Antarctic snow firn.  Not optimal, but
1842           ! probably better than just a constant.
1843           DF1N = 0.32333 - ( 0.10073 * ZMD )
1844        ENDIF
1845 ! ----------------------------------------------------------------------
1846 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
1847 ! ----------------------------------------------------------------------
1848          IF (K /= NSOIL) THEN
1849             DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
1850
1851 ! ----------------------------------------------------------------------
1852 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
1853 ! ----------------------------------------------------------------------
1854             DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
1855             DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
1856             CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K))*HCPCT)
1857
1858 ! ----------------------------------------------------------------------
1859 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
1860 ! ----------------------------------------------------------------------
1861          ELSE
1862
1863 ! ----------------------------------------------------------------------
1864 ! SET MATRIX COEF, CI TO ZERO.
1865 ! ----------------------------------------------------------------------
1866             DTSDZ2 = (STC (K) - TBOT)/ (.5 * (ZSOIL (K -1) + ZSOIL (K)) &
1867                      - ZBOT)
1868             CI (K) = 0.
1869 ! ----------------------------------------------------------------------
1870 ! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
1871 ! ----------------------------------------------------------------------
1872          END IF
1873          DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
1874
1875 ! ----------------------------------------------------------------------
1876 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
1877 ! ----------------------------------------------------------------------
1878          RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
1879          AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
1880
1881 ! ----------------------------------------------------------------------
1882 ! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
1883 ! ----------------------------------------------------------------------
1884          BI (K) = - (AI (K) + CI (K))
1885          DF1K = DF1N
1886          DTSDZ = DTSDZ2
1887          DDZ = DDZ2
1888       END DO
1889 ! ----------------------------------------------------------------------
1890   END SUBROUTINE HRTICE
1891 ! ----------------------------------------------------------------------
1892
1893       SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
1894
1895 ! ----------------------------------------------------------------------
1896 ! SUBROUTINE HSTEP
1897 ! ----------------------------------------------------------------------
1898 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
1899 ! ----------------------------------------------------------------------
1900       IMPLICIT NONE
1901       INTEGER, INTENT(IN)  :: NSOIL
1902       INTEGER              :: K
1903
1904       REAL, DIMENSION(1:NSOIL), INTENT(IN):: STCIN
1905       REAL, DIMENSION(1:NSOIL), INTENT(OUT):: STCOUT
1906       REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: RHSTS
1907       REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: AI,BI,CI
1908       REAL, DIMENSION(1:NSOIL) :: RHSTSin
1909       REAL, DIMENSION(1:NSOIL) :: CIin
1910       REAL                 :: DT
1911
1912 ! ----------------------------------------------------------------------
1913 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
1914 ! ----------------------------------------------------------------------
1915       DO K = 1,NSOIL
1916          RHSTS (K) = RHSTS (K) * DT
1917          AI (K) = AI (K) * DT
1918          BI (K) = 1. + BI (K) * DT
1919          CI (K) = CI (K) * DT
1920       END DO
1921 ! ----------------------------------------------------------------------
1922 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
1923 ! ----------------------------------------------------------------------
1924       DO K = 1,NSOIL
1925          RHSTSin (K) = RHSTS (K)
1926       END DO
1927       DO K = 1,NSOIL
1928          CIin (K) = CI (K)
1929       END DO
1930 ! ----------------------------------------------------------------------
1931 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
1932 ! ----------------------------------------------------------------------
1933       CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
1934 ! ----------------------------------------------------------------------
1935 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
1936 ! ----------------------------------------------------------------------
1937       DO K = 1,NSOIL
1938          STCOUT (K) = STCIN (K) + CI (K)
1939       END DO
1940 ! ----------------------------------------------------------------------
1941   END SUBROUTINE HSTEP
1942 ! ----------------------------------------------------------------------
1943
1944       SUBROUTINE NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT,                 &
1945                          SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC,      &
1946                          SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI,    &
1947                          SSOIL,                                         &
1948                          STC,EPSCA,BEXP,PC,RCH,RR,CFACTR,               &
1949                          SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL,           &
1950                          DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2,         &
1951                          RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS,        &
1952                          QUARTZ,FXEXP,CSOIL,                            &
1953                          BETA,DRIP,DEW,FLX1,FLX3,VEGTYP,ISURBAN)
1954
1955 ! ----------------------------------------------------------------------
1956 ! SUBROUTINE NOPAC
1957 ! ----------------------------------------------------------------------
1958 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
1959 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
1960 ! PRESENT.
1961 ! ----------------------------------------------------------------------
1962       IMPLICIT NONE
1963
1964       INTEGER, INTENT(IN)  :: ICE, NROOT,NSOIL,VEGTYP
1965       INTEGER, INTENT(IN)  :: ISURBAN
1966       INTEGER              :: K
1967
1968       REAL, INTENT(IN)     :: BEXP,CFACTR, CMCMAX,CSOIL,DKSAT,DT,DWSAT, &
1969                               EPSCA,ETP,FDOWN,F1,FXEXP,FRZFACT,KDT,PC,  &
1970                               PRCP,PSISAT,Q2,QUARTZ,RCH,RR,SBETA,SFCTMP,&
1971                               SHDFAC,SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, &
1972                               T24,TBOT,TH2,ZBOT,EMISSI
1973       REAL, INTENT(INOUT)  :: CMC,BETA,T1
1974       REAL, INTENT(OUT)    :: DEW,DRIP,EC,EDIR,ETA,ETT,FLX1,FLX3,       &
1975                               RUNOFF1,RUNOFF2,RUNOFF3,SSOIL
1976       REAL, DIMENSION(1:NSOIL),INTENT(IN)     :: RTDIS,ZSOIL
1977       REAL, DIMENSION(1:NSOIL),INTENT(OUT)    :: ET
1978       REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
1979       REAL, DIMENSION(1:NSOIL) :: ET1
1980       REAL                 :: EC1,EDIR1,ETT1,DF1,ETA1,ETP1,PRCP1,YY,    &
1981                               YYNUM,ZZ1
1982
1983 ! ----------------------------------------------------------------------
1984 ! EXECUTABLE CODE BEGINS HERE:
1985 ! CONVERT ETP Fnd PRCP FROM KG M-2 S-1 TO M S-1 AND INITIALIZE DEW.
1986 ! ----------------------------------------------------------------------
1987       PRCP1 = PRCP * 0.001
1988       ETP1 = ETP * 0.001
1989       DEW = 0.0
1990 ! ----------------------------------------------------------------------
1991 ! INITIALIZE EVAP TERMS.
1992 ! ----------------------------------------------------------------------
1993       EDIR = 0.
1994       EDIR1 = 0.
1995       EC1 = 0.
1996       EC = 0.
1997       DO K = 1,NSOIL
1998         ET(K) = 0.
1999         ET1(K) = 0.
2000       END DO
2001       ETT = 0.
2002       ETT1 = 0.
2003
2004       IF (ETP > 0.0) THEN
2005          CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,                  &
2006                       SH2O,                                             &
2007                       SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,                &
2008                       SMCREF,SHDFAC,CMCMAX,                             &
2009                       SMCDRY,CFACTR,                                    &
2010                        EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2011          CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL,                      &
2012                       SH2O,SLOPE,KDT,FRZFACT,                           &
2013                       SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT,                   &
2014                       SHDFAC,CMCMAX,                                    &
2015                       RUNOFF1,RUNOFF2,RUNOFF3,                          &
2016                       EDIR1,EC1,ET1,                                    &
2017                       DRIP)
2018
2019 ! ----------------------------------------------------------------------
2020 ! CONVERT MODELED EVAPOTRANSPIRATION FROM  M S-1  TO  KG M-2 S-1.
2021 ! ----------------------------------------------------------------------
2022
2023          ETA = ETA1 * 1000.0
2024
2025 ! ----------------------------------------------------------------------
2026 ! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
2027 ! ETP1 TO ZERO).
2028 ! ----------------------------------------------------------------------
2029       ELSE
2030          DEW = - ETP1
2031
2032 ! ----------------------------------------------------------------------
2033 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
2034 ! ----------------------------------------------------------------------
2035
2036          PRCP1 = PRCP1+ DEW
2037          CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL,                      &
2038                       SH2O,SLOPE,KDT,FRZFACT,                           &
2039                       SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT,                   &
2040                       SHDFAC,CMCMAX,                                    &
2041                       RUNOFF1,RUNOFF2,RUNOFF3,                          &
2042                       EDIR1,EC1,ET1,                                    &
2043                       DRIP)
2044
2045 ! ----------------------------------------------------------------------
2046 ! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
2047 ! ----------------------------------------------------------------------
2048 !         ETA = ETA1 * 1000.0
2049       END IF
2050
2051 ! ----------------------------------------------------------------------
2052 ! BASED ON ETP AND E VALUES, DETERMINE BETA
2053 ! ----------------------------------------------------------------------
2054
2055       IF ( ETP <= 0.0 ) THEN
2056          BETA = 0.0
2057          ETA = ETP
2058          IF ( ETP < 0.0 ) THEN
2059             BETA = 1.0
2060          END IF
2061       ELSE
2062          BETA = ETA / ETP
2063       END IF
2064
2065 ! ----------------------------------------------------------------------
2066 ! CONVERT MODELED EVAPOTRANSPIRATION COMPONENTS 'M S-1' TO 'KG M-2 S-1'.
2067 ! ----------------------------------------------------------------------
2068       EDIR = EDIR1*1000.
2069       EC = EC1*1000.
2070       DO K = 1,NSOIL
2071         ET(K) = ET1(K)*1000.
2072       END DO
2073       ETT = ETT1*1000.
2074
2075 ! ----------------------------------------------------------------------
2076 ! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
2077 ! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
2078 ! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
2079 ! ----------------------------------------------------------------------
2080
2081       CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
2082
2083 !urban
2084       IF ( VEGTYP == ISURBAN ) DF1=3.24
2085 !
2086
2087 ! ----------------------------------------------------------------------
2088 ! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
2089 ! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
2090 ! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
2091 ! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
2092 ! ROUTINE SFLX)
2093 ! ----------------------------------------------------------------------
2094       DF1 = DF1 * EXP (SBETA * SHDFAC)
2095 ! ----------------------------------------------------------------------
2096 ! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
2097 ! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
2098 ! ----------------------------------------------------------------------
2099       YYNUM = FDOWN - EMISSI*SIGMA * T24
2100       YY = SFCTMP + (YYNUM / RCH + TH2- SFCTMP - BETA * EPSCA) / RR
2101
2102       ZZ1 = DF1 / ( -0.5 * ZSOIL (1) * RCH * RR ) + 1.0
2103
2104 !urban
2105       CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL,        &
2106                   TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE,         &
2107                   QUARTZ,CSOIL,VEGTYP,ISURBAN)
2108
2109 ! ----------------------------------------------------------------------
2110 ! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
2111 ! THEY ARE NOT USED HERE IN SNOPAC.  FLX2 (FREEZING RAIN HEAT FLUX) WAS
2112 ! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
2113 ! ----------------------------------------------------------------------
2114       FLX1 = 0.0
2115       FLX3 = 0.0
2116
2117 ! ----------------------------------------------------------------------
2118   END SUBROUTINE NOPAC
2119 ! ----------------------------------------------------------------------
2120
2121       SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
2122      &                   Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA,       &
2123      &              DQSDT2,FLX2,EMISSI_IN,SNEQV,T1,ICE,SNCOVR)
2124
2125 ! ----------------------------------------------------------------------
2126 ! SUBROUTINE PENMAN
2127 ! ----------------------------------------------------------------------
2128 ! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT.  VARIOUS
2129 ! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
2130 ! CALLING ROUTINE FOR LATER USE.
2131 ! ----------------------------------------------------------------------
2132       IMPLICIT NONE
2133       LOGICAL, INTENT(IN)     :: SNOWNG, FRZGRA
2134       REAL, INTENT(IN)        :: CH, DQSDT2,FDOWN,PRCP,                 &
2135                                  Q2, Q2SAT,SSOIL, SFCPRS, SFCTMP,       &
2136                                  T2V, TH2,EMISSI_IN,SNEQV
2137       REAL, INTENT(IN)        :: T1 , SNCOVR
2138 !
2139 ! kmh 09/13/2006
2140       INTEGER, INTENT(IN)     :: ICE
2141 ! kmh 09/03/2006
2142 !
2143       REAL, INTENT(OUT)       :: EPSCA,ETP,FLX2,RCH,RR,T24
2144       REAL                    :: A, DELTA, FNET,RAD,RHO,EMISSI,ELCP1,LVS
2145
2146       REAL, PARAMETER      :: ELCP = 2.4888E+3, LSUBC = 2.501000E+6,CP = 1004.6
2147       REAL, PARAMETER      :: LSUBS = 2.83E+6
2148
2149 ! ----------------------------------------------------------------------
2150 ! EXECUTABLE CODE BEGINS HERE:
2151 ! ----------------------------------------------------------------------
2152 ! ----------------------------------------------------------------------
2153 ! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
2154 ! ----------------------------------------------------------------------
2155         EMISSI=EMISSI_IN
2156         IF (ICE==0) THEN
2157            ELCP1  = (1.0-SNCOVR)*ELCP  + SNCOVR*ELCP*LSUBS/LSUBC
2158            LVS    = (1.0-SNCOVR)*LSUBC + SNCOVR*LSUBS
2159         ELSE
2160            IF ( T1 > 273.15 ) THEN
2161               ELCP1=ELCP
2162               LVS=LSUBC
2163            ELSE
2164               ELCP1  = ELCP*LSUBS/LSUBC
2165               LVS    = LSUBS
2166            ENDIF
2167         ENDIF
2168
2169       FLX2 = 0.0
2170 !      DELTA = ELCP * DQSDT2
2171       DELTA = ELCP1 * DQSDT2
2172       T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
2173 !      RR = T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
2174       RR = EMISSI*T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
2175       RHO = SFCPRS / (RD * T2V)
2176
2177 ! ----------------------------------------------------------------------
2178 ! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
2179 ! EFFECTS CAUSED BY FALLING PRECIPITATION.
2180 ! ----------------------------------------------------------------------
2181       RCH = RHO * CP * CH
2182       IF (.NOT. SNOWNG) THEN
2183          IF (PRCP >  0.0) RR = RR + CPH2O * PRCP / RCH
2184       ELSE
2185          RR = RR + CPICE * PRCP / RCH
2186       END IF
2187
2188 ! ----------------------------------------------------------------------
2189 ! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
2190 ! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
2191 ! ----------------------------------------------------------------------
2192 !      FNET = FDOWN - SIGMA * T24- SSOIL
2193       FNET = FDOWN -  EMISSI*SIGMA * T24- SSOIL
2194       IF (FRZGRA) THEN
2195          FLX2 = - LSUBF * PRCP
2196          FNET = FNET - FLX2
2197 ! ----------------------------------------------------------------------
2198 ! FINISH PENMAN EQUATION CALCULATIONS.
2199 ! ----------------------------------------------------------------------
2200       END IF
2201       RAD = FNET / RCH + TH2- SFCTMP
2202 !      A = ELCP * (Q2SAT - Q2)
2203       A = ELCP1 * (Q2SAT - Q2)
2204       EPSCA = (A * RR + RAD * DELTA) / (DELTA + RR)
2205 !      ETP = EPSCA * RCH / LSUBC
2206       ETP = EPSCA * RCH / LVS
2207
2208 ! ----------------------------------------------------------------------
2209   END SUBROUTINE PENMAN
2210 ! ----------------------------------------------------------------------
2211
2212       SUBROUTINE REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX,      &
2213                          TOPT,                                             &
2214                          REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,  &
2215                          PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT,          &
2216                          SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP,      &
2217                          RTDIS,SLDPTH,ZSOIL, NROOT,NSOIL,CZIL,             &
2218                          LAIMIN, LAIMAX, EMISSMIN, EMISSMAX, ALBEDOMIN,    &
2219                          ALBEDOMAX, Z0MIN, Z0MAX, CSOIL, PTU, LLANDUSE,    &
2220                          LSOIL, LOCAL,LVCOEF)
2221
2222       IMPLICIT NONE
2223 ! ----------------------------------------------------------------------
2224 ! Internally set (default valuess)
2225 ! all soil and vegetation parameters required for the execusion oF
2226 ! the Noah lsm are defined in VEGPARM.TBL, SOILPARM.TB, and GENPARM.TBL.
2227 ! ----------------------------------------------------------------------
2228 !     Vegetation parameters:
2229 !             ALBBRD: SFC background snow-free albedo
2230 !             CMXTBL: MAX CNPY Capacity
2231 !              Z0BRD: Background roughness length
2232 !             SHDFAC: Green vegetation fraction
2233 !              NROOT: Rooting depth
2234 !              RSMIN: Mimimum stomatal resistance
2235 !              RSMAX: Max. stomatal resistance
2236 !                RGL: Parameters used in radiation stress function
2237 !                 HS: Parameter used in vapor pressure deficit functio
2238 !               TOPT: Optimum transpiration air temperature.
2239 !             CMCMAX: Maximum canopy water capacity
2240 !             CFACTR: Parameter used in the canopy inteception calculation
2241 !               SNUP: Threshold snow depth (in water equivalent m) that
2242 !                     implies 100 percent snow cover
2243 !                LAI: Leaf area index
2244 !
2245 ! ----------------------------------------------------------------------
2246 !      Soil parameters:
2247 !        SMCMAX: MAX soil moisture content (porosity)
2248 !        SMCREF: Reference soil moisture  (field capacity)
2249 !        SMCWLT: Wilting point soil moisture
2250 !        SMCWLT: Air dry soil moist content limits
2251 !       SSATPSI: SAT (saturation) soil potential
2252 !         DKSAT: SAT soil conductivity
2253 !          BEXP: B parameter
2254 !        SSATDW: SAT soil diffusivity
2255 !           F1: Soil thermal diffusivity/conductivity coef.
2256 !        QUARTZ: Soil quartz content
2257 !  Modified by F. Chen (12/22/97)  to use the STATSGO soil map
2258 !  Modified By F. Chen (01/22/00)  to include PLaya, Lava, and White San
2259 !  Modified By F. Chen (08/05/02)  to include additional parameters for the Noah
2260 ! NOTE: SATDW = BB*SATDK*(SATPSI/MAXSMC)
2261 !         F11 = ALOG10(SATPSI) + BB*ALOG10(MAXSMC) + 2.0
2262 !       REFSMC1=MAXSMC*(5.79E-9/SATDK)**(1/(2*BB+3)) 5.79E-9 m/s= 0.5 mm
2263 !       REFSMC=REFSMC1+1./3.(MAXSMC-REFSMC1)
2264 !       WLTSMC1=MAXSMC*(200./SATPSI)**(-1./BB)    (Wetzel and Chang, 198
2265 !       WLTSMC=WLTSMC1-0.5*WLTSMC1
2266 ! Note: the values for playa is set for it to have a thermal conductivit
2267 ! as sand and to have a hydrulic conductivity as clay
2268 !
2269 ! ----------------------------------------------------------------------
2270 ! Class parameter 'SLOPETYP' was included to estimate linear reservoir
2271 ! coefficient 'SLOPE' to the baseflow runoff out of the bottom layer.
2272 ! lowest class (slopetyp=0) means highest slope parameter = 1.
2273 ! definition of slopetyp from 'zobler' slope type:
2274 ! slope class  percent slope
2275 ! 1            0-8
2276 ! 2            8-30
2277 ! 3            > 30
2278 ! 4            0-30
2279 ! 5            0-8 & > 30
2280 ! 6            8-30 & > 30
2281 ! 7            0-8, 8-30, > 30
2282 ! 9            GLACIAL ICE
2283 ! BLANK        OCEAN/SEA
2284 !       SLOPE_DATA: linear reservoir coefficient
2285 !       SBETA_DATA: parameter used to caluculate vegetation effect on soil heat
2286 !       FXEXP_DAT:  soil evaporation exponent used in DEVAP
2287 !       CSOIL_DATA: soil heat capacity [J M-3 K-1]
2288 !       SALP_DATA: shape parameter of  distribution function of snow cover
2289 !       REFDK_DATA and REFKDT_DATA: parameters in the surface runoff parameteriz
2290 !       FRZK_DATA: frozen ground parameter
2291 !       ZBOT_DATA: depth[M] of lower boundary soil temperature
2292 !       CZIL_DATA: calculate roughness length of heat
2293 !       SMLOW_DATA and MHIGH_DATA: two soil moisture wilt, soil moisture referen
2294 !                 parameters
2295 ! Set maximum number of soil-, veg-, and slopetyp in data statement.
2296 ! ----------------------------------------------------------------------
2297       INTEGER, PARAMETER     :: MAX_SLOPETYP=30,MAX_SOILTYP=30,MAX_VEGTYP=30
2298       LOGICAL                :: LOCAL
2299       CHARACTER (LEN=256), INTENT(IN)::  LLANDUSE, LSOIL
2300
2301 ! Veg parameters
2302       INTEGER, INTENT(IN)    :: VEGTYP
2303       INTEGER, INTENT(OUT)   :: NROOT
2304       REAL, INTENT(OUT)      :: HS,RSMIN,RGL,SHDFAC,SNUP,                   &
2305                                 CMCMAX,RSMAX,TOPT,                          &
2306                                 EMISSMIN,  EMISSMAX,                        &
2307                                 LAIMIN,    LAIMAX,                          &
2308                                 Z0MIN,     Z0MAX,                           &
2309                                 ALBEDOMIN, ALBEDOMAX
2310 ! Soil parameters
2311       INTEGER, INTENT(IN)    :: SOILTYP
2312       REAL, INTENT(OUT)      :: BEXP,DKSAT,DWSAT,F1,QUARTZ,SMCDRY,          &
2313                                 SMCMAX,SMCREF,SMCWLT,PSISAT
2314 ! General parameters
2315       INTEGER, INTENT(IN)    :: SLOPETYP,NSOIL
2316       INTEGER                :: I
2317
2318       REAL,    INTENT(OUT)   :: SLOPE,CZIL,SBETA,FXEXP,                     &
2319                                 CSOIL,SALP,FRZX,KDT,CFACTR,      &
2320                                 ZBOT,REFKDT,PTU
2321       REAL,    INTENT(OUT)   :: LVCOEF
2322       REAL,DIMENSION(1:NSOIL),INTENT(IN) :: SLDPTH,ZSOIL
2323       REAL,DIMENSION(1:NSOIL),INTENT(OUT):: RTDIS
2324       REAL                   :: FRZFACT,FRZK,REFDK
2325
2326 !      SAVE
2327 ! ----------------------------------------------------------------------
2328 !
2329                IF (SOILTYP .gt. SLCATS) THEN
2330                         CALL wrf_error_fatal ( 'Warning: too many input soil types' )
2331                END IF
2332                IF (VEGTYP .gt. LUCATS) THEN
2333                      CALL wrf_error_fatal ( 'Warning: too many input landuse types' )
2334                END IF
2335                IF (SLOPETYP .gt. SLPCATS) THEN
2336                      CALL wrf_error_fatal ( 'Warning: too many input slope types' )
2337                END IF
2338
2339 ! ----------------------------------------------------------------------
2340 !  SET-UP SOIL PARAMETERS
2341 ! ----------------------------------------------------------------------
2342       CSOIL = CSOIL_DATA
2343       BEXP = BB (SOILTYP)
2344       DKSAT = SATDK (SOILTYP)
2345       DWSAT = SATDW (SOILTYP)
2346       F1 = F11 (SOILTYP)
2347       PSISAT = SATPSI (SOILTYP)
2348       QUARTZ = QTZ (SOILTYP)
2349       SMCDRY = DRYSMC (SOILTYP)
2350       SMCMAX = MAXSMC (SOILTYP)
2351       SMCREF = REFSMC (SOILTYP)
2352       SMCWLT = WLTSMC (SOILTYP)
2353 ! ----------------------------------------------------------------------
2354 ! Set-up universal parameters (not dependent on SOILTYP, VEGTYP or
2355 ! SLOPETYP)
2356 ! ----------------------------------------------------------------------
2357       ZBOT = ZBOT_DATA
2358       SALP = SALP_DATA
2359       SBETA = SBETA_DATA
2360       REFDK = REFDK_DATA
2361       FRZK = FRZK_DATA
2362       FXEXP = FXEXP_DATA
2363       REFKDT = REFKDT_DATA
2364       PTU = 0.    ! (not used yet) to satisify intent(out)
2365       KDT = REFKDT * DKSAT / REFDK
2366       CZIL = CZIL_DATA
2367       SLOPE = SLOPE_DATA (SLOPETYP)
2368       LVCOEF = LVCOEF_DATA
2369
2370 ! ----------------------------------------------------------------------
2371 ! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
2372 ! ----------------------------------------------------------------------
2373       FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
2374       FRZX = FRZK * FRZFACT
2375
2376 ! ----------------------------------------------------------------------
2377 ! SET-UP VEGETATION PARAMETERS
2378 ! ----------------------------------------------------------------------
2379       TOPT = TOPT_DATA
2380       CMCMAX = CMCMAX_DATA
2381       CFACTR = CFACTR_DATA
2382       RSMAX = RSMAX_DATA
2383       NROOT = NROTBL (VEGTYP)
2384       SNUP = SNUPTBL (VEGTYP)
2385       RSMIN = RSTBL (VEGTYP)
2386       RGL = RGLTBL (VEGTYP)
2387       HS = HSTBL (VEGTYP)
2388       EMISSMIN  = EMISSMINTBL  (VEGTYP)
2389       EMISSMAX  = EMISSMAXTBL  (VEGTYP)
2390       LAIMIN    = LAIMINTBL    (VEGTYP)
2391       LAIMAX    = LAIMAXTBL    (VEGTYP)
2392       Z0MIN     = Z0MINTBL     (VEGTYP)
2393       Z0MAX     = Z0MAXTBL     (VEGTYP)
2394       ALBEDOMIN = ALBEDOMINTBL (VEGTYP)
2395       ALBEDOMAX = ALBEDOMAXTBL (VEGTYP)
2396
2397                IF (VEGTYP .eq. BARE) SHDFAC = 0.0
2398                IF (NROOT .gt. NSOIL) THEN
2399                   WRITE (err_message,*) 'Error: too many root layers ',  &
2400                                                  NSOIL,NROOT
2401                   CALL wrf_error_fatal ( err_message )
2402 ! ----------------------------------------------------------------------
2403 ! CALCULATE ROOT DISTRIBUTION.  PRESENT VERSION ASSUMES UNIFORM
2404 ! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
2405 ! ----------------------------------------------------------------------
2406                END IF
2407                DO I = 1,NROOT
2408                   RTDIS (I) = - SLDPTH (I)/ ZSOIL (NROOT)
2409 ! ----------------------------------------------------------------------
2410 !  SET-UP SLOPE PARAMETER
2411 ! ----------------------------------------------------------------------
2412                END DO
2413
2414 !        print*,'end of PRMRED'
2415 !       print*,'VEGTYP',VEGTYP,'SOILTYP',SOILTYP,'SLOPETYP',SLOPETYP,    &
2416 !    & 'CFACTR',CFACTR,'CMCMAX',CMCMAX,'RSMAX',RSMAX,'TOPT',TOPT,        &
2417 !    & 'REFKDT',REFKDT,'KDT',KDT,'SBETA',SBETA, 'SHDFAC',SHDFAC,         &
2418 !    &  'RSMIN',RSMIN,'RGL',RGL,'HS',HS,'ZBOT',ZBOT,'FRZX',FRZX,         &
2419 !    &  'PSISAT',PSISAT,'SLOPE',SLOPE,'SNUP',SNUP,'SALP',SALP,'BEXP',    &
2420 !    &   BEXP,                                                           &
2421 !    &  'DKSAT',DKSAT,'DWSAT',DWSAT,                                     &
2422 !    &  'SMCMAX',SMCMAX,'SMCWLT',SMCWLT,'SMCREF',SMCREF,'SMCDRY',SMCDRY, &
2423 !    &  'F1',F1,'QUARTZ',QUARTZ,'FXEXP',FXEXP,                           &
2424 !    &  'RTDIS',RTDIS,'SLDPTH',SLDPTH,'ZSOIL',ZSOIL, 'NROOT',NROOT,      &
2425 !    &  'NSOIL',NSOIL,'Z0',Z0,'CZIL',CZIL,'LAI',LAI,                     &
2426 !    &  'CSOIL',CSOIL,'PTU',PTU,                                         &
2427 !    &  'LOCAL', LOCAL
2428
2429       END  SUBROUTINE REDPRM
2430
2431       SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
2432
2433 ! ----------------------------------------------------------------------
2434 ! SUBROUTINE ROSR12
2435 ! ----------------------------------------------------------------------
2436 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
2437 ! ###                                            ### ###  ###   ###  ###
2438 ! #B(1), C(1),  0  ,  0  ,  0  ,   . . .  ,    0   # #      #   #      #
2439 ! #A(2), B(2), C(2),  0  ,  0  ,   . . .  ,    0   # #      #   #      #
2440 ! # 0  , A(3), B(3), C(3),  0  ,   . . .  ,    0   # #      #   # D(3) #
2441 ! # 0  ,  0  , A(4), B(4), C(4),   . . .  ,    0   # # P(4) #   # D(4) #
2442 ! # 0  ,  0  ,  0  , A(5), B(5),   . . .  ,    0   # # P(5) #   # D(5) #
2443 ! # .                                          .   # #  .   # = #   .  #
2444 ! # .                                          .   # #  .   #   #   .  #
2445 ! # .                                          .   # #  .   #   #   .  #
2446 ! # 0  , . . . , 0 , A(M-2), B(M-2), C(M-2),   0   # #P(M-2)#   #D(M-2)#
2447 ! # 0  , . . . , 0 ,   0   , A(M-1), B(M-1), C(M-1)# #P(M-1)#   #D(M-1)#
2448 ! # 0  , . . . , 0 ,   0   ,   0   ,  A(M) ,  B(M) # # P(M) #   # D(M) #
2449 ! ###                                            ### ###  ###   ###  ###
2450 ! ----------------------------------------------------------------------
2451       IMPLICIT NONE
2452
2453       INTEGER, INTENT(IN)   :: NSOIL
2454       INTEGER               :: K, KK
2455
2456       REAL, DIMENSION(1:NSOIL), INTENT(IN):: A, B, D
2457       REAL, DIMENSION(1:NSOIL),INTENT(INOUT):: C,P,DELTA
2458
2459 ! ----------------------------------------------------------------------
2460 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
2461 ! ----------------------------------------------------------------------
2462       C (NSOIL) = 0.0
2463       P (1) = - C (1) / B (1)
2464 ! ----------------------------------------------------------------------
2465 ! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
2466 ! ----------------------------------------------------------------------
2467
2468 ! ----------------------------------------------------------------------
2469 ! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
2470 ! ----------------------------------------------------------------------
2471       DELTA (1) = D (1) / B (1)
2472       DO K = 2,NSOIL
2473          P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
2474          DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
2475                     * P (K -1)))
2476       END DO
2477 ! ----------------------------------------------------------------------
2478 ! SET P TO DELTA FOR LOWEST SOIL LAYER
2479 ! ----------------------------------------------------------------------
2480       P (NSOIL) = DELTA (NSOIL)
2481
2482 ! ----------------------------------------------------------------------
2483 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
2484 ! ----------------------------------------------------------------------
2485       DO K = 2,NSOIL
2486          KK = NSOIL - K + 1
2487          P (KK) = P (KK) * P (KK +1) + DELTA (KK)
2488       END DO
2489 ! ----------------------------------------------------------------------
2490   END SUBROUTINE ROSR12
2491 ! ----------------------------------------------------------------------
2492
2493
2494       SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL,  &
2495                          TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE,  &
2496                          QUARTZ,CSOIL,VEGTYP,ISURBAN)
2497
2498 ! ----------------------------------------------------------------------
2499 ! SUBROUTINE SHFLX
2500 ! ----------------------------------------------------------------------
2501 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
2502 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
2503 ! ON THE TEMPERATURE.
2504 ! ----------------------------------------------------------------------
2505       IMPLICIT NONE
2506
2507       INTEGER, INTENT(IN)   :: ICE, NSOIL, VEGTYP, ISURBAN
2508       INTEGER               :: I
2509
2510       REAL, INTENT(IN)      :: BEXP,CSOIL,DF1,DT,F1,PSISAT,QUARTZ,     &
2511                                SMCMAX, SMCWLT, TBOT,YY, ZBOT,ZZ1
2512       REAL, INTENT(INOUT)   :: T1
2513       REAL, INTENT(OUT)     :: SSOIL
2514       REAL, DIMENSION(1:NSOIL), INTENT(IN)    :: SMC,ZSOIL
2515       REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O
2516       REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
2517       REAL, DIMENSION(1:NSOIL)             :: AI, BI, CI, STCF,RHSTS
2518       REAL, PARAMETER       :: T0 = 273.15
2519
2520 ! ----------------------------------------------------------------------
2521 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
2522 ! ----------------------------------------------------------------------
2523
2524 ! ----------------------------------------------------------------------
2525 ! SEA-ICE CASE, GLACIAL ICE CASE
2526 ! ----------------------------------------------------------------------
2527       IF ( ICE /= 0 ) THEN
2528
2529          CALL HRTICE (RHSTS,STC,TBOT,ICE,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2530
2531          CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
2532
2533 ! ----------------------------------------------------------------------
2534 ! LAND-MASS CASE
2535 ! ----------------------------------------------------------------------
2536       ELSE
2537          CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT,        &
2538                     ZBOT,PSISAT,SH2O,DT,                                &
2539                     BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP,ISURBAN)
2540
2541          CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
2542       END IF
2543       DO I = 1,NSOIL
2544          STC (I) = STCF (I)
2545       END DO
2546
2547 ! ----------------------------------------------------------------------
2548 ! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
2549 ! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
2550 ! PROFILE ABOVE.  (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
2551 ! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
2552 ! DIFFERENTLY IN ROUTINE SNOPAC)
2553 ! ----------------------------------------------------------------------
2554 ! ----------------------------------------------------------------------
2555 ! CALCULATE SURFACE SOIL HEAT FLUX
2556 ! ----------------------------------------------------------------------
2557       T1 = (YY + (ZZ1- 1.0) * STC (1)) / ZZ1
2558       SSOIL = DF1 * (STC (1) - T1) / (0.5 * ZSOIL (1))
2559
2560 ! ----------------------------------------------------------------------
2561   END SUBROUTINE SHFLX
2562 ! ----------------------------------------------------------------------
2563
2564       SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL,                   &
2565      &                   SH2O,SLOPE,KDT,FRZFACT,                        &
2566      &                   SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT,                &
2567      &                   SHDFAC,CMCMAX,                                 &
2568      &                   RUNOFF1,RUNOFF2,RUNOFF3,                       &
2569      &                   EDIR,EC,ET,                                    &
2570      &                   DRIP)
2571
2572 ! ----------------------------------------------------------------------
2573 ! SUBROUTINE SMFLX
2574 ! ----------------------------------------------------------------------
2575 ! CALCULATE SOIL MOISTURE FLUX.  THE SOIL MOISTURE CONTENT (SMC - A PER
2576 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
2577 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
2578 ! FROZEN GROUND VERSION:  NEW STATES ADDED: SH2O, AND FROZEN GROUND
2579 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
2580 ! ----------------------------------------------------------------------
2581       IMPLICIT NONE
2582
2583       INTEGER, INTENT(IN)   :: NSOIL
2584       INTEGER               :: I,K
2585
2586       REAL, INTENT(IN)      :: BEXP, CMCMAX, DKSAT,DWSAT, DT, EC, EDIR,  &
2587                                KDT, PRCP1, SHDFAC, SLOPE, SMCMAX, SMCWLT
2588       REAL, INTENT(OUT)                      :: DRIP, RUNOFF1, RUNOFF2, RUNOFF3
2589       REAL, INTENT(INOUT)   :: CMC
2590       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: ET,ZSOIL
2591       REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SMC, SH2O
2592       REAL, DIMENSION(1:NSOIL)             :: AI, BI, CI, STCF,RHSTS, RHSTT, &
2593                                               SICE, SH2OA, SH2OFG
2594       REAL                  :: DUMMY, EXCESS,FRZFACT,PCPDRP,RHSCT,TRHSCT
2595       REAL :: FAC2
2596       REAL :: FLIMIT
2597
2598 ! ----------------------------------------------------------------------
2599 ! EXECUTABLE CODE BEGINS HERE.
2600 ! ----------------------------------------------------------------------
2601 ! ----------------------------------------------------------------------
2602 ! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
2603 ! ----------------------------------------------------------------------
2604       DUMMY = 0.
2605
2606 ! ----------------------------------------------------------------------
2607 ! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
2608 ! CMC.  IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
2609 ! FALL TO THE GRND.
2610 ! ----------------------------------------------------------------------
2611       RHSCT = SHDFAC * PRCP1- EC
2612       DRIP = 0.
2613       TRHSCT = DT * RHSCT
2614       EXCESS = CMC + TRHSCT
2615
2616 ! ----------------------------------------------------------------------
2617 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
2618 ! SOIL
2619 ! ----------------------------------------------------------------------
2620       IF (EXCESS > CMCMAX) DRIP = EXCESS - CMCMAX
2621       PCPDRP = (1. - SHDFAC) * PRCP1+ DRIP / DT
2622
2623 ! ----------------------------------------------------------------------
2624 ! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT and SSTEP
2625 !
2626       DO I = 1,NSOIL
2627          SICE (I) = SMC (I) - SH2O (I)
2628       END DO
2629 ! ----------------------------------------------------------------------
2630 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
2631 ! TENDENCY EQUATIONS.
2632 ! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
2633 !   (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
2634 !    EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
2635 !    THE FIRST SOIL LAYER)
2636 ! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
2637 !   TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
2638 !   OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
2639 !   PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
2640 !   SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
2641 !   OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
2642 !   DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
2643 !   SOIL MOISTURE STATE
2644 ! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
2645 !   TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
2646 !   OF SECTION 2 OF KALNAY AND KANAMITSU
2647 ! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
2648 ! ----------------------------------------------------------------------
2649 !  According to Dr. Ken Mitchell's suggestion, add the second contraint
2650 !  to remove numerical instability of runoff and soil moisture
2651 !  FLIMIT is a limit value for FAC2
2652       FAC2=0.0
2653       DO I=1,NSOIL
2654          FAC2=MAX(FAC2,SH2O(I)/SMCMAX)
2655       ENDDO
2656       CALL FAC2MIT(SMCMAX,FLIMIT)
2657
2658 ! ----------------------------------------------------------------------
2659 ! FROZEN GROUND VERSION:
2660 ! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR.  SH2O & SICE STATES
2661 ! INC&UDED IN SSTEP SUBR.  FROZEN GROUND CORRECTION FACTOR, FRZFACT
2662 ! ADDED.  ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
2663 ! ----------------------------------------------------------------------
2664       IF ( ( (PCPDRP * DT) > (0.0001*1000.0* (- ZSOIL (1))* SMCMAX) )   &
2665            .OR. (FAC2 > FLIMIT) ) THEN
2666          CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL,          &
2667                     DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1,                    &
2668                     RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2669          CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX,     &
2670                         CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2671          DO K = 1,NSOIL
2672             SH2OA (K) = (SH2O (K) + SH2OFG (K)) * 0.5
2673          END DO
2674          CALL SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL,         &
2675                     DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1,                    &
2676                     RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2677          CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX,         &
2678                         CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2679
2680       ELSE
2681          CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL,          &
2682                     DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1,                    &
2683                       RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
2684          CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX,         &
2685                         CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
2686 !      RUNOF = RUNOFF
2687
2688       END IF
2689
2690 ! ----------------------------------------------------------------------
2691   END SUBROUTINE SMFLX
2692 ! ----------------------------------------------------------------------
2693
2694
2695       SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
2696
2697 ! ----------------------------------------------------------------------
2698 ! SUBROUTINE SNFRAC
2699 ! ----------------------------------------------------------------------
2700 ! CALCULATE SNOW FRACTION (0 -> 1)
2701 ! SNEQV   SNOW WATER EQUIVALENT (M)
2702 ! SNUP    THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
2703 ! SALP    TUNING PARAMETER
2704 ! SNCOVR  FRACTIONAL SNOW COVER
2705 ! ----------------------------------------------------------------------
2706       IMPLICIT NONE
2707
2708       REAL, INTENT(IN)     :: SNEQV,SNUP,SALP,SNOWH
2709       REAL, INTENT(OUT)    :: SNCOVR
2710       REAL                 :: RSNOW, Z0N
2711
2712 ! ----------------------------------------------------------------------
2713 ! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
2714 ! REDPRM) ABOVE WHICH SNOCVR=1.
2715 ! ----------------------------------------------------------------------
2716       IF (SNEQV < SNUP) THEN
2717          RSNOW = SNEQV / SNUP
2718          SNCOVR = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
2719       ELSE
2720          SNCOVR = 1.0
2721       END IF
2722
2723 !     FORMULATION OF DICKINSON ET AL. 1986
2724 !     Z0N = 0.035
2725
2726 !        SNCOVR=SNOWH/(SNOWH + 5*Z0N)
2727
2728 !     FORMULATION OF MARSHALL ET AL. 1994
2729 !        SNCOVR=SNEQV/(SNEQV + 2*Z0N)
2730
2731 ! ----------------------------------------------------------------------
2732   END SUBROUTINE SNFRAC
2733 ! ----------------------------------------------------------------------
2734
2735       SUBROUTINE SNKSRC (TSNSR,TAVG,SMC,SH2O,ZSOIL,NSOIL,               &
2736      &                      SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2737 ! ----------------------------------------------------------------------
2738 ! SUBROUTINE SNKSRC
2739 ! ----------------------------------------------------------------------
2740 ! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
2741 ! AVAILABLE LIQUED WATER.
2742 ! ----------------------------------------------------------------------
2743       IMPLICIT NONE
2744
2745       INTEGER, INTENT(IN)   :: K,NSOIL
2746       REAL, INTENT(IN)      :: BEXP, DT, PSISAT, QTOT, SMC, SMCMAX,    &
2747                                TAVG
2748       REAL, INTENT(INOUT)   :: SH2O
2749
2750       REAL, DIMENSION(1:NSOIL), INTENT(IN):: ZSOIL
2751
2752       REAL                  :: DF, DZ, DZH, FREE, TSNSR,               &
2753                                TDN, TM, TUP, TZ, X0, XDN, XH2O, XUP
2754
2755       REAL, PARAMETER       :: DH2O = 1.0000E3, HLICE = 3.3350E5,      &
2756                                T0 = 2.7315E2
2757
2758       IF (K == 1) THEN
2759          DZ = - ZSOIL (1)
2760       ELSE
2761          DZ = ZSOIL (K -1) - ZSOIL (K)
2762       END IF
2763 ! ----------------------------------------------------------------------
2764 ! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
2765 ! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
2766 ! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
2767 ! 104, PG 19573).  (ASIDE:  LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
2768 ! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
2769 ! ----------------------------------------------------------------------
2770 !      FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
2771
2772 ! ----------------------------------------------------------------------
2773 ! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
2774 ! VOL. 104, PG 19573.)  THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
2775 ! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
2776 ! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
2777 ! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
2778 ! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
2779 ! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
2780 ! ----------------------------------------------------------------------
2781       CALL FRH2O (FREE,TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
2782
2783 ! ----------------------------------------------------------------------
2784 ! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
2785 ! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
2786 ! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
2787 ! ----------------------------------------------------------------------
2788       XH2O = SH2O + QTOT * DT / (DH2O * HLICE * DZ)
2789       IF ( XH2O < SH2O .AND. XH2O < FREE) THEN
2790          IF ( FREE > SH2O ) THEN
2791             XH2O = SH2O
2792          ELSE
2793             XH2O = FREE
2794          END IF
2795       END IF
2796 ! ----------------------------------------------------------------------
2797 ! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
2798 ! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
2799 ! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
2800 ! ----------------------------------------------------------------------
2801       IF ( XH2O > SH2O .AND. XH2O > FREE ) THEN
2802          IF ( FREE < SH2O ) THEN
2803             XH2O = SH2O
2804          ELSE
2805             XH2O = FREE
2806          END IF
2807       END IF
2808
2809 ! ----------------------------------------------------------------------
2810 ! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
2811 ! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
2812 ! ----------------------------------------------------------------------
2813 !      SNKSRC = -DH2O*HLICE*DZ*(XH2O-SH2O)/DT
2814       IF (XH2O < 0.) XH2O = 0.
2815       IF (XH2O > SMC) XH2O = SMC
2816       TSNSR = - DH2O * HLICE * DZ * (XH2O - SH2O)/ DT
2817       SH2O = XH2O
2818
2819 ! ----------------------------------------------------------------------
2820   END SUBROUTINE SNKSRC
2821 ! ----------------------------------------------------------------------
2822
2823       SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT,   &
2824                           SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,            &
2825                           SBETA,DF1,                                    &
2826                           Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA,&
2827                          SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
2828                           SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT,          &
2829                           ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1,   &
2830                           RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT,  &
2831                           ICE,RTDIS,QUARTZ,FXEXP,CSOIL,                 &
2832                           BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI,&
2833                           RIBB,SOLDN,                                   &
2834                           ISURBAN,                                      &
2835
2836                           VEGTYP)
2837
2838 ! ----------------------------------------------------------------------
2839 ! SUBROUTINE SNOPAC
2840 ! ----------------------------------------------------------------------
2841 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
2842 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
2843 ! PRESENT.
2844 ! ----------------------------------------------------------------------
2845       IMPLICIT NONE
2846
2847       INTEGER, INTENT(IN)   :: ICE, NROOT, NSOIL,VEGTYP
2848       INTEGER, INTENT(IN)   :: ISURBAN
2849       INTEGER               :: K
2850 !
2851 ! kmh 09/03/2006 add IT16 for surface temperature iteration
2852 !
2853       INTEGER               :: IT16
2854       LOGICAL, INTENT(IN)   :: SNOWNG
2855       REAL, INTENT(IN)      :: BEXP,CFACTR, CMCMAX,CSOIL,DF1,DKSAT,     &
2856                                DT,DWSAT, EPSCA,FDOWN,F1,FXEXP,          &
2857                                FRZFACT,KDT,PC, PRCP,PSISAT,Q2,QUARTZ,   &
2858                                RCH,RR,SBETA,SFCPRS, SFCTMP, SHDFAC,     &
2859                                SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, T24,  &
2860                                TBOT,TH2,ZBOT,EMISSI,SOLDN
2861       REAL, INTENT(INOUT)   :: CMC, BETA, ESD,FLX2,PRCPF,SNOWH,SNCOVR,  &
2862                                SNDENS, T1, RIBB, ETP
2863       REAL, INTENT(OUT)     :: DEW,DRIP,EC,EDIR, ETNS, ESNOW,ETT,       &
2864                                FLX1,FLX3, RUNOFF1,RUNOFF2,RUNOFF3,      &
2865                                SSOIL,SNOMLT
2866       REAL, DIMENSION(1:NSOIL),INTENT(IN)     :: RTDIS,ZSOIL
2867       REAL, DIMENSION(1:NSOIL),INTENT(OUT)    :: ET
2868       REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
2869       REAL, DIMENSION(1:NSOIL) :: ET1
2870       REAL                  :: DENOM,DSOIL,DTOT,EC1,EDIR1,ESDFLX,ETA,   &
2871                                ETT1, ESNOW1, ESNOW2, ETA1,ETP1,ETP2,    &
2872                                ETP3, ETNS1, ETANRG, ETAX, EX, FLX3X,    &
2873                                FRCSNO,FRCSOI, PRCP1, QSAT,RSNOW, SEH,   &
2874                                SNCOND,SSOIL1, T11,T12, T12A, T12AX,     &
2875                                T12B, T14, YY, ZZ1
2876 !                               T12B, T14, YY, ZZ1,EMISSI_S
2877 !
2878 ! kmh 01/11/2007 add T15, T16, and DTOT2 for SFC T iteration and snow heat flux
2879 !
2880       REAL                  :: T15, T16, DTOT2
2881       REAL, PARAMETER       :: ESDMIN = 1.E-6, LSUBC = 2.501000E+6,     &
2882                                LSUBS = 2.83E+6, TFREEZ = 273.15,        &
2883                                SNOEXP = 2.0
2884
2885 ! ----------------------------------------------------------------------
2886 ! EXECUTABLE CODE BEGINS HERE:
2887 ! ----------------------------------------------------------------------
2888 ! IF SEA-ICE (ICE=1) OR GLACIAL-ICE (ICE=-1), SNOWCOVER FRACTION = 1.0,
2889 ! AND SUBLIMATION IS AT THE POTENTIAL RATE.
2890 ! FOR NON-GLACIAL LAND (ICE=0), IF SNOWCOVER FRACTION < 1.0, TOTAL
2891 ! EVAPORATION < POTENTIAL DUE TO NON-POTENTIAL CONTRIBUTION FROM
2892 ! NON-SNOW COVERED FRACTION.
2893 ! ----------------------------------------------------------------------
2894 ! INITIALIZE EVAP TERMS.
2895 ! ----------------------------------------------------------------------
2896 ! conversions:
2897 ! ESNOW [KG M-2 S-1]
2898 ! ESDFLX [KG M-2 S-1] .le. ESNOW
2899 ! ESNOW1 [M S-1]
2900 ! ESNOW2 [M]
2901 ! ETP [KG M-2 S-1]
2902 ! ETP1 [M S-1]
2903 ! ETP2 [M]
2904 ! ----------------------------------------------------------------------
2905       DEW = 0.
2906       EDIR = 0.
2907       EDIR1 = 0.
2908       EC1 = 0.
2909       EC = 0.
2910 !      EMISSI_S=0.95    ! For snow
2911
2912       DO K = 1,NSOIL
2913          ET (K) = 0.
2914          ET1 (K) = 0.
2915       END DO
2916       ETT = 0.
2917       ETT1 = 0.
2918       ETNS = 0.
2919       ETNS1 = 0.
2920       ESNOW = 0.
2921       ESNOW1 = 0.
2922       ESNOW2 = 0.
2923
2924 ! ----------------------------------------------------------------------
2925 ! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO ETP1 IN M S-1
2926 ! ----------------------------------------------------------------------
2927       PRCP1 = PRCPF *0.001
2928 ! ----------------------------------------------------------------------
2929 ! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
2930 ! ----------------------------------------------------------------------
2931       BETA = 1.0
2932       IF (ETP <= 0.0) THEN
2933          IF ( ( RIBB >= 0.1 ) .AND. ( FDOWN > 150.0 ) ) THEN
2934             ETP=(MIN(ETP*(1.0-RIBB),0.)*SNCOVR/0.980 + ETP*(0.980-SNCOVR))/0.980
2935          ENDIF
2936         IF(ETP == 0.) BETA = 0.0
2937         ETP1 = ETP * 0.001
2938          DEW = -ETP1
2939          ESNOW2 = ETP1*DT
2940          ETANRG = ETP*((1.-SNCOVR)*LSUBC + SNCOVR*LSUBS)
2941       ELSE
2942          ETP1 = ETP * 0.001
2943          IF ( ICE /= 0 ) THEN
2944             ! SEA-ICE AND GLACIAL-ICE CASE
2945             ESNOW = ETP
2946             ESNOW1 = ESNOW*0.001
2947             ESNOW2 = ESNOW1*DT
2948             ETANRG = ESNOW*LSUBS
2949          ELSE IF ( ICE ==  0) THEN
2950             ! NON-GLACIAL LAND CASE
2951              IF (SNCOVR <  1.) THEN
2952                CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,           &
2953                             SH2O,                                       &
2954                             SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,          &
2955                             SMCREF,SHDFAC,CMCMAX,                       &
2956                             SMCDRY,CFACTR,                              &
2957                             EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,   &
2958                             FXEXP)
2959 ! ----------------------------------------------------------------------------
2960                EDIR1 = EDIR1* (1. - SNCOVR)
2961                EC1 = EC1* (1. - SNCOVR)
2962                DO K = 1,NSOIL
2963                   ET1 (K) = ET1 (K)* (1. - SNCOVR)
2964                END DO
2965                ETT1 = ETT1*(1.-SNCOVR)
2966 !               ETNS1 = EDIR1+ EC1+ ETT1
2967                ETNS1 = ETNS1*(1.-SNCOVR)
2968 ! ----------------------------------------------------------------------------
2969                EDIR = EDIR1*1000.
2970                EC = EC1*1000.
2971                DO K = 1,NSOIL
2972                   ET (K) = ET1 (K)*1000.
2973                END DO
2974                ETT = ETT1*1000.
2975                ETNS = ETNS1*1000.
2976 ! ----------------------------------------------------------------------
2977
2978 !   end IF (SNCOVR .lt. 1.)
2979             END IF
2980 !   end IF (ICE .ne. 1)
2981          END IF
2982          ESNOW = ETP*SNCOVR
2983          ESNOW1 = ESNOW*0.001
2984          ESNOW2 = ESNOW1*DT
2985          ETANRG = ESNOW*LSUBS + ETNS*LSUBC
2986 !   end IF (ETP .le. 0.0)
2987       END IF
2988
2989 ! ----------------------------------------------------------------------
2990 ! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
2991 ! ACCUMULATING PRECIP.  NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
2992 ! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1).  ASSUMES TEMPERATURE OF THE
2993 ! SNOWFALL STRIKING THE GROUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
2994 ! ----------------------------------------------------------------------
2995       FLX1 = 0.0
2996       IF (SNOWNG) THEN
2997          FLX1 = CPICE * PRCP * (T1- SFCTMP)
2998       ELSE
2999          IF (PRCP >  0.0) FLX1 = CPH2O * PRCP * (T1- SFCTMP)
3000 ! ----------------------------------------------------------------------
3001 ! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
3002 ! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
3003 ! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
3004 ! FLUXES.  FLX1 FROM ABOVE, FLX2 BROUGHT IN VIA COMMOM BLOCK RITE.
3005 ! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
3006 ! PENMAN.
3007 ! ----------------------------------------------------------------------
3008       END IF
3009       DSOIL = - (0.5 * ZSOIL (1))
3010       DTOT = SNOWH + DSOIL
3011       DENOM = 1.0+ DF1 / (DTOT * RR * RCH)
3012 ! surface emissivity weighted by snow cover fraction
3013 !      T12A = ( (FDOWN - FLX1 - FLX2 -                                   &
3014 !     &       ((SNCOVR*EMISSI_S)+EMISSI*(1.0-SNCOVR))*SIGMA *T24)/RCH    &
3015 !     &       + TH2 - SFCTMP - ETANRG/RCH ) / RR
3016       T12A = ( (FDOWN - FLX1- FLX2- EMISSI * SIGMA * T24)/ RCH                    &
3017                 + TH2- SFCTMP - ETANRG / RCH ) / RR
3018
3019       T12B = DF1 * STC (1) / (DTOT * RR * RCH)
3020
3021 ! ----------------------------------------------------------------------
3022 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
3023 ! MELT WILL OCCUR.  SET THE SKIN TEMP TO THIS EFFECTIVE TEMP.  REDUCE
3024 ! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
3025 ! DEPENDING ON SIGN OF ETP.
3026 ! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
3027 ! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
3028 ! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
3029 ! TO ZERO.
3030 ! ----------------------------------------------------------------------
3031 ! SUB-FREEZING BLOCK
3032 ! ----------------------------------------------------------------------
3033       T12 = (SFCTMP + T12A + T12B) / DENOM
3034       IF (T12 <=  TFREEZ) THEN
3035          T1 = T12
3036          SSOIL = DF1 * (T1- STC (1)) / DTOT
3037 !        ESD = MAX (0.0, ESD- ETP2)
3038          ESD = MAX(0.0, ESD-ESNOW2)
3039          FLX3 = 0.0
3040          EX = 0.0
3041
3042          SNOMLT = 0.0
3043 ! ----------------------------------------------------------------------
3044 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
3045 ! WILL OCCUR.  CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT.  REVISE THE
3046 ! EFFECTIVE SNOW DEPTH.  REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
3047 ! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
3048 ! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
3049 ! EX FOR USE IN SMFLX.  ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
3050 ! CALCULATE QSAT VALID AT FREEZING POINT.  NOTE THAT ESAT (SATURATION
3051 ! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
3052 ! POINT.  NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
3053 ! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
3054 ! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
3055 ! ----------------------------------------------------------------------
3056 ! ABOVE FREEZING BLOCK
3057 ! ----------------------------------------------------------------------
3058       ELSE
3059          T1 = TFREEZ * SNCOVR ** SNOEXP + T12 * (1.0- SNCOVR ** SNOEXP)
3060          BETA = 1.0
3061
3062 ! ----------------------------------------------------------------------
3063 ! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
3064 ! BETA<1
3065 ! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
3066 ! ----------------------------------------------------------------------
3067 !  kmh  12/15/2005 modify SSOIL
3068 !  kmh  09/03/2006 modify DTOT
3069          if ( DTOT .GT. 2.0*DSOIL ) then
3070            DTOT = 2.0*DSOIL
3071          endif
3072          SSOIL = DF1 * (T1- STC (1)) / DTOT
3073          IF (ESD-ESNOW2 <= ESDMIN) THEN
3074             ESD = 0.0
3075             EX = 0.0
3076             SNOMLT = 0.0
3077             FLX3 = 0.0
3078 ! ----------------------------------------------------------------------
3079 ! SUBLIMATION LESS THAN DEPTH OF SNOWPACK
3080 ! SNOWPACK (ESD) REDUCED BY ESNOW2 (DEPTH OF SUBLIMATED SNOW)
3081 ! ----------------------------------------------------------------------
3082          ELSE
3083             ESD = ESD-ESNOW2
3084             ETP3 = ETP * LSUBC
3085             SEH = RCH * (T1- TH2)
3086             T14 = T1* T1
3087             T14 = T14* T14
3088 !           FLX3 = FDOWN - FLX1 - FLX2 -                                 &
3089 !                  ((SNCOVR*EMISSI_S)+EMISSI*(1-SNCOVR))*SIGMA*T14 -     &
3090 !                    SSOIL - SEH - ETANRG
3091             FLX3 = FDOWN - FLX1- FLX2- EMISSI*SIGMA * T14- SSOIL - SEH - ETANRG
3092             IF (FLX3 <= 0.0) FLX3 = 0.0
3093 ! ----------------------------------------------------------------------
3094 ! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
3095 ! ----------------------------------------------------------------------
3096             EX = FLX3*0.001/ LSUBF
3097
3098 ! ----------------------------------------------------------------------
3099 ! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
3100 ! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
3101 ! ----------------------------------------------------------------------
3102             SNOMLT = EX * DT
3103             IF (ESD- SNOMLT >=  ESDMIN) THEN
3104                ESD = ESD- SNOMLT
3105 ! ----------------------------------------------------------------------
3106 ! SNOWMELT EXCEEDS SNOW DEPTH
3107 ! ----------------------------------------------------------------------
3108             ELSE
3109                EX = ESD / DT
3110                FLX3 = EX *1000.0* LSUBF
3111                SNOMLT = ESD
3112
3113                ESD = 0.0
3114 ! ----------------------------------------------------------------------
3115 ! END OF 'ESD .LE. ETP2' IF-BLOCK
3116 ! ----------------------------------------------------------------------
3117             END IF
3118          END IF
3119
3120 ! ----------------------------------------------------------------------
3121 ! END OF 'T12 .LE. TFREEZ' IF-BLOCK
3122 ! ----------------------------------------------------------------------
3123 ! ----------------------------------------------------------------------
3124 ! IF NON-GLACIAL LAND, ADD SNOWMELT RATE (EX) TO PRECIP RATE TO BE USED
3125 ! IN SUBROUTINE SMFLX (SOIL MOISTURE EVOLUTION) VIA INFILTRATION.
3126 !
3127 ! FOR SEA-ICE AND GLACIAL-ICE, THE SNOWMELT WILL BE ADDED TO SUBSURFACE
3128 ! RUNOFF/BASEFLOW LATER NEAR THE END OF SFLX (AFTER RETURN FROM CALL TO
3129 ! SUBROUTINE SNOPAC)
3130 ! ----------------------------------------------------------------------
3131          IF (ICE == 0) PRCP1 = PRCP1+ EX
3132
3133 ! ----------------------------------------------------------------------
3134 ! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
3135 ! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
3136 ! (BELOW).
3137 ! IF SEAICE (ICE==1) SKIP CALL TO SMFLX.
3138 ! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES FOR NON-GLACIAL LAND.
3139 ! IF SEA-ICE (ICE==1) OR GLACIAL-ICE (ICE==-1), SKIP CALL TO SMFLX,
3140 ! SINCE NO SOIL MEDIUM FOR SEA-ICE OR GLACIAL-ICE.
3141 ! ----------------------------------------------------------------------
3142       END IF
3143       IF (ICE == 0) THEN
3144          CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL,                      &
3145                       SH2O,SLOPE,KDT,FRZFACT,                           &
3146                       SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT,                   &
3147                       SHDFAC,CMCMAX,                                    &
3148                       RUNOFF1,RUNOFF2,RUNOFF3,                          &
3149                       EDIR1,EC1,ET1,                                    &
3150                       DRIP)
3151 ! ----------------------------------------------------------------------
3152 ! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
3153 ! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
3154 ! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
3155 ! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
3156 ! SNOW TOP SURFACE.  T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
3157 ! SKIN TEMP VALUE AS REVISED BY SHFLX.
3158 ! ----------------------------------------------------------------------
3159       END IF
3160       ZZ1 = 1.0
3161       YY = STC (1) -0.5* SSOIL * ZSOIL (1)* ZZ1/ DF1
3162
3163 ! ----------------------------------------------------------------------
3164 ! SHFLX WILL CALC/UPDATE THE SOIL TEMPS.  NOTE:  THE SUB-SFC HEAT FLUX
3165 ! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
3166 ! USED  IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
3167 ! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
3168 ! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
3169 ! ----------------------------------------------------------------------
3170       T11 = T1
3171       CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL,      &
3172                    TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE,        &
3173                    QUARTZ,CSOIL,VEGTYP,ISURBAN)
3174
3175 ! ----------------------------------------------------------------------
3176 ! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION.  YY IS
3177 ! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
3178 ! ----------------------------------------------------------------------
3179       IF (ICE == 0) THEN
3180          ! NON-GLACIAL LAND
3181          IF (ESD >  0.) THEN
3182             CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
3183          ELSE
3184             ESD = 0.
3185             SNOWH = 0.
3186             SNDENS = 0.
3187             SNCOND = 1.
3188             SNCOVR = 0.
3189          END IF
3190       ELSEIF (ICE == 1) THEN
3191          ! SEA-ICE
3192          IF (ESD .GE. 0.01) THEN
3193             CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
3194         ELSE
3195            ESD = 0.01
3196            SNOWH = 0.05
3197            !KWM???? SNDENS =
3198            !KWM???? SNCOND =
3199            SNCOVR = 1.0
3200         ENDIF
3201       ELSEIF (ICE == -1) THEN
3202          ! GLACIAL-ICE
3203          IF (ESD .GE. 0.10) THEN
3204             CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
3205          ELSE
3206             ESD = 0.10
3207             SNOWH = 0.50
3208             !KWM???? SNDENS =
3209             !KWM???? SNCOND =
3210             SNCOVR = 1.0
3211         ENDIF
3212       ENDIF
3213 ! ----------------------------------------------------------------------
3214   END SUBROUTINE SNOPAC
3215 ! ----------------------------------------------------------------------
3216
3217
3218       SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
3219
3220 ! ----------------------------------------------------------------------
3221 ! SUBROUTINE SNOWPACK
3222 ! ----------------------------------------------------------------------
3223 ! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
3224 ! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
3225 ! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
3226 ! KOREN, 03/25/95.
3227 ! ----------------------------------------------------------------------
3228 ! ESD     WATER EQUIVALENT OF SNOW (M)
3229 ! DTSEC   TIME STEP (SEC)
3230 ! SNOWH   SNOW DEPTH (M)
3231 ! SNDENS  SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
3232 ! TSNOW   SNOW SURFACE TEMPERATURE (K)
3233 ! TSOIL   SOIL SURFACE TEMPERATURE (K)
3234
3235 ! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
3236 ! ----------------------------------------------------------------------
3237       IMPLICIT NONE
3238
3239       INTEGER                :: IPOL, J
3240       REAL, INTENT(IN)       :: ESD, DTSEC,TSNOW,TSOIL
3241       REAL, INTENT(INOUT)    :: SNOWH, SNDENS
3242       REAL                   :: BFAC,DSX,DTHR,DW,SNOWHC,PEXP,           &
3243                                 TAVGC,TSNOWC,TSOILC,ESDC,ESDCX
3244       REAL, PARAMETER        :: C1 = 0.01, C2 = 21.0, G = 9.81,         &
3245                                 KN = 4000.0
3246 ! ----------------------------------------------------------------------
3247 ! CONVERSION INTO SIMULATION UNITS
3248 ! ----------------------------------------------------------------------
3249       SNOWHC = SNOWH *100.
3250       ESDC = ESD *100.
3251       DTHR = DTSEC /3600.
3252       TSNOWC = TSNOW -273.15
3253       TSOILC = TSOIL -273.15
3254
3255 ! ----------------------------------------------------------------------
3256 ! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
3257 ! ----------------------------------------------------------------------
3258 ! ----------------------------------------------------------------------
3259 ! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
3260 !  SNDENS=DS0*(EXP(BFAC*ESD)-1.)/(BFAC*ESD)
3261 !  BFAC=DTHR*C1*EXP(0.08*TAVGC-C2*DS0)
3262 ! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
3263 ! NUMERICALLY BELOW:
3264 !   C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
3265 !   C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
3266 ! ----------------------------------------------------------------------
3267       TAVGC = 0.5* (TSNOWC + TSOILC)
3268       IF (ESDC >  1.E-2) THEN
3269          ESDCX = ESDC
3270       ELSE
3271          ESDCX = 1.E-2
3272       END IF
3273
3274 !      DSX = SNDENS*((DEXP(BFAC*ESDC)-1.)/(BFAC*ESDC))
3275 ! ----------------------------------------------------------------------
3276 ! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
3277 ! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
3278 ! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
3279 ! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
3280 ! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
3281 ! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
3282 ! POLYNOMIAL EXPANSION.
3283
3284 ! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
3285 ! IS GOVERNED BY ITERATION LIMIT "IPOL".
3286 !      IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
3287 !            PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
3288 !       IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
3289 !       IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
3290 !       IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
3291 ! ----------------------------------------------------------------------
3292       BFAC = DTHR * C1* EXP (0.08* TAVGC - C2* SNDENS)
3293       IPOL = 4
3294       PEXP = 0.
3295 !        PEXP = (1. + PEXP)*BFAC*ESDC/REAL(J+1)
3296       DO J = IPOL,1, -1
3297          PEXP = (1. + PEXP)* BFAC * ESDCX / REAL (J +1)
3298       END DO
3299
3300       PEXP = PEXP + 1.
3301 ! ----------------------------------------------------------------------
3302 ! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
3303 ! ----------------------------------------------------------------------
3304 !     END OF KOREAN FORMULATION
3305
3306 !     BASE FORMULATION (COGLEY ET AL., 1990)
3307 !     CONVERT DENSITY FROM G/CM3 TO KG/M3
3308 !       DSM=SNDENS*1000.0
3309
3310 !       DSX=DSM+DTSEC*0.5*DSM*G*ESD/
3311 !    &      (1E7*EXP(-0.02*DSM+KN/(TAVGC+273.16)-14.643))
3312
3313 !  &   CONVERT DENSITY FROM KG/M3 TO G/CM3
3314 !       DSX=DSX/1000.0
3315
3316 !     END OF COGLEY ET AL. FORMULATION
3317
3318 ! ----------------------------------------------------------------------
3319 ! SET UPPER/LOWER LIMIT ON SNOW DENSITY
3320 ! ----------------------------------------------------------------------
3321       DSX = SNDENS * (PEXP)
3322       IF (DSX > 0.40) DSX = 0.40
3323       IF (DSX < 0.05) DSX = 0.05
3324 ! ----------------------------------------------------------------------
3325 ! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
3326 ! SNOWMELT.  ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
3327 ! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
3328 ! ----------------------------------------------------------------------
3329       SNDENS = DSX
3330       IF (TSNOWC >=  0.) THEN
3331          DW = 0.13* DTHR /24.
3332          SNDENS = SNDENS * (1. - DW) + DW
3333          IF (SNDENS >=  0.40) SNDENS = 0.40
3334 ! ----------------------------------------------------------------------
3335 ! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
3336 ! CHANGE SNOW DEPTH UNITS TO METERS
3337 ! ----------------------------------------------------------------------
3338       END IF
3339       SNOWHC = ESDC / SNDENS
3340       SNOWH = SNOWHC *0.01
3341
3342 ! ----------------------------------------------------------------------
3343   END SUBROUTINE SNOWPACK
3344 ! ----------------------------------------------------------------------
3345
3346       SUBROUTINE SNOWZ0 (SNCOVR,Z0, Z0BRD, SNOWH)
3347
3348 ! ----------------------------------------------------------------------
3349 ! SUBROUTINE SNOWZ0
3350 ! ----------------------------------------------------------------------
3351 ! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
3352 ! SNCOVR  FRACTIONAL SNOW COVER
3353 ! Z0      ROUGHNESS LENGTH (m)
3354 ! Z0S     SNOW ROUGHNESS LENGTH:=0.001 (m)
3355 ! ----------------------------------------------------------------------
3356       IMPLICIT NONE
3357       REAL, INTENT(IN)        :: SNCOVR, Z0BRD
3358       REAL, INTENT(OUT)       :: Z0
3359       REAL, PARAMETER         :: Z0S=0.001
3360       REAL, INTENT(IN)        :: SNOWH
3361       REAL BURIAL, Z0EFF
3362
3363 !m      Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0S
3364       BURIAL = 7.0*Z0BRD - SNOWH
3365       IF(BURIAL.LE.0.0007) THEN
3366         Z0EFF = Z0S
3367       ELSE
3368         Z0EFF = BURIAL/7.0
3369       ENDIF
3370
3371       Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0EFF
3372
3373 ! ----------------------------------------------------------------------
3374   END SUBROUTINE SNOWZ0
3375 ! ----------------------------------------------------------------------
3376
3377
3378       SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
3379
3380 ! ----------------------------------------------------------------------
3381 ! SUBROUTINE SNOW_NEW
3382 ! ----------------------------------------------------------------------
3383 ! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
3384 ! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
3385
3386 ! TEMP    AIR TEMPERATURE (K)
3387 ! NEWSN   NEW SNOWFALL (M)
3388 ! SNOWH   SNOW DEPTH (M)
3389 ! SNDENS  SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
3390 ! ----------------------------------------------------------------------
3391       IMPLICIT NONE
3392       REAL, INTENT(IN)        :: NEWSN, TEMP
3393       REAL, INTENT(INOUT)     :: SNDENS, SNOWH
3394       REAL                    :: DSNEW, HNEWC, SNOWHC,NEWSNC,TEMPC
3395
3396 ! ----------------------------------------------------------------------
3397 ! CONVERSION INTO SIMULATION UNITS
3398 ! ----------------------------------------------------------------------
3399       SNOWHC = SNOWH *100.
3400       NEWSNC = NEWSN *100.
3401
3402 ! ----------------------------------------------------------------------
3403 ! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
3404 ! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
3405 ! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
3406 ! VEMADOLEN, SWEDEN, 1980, 172-177PP.
3407 !-----------------------------------------------------------------------
3408       TEMPC = TEMP -273.15
3409       IF (TEMPC <=  -15.) THEN
3410          DSNEW = 0.05
3411       ELSE
3412          DSNEW = 0.05+0.0017* (TEMPC +15.)**1.5
3413       END IF
3414 ! ----------------------------------------------------------------------
3415 ! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
3416 ! ----------------------------------------------------------------------
3417       HNEWC = NEWSNC / DSNEW
3418       IF (SNOWHC + HNEWC .LT. 1.0E-3) THEN
3419          SNDENS = MAX(DSNEW,SNDENS)
3420       ELSE
3421          SNDENS = (SNOWHC * SNDENS + HNEWC * DSNEW)/ (SNOWHC + HNEWC)
3422       ENDIF
3423       SNOWHC = SNOWHC + HNEWC
3424       SNOWH = SNOWHC *0.01
3425
3426 ! ----------------------------------------------------------------------
3427   END SUBROUTINE SNOW_NEW
3428 ! ----------------------------------------------------------------------
3429
3430       SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,            &
3431                        ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1,           &
3432                        RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
3433
3434 ! ----------------------------------------------------------------------
3435 ! SUBROUTINE SRT
3436 ! ----------------------------------------------------------------------
3437 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
3438 ! WATER DIFFUSION EQUATION.  ALSO TO COMPUTE ( PREPARE ) THE MATRIX
3439 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
3440 ! ----------------------------------------------------------------------
3441       IMPLICIT NONE
3442       INTEGER, INTENT(IN)       :: NSOIL
3443       INTEGER                   :: IALP1, IOHINF, J, JJ,  K, KS
3444       REAL, INTENT(IN)          :: BEXP, DKSAT, DT, DWSAT, EDIR, FRZX,  &
3445                                    KDT, PCPDRP, SLOPE, SMCMAX, SMCWLT
3446       REAL, INTENT(OUT)         :: RUNOFF1, RUNOFF2
3447       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: ET, SH2O, SH2OA, SICE,  &
3448                                                 ZSOIL
3449       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: RHSTT
3450       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: AI, BI, CI
3451       REAL, DIMENSION(1:NSOIL)  :: DMAX
3452       REAL                      :: ACRT, DD, DDT, DDZ, DDZ2, DENOM,     &
3453                                    DENOM2,DICE, DSMDZ, DSMDZ2, DT1,     &
3454                                    FCR,INFMAX,MXSMC,MXSMC2,NUMER,PDDUM, &
3455                                    PX, SICEMAX,SLOPX, SMCAV, SSTT,      &
3456                                    SUM, VAL, WCND, WCND2, WDF, WDF2
3457       INTEGER, PARAMETER        :: CVFRZ = 3
3458
3459 ! ----------------------------------------------------------------------
3460 ! FROZEN GROUND VERSION:
3461 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
3462 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
3463 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT.  BASED
3464 ! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
3465 ! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM.  THAT IS
3466 ! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
3467 ! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
3468 ! ----------------------------------------------------------------------
3469
3470 ! ----------------------------------------------------------------------
3471 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF.  INCLUDE THE
3472 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
3473 ! MODIFIED BY Q DUAN
3474 ! ----------------------------------------------------------------------
3475 ! ----------------------------------------------------------------------
3476 ! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
3477 ! LAYERS.
3478 ! ----------------------------------------------------------------------
3479       IOHINF = 1
3480       SICEMAX = 0.0
3481       DO KS = 1,NSOIL
3482          IF (SICE (KS) >  SICEMAX) SICEMAX = SICE (KS)
3483 ! ----------------------------------------------------------------------
3484 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
3485 ! ----------------------------------------------------------------------
3486       END DO
3487       PDDUM = PCPDRP
3488       RUNOFF1 = 0.0
3489
3490 ! ----------------------------------------------------------------------
3491 ! MODIFIED BY Q. DUAN, 5/16/94
3492 ! ----------------------------------------------------------------------
3493 !        IF (IOHINF == 1) THEN
3494
3495       IF (PCPDRP /=  0.0) THEN
3496          DT1 = DT /86400.
3497          SMCAV = SMCMAX - SMCWLT
3498
3499 ! ----------------------------------------------------------------------
3500 ! FROZEN GROUND VERSION:
3501 ! ----------------------------------------------------------------------
3502          DMAX (1)= - ZSOIL (1)* SMCAV
3503
3504          DICE = - ZSOIL (1) * SICE (1)
3505          DMAX (1)= DMAX (1)* (1.0- (SH2OA (1) + SICE (1) - SMCWLT)/      &
3506                     SMCAV)
3507
3508          DD = DMAX (1)
3509
3510 ! ----------------------------------------------------------------------
3511 ! FROZEN GROUND VERSION:
3512 ! ----------------------------------------------------------------------
3513          DO KS = 2,NSOIL
3514
3515             DICE = DICE+ ( ZSOIL (KS -1) - ZSOIL (KS) ) * SICE (KS)
3516             DMAX (KS) = (ZSOIL (KS -1) - ZSOIL (KS))* SMCAV
3517             DMAX (KS) = DMAX (KS)* (1.0- (SH2OA (KS) + SICE (KS)        &
3518                         - SMCWLT)/ SMCAV)
3519             DD = DD+ DMAX (KS)
3520 ! ----------------------------------------------------------------------
3521 ! VAL = (1.-EXP(-KDT*SQRT(DT1)))
3522 ! IN BELOW, REMOVE THE SQRT IN ABOVE
3523 ! ----------------------------------------------------------------------
3524          END DO
3525          VAL = (1. - EXP ( - KDT * DT1))
3526          DDT = DD * VAL
3527          PX = PCPDRP * DT
3528          IF (PX <  0.0) PX = 0.0
3529
3530 ! ----------------------------------------------------------------------
3531 ! FROZEN GROUND VERSION:
3532 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
3533 ! ----------------------------------------------------------------------
3534          INFMAX = (PX * (DDT / (PX + DDT)))/ DT
3535          FCR = 1.
3536          IF (DICE >  1.E-2) THEN
3537             ACRT = CVFRZ * FRZX / DICE
3538             SUM = 1.
3539             IALP1 = CVFRZ - 1
3540             DO J = 1,IALP1
3541                K = 1
3542                DO JJ = J +1,IALP1
3543                   K = K * JJ
3544                END DO
3545                SUM = SUM + (ACRT ** ( CVFRZ - J)) / FLOAT (K)
3546             END DO
3547             FCR = 1. - EXP ( - ACRT) * SUM
3548          END IF
3549
3550 ! ----------------------------------------------------------------------
3551 ! CORRECTION OF INFILTRATION LIMITATION:
3552 ! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
3553 ! HYDROLIC CONDUCTIVITY
3554 ! ----------------------------------------------------------------------
3555 !         MXSMC = MAX ( SH2OA(1), SH2OA(2) )
3556          INFMAX = INFMAX * FCR
3557
3558          MXSMC = SH2OA (1)
3559          CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT,           &
3560                          SICEMAX)
3561          INFMAX = MAX (INFMAX,WCND)
3562
3563          INFMAX = MIN (INFMAX,PX/DT)
3564          IF (PCPDRP >  INFMAX) THEN
3565             RUNOFF1 = PCPDRP - INFMAX
3566             PDDUM = INFMAX
3567          END IF
3568 ! ----------------------------------------------------------------------
3569 ! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
3570 ! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
3571 ! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
3572 ! ----------------------------------------------------------------------
3573       END IF
3574
3575       MXSMC = SH2OA (1)
3576       CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT,              &
3577                     SICEMAX)
3578 ! ----------------------------------------------------------------------
3579 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
3580 ! ----------------------------------------------------------------------
3581       DDZ = 1. / ( - .5 * ZSOIL (2) )
3582       AI (1) = 0.0
3583       BI (1) = WDF * DDZ / ( - ZSOIL (1) )
3584
3585 ! ----------------------------------------------------------------------
3586 ! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
3587 ! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
3588 ! ----------------------------------------------------------------------
3589       CI (1) = - BI (1)
3590       DSMDZ = ( SH2O (1) - SH2O (2) ) / ( - .5 * ZSOIL (2) )
3591       RHSTT (1) = (WDF * DSMDZ + WCND- PDDUM + EDIR + ET (1))/ ZSOIL (1)
3592
3593 ! ----------------------------------------------------------------------
3594 ! INITIALIZE DDZ2
3595 ! ----------------------------------------------------------------------
3596       SSTT = WDF * DSMDZ + WCND+ EDIR + ET (1)
3597
3598 ! ----------------------------------------------------------------------
3599 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
3600 ! ----------------------------------------------------------------------
3601       DDZ2 = 0.0
3602       DO K = 2,NSOIL
3603          DENOM2 = (ZSOIL (K -1) - ZSOIL (K))
3604          IF (K /= NSOIL) THEN
3605
3606 ! ----------------------------------------------------------------------
3607 ! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
3608 ! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
3609 ! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
3610 ! ----------------------------------------------------------------------
3611             SLOPX = 1.
3612
3613             MXSMC2 = SH2OA (K)
3614             CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT,     &
3615                           SICEMAX)
3616 ! -----------------------------------------------------------------------
3617 ! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
3618 ! ----------------------------------------------------------------------
3619             DENOM = (ZSOIL (K -1) - ZSOIL (K +1))
3620
3621 ! ----------------------------------------------------------------------
3622 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
3623 ! ----------------------------------------------------------------------
3624             DSMDZ2 = (SH2O (K) - SH2O (K +1)) / (DENOM * 0.5)
3625             DDZ2 = 2.0 / DENOM
3626             CI (K) = - WDF2 * DDZ2 / DENOM2
3627
3628          ELSE
3629 ! ----------------------------------------------------------------------
3630 ! SLOPE OF BOTTOM LAYER IS INTRODUCED
3631 ! ----------------------------------------------------------------------
3632
3633 ! ----------------------------------------------------------------------
3634 ! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
3635 ! THIS LAYER
3636 ! ----------------------------------------------------------------------
3637             SLOPX = SLOPE
3638           CALL WDFCND (WDF2,WCND2,SH2OA (NSOIL),SMCMAX,BEXP,DKSAT,DWSAT,     &
3639                             SICEMAX)
3640
3641 ! ----------------------------------------------------------------------
3642 ! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
3643 ! ----------------------------------------------------------------------
3644
3645 ! ----------------------------------------------------------------------
3646 ! SET MATRIX COEF CI TO ZERO
3647 ! ----------------------------------------------------------------------
3648             DSMDZ2 = 0.0
3649             CI (K) = 0.0
3650 ! ----------------------------------------------------------------------
3651 ! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
3652 ! ----------------------------------------------------------------------
3653          END IF
3654          NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2- (WDF * DSMDZ)         &
3655                  - WCND+ ET (K)
3656
3657 ! ----------------------------------------------------------------------
3658 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
3659 ! ----------------------------------------------------------------------
3660          RHSTT (K) = NUMER / ( - DENOM2)
3661          AI (K) = - WDF * DDZ / DENOM2
3662
3663 ! ----------------------------------------------------------------------
3664 ! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
3665 ! RUNOFF2:  SUB-SURFACE OR BASEFLOW RUNOFF
3666 ! ----------------------------------------------------------------------
3667          BI (K) = - ( AI (K) + CI (K) )
3668          IF (K .eq. NSOIL) THEN
3669             RUNOFF2 = SLOPX * WCND2
3670          END IF
3671          IF (K .ne. NSOIL) THEN
3672             WDF = WDF2
3673             WCND = WCND2
3674             DSMDZ = DSMDZ2
3675             DDZ = DDZ2
3676          END IF
3677       END DO
3678 ! ----------------------------------------------------------------------
3679   END SUBROUTINE SRT
3680 ! ----------------------------------------------------------------------
3681
3682       SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT,              &
3683                         NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,     &
3684                         AI,BI,CI)
3685
3686 ! ----------------------------------------------------------------------
3687 ! SUBROUTINE SSTEP
3688 ! ----------------------------------------------------------------------
3689 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
3690 ! CONTENT VALUES.
3691 ! ----------------------------------------------------------------------
3692       IMPLICIT NONE
3693       INTEGER, INTENT(IN)       :: NSOIL
3694       INTEGER                   :: I, K, KK11
3695
3696       REAL, INTENT(IN)          :: CMCMAX, DT, SMCMAX
3697       REAL, INTENT(OUT)         :: RUNOFF3
3698       REAL, INTENT(INOUT)       :: CMC
3699       REAL, DIMENSION(1:NSOIL), INTENT(IN)     :: SH2OIN, SICE, ZSOIL
3700       REAL, DIMENSION(1:NSOIL), INTENT(OUT)    :: SH2OOUT
3701       REAL, DIMENSION(1:NSOIL), INTENT(INOUT)  :: RHSTT, SMC
3702       REAL, DIMENSION(1:NSOIL), INTENT(INOUT)  :: AI, BI, CI
3703       REAL, DIMENSION(1:NSOIL)  :: RHSTTin
3704       REAL, DIMENSION(1:NSOIL)  :: CIin
3705       REAL                      :: DDZ, RHSCT, STOT, WPLUS
3706
3707 ! ----------------------------------------------------------------------
3708 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
3709 ! TRI-DIAGONAL MATRIX ROUTINE.
3710 ! ----------------------------------------------------------------------
3711       DO K = 1,NSOIL
3712          RHSTT (K) = RHSTT (K) * DT
3713          AI (K) = AI (K) * DT
3714          BI (K) = 1. + BI (K) * DT
3715          CI (K) = CI (K) * DT
3716       END DO
3717 ! ----------------------------------------------------------------------
3718 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
3719 ! ----------------------------------------------------------------------
3720       DO K = 1,NSOIL
3721          RHSTTin (K) = RHSTT (K)
3722       END DO
3723       DO K = 1,NSOIL
3724          CIin (K) = CI (K)
3725       END DO
3726 ! ----------------------------------------------------------------------
3727 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
3728 ! ----------------------------------------------------------------------
3729       CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
3730 ! ----------------------------------------------------------------------
3731 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
3732 ! NEW VALUE.  MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
3733 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
3734 ! ----------------------------------------------------------------------
3735       WPLUS = 0.0
3736       RUNOFF3 = 0.
3737
3738       DDZ = - ZSOIL (1)
3739       DO K = 1,NSOIL
3740          IF (K /= 1) DDZ = ZSOIL (K - 1) - ZSOIL (K)
3741          SH2OOUT (K) = SH2OIN (K) + CI (K) + WPLUS / DDZ
3742          STOT = SH2OOUT (K) + SICE (K)
3743          IF (STOT > SMCMAX) THEN
3744             IF (K .eq. 1) THEN
3745                DDZ = - ZSOIL (1)
3746             ELSE
3747                KK11 = K - 1
3748                DDZ = - ZSOIL (K) + ZSOIL (KK11)
3749             END IF
3750             WPLUS = (STOT - SMCMAX) * DDZ
3751          ELSE
3752             WPLUS = 0.
3753          END IF
3754          SMC (K) = MAX ( MIN (STOT,SMCMAX),0.02 )
3755          SH2OOUT (K) = MAX ( (SMC (K) - SICE (K)),0.0)
3756       END DO
3757
3758 ! ----------------------------------------------------------------------
3759 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC).  CONVERT RHSCT TO
3760 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
3761 ! ----------------------------------------------------------------------
3762       RUNOFF3 = WPLUS
3763       CMC = CMC + DT * RHSCT
3764       IF (CMC < 1.E-20) CMC = 0.0
3765       CMC = MIN (CMC,CMCMAX)
3766
3767 ! ----------------------------------------------------------------------
3768   END SUBROUTINE SSTEP
3769 ! ----------------------------------------------------------------------
3770
3771       SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
3772
3773 ! ----------------------------------------------------------------------
3774 ! SUBROUTINE TBND
3775 ! ----------------------------------------------------------------------
3776 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
3777 ! THE MIDDLE LAYER TEMPERATURES
3778 ! ----------------------------------------------------------------------
3779       IMPLICIT NONE
3780       INTEGER, INTENT(IN)       :: NSOIL
3781       INTEGER                   :: K
3782       REAL, INTENT(IN)          :: TB, TU, ZBOT
3783       REAL, INTENT(OUT)         :: TBND1
3784       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: ZSOIL
3785       REAL                      :: ZB, ZUP
3786       REAL, PARAMETER           :: T0 = 273.15
3787
3788 ! ----------------------------------------------------------------------
3789 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
3790 ! ----------------------------------------------------------------------
3791      IF (K == 1) THEN
3792          ZUP = 0.
3793       ELSE
3794          ZUP = ZSOIL (K -1)
3795       END IF
3796 ! ----------------------------------------------------------------------
3797 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
3798 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
3799 ! ----------------------------------------------------------------------
3800       IF (K ==  NSOIL) THEN
3801          ZB = 2.* ZBOT - ZSOIL (K)
3802       ELSE
3803          ZB = ZSOIL (K +1)
3804       END IF
3805 ! ----------------------------------------------------------------------
3806 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
3807 ! ----------------------------------------------------------------------
3808
3809       TBND1 = TU + (TB - TU)* (ZUP - ZSOIL (K))/ (ZUP - ZB)
3810 ! ----------------------------------------------------------------------
3811   END SUBROUTINE TBND
3812 ! ----------------------------------------------------------------------
3813
3814
3815       SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
3816
3817 ! ----------------------------------------------------------------------
3818 ! SUBROUTINE TDFCND
3819 ! ----------------------------------------------------------------------
3820 ! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
3821 ! POINT AND TIME.
3822 ! ----------------------------------------------------------------------
3823 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
3824 ! June 2001 CHANGES: FROZEN SOIL CONDITION.
3825 ! ----------------------------------------------------------------------
3826       IMPLICIT NONE
3827       REAL, INTENT(IN)          :: QZ,  SMC, SMCMAX, SH2O
3828       REAL, INTENT(OUT)         :: DF
3829       REAL                      :: AKE, GAMMD, THKDRY, THKICE, THKO,    &
3830                                    THKQTZ,THKSAT,THKS,THKW,SATRATIO,XU, &
3831                                    XUNFROZ
3832
3833 ! ----------------------------------------------------------------------
3834 ! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
3835 !      DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
3836 !     &             0.35, 0.60, 0.40, 0.82/
3837 ! ----------------------------------------------------------------------
3838 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
3839 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
3840 ! ----------------------------------------------------------------------
3841 !  THKW ......WATER THERMAL CONDUCTIVITY
3842 !  THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
3843 !  THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
3844 !  THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
3845 !  THKICE ....ICE THERMAL CONDUCTIVITY
3846 !  SMCMAX ....POROSITY (= SMCMAX)
3847 !  QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
3848 ! ----------------------------------------------------------------------
3849 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
3850
3851 !                                  PABLO GRUNMANN, 08/17/98
3852 ! REFS.:
3853 !      FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
3854 !              AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
3855 !      JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
3856 !              UNIVERSITY OF TRONDHEIM,
3857 !      PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
3858 !              CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
3859 !              AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
3860 !              VOL. 55, PP. 1209-1224.
3861 ! ----------------------------------------------------------------------
3862 ! NEEDS PARAMETERS
3863 ! POROSITY(SOIL TYPE):
3864 !      POROS = SMCMAX
3865 ! SATURATION RATIO:
3866 ! PARAMETERS  W/(M.K)
3867       SATRATIO = SMC / SMCMAX
3868 ! ICE CONDUCTIVITY:
3869       THKICE = 2.2
3870 ! WATER CONDUCTIVITY:
3871       THKW = 0.57
3872 ! THERMAL CONDUCTIVITY OF "OTHER" SOIL COMPONENTS
3873 !      IF (QZ .LE. 0.2) THKO = 3.0
3874       THKO = 2.0
3875 ! QUARTZ' CONDUCTIVITY
3876       THKQTZ = 7.7
3877 ! SOLIDS' CONDUCTIVITY
3878       THKS = (THKQTZ ** QZ)* (THKO ** (1. - QZ))
3879
3880 ! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
3881       XUNFROZ = SH2O / SMC
3882 ! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
3883       XU = XUNFROZ * SMCMAX
3884
3885 ! SATURATED THERMAL CONDUCTIVITY
3886       THKSAT = THKS ** (1. - SMCMAX)* THKICE ** (SMCMAX - XU)* THKW **   &
3887          (XU)
3888
3889 ! DRY DENSITY IN KG/M3
3890       GAMMD = (1. - SMCMAX)*2700.
3891
3892 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
3893       THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD)
3894 ! FROZEN
3895       IF ( (SH2O + 0.0005) <  SMC ) THEN
3896          AKE = SATRATIO
3897 ! UNFROZEN
3898 ! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
3899       ELSE
3900
3901 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
3902 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
3903 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
3904
3905          IF ( SATRATIO >  0.1 ) THEN
3906
3907             AKE = LOG10 (SATRATIO) + 1.0
3908
3909 ! USE K = KDRY
3910          ELSE
3911
3912             AKE = 0.0
3913          END IF
3914 !  THERMAL CONDUCTIVITY
3915
3916       END IF
3917
3918       DF = AKE * (THKSAT - THKDRY) + THKDRY
3919 ! ----------------------------------------------------------------------
3920   END SUBROUTINE TDFCND
3921 ! ----------------------------------------------------------------------
3922
3923       SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
3924
3925 ! ----------------------------------------------------------------------
3926 ! SUBROUTINE TMPAVG
3927 ! ----------------------------------------------------------------------
3928 ! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
3929 ! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
3930 ! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
3931 ! LAYER.  TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
3932 ! ----------------------------------------------------------------------
3933       IMPLICIT NONE
3934       INTEGER  K
3935
3936       INTEGER  NSOIL
3937       REAL     DZ
3938       REAL     DZH
3939       REAL     T0
3940       REAL     TAVG
3941       REAL     TDN
3942       REAL     TM
3943       REAL     TUP
3944       REAL     X0
3945       REAL     XDN
3946       REAL     XUP
3947
3948       REAL     ZSOIL (NSOIL)
3949
3950 ! ----------------------------------------------------------------------
3951       PARAMETER (T0 = 2.7315E2)
3952       IF (K .eq. 1) THEN
3953          DZ = - ZSOIL (1)
3954       ELSE
3955          DZ = ZSOIL (K -1) - ZSOIL (K)
3956       END IF
3957
3958       DZH = DZ *0.5
3959       IF (TUP .lt. T0) THEN
3960          IF (TM .lt. T0) THEN
3961 ! ----------------------------------------------------------------------
3962 ! TUP, TM, TDN < T0
3963 ! ----------------------------------------------------------------------
3964             IF (TDN .lt. T0) THEN
3965                TAVG = (TUP + 2.0* TM + TDN)/ 4.0
3966 ! ----------------------------------------------------------------------
3967 ! TUP & TM < T0,  TDN .ge. T0
3968 ! ----------------------------------------------------------------------
3969             ELSE
3970                X0 = (T0- TM) * DZH / (TDN - TM)
3971                TAVG = 0.5 * (TUP * DZH + TM * (DZH + X0) + T0* (        &
3972      &               2.* DZH - X0)) / DZ
3973             END IF
3974          ELSE
3975 ! ----------------------------------------------------------------------
3976 ! TUP < T0, TM .ge. T0, TDN < T0
3977 ! ----------------------------------------------------------------------
3978             IF (TDN .lt. T0) THEN
3979                XUP = (T0- TUP) * DZH / (TM - TUP)
3980                XDN = DZH - (T0- TM) * DZH / (TDN - TM)
3981                TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP - XDN)       &
3982      &                + TDN * XDN) / DZ
3983 ! ----------------------------------------------------------------------
3984 ! TUP < T0, TM .ge. T0, TDN .ge. T0
3985 ! ----------------------------------------------------------------------
3986             ELSE
3987                XUP = (T0- TUP) * DZH / (TM - TUP)
3988                TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP)) / DZ
3989             END IF
3990          END IF
3991       ELSE
3992          IF (TM .lt. T0) THEN
3993 ! ----------------------------------------------------------------------
3994 ! TUP .ge. T0, TM < T0, TDN < T0
3995 ! ----------------------------------------------------------------------
3996             IF (TDN .lt. T0) THEN
3997                XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
3998                TAVG = 0.5 * (T0* (DZ - XUP) + TM * (DZH + XUP)          &
3999      &                + TDN * DZH) / DZ
4000 ! ----------------------------------------------------------------------
4001 ! TUP .ge. T0, TM < T0, TDN .ge. T0
4002 ! ----------------------------------------------------------------------
4003             ELSE
4004                XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4005                XDN = (T0- TM) * DZH / (TDN - TM)
4006                TAVG = 0.5 * (T0* (2.* DZ - XUP - XDN) + TM *            &
4007      & (XUP + XDN)) / DZ
4008             END IF
4009          ELSE
4010 ! ----------------------------------------------------------------------
4011 ! TUP .ge. T0, TM .ge. T0, TDN < T0
4012 ! ----------------------------------------------------------------------
4013             IF (TDN .lt. T0) THEN
4014                XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4015                TAVG = (T0* (DZ - XDN) +0.5* (T0+ TDN)* XDN) / DZ
4016 ! ----------------------------------------------------------------------
4017 ! TUP .ge. T0, TM .ge. T0, TDN .ge. T0
4018 ! ----------------------------------------------------------------------
4019             ELSE
4020                TAVG = (TUP + 2.0* TM + TDN) / 4.0
4021             END IF
4022          END IF
4023       END IF
4024 ! ----------------------------------------------------------------------
4025   END SUBROUTINE TMPAVG
4026 ! ----------------------------------------------------------------------
4027
4028       SUBROUTINE TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT,     &
4029      &                      CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,    &
4030      &                      RTDIS)
4031
4032 ! ----------------------------------------------------------------------
4033 ! SUBROUTINE TRANSP
4034 ! ----------------------------------------------------------------------
4035 ! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
4036 ! ----------------------------------------------------------------------
4037       IMPLICIT NONE
4038       INTEGER  I
4039       INTEGER  K
4040       INTEGER  NSOIL
4041
4042       INTEGER  NROOT
4043       REAL     CFACTR
4044       REAL     CMC
4045       REAL     CMCMAX
4046       REAL     DENOM
4047       REAL     ET (NSOIL)
4048       REAL     ETP1
4049       REAL     ETP1A
4050 !.....REAL PART(NSOIL)
4051       REAL     GX (NROOT)
4052       REAL     PC
4053       REAL     Q2
4054       REAL     RTDIS (NSOIL)
4055       REAL     RTX
4056       REAL     SFCTMP
4057       REAL     SGX
4058       REAL     SHDFAC
4059       REAL     SMC (NSOIL)
4060       REAL     SMCREF
4061       REAL     SMCWLT
4062
4063 ! ----------------------------------------------------------------------
4064 ! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
4065 ! ----------------------------------------------------------------------
4066       REAL     ZSOIL (NSOIL)
4067       DO K = 1,NSOIL
4068          ET (K) = 0.
4069 ! ----------------------------------------------------------------------
4070 ! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
4071 ! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
4072 ! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
4073 ! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
4074 ! TOTAL ETP1A.
4075 ! ----------------------------------------------------------------------
4076       END DO
4077       IF (CMC .ne. 0.0) THEN
4078          ETP1A = SHDFAC * PC * ETP1 * (1.0- (CMC / CMCMAX) ** CFACTR)
4079       ELSE
4080          ETP1A = SHDFAC * PC * ETP1
4081       END IF
4082       SGX = 0.0
4083       DO I = 1,NROOT
4084          GX (I) = ( SMC (I) - SMCWLT ) / ( SMCREF - SMCWLT )
4085          GX (I) = MAX ( MIN ( GX (I), 1. ), 0. )
4086          SGX = SGX + GX (I)
4087       END DO
4088
4089       SGX = SGX / NROOT
4090       DENOM = 0.
4091       DO I = 1,NROOT
4092          RTX = RTDIS (I) + GX (I) - SGX
4093          GX (I) = GX (I) * MAX ( RTX, 0. )
4094          DENOM = DENOM + GX (I)
4095       END DO
4096
4097       IF (DENOM .le. 0.0) DENOM = 1.
4098       DO I = 1,NROOT
4099          ET (I) = ETP1A * GX (I) / DENOM
4100 ! ----------------------------------------------------------------------
4101 ! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
4102 ! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
4103 ! ----------------------------------------------------------------------
4104 !      ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
4105 !      ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
4106 ! ----------------------------------------------------------------------
4107 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
4108 ! ----------------------------------------------------------------------
4109 !      ET(1) = RTDIS(1) * ETP1A
4110 !      ET(1) = ETP1A * PART(1)
4111 ! ----------------------------------------------------------------------
4112 ! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
4113 ! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
4114 ! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
4115 ! ----------------------------------------------------------------------
4116 !      DO K = 2,NROOT
4117 !        GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
4118 !        GX = MAX ( MIN ( GX, 1. ), 0. )
4119 ! TEST CANOPY RESISTANCE
4120 !        GX = 1.0
4121 !        ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A
4122 !        ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
4123 ! ----------------------------------------------------------------------
4124 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
4125 ! ----------------------------------------------------------------------
4126 !        ET(K) = RTDIS(K) * ETP1A
4127 !        ET(K) = ETP1A*PART(K)
4128 !      END DO
4129       END DO
4130 ! ----------------------------------------------------------------------
4131   END SUBROUTINE TRANSP
4132 ! ----------------------------------------------------------------------
4133
4134       SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT,          &
4135      &                      SICEMAX)
4136
4137 ! ----------------------------------------------------------------------
4138 ! SUBROUTINE WDFCND
4139 ! ----------------------------------------------------------------------
4140 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
4141 ! ----------------------------------------------------------------------
4142       IMPLICIT NONE
4143       REAL     BEXP
4144       REAL     DKSAT
4145       REAL     DWSAT
4146       REAL     EXPON
4147       REAL     FACTR1
4148       REAL     FACTR2
4149       REAL     SICEMAX
4150       REAL     SMC
4151       REAL     SMCMAX
4152       REAL     VKwgt
4153       REAL     WCND
4154
4155 ! ----------------------------------------------------------------------
4156 !     CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
4157 ! ----------------------------------------------------------------------
4158       REAL     WDF
4159       FACTR1 = 0.05 / SMCMAX
4160
4161 ! ----------------------------------------------------------------------
4162 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
4163 ! ----------------------------------------------------------------------
4164       FACTR2 = SMC / SMCMAX
4165       FACTR1 = MIN(FACTR1,FACTR2)
4166       EXPON = BEXP + 2.0
4167
4168 ! ----------------------------------------------------------------------
4169 ! FROZEN SOIL HYDRAULIC DIFFUSIVITY.  VERY SENSITIVE TO THE VERTICAL
4170 ! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
4171 ! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
4172 ! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
4173 ! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
4174 ! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
4175 ! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
4176 ! --
4177 ! VERSION D_10CM: ........  FACTR1 = 0.2/SMCMAX
4178 ! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
4179 ! ----------------------------------------------------------------------
4180       WDF = DWSAT * FACTR2 ** EXPON
4181       IF (SICEMAX .gt. 0.0) THEN
4182          VKWGT = 1./ (1. + (500.* SICEMAX)**3.)
4183          WDF = VKWGT * WDF + (1. - VKWGT)* DWSAT * FACTR1** EXPON
4184 ! ----------------------------------------------------------------------
4185 ! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
4186 ! ----------------------------------------------------------------------
4187       END IF
4188       EXPON = (2.0 * BEXP) + 3.0
4189       WCND = DKSAT * FACTR2 ** EXPON
4190
4191 ! ----------------------------------------------------------------------
4192   END SUBROUTINE WDFCND
4193 ! ----------------------------------------------------------------------
4194
4195       SUBROUTINE SFCDIF_off (ZLM,Z0,THZ0,THLM,SFCSPD,CZIL,AKMS,AKHS)
4196
4197 ! ----------------------------------------------------------------------
4198 ! SUBROUTINE SFCDIF (renamed SFCDIF_off to avoid clash with Eta PBL)
4199 ! ----------------------------------------------------------------------
4200 ! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
4201 ! SEE CHEN ET AL (1997, BLM)
4202 ! ----------------------------------------------------------------------
4203
4204       IMPLICIT NONE
4205       REAL     WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
4206       REAL     PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL,     &
4207      & SQVISC
4208       REAL     RIC, RRIC, FHNEU, RFC, RFAC, ZZ, PSLMU, PSLMS, PSLHU,     &
4209      & PSLHS
4210       REAL     XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
4211       REAL     SFCSPD, CZIL, AKMS, AKHS, ZILFC, ZU, ZT, RDZ, CXCH
4212       REAL     DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
4213       REAL     RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
4214 !CC   ......REAL ZTFC
4215
4216       REAL     XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN,  &
4217      &         RLMA
4218
4219       INTEGER  ITRMX, ILECH, ITR
4220       PARAMETER                                                         &
4221      &        (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40,      &
4222      &         EXCM = 0.001                                             &
4223      &        ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG          &
4224      &                  ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05,         &
4225      &                   PIHF = 3.14159265/2.)
4226       PARAMETER                                                         &
4227      &         (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8  &
4228      &         ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0                    &
4229      &          ,SQVISC = 258.2)
4230       PARAMETER                                                         &
4231      &       (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191       &
4232      &        ,RFAC = RIC / (FHNEU * RFC * RFC))
4233
4234 ! ----------------------------------------------------------------------
4235 ! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
4236 ! ----------------------------------------------------------------------
4237 ! LECH'S SURFACE FUNCTIONS
4238 ! ----------------------------------------------------------------------
4239       PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
4240       PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
4241       PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
4242
4243 ! ----------------------------------------------------------------------
4244 ! PAULSON'S SURFACE FUNCTIONS
4245 ! ----------------------------------------------------------------------
4246       PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
4247       PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5)   &
4248      &        +2.* ATAN (XX)                                            &
4249      &- PIHF
4250       PSPMS (YY)= 5.* YY
4251       PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
4252
4253 ! ----------------------------------------------------------------------
4254 ! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
4255 ! OVER SOLID SURFACE (LAND, SEA-ICE).
4256 ! ----------------------------------------------------------------------
4257       PSPHS (YY)= 5.* YY
4258
4259 ! ----------------------------------------------------------------------
4260 !     ZTFC: RATIO OF ZOH/ZOM  LESS OR EQUAL THAN 1
4261 !     C......ZTFC=0.1
4262 !     CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
4263 ! ----------------------------------------------------------------------
4264       ILECH = 0
4265
4266 ! ----------------------------------------------------------------------
4267       ZILFC = - CZIL * VKRM * SQVISC
4268 !     C.......ZT=Z0*ZTFC
4269       ZU = Z0
4270       RDZ = 1./ ZLM
4271       CXCH = EXCM * RDZ
4272       DTHV = THLM - THZ0
4273
4274 ! ----------------------------------------------------------------------
4275 ! BELJARS CORRECTION OF USTAR
4276 ! ----------------------------------------------------------------------
4277       DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
4278 !cc   If statements to avoid TANGENT LINEAR problems near zero
4279       BTGH = BTG * HPBL
4280       IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
4281          WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
4282       ELSE
4283          WSTAR2 = 0.0
4284       END IF
4285
4286 ! ----------------------------------------------------------------------
4287 ! ZILITINKEVITCH APPROACH FOR ZT
4288 ! ----------------------------------------------------------------------
4289       USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
4290
4291 ! ----------------------------------------------------------------------
4292       ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
4293       ZSLU = ZLM + ZU
4294 !     PRINT*,'ZSLT=',ZSLT
4295 !     PRINT*,'ZLM=',ZLM
4296 !     PRINT*,'ZT=',ZT
4297
4298       ZSLT = ZLM + ZT
4299       RLOGU = log (ZSLU / ZU)
4300
4301       RLOGT = log (ZSLT / ZT)
4302 !     PRINT*,'RLMO=',RLMO
4303 !     PRINT*,'ELFC=',ELFC
4304 !     PRINT*,'AKHS=',AKHS
4305 !     PRINT*,'DTHV=',DTHV
4306 !     PRINT*,'USTAR=',USTAR
4307
4308       RLMO = ELFC * AKHS * DTHV / USTAR **3
4309 ! ----------------------------------------------------------------------
4310 ! 1./MONIN-OBUKKHOV LENGTH-SCALE
4311 ! ----------------------------------------------------------------------
4312       DO ITR = 1,ITRMX
4313          ZETALT = MAX (ZSLT * RLMO,ZTMIN)
4314          RLMO = ZETALT / ZSLT
4315          ZETALU = ZSLU * RLMO
4316          ZETAU = ZU * RLMO
4317
4318          ZETAT = ZT * RLMO
4319          IF (ILECH .eq. 0) THEN
4320             IF (RLMO .lt. 0.)THEN
4321                XLU4 = 1. -16.* ZETALU
4322                XLT4 = 1. -16.* ZETALT
4323                XU4 = 1. -16.* ZETAU
4324
4325                XT4 = 1. -16.* ZETAT
4326                XLU = SQRT (SQRT (XLU4))
4327                XLT = SQRT (SQRT (XLT4))
4328                XU = SQRT (SQRT (XU4))
4329
4330                XT = SQRT (SQRT (XT4))
4331 !     PRINT*,'-----------1------------'
4332 !     PRINT*,'PSMZ=',PSMZ
4333 !     PRINT*,'PSPMU(ZETAU)=',PSPMU(ZETAU)
4334 !     PRINT*,'XU=',XU
4335 !     PRINT*,'------------------------'
4336                PSMZ = PSPMU (XU)
4337                SIMM = PSPMU (XLU) - PSMZ + RLOGU
4338                PSHZ = PSPHU (XT)
4339                SIMH = PSPHU (XLT) - PSHZ + RLOGT
4340             ELSE
4341                ZETALU = MIN (ZETALU,ZTMAX)
4342                ZETALT = MIN (ZETALT,ZTMAX)
4343 !     PRINT*,'-----------2------------'
4344 !     PRINT*,'PSMZ=',PSMZ
4345 !     PRINT*,'PSPMS(ZETAU)=',PSPMS(ZETAU)
4346 !     PRINT*,'ZETAU=',ZETAU
4347 !     PRINT*,'------------------------'
4348                PSMZ = PSPMS (ZETAU)
4349                SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
4350                PSHZ = PSPHS (ZETAT)
4351                SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
4352             END IF
4353 ! ----------------------------------------------------------------------
4354 ! LECH'S FUNCTIONS
4355 ! ----------------------------------------------------------------------
4356          ELSE
4357             IF (RLMO .lt. 0.)THEN
4358 !     PRINT*,'-----------3------------'
4359 !     PRINT*,'PSMZ=',PSMZ
4360 !     PRINT*,'PSLMU(ZETAU)=',PSLMU(ZETAU)
4361 !     PRINT*,'ZETAU=',ZETAU
4362 !     PRINT*,'------------------------'
4363                PSMZ = PSLMU (ZETAU)
4364                SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
4365                PSHZ = PSLHU (ZETAT)
4366                SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
4367             ELSE
4368                ZETALU = MIN (ZETALU,ZTMAX)
4369
4370                ZETALT = MIN (ZETALT,ZTMAX)
4371 !     PRINT*,'-----------4------------'
4372 !     PRINT*,'PSMZ=',PSMZ
4373 !     PRINT*,'PSLMS(ZETAU)=',PSLMS(ZETAU)
4374 !     PRINT*,'ZETAU=',ZETAU
4375 !     PRINT*,'------------------------'
4376                PSMZ = PSLMS (ZETAU)
4377                SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
4378                PSHZ = PSLHS (ZETAT)
4379                SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
4380             END IF
4381 ! ----------------------------------------------------------------------
4382 ! BELJAARS CORRECTION FOR USTAR
4383 ! ----------------------------------------------------------------------
4384          END IF
4385
4386 ! ----------------------------------------------------------------------
4387 ! ZILITINKEVITCH FIX FOR ZT
4388 ! ----------------------------------------------------------------------
4389          USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
4390
4391          ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
4392          ZSLT = ZLM + ZT
4393 !-----------------------------------------------------------------------
4394          RLOGT = log (ZSLT / ZT)
4395          USTARK = USTAR * VKRM
4396          AKMS = MAX (USTARK / SIMM,CXCH)
4397 !-----------------------------------------------------------------------
4398 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
4399 !-----------------------------------------------------------------------
4400          AKHS = MAX (USTARK / SIMH,CXCH)
4401          IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
4402             WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
4403          ELSE
4404             WSTAR2 = 0.0
4405          END IF
4406 !-----------------------------------------------------------------------
4407          RLMN = ELFC * AKHS * DTHV / USTAR **3
4408 !-----------------------------------------------------------------------
4409 !     IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT)    GO TO 110
4410 !-----------------------------------------------------------------------
4411          RLMA = RLMO * WOLD+ RLMN * WNEW
4412 !-----------------------------------------------------------------------
4413          RLMO = RLMA
4414 !     PRINT*,'----------------------------'
4415 !     PRINT*,'SFCDIF OUTPUT !  ! ! ! ! ! ! ! !  !   !    !'
4416
4417 !     PRINT*,'ZLM=',ZLM
4418 !     PRINT*,'Z0=',Z0
4419 !     PRINT*,'THZ0=',THZ0
4420 !     PRINT*,'THLM=',THLM
4421 !     PRINT*,'SFCSPD=',SFCSPD
4422 !     PRINT*,'CZIL=',CZIL
4423 !     PRINT*,'AKMS=',AKMS
4424 !     PRINT*,'AKHS=',AKHS
4425 !     PRINT*,'----------------------------'
4426
4427       END DO
4428 ! ----------------------------------------------------------------------
4429   END SUBROUTINE SFCDIF_off
4430 ! ----------------------------------------------------------------------
4431
4432 END MODULE module_sf_noahlsm