phys/module_sf_noah_seaice.F

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