phys/module_sf_noahlsm.F

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