phys/module_sf_noahlsm_glacial_only.F

   1 MODULE module_sf_noahlsm_glacial_only
   2 #if defined(mpas)
   3 use mpas_atmphys_constants
   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
   8 use module_wrf_error
   9 #define FATAL_ERROR(M) call wrf_error_fatal( M )
  10 #endif
  11
  12   USE module_sf_noahlsm, ONLY : RD, SIGMA, CPH2O, CPICE, LSUBF, EMISSI_S, ROSR12
  13   USE module_sf_noahlsm, ONLY : LVCOEF_DATA
  14
  15   PRIVATE :: ALCALC
  16   PRIVATE :: CSNOW
  17   PRIVATE :: HRTICE
  18   PRIVATE :: HSTEP
  19   PRIVATE :: PENMAN
  20   PRIVATE :: SHFLX
  21   PRIVATE :: SNOPAC
  22   PRIVATE :: SNOWPACK
  23   PRIVATE :: SNOWZ0
  24   PRIVATE :: SNOW_NEW
  25
  26   integer, private :: iloc, jloc
  27 !$omp threadprivate(iloc, jloc)
  28
  29 CONTAINS
  30
  31   SUBROUTINE SFLX_GLACIAL (IILOC,JJLOC,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH,     &    !C
  32        &                   LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2,           &    !F
  33        &                   TH2,Q2SAT,DQSDT2,                            &    !I
  34        &                   ALB, SNOALB,TBOT, Z0BRD, Z0, EMISSI, EMBRD,  &    !S
  35        &                   T1,STC,SNOWH,SNEQV,ALBEDO,CH,                &    !H
  36 ! ----------------------------------------------------------------------
  37 ! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
  38 ! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
  39 ! MODEL).  OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
  40 ! ----------------------------------------------------------------------
  41        &                   ETA,SHEAT, ETA_KINEMATIC,FDOWN,              &    !O
  42        &                   ESNOW,DEW,                                   &    !O
  43        &                   ETP,SSOIL,                                   &    !O
  44        &                   FLX1,FLX2,FLX3,                              &    !O
  45        &                   SNOMLT,SNCOVR,                               &    !O
  46        &                   RUNOFF1,                                     &    !O
  47        &                   Q1,                                          &    !D
  48        &                   SNOTIME1,                                    &
  49        &                   RIBB)
  50 ! ----------------------------------------------------------------------
  51 ! SUB-DRIVER FOR "Noah LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
  52 ! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE ICE TEMPERATURE, SKIN
  53 ! TEMPERATURE, SNOWPACK WATER CONTENT, SNOWDEPTH, AND ALL TERMS OF THE
  54 ! SURFACE ENERGY BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF
  55 ! DOWNWARD RADIATION AND PRECIP)
  56 ! ----------------------------------------------------------------------
  57 ! SFLX ARGUMENT LIST KEY:
  58 ! ----------------------------------------------------------------------
  59 !  C  CONFIGURATION INFORMATION
  60 !  F  FORCING DATA
  61 !  I  OTHER (INPUT) FORCING DATA
  62 !  S  SURFACE CHARACTERISTICS
  63 !  H  HISTORY (STATE) VARIABLES
  64 !  O  OUTPUT VARIABLES
  65 !  D  DIAGNOSTIC OUTPUT
  66 ! ----------------------------------------------------------------------
  67 ! 1. CONFIGURATION INFORMATION (C):
  68 ! ----------------------------------------------------------------------
  69 !   DT         TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
  70 !                1800 SECS OR LESS)
  71 !   ZLVL       HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
  72 !   NSOIL      NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
  73 !                PARAMETER NSOLD SET BELOW)
  74 !   SLDPTH     THE THICKNESS OF EACH SOIL LAYER (M)
  75 ! ----------------------------------------------------------------------
  76 ! 3. FORCING DATA (F):
  77 ! ----------------------------------------------------------------------
  78 !   LWDN       LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
  79 !   SOLNET     NET DOWNWARD SOLAR RADIATION ((W M-2; POSITIVE)
  80 !   SFCPRS     PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
  81 !   PRCP       PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
  82 !   SFCTMP     AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
  83 !   TH2        AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
  84 !   Q2         MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
  85 !   FFROZP     FRACTION OF FROZEN PRECIPITATION
  86 ! ----------------------------------------------------------------------
  87 ! 4. OTHER FORCING (INPUT) DATA (I):
  88 ! ----------------------------------------------------------------------
  89 !   Q2SAT      SAT SPECIFIC HUMIDITY AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
  90 !   DQSDT2     SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
  91 !                (KG KG-1 K-1)
  92 ! ----------------------------------------------------------------------
  93 ! 5. CANOPY/SOIL CHARACTERISTICS (S):
  94 ! ----------------------------------------------------------------------
  95 !   ALB        BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
  96 !                DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
  97 !                MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
  98 !                INCLUDE DIURNAL SUN ANGLE EFFECT)
  99 !   SNOALB     UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
 100 !                ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
 101 !   TBOT       BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
 102 !                TEMPERATURE)
 103 !   Z0BRD      Background fixed roughness length (M)
 104 !   Z0         Time varying roughness length (M) as function of snow depth
 105 !   EMBRD      Background surface emissivity (between 0 and 1)
 106 !   EMISSI     Surface emissivity (between 0 and 1)
 107 ! ----------------------------------------------------------------------
 108 ! 6. HISTORY (STATE) VARIABLES (H):
 109 ! ----------------------------------------------------------------------
 110 !  T1          GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
 111 !  STC(NSOIL)  SOIL TEMP (K)
 112 !  SNOWH       ACTUAL SNOW DEPTH (M)
 113 !  SNEQV       LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
 114 !                NOTE: SNOW DENSITY = SNEQV/SNOWH
 115 !  ALBEDO      SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
 116 !                =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
 117 !                =FCT(MSNOALB,ALB,SHDFAC,SHDMIN) WHEN SNEQV>0
 118 !  CH          SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
 119 !                (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
 120 !                IT HAS BEEN MULTIPLIED BY WIND SPEED.
 121 ! ----------------------------------------------------------------------
 122 ! 7. OUTPUT (O):
 123 ! ----------------------------------------------------------------------
 124 ! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
 125 ! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL.  FOR THIS APPLICATION,
 126 ! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
 127 ! NECESSARY.  OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
 128 !   ETA        ACTUAL LATENT HEAT FLUX (W m-2: NEGATIVE, IF UP FROM
 129 !              SURFACE)
 130 !  ETA_KINEMATIC atctual latent heat flux in Kg m-2 s-1
 131 !   SHEAT      SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
 132 !              SURFACE)
 133 !   FDOWN      Radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
 134 ! ----------------------------------------------------------------------
 135 !   ESNOW      SUBLIMATION FROM (OR DEPOSITION TO IF <0) SNOWPACK
 136 !                (W m-2)
 137 !   DEW        DEWFALL (OR FROSTFALL FOR T<273.15) (M)
 138 ! ----------------------------------------------------------------------
 139 !   ETP        POTENTIAL EVAPORATION (W m-2)
 140 !   SSOIL      SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
 141 ! ----------------------------------------------------------------------
 142 !   FLX1       PRECIP-SNOW SFC (W M-2)
 143 !   FLX2       FREEZING RAIN LATENT HEAT FLUX (W M-2)
 144 !   FLX3       PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
 145 ! ----------------------------------------------------------------------
 146 !   SNOMLT     SNOW MELT (M) (WATER EQUIVALENT)
 147 !   SNCOVR     FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
 148 ! ----------------------------------------------------------------------
 149 !   RUNOFF1    SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
 150 ! ----------------------------------------------------------------------
 151 ! 8. DIAGNOSTIC OUTPUT (D):
 152 ! ----------------------------------------------------------------------
 153 !   Q1         Effective mixing ratio at surface (kg kg-1), used for
 154 !              diagnosing the mixing ratio at 2 meter for coupled model
 155 !  Documentation for SNOTIME1 and SNOABL2 ?????
 156 !  What categories of arguments do these variables fall into ????
 157 !  Documentation for RIBB ?????
 158 !  What category of argument does RIBB fall into ?????
 159 ! ----------------------------------------------------------------------
 160
 161     IMPLICIT NONE
 162 ! ----------------------------------------------------------------------
 163     integer, intent(in) ::  iiloc, jjloc
 164     INTEGER, INTENT(IN) ::  ISICE
 165 ! ----------------------------------------------------------------------
 166     LOGICAL             ::  FRZGRA, SNOWNG
 167
 168 ! ----------------------------------------------------------------------
 169 ! 1. CONFIGURATION INFORMATION (C):
 170 ! ----------------------------------------------------------------------
 171     INTEGER, INTENT(IN) ::  NSOIL
 172     INTEGER             ::  KZ
 173
 174 ! ----------------------------------------------------------------------
 175 ! 2. LOGICAL:
 176 ! ----------------------------------------------------------------------
 177
 178     REAL, INTENT(IN)   :: DT,DQSDT2,LWDN,PRCP,     &
 179          &                Q2,Q2SAT,SFCPRS,SFCTMP, SNOALB,          &
 180          &                SOLNET,TBOT,TH2,ZLVL,FFROZP
 181     REAL, INTENT(OUT)  :: EMBRD, ALBEDO
 182     REAL, INTENT(INOUT):: CH,SNEQV,SNCOVR,SNOWH,T1,Z0BRD,EMISSI,ALB
 183     REAL, INTENT(INOUT):: SNOTIME1
 184     REAL, INTENT(INOUT):: RIBB
 185     REAL, DIMENSION(1:NSOIL), INTENT(IN)    :: SLDPTH
 186     REAL, DIMENSION(1:NSOIL), INTENT(INOUT) ::  STC
 187     REAL, DIMENSION(1:NSOIL) ::   ZSOIL
 188
 189     REAL,INTENT(OUT)   :: ETA_KINEMATIC,DEW,ESNOW,ETA,  &
 190          &                ETP,FLX1,FLX2,FLX3,SHEAT,RUNOFF1,    &
 191          &                SSOIL,SNOMLT,FDOWN,Q1
 192     REAL               :: DF1,DSOIL,DTOT,FRCSNO,FRCSOI,          &
 193          &                PRCP1,RCH,RR,RSNOW,SNDENS,SNCOND,SN_NEW,     &
 194          &                T1V,T24,T2V,TH2V,TSNOW,Z0,PRCPF,RHO
 195
 196 ! ----------------------------------------------------------------------
 197 ! DECLARATIONS - PARAMETERS
 198 ! ----------------------------------------------------------------------
 199     REAL, PARAMETER :: TFREEZ = 273.15
 200     REAL, PARAMETER :: LVH2O = 2.501E+6
 201     REAL, PARAMETER :: LSUBS = 2.83E+6
 202     REAL, PARAMETER :: R = 287.04
 203
 204 ! ----------------------------------------------------------------------
 205     iloc = iiloc
 206     jloc = jjloc
 207 ! ----------------------------------------------------------------------
 208     ZSOIL (1) = - SLDPTH (1)
 209     DO KZ = 2,NSOIL
 210        ZSOIL (KZ) = - SLDPTH (KZ) + ZSOIL (KZ -1)
 211     END DO
 212
 213 ! ----------------------------------------------------------------------
 214 ! IF S.W.E. (SNEQV) BELOW THRESHOLD LOWER BOUND (0.10 M FOR GLACIAL
 215 ! ICE), THEN SET AT LOWER BOUND
 216 ! ----------------------------------------------------------------------
 217     IF ( SNEQV < 0.10 ) THEN
 218        SNEQV = 0.10
 219        SNOWH = 0.50
 220     ENDIF
 221 ! ----------------------------------------------------------------------
 222 ! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
 223 ! SNOW THERMAL CONDUCTIVITY "SNCOND"
 224 ! ----------------------------------------------------------------------
 225     SNDENS = SNEQV / SNOWH
 226     IF(SNDENS > 1.0) THEN
 227        FATAL_ERROR( 'Physical snow depth is less than snow water equiv.' )
 228     ENDIF
 229
 230     CALL CSNOW (SNCOND,SNDENS)
 231 ! ----------------------------------------------------------------------
 232 ! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
 233 ! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
 234 ! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
 235 ! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
 236 ! ----------------------------------------------------------------------
 237
 238     SNOWNG = .FALSE.
 239     FRZGRA = .FALSE.
 240     IF (PRCP > 0.0) THEN
 241 ! ----------------------------------------------------------------------
 242 ! Snow defined when fraction of frozen precip (FFROZP) > 0.5,
 243 ! passed in from model microphysics.
 244 ! ----------------------------------------------------------------------
 245        IF (FFROZP .GT. 0.5) THEN
 246           SNOWNG = .TRUE.
 247        ELSE
 248           IF (T1 <= TFREEZ) FRZGRA = .TRUE.
 249        END IF
 250     END IF
 251 ! ----------------------------------------------------------------------
 252 ! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
 253 ! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
 254 ! IT TO THE EXISTING SNOWPACK.
 255 ! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
 256 ! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
 257 ! ----------------------------------------------------------------------
 258     IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
 259        SN_NEW = PRCP * DT * 0.001
 260        SNEQV = SNEQV + SN_NEW
 261        PRCPF = 0.0
 262
 263 ! ----------------------------------------------------------------------
 264 ! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
 265 ! UPDATE SNOW THERMAL CONDUCTIVITY
 266 ! ----------------------------------------------------------------------
 267        CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
 268
 269 ! ----------------------------------------------------------------------
 270 ! kmh 09/04/2006 set Snow Density at 0.2 g/cm**3
 271 ! for "cold permanent ice" or new "dry" snow
 272 ! if soil temperature less than 268.15 K, treat as typical
 273 ! Antarctic/Greenland snow firn
 274 ! ----------------------------------------------------------------------
 275        IF ( SNCOVR .GT. 0.99 ) THEN
 276           IF ( STC(1) .LT. (TFREEZ - 5.) ) SNDENS = 0.2
 277           IF ( SNOWNG .AND. (T1.LT.273.) .AND. (SFCTMP.LT.273.) ) SNDENS=0.2
 278        ENDIF
 279
 280        CALL CSNOW (SNCOND,SNDENS)
 281
 282 ! ----------------------------------------------------------------------
 283 ! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
 284 ! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL
 285 ! ----------------------------------------------------------------------
 286     ELSE
 287        PRCPF = PRCP
 288     ENDIF
 289
 290 ! ----------------------------------------------------------------------
 291 ! DETERMINE SNOW FRACTIONAL COVERAGE.
 292 !      KWM:  Set SNCOVR to 1.0 because SNUP is set small in VEGPARM.TBL,
 293 !      and SNEQV is at least 0.1 (as set above)
 294 ! ----------------------------------------------------------------------
 295     SNCOVR = 1.0
 296
 297 ! ----------------------------------------------------------------------
 298 ! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
 299 ! ----------------------------------------------------------------------
 300
 301     CALL ALCALC (ALB,SNOALB,EMBRD,T1,ALBEDO,EMISSI,   &
 302          &       DT,SNOWNG,SNOTIME1)
 303
 304 ! ----------------------------------------------------------------------
 305 ! THERMAL CONDUCTIVITY
 306 ! ----------------------------------------------------------------------
 307     DF1 = SNCOND
 308
 309     DSOIL = - (0.5 * ZSOIL (1))
 310     DTOT = SNOWH + DSOIL
 311     FRCSNO = SNOWH / DTOT
 312
 313 ! 1. HARMONIC MEAN (SERIES FLOW)
 314 !        DF1 = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
 315     FRCSOI = DSOIL / DTOT
 316
 317 ! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
 318 !        DF1 = (SNCOND**FRCSNO)*(DF1**FRCSOI)
 319     DF1 = FRCSNO * SNCOND + FRCSOI * DF1
 320
 321 ! ----------------------------------------------------------------------
 322 ! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
 323 ! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
 324 ! MID-LAYER SOIL TEMPERATURE
 325 ! ----------------------------------------------------------------------
 326     IF ( DTOT .GT. 2.*DSOIL ) then
 327        DTOT = 2.*DSOIL
 328     ENDIF
 329     SSOIL = DF1 * ( T1 - STC(1) ) / DTOT
 330
 331 ! ----------------------------------------------------------------------
 332 ! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
 333 ! THE PREVIOUS TIMESTEP.
 334 ! ----------------------------------------------------------------------
 335
 336     CALL SNOWZ0 (Z0,Z0BRD,SNOWH)
 337
 338 ! ----------------------------------------------------------------------
 339 ! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
 340 ! PENMAN EP SUBROUTINE THAT FOLLOWS
 341 ! ----------------------------------------------------------------------
 342
 343     FDOWN = SOLNET + LWDN
 344
 345 ! ----------------------------------------------------------------------
 346 ! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
 347 ! PENMAN.
 348 ! ----------------------------------------------------------------------
 349
 350     T2V = SFCTMP * (1.0+ 0.61 * Q2 )
 351     RHO = SFCPRS / (RD * T2V)
 352     RCH = RHO * 1004.6 * CH
 353     T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
 354
 355 ! ----------------------------------------------------------------------
 356 ! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
 357 ! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
 358 ! CALCULATIONS.
 359 ! ----------------------------------------------------------------------
 360
 361     ! PENMAN returns ETP, FLX2, and RR
 362     CALL PENMAN (SFCTMP,SFCPRS,CH,TH2,PRCP,FDOWN,T24,SSOIL,     &
 363          &       Q2,Q2SAT,ETP,RCH,RR,SNOWNG,FRZGRA,             &
 364          &       DQSDT2,FLX2,EMISSI,T1)
 365
 366     CALL SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,NSOIL,DT,DF1,        &
 367          &       Q2,T1,SFCTMP,T24,TH2,FDOWN,SSOIL,STC,          &
 368          &       SFCPRS,RCH,RR,SNEQV,SNDENS,SNOWH,ZSOIL,TBOT,   &
 369          &       SNOMLT,DEW,FLX1,FLX2,FLX3,ESNOW,EMISSI,RIBB)
 370
 371 !   ETA_KINEMATIC =  ESNOW
 372     ETA_KINEMATIC =  ETP
 373
 374 ! ----------------------------------------------------------------------
 375 ! Effective mixing ratio at grnd level (skin)
 376 ! ----------------------------------------------------------------------
 377     Q1=Q2+ETA_KINEMATIC*CP/RCH
 378
 379 ! ----------------------------------------------------------------------
 380 ! DETERMINE SENSIBLE HEAT (H) IN ENERGY UNITS (W M-2)
 381 ! ----------------------------------------------------------------------
 382     SHEAT = - (CH * CP * SFCPRS)/ (R * T2V) * ( TH2- T1 )
 383
 384 ! ----------------------------------------------------------------------
 385 ! CONVERT EVAP TERMS FROM KINEMATIC (KG M-2 S-1) TO ENERGY UNITS (W M-2)
 386 ! ----------------------------------------------------------------------
 387     ESNOW = ESNOW * LSUBS
 388     ETP   = ETP   * LSUBS
 389     IF (ETP .GT. 0.) THEN
 390        ETA = ESNOW
 391     ELSE
 392        ETA = ETP
 393     ENDIF
 394
 395 ! ----------------------------------------------------------------------
 396 ! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
 397 !   SSOIL>0: WARM THE SURFACE  (NIGHT TIME)
 398 !   SSOIL<0: COOL THE SURFACE  (DAY TIME)
 399 ! ----------------------------------------------------------------------
 400     SSOIL = -1.0* SSOIL
 401
 402 ! ----------------------------------------------------------------------
 403 ! FOR THE CASE OF GLACIAL-ICE, ADD ANY SNOWMELT DIRECTLY TO SURFACE
 404 ! RUNOFF (RUNOFF1) SINCE THERE IS NO SOIL MEDIUM
 405 ! ----------------------------------------------------------------------
 406     RUNOFF1 = SNOMLT / DT
 407
 408 ! ----------------------------------------------------------------------
 409   END SUBROUTINE SFLX_GLACIAL
 410 ! ----------------------------------------------------------------------
 411
 412   SUBROUTINE ALCALC (ALB,SNOALB,EMBRD,TSNOW,ALBEDO,EMISSI,   &
 413        &             DT,SNOWNG,SNOTIME1)
 414
 415 ! ----------------------------------------------------------------------
 416 ! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
 417 !   ALB     SNOWFREE ALBEDO
 418 !   SNOALB  MAXIMUM (DEEP) SNOW ALBEDO
 419 !   ALBEDO  SURFACE ALBEDO INCLUDING SNOW EFFECT
 420 !   TSNOW   SNOW SURFACE TEMPERATURE (K)
 421 ! ----------------------------------------------------------------------
 422     IMPLICIT NONE
 423
 424 ! ----------------------------------------------------------------------
 425 ! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
 426 ! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
 427 ! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
 428 ! (1985, JCAM, VOL 24, 402-411)
 429 ! ----------------------------------------------------------------------
 430     REAL,    INTENT(IN)    :: ALB, SNOALB, EMBRD, TSNOW
 431     REAL,    INTENT(IN)    :: DT
 432     LOGICAL, INTENT(IN)    :: SNOWNG
 433     REAL,    INTENT(INOUT) :: SNOTIME1
 434     REAL,    INTENT(OUT)   :: ALBEDO, EMISSI
 435     REAL                   :: SNOALB2
 436     REAL                   :: TM,SNOALB1
 437     REAL,    PARAMETER     :: SNACCA=0.94,SNACCB=0.58,SNTHWA=0.82,SNTHWB=0.46
 438 ! turn off vegetation effect
 439 !      ALBEDO = ALB + (1.0- (SHDFAC - SHDMIN))* SNCOVR * (SNOALB - ALB)
 440 !      ALBEDO = (1.0-SNCOVR)*ALB + SNCOVR*SNOALB !this is equivalent to below
 441     ALBEDO = ALB + (SNOALB-ALB)
 442     EMISSI = EMBRD + (EMISSI_S - EMBRD)
 443
 444 !     BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
 445 !          IF (TSNOW.LE.263.16) THEN
 446 !            ALBEDO=SNOALB
 447 !          ELSE
 448 !            IF (TSNOW.LT.273.16) THEN
 449 !              TM=0.1*(TSNOW-263.16)
 450 !              SNOALB1=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
 451 !            ELSE
 452 !              SNOALB1=0.67
 453 !             IF(SNCOVR.GT.0.95) SNOALB1= 0.6
 454 !             SNOALB1 = ALB + SNCOVR*(SNOALB-ALB)
 455 !            ENDIF
 456 !          ENDIF
 457 !            ALBEDO = ALB + SNCOVR*(SNOALB1-ALB)
 458
 459 !     ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
 460 !          SNOALB1 = SNOALB+COEF*(0.85-SNOALB)
 461 !          SNOALB2=SNOALB1
 462 !!m          LSTSNW=LSTSNW+1
 463 !          SNOTIME1 = SNOTIME1 + DT
 464 !          IF (SNOWNG) THEN
 465 !             SNOALB2=SNOALB
 466 !!m             LSTSNW=0
 467 !             SNOTIME1 = 0.0
 468 !          ELSE
 469 !            IF (TSNOW.LT.273.16) THEN
 470 !!              SNOALB2=SNOALB-0.008*LSTSNW*DT/86400
 471 !!m              SNOALB2=SNOALB-0.008*SNOTIME1/86400
 472 !              SNOALB2=(SNOALB2-0.65)*EXP(-0.05*DT/3600)+0.65
 473 !!              SNOALB2=(ALBEDO-0.65)*EXP(-0.01*DT/3600)+0.65
 474 !            ELSE
 475 !              SNOALB2=(SNOALB2-0.5)*EXP(-0.0005*DT/3600)+0.5
 476 !!              SNOALB2=(SNOALB-0.5)*EXP(-0.24*LSTSNW*DT/86400)+0.5
 477 !!m              SNOALB2=(SNOALB-0.5)*EXP(-0.24*SNOTIME1/86400)+0.5
 478 !            ENDIF
 479 !          ENDIF
 480 !
 481 !!               print*,'SNOALB2',SNOALB2,'ALBEDO',ALBEDO,'DT',DT
 482 !          ALBEDO = ALB + SNCOVR*(SNOALB2-ALB)
 483 !          IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
 484 !!m          LSTSNW1=LSTSNW
 485 !!          SNOTIME = SNOTIME1
 486
 487 ! formulation by Livneh
 488 ! ----------------------------------------------------------------------
 489 ! SNOALB IS CONSIDERED AS THE MAXIMUM SNOW ALBEDO FOR NEW SNOW, AT
 490 ! A VALUE OF 85%. SNOW ALBEDO CURVE DEFAULTS ARE FROM BRAS P.263. SHOULD
 491 ! NOT BE CHANGED EXCEPT FOR SERIOUS PROBLEMS WITH SNOW MELT.
 492 ! TO IMPLEMENT ACCUMULATIN PARAMETERS, SNACCA AND SNACCB, ASSERT THAT IT
 493 ! IS INDEED ACCUMULATION SEASON. I.E. THAT SNOW SURFACE TEMP IS BELOW
 494 ! ZERO AND THE DATE FALLS BETWEEN OCTOBER AND FEBRUARY
 495 ! ----------------------------------------------------------------------
 496     SNOALB1 = SNOALB+LVCOEF_DATA*(0.85-SNOALB)
 497     SNOALB2=SNOALB1
 498 ! ---------------- Initial LSTSNW --------------------------------------
 499     IF (SNOWNG) THEN
 500        SNOTIME1 = 0.
 501     ELSE
 502        SNOTIME1=SNOTIME1+DT
 503 !               IF (TSNOW.LT.273.16) THEN
 504        SNOALB2=SNOALB1*(SNACCA**((SNOTIME1/86400.0)**SNACCB))
 505 !               ELSE
 506 !                  SNOALB2 =SNOALB1*(SNTHWA**((SNOTIME1/86400.0)**SNTHWB))
 507 !               ENDIF
 508     ENDIF
 509
 510     SNOALB2 = MAX ( SNOALB2, ALB )
 511     ALBEDO = ALB + (SNOALB2-ALB)
 512     IF (ALBEDO .GT. SNOALB2) ALBEDO=SNOALB2
 513
 514 !          IF (TSNOW.LT.273.16) THEN
 515 !            ALBEDO=SNOALB-0.008*DT/86400
 516 !          ELSE
 517 !            ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
 518 !          ENDIF
 519
 520 !      IF (ALBEDO > SNOALB) ALBEDO = SNOALB
 521
 522 ! ----------------------------------------------------------------------
 523   END SUBROUTINE ALCALC
 524 ! ----------------------------------------------------------------------
 525
 526   SUBROUTINE CSNOW (SNCOND,DSNOW)
 527
 528 ! ----------------------------------------------------------------------
 529 ! CALCULATE SNOW TERMAL CONDUCTIVITY
 530 ! ----------------------------------------------------------------------
 531     IMPLICIT NONE
 532     REAL, INTENT(IN)  :: DSNOW
 533     REAL, INTENT(OUT) :: SNCOND
 534     REAL              :: C
 535     REAL, PARAMETER   :: UNIT = 0.11631
 536
 537 ! ----------------------------------------------------------------------
 538 ! SNCOND IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
 539 ! CSNOW IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
 540 ! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
 541 ! ----------------------------------------------------------------------
 542     C = 0.328*10** (2.25* DSNOW)
 543 !      CSNOW=UNIT*C
 544
 545 ! ----------------------------------------------------------------------
 546 ! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
 547 ! ----------------------------------------------------------------------
 548 !      SNCOND=0.0293*(1.+100.*DSNOW**2)
 549 !      CSNOW=0.0293*(1.+100.*DSNOW**2)
 550
 551 ! ----------------------------------------------------------------------
 552 ! E. ANDERSEN FROM FLERCHINGER
 553 ! ----------------------------------------------------------------------
 554 !      SNCOND=0.021+2.51*DSNOW**2
 555 !      CSNOW=0.021+2.51*DSNOW**2
 556
 557 !      SNCOND = UNIT * C
 558 ! double snow thermal conductivity
 559     SNCOND = 2.0 * UNIT * C
 560
 561 ! ----------------------------------------------------------------------
 562   END SUBROUTINE CSNOW
 563 ! ----------------------------------------------------------------------
 564
 565   SUBROUTINE HRTICE (RHSTS,STC,TBOT,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
 566
 567 ! ----------------------------------------------------------------------
 568 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
 569 ! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE (ICE=1) OR GLACIAL
 570 ! ICE (ICE=-1). COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE
 571 ! TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
 572 !
 573 ! (NOTE:  THIS SUBROUTINE ONLY CALLED FOR SEA-ICE OR GLACIAL ICE, BUT
 574 ! NOT FOR NON-GLACIAL LAND (ICE = 0).
 575 ! ----------------------------------------------------------------------
 576     IMPLICIT NONE
 577
 578
 579     INTEGER,                  INTENT(IN)  :: NSOIL
 580     REAL,                     INTENT(IN)  :: DF1,YY,ZZ1
 581     REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI,CI
 582     REAL, DIMENSION(1:NSOIL), INTENT(IN)  :: STC, ZSOIL
 583     REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTS
 584     REAL,                     INTENT(IN)  :: TBOT
 585     INTEGER                :: K
 586     REAL                   :: DDZ,DDZ2,DENOM,DTSDZ,DTSDZ2,SSOIL,HCPCT
 587     REAL                   :: DF1K,DF1N
 588     REAL                   :: ZMD
 589     REAL, PARAMETER        :: ZBOT = -25.0
 590
 591 ! ----------------------------------------------------------------------
 592 ! SET A NOMINAL UNIVERSAL VALUE OF GLACIAL-ICE SPECIFIC HEAT CAPACITY,
 593 !   HCPCT = 2100.0*900.0 = 1.89000E+6 (SOURCE:  BOB GRUMBINE, 2005)
 594 !   TBOT PASSED IN AS ARGUMENT, VALUE FROM GLOBAL DATA SET
 595     !
 596     ! A least-squares fit for the four points provided by
 597     ! Keith Hines for the Yen (1981) values for Antarctic
 598     ! snow firn.
 599     !
 600     HCPCT = 1.E6 * (0.8194 - 0.1309*0.5*ZSOIL(1))
 601     DF1K = DF1
 602
 603 ! ----------------------------------------------------------------------
 604 ! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
 605 ! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
 606 ! ----------------------------------------------------------------------
 607 ! SET ICE PACK DEPTH.  USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
 608 ! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK).  ASSUME ICE
 609 ! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
 610 ! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
 611 ! ----------------------------------------------------------------------
 612 ! ----------------------------------------------------------------------
 613 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
 614 ! ----------------------------------------------------------------------
 615     DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
 616     AI (1) = 0.0
 617     CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
 618
 619 ! ----------------------------------------------------------------------
 620 ! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
 621 ! RECALC/ADJUST THE SOIL HEAT FLUX.  USE THE GRADIENT AND FLUX TO CALC
 622 ! RHSTS FOR THE TOP SOIL LAYER.
 623 ! ----------------------------------------------------------------------
 624     BI (1) = - CI (1) + DF1/ (0.5 * ZSOIL (1) * ZSOIL (1) * HCPCT *    &
 625          &   ZZ1)
 626     DTSDZ = ( STC (1) - STC (2) ) / ( -0.5 * ZSOIL (2) )
 627     SSOIL = DF1 * ( STC (1) - YY ) / ( 0.5 * ZSOIL (1) * ZZ1 )
 628
 629 ! ----------------------------------------------------------------------
 630 ! INITIALIZE DDZ2
 631 ! ----------------------------------------------------------------------
 632     RHSTS (1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL (1) * HCPCT )
 633
 634 ! ----------------------------------------------------------------------
 635 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
 636 ! ----------------------------------------------------------------------
 637     DDZ2 = 0.0
 638     DF1K = DF1
 639     DF1N = DF1
 640     DO K = 2,NSOIL
 641
 642        ZMD = 0.5 * (ZSOIL(K)+ZSOIL(K-1))
 643        ! For the land-ice case
 644 ! kmh 09/03/2006 use Yen (1981)'s values for Antarctic snow firn
 645 !         IF ( K .eq. 2 ) HCPCT = 0.855108E6
 646 !         IF ( K .eq. 3 ) HCPCT = 0.922906E6
 647 !         IF ( K .eq. 4 ) HCPCT = 1.009986E6
 648
 649        ! Least squares fit to the four points supplied by Keith Hines
 650        ! from Yen (1981) for Antarctic snow firn.  Not optimal, but
 651        ! probably better than just a constant.
 652        HCPCT = 1.E6 * ( 0.8194 - 0.1309*ZMD )
 653
 654 !         IF ( K .eq. 2 ) DF1N = 0.345356
 655 !         IF ( K .eq. 3 ) DF1N = 0.398777
 656 !         IF ( K .eq. 4 ) DF1N = 0.472653
 657
 658        ! Least squares fit to the three points supplied by Keith Hines
 659        ! from Yen (1981) for Antarctic snow firn.  Not optimal, but
 660        ! probably better than just a constant.
 661        DF1N = 0.32333 - ( 0.10073 * ZMD )
 662 ! ----------------------------------------------------------------------
 663 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
 664 ! ----------------------------------------------------------------------
 665        IF (K /= NSOIL) THEN
 666           DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
 667
 668 ! ----------------------------------------------------------------------
 669 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
 670 ! ----------------------------------------------------------------------
 671           DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
 672           DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
 673           CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K))*HCPCT)
 674
 675 ! ----------------------------------------------------------------------
 676 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
 677 ! ----------------------------------------------------------------------
 678        ELSE
 679
 680 ! ----------------------------------------------------------------------
 681 ! SET MATRIX COEF, CI TO ZERO.
 682 ! ----------------------------------------------------------------------
 683           DTSDZ2 = (STC (K) - TBOT)/ (.5 * (ZSOIL (K -1) + ZSOIL (K)) &
 684                &   - ZBOT)
 685           CI (K) = 0.
 686 ! ----------------------------------------------------------------------
 687 ! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
 688 ! ----------------------------------------------------------------------
 689        END IF
 690        DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
 691
 692 ! ----------------------------------------------------------------------
 693 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
 694 ! ----------------------------------------------------------------------
 695        RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
 696        AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
 697
 698 ! ----------------------------------------------------------------------
 699 ! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
 700 ! ----------------------------------------------------------------------
 701        BI (K) = - (AI (K) + CI (K))
 702        DF1K = DF1N
 703        DTSDZ = DTSDZ2
 704        DDZ = DDZ2
 705     END DO
 706 ! ----------------------------------------------------------------------
 707   END SUBROUTINE HRTICE
 708 ! ----------------------------------------------------------------------
 709
 710   SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
 711
 712 ! ----------------------------------------------------------------------
 713 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
 714 ! ----------------------------------------------------------------------
 715     IMPLICIT NONE
 716     INTEGER,                  INTENT(IN)    :: NSOIL
 717     REAL, DIMENSION(1:NSOIL), INTENT(IN)    :: STCIN
 718     REAL, DIMENSION(1:NSOIL), INTENT(OUT)   :: STCOUT
 719     REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: RHSTS
 720     REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: AI,BI,CI
 721     REAL, DIMENSION(1:NSOIL) :: RHSTSin
 722     REAL, DIMENSION(1:NSOIL) :: CIin
 723     REAL                     :: DT
 724     INTEGER                  :: K
 725
 726 ! ----------------------------------------------------------------------
 727 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
 728 ! ----------------------------------------------------------------------
 729     DO K = 1,NSOIL
 730        RHSTS (K) = RHSTS (K) * DT
 731        AI (K) = AI (K) * DT
 732        BI (K) = 1. + BI (K) * DT
 733        CI (K) = CI (K) * DT
 734     END DO
 735 ! ----------------------------------------------------------------------
 736 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
 737 ! ----------------------------------------------------------------------
 738     DO K = 1,NSOIL
 739        RHSTSin (K) = RHSTS (K)
 740     END DO
 741     DO K = 1,NSOIL
 742        CIin (K) = CI (K)
 743     END DO
 744 ! ----------------------------------------------------------------------
 745 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
 746 ! ----------------------------------------------------------------------
 747     CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
 748 ! ----------------------------------------------------------------------
 749 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
 750 ! ----------------------------------------------------------------------
 751     DO K = 1,NSOIL
 752        STCOUT (K) = STCIN (K) + CI (K)
 753     END DO
 754 ! ----------------------------------------------------------------------
 755   END SUBROUTINE HSTEP
 756 ! ----------------------------------------------------------------------
 757
 758   SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,TH2,PRCP,FDOWN,T24,SSOIL, &
 759        &             Q2,Q2SAT,ETP,RCH,RR,SNOWNG,FRZGRA,       &
 760        &             DQSDT2,FLX2,EMISSI,T1)
 761
 762 ! ----------------------------------------------------------------------
 763 ! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT.  VARIOUS
 764 ! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
 765 ! CALLING ROUTINE FOR LATER USE.
 766 ! ----------------------------------------------------------------------
 767     IMPLICIT NONE
 768     LOGICAL, INTENT(IN)  :: SNOWNG, FRZGRA
 769     REAL,    INTENT(IN)  :: CH, DQSDT2,FDOWN,PRCP,Q2,Q2SAT,SSOIL,SFCPRS, &
 770          &                  SFCTMP,TH2,EMISSI,T1,RCH,T24
 771     REAL,    INTENT(OUT) :: ETP,FLX2,RR
 772
 773     REAL                 :: A, DELTA, FNET,RAD,ELCP1,LVS,EPSCA
 774
 775     REAL, PARAMETER      :: ELCP = 2.4888E+3, LSUBC = 2.501000E+6
 776     REAL, PARAMETER      :: LSUBS = 2.83E+6
 777
 778 ! ----------------------------------------------------------------------
 779 ! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
 780 ! ----------------------------------------------------------------------
 781     IF ( T1 > 273.15 ) THEN
 782        ELCP1 = ELCP
 783        LVS   = LSUBC
 784     ELSE
 785        ELCP1 = ELCP*LSUBS/LSUBC
 786        LVS   = LSUBS
 787     ENDIF
 788     DELTA = ELCP1 * DQSDT2
 789     A = ELCP1 * (Q2SAT - Q2)
 790     RR = EMISSI*T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
 791
 792 ! ----------------------------------------------------------------------
 793 ! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
 794 ! EFFECTS CAUSED BY FALLING PRECIPITATION.
 795 ! ----------------------------------------------------------------------
 796     IF (.NOT. SNOWNG) THEN
 797        IF (PRCP >  0.0) RR = RR + CPH2O * PRCP / RCH
 798     ELSE
 799        RR = RR + CPICE * PRCP / RCH
 800     END IF
 801
 802 ! ----------------------------------------------------------------------
 803 ! INCLUDE THE LATENT HEAT EFFECTS OF FREEZING RAIN CONVERTING TO ICE ON
 804 ! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
 805 ! ----------------------------------------------------------------------
 806     IF (FRZGRA) THEN
 807        FLX2 = - LSUBF * PRCP
 808     ELSE
 809        FLX2 = 0.0
 810     ENDIF
 811     FNET = FDOWN -  ( EMISSI * SIGMA * T24 ) - SSOIL - FLX2
 812
 813 ! ----------------------------------------------------------------------
 814 ! FINISH PENMAN EQUATION CALCULATIONS.
 815 ! ----------------------------------------------------------------------
 816     RAD = FNET / RCH + TH2 - SFCTMP
 817     EPSCA = (A * RR + RAD * DELTA) / (DELTA + RR)
 818     ETP = EPSCA * RCH / LVS
 819
 820 ! ----------------------------------------------------------------------
 821   END SUBROUTINE PENMAN
 822 ! ----------------------------------------------------------------------
 823
 824   SUBROUTINE SHFLX (STC,NSOIL,DT,YY,ZZ1,ZSOIL,TBOT,DF1)
 825 ! ----------------------------------------------------------------------
 826 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
 827 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
 828 ! ON THE TEMPERATURE.
 829 ! ----------------------------------------------------------------------
 830     IMPLICIT NONE
 831
 832     INTEGER,                  INTENT(IN)    :: NSOIL
 833     REAL,                     INTENT(IN)    :: DF1,DT,TBOT,YY, ZZ1
 834     REAL, DIMENSION(1:NSOIL), INTENT(IN)    :: ZSOIL
 835     REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
 836
 837     REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS
 838     INTEGER                  :: I
 839     REAL, PARAMETER          :: T0 = 273.15
 840
 841 ! ----------------------------------------------------------------------
 842 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
 843 ! ----------------------------------------------------------------------
 844
 845     CALL HRTICE (RHSTS,STC,TBOT, NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
 846
 847     CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
 848
 849     DO I = 1,NSOIL
 850        STC (I) = STCF (I)
 851     END DO
 852 ! ----------------------------------------------------------------------
 853   END SUBROUTINE SHFLX
 854 ! ----------------------------------------------------------------------
 855
 856   SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,NSOIL,DT,DF1,           &
 857        &             Q2,T1,SFCTMP,T24,TH2,FDOWN,SSOIL,STC,             &
 858        &             SFCPRS,RCH,RR,SNEQV,SNDENS,SNOWH,ZSOIL,TBOT,      &
 859        &             SNOMLT,DEW,FLX1,FLX2,FLX3,ESNOW,EMISSI,RIBB)
 860
 861 ! ----------------------------------------------------------------------
 862 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
 863 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
 864 ! PRESENT.
 865 ! ----------------------------------------------------------------------
 866     IMPLICIT NONE
 867
 868     INTEGER, INTENT(IN)   :: NSOIL
 869     LOGICAL, INTENT(IN)   :: SNOWNG
 870     REAL,    INTENT(IN)   :: DF1,DT,FDOWN,PRCP,Q2,RCH,RR,SFCPRS,SFCTMP, &
 871          &                   T24,TBOT,TH2,EMISSI
 872     REAL, INTENT(INOUT)   :: SNEQV,FLX2,PRCPF,SNOWH,SNDENS,T1,RIBB,ETP
 873     REAL, INTENT(OUT)     :: DEW,ESNOW,FLX1,FLX3,SSOIL,SNOMLT
 874     REAL, DIMENSION(1:NSOIL),INTENT(IN)     :: ZSOIL
 875     REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
 876     REAL, DIMENSION(1:NSOIL) :: ET1
 877     INTEGER               :: K
 878     REAL                  :: DENOM,DSOIL,DTOT,ESDFLX,ETA,        &
 879          &                   ESNOW1,ESNOW2,ETA1,ETP1,ETP2,       &
 880          &                   ETP3,ETANRG,EX,                     &
 881          &                   FRCSNO,FRCSOI,PRCP1,QSAT,RSNOW,SEH, &
 882          &                   SNCOND,T12,T12A,T12B,T14,YY,ZZ1
 883
 884     REAL, PARAMETER       :: ESDMIN = 1.E-6, LSUBC = 2.501000E+6,     &
 885          &                   LSUBS = 2.83E+6, TFREEZ = 273.15,        &
 886          &                   SNOEXP = 2.0
 887
 888 ! ----------------------------------------------------------------------
 889 ! FOR GLACIAL-ICE, SNOWCOVER FRACTION = 1.0, AND SUBLIMATION IS AT THE
 890 ! POTENTIAL RATE.
 891 ! ----------------------------------------------------------------------
 892 ! INITIALIZE EVAP TERMS.
 893 ! ----------------------------------------------------------------------
 894 ! conversions:
 895 ! ESNOW [KG M-2 S-1]
 896 ! ESDFLX [KG M-2 S-1] .le. ESNOW
 897 ! ESNOW1 [M S-1]
 898 ! ESNOW2 [M]
 899 ! ETP [KG M-2 S-1]
 900 ! ETP1 [M S-1]
 901 ! ETP2 [M]
 902 ! ----------------------------------------------------------------------
 903     SNOMLT = 0.0
 904     DEW = 0.
 905     ESNOW = 0.
 906     ESNOW1 = 0.
 907     ESNOW2 = 0.
 908
 909 ! ----------------------------------------------------------------------
 910 ! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO ETP1 IN M S-1
 911 ! ----------------------------------------------------------------------
 912     PRCP1 = PRCPF *0.001
 913 ! ----------------------------------------------------------------------
 914 ! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
 915 ! ----------------------------------------------------------------------
 916     IF (ETP <= 0.0) THEN
 917        IF ( ( RIBB >= 0.1 ) .AND. ( FDOWN > 150.0 ) ) THEN
 918           ETP=(MIN(ETP*(1.0-RIBB),0.)/0.980 + ETP*(0.980-1.0))/0.980
 919        ENDIF
 920        ETP1 = ETP * 0.001
 921        DEW = -ETP1
 922        ESNOW2 = ETP1*DT
 923        ETANRG = ETP*LSUBS
 924     ELSE
 925        ETP1 = ETP * 0.001
 926        ESNOW = ETP
 927        ESNOW1 = ESNOW*0.001
 928        ESNOW2 = ESNOW1*DT
 929        ETANRG = ESNOW*LSUBS
 930     END IF
 931
 932 ! ----------------------------------------------------------------------
 933 ! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
 934 ! ACCUMULATING PRECIP.  NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
 935 ! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1).  ASSUMES TEMPERATURE OF THE
 936 ! SNOWFALL STRIKING THE GROUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
 937 ! ----------------------------------------------------------------------
 938     FLX1 = 0.0
 939     IF (SNOWNG) THEN
 940        FLX1 = CPICE * PRCP * (T1- SFCTMP)
 941     ELSE
 942        IF (PRCP >  0.0) FLX1 = CPH2O * PRCP * (T1- SFCTMP)
 943     END IF
 944 ! ----------------------------------------------------------------------
 945 ! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
 946 ! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
 947 ! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
 948 ! FLUXES.  FLX1 FROM ABOVE, FLX2 BROUGHT IN VIA COMMOM BLOCK RITE.
 949 ! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
 950 ! PENMAN.
 951 ! ----------------------------------------------------------------------
 952     DSOIL = - (0.5 * ZSOIL (1))
 953     DTOT = SNOWH + DSOIL
 954     DENOM = 1.0+ DF1 / (DTOT * RR * RCH)
 955     T12A = ( (FDOWN - FLX1- FLX2- EMISSI * SIGMA * T24)/ RCH                    &
 956          + TH2- SFCTMP - ETANRG / RCH ) / RR
 957     T12B = DF1 * STC (1) / (DTOT * RR * RCH)
 958
 959     T12 = (SFCTMP + T12A + T12B) / DENOM
 960     IF (T12 <=  TFREEZ) THEN
 961 ! ----------------------------------------------------------------------
 962 ! SUB-FREEZING BLOCK
 963 ! ----------------------------------------------------------------------
 964 ! ----------------------------------------------------------------------
 965 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
 966 ! MELT WILL OCCUR.  SET THE SKIN TEMP TO THIS EFFECTIVE TEMP.  REDUCE
 967 ! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
 968 ! DEPENDING ON SIGN OF ETP.
 969 ! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
 970 ! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
 971 ! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
 972 ! TO ZERO.
 973 ! ----------------------------------------------------------------------
 974        T1 = T12
 975        SSOIL = DF1 * (T1- STC (1)) / DTOT
 976        SNEQV = MAX(0.0, SNEQV-ESNOW2)
 977        FLX3 = 0.0
 978        EX = 0.0
 979        SNOMLT = 0.0
 980     ELSE
 981 ! ----------------------------------------------------------------------
 982 ! ABOVE FREEZING BLOCK
 983 ! ----------------------------------------------------------------------
 984 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
 985 ! WILL OCCUR.  CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT.  REVISE THE
 986 ! EFFECTIVE SNOW DEPTH.  REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
 987 ! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
 988 ! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
 989 ! EX FOR USE IN SMFLX.  ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
 990 ! CALCULATE QSAT VALID AT FREEZING POINT.  NOTE THAT ESAT (SATURATION
 991 ! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
 992 ! POINT.  NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
 993 ! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
 994 ! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
 995 ! ----------------------------------------------------------------------
 996        T1 = TFREEZ
 997        IF ( DTOT .GT. 2.0*DSOIL ) THEN
 998           DTOT = 2.0*DSOIL
 999        ENDIF
1000        SSOIL = DF1 * (T1- STC (1)) / DTOT
1001        IF (SNEQV-ESNOW2 <= ESDMIN) THEN
1002           SNEQV = 0.0
1003           EX = 0.0
1004           SNOMLT = 0.0
1005           FLX3 = 0.0
1006 ! ----------------------------------------------------------------------
1007 ! SUBLIMATION LESS THAN DEPTH OF SNOWPACK
1008 ! SNOWPACK (SNEQV) REDUCED BY ESNOW2 (DEPTH OF SUBLIMATED SNOW)
1009 ! ----------------------------------------------------------------------
1010        ELSE
1011           SNEQV = SNEQV-ESNOW2
1012           ETP3 = ETP * LSUBC
1013           SEH = RCH * (T1- TH2)
1014           T14 = ( T1 * T1 ) * ( T1 * T1 )
1015           FLX3 = FDOWN - FLX1- FLX2- EMISSI*SIGMA * T14- SSOIL - SEH - ETANRG
1016           IF (FLX3 <= 0.0) FLX3 = 0.0
1017           EX = FLX3*0.001/ LSUBF
1018           SNOMLT = EX * DT
1019 ! ----------------------------------------------------------------------
1020 ! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
1021 ! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
1022 ! ----------------------------------------------------------------------
1023           IF (SNEQV- SNOMLT >=  ESDMIN) THEN
1024              SNEQV = SNEQV- SNOMLT
1025           ELSE
1026 ! ----------------------------------------------------------------------
1027 ! SNOWMELT EXCEEDS SNOW DEPTH
1028 ! ----------------------------------------------------------------------
1029              EX = SNEQV / DT
1030              FLX3 = EX *1000.0* LSUBF
1031              SNOMLT = SNEQV
1032
1033              SNEQV = 0.0
1034           ENDIF
1035        ENDIF
1036
1037 ! ----------------------------------------------------------------------
1038 ! FOR GLACIAL ICE, THE SNOWMELT WILL BE ADDED TO SUBSURFACE
1039 ! RUNOFF/BASEFLOW LATER NEAR THE END OF SFLX (AFTER RETURN FROM CALL TO
1040 ! SUBROUTINE SNOPAC)
1041 ! ----------------------------------------------------------------------
1042
1043     ENDIF
1044
1045 ! ----------------------------------------------------------------------
1046 ! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
1047 ! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
1048 ! MATCHES THAT ALREADY COMPUTED FOR BELOW THE SNOWPACK, THUS THE SFC
1049 ! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
1050 ! SNOW TOP SURFACE.
1051 ! ----------------------------------------------------------------------
1052     ZZ1 = 1.0
1053     YY = STC (1) -0.5* SSOIL * ZSOIL (1)* ZZ1/ DF1
1054
1055 ! ----------------------------------------------------------------------
1056 ! SHFLX WILL CALC/UPDATE THE SOIL TEMPS.
1057 ! ----------------------------------------------------------------------
1058     CALL SHFLX (STC,NSOIL,DT,YY,ZZ1,ZSOIL,TBOT,DF1)
1059
1060 ! ----------------------------------------------------------------------
1061 ! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION.  YY IS
1062 ! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
1063 ! ----------------------------------------------------------------------
1064     IF (SNEQV .GE. 0.10) THEN
1065        CALL SNOWPACK (SNEQV,DT,SNOWH,SNDENS,T1,YY)
1066     ELSE
1067        SNEQV = 0.10
1068        SNOWH = 0.50
1069 !KWM???? SNDENS =
1070 !KWM???? SNCOND =
1071     ENDIF
1072 ! ----------------------------------------------------------------------
1073   END SUBROUTINE SNOPAC
1074 ! ----------------------------------------------------------------------
1075
1076   SUBROUTINE SNOWPACK (SNEQV,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
1077
1078 ! ----------------------------------------------------------------------
1079 ! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
1080 ! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
1081 ! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
1082 ! KOREN, 03/25/95.
1083 ! ----------------------------------------------------------------------
1084 ! SNEQV     WATER EQUIVALENT OF SNOW (M)
1085 ! DTSEC   TIME STEP (SEC)
1086 ! SNOWH   SNOW DEPTH (M)
1087 ! SNDENS  SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
1088 ! TSNOW   SNOW SURFACE TEMPERATURE (K)
1089 ! TSOIL   SOIL SURFACE TEMPERATURE (K)
1090
1091 ! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
1092 ! ----------------------------------------------------------------------
1093       IMPLICIT NONE
1094
1095       INTEGER                :: IPOL, J
1096       REAL, INTENT(IN)       :: SNEQV, DTSEC,TSNOW,TSOIL
1097       REAL, INTENT(INOUT)    :: SNOWH, SNDENS
1098       REAL                   :: BFAC,DSX,DTHR,DW,SNOWHC,PEXP,           &
1099                                 TAVGC,TSNOWC,TSOILC,ESDC,ESDCX
1100       REAL, PARAMETER        :: C1 = 0.01, C2 = 21.0, G = 9.81,         &
1101                                 KN = 4000.0
1102 ! ----------------------------------------------------------------------
1103 ! CONVERSION INTO SIMULATION UNITS
1104 ! ----------------------------------------------------------------------
1105       SNOWHC = SNOWH *100.
1106       ESDC   = SNEQV *100.
1107       DTHR   = DTSEC /3600.
1108       TSNOWC = TSNOW -273.15
1109       TSOILC = TSOIL -273.15
1110
1111 ! ----------------------------------------------------------------------
1112 ! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
1113 ! ----------------------------------------------------------------------
1114 ! ----------------------------------------------------------------------
1115 ! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
1116 !  SNDENS=DS0*(EXP(BFAC*SNEQV)-1.)/(BFAC*SNEQV)
1117 !  BFAC=DTHR*C1*EXP(0.08*TAVGC-C2*DS0)
1118 ! NOTE: BFAC*SNEQV IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
1119 ! NUMERICALLY BELOW:
1120 !   C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
1121 !   C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
1122 ! ----------------------------------------------------------------------
1123       TAVGC = 0.5* (TSNOWC + TSOILC)
1124       IF (ESDC >  1.E-2) THEN
1125          ESDCX = ESDC
1126       ELSE
1127          ESDCX = 1.E-2
1128       END IF
1129
1130 !      DSX = SNDENS*((DEXP(BFAC*ESDC)-1.)/(BFAC*ESDC))
1131 ! ----------------------------------------------------------------------
1132 ! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
1133 ! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
1134 ! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
1135 ! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
1136 ! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
1137 ! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
1138 ! POLYNOMIAL EXPANSION.
1139
1140 ! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
1141 ! IS GOVERNED BY ITERATION LIMIT "IPOL".
1142 !      IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
1143 !            PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
1144 !       IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
1145 !       IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
1146 !       IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
1147 ! ----------------------------------------------------------------------
1148       BFAC = DTHR * C1* EXP (0.08* TAVGC - C2* SNDENS)
1149       IPOL = 4
1150       PEXP = 0.
1151 !        PEXP = (1. + PEXP)*BFAC*ESDC/REAL(J+1)
1152       DO J = IPOL,1, -1
1153          PEXP = (1. + PEXP)* BFAC * ESDCX / REAL (J +1)
1154       END DO
1155
1156       PEXP = PEXP + 1.
1157 ! ----------------------------------------------------------------------
1158 ! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
1159 ! ----------------------------------------------------------------------
1160 !     END OF KOREAN FORMULATION
1161
1162 !     BASE FORMULATION (COGLEY ET AL., 1990)
1163 !     CONVERT DENSITY FROM G/CM3 TO KG/M3
1164 !       DSM=SNDENS*1000.0
1165
1166 !       DSX=DSM+DTSEC*0.5*DSM*G*SNEQV/
1167 !    &      (1E7*EXP(-0.02*DSM+KN/(TAVGC+273.16)-14.643))
1168
1169 !  &   CONVERT DENSITY FROM KG/M3 TO G/CM3
1170 !       DSX=DSX/1000.0
1171
1172 !     END OF COGLEY ET AL. FORMULATION
1173
1174 ! ----------------------------------------------------------------------
1175 ! SET UPPER/LOWER LIMIT ON SNOW DENSITY
1176 ! ----------------------------------------------------------------------
1177       DSX = SNDENS * (PEXP)
1178       IF (DSX > 0.40) DSX = 0.40
1179       IF (DSX < 0.05) DSX = 0.05
1180 ! ----------------------------------------------------------------------
1181 ! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
1182 ! SNOWMELT.  ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
1183 ! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
1184 ! ----------------------------------------------------------------------
1185       SNDENS = DSX
1186       IF (TSNOWC >=  0.) THEN
1187          DW = 0.13* DTHR /24.
1188          SNDENS = SNDENS * (1. - DW) + DW
1189          IF (SNDENS >=  0.40) SNDENS = 0.40
1190 ! ----------------------------------------------------------------------
1191 ! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
1192 ! CHANGE SNOW DEPTH UNITS TO METERS
1193 ! ----------------------------------------------------------------------
1194       END IF
1195       SNOWHC = ESDC / SNDENS
1196       SNOWH = SNOWHC * 0.01
1197
1198 ! ----------------------------------------------------------------------
1199   END SUBROUTINE SNOWPACK
1200 ! ----------------------------------------------------------------------
1201
1202   SUBROUTINE SNOWZ0 (Z0, Z0BRD, SNOWH)
1203 ! ----------------------------------------------------------------------
1204 ! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
1205 ! Z0      ROUGHNESS LENGTH (m)
1206 ! Z0S     SNOW ROUGHNESS LENGTH:=0.001 (m)
1207 ! ----------------------------------------------------------------------
1208     IMPLICIT NONE
1209     REAL, INTENT(IN)        :: Z0BRD
1210     REAL, INTENT(OUT)       :: Z0
1211     REAL, PARAMETER         :: Z0S=0.001
1212     REAL, INTENT(IN)        :: SNOWH
1213     REAL                    :: BURIAL
1214     REAL                    :: Z0EFF
1215
1216     BURIAL = 7.0*Z0BRD - SNOWH
1217     IF(BURIAL.LE.0.0007) THEN
1218        Z0EFF = Z0S
1219     ELSE
1220        Z0EFF = BURIAL/7.0
1221     ENDIF
1222
1223     Z0 = Z0EFF
1224
1225 ! ----------------------------------------------------------------------
1226   END SUBROUTINE SNOWZ0
1227 ! ----------------------------------------------------------------------
1228
1229   SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
1230
1231 ! ----------------------------------------------------------------------
1232 ! CALCULATE SNOW DEPTH AND DENSITY TO ACCOUNT FOR THE NEW SNOWFALL.
1233 ! UPDATED VALUES OF SNOW DEPTH AND DENSITY ARE RETURNED.
1234
1235 ! TEMP    AIR TEMPERATURE (K)
1236 ! NEWSN   NEW SNOWFALL (M)
1237 ! SNOWH   SNOW DEPTH (M)
1238 ! SNDENS  SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
1239 ! ----------------------------------------------------------------------
1240     IMPLICIT NONE
1241     REAL, INTENT(IN)        :: NEWSN, TEMP
1242     REAL, INTENT(INOUT)     :: SNDENS, SNOWH
1243     REAL                    :: DSNEW, HNEWC, SNOWHC,NEWSNC,TEMPC
1244
1245 ! ----------------------------------------------------------------------
1246 ! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
1247 ! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
1248 ! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
1249 ! VEMADOLEN, SWEDEN, 1980, 172-177PP.
1250 !-----------------------------------------------------------------------
1251     TEMPC = TEMP - 273.15
1252     IF ( TEMPC <=  -15. ) THEN
1253        DSNEW = 0.05
1254     ELSE
1255        DSNEW = 0.05 + 0.0017 * ( TEMPC + 15. ) ** 1.5
1256     ENDIF
1257
1258 ! ----------------------------------------------------------------------
1259 ! CONVERSION INTO SIMULATION UNITS
1260 ! ----------------------------------------------------------------------
1261     SNOWHC = SNOWH * 100.
1262     NEWSNC = NEWSN * 100.
1263
1264 ! ----------------------------------------------------------------------
1265 ! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
1266 ! ----------------------------------------------------------------------
1267     HNEWC = NEWSNC / DSNEW
1268     IF ( SNOWHC + HNEWC < 1.0E-3 ) THEN
1269        SNDENS = MAX ( DSNEW , SNDENS )
1270     ELSE
1271        SNDENS = ( SNOWHC * SNDENS + HNEWC * DSNEW ) / ( SNOWHC + HNEWC )
1272     ENDIF
1273     SNOWHC = SNOWHC + HNEWC
1274     SNOWH = SNOWHC * 0.01
1275
1276 ! ----------------------------------------------------------------------
1277   END SUBROUTINE SNOW_NEW
1278 ! ----------------------------------------------------------------------
1279
1280 END MODULE module_sf_noahlsm_glacial_only