wrfv2_fire/phys/module_sf_noahlsm.F

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