MODULE NOAHMP_GLACIER_GLOBALS 2
implicit none
! ==================================================================================================
!------------------------------------------------------------------------------------------!
! Physical Constants: !
!------------------------------------------------------------------------------------------!
REAL, PARAMETER :: GRAV = 9.80616 !acceleration due to gravity (m/s2)
REAL, PARAMETER :: SB = 5.67E-08 !Stefan-Boltzmann constant (w/m2/k4)
REAL, PARAMETER :: VKC = 0.40 !von Karman constant
REAL, PARAMETER :: TFRZ = 273.16 !freezing/melting point (k)
REAL, PARAMETER :: HSUB = 2.8440E06 !latent heat of sublimation (j/kg)
REAL, PARAMETER :: HVAP = 2.5104E06 !latent heat of vaporization (j/kg)
REAL, PARAMETER :: HFUS = 0.3336E06 !latent heat of fusion (j/kg)
REAL, PARAMETER :: CWAT = 4.188E06 !specific heat capacity of water (j/m3/k)
REAL, PARAMETER :: CICE = 2.094E06 !specific heat capacity of ice (j/m3/k)
REAL, PARAMETER :: CPAIR = 1004.64 !heat capacity dry air at const pres (j/kg/k)
REAL, PARAMETER :: TKWAT = 0.6 !thermal conductivity of water (w/m/k)
REAL, PARAMETER :: TKICE = 2.2 !thermal conductivity of ice (w/m/k)
REAL, PARAMETER :: TKAIR = 0.023 !thermal conductivity of air (w/m/k)
REAL, PARAMETER :: RAIR = 287.04 !gas constant for dry air (j/kg/k)
REAL, PARAMETER :: RW = 461.269 !gas constant for water vapor (j/kg/k)
REAL, PARAMETER :: DENH2O = 1000. !density of water (kg/m3)
REAL, PARAMETER :: DENICE = 917. !density of ice (kg/m3)
! =====================================options for different schemes================================
! options for dynamic vegetation:
! 1 -> off (use table LAI; use FVEG = SHDFAC from input)
! 2 -> on (together with OPT_CRS = 1)
! 3 -> off (use table LAI; calculate FVEG)
! 4 -> off (use table LAI; use maximum vegetation fraction)
INTEGER :: DVEG != 2 !
! options for canopy stomatal resistance
! 1-> Ball-Berry; 2->Jarvis
INTEGER :: OPT_CRS != 1 !(must 1 when DVEG = 2)
! options for soil moisture factor for stomatal resistance
! 1-> Noah (soil moisture)
! 2-> CLM (matric potential)
! 3-> SSiB (matric potential)
INTEGER :: OPT_BTR != 1 !(suggested 1)
! options for runoff and groundwater
! 1 -> TOPMODEL with groundwater (Niu et al. 2007 JGR) ;
! 2 -> TOPMODEL with an equilibrium water table (Niu et al. 2005 JGR) ;
! 3 -> original surface and subsurface runoff (free drainage)
! 4 -> BATS surface and subsurface runoff (free drainage)
INTEGER :: OPT_RUN != 1 !(suggested 1)
! options for surface layer drag coeff (CH & CM)
! 1->M-O ; 2->original Noah (Chen97); 3->MYJ consistent; 4->YSU consistent.
INTEGER :: OPT_SFC != 1 !(1 or 2 or 3 or 4)
! options for supercooled liquid water (or ice fraction)
! 1-> no iteration (Niu and Yang, 2006 JHM); 2: Koren's iteration
INTEGER :: OPT_FRZ != 1 !(1 or 2)
! options for frozen soil permeability
! 1 -> linear effects, more permeable (Niu and Yang, 2006, JHM)
! 2 -> nonlinear effects, less permeable (old)
INTEGER :: OPT_INF != 1 !(suggested 1)
! options for radiation transfer
! 1 -> modified two-stream (gap = F(solar angle, 3D structure ...)<1-FVEG)
! 2 -> two-stream applied to grid-cell (gap = 0)
! 3 -> two-stream applied to vegetated fraction (gap=1-FVEG)
INTEGER :: OPT_RAD != 1 !(suggested 1)
! options for ground snow surface albedo
! 1-> BATS; 2 -> CLASS
INTEGER :: OPT_ALB != 2 !(suggested 2)
! options for partitioning precipitation into rainfall & snowfall
! 1 -> Jordan (1991); 2 -> BATS: when SFCTMP<TFRZ+2.2 ; 3-> SFCTMP<TFRZ
INTEGER :: OPT_SNF != 1 !(suggested 1)
! options for lower boundary condition of soil temperature
! 1 -> zero heat flux from bottom (ZBOT and TBOT not used)
! 2 -> TBOT at ZBOT (8m) read from a file (original Noah)
INTEGER :: OPT_TBOT != 2 !(suggested 2)
! options for snow/soil temperature time scheme (only layer 1)
! 1 -> semi-implicit; 2 -> full implicit (original Noah)
INTEGER :: OPT_STC != 1 !(suggested 1)
! adjustable parameters for snow processes
REAL, PARAMETER :: Z0SNO = 0.002 !snow surface roughness length (m) (0.002)
REAL, PARAMETER :: SSI = 0.03 !liquid water holding capacity for snowpack (m3/m3) (0.03)
REAL, PARAMETER :: SWEMX = 1.00 !new snow mass to fully cover old snow (mm)
!equivalent to 10mm depth (density = 100 kg/m3)
!------------------------------------------------------------------------------------------!
END MODULE NOAHMP_GLACIER_GLOBALS
!------------------------------------------------------------------------------------------!
MODULE NOAHMP_GLACIER_ROUTINES 1,1
USE NOAHMP_GLACIER_GLOBALS
IMPLICIT NONE
public :: NOAHMP_OPTIONS_GLACIER
public :: NOAHMP_GLACIER
private :: ATM_GLACIER
private :: ENERGY_GLACIER
private :: THERMOPROP_GLACIER
private :: CSNOW_GLACIER
private :: RADIATION_GLACIER
private :: SNOW_AGE_GLACIER
private :: SNOWALB_BATS_GLACIER
private :: SNOWALB_CLASS_GLACIER
private :: GLACIER_FLUX
private :: SFCDIF1_GLACIER
private :: TSNOSOI_GLACIER
private :: HRT_GLACIER
private :: HSTEP_GLACIER
private :: ROSR12_GLACIER
private :: PHASECHANGE_GLACIER
private :: WATER_GLACIER
private :: SNOWWATER_GLACIER
private :: SNOWFALL_GLACIER
private :: COMBINE_GLACIER
private :: DIVIDE_GLACIER
private :: COMBO_GLACIER
private :: COMPACT_GLACIER
private :: SNOWH2O_GLACIER
private :: ERROR_GLACIER
contains
!
! ==================================================================================================
SUBROUTINE NOAHMP_GLACIER (& 1,4
ILOC ,JLOC ,COSZ ,NSNOW ,NSOIL ,DT , & ! IN : Time/Space/Model-related
SFCTMP ,SFCPRS ,UU ,VV ,Q2 ,SOLDN , & ! IN : Forcing
PRCP ,LWDN ,TBOT ,ZLVL ,FICEOLD ,ZSOIL , & ! IN : Forcing
QSNOW ,SNEQVO ,ALBOLD ,CM ,CH ,ISNOW , & ! IN/OUT :
SNEQV ,SMC ,ZSNSO ,SNOWH ,SNICE ,SNLIQ , & ! IN/OUT :
TG ,STC ,SH2O ,TAUSS ,QSFC , & ! IN/OUT :
FSA ,FSR ,FIRA ,FSH ,FGEV ,SSOIL , & ! OUT :
TRAD ,EDIR ,RUNSRF ,RUNSUB ,SAG ,ALBEDO , & ! OUT :
QSNBOT ,PONDING ,PONDING1,PONDING2,T2M ,Q2E , & ! OUT :
EMISSI, FPICE, CH2B) ! OUT :
! --------------------------------------------------------------------------------------------------
! Initial code: Guo-Yue Niu, Oct. 2007
! Modified to glacier: Michael Barlage, June 2012
! --------------------------------------------------------------------------------------------------
implicit none
! --------------------------------------------------------------------------------------------------
! input
INTEGER , INTENT(IN) :: ILOC !grid index
INTEGER , INTENT(IN) :: JLOC !grid index
REAL , INTENT(IN) :: COSZ !cosine solar zenith angle [0-1]
INTEGER , INTENT(IN) :: NSNOW !maximum no. of snow layers
INTEGER , INTENT(IN) :: NSOIL !no. of soil layers
REAL , INTENT(IN) :: DT !time step [sec]
REAL , INTENT(IN) :: SFCTMP !surface air temperature [K]
REAL , INTENT(IN) :: SFCPRS !pressure (pa)
REAL , INTENT(IN) :: UU !wind speed in eastward dir (m/s)
REAL , INTENT(IN) :: VV !wind speed in northward dir (m/s)
REAL , INTENT(IN) :: Q2 !mixing ratio (kg/kg) lowest model layer
REAL , INTENT(IN) :: SOLDN !downward shortwave radiation (w/m2)
REAL , INTENT(IN) :: PRCP !precipitation rate (kg m-2 s-1)
REAL , INTENT(IN) :: LWDN !downward longwave radiation (w/m2)
REAL , INTENT(IN) :: TBOT !bottom condition for soil temp. [K]
REAL , INTENT(IN) :: ZLVL !reference height (m)
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: FICEOLD!ice fraction at last timestep
REAL, DIMENSION( 1:NSOIL), INTENT(IN) :: ZSOIL !layer-bottom depth from soil surf (m)
! input/output : need arbitary intial values
REAL , INTENT(INOUT) :: QSNOW !snowfall [mm/s]
REAL , INTENT(INOUT) :: SNEQVO !snow mass at last time step (mm)
REAL , INTENT(INOUT) :: ALBOLD !snow albedo at last time step (CLASS type)
REAL , INTENT(INOUT) :: CM !momentum drag coefficient
REAL , INTENT(INOUT) :: CH !sensible heat exchange coefficient
! prognostic variables
INTEGER , INTENT(INOUT) :: ISNOW !actual no. of snow layers [-]
REAL , INTENT(INOUT) :: SNEQV !snow water eqv. [mm]
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SMC !soil moisture (ice + liq.) [m3/m3]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: ZSNSO !layer-bottom depth from snow surf [m]
REAL , INTENT(INOUT) :: SNOWH !snow height [m]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ !snow layer liquid water [mm]
REAL , INTENT(INOUT) :: TG !ground temperature (k)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow/soil temperature [k]
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O !liquid soil moisture [m3/m3]
REAL , INTENT(INOUT) :: TAUSS !non-dimensional snow age
REAL , INTENT(INOUT) :: QSFC !mixing ratio at lowest model layer
! output
REAL , INTENT(OUT) :: FSA !total absorbed solar radiation (w/m2)
REAL , INTENT(OUT) :: FSR !total reflected solar radiation (w/m2)
REAL , INTENT(OUT) :: FIRA !total net LW rad (w/m2) [+ to atm]
REAL , INTENT(OUT) :: FSH !total sensible heat (w/m2) [+ to atm]
REAL , INTENT(OUT) :: FGEV !ground evap heat (w/m2) [+ to atm]
REAL , INTENT(OUT) :: SSOIL !ground heat flux (w/m2) [+ to soil]
REAL , INTENT(OUT) :: TRAD !surface radiative temperature (k)
REAL , INTENT(OUT) :: EDIR !soil surface evaporation rate (mm/s]
REAL , INTENT(OUT) :: RUNSRF !surface runoff [mm/s]
REAL , INTENT(OUT) :: RUNSUB !baseflow (saturation excess) [mm/s]
REAL , INTENT(OUT) :: SAG !solar rad absorbed by ground (w/m2)
REAL , INTENT(OUT) :: ALBEDO !surface albedo [-]
REAL , INTENT(OUT) :: QSNBOT !snowmelt [mm/s]
REAL , INTENT(OUT) :: PONDING!surface ponding [mm]
REAL , INTENT(OUT) :: PONDING1!surface ponding [mm]
REAL , INTENT(OUT) :: PONDING2!surface ponding [mm]
REAL , INTENT(OUT) :: T2M !2-m air temperature over bare ground part [k]
REAL , INTENT(OUT) :: Q2E
REAL , INTENT(OUT) :: EMISSI
REAL , INTENT(OUT) :: FPICE
REAL , INTENT(OUT) :: CH2B
! local
INTEGER :: IZ !do-loop index
INTEGER, DIMENSION(-NSNOW+1:NSOIL) :: IMELT !phase change index [1-melt; 2-freeze]
REAL :: RHOAIR !density air (kg/m3)
REAL, DIMENSION(-NSNOW+1:NSOIL) :: DZSNSO !snow/soil layer thickness [m]
REAL :: THAIR !potential temperature (k)
REAL :: QAIR !specific humidity (kg/kg) (q2/(1+q2))
REAL :: EAIR !vapor pressure air (pa)
REAL, DIMENSION( 1: 2) :: SOLAD !incoming direct solar rad (w/m2)
REAL, DIMENSION( 1: 2) :: SOLAI !incoming diffuse solar rad (w/m2)
REAL, DIMENSION( 1:NSOIL) :: SICE !soil ice content (m3/m3)
REAL, DIMENSION(-NSNOW+1: 0) :: SNICEV !partial volume ice of snow [m3/m3]
REAL, DIMENSION(-NSNOW+1: 0) :: SNLIQV !partial volume liq of snow [m3/m3]
REAL, DIMENSION(-NSNOW+1: 0) :: EPORE !effective porosity [m3/m3]
REAL :: QDEW !ground surface dew rate [mm/s]
REAL :: QVAP !ground surface evap. rate [mm/s]
REAL :: LATHEA !latent heat [j/kg]
REAL :: QMELT !internal pack melt
REAL :: SWDOWN !downward solar [w/m2]
REAL :: BEG_WB !beginning water for error check
REAL :: ZBOT = -8.0
! --------------------------------------------------------------------------------------------------
! re-process atmospheric forcing
CALL ATM_GLACIER (SFCPRS ,SFCTMP ,Q2 ,SOLDN ,COSZ ,THAIR , &
QAIR ,EAIR ,RHOAIR ,SOLAD ,SOLAI ,SWDOWN )
BEG_WB = SNEQV
! snow/soil layer thickness (m); interface depth: ZSNSO < 0; layer thickness DZSNSO > 0
DO IZ = ISNOW+1, NSOIL
IF(IZ == ISNOW+1) THEN
DZSNSO(IZ) = - ZSNSO(IZ)
ELSE
DZSNSO(IZ) = ZSNSO(IZ-1) - ZSNSO(IZ)
END IF
END DO
! compute energy budget (momentum & energy fluxes and phase changes)
CALL ENERGY_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT ,QSNOW ,RHOAIR , & !in
EAIR ,SFCPRS ,QAIR ,SFCTMP ,LWDN ,UU , & !in
VV ,SOLAD ,SOLAI ,COSZ ,ZLVL , & !in
TBOT ,ZBOT ,ZSNSO ,DZSNSO , & !in
TG ,STC ,SNOWH ,SNEQV ,SNEQVO ,SH2O , & !inout
SMC ,SNICE ,SNLIQ ,ALBOLD ,CM ,CH , & !inout
TAUSS ,QSFC , & !inout
IMELT ,SNICEV ,SNLIQV ,EPORE ,QMELT ,PONDING, & !out
SAG ,FSA ,FSR ,FIRA ,FSH ,FGEV , & !out
TRAD ,T2M ,SSOIL ,LATHEA ,Q2E ,EMISSI, CH2B ) !out
SICE = MAX(0.0, SMC - SH2O)
SNEQVO = SNEQV
QVAP = MAX( FGEV/LATHEA, 0.) ! positive part of fgev [mm/s] > 0
QDEW = ABS( MIN(FGEV/LATHEA, 0.)) ! negative part of fgev [mm/s] > 0
EDIR = QVAP - QDEW
! compute water budgets (water storages, ET components, and runoff)
CALL WATER_GLACIER (NSNOW ,NSOIL ,IMELT ,DT ,PRCP ,SFCTMP , & !in
QVAP ,QDEW ,FICEOLD,ZSOIL , & !in
ISNOW ,SNOWH ,SNEQV ,SNICE ,SNLIQ ,STC , & !inout
DZSNSO ,SH2O ,SICE ,PONDING,ZSNSO , & !inout
RUNSRF ,RUNSUB ,QSNOW ,PONDING1 ,PONDING2,QSNBOT,FPICE) !out
! water and energy balance check
CALL ERROR_GLACIER (ILOC ,JLOC ,SWDOWN ,FSA ,FSR ,FIRA , &
FSH ,FGEV ,SSOIL ,SAG ,PRCP ,EDIR , &
RUNSRF ,RUNSUB ,SNEQV ,DT ,BEG_WB )
IF(SWDOWN.NE.0.) THEN
ALBEDO = FSR / SWDOWN
ELSE
ALBEDO = -999.9
END IF
END SUBROUTINE NOAHMP_GLACIER
! ==================================================================================================
SUBROUTINE ATM_GLACIER (SFCPRS ,SFCTMP ,Q2 ,SOLDN ,COSZ ,THAIR , & 1
QAIR ,EAIR ,RHOAIR ,SOLAD ,SOLAI , &
SWDOWN )
! --------------------------------------------------------------------------------------------------
! re-process atmospheric forcing
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
! inputs
REAL , INTENT(IN) :: SFCPRS !pressure (pa)
REAL , INTENT(IN) :: SFCTMP !surface air temperature [k]
REAL , INTENT(IN) :: Q2 !mixing ratio (kg/kg)
REAL , INTENT(IN) :: SOLDN !downward shortwave radiation (w/m2)
REAL , INTENT(IN) :: COSZ !cosine solar zenith angle [0-1]
! outputs
REAL , INTENT(OUT) :: THAIR !potential temperature (k)
REAL , INTENT(OUT) :: QAIR !specific humidity (kg/kg) (q2/(1+q2))
REAL , INTENT(OUT) :: EAIR !vapor pressure air (pa)
REAL, DIMENSION( 1: 2), INTENT(OUT) :: SOLAD !incoming direct solar radiation (w/m2)
REAL, DIMENSION( 1: 2), INTENT(OUT) :: SOLAI !incoming diffuse solar radiation (w/m2)
REAL , INTENT(OUT) :: RHOAIR !density air (kg/m3)
REAL , INTENT(OUT) :: SWDOWN !downward solar filtered by sun angle [w/m2]
!locals
REAL :: PAIR !atm bottom level pressure (pa)
! --------------------------------------------------------------------------------------------------
PAIR = SFCPRS ! atm bottom level pressure (pa)
THAIR = SFCTMP * (SFCPRS/PAIR)**(RAIR/CPAIR)
QAIR = Q2 / (1.0+Q2) ! mixing ratio to specific humidity [kg/kg]
EAIR = QAIR*SFCPRS / (0.622+0.378*QAIR)
RHOAIR = (SFCPRS-0.378*EAIR) / (RAIR*SFCTMP)
IF(COSZ <= 0.) THEN
SWDOWN = 0.
ELSE
SWDOWN = SOLDN
END IF
SOLAD(1) = SWDOWN*0.7*0.5 ! direct vis
SOLAD(2) = SWDOWN*0.7*0.5 ! direct nir
SOLAI(1) = SWDOWN*0.3*0.5 ! diffuse vis
SOLAI(2) = SWDOWN*0.3*0.5 ! diffuse nir
END SUBROUTINE ATM_GLACIER
! ==================================================================================================
! --------------------------------------------------------------------------------------------------
SUBROUTINE ENERGY_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT ,QSNOW ,RHOAIR , & !in 1,6
EAIR ,SFCPRS ,QAIR ,SFCTMP ,LWDN ,UU , & !in
VV ,SOLAD ,SOLAI ,COSZ ,ZREF , & !in
TBOT ,ZBOT ,ZSNSO ,DZSNSO , & !in
TG ,STC ,SNOWH ,SNEQV ,SNEQVO ,SH2O , & !inout
SMC ,SNICE ,SNLIQ ,ALBOLD ,CM ,CH , & !inout
TAUSS ,QSFC , & !inout
IMELT ,SNICEV ,SNLIQV ,EPORE ,QMELT ,PONDING, & !out
SAG ,FSA ,FSR ,FIRA ,FSH ,FGEV , & !out
TRAD ,T2M ,SSOIL ,LATHEA ,Q2E ,EMISSI, CH2B ) !out
! --------------------------------------------------------------------------------------------------
! --------------------------------------------------------------------------------------------------
! USE NOAHMP_VEG_PARAMETERS
! USE NOAHMP_RAD_PARAMETERS
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
! inputs
INTEGER , INTENT(IN) :: NSNOW !maximum no. of snow layers
INTEGER , INTENT(IN) :: NSOIL !number of soil layers
INTEGER , INTENT(IN) :: ISNOW !actual no. of snow layers
REAL , INTENT(IN) :: DT !time step [sec]
REAL , INTENT(IN) :: QSNOW !snowfall on the ground (mm/s)
REAL , INTENT(IN) :: RHOAIR !density air (kg/m3)
REAL , INTENT(IN) :: EAIR !vapor pressure air (pa)
REAL , INTENT(IN) :: SFCPRS !pressure (pa)
REAL , INTENT(IN) :: QAIR !specific humidity (kg/kg)
REAL , INTENT(IN) :: SFCTMP !air temperature (k)
REAL , INTENT(IN) :: LWDN !downward longwave radiation (w/m2)
REAL , INTENT(IN) :: UU !wind speed in e-w dir (m/s)
REAL , INTENT(IN) :: VV !wind speed in n-s dir (m/s)
REAL , DIMENSION( 1: 2), INTENT(IN) :: SOLAD !incoming direct solar rad. (w/m2)
REAL , DIMENSION( 1: 2), INTENT(IN) :: SOLAI !incoming diffuse solar rad. (w/m2)
REAL , INTENT(IN) :: COSZ !cosine solar zenith angle (0-1)
REAL , INTENT(IN) :: ZREF !reference height (m)
REAL , INTENT(IN) :: TBOT !bottom condition for soil temp. (k)
REAL , INTENT(IN) :: ZBOT !depth for TBOT [m]
REAL , DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: ZSNSO !layer-bottom depth from snow surf [m]
REAL , DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO !depth of snow & soil layer-bottom [m]
! input & output
REAL , INTENT(INOUT) :: TG !ground temperature (k)
REAL , DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow/soil temperature [k]
REAL , INTENT(INOUT) :: SNOWH !snow height [m]
REAL , INTENT(INOUT) :: SNEQV !snow mass (mm)
REAL , INTENT(INOUT) :: SNEQVO !snow mass at last time step (mm)
REAL , DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O !liquid soil moisture [m3/m3]
REAL , DIMENSION( 1:NSOIL), INTENT(INOUT) :: SMC !soil moisture (ice + liq.) [m3/m3]
REAL , DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE !snow ice mass (kg/m2)
REAL , DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ !snow liq mass (kg/m2)
REAL , INTENT(INOUT) :: ALBOLD !snow albedo at last time step(CLASS type)
REAL , INTENT(INOUT) :: CM !momentum drag coefficient
REAL , INTENT(INOUT) :: CH !sensible heat exchange coefficient
REAL , INTENT(INOUT) :: TAUSS !snow aging factor
REAL , INTENT(INOUT) :: QSFC !mixing ratio at lowest model layer
! outputs
INTEGER, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: IMELT !phase change index [1-melt; 2-freeze]
REAL , DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNICEV !partial volume ice [m3/m3]
REAL , DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNLIQV !partial volume liq. water [m3/m3]
REAL , DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: EPORE !effective porosity [m3/m3]
REAL , INTENT(OUT) :: QMELT !snowmelt [mm/s]
REAL , INTENT(OUT) :: PONDING!pounding at ground [mm]
REAL , INTENT(OUT) :: SAG !solar rad. absorbed by ground (w/m2)
REAL , INTENT(OUT) :: FSA !tot. absorbed solar radiation (w/m2)
REAL , INTENT(OUT) :: FSR !tot. reflected solar radiation (w/m2)
REAL , INTENT(OUT) :: FIRA !total net LW. rad (w/m2) [+ to atm]
REAL , INTENT(OUT) :: FSH !total sensible heat (w/m2) [+ to atm]
REAL , INTENT(OUT) :: FGEV !ground evaporation (w/m2) [+ to atm]
REAL , INTENT(OUT) :: TRAD !radiative temperature (k)
REAL , INTENT(OUT) :: T2M !2 m height air temperature (k)
REAL , INTENT(OUT) :: SSOIL !ground heat flux (w/m2) [+ to soil]
REAL , INTENT(OUT) :: LATHEA !latent heat vap./sublimation (j/kg)
REAL , INTENT(OUT) :: Q2E
REAL , INTENT(OUT) :: EMISSI
REAL , INTENT(OUT) :: CH2B !sensible heat conductance, canopy air to ZLVL air (m/s)
! local
REAL :: UR !wind speed at height ZLVL (m/s)
REAL :: ZLVL !reference height (m)
REAL :: RSURF !ground surface resistance (s/m)
REAL :: ZPD !zero plane displacement (m)
REAL :: Z0MG !z0 momentum, ground (m)
REAL :: EMG !ground emissivity
REAL :: FIRE !emitted IR (w/m2)
REAL, DIMENSION(-NSNOW+1:NSOIL) :: FACT !temporary used in phase change
REAL, DIMENSION(-NSNOW+1:NSOIL) :: DF !thermal conductivity [w/m/k]
REAL, DIMENSION(-NSNOW+1:NSOIL) :: HCPCT !heat capacity [j/m3/k]
REAL :: GAMMA !psychrometric constant (pa/k)
REAL :: RHSUR !raltive humidity in surface soil/snow air space (-)
! ---------------------------------------------------------------------------------------------------
! wind speed at reference height: ur >= 1
UR = MAX( SQRT(UU**2.+VV**2.), 1. )
! roughness length and displacement height
Z0MG = Z0SNO
ZPD = SNOWH
ZLVL = ZPD + ZREF
! Thermal properties of soil, snow, lake, and frozen soil
CALL THERMOPROP_GLACIER (NSOIL ,NSNOW ,ISNOW ,DZSNSO , & !in
DT ,SNOWH ,SNICE ,SNLIQ , & !in
DF ,HCPCT ,SNICEV ,SNLIQV ,EPORE , & !out
FACT ) !out
! Solar radiation: absorbed & reflected by the ground
CALL RADIATION_GLACIER (DT ,TG ,SNEQVO ,SNEQV ,COSZ , & !in
QSNOW ,SOLAD ,SOLAI , & !in
ALBOLD ,TAUSS , & !inout
SAG ,FSR ,FSA) !out
! vegetation and ground emissivity
EMG = 0.98
! soil surface resistance for ground evap.
RHSUR = 1.0
RSURF = 1.0
! set psychrometric constant
IF (SFCTMP .GT. TFRZ) THEN
LATHEA = HVAP
ELSE
LATHEA = HSUB
END IF
GAMMA = CPAIR*SFCPRS/(0.622*LATHEA)
! Surface temperatures of the ground and energy fluxes
CALL GLACIER_FLUX (NSOIL ,NSNOW ,EMG ,ISNOW ,DF ,DZSNSO ,Z0MG , & !in
ZLVL ,ZPD ,QAIR ,SFCTMP ,RHOAIR ,SFCPRS , & !in
UR ,GAMMA ,RSURF ,LWDN ,RHSUR , & !in
EAIR ,STC ,SAG ,SNOWH ,LATHEA , & !in
CM ,CH ,TG ,QSFC , & !inout
FIRA ,FSH ,FGEV ,SSOIL , & !out
T2M ,Q2E ,CH2B) !out
!energy balance at surface: SAG=(IRB+SHB+EVB+GHB)
FIRE = LWDN + FIRA
IF(FIRE <=0.) call wrf_error_fatal("STOP in Noah-MP: emitted longwave <0")
! Compute a net emissivity
EMISSI = EMG
! When we're computing a TRAD, subtract from the emitted IR the
! reflected portion of the incoming LWDN, so we're just
! considering the IR originating in the canopy/ground system.
TRAD = ( ( FIRE - (1-EMISSI)*LWDN ) / (EMISSI*SB) ) ** 0.25
! 3L snow & 4L soil temperatures
CALL TSNOSOI_GLACIER (NSOIL ,NSNOW ,ISNOW ,DT ,TBOT , & !in
SSOIL ,SNOWH ,ZBOT ,ZSNSO ,DF , & !in
HCPCT , & !in
STC ) !inout
! adjusting snow surface temperature
IF(OPT_STC == 2) THEN
IF (SNOWH > 0.05 .AND. TG > TFRZ) TG = TFRZ
END IF
! Energy released or consumed by snow & frozen soil
CALL PHASECHANGE_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT ,FACT , & !in
DZSNSO , & !in
STC ,SNICE ,SNLIQ ,SNEQV ,SNOWH , & !inout
SMC ,SH2O , & !inout
QMELT ,IMELT ,PONDING ) !out
END SUBROUTINE ENERGY_GLACIER
! ==================================================================================================
SUBROUTINE THERMOPROP_GLACIER (NSOIL ,NSNOW ,ISNOW ,DZSNSO , & !in 1,1
DT ,SNOWH ,SNICE ,SNLIQ , & !in
DF ,HCPCT ,SNICEV ,SNLIQV ,EPORE , & !out
FACT ) !out
! -------------------------------------------------------------------------------------------------
! -------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
! inputs
INTEGER , INTENT(IN) :: NSOIL !number of soil layers
INTEGER , INTENT(IN) :: NSNOW !maximum no. of snow layers
INTEGER , INTENT(IN) :: ISNOW !actual no. of snow layers
REAL , INTENT(IN) :: DT !time step [s]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNICE !snow ice mass (kg/m2)
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNLIQ !snow liq mass (kg/m2)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO !thickness of snow/soil layers [m]
REAL , INTENT(IN) :: SNOWH !snow height [m]
! outputs
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: DF !thermal conductivity [w/m/k]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: HCPCT !heat capacity [j/m3/k]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNICEV !partial volume of ice [m3/m3]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNLIQV !partial volume of liquid water [m3/m3]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: EPORE !effective porosity [m3/m3]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: FACT !computing energy for phase change
! --------------------------------------------------------------------------------------------------
! locals
INTEGER :: IZ, IZ2
REAL, DIMENSION(-NSNOW+1: 0) :: CVSNO !volumetric specific heat (j/m3/k)
REAL, DIMENSION(-NSNOW+1: 0) :: TKSNO !snow thermal conductivity (j/m3/k)
REAL :: ZMID !mid-point soil depth
! --------------------------------------------------------------------------------------------------
! compute snow thermal conductivity and heat capacity
CALL CSNOW_GLACIER (ISNOW ,NSNOW ,NSOIL ,SNICE ,SNLIQ ,DZSNSO , & !in
TKSNO ,CVSNO ,SNICEV ,SNLIQV ,EPORE ) !out
DO IZ = ISNOW+1, 0
DF (IZ) = TKSNO(IZ)
HCPCT(IZ) = CVSNO(IZ)
END DO
! compute soil thermal properties (using Noah glacial ice approximations)
DO IZ = 1, NSOIL
ZMID = 0.5 * (DZSNSO(IZ))
DO IZ2 = 1, IZ-1
ZMID = ZMID + DZSNSO(IZ2)
END DO
HCPCT(IZ) = 1.E6 * ( 0.8194 - 0.1309*ZMID )
DF(IZ) = 0.32333 + ( 0.10073 * ZMID )
END DO
! combine a temporary variable used for melting/freezing of snow and frozen soil
DO IZ = ISNOW+1,NSOIL
FACT(IZ) = DT/(HCPCT(IZ)*DZSNSO(IZ))
END DO
! snow/soil interface
IF(ISNOW == 0) THEN
DF(1) = (DF(1)*DZSNSO(1)+0.35*SNOWH) / (SNOWH +DZSNSO(1))
ELSE
DF(1) = (DF(1)*DZSNSO(1)+DF(0)*DZSNSO(0)) / (DZSNSO(0)+DZSNSO(1))
END IF
END SUBROUTINE THERMOPROP_GLACIER
! ==================================================================================================
! --------------------------------------------------------------------------------------------------
SUBROUTINE CSNOW_GLACIER (ISNOW ,NSNOW ,NSOIL ,SNICE ,SNLIQ ,DZSNSO , & !in 1
TKSNO ,CVSNO ,SNICEV ,SNLIQV ,EPORE ) !out
! --------------------------------------------------------------------------------------------------
! Snow bulk density,volumetric capacity, and thermal conductivity
!---------------------------------------------------------------------------------------------------
IMPLICIT NONE
!---------------------------------------------------------------------------------------------------
! inputs
INTEGER, INTENT(IN) :: ISNOW !number of snow layers (-)
INTEGER , INTENT(IN) :: NSNOW !maximum no. of snow layers
INTEGER , INTENT(IN) :: NSOIL !number of soil layers
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNICE !snow ice mass (kg/m2)
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNLIQ !snow liq mass (kg/m2)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO !snow/soil layer thickness [m]
! outputs
REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: CVSNO !volumetric specific heat (j/m3/k)
REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: TKSNO !thermal conductivity (w/m/k)
REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNICEV !partial volume of ice [m3/m3]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNLIQV !partial volume of liquid water [m3/m3]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: EPORE !effective porosity [m3/m3]
! locals
INTEGER :: IZ
REAL, DIMENSION(-NSNOW+1: 0) :: BDSNOI !bulk density of snow(kg/m3)
!---------------------------------------------------------------------------------------------------
! thermal capacity of snow
DO IZ = ISNOW+1, 0
SNICEV(IZ) = MIN(1., SNICE(IZ)/(DZSNSO(IZ)*DENICE) )
EPORE(IZ) = 1. - SNICEV(IZ)
SNLIQV(IZ) = MIN(EPORE(IZ),SNLIQ(IZ)/(DZSNSO(IZ)*DENH2O))
ENDDO
DO IZ = ISNOW+1, 0
BDSNOI(IZ) = (SNICE(IZ)+SNLIQ(IZ))/DZSNSO(IZ)
CVSNO(IZ) = CICE*SNICEV(IZ)+CWAT*SNLIQV(IZ)
! CVSNO(IZ) = 0.525E06 ! constant
enddo
! thermal conductivity of snow
DO IZ = ISNOW+1, 0
TKSNO(IZ) = 3.2217E-6*BDSNOI(IZ)**2. ! Stieglitz(yen,1965)
! TKSNO(IZ) = 2E-2+2.5E-6*BDSNOI(IZ)*BDSNOI(IZ) ! Anderson, 1976
! TKSNO(IZ) = 0.35 ! constant
! TKSNO(IZ) = 2.576E-6*BDSNOI(IZ)**2. + 0.074 ! Verseghy (1991)
! TKSNO(IZ) = 2.22*(BDSNOI(IZ)/1000.)**1.88 ! Douvill(Yen, 1981)
ENDDO
END SUBROUTINE CSNOW_GLACIER
!===================================================================================================
SUBROUTINE RADIATION_GLACIER (DT ,TG ,SNEQVO ,SNEQV ,COSZ , & !in 1,3
QSNOW ,SOLAD ,SOLAI , & !in
ALBOLD ,TAUSS , & !inout
SAG ,FSR ,FSA) !out
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
! input
REAL, INTENT(IN) :: DT !time step [s]
REAL, INTENT(IN) :: TG !ground temperature (k)
REAL, INTENT(IN) :: SNEQVO !snow mass at last time step(mm)
REAL, INTENT(IN) :: SNEQV !snow mass (mm)
REAL, INTENT(IN) :: COSZ !cosine solar zenith angle (0-1)
REAL, INTENT(IN) :: QSNOW !snowfall (mm/s)
REAL, DIMENSION(1:2) , INTENT(IN) :: SOLAD !incoming direct solar radiation (w/m2)
REAL, DIMENSION(1:2) , INTENT(IN) :: SOLAI !incoming diffuse solar radiation (w/m2)
! inout
REAL, INTENT(INOUT) :: ALBOLD !snow albedo at last time step (CLASS type)
REAL, INTENT(INOUT) :: TAUSS !non-dimensional snow age
! output
REAL, INTENT(OUT) :: SAG !solar radiation absorbed by ground (w/m2)
REAL, INTENT(OUT) :: FSR !total reflected solar radiation (w/m2)
REAL, INTENT(OUT) :: FSA !total absorbed solar radiation (w/m2)
! local
INTEGER :: IB !number of radiation bands
INTEGER :: NBAND !number of radiation bands
REAL :: FAGE !snow age function (0 - new snow)
REAL, DIMENSION(1:2) :: ALBSND !snow albedo (direct)
REAL, DIMENSION(1:2) :: ALBSNI !snow albedo (diffuse)
REAL :: ALB !current CLASS albedo
REAL :: ABS !temporary absorbed rad
REAL :: REF !temporary reflected rad
REAL :: FSNO !snow-cover fraction, = 1 if any snow
REAL, DIMENSION(1:2) :: ALBICE !albedo land ice: 1=vis, 2=nir
REAL,PARAMETER :: MPE = 1.E-6
! --------------------------------------------------------------------------------------------------
NBAND = 2
ALBSND = 0.0
ALBSNI = 0.0
ALBICE(1) = 0.80 !albedo land ice: 1=vis, 2=nir
ALBICE(2) = 0.55
! snow age
CALL SNOW_AGE_GLACIER (DT,TG,SNEQVO,SNEQV,TAUSS,FAGE)
! snow albedos: age even when sun is not present
IF(OPT_ALB == 1) &
CALL SNOWALB_BATS_GLACIER (NBAND,COSZ,FAGE,ALBSND,ALBSNI)
IF(OPT_ALB == 2) THEN
CALL SNOWALB_CLASS_GLACIER(NBAND,QSNOW,DT,ALB,ALBOLD,ALBSND,ALBSNI)
ALBOLD = ALB
END IF
! zero summed solar fluxes
SAG = 0.
FSA = 0.
FSR = 0.
FSNO = 0.0
IF(SNEQV > 0.0) FSNO = 1.0
! loop over nband wavebands
DO IB = 1, NBAND
ALBSND(IB) = ALBICE(IB)*(1.-FSNO) + ALBSND(IB)*FSNO
ALBSNI(IB) = ALBICE(IB)*(1.-FSNO) + ALBSNI(IB)*FSNO
! solar radiation absorbed by ground surface
ABS = SOLAD(IB)*(1.-ALBSND(IB)) + SOLAI(IB)*(1.-ALBSNI(IB))
SAG = SAG + ABS
FSA = FSA + ABS
REF = SOLAD(IB)*ALBSND(IB) + SOLAI(IB)*ALBSNI(IB)
FSR = FSR + REF
END DO
END SUBROUTINE RADIATION_GLACIER
! ==================================================================================================
SUBROUTINE SNOW_AGE_GLACIER (DT,TG,SNEQVO,SNEQV,TAUSS,FAGE) 1
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! ------------------------ code history ------------------------------------------------------------
! from BATS
! ------------------------ input/output variables --------------------------------------------------
!input
REAL, INTENT(IN) :: DT !main time step (s)
REAL, INTENT(IN) :: TG !ground temperature (k)
REAL, INTENT(IN) :: SNEQVO !snow mass at last time step(mm)
REAL, INTENT(IN) :: SNEQV !snow water per unit ground area (mm)
! inout
REAL, INTENT(INOUT) :: TAUSS !non-dimensional snow age
!output
REAL, INTENT(OUT) :: FAGE !snow age
!local
REAL :: TAGE !total aging effects
REAL :: AGE1 !effects of grain growth due to vapor diffusion
REAL :: AGE2 !effects of grain growth at freezing of melt water
REAL :: AGE3 !effects of soot
REAL :: DELA !temporary variable
REAL :: SGE !temporary variable
REAL :: DELS !temporary variable
REAL :: DELA0 !temporary variable
REAL :: ARG !temporary variable
! See Yang et al. (1997) J.of Climate for detail.
!---------------------------------------------------------------------------------------------------
IF(SNEQV.LE.0.0) THEN
TAUSS = 0.
ELSE IF (SNEQV.GT.800.) THEN
TAUSS = 0.
ELSE
TAUSS = 0.
DELA0 = 1.E-6*DT
ARG = 5.E3*(1./TFRZ-1./TG)
AGE1 = EXP(ARG)
AGE2 = EXP(AMIN1(0.,10.*ARG))
AGE3 = 0.3
TAGE = AGE1+AGE2+AGE3
DELA = DELA0*TAGE
DELS = AMAX1(0.0,SNEQV-SNEQVO) / SWEMX
SGE = (TAUSS+DELA)*(1.0-DELS)
TAUSS = AMAX1(0.,SGE)
ENDIF
FAGE= TAUSS/(TAUSS+1.)
END SUBROUTINE SNOW_AGE_GLACIER
! ==================================================================================================
! --------------------------------------------------------------------------------------------------
SUBROUTINE SNOWALB_BATS_GLACIER (NBAND,COSZ,FAGE,ALBSND,ALBSNI) 1
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
! input
INTEGER,INTENT(IN) :: NBAND !number of waveband classes
REAL,INTENT(IN) :: COSZ !cosine solar zenith angle
REAL,INTENT(IN) :: FAGE !snow age correction
! output
REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSND !snow albedo for direct(1=vis, 2=nir)
REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSNI !snow albedo for diffuse
! ---------------------------------------------------------------------------------------------
REAL :: FZEN !zenith angle correction
REAL :: CF1 !temperary variable
REAL :: SL2 !2.*SL
REAL :: SL1 !1/SL
REAL :: SL !adjustable parameter
REAL, PARAMETER :: C1 = 0.2 !default in BATS
REAL, PARAMETER :: C2 = 0.5 !default in BATS
! REAL, PARAMETER :: C1 = 0.2 * 2. ! double the default to match Sleepers River's
! REAL, PARAMETER :: C2 = 0.5 * 2. ! snow surface albedo (double aging effects)
! ---------------------------------------------------------------------------------------------
! zero albedos for all points
ALBSND(1: NBAND) = 0.
ALBSNI(1: NBAND) = 0.
! when cosz > 0
SL=2.0
SL1=1./SL
SL2=2.*SL
CF1=((1.+SL1)/(1.+SL2*COSZ)-SL1)
FZEN=AMAX1(CF1,0.)
ALBSNI(1)=0.95*(1.-C1*FAGE)
ALBSNI(2)=0.65*(1.-C2*FAGE)
ALBSND(1)=ALBSNI(1)+0.4*FZEN*(1.-ALBSNI(1)) ! vis direct
ALBSND(2)=ALBSNI(2)+0.4*FZEN*(1.-ALBSNI(2)) ! nir direct
END SUBROUTINE SNOWALB_BATS_GLACIER
! ==================================================================================================
! --------------------------------------------------------------------------------------------------
SUBROUTINE SNOWALB_CLASS_GLACIER (NBAND,QSNOW,DT,ALB,ALBOLD,ALBSND,ALBSNI) 1
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
! input
INTEGER,INTENT(IN) :: NBAND !number of waveband classes
REAL,INTENT(IN) :: QSNOW !snowfall (mm/s)
REAL,INTENT(IN) :: DT !time step (sec)
REAL,INTENT(IN) :: ALBOLD !snow albedo at last time step
! in & out
REAL, INTENT(INOUT) :: ALB !
! output
REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSND !snow albedo for direct(1=vis, 2=nir)
REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSNI !snow albedo for diffuse
! ---------------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------------
! zero albedos for all points
ALBSND(1: NBAND) = 0.
ALBSNI(1: NBAND) = 0.
! when cosz > 0
ALB = 0.55 + (ALBOLD-0.55) * EXP(-0.01*DT/3600.)
! 1 mm fresh snow(SWE) -- 10mm snow depth, assumed the fresh snow density 100kg/m3
! here assume 1cm snow depth will fully cover the old snow
IF (QSNOW > 0.) then
ALB = ALB + MIN(QSNOW*DT,SWEMX) * (0.84-ALB)/(SWEMX)
ENDIF
ALBSNI(1)= ALB ! vis diffuse
ALBSNI(2)= ALB ! nir diffuse
ALBSND(1)= ALB ! vis direct
ALBSND(2)= ALB ! nir direct
END SUBROUTINE SNOWALB_CLASS_GLACIER
! ==================================================================================================
SUBROUTINE GLACIER_FLUX (NSOIL ,NSNOW ,EMG ,ISNOW ,DF ,DZSNSO ,Z0M , & !in 1,3
ZLVL ,ZPD ,QAIR ,SFCTMP ,RHOAIR ,SFCPRS , & !in
UR ,GAMMA ,RSURF ,LWDN ,RHSUR , & !in
EAIR ,STC ,SAG ,SNOWH ,LATHEA , & !in
CM ,CH ,TGB ,QSFC , & !inout
IRB ,SHB ,EVB ,GHB , & !out
T2MB ,Q2B ,EHB2) !out
! --------------------------------------------------------------------------------------------------
! use newton-raphson iteration to solve ground (tg) temperature
! that balances the surface energy budgets for glacier.
! bare soil:
! -SAB + IRB[TG] + SHB[TG] + EVB[TG] + GHB[TG] = 0
! ----------------------------------------------------------------------
! USE MODULE_MODEL_CONSTANTS
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers
INTEGER, INTENT(IN) :: NSOIL !number of soil layers
REAL, INTENT(IN) :: EMG !ground emissivity
INTEGER, INTENT(IN) :: ISNOW !actual no. of snow layers
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DF !thermal conductivity of snow/soil (w/m/k)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO !thickness of snow/soil layers (m)
REAL, INTENT(IN) :: Z0M !roughness length, momentum, ground (m)
REAL, INTENT(IN) :: ZLVL !reference height (m)
REAL, INTENT(IN) :: ZPD !zero plane displacement (m)
REAL, INTENT(IN) :: QAIR !specific humidity at height zlvl (kg/kg)
REAL, INTENT(IN) :: SFCTMP !air temperature at reference height (k)
REAL, INTENT(IN) :: RHOAIR !density air (kg/m3)
REAL, INTENT(IN) :: SFCPRS !density air (kg/m3)
REAL, INTENT(IN) :: UR !wind speed at height zlvl (m/s)
REAL, INTENT(IN) :: GAMMA !psychrometric constant (pa/k)
REAL, INTENT(IN) :: RSURF !ground surface resistance (s/m)
REAL, INTENT(IN) :: LWDN !atmospheric longwave radiation (w/m2)
REAL, INTENT(IN) :: RHSUR !raltive humidity in surface soil/snow air space (-)
REAL, INTENT(IN) :: EAIR !vapor pressure air at height (pa)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: STC !soil/snow temperature (k)
REAL, INTENT(IN) :: SAG !solar radiation absorbed by ground (w/m2)
REAL, INTENT(IN) :: SNOWH !actual snow depth [m]
REAL, INTENT(IN) :: LATHEA !latent heat of vaporization/subli (j/kg)
! input/output
REAL, INTENT(INOUT) :: CM !momentum drag coefficient
REAL, INTENT(INOUT) :: CH !sensible heat exchange coefficient
REAL, INTENT(INOUT) :: TGB !ground temperature (k)
REAL, INTENT(INOUT) :: QSFC !mixing ratio at lowest model layer
! output
! -SAB + IRB[TG] + SHB[TG] + EVB[TG] + GHB[TG] = 0
REAL, INTENT(OUT) :: IRB !net longwave rad (w/m2) [+ to atm]
REAL, INTENT(OUT) :: SHB !sensible heat flux (w/m2) [+ to atm]
REAL, INTENT(OUT) :: EVB !latent heat flux (w/m2) [+ to atm]
REAL, INTENT(OUT) :: GHB !ground heat flux (w/m2) [+ to soil]
REAL, INTENT(OUT) :: T2MB !2 m height air temperature (k)
REAL, INTENT(OUT) :: Q2B !bare ground heat conductance
REAL, INTENT(OUT) :: EHB2 !sensible heat conductance for diagnostics
! local variables
INTEGER :: NITERB !number of iterations for surface temperature
REAL :: MPE !prevents overflow error if division by zero
REAL :: DTG !change in tg, last iteration (k)
INTEGER :: MOZSGN !number of times MOZ changes sign
REAL :: MOZOLD !Monin-Obukhov stability parameter from prior iteration
REAL :: H !temporary sensible heat flux (w/m2)
REAL :: FV !friction velocity (m/s)
REAL :: CIR !coefficients for ir as function of ts**4
REAL :: CGH !coefficients for st as function of ts
REAL :: CSH !coefficients for sh as function of ts
REAL :: CEV !coefficients for ev as function of esat[ts]
REAL :: CQ2B !
INTEGER :: ITER !iteration index
REAL :: Z0H !roughness length, sensible heat, ground (m)
REAL :: MOZ !Monin-Obukhov stability parameter
REAL :: FM !momentum stability correction, weighted by prior iters
REAL :: FH !sen heat stability correction, weighted by prior iters
REAL :: RAMB !aerodynamic resistance for momentum (s/m)
REAL :: RAHB !aerodynamic resistance for sensible heat (s/m)
REAL :: RAWB !aerodynamic resistance for water vapor (s/m)
REAL :: ESTG !saturation vapor pressure at tg (pa)
REAL :: DESTG !d(es)/dt at tg (pa/K)
REAL :: ESATW !es for water
REAL :: ESATI !es for ice
REAL :: DSATW !d(es)/dt at tg (pa/K) for water
REAL :: DSATI !d(es)/dt at tg (pa/K) for ice
REAL :: A !temporary calculation
REAL :: B !temporary calculation
REAL :: T, TDC !Kelvin to degree Celsius with limit -50 to +50
TDC(T) = MIN( 50., MAX(-50.,(T-TFRZ)) )
! -----------------------------------------------------------------
! initialization variables that do not depend on stability iteration
! -----------------------------------------------------------------
NITERB = 5
MPE = 1E-6
DTG = 0.
MOZSGN = 0
MOZOLD = 0.
H = 0.
FV = 0.1
CIR = EMG*SB
CGH = 2.*DF(ISNOW+1)/DZSNSO(ISNOW+1)
! -----------------------------------------------------------------
loop3: DO ITER = 1, NITERB ! begin stability iteration
Z0H = Z0M
! For now, only allow SFCDIF1 until others can be fixed
CALL SFCDIF1_GLACIER(ITER ,ZLVL ,ZPD ,Z0H ,Z0M , & !in
QAIR ,SFCTMP ,H ,RHOAIR ,MPE ,UR , & !in
& MOZ ,MOZSGN ,FM ,FH , & !inout
& FV ,CM ,CH ) !out
RAMB = MAX(1.,1./(CM*UR))
RAHB = MAX(1.,1./(CH*UR))
RAWB = RAHB
! es and d(es)/dt evaluated at tg
T = TDC(TGB)
CALL ESAT(T, ESATW, ESATI, DSATW, DSATI)
IF (T .GT. 0.) THEN
ESTG = ESATW
DESTG = DSATW
ELSE
ESTG = ESATI
DESTG = DSATI
END IF
CSH = RHOAIR*CPAIR/RAHB
CEV = RHOAIR*CPAIR/GAMMA/(RSURF+RAWB)
! surface fluxes and dtg
IRB = CIR * TGB**4 - EMG*LWDN
SHB = CSH * (TGB - SFCTMP )
EVB = CEV * (ESTG*RHSUR - EAIR )
GHB = CGH * (TGB - STC(ISNOW+1))
B = SAG-IRB-SHB-EVB-GHB
A = 4.*CIR*TGB**3 + CSH + CEV*DESTG + CGH
DTG = B/A
IRB = IRB + 4.*CIR*TGB**3*DTG
SHB = SHB + CSH*DTG
EVB = EVB + CEV*DESTG*DTG
GHB = GHB + CGH*DTG
! update ground surface temperature
TGB = TGB + DTG
! for M-O length
H = CSH * (TGB - SFCTMP)
T = TDC(TGB)
CALL ESAT(T, ESATW, ESATI, DSATW, DSATI)
IF (T .GT. 0.) THEN
ESTG = ESATW
ELSE
ESTG = ESATI
END IF
QSFC = 0.622*(ESTG*RHSUR)/(SFCPRS-0.378*(ESTG*RHSUR))
END DO loop3 ! end stability iteration
! -----------------------------------------------------------------
! if snow on ground and TG > TFRZ: reset TG = TFRZ. reevaluate ground fluxes.
IF(OPT_STC == 1) THEN
IF (SNOWH > 0.05 .AND. TGB > TFRZ) THEN
TGB = TFRZ
IRB = CIR * TGB**4 - EMG*LWDN
SHB = CSH * (TGB - SFCTMP)
EVB = CEV * (ESTG*RHSUR - EAIR ) !ESTG reevaluate ?
GHB = SAG - (IRB+SHB+EVB)
END IF
END IF
! 2m air temperature
EHB2 = FV*VKC/LOG((2.+Z0H)/Z0H)
CQ2B = EHB2
IF (EHB2.lt.1.E-5 ) THEN
T2MB = TGB
Q2B = QSFC
ELSE
T2MB = TGB - SHB/(RHOAIR*CPAIR*FV) * 1./VKC * LOG((2.+Z0H)/Z0H)
Q2B = QSFC - EVB/(LATHEA*RHOAIR)*(1./CQ2B + RSURF)
ENDIF
! update CH
CH = 1./RAHB
END SUBROUTINE GLACIER_FLUX
! ==================================================================================================
SUBROUTINE ESAT(T, ESW, ESI, DESW, DESI) 7
!---------------------------------------------------------------------------------------------------
! use polynomials to calculate saturation vapor pressure and derivative with
! respect to temperature: over water when t > 0 c and over ice when t <= 0 c
IMPLICIT NONE
!---------------------------------------------------------------------------------------------------
! in
REAL, intent(in) :: T !temperature
!out
REAL, intent(out) :: ESW !saturation vapor pressure over water (pa)
REAL, intent(out) :: ESI !saturation vapor pressure over ice (pa)
REAL, intent(out) :: DESW !d(esat)/dt over water (pa/K)
REAL, intent(out) :: DESI !d(esat)/dt over ice (pa/K)
! local
REAL :: A0,A1,A2,A3,A4,A5,A6 !coefficients for esat over water
REAL :: B0,B1,B2,B3,B4,B5,B6 !coefficients for esat over ice
REAL :: C0,C1,C2,C3,C4,C5,C6 !coefficients for dsat over water
REAL :: D0,D1,D2,D3,D4,D5,D6 !coefficients for dsat over ice
PARAMETER (A0=6.107799961 , A1=4.436518521E-01, &
A2=1.428945805E-02, A3=2.650648471E-04, &
A4=3.031240396E-06, A5=2.034080948E-08, &
A6=6.136820929E-11)
PARAMETER (B0=6.109177956 , B1=5.034698970E-01, &
B2=1.886013408E-02, B3=4.176223716E-04, &
B4=5.824720280E-06, B5=4.838803174E-08, &
B6=1.838826904E-10)
PARAMETER (C0= 4.438099984E-01, C1=2.857002636E-02, &
C2= 7.938054040E-04, C3=1.215215065E-05, &
C4= 1.036561403E-07, C5=3.532421810e-10, &
C6=-7.090244804E-13)
PARAMETER (D0=5.030305237E-01, D1=3.773255020E-02, &
D2=1.267995369E-03, D3=2.477563108E-05, &
D4=3.005693132E-07, D5=2.158542548E-09, &
D6=7.131097725E-12)
ESW = 100.*(A0+T*(A1+T*(A2+T*(A3+T*(A4+T*(A5+T*A6))))))
ESI = 100.*(B0+T*(B1+T*(B2+T*(B3+T*(B4+T*(B5+T*B6))))))
DESW = 100.*(C0+T*(C1+T*(C2+T*(C3+T*(C4+T*(C5+T*C6))))))
DESI = 100.*(D0+T*(D1+T*(D2+T*(D3+T*(D4+T*(D5+T*D6))))))
END SUBROUTINE ESAT
! ==================================================================================================
SUBROUTINE SFCDIF1_GLACIER(ITER ,ZLVL ,ZPD ,Z0H ,Z0M , & !in 1,1
QAIR ,SFCTMP ,H ,RHOAIR ,MPE ,UR , & !in
& MOZ ,MOZSGN ,FM ,FH , & !inout
& FV ,CM ,CH ) !out
! -------------------------------------------------------------------------------------------------
! computing surface drag coefficient CM for momentum and CH for heat
! -------------------------------------------------------------------------------------------------
IMPLICIT NONE
! -------------------------------------------------------------------------------------------------
! inputs
INTEGER, INTENT(IN) :: ITER !iteration index
REAL, INTENT(IN) :: ZLVL !reference height (m)
REAL, INTENT(IN) :: ZPD !zero plane displacement (m)
REAL, INTENT(IN) :: Z0H !roughness length, sensible heat, ground (m)
REAL, INTENT(IN) :: Z0M !roughness length, momentum, ground (m)
REAL, INTENT(IN) :: QAIR !specific humidity at reference height (kg/kg)
REAL, INTENT(IN) :: SFCTMP !temperature at reference height (k)
REAL, INTENT(IN) :: H !sensible heat flux (w/m2) [+ to atm]
REAL, INTENT(IN) :: RHOAIR !density air (kg/m**3)
REAL, INTENT(IN) :: MPE !prevents overflow error if division by zero
REAL, INTENT(IN) :: UR !wind speed (m/s)
! in & out
REAL, INTENT(INOUT) :: MOZ !Monin-Obukhov stability (z/L)
INTEGER, INTENT(INOUT) :: MOZSGN !number of times moz changes sign
REAL, INTENT(INOUT) :: FM !momentum stability correction, weighted by prior iters
REAL, INTENT(INOUT) :: FH !sen heat stability correction, weighted by prior iters
! outputs
REAL, INTENT(OUT) :: FV !friction velocity (m/s)
REAL, INTENT(OUT) :: CM !drag coefficient for momentum
REAL, INTENT(OUT) :: CH !drag coefficient for heat
! locals
REAL :: MOZOLD !Monin-Obukhov stability parameter from prior iteration
REAL :: TMPCM !temporary calculation for CM
REAL :: TMPCH !temporary calculation for CH
REAL :: MOL !Monin-Obukhov length (m)
REAL :: TVIR !temporary virtual temperature (k)
REAL :: TMP1,TMP2,TMP3 !temporary calculation
REAL :: FMNEW !stability correction factor, momentum, for current moz
REAL :: FHNEW !stability correction factor, sen heat, for current moz
REAL :: CMFM, CHFH
! -------------------------------------------------------------------------------------------------
! Monin-Obukhov stability parameter moz for next iteration
MOZOLD = MOZ
IF(ZLVL <= ZPD) THEN
write(*,*) 'critical problem: ZLVL <= ZPD; model stops'
call wrf_error_fatal("STOP in Noah-MP")
ENDIF
TMPCM = LOG((ZLVL-ZPD) / Z0M)
TMPCH = LOG((ZLVL-ZPD) / Z0H)
IF(ITER == 1) THEN
FV = 0.0
MOZ = 0.0
MOL = 0.0
ELSE
TVIR = (1. + 0.61*QAIR) * SFCTMP
TMP1 = VKC * (GRAV/TVIR) * H/(RHOAIR*CPAIR)
IF (ABS(TMP1) .LE. MPE) TMP1 = MPE
MOL = -1. * FV**3 / TMP1
MOZ = MIN( (ZLVL-ZPD)/MOL, 1.)
ENDIF
! accumulate number of times moz changes sign.
IF (MOZOLD*MOZ .LT. 0.) MOZSGN = MOZSGN+1
IF (MOZSGN .GE. 2) THEN
MOZ = 0.
FM = 0.
FH = 0.
ENDIF
! evaluate stability-dependent variables using moz from prior iteration
IF (MOZ .LT. 0.) THEN
TMP1 = (1. - 16.*MOZ)**0.25
TMP2 = LOG((1.+TMP1*TMP1)/2.)
TMP3 = LOG((1.+TMP1)/2.)
FMNEW = 2.*TMP3 + TMP2 - 2.*ATAN(TMP1) + 1.5707963
FHNEW = 2*TMP2
ELSE
FMNEW = -5.*MOZ
FHNEW = FMNEW
ENDIF
! except for first iteration, weight stability factors for previous
! iteration to help avoid flip-flops from one iteration to the next
IF (ITER == 1) THEN
FM = FMNEW
FH = FHNEW
ELSE
FM = 0.5 * (FM+FMNEW)
FH = 0.5 * (FH+FHNEW)
ENDIF
! exchange coefficients
CMFM = TMPCM-FM
CHFH = TMPCH-FH
IF(ABS(CMFM) <= MPE) CMFM = MPE
IF(ABS(CHFH) <= MPE) CHFH = MPE
CM = VKC*VKC/(CMFM*CMFM)
CH = VKC*VKC/(CMFM*CHFH)
! friction velocity
FV = UR * SQRT(CM)
END SUBROUTINE SFCDIF1_GLACIER
! ==================================================================================================
SUBROUTINE TSNOSOI_GLACIER (NSOIL ,NSNOW ,ISNOW ,DT ,TBOT , & !in 1,2
SSOIL ,SNOWH ,ZBOT ,ZSNSO ,DF , & !in
HCPCT , & !in
STC ) !inout
! --------------------------------------------------------------------------------------------------
! Compute snow (up to 3L) and soil (4L) temperature. Note that snow temperatures
! during melting season may exceed melting point (TFRZ) but later in PHASECHANGE
! subroutine the snow temperatures are reset to TFRZ for melting snow.
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
!input
INTEGER, INTENT(IN) :: NSOIL !no of soil layers (4)
INTEGER, INTENT(IN) :: NSNOW !maximum no of snow layers (3)
INTEGER, INTENT(IN) :: ISNOW !actual no of snow layers
REAL, INTENT(IN) :: DT !time step (s)
REAL, INTENT(IN) :: TBOT !
REAL, INTENT(IN) :: SSOIL !ground heat flux (w/m2)
REAL, INTENT(IN) :: SNOWH !snow depth (m)
REAL, INTENT(IN) :: ZBOT !from soil surface (m)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: ZSNSO !layer-bot. depth from snow surf.(m)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DF !thermal conductivity
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: HCPCT !heat capacity (J/m3/k)
!input and output
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC
!local
INTEGER :: IZ
REAL :: ZBOTSNO !ZBOT from snow surface
REAL, DIMENSION(-NSNOW+1:NSOIL) :: AI, BI, CI, RHSTS
REAL :: EFLXB !energy influx from soil bottom (w/m2)
REAL, DIMENSION(-NSNOW+1:NSOIL) :: PHI !light through water (w/m2)
! ----------------------------------------------------------------------
! prescribe solar penetration into ice/snow
PHI(ISNOW+1:NSOIL) = 0.
! adjust ZBOT from soil surface to ZBOTSNO from snow surface
ZBOTSNO = ZBOT - SNOWH !from snow surface
! compute ice temperatures
CALL HRT_GLACIER (NSNOW ,NSOIL ,ISNOW ,ZSNSO , &
STC ,TBOT ,ZBOTSNO ,DF , &
HCPCT ,SSOIL ,PHI , &
AI ,BI ,CI ,RHSTS , &
EFLXB )
CALL HSTEP_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT , &
AI ,BI ,CI ,RHSTS , &
STC )
END SUBROUTINE TSNOSOI_GLACIER
! ==================================================================================================
! ----------------------------------------------------------------------
SUBROUTINE HRT_GLACIER (NSNOW ,NSOIL ,ISNOW ,ZSNSO , & !in 1
STC ,TBOT ,ZBOT ,DF , & !in
HCPCT ,SSOIL ,PHI , & !in
AI ,BI ,CI ,RHSTS , & !out
BOTFLX ) !out
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
! calculate the right hand side of the time tendency term of the soil
! thermal diffusion equation. also to compute ( prepare ) the matrix
! coefficients for the tri-diagonal matrix of the implicit time scheme.
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSOIL !no of soil layers (4)
INTEGER, INTENT(IN) :: NSNOW !maximum no of snow layers (3)
INTEGER, INTENT(IN) :: ISNOW !actual no of snow layers
REAL, INTENT(IN) :: TBOT !bottom soil temp. at ZBOT (k)
REAL, INTENT(IN) :: ZBOT !depth of lower boundary condition (m)
!from soil surface not snow surface
REAL, INTENT(IN) :: SSOIL !ground heat flux (w/m2)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: ZSNSO !depth of layer-bottom of snow/soil (m)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: STC !snow/soil temperature (k)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DF !thermal conductivity [w/m/k]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: HCPCT !heat capacity [j/m3/k]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: PHI !light through water (w/m2)
! output
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: RHSTS !right-hand side of the matrix
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: AI !left-hand side coefficient
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: BI !left-hand side coefficient
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: CI !left-hand side coefficient
REAL, INTENT(OUT) :: BOTFLX !energy influx from soil bottom (w/m2)
! local
INTEGER :: K
REAL, DIMENSION(-NSNOW+1:NSOIL) :: DDZ
REAL, DIMENSION(-NSNOW+1:NSOIL) :: DENOM
REAL, DIMENSION(-NSNOW+1:NSOIL) :: DTSDZ
REAL, DIMENSION(-NSNOW+1:NSOIL) :: EFLUX
REAL :: TEMP1
! ----------------------------------------------------------------------
DO K = ISNOW+1, NSOIL
IF (K == ISNOW+1) THEN
DENOM(K) = - ZSNSO(K) * HCPCT(K)
TEMP1 = - ZSNSO(K+1)
DDZ(K) = 2.0 / TEMP1
DTSDZ(K) = 2.0 * (STC(K) - STC(K+1)) / TEMP1
EFLUX(K) = DF(K) * DTSDZ(K) - SSOIL - PHI(K)
ELSE IF (K < NSOIL) THEN
DENOM(K) = (ZSNSO(K-1) - ZSNSO(K)) * HCPCT(K)
TEMP1 = ZSNSO(K-1) - ZSNSO(K+1)
DDZ(K) = 2.0 / TEMP1
DTSDZ(K) = 2.0 * (STC(K) - STC(K+1)) / TEMP1
EFLUX(K) = (DF(K)*DTSDZ(K) - DF(K-1)*DTSDZ(K-1)) - PHI(K)
ELSE IF (K == NSOIL) THEN
DENOM(K) = (ZSNSO(K-1) - ZSNSO(K)) * HCPCT(K)
TEMP1 = ZSNSO(K-1) - ZSNSO(K)
IF(OPT_TBOT == 1) THEN
BOTFLX = 0.
END IF
IF(OPT_TBOT == 2) THEN
DTSDZ(K) = (STC(K) - TBOT) / ( 0.5*(ZSNSO(K-1)+ZSNSO(K)) - ZBOT)
BOTFLX = -DF(K) * DTSDZ(K)
END IF
EFLUX(K) = (-BOTFLX - DF(K-1)*DTSDZ(K-1) ) - PHI(K)
END IF
END DO
DO K = ISNOW+1, NSOIL
IF (K == ISNOW+1) THEN
AI(K) = 0.0
CI(K) = - DF(K) * DDZ(K) / DENOM(K)
IF (OPT_STC == 1) THEN
BI(K) = - CI(K)
END IF
IF (OPT_STC == 2) THEN
BI(K) = - CI(K) + DF(K)/(0.5*ZSNSO(K)*ZSNSO(K)*HCPCT(K))
END IF
ELSE IF (K < NSOIL) THEN
AI(K) = - DF(K-1) * DDZ(K-1) / DENOM(K)
CI(K) = - DF(K ) * DDZ(K ) / DENOM(K)
BI(K) = - (AI(K) + CI (K))
ELSE IF (K == NSOIL) THEN
AI(K) = - DF(K-1) * DDZ(K-1) / DENOM(K)
CI(K) = 0.0
BI(K) = - (AI(K) + CI(K))
END IF
RHSTS(K) = EFLUX(K)/ (-DENOM(K))
END DO
END SUBROUTINE HRT_GLACIER
! ==================================================================================================
! ----------------------------------------------------------------------
SUBROUTINE HSTEP_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT , & !in 1,1
AI ,BI ,CI ,RHSTS , & !inout
STC ) !inout
! ----------------------------------------------------------------------
! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSOIL
INTEGER, INTENT(IN) :: NSNOW
INTEGER, INTENT(IN) :: ISNOW
REAL, INTENT(IN) :: DT
! output & input
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: AI
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: BI
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: CI
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: RHSTS
! local
INTEGER :: K
REAL, DIMENSION(-NSNOW+1:NSOIL) :: RHSTSIN
REAL, DIMENSION(-NSNOW+1:NSOIL) :: CIIN
! ----------------------------------------------------------------------
DO K = ISNOW+1,NSOIL
RHSTS(K) = RHSTS(K) * DT
AI(K) = AI(K) * DT
BI(K) = 1. + BI(K) * DT
CI(K) = CI(K) * DT
END DO
! copy values for input variables before call to rosr12
DO K = ISNOW+1,NSOIL
RHSTSIN(K) = RHSTS(K)
CIIN(K) = CI(K)
END DO
! solve the tri-diagonal matrix equation
CALL ROSR12_GLACIER (CI,AI,BI,CIIN,RHSTSIN,RHSTS,ISNOW+1,NSOIL,NSNOW)
! update snow & soil temperature
DO K = ISNOW+1,NSOIL
STC (K) = STC (K) + CI (K)
END DO
END SUBROUTINE HSTEP_GLACIER
! ==================================================================================================
SUBROUTINE ROSR12_GLACIER (P,A,B,C,D,DELTA,NTOP,NSOIL,NSNOW) 1
! ----------------------------------------------------------------------
! SUBROUTINE ROSR12
! ----------------------------------------------------------------------
! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
! ### ### ### ### ### ###
! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
! # . . # # . # = # . #
! # . . # # . # # . #
! # . . # # . # # . #
! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
! ### ### ### ### ### ###
! ----------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: NTOP
INTEGER, INTENT(IN) :: NSOIL,NSNOW
INTEGER :: K, KK
REAL, DIMENSION(-NSNOW+1:NSOIL),INTENT(IN):: A, B, D
REAL, DIMENSION(-NSNOW+1:NSOIL),INTENT(INOUT):: C,P,DELTA
! ----------------------------------------------------------------------
! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
! ----------------------------------------------------------------------
C (NSOIL) = 0.0
P (NTOP) = - C (NTOP) / B (NTOP)
! ----------------------------------------------------------------------
! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
! ----------------------------------------------------------------------
DELTA (NTOP) = D (NTOP) / B (NTOP)
! ----------------------------------------------------------------------
! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
! ----------------------------------------------------------------------
DO K = NTOP+1,NSOIL
P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
* P (K -1)))
END DO
! ----------------------------------------------------------------------
! SET P TO DELTA FOR LOWEST SOIL LAYER
! ----------------------------------------------------------------------
P (NSOIL) = DELTA (NSOIL)
! ----------------------------------------------------------------------
! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
! ----------------------------------------------------------------------
DO K = NTOP+1,NSOIL
KK = NSOIL - K + (NTOP-1) + 1
P (KK) = P (KK) * P (KK +1) + DELTA (KK)
END DO
! ----------------------------------------------------------------------
END SUBROUTINE ROSR12_GLACIER
! ----------------------------------------------------------------------
! ==================================================================================================
SUBROUTINE PHASECHANGE_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT ,FACT , & !in 1
DZSNSO , & !in
STC ,SNICE ,SNLIQ ,SNEQV ,SNOWH , & !inout
SMC ,SH2O , & !inout
QMELT ,IMELT ,PONDING ) !out
! ----------------------------------------------------------------------
! melting/freezing of snow water and soil water
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! inputs
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers [=3]
INTEGER, INTENT(IN) :: NSOIL !No. of soil layers [=4]
INTEGER, INTENT(IN) :: ISNOW !actual no. of snow layers [<=3]
REAL, INTENT(IN) :: DT !land model time step (sec)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: FACT !temporary
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO !snow/soil layer thickness [m]
! inputs/outputs
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow/soil layer temperature [k]
REAL, DIMENSION(-NSNOW+1:0) , INTENT(INOUT) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1:0) , INTENT(INOUT) :: SNLIQ !snow layer liquid water [mm]
REAL, INTENT(INOUT) :: SNEQV
REAL, INTENT(INOUT) :: SNOWH
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O !soil liquid water [m3/m3]
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SMC !total soil water [m3/m3]
! outputs
REAL, INTENT(OUT) :: QMELT !snowmelt rate [mm/s]
INTEGER, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: IMELT !phase change index
REAL, INTENT(OUT) :: PONDING!snowmelt when snow has no layer [mm]
! local
INTEGER :: J !do loop index
REAL, DIMENSION(-NSNOW+1:NSOIL) :: HM !energy residual [w/m2]
REAL, DIMENSION(-NSNOW+1:NSOIL) :: XM !melting or freezing water [kg/m2]
REAL, DIMENSION(-NSNOW+1:NSOIL) :: WMASS0
REAL, DIMENSION(-NSNOW+1:NSOIL) :: WICE0
REAL, DIMENSION(-NSNOW+1:NSOIL) :: WLIQ0
REAL, DIMENSION(-NSNOW+1:NSOIL) :: MICE !soil/snow ice mass [mm]
REAL, DIMENSION(-NSNOW+1:NSOIL) :: MLIQ !soil/snow liquid water mass [mm]
REAL :: HEATR !energy residual or loss after melting/freezing
REAL :: TEMP1 !temporary variables [kg/m2]
REAL :: PROPOR
REAL :: XMF !total latent heat of phase change
! ----------------------------------------------------------------------
! Initialization
QMELT = 0.
PONDING = 0.
XMF = 0.
DO J = ISNOW+1,0 ! all snow layers
MICE(J) = SNICE(J)
MLIQ(J) = SNLIQ(J)
END DO
DO J = 1, NSOIL ! all soil layers
MLIQ(J) = SH2O(J) * DZSNSO(J) * 1000.
MICE(J) = (SMC(J) - SH2O(J)) * DZSNSO(J) * 1000.
END DO
DO J = ISNOW+1,NSOIL ! all layers
IMELT(J) = 0
HM(J) = 0.
XM(J) = 0.
WICE0(J) = MICE(J)
WLIQ0(J) = MLIQ(J)
WMASS0(J) = MICE(J) + MLIQ(J)
ENDDO
DO J = ISNOW+1,NSOIL
IF (MICE(J) > 0. .AND. STC(J) >= TFRZ) THEN ! melting
IMELT(J) = 1
ENDIF
IF (MLIQ(J) > 0. .AND. STC(J) < TFRZ) THEN ! freezing
IMELT(J) = 2
ENDIF
! If snow exists, but its thickness is not enough to create a layer
IF (ISNOW == 0 .AND. SNEQV > 0. .AND. J == 1) THEN
IF (STC(J) >= TFRZ) THEN
IMELT(J) = 1
ENDIF
ENDIF
ENDDO
! Calculate the energy surplus and loss for melting and freezing
DO J = ISNOW+1,NSOIL
IF (IMELT(J) > 0) THEN
HM(J) = (STC(J)-TFRZ)/FACT(J)
STC(J) = TFRZ
ENDIF
IF (IMELT(J) == 1 .AND. HM(J) < 0.) THEN
HM(J) = 0.
IMELT(J) = 0
ENDIF
IF (IMELT(J) == 2 .AND. HM(J) > 0.) THEN
HM(J) = 0.
IMELT(J) = 0
ENDIF
XM(J) = HM(J)*DT/HFUS
ENDDO
! The rate of melting and freezing for snow without a layer, needs more work.
IF (ISNOW == 0 .AND. SNEQV > 0. .AND. XM(1) > 0.) THEN
TEMP1 = SNEQV
SNEQV = MAX(0.,TEMP1-XM(1))
PROPOR = SNEQV/TEMP1
SNOWH = MAX(0.,PROPOR * SNOWH)
HEATR = HM(1) - HFUS*(TEMP1-SNEQV)/DT
IF (HEATR > 0.) THEN
XM(1) = HEATR*DT/HFUS
HM(1) = HEATR
ELSE
XM(1) = 0.
HM(1) = 0.
ENDIF
QMELT = MAX(0.,(TEMP1-SNEQV))/DT
XMF = HFUS*QMELT
PONDING = TEMP1-SNEQV
ENDIF
! The rate of melting and freezing for snow and soil
DO J = ISNOW+1,NSOIL
IF (IMELT(J) > 0 .AND. ABS(HM(J)) > 0.) THEN
HEATR = 0.
IF (XM(J) > 0.) THEN
MICE(J) = MAX(0., WICE0(J)-XM(J))
HEATR = HM(J) - HFUS*(WICE0(J)-MICE(J))/DT
ELSE IF (XM(J) < 0.) THEN
MICE(J) = MIN(WMASS0(J), WICE0(J)-XM(J))
HEATR = HM(J) - HFUS*(WICE0(J)-MICE(J))/DT
ENDIF
MLIQ(J) = MAX(0.,WMASS0(J)-MICE(J))
IF (ABS(HEATR) > 0.) THEN
STC(J) = STC(J) + FACT(J)*HEATR
IF (J <= 0) THEN ! snow
IF (MLIQ(J)*MICE(J)>0.) STC(J) = TFRZ
END IF
ENDIF
XMF = XMF + HFUS * (WICE0(J)-MICE(J))/DT
IF (J < 1) THEN
QMELT = QMELT + MAX(0.,(WICE0(J)-MICE(J)))/DT
ENDIF
ENDIF
ENDDO
DO J = ISNOW+1,0 ! snow
SNLIQ(J) = MLIQ(J)
SNICE(J) = MICE(J)
END DO
DO J = 1, NSOIL ! soil
SH2O(J) = MLIQ(J) / (1000. * DZSNSO(J))
SMC(J) = (MLIQ(J) + MICE(J)) / (1000. * DZSNSO(J))
END DO
END SUBROUTINE PHASECHANGE_GLACIER
! ==================================================================================================
SUBROUTINE WATER_GLACIER (NSNOW ,NSOIL ,IMELT ,DT ,PRCP ,SFCTMP , & !in 1,1
QVAP ,QDEW ,FICEOLD,ZSOIL , & !in
ISNOW ,SNOWH ,SNEQV ,SNICE ,SNLIQ ,STC , & !inout
DZSNSO ,SH2O ,SICE ,PONDING,ZSNSO , & !inout
RUNSRF ,RUNSUB ,QSNOW ,PONDING1 ,PONDING2,QSNBOT,FPICE ) !out
! ----------------------------------------------------------------------
! Code history:
! Initial code: Guo-Yue Niu, Oct. 2007
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers
INTEGER, INTENT(IN) :: NSOIL !no. of soil layers
INTEGER, DIMENSION(-NSNOW+1:0) , INTENT(IN) :: IMELT !melting state index [1-melt; 2-freeze]
REAL, INTENT(IN) :: DT !main time step (s)
REAL, INTENT(IN) :: PRCP !precipitation (mm/s)
REAL, INTENT(IN) :: SFCTMP !surface air temperature [k]
REAL, INTENT(IN) :: QVAP !soil surface evaporation rate[mm/s]
REAL, INTENT(IN) :: QDEW !soil surface dew rate[mm/s]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: FICEOLD !ice fraction at last timestep
REAL, DIMENSION( 1:NSOIL), INTENT(IN) :: ZSOIL !layer-bottom depth from soil surf (m)
! input/output
INTEGER, INTENT(INOUT) :: ISNOW !actual no. of snow layers
REAL, INTENT(INOUT) :: SNOWH !snow height [m]
REAL, INTENT(INOUT) :: SNEQV !snow water eqv. [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ !snow layer liquid water [mm]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow/soil layer temperature [k]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: DZSNSO !snow/soil layer thickness [m]
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O !soil liquid water content [m3/m3]
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SICE !soil ice content [m3/m3]
REAL , INTENT(INOUT) :: PONDING ![mm]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: ZSNSO !layer-bottom depth from snow surf [m]
! output
REAL, INTENT(OUT) :: RUNSRF !surface runoff [mm/s]
REAL, INTENT(OUT) :: RUNSUB !baseflow (sturation excess) [mm/s]
REAL, INTENT(OUT) :: QSNOW !snow at ground srf (mm/s) [+]
REAL, INTENT(OUT) :: PONDING1
REAL, INTENT(OUT) :: PONDING2
REAL, INTENT(OUT) :: QSNBOT !melting water out of snow bottom [mm/s]
REAL, INTENT(OUT) :: FPICE !precipitation frozen fraction
! local
REAL :: QRAIN !rain at ground srf (mm) [+]
REAL :: QSEVA !soil surface evap rate [mm/s]
REAL :: QSDEW !soil surface dew rate [mm/s]
REAL :: QSNFRO !snow surface frost rate[mm/s]
REAL :: QSNSUB !snow surface sublimation rate [mm/s]
REAL :: SNOWHIN !snow depth increasing rate (m/s)
REAL :: SNOFLOW !glacier flow [mm/s]
REAL :: BDFALL !density of new snow (mm water/m snow)
! ----------------------------------------------------------------------
! initialize
SNOFLOW = 0.
RUNSUB = 0.
RUNSRF = 0.
! --------------------------------------------------------------------
! partition precipitation into rain and snow (from CANWATER)
! Jordan (1991)
IF(OPT_SNF == 1) THEN
IF(SFCTMP > TFRZ+2.5)THEN
FPICE = 0.
ELSE
IF(SFCTMP <= TFRZ+0.5)THEN
FPICE = 1.0
ELSE IF(SFCTMP <= TFRZ+2.)THEN
FPICE = 1.-(-54.632 + 0.2*SFCTMP)
ELSE
FPICE = 0.6
ENDIF
ENDIF
ENDIF
IF(OPT_SNF == 2) THEN
IF(SFCTMP >= TFRZ+2.2) THEN
FPICE = 0.
ELSE
FPICE = 1.0
ENDIF
ENDIF
IF(OPT_SNF == 3) THEN
IF(SFCTMP >= TFRZ) THEN
FPICE = 0.
ELSE
FPICE = 1.0
ENDIF
ENDIF
! print*, 'fpice: ',fpice
! Hedstrom NR and JW Pomeroy (1998), Hydrol. Processes, 12, 1611-1625
! fresh snow density
BDFALL = MAX(120.,67.92+51.25*EXP((SFCTMP-TFRZ)/2.59))
QRAIN = PRCP * (1.-FPICE)
QSNOW = PRCP * FPICE
SNOWHIN = QSNOW/BDFALL
! print *, 'qrain, qsnow',qrain,qsnow,qrain*dt,qsnow*dt
! sublimation, frost, evaporation, and dew
QSNSUB = 0.
IF (SNEQV > 0.) THEN
QSNSUB = MIN(QVAP, SNEQV/DT)
ENDIF
QSEVA = QVAP-QSNSUB
QSNFRO = 0.
IF (SNEQV > 0.) THEN
QSNFRO = QDEW
ENDIF
QSDEW = QDEW - QSNFRO
! print *, 'qvap',qvap,qvap*dt
! print *, 'qsnsub',qsnsub,qsnsub*dt
! print *, 'qseva',qseva,qseva*dt
! print *, 'qsnfro',qsnfro,qsnfro*dt
! print *, 'qdew',qdew,qdew*dt
! print *, 'qsdew',qsdew,qsdew*dt
CALL SNOWWATER_GLACIER (NSNOW ,NSOIL ,IMELT ,DT ,SFCTMP , & !in
SNOWHIN,QSNOW ,QSNFRO ,QSNSUB ,QRAIN , & !in
FICEOLD,ZSOIL , & !in
ISNOW ,SNOWH ,SNEQV ,SNICE ,SNLIQ , & !inout
SH2O ,SICE ,STC ,DZSNSO ,ZSNSO , & !inout
QSNBOT ,SNOFLOW,PONDING1 ,PONDING2) !out
!PONDING: melting water from snow when there is no layer
RUNSRF = (PONDING+PONDING1+PONDING2)/DT
IF(ISNOW == 0) THEN
RUNSRF = RUNSRF+(QSNBOT + QSDEW + QRAIN)
ELSE
RUNSRF = RUNSRF+(QSNBOT + QSDEW)
ENDIF
RUNSUB = SNOFLOW !mm/s allow excess snow to disappear as runsub
END SUBROUTINE WATER_GLACIER
! ==================================================================================================
! ----------------------------------------------------------------------
SUBROUTINE SNOWWATER_GLACIER (NSNOW ,NSOIL ,IMELT ,DT ,SFCTMP , & !in 1,5
SNOWHIN,QSNOW ,QSNFRO ,QSNSUB ,QRAIN , & !in
FICEOLD,ZSOIL , & !in
ISNOW ,SNOWH ,SNEQV ,SNICE ,SNLIQ , & !inout
SH2O ,SICE ,STC ,DZSNSO ,ZSNSO , & !inout
QSNBOT ,SNOFLOW,PONDING1 ,PONDING2) !out
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers
INTEGER, INTENT(IN) :: NSOIL !no. of soil layers
INTEGER, DIMENSION(-NSNOW+1:0) , INTENT(IN) :: IMELT !melting state index [0-no melt;1-melt]
REAL, INTENT(IN) :: DT !time step (s)
REAL, INTENT(IN) :: SFCTMP !surface air temperature [k]
REAL, INTENT(IN) :: SNOWHIN!snow depth increasing rate (m/s)
REAL, INTENT(IN) :: QSNOW !snow at ground srf (mm/s) [+]
REAL, INTENT(IN) :: QSNFRO !snow surface frost rate[mm/s]
REAL, INTENT(IN) :: QSNSUB !snow surface sublimation rate[mm/s]
REAL, INTENT(IN) :: QRAIN !snow surface rain rate[mm/s]
REAL, DIMENSION(-NSNOW+1:0) , INTENT(IN) :: FICEOLD!ice fraction at last timestep
REAL, DIMENSION( 1:NSOIL), INTENT(IN) :: ZSOIL !layer-bottom depth from soil surf (m)
! input & output
INTEGER, INTENT(INOUT) :: ISNOW !actual no. of snow layers
REAL, INTENT(INOUT) :: SNOWH !snow height [m]
REAL, INTENT(INOUT) :: SNEQV !snow water eqv. [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ !snow layer liquid water [mm]
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O !soil liquid moisture (m3/m3)
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SICE !soil ice moisture (m3/m3)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow layer temperature [k]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: DZSNSO !snow/soil layer thickness [m]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: ZSNSO !layer-bottom depth from snow surf [m]
! output
REAL, INTENT(OUT) :: QSNBOT !melting water out of snow bottom [mm/s]
REAL, INTENT(OUT) :: SNOFLOW!glacier flow [mm]
REAL, INTENT(OUT) :: PONDING1
REAL, INTENT(OUT) :: PONDING2
! local
INTEGER :: IZ
REAL :: BDSNOW !bulk density of snow (kg/m3)
! ----------------------------------------------------------------------
SNOFLOW = 0.0
PONDING1 = 0.0
PONDING2 = 0.0
CALL SNOWFALL_GLACIER (NSOIL ,NSNOW ,DT ,QSNOW ,SNOWHIN, & !in
SFCTMP , & !in
ISNOW ,SNOWH ,DZSNSO ,STC ,SNICE , & !inout
SNLIQ ,SNEQV ) !inout
IF(ISNOW < 0) THEN !WHEN MORE THAN ONE LAYER
CALL COMPACT_GLACIER (NSNOW ,NSOIL ,DT ,STC ,SNICE , & !in
SNLIQ ,IMELT ,FICEOLD, & !in
ISNOW ,DZSNSO ) !inout
CALL COMBINE_GLACIER (NSNOW ,NSOIL , & !in
ISNOW ,SH2O ,STC ,SNICE ,SNLIQ , & !inout
DZSNSO ,SICE ,SNOWH ,SNEQV , & !inout
PONDING1 ,PONDING2) !out
CALL DIVIDE_GLACIER (NSNOW ,NSOIL , & !in
ISNOW ,STC ,SNICE ,SNLIQ ,DZSNSO ) !inout
END IF
!SET EMPTY SNOW LAYERS TO ZERO
DO IZ = -NSNOW+1, ISNOW
SNICE(IZ) = 0.
SNLIQ(IZ) = 0.
STC(IZ) = 0.
DZSNSO(IZ)= 0.
ZSNSO(IZ) = 0.
ENDDO
CALL SNOWH2O_GLACIER (NSNOW ,NSOIL ,DT ,QSNFRO ,QSNSUB , & !in
QRAIN , & !in
ISNOW ,DZSNSO ,SNOWH ,SNEQV ,SNICE , & !inout
SNLIQ ,SH2O ,SICE ,STC , & !inout
PONDING1 ,PONDING2 , & !inout
QSNBOT ) !out
!to obtain equilibrium state of snow in glacier region
IF(SNEQV > 2000.) THEN ! 2000 mm -> maximum water depth
BDSNOW = SNICE(0) / DZSNSO(0)
SNOFLOW = (SNEQV - 2000.)
SNICE(0) = SNICE(0) - SNOFLOW
DZSNSO(0) = DZSNSO(0) - SNOFLOW/BDSNOW
SNOFLOW = SNOFLOW / DT
END IF
! sum up snow mass for layered snow
IF(ISNOW /= 0) THEN
SNEQV = 0.
DO IZ = ISNOW+1,0
SNEQV = SNEQV + SNICE(IZ) + SNLIQ(IZ)
ENDDO
END IF
! Reset ZSNSO and layer thinkness DZSNSO
DO IZ = ISNOW+1, 0
DZSNSO(IZ) = -DZSNSO(IZ)
END DO
DZSNSO(1) = ZSOIL(1)
DO IZ = 2,NSOIL
DZSNSO(IZ) = (ZSOIL(IZ) - ZSOIL(IZ-1))
END DO
ZSNSO(ISNOW+1) = DZSNSO(ISNOW+1)
DO IZ = ISNOW+2 ,NSOIL
ZSNSO(IZ) = ZSNSO(IZ-1) + DZSNSO(IZ)
ENDDO
DO IZ = ISNOW+1 ,NSOIL
DZSNSO(IZ) = -DZSNSO(IZ)
END DO
END SUBROUTINE SNOWWATER_GLACIER
! ==================================================================================================
SUBROUTINE SNOWFALL_GLACIER (NSOIL ,NSNOW ,DT ,QSNOW ,SNOWHIN , & !in 1
SFCTMP , & !in
ISNOW ,SNOWH ,DZSNSO ,STC ,SNICE , & !inout
SNLIQ ,SNEQV ) !inout
! ----------------------------------------------------------------------
! snow depth and density to account for the new snowfall.
! new values of snow depth & density returned.
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSOIL !no. of soil layers
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers
REAL, INTENT(IN) :: DT !main time step (s)
REAL, INTENT(IN) :: QSNOW !snow at ground srf (mm/s) [+]
REAL, INTENT(IN) :: SNOWHIN!snow depth increasing rate (m/s)
REAL, INTENT(IN) :: SFCTMP !surface air temperature [k]
! input and output
INTEGER, INTENT(INOUT) :: ISNOW !actual no. of snow layers
REAL, INTENT(INOUT) :: SNOWH !snow depth [m]
REAL, INTENT(INOUT) :: SNEQV !swow water equivalent [m]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: DZSNSO !thickness of snow/soil layers (m)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow layer temperature [k]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ !snow layer liquid water [mm]
! local
INTEGER :: NEWNODE ! 0-no new layers, 1-creating new layers
! ----------------------------------------------------------------------
NEWNODE = 0
! shallow snow / no layer
IF(ISNOW == 0 .and. QSNOW > 0.) THEN
SNOWH = SNOWH + SNOWHIN * DT
SNEQV = SNEQV + QSNOW * DT
END IF
! creating a new layer
IF(ISNOW == 0 .AND. QSNOW>0. .AND. SNOWH >= 0.05) THEN
ISNOW = -1
NEWNODE = 1
DZSNSO(0)= SNOWH
SNOWH = 0.
STC(0) = MIN(273.16, SFCTMP) ! temporary setup
SNICE(0) = SNEQV
SNLIQ(0) = 0.
END IF
! snow with layers
IF(ISNOW < 0 .AND. NEWNODE == 0 .AND. QSNOW > 0.) then
SNICE(ISNOW+1) = SNICE(ISNOW+1) + QSNOW * DT
DZSNSO(ISNOW+1) = DZSNSO(ISNOW+1) + SNOWHIN * DT
ENDIF
! ----------------------------------------------------------------------
END SUBROUTINE SNOWFALL_GLACIER
! ==================================================================================================
! ----------------------------------------------------------------------
SUBROUTINE COMPACT_GLACIER (NSNOW ,NSOIL ,DT ,STC ,SNICE , & !in 1
SNLIQ ,IMELT ,FICEOLD, & !in
ISNOW ,DZSNSO ) !inout
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSOIL !no. of soil layers [ =4]
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers [ =3]
INTEGER, DIMENSION(-NSNOW+1:0) , INTENT(IN) :: IMELT !melting state index [0-no melt;1-melt]
REAL, INTENT(IN) :: DT !time step (sec)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: STC !snow layer temperature [k]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNLIQ !snow layer liquid water [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: FICEOLD!ice fraction at last timestep
! input and output
INTEGER, INTENT(INOUT) :: ISNOW ! actual no. of snow layers
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: DZSNSO ! snow layer thickness [m]
! local
REAL, PARAMETER :: C2 = 21.e-3 ![m3/kg] ! default 21.e-3
REAL, PARAMETER :: C3 = 2.5e-6 ![1/s]
REAL, PARAMETER :: C4 = 0.04 ![1/k]
REAL, PARAMETER :: C5 = 2.0 !
REAL, PARAMETER :: DM = 100.0 !upper Limit on destructive metamorphism compaction [kg/m3]
REAL, PARAMETER :: ETA0 = 0.8e+6 !viscosity coefficient [kg-s/m2]
!according to Anderson, it is between 0.52e6~1.38e6
REAL :: BURDEN !pressure of overlying snow [kg/m2]
REAL :: DDZ1 !rate of settling of snow pack due to destructive metamorphism.
REAL :: DDZ2 !rate of compaction of snow pack due to overburden.
REAL :: DDZ3 !rate of compaction of snow pack due to melt [1/s]
REAL :: DEXPF !EXPF=exp(-c4*(273.15-STC)).
REAL :: TD !STC - TFRZ [K]
REAL :: PDZDTC !nodal rate of change in fractional-thickness due to compaction [fraction/s]
REAL :: VOID !void (1 - SNICE - SNLIQ)
REAL :: WX !water mass (ice + liquid) [kg/m2]
REAL :: BI !partial density of ice [kg/m3]
REAL, DIMENSION(-NSNOW+1:0) :: FICE !fraction of ice at current time step
INTEGER :: J
! ----------------------------------------------------------------------
BURDEN = 0.0
DO J = ISNOW+1, 0
WX = SNICE(J) + SNLIQ(J)
FICE(J) = SNICE(J) / WX
VOID = 1. - (SNICE(J)/DENICE + SNLIQ(J)/DENH2O) / DZSNSO(J)
! Allow compaction only for non-saturated node and higher ice lens node.
IF (VOID > 0.001 .AND. SNICE(J) > 0.1) THEN
BI = SNICE(J) / DZSNSO(J)
TD = MAX(0.,TFRZ-STC(J))
DEXPF = EXP(-C4*TD)
! Settling as a result of destructive metamorphism
DDZ1 = -C3*DEXPF
IF (BI > DM) DDZ1 = DDZ1*EXP(-46.0E-3*(BI-DM))
! Liquid water term
IF (SNLIQ(J) > 0.01*DZSNSO(J)) DDZ1=DDZ1*C5
! Compaction due to overburden
DDZ2 = -(BURDEN+0.5*WX)*EXP(-0.08*TD-C2*BI)/ETA0 ! 0.5*WX -> self-burden
! Compaction occurring during melt
IF (IMELT(J) == 1) THEN
DDZ3 = MAX(0.,(FICEOLD(J) - FICE(J))/MAX(1.E-6,FICEOLD(J)))
DDZ3 = - DDZ3/DT ! sometimes too large
ELSE
DDZ3 = 0.
END IF
! Time rate of fractional change in DZ (units of s-1)
PDZDTC = (DDZ1 + DDZ2 + DDZ3)*DT
PDZDTC = MAX(-0.5,PDZDTC)
! The change in DZ due to compaction
DZSNSO(J) = DZSNSO(J)*(1.+PDZDTC)
END IF
! Pressure of overlying snow
BURDEN = BURDEN + WX
END DO
END SUBROUTINE COMPACT_GLACIER
! ==================================================================================================
SUBROUTINE COMBINE_GLACIER (NSNOW ,NSOIL , & !in 2,1
ISNOW ,SH2O ,STC ,SNICE ,SNLIQ , & !inout
DZSNSO ,SICE ,SNOWH ,SNEQV , & !inout
PONDING1 ,PONDING2) !inout
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers
INTEGER, INTENT(IN) :: NSOIL !no. of soil layers
! input and output
INTEGER, INTENT(INOUT) :: ISNOW !actual no. of snow layers
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O !soil liquid moisture (m3/m3)
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SICE !soil ice moisture (m3/m3)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow layer temperature [k]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ !snow layer liquid water [mm]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: DZSNSO!snow layer depth [m]
REAL, INTENT(INOUT) :: SNEQV !snow water equivalent [m]
REAL, INTENT(INOUT) :: SNOWH !snow depth [m]
REAL, INTENT(INOUT) :: PONDING1
REAL, INTENT(INOUT) :: PONDING2
! local variables:
INTEGER :: I,J,K,L ! node indices
INTEGER :: ISNOW_OLD ! number of top snow layer
INTEGER :: MSSI ! node index
INTEGER :: NEIBOR ! adjacent node selected for combination
REAL :: ZWICE ! total ice mass in snow
REAL :: ZWLIQ ! total liquid water in snow
REAL :: DZMIN(3) ! minimum of top snow layer
DATA DZMIN /0.045, 0.05, 0.2/
!-----------------------------------------------------------------------
ISNOW_OLD = ISNOW
DO J = ISNOW_OLD+1,0
IF (SNICE(J) <= .1) THEN
IF(J /= 0) THEN
SNLIQ(J+1) = SNLIQ(J+1) + SNLIQ(J)
SNICE(J+1) = SNICE(J+1) + SNICE(J)
ELSE
IF (ISNOW_OLD < -1) THEN
SNLIQ(J-1) = SNLIQ(J-1) + SNLIQ(J)
SNICE(J-1) = SNICE(J-1) + SNICE(J)
ELSE
PONDING1 = PONDING1 + SNLIQ(J) ! ISNOW WILL GET SET TO ZERO BELOW
SNEQV = SNICE(J) ! PONDING WILL GET ADDED TO PONDING FROM
SNOWH = DZSNSO(J) ! PHASECHANGE WHICH SHOULD BE ZERO HERE
SNLIQ(J) = 0.0 ! BECAUSE THERE IT WAS ONLY CALCULATED
SNICE(J) = 0.0 ! FOR THIN SNOW
DZSNSO(J) = 0.0
ENDIF
! SH2O(1) = SH2O(1)+SNLIQ(J)/(DZSNSO(1)*1000.)
! SICE(1) = SICE(1)+SNICE(J)/(DZSNSO(1)*1000.)
ENDIF
! shift all elements above this down by one.
IF (J > ISNOW+1 .AND. ISNOW < -1) THEN
DO I = J, ISNOW+2, -1
STC(I) = STC(I-1)
SNLIQ(I) = SNLIQ(I-1)
SNICE(I) = SNICE(I-1)
DZSNSO(I)= DZSNSO(I-1)
END DO
END IF
ISNOW = ISNOW + 1
END IF
END DO
! to conserve water in case of too large surface sublimation
IF(SICE(1) < 0.) THEN
SH2O(1) = SH2O(1) + SICE(1)
SICE(1) = 0.
END IF
SNEQV = 0.
SNOWH = 0.
ZWICE = 0.
ZWLIQ = 0.
DO J = ISNOW+1,0
SNEQV = SNEQV + SNICE(J) + SNLIQ(J)
SNOWH = SNOWH + DZSNSO(J)
ZWICE = ZWICE + SNICE(J)
ZWLIQ = ZWLIQ + SNLIQ(J)
END DO
! check the snow depth - all snow gone
! the liquid water assumes ponding on soil surface.
IF (SNOWH < 0.05 .AND. ISNOW < 0 ) THEN
ISNOW = 0
SNEQV = ZWICE
PONDING2 = PONDING2 + ZWLIQ ! LIMIT OF ISNOW < 0 MEANS INPUT PONDING
IF(SNEQV <= 0.) SNOWH = 0. ! SHOULD BE ZERO; SEE ABOVE
END IF
! IF (SNOWH < 0.05 ) THEN
! ISNOW = 0
! SNEQV = ZWICE
! SH2O(1) = SH2O(1) + ZWLIQ / (DZSNSO(1) * 1000.)
! IF(SNEQV <= 0.) SNOWH = 0.
! END IF
! check the snow depth - snow layers combined
IF (ISNOW < -1) THEN
ISNOW_OLD = ISNOW
MSSI = 1
DO I = ISNOW_OLD+1,0
IF (DZSNSO(I) < DZMIN(MSSI)) THEN
IF (I == ISNOW+1) THEN
NEIBOR = I + 1
ELSE IF (I == 0) THEN
NEIBOR = I - 1
ELSE
NEIBOR = I + 1
IF ((DZSNSO(I-1)+DZSNSO(I)) < (DZSNSO(I+1)+DZSNSO(I))) NEIBOR = I-1
END IF
! Node l and j are combined and stored as node j.
IF (NEIBOR > I) THEN
J = NEIBOR
L = I
ELSE
J = I
L = NEIBOR
END IF
CALL COMBO_GLACIER (DZSNSO(J), SNLIQ(J), SNICE(J), &
STC(J), DZSNSO(L), SNLIQ(L), SNICE(L), STC(L) )
! Now shift all elements above this down one.
IF (J-1 > ISNOW+1) THEN
DO K = J-1, ISNOW+2, -1
STC(K) = STC(K-1)
SNICE(K) = SNICE(K-1)
SNLIQ(K) = SNLIQ(K-1)
DZSNSO(K) = DZSNSO(K-1)
END DO
END IF
! Decrease the number of snow layers
ISNOW = ISNOW + 1
IF (ISNOW >= -1) EXIT
ELSE
! The layer thickness is greater than the prescribed minimum value
MSSI = MSSI + 1
END IF
END DO
END IF
END SUBROUTINE COMBINE_GLACIER
! ==================================================================================================
! ----------------------------------------------------------------------
SUBROUTINE COMBO_GLACIER(DZ, WLIQ, WICE, T, DZ2, WLIQ2, WICE2, T2) 3
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! ----------------------------------------------------------------------s
! input
REAL, INTENT(IN) :: DZ2 !nodal thickness of 2 elements being combined [m]
REAL, INTENT(IN) :: WLIQ2 !liquid water of element 2 [kg/m2]
REAL, INTENT(IN) :: WICE2 !ice of element 2 [kg/m2]
REAL, INTENT(IN) :: T2 !nodal temperature of element 2 [k]
REAL, INTENT(INOUT) :: DZ !nodal thickness of 1 elements being combined [m]
REAL, INTENT(INOUT) :: WLIQ !liquid water of element 1
REAL, INTENT(INOUT) :: WICE !ice of element 1 [kg/m2]
REAL, INTENT(INOUT) :: T !node temperature of element 1 [k]
! local
REAL :: DZC !total thickness of nodes 1 and 2 (DZC=DZ+DZ2).
REAL :: WLIQC !combined liquid water [kg/m2]
REAL :: WICEC !combined ice [kg/m2]
REAL :: TC !combined node temperature [k]
REAL :: H !enthalpy of element 1 [J/m2]
REAL :: H2 !enthalpy of element 2 [J/m2]
REAL :: HC !temporary
!-----------------------------------------------------------------------
DZC = DZ+DZ2
WICEC = (WICE+WICE2)
WLIQC = (WLIQ+WLIQ2)
H = (CICE*WICE+CWAT*WLIQ) * (T-TFRZ)+HFUS*WLIQ
H2= (CICE*WICE2+CWAT*WLIQ2) * (T2-TFRZ)+HFUS*WLIQ2
HC = H + H2
IF(HC < 0.)THEN
TC = TFRZ + HC/(CICE*WICEC + CWAT*WLIQC)
ELSE IF (HC.LE.HFUS*WLIQC) THEN
TC = TFRZ
ELSE
TC = TFRZ + (HC - HFUS*WLIQC) / (CICE*WICEC + CWAT*WLIQC)
END IF
DZ = DZC
WICE = WICEC
WLIQ = WLIQC
T = TC
END SUBROUTINE COMBO_GLACIER
! ==================================================================================================
SUBROUTINE DIVIDE_GLACIER (NSNOW ,NSOIL , & !in 1,2
ISNOW ,STC ,SNICE ,SNLIQ ,DZSNSO ) !inout
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers [ =3]
INTEGER, INTENT(IN) :: NSOIL !no. of soil layers [ =4]
! input and output
INTEGER , INTENT(INOUT) :: ISNOW !actual no. of snow layers
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow layer temperature [k]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ !snow layer liquid water [mm]
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: DZSNSO!snow layer depth [m]
! local variables:
INTEGER :: J !indices
INTEGER :: MSNO !number of layer (top) to MSNO (bot)
REAL :: DRR !thickness of the combined [m]
REAL, DIMENSION( 1:NSNOW) :: DZ !snow layer thickness [m]
REAL, DIMENSION( 1:NSNOW) :: SWICE !partial volume of ice [m3/m3]
REAL, DIMENSION( 1:NSNOW) :: SWLIQ !partial volume of liquid water [m3/m3]
REAL, DIMENSION( 1:NSNOW) :: TSNO !node temperature [k]
REAL :: ZWICE !temporary
REAL :: ZWLIQ !temporary
REAL :: PROPOR!temporary
REAL :: DTDZ !temporary
! ----------------------------------------------------------------------
DO J = 1,NSNOW
IF (J <= ABS(ISNOW)) THEN
DZ(J) = DZSNSO(J+ISNOW)
SWICE(J) = SNICE(J+ISNOW)
SWLIQ(J) = SNLIQ(J+ISNOW)
TSNO(J) = STC(J+ISNOW)
END IF
END DO
MSNO = ABS(ISNOW)
IF (MSNO == 1) THEN
! Specify a new snow layer
IF (DZ(1) > 0.05) THEN
MSNO = 2
DZ(1) = DZ(1)/2.
SWICE(1) = SWICE(1)/2.
SWLIQ(1) = SWLIQ(1)/2.
DZ(2) = DZ(1)
SWICE(2) = SWICE(1)
SWLIQ(2) = SWLIQ(1)
TSNO(2) = TSNO(1)
END IF
END IF
IF (MSNO > 1) THEN
IF (DZ(1) > 0.05) THEN
DRR = DZ(1) - 0.05
PROPOR = DRR/DZ(1)
ZWICE = PROPOR*SWICE(1)
ZWLIQ = PROPOR*SWLIQ(1)
PROPOR = 0.05/DZ(1)
SWICE(1) = PROPOR*SWICE(1)
SWLIQ(1) = PROPOR*SWLIQ(1)
DZ(1) = 0.05
CALL COMBO_GLACIER (DZ(2), SWLIQ(2), SWICE(2), TSNO(2), DRR, &
ZWLIQ, ZWICE, TSNO(1))
! subdivide a new layer
IF (MSNO <= 2 .AND. DZ(2) > 0.10) THEN
MSNO = 3
DTDZ = (TSNO(1) - TSNO(2))/((DZ(1)+DZ(2))/2.)
DZ(2) = DZ(2)/2.
SWICE(2) = SWICE(2)/2.
SWLIQ(2) = SWLIQ(2)/2.
DZ(3) = DZ(2)
SWICE(3) = SWICE(2)
SWLIQ(3) = SWLIQ(2)
TSNO(3) = TSNO(2) - DTDZ*DZ(2)/2.
IF (TSNO(3) >= TFRZ) THEN
TSNO(3) = TSNO(2)
ELSE
TSNO(2) = TSNO(2) + DTDZ*DZ(2)/2.
ENDIF
END IF
END IF
END IF
IF (MSNO > 2) THEN
IF (DZ(2) > 0.2) THEN
DRR = DZ(2) - 0.2
PROPOR = DRR/DZ(2)
ZWICE = PROPOR*SWICE(2)
ZWLIQ = PROPOR*SWLIQ(2)
PROPOR = 0.2/DZ(2)
SWICE(2) = PROPOR*SWICE(2)
SWLIQ(2) = PROPOR*SWLIQ(2)
DZ(2) = 0.2
CALL COMBO_GLACIER (DZ(3), SWLIQ(3), SWICE(3), TSNO(3), DRR, &
ZWLIQ, ZWICE, TSNO(2))
END IF
END IF
ISNOW = -MSNO
DO J = ISNOW+1,0
DZSNSO(J) = DZ(J-ISNOW)
SNICE(J) = SWICE(J-ISNOW)
SNLIQ(J) = SWLIQ(J-ISNOW)
STC(J) = TSNO(J-ISNOW)
END DO
! DO J = ISNOW+1,NSOIL
! WRITE(*,'(I5,7F10.3)') J, DZSNSO(J), SNICE(J), SNLIQ(J),STC(J)
! END DO
END SUBROUTINE DIVIDE_GLACIER
! ==================================================================================================
SUBROUTINE SNOWH2O_GLACIER (NSNOW ,NSOIL ,DT ,QSNFRO ,QSNSUB , & !in 1,1
QRAIN , & !in
ISNOW ,DZSNSO ,SNOWH ,SNEQV ,SNICE , & !inout
SNLIQ ,SH2O ,SICE ,STC , & !inout
PONDING1 ,PONDING2 , & !inout
QSNBOT ) !out
! ----------------------------------------------------------------------
! Renew the mass of ice lens (SNICE) and liquid (SNLIQ) of the
! surface snow layer resulting from sublimation (frost) / evaporation (dew)
! ----------------------------------------------------------------------
IMPLICIT NONE
! ----------------------------------------------------------------------
! input
INTEGER, INTENT(IN) :: NSNOW !maximum no. of snow layers[=3]
INTEGER, INTENT(IN) :: NSOIL !No. of soil layers[=4]
REAL, INTENT(IN) :: DT !time step
REAL, INTENT(IN) :: QSNFRO !snow surface frost rate[mm/s]
REAL, INTENT(IN) :: QSNSUB !snow surface sublimation rate[mm/s]
REAL, INTENT(IN) :: QRAIN !snow surface rain rate[mm/s]
! output
REAL, INTENT(OUT) :: QSNBOT !melting water out of snow bottom [mm/s]
! input and output
INTEGER, INTENT(INOUT) :: ISNOW !actual no. of snow layers
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: DZSNSO ! snow layer depth [m]
REAL, INTENT(INOUT) :: SNOWH !snow height [m]
REAL, INTENT(INOUT) :: SNEQV !snow water eqv. [mm]
REAL, DIMENSION(-NSNOW+1:0), INTENT(INOUT) :: SNICE !snow layer ice [mm]
REAL, DIMENSION(-NSNOW+1:0), INTENT(INOUT) :: SNLIQ !snow layer liquid water [mm]
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O !soil liquid moisture (m3/m3)
REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SICE !soil ice moisture (m3/m3)
REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC !snow layer temperature [k]
REAL, INTENT(INOUT) :: PONDING1
REAL, INTENT(INOUT) :: PONDING2
! local variables:
INTEGER :: J !do loop/array indices
REAL :: QIN !water flow into the element (mm/s)
REAL :: QOUT !water flow out of the element (mm/s)
REAL :: WGDIF !ice mass after minus sublimation
REAL, DIMENSION(-NSNOW+1:0) :: VOL_LIQ !partial volume of liquid water in layer
REAL, DIMENSION(-NSNOW+1:0) :: VOL_ICE !partial volume of ice lens in layer
REAL, DIMENSION(-NSNOW+1:0) :: EPORE !effective porosity = porosity - VOL_ICE
REAL :: PROPOR, TEMP
! ----------------------------------------------------------------------
!for the case when SNEQV becomes '0' after 'COMBINE'
IF(SNEQV == 0.) THEN
SH2O(1) = SH2O(1) + (QSNFRO-QSNSUB)*DT/(DZSNSO(1)*1000.)
END IF
! for shallow snow without a layer
! snow surface sublimation may be larger than existing snow mass. To conserve water,
! excessive sublimation is used to reduce soil water. Smaller time steps would tend
! to aviod this problem.
IF(ISNOW == 0 .and. SNEQV > 0.) THEN
TEMP = SNEQV
SNEQV = SNEQV - QSNSUB*DT + QSNFRO*DT
PROPOR = SNEQV/TEMP
SNOWH = MAX(0.,PROPOR * SNOWH)
IF(SNEQV < 0.) THEN
SICE(1) = SICE(1) + SNEQV/(DZSNSO(1)*1000.)
SNEQV = 0.
END IF
IF(SICE(1) < 0.) THEN
SH2O(1) = SH2O(1) + SICE(1)
SICE(1) = 0.
END IF
END IF
IF(SNOWH <= 1.E-8) SNOWH = 0.0
IF(SNEQV <= 1.E-6) SNEQV = 0.0
! for deep snow
IF ( ISNOW < 0 ) THEN !KWM added this IF statement to prevent out-of-bounds array references
WGDIF = SNICE(ISNOW+1) - QSNSUB*DT + QSNFRO*DT
SNICE(ISNOW+1) = WGDIF
IF (WGDIF < 1.e-6 .and. ISNOW <0) THEN
CALL COMBINE_GLACIER (NSNOW ,NSOIL , & !in
ISNOW ,SH2O ,STC ,SNICE ,SNLIQ , & !inout
DZSNSO ,SICE ,SNOWH ,SNEQV , & !inout
PONDING1, PONDING2 ) !inout
ENDIF
!KWM: Subroutine COMBINE can change ISNOW to make it 0 again?
IF ( ISNOW < 0 ) THEN !KWM added this IF statement to prevent out-of-bounds array references
SNLIQ(ISNOW+1) = SNLIQ(ISNOW+1) + QRAIN * DT
SNLIQ(ISNOW+1) = MAX(0., SNLIQ(ISNOW+1))
ENDIF
ENDIF !KWM -- Can the ENDIF be moved toward the end of the subroutine (Just set QSNBOT=0)?
! Porosity and partial volume
!KWM Looks to me like loop index / IF test can be simplified.
DO J = -NSNOW+1, 0
IF (J >= ISNOW+1) THEN
VOL_ICE(J) = MIN(1., SNICE(J)/(DZSNSO(J)*DENICE))
EPORE(J) = 1. - VOL_ICE(J)
VOL_LIQ(J) = MIN(EPORE(J),SNLIQ(J)/(DZSNSO(J)*DENH2O))
END IF
END DO
QIN = 0.
QOUT = 0.
!KWM Looks to me like loop index / IF test can be simplified.
DO J = -NSNOW+1, 0
IF (J >= ISNOW+1) THEN
SNLIQ(J) = SNLIQ(J) + QIN
IF (J <= -1) THEN
IF (EPORE(J) < 0.05 .OR. EPORE(J+1) < 0.05) THEN
QOUT = 0.
ELSE
QOUT = MAX(0.,(VOL_LIQ(J)-SSI*EPORE(J))*DZSNSO(J))
QOUT = MIN(QOUT,(1.-VOL_ICE(J+1)-VOL_LIQ(J+1))*DZSNSO(J+1))
END IF
ELSE
QOUT = MAX(0.,(VOL_LIQ(J) - SSI*EPORE(J))*DZSNSO(J))
END IF
QOUT = QOUT*1000.
SNLIQ(J) = SNLIQ(J) - QOUT
QIN = QOUT
END IF
END DO
! Liquid water from snow bottom to soil
QSNBOT = QOUT / DT ! mm/s
END SUBROUTINE SNOWH2O_GLACIER
! ********************* end of water subroutines ******************************************
! ==================================================================================================
SUBROUTINE ERROR_GLACIER (ILOC ,JLOC ,SWDOWN ,FSA ,FSR ,FIRA , & 1,11
FSH ,FGEV ,SSOIL ,SAG ,PRCP ,EDIR , &
RUNSRF ,RUNSUB ,SNEQV ,DT ,BEG_WB )
! --------------------------------------------------------------------------------------------------
! check surface energy balance and water balance
! --------------------------------------------------------------------------------------------------
IMPLICIT NONE
! --------------------------------------------------------------------------------------------------
! inputs
INTEGER , INTENT(IN) :: ILOC !grid index
INTEGER , INTENT(IN) :: JLOC !grid index
REAL , INTENT(IN) :: SWDOWN !downward solar filtered by sun angle [w/m2]
REAL , INTENT(IN) :: FSA !total absorbed solar radiation (w/m2)
REAL , INTENT(IN) :: FSR !total reflected solar radiation (w/m2)
REAL , INTENT(IN) :: FIRA !total net longwave rad (w/m2) [+ to atm]
REAL , INTENT(IN) :: FSH !total sensible heat (w/m2) [+ to atm]
REAL , INTENT(IN) :: FGEV !ground evaporation heat (w/m2) [+ to atm]
REAL , INTENT(IN) :: SSOIL !ground heat flux (w/m2) [+ to soil]
REAL , INTENT(IN) :: SAG
REAL , INTENT(IN) :: PRCP !precipitation rate (kg m-2 s-1)
REAL , INTENT(IN) :: EDIR !soil surface evaporation rate[mm/s]
REAL , INTENT(IN) :: RUNSRF !surface runoff [mm/s]
REAL , INTENT(IN) :: RUNSUB !baseflow (saturation excess) [mm/s]
REAL , INTENT(IN) :: SNEQV !snow water eqv. [mm]
REAL , INTENT(IN) :: DT !time step [sec]
REAL , INTENT(IN) :: BEG_WB !water storage at begin of a timesetp [mm]
REAL :: END_WB !water storage at end of a timestep [mm]
REAL :: ERRWAT !error in water balance [mm/timestep]
REAL :: ERRENG !error in surface energy balance [w/m2]
REAL :: ERRSW !error in shortwave radiation balance [w/m2]
CHARACTER(len=256) :: message
! --------------------------------------------------------------------------------------------------
ERRSW = SWDOWN - (FSA + FSR)
IF (ERRSW > 0.01) THEN ! w/m2
WRITE(*,*) "SAG =",SAG
WRITE(*,*) "FSA =",FSA
WRITE(*,*) "FSR =",FSR
WRITE(message,*) 'ERRSW =',ERRSW
call wrf_message(trim(message))
call wrf_error_fatal("Radiation budget problem in NOAHMP GLACIER")
END IF
ERRENG = SAG-(FIRA+FSH+FGEV+SSOIL)
IF(ERRENG > 0.01) THEN
write(message,*) 'ERRENG =',ERRENG
call wrf_message(trim(message))
WRITE(message,'(i6,1x,i6,1x,5F10.4)')ILOC,JLOC,SAG,FIRA,FSH,FGEV,SSOIL
call wrf_message(trim(message))
call wrf_error_fatal("Energy budget problem in NOAHMP GLACIER")
END IF
END_WB = SNEQV
ERRWAT = END_WB-BEG_WB-(PRCP-EDIR-RUNSRF-RUNSUB)*DT
IF(ABS(ERRWAT) > 0.1) THEN
if (ERRWAT > 0) then
call wrf_message ('The model is gaining water (ERRWAT is positive)')
else
call wrf_message('The model is losing water (ERRWAT is negative)')
endif
write(message, *) 'ERRWAT =',ERRWAT, "kg m{-2} timestep{-1}"
call wrf_message(trim(message))
WRITE(message,'(" I J END_WB BEG_WB PRCP EDIR RUNSRF RUNSUB")')
call wrf_message(trim(message))
WRITE(message,'(i6,1x,i6,1x,2f15.3,4f11.5)')ILOC,JLOC,END_WB,BEG_WB,PRCP*DT,&
EDIR*DT,RUNSRF*DT,RUNSUB*DT
call wrf_message(trim(message))
call wrf_error_fatal("Water budget problem in NOAHMP GLACIER")
END IF
END SUBROUTINE ERROR_GLACIER
! ==================================================================================================
SUBROUTINE NOAHMP_OPTIONS_GLACIER(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz , & ,2
iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc )
IMPLICIT NONE
INTEGER, INTENT(IN) :: idveg !dynamic vegetation (1 -> off ; 2 -> on) with opt_crs = 1
INTEGER, INTENT(IN) :: iopt_crs !canopy stomatal resistance (1-> Ball-Berry; 2->Jarvis)
INTEGER, INTENT(IN) :: iopt_btr !soil moisture factor for stomatal resistance (1-> Noah; 2-> CLM; 3-> SSiB)
INTEGER, INTENT(IN) :: iopt_run !runoff and groundwater (1->SIMGM; 2->SIMTOP; 3->Schaake96; 4->BATS)
INTEGER, INTENT(IN) :: iopt_sfc !surface layer drag coeff (CH & CM) (1->M-O; 2->Chen97)
INTEGER, INTENT(IN) :: iopt_frz !supercooled liquid water (1-> NY06; 2->Koren99)
INTEGER, INTENT(IN) :: iopt_inf !frozen soil permeability (1-> NY06; 2->Koren99)
INTEGER, INTENT(IN) :: iopt_rad !radiation transfer (1->gap=F(3D,cosz); 2->gap=0; 3->gap=1-Fveg)
INTEGER, INTENT(IN) :: iopt_alb !snow surface albedo (1->BATS; 2->CLASS)
INTEGER, INTENT(IN) :: iopt_snf !rainfall & snowfall (1-Jordan91; 2->BATS; 3->Noah)
INTEGER, INTENT(IN) :: iopt_tbot !lower boundary of soil temperature (1->zero-flux; 2->Noah)
INTEGER, INTENT(IN) :: iopt_stc !snow/soil temperature time scheme (only layer 1)
! 1 -> semi-implicit; 2 -> full implicit (original Noah)
! -------------------------------------------------------------------------------------------------
dveg = idveg
opt_crs = iopt_crs
opt_btr = iopt_btr
opt_run = iopt_run
opt_sfc = iopt_sfc
opt_frz = iopt_frz
opt_inf = iopt_inf
opt_rad = iopt_rad
opt_alb = iopt_alb
opt_snf = iopt_snf
opt_tbot = iopt_tbot
opt_stc = iopt_stc
end subroutine noahmp_options_glacier
END MODULE NOAHMP_GLACIER_ROUTINES
! ==================================================================================================
MODULE MODULE_SF_NOAHMP_GLACIER 1
USE NOAHMP_GLACIER_ROUTINES
USE NOAHMP_GLACIER_GLOBALS
END MODULE MODULE_SF_NOAHMP_GLACIER