CLASSIC
Canadian Land Surface Scheme including Biogeochemical Cycles
Functions/Subroutines
energyBudgetPrep.f90 File Reference

Initializes subarea variables and calculate various parameters for surface energy budget calculations. More...

Functions/Subroutines

subroutine energybudgetprep (THLIQC, THLIQG, THICEC, THICEG, TBARC, TBARG, TBARCS, TBARGS, HCPC, HCPG, TCTOPC, TCBOTC, TCTOPG, TCBOTG, HCPSCS, HCPSGS, TCSNOW, TSNOCS, TSNOGS, WSNOCS, WSNOGS, RHOSCS, RHOSGS, TCANO, TCANS, CEVAP, IEVAP, TBAR1P, WTABLE, ZERO, EVAPC, EVAPCG, EVAPG, EVAPCS, EVPCSG, EVAPGS, GSNOWC, GSNOWG, GZEROC, GZEROG, GZROCS, GZROGS, QMELTC, QMELTG, EVAP, GSNOW, TPONDC, TPONDG, TPNDCS, TPNDGS, QSENSC, QSENSG, QEVAPC, QEVAPG, TACCO, QACCO, TACCS, QACCS, ILMOX, UEX, HBLX, ILMO, UE, HBL, ST, SU, SV, SQ, SRH, CDH, CDM, QSENS, QEVAP, QLWAVG, FSGV, FSGS, FSGG, FLGV, FLGS, FLGG, HFSC, HFSS, HFSG, HEVC, HEVS, HEVG, HMFC, HMFN, QFCF, QFCL, EVPPOT, ACOND, DRAG, THLIQ, THICE, TBAR, ZPOND, TPOND, THPOR, THLMIN, THLRET, THFC, HCPS, TCS, TA, RHOSNO, TSNOW, ZSNOW, WSNOW, TCAN, FC, FCS, DELZ, DELZW, ZBOTW, ISAND, ILG, IL1, IL2, JL, IG, FVEG, TCSATU, TCSATF, FTEMP, FTEMPX, FVAP, FVAPX, RIB, RIBX)
 

Detailed Description

Initializes subarea variables and calculate various parameters for surface energy budget calculations.

Author
D. Verseghy, M. Lazare, A. Wu, Y. Delage, J. P. Paquin, R. Harvey, J. Melton, G. Meyer, E. Humphreys

Function/Subroutine Documentation

◆ energybudgetprep()

subroutine energybudgetprep ( real, dimension(ilg,ig), intent(inout)  THLIQC,
real, dimension(ilg,ig), intent(inout)  THLIQG,
real, dimension(ilg,ig), intent(inout)  THICEC,
real, dimension(ilg,ig), intent(inout)  THICEG,
real, dimension (ilg,ig), intent(out)  TBARC,
real, dimension (ilg,ig), intent(out)  TBARG,
real, dimension(ilg,ig), intent(out)  TBARCS,
real, dimension(ilg,ig), intent(out)  TBARGS,
real, dimension (ilg,ig), intent(out)  HCPC,
real, dimension (ilg,ig), intent(inout)  HCPG,
real, dimension(ilg,ig), intent(inout)  TCTOPC,
real, dimension(ilg,ig), intent(out)  TCBOTC,
real, dimension(ilg,ig), intent(out)  TCTOPG,
real, dimension(ilg,ig), intent(out)  TCBOTG,
real, dimension(ilg), intent(inout)  HCPSCS,
real, dimension(ilg), intent(out)  HCPSGS,
real, dimension(ilg), intent(out)  TCSNOW,
real, dimension(ilg), intent(out)  TSNOCS,
real, dimension(ilg), intent(out)  TSNOGS,
real, dimension(ilg), intent(out)  WSNOCS,
real, dimension(ilg), intent(out)  WSNOGS,
real, dimension(ilg), intent(out)  RHOSCS,
real, dimension(ilg), intent(out)  RHOSGS,
real, dimension (ilg), intent(out)  TCANO,
real, dimension (ilg), intent(out)  TCANS,
real, dimension (ilg), intent(out)  CEVAP,
integer, dimension (ilg), intent(out)  IEVAP,
real, dimension(ilg), intent(out)  TBAR1P,
real, dimension(ilg), intent(out)  WTABLE,
real, dimension (ilg), intent(out)  ZERO,
real, dimension (ilg), intent(out)  EVAPC,
real, dimension(ilg), intent(out)  EVAPCG,
real, dimension (ilg), intent(out)  EVAPG,
real, dimension(ilg), intent(out)  EVAPCS,
real, dimension(ilg), intent(out)  EVPCSG,
real, dimension(ilg), intent(out)  EVAPGS,
real, dimension(ilg), intent(out)  GSNOWC,
real, dimension(ilg), intent(out)  GSNOWG,
real, dimension(ilg), intent(out)  GZEROC,
real, dimension(ilg), intent(out)  GZEROG,
real, dimension(ilg), intent(out)  GZROCS,
real, dimension(ilg), intent(out)  GZROGS,
real, dimension(ilg), intent(out)  QMELTC,
real, dimension(ilg), intent(out)  QMELTG,
real, dimension (ilg), intent(out)  EVAP,
real, dimension (ilg), intent(out)  GSNOW,
real, dimension(ilg), intent(out)  TPONDC,
real, dimension(ilg), intent(out)  TPONDG,
real, dimension(ilg), intent(out)  TPNDCS,
real, dimension(ilg), intent(out)  TPNDGS,
real, dimension(ilg), intent(out)  QSENSC,
real, dimension(ilg), intent(out)  QSENSG,
real, dimension(ilg), intent(out)  QEVAPC,
real, dimension(ilg), intent(out)  QEVAPG,
real, dimension (ilg), intent(out)  TACCO,
real, dimension (ilg), intent(out)  QACCO,
real, dimension (ilg), intent(out)  TACCS,
real, dimension (ilg), intent(out)  QACCS,
real, dimension (ilg), intent(out)  ILMOX,
real, dimension (ilg), intent(out)  UEX,
real, dimension (ilg), intent(out)  HBLX,
real, dimension (ilg), intent(out)  ILMO,
real, dimension (ilg), intent(out)  UE,
real, dimension (ilg), intent(out)  HBL,
real, dimension (ilg), intent(out)  ST,
real, dimension (ilg), intent(out)  SU,
real, dimension (ilg), intent(out)  SV,
real, dimension (ilg), intent(out)  SQ,
real, dimension (ilg), intent(out)  SRH,
real, dimension (ilg), intent(out)  CDH,
real, dimension (ilg), intent(out)  CDM,
real, dimension (ilg), intent(out)  QSENS,
real, dimension (ilg), intent(out)  QEVAP,
real, dimension(ilg), intent(out)  QLWAVG,
real, dimension (ilg), intent(out)  FSGV,
real, dimension (ilg), intent(out)  FSGS,
real, dimension (ilg), intent(out)  FSGG,
real, dimension (ilg), intent(out)  FLGV,
real, dimension (ilg), intent(out)  FLGS,
real, dimension (ilg), intent(out)  FLGG,
real, dimension (ilg), intent(out)  HFSC,
real, dimension (ilg), intent(out)  HFSS,
real, dimension (ilg), intent(out)  HFSG,
real, dimension (ilg), intent(out)  HEVC,
real, dimension (ilg), intent(out)  HEVS,
real, dimension (ilg), intent(out)  HEVG,
real, dimension (ilg), intent(out)  HMFC,
real, dimension (ilg), intent(out)  HMFN,
real, dimension(ilg), intent(out)  QFCF,
real, dimension(ilg), intent(out)  QFCL,
real, dimension(ilg), intent(out)  EVPPOT,
real, dimension (ilg), intent(out)  ACOND,
real, dimension (ilg), intent(out)  DRAG,
real, dimension (ilg,ig), intent(in)  THLIQ,
real, dimension (ilg,ig), intent(in)  THICE,
real, dimension (ilg,ig), intent(in)  TBAR,
real, dimension (ilg), intent(in)  ZPOND,
real, dimension (ilg), intent(in)  TPOND,
real, dimension(ilg,ig), intent(in)  THPOR,
real, dimension(ilg,ig), intent(in)  THLMIN,
real, dimension(ilg,ig), intent(in)  THLRET,
real, dimension (ilg,ig), intent(in)  THFC,
real, dimension (ilg,ig), intent(in)  HCPS,
real, dimension (ilg,ig), intent(in)  TCS,
real, dimension (ilg), intent(in)  TA,
real, dimension(ilg), intent(in)  RHOSNO,
real, dimension (ilg), intent(in)  TSNOW,
real, dimension (ilg), intent(in)  ZSNOW,
real, dimension (ilg), intent(in)  WSNOW,
real, dimension (ilg), intent(in)  TCAN,
real, dimension (ilg), intent(in)  FC,
real, dimension (ilg), intent(in)  FCS,
real, dimension(ig), intent(in)  DELZ,
real, dimension(ilg,ig), intent(in)  DELZW,
real, dimension(ilg,ig), intent(in)  ZBOTW,
integer, dimension (ilg,ig), intent(in)  ISAND,
integer, intent(in)  ILG,
integer, intent(in)  IL1,
integer, intent(in)  IL2,
integer, intent(in)  JL,
integer, intent(in)  IG,
real, dimension(ilg), intent(inout)  FVEG,
real, dimension(ilg), intent(inout)  TCSATU,
real, dimension(ilg), intent(inout)  TCSATF,
real, dimension(ilg), intent(out)  FTEMP,
real, dimension(ilg), intent(out)  FTEMPX,
real, dimension(ilg), intent(out)  FVAP,
real, dimension(ilg), intent(out)  FVAPX,
real, dimension(ilg), intent(out)  RIB,
real, dimension(ilg), intent(out)  RIBX 
)
Parameters
[out]tbarcSubarea temperatures of soil layers [C]
[out]tbargSubarea temperatures of soil layers [C]
[out]tbarcsSubarea temperatures of soil layers [C]
[out]tbargsSubarea temperatures of soil layers [C]
[in,out]thliqcLiquid water content of soil layers under vegetation \([m^3 m^{-3}]\)
[in,out]thliqgLiquid water content of soil layers in bare areas \([m^3 m^{-3}]\)
[in,out]thicecFrozen water content of soil layers under vegetation \([m^3 m^{-3}]\)
[in,out]thicegFrozen water content of soil layers in bare areas \([m^3 m^{-3}]\)
[out]hcpcHeat capacity of soil layers under vegetation \([J m^{-3} K^{-1}] (C_g)\)
[in,out]hcpgHeat capacity of soil layers in bare areas \([J m^{-3} K^{1}] (Cg)\)
[in,out]tctopcThermal conductivity of soil at top of layer in canopy-covered subareas \([W m^{-1} K^{-1}] (\lambda)\)
[out]tcbotcThermal conductivity of soil at bottom of layer in canopy-covered subareas \([W m^{-1} K^{-1}] (\lambda)\)
[out]tctopgThermal conductivity of soil at top of layer in bare ground subareas \([W m^{-1} K^{-1}] (\lambda)\)
[out]tcbotgThermal conductivity of soil at bottom of layer in bare ground subareas \([W m^{-1} K^{-1}] (\lambda)\)
[in,out]hcpscsHeat capacity of snow pack under vegetation canopy \([J m^{-3} K^1] (C_s)\)
[out]hcpsgsHeat capacity of snow pack in bare areas \([J m^{-3} K^{1}] (C_s)\)
[out]tcsnowThermal conductivity of snow \([W m^{-1} K^{-1}]\)
[out]tsnogsTemperature of snow pack in bare areas [K]
[out]tsnocsTemperature of snow pack under vegetation canopy [K]
[out]wsnocsLiquid water content of snow pack under vegetation \([kg m^{-2}]\)
[out]wsnogsLiquid water content of snow pack in bare areas \([kg m^{-2}]\)
[out]rhoscsDensity of snow pack under vegetation canopy \([kg m^{-3}]\)
[out]rhosgsDensity of snow pack in bare areas \([kg m^{-3}]\)
[out]tcanoTemperature of canopy over ground [K]
[out]tcansTemperature of canopy over snow [K]
[out]cevapSoil evaporation efficiency coefficient \([ ] (\beta)\)
[out]tbar1pLumped temperature of ponded water and first soil layer [K]
[out]wtableDepth of water table in soil \([m] (z_{wt})\)
[out]zeroDummy vector containing all zeros
[out]tpondcSubarea temperature of surface ponded water [C]
[out]tpondgSubarea temperature of surface ponded water [C]
[out]tpndcsSubarea temperature of surface ponded water [C]
[out]tpndgsSubarea temperature of surface ponded water [C]
[out]ievapFlag indicating whether soil evaporation is occurring or not
[out]evapcEvaporation from vegetation over ground \([m s^{-1}]\)
[out]evapcgEvaporation from ground under vegetation \([m s^{-1}]\)
[out]evapgEvaporation from bare ground \([m s^{-1}]\)
[out]evapcsEvaporation from vegetation over snow \([m s^{-1}]\)
[out]evpcsgEvaporation from snow under vegetation \([m s^{-1}]\)
[out]evapgsEvaporation from snow on bare ground \([m s^{-1}]\)
[out]gsnowcHeat flux at top of snow pack under canopy \([W m^{-2}]\)
[out]gsnowgHeat flux at top of snow pack over bare ground \([W m^{-2}]\)
[out]gzerocSubarea heat flux at soil surface \([W m^{-2}]\)
[out]gzerogSubarea heat flux at soil surface \([W m^{-2}]\)
[out]gzrocsSubarea heat flux at soil surface \([W m^{-2}]\)
[out]gzrogsSubarea heat flux at soil surface \([W m^{-2}]\)
[out]qmeltcHeat to be used for melting snow under canopy \([W m^{-2}]\)
[out]qmeltgHeat to be used for melting snow on bare ground \([W m^{-2}]\)
[out]evapDiagnosed total surface water vapour flux over modelled area \([kg m^{-2} s^{-1}]\)
[out]qsenscSensible heat flux from vegetation canopy over subarea \([W m^{-2}]\)
[out]qsensgSensible heat flux from ground over subarea \([W m^{-2}]\)
[out]qevapcLatent heat flux from vegetation canopy over subarea \([W m^{-2}]\)
[out]qevapgLatent heat flux from ground over subarea \([W m^{-2}]\)
[out]taccoTemperature of air within vegetation canopy space over bare ground [K]
[out]qaccoSpecific humidity of air within vegetation canopy space over bare ground \([kg kg^{-1}]\)
[out]taccsTemperature of air within vegetation canopy space over snow [K]
[out]qaccsSpecific humidity of air within vegetation canopy space over snow \([kg kg^{-1}]\)
[out]gsnowDiagnostic heat flux at snow surface for use in CCCma black carbon deposition scheme \([W m^{-2}]\)
[out]stDiagnosed screen-level air temperature [K]
[out]suDiagnosed anemometer-level zonal wind \([m s^{-1}]\)
[out]svDiagnosed anemometer-level meridional wind \([m s^{-1}]\)
[out]sqDiagnosed screen-level specific humidity \([kg kg^{-1}]\)
[out]srhDiagnosed screen-level relative humidity [%]
[out]cdhSurface drag coefficient for heat [ ]
[out]cdmSurface drag coefficient for momentum [ ]
[out]qsensDiagnosed total surface sensible heat flux over modelled area \([W m^{-2}]\)
[out]qevapDiagnosed total surface latent heat flux over modelled area \([W m^{-2}]\)
[out]qlwavgUpwelling longwave radiation over modelled area \([W m^{-2}]\)
[out]fsgvDiagnosed net shortwave radiation on vegetation canopy \([W m^{-2}]\)
[out]fsgsDiagnosed net shortwave radiation at snow surface \([W m^{-2}]\)
[out]fsggDiagnosed net shortwave radiation at soil surface \([W m^{-2}]\)
[out]flgvDiagnosed net longwave radiation on vegetation canopy \([W m^{-2}]\)
[out]flgsDiagnosed net longwave radiation at snow surface \([W m^{-2}]\)
[out]flggDiagnosed net longwave radiation at soil surface \([W m^{-2}]\)
[out]hfscDiagnosed sensible heat flux on vegetation canopy \([W m^{-2}]\)
[out]hfssDiagnosed sensible heat flux at snow surface \([W m^{-2}]\)
[out]hfsgDiagnosed sensible heat flux at soil surface \([W m^{-2}]\)
[out]hevcDiagnosed latent heat flux on vegetation canopy \([W m^{-2}]\)
[out]hevsDiagnosed latent heat flux at snow surface \([W m^{-2}]\)
[out]hevgDiagnosed latent heat flux at soil surface \([W m^{-2}]\)
[out]hmfcDiagnosed energy associated with phase change of water on vegetation \([W m^{-2}]\)
[out]hmfnDiagnosed energy associated with phase change of water in snow pack \([W m^{-2}]\)
[out]evppotDiagnosed potential evapotranspiration \([kg m^{-2} s^{-1}]\)
[out]acondDiagnosed product of drag coefficient and wind speed over modelled area \([m s^{-1}]\)
[out]dragSurface drag coefficient under neutral stability [ ]
[out]ilmoSurface drag coefficient under neutral stability [ ]
[out]ueFriction velocity of air \([m s^{-1}]\)
[out]hblHeight of the atmospheric boundary layer [m]
[out]ilmoxInverse of Monin-Obukhov roughness length over each subarea \([m^{-1}]\)
[out]uexFriction velocity of air over each subarea \([m s^{-1}]\)
[out]hblxHeight of the atmospheric boundary layer over each subarea [m]
[in]thliqVolumetric liquid water content of soil layers \([m^3 m^{-3}] (\theta_l)\)
[in]thiceVolumetric frozen water content of soil layers \([m^3 m^{-3}] (\theta_i)\)
[in]tbarTemperature of soil layers [K]
[in]zpondDepth of ponded water on surface [m]
[in]tpondTemperature of ponded water [K]
[in]taAir temperature at reference height [K]
[in]rhosnoDensity of snow \([kg m^{-3}] (\rho_s)\)
[in]tsnowSnowpack temperature [K]
[in]zsnowDepth of snow pack \([m] (z_s)\)
[in]wsnowLiquid water content of snow pack \([kg m^{-2}] (w_s)\)
[in]tcanVegetation canopy temperature [K]
[in]fcFractional coverage of canopy over bare ground for modelled area [ ]
[in]fcsFractional coverage of canopy over snow for modelled area [ ]
[in]thporPore volume in soil layer \([m^3 m^{-3}] (\theta_p)\)
[in]thlminResidual soil liquid water content remaining after freezing or evaporation \([m^3 m^{-3}]\)
[in]thlretLiquid water retention capacity for organic soil \([m^3 m^{-3}] (\theta_{ret})\)
[in]thfcField capacity \([m^3 m^{-3}] (theta_{fc})\)
[in]hcpsHeat capacity of soil material \([J m^{-3} K^{-1}] (C_m)\)
[in]tcsThermal conductivity of soil particles \([W m^{-1} K^{-1}] (\theta_s)\)
[in]delzwPermeable thickness of soil layer \([m] (\Delta z_w)\)
[in]zbotwDepth to permeable bottom of soil layer \([m] (z_{b,w})\)
[in]delzOverall thickness of soil layer [m]
[in]isandSand content flag

In the first two loops, various subarea arrays and internal energyBudgetDriver variables are initialized. The initial temperatures of the vegetation canopy above snow and above bare ground (TCANS and TCANO) are set to the temperature of the vegetation over the whole modelled area (TCAN) if TCAN is not effectively 0 K (the value it is assigned if vegetation is not present). Otherwise, the canopy temperatures are initialized to the air temperature TA.

In loop 200 the soil surface evaporation flag IEVAP and the evaporation efficiency coefficient CEVAP are assigned. If the liquid water content of the first soil layer is effectively equal to the minimum water content THLMIN, IEVAP and CEVAP are set to zero. If the liquid water content of the first soil layer is greater than the soil porosity (THPOR), IEVAP and CEVAP are set to unity. Otherwise, IEVAP is set to 1 and CEVAP (or \(\beta\) as it is typically symbolized in the literature) is calculated using a relation presented by Merlin et al. (2011) [Merlin2011-xy:]

\(\beta = 0.25 [1 – cos(\theta_l \pi / \theta_{p})]^2\)

where \(\theta_l\) is the liquid water content of the first soil layer and \(\theta_{p}\) is its porosity. We follow Merlin et al. in using a 'P' value (their terminology for exponent term) of 2.

In loop 300 the volumetric heat capacities Cg of the soil layers under a bare surface (HCPG) and under vegetation (HCPC) are calculated, from their respective liquid and frozen water contents \(\theta_l\) and \(\theta_i\):

\(C_g = C_w \theta_l + C_i \theta_i + C_m(1 - \theta_p)\)

where \(C_m\) is the heat capacity of the soil matter and \(\theta_p\) is the pore volume. (The heat capacity of air is neglected.)

In loop 400, the thermal properties of the snow pack under the vegetation canopy and over bare soil are assigned on the basis of the properties of the snow pack over the whole modelled area. The heat capacity of the snow pack Cs is calculated from the volume fractions of snow particles and liquid water in the pack. The former is obtained from the ratio of the densities of the snow pack and ice, and the latter from the ratio of the liquid water content, normalized by the snow depth, and the density of water:

\(C_s = C_i [\rho_s /\rho_i ] + C_w w_s /[\rho_w z_s]\)

The thermal conductivity of snow \(\lambda_s\) is obtained from the snow density using an empirical relationship derived by Sturm et al. (1997):

\(\lambda_s = 3.233 x 10^{-6} \rho_s^2 – 1.01 x 10^{-3} \rho_s + 0.138 \rho_s \geq 156.0 \)

\(\lambda_s = 0.234 x 10^{-3} \rho_s + 0.023 \rho_s < 156.0 \)

In loop 500, the thermal conductivities of the soil layers are assigned. If the ISAND flag for the first soil layer is -4 (indicating glacier or ice sheet), or if the ISAND flag is -3 (indicating rock), then literature values for glacier ice or sand particles respectively are assigned. If the ISAND flag is equal to -2, indicating organic soil, the depth of the water table \(z_{wt}\) is first calculated. This is taken to lie within the first layer, counting from the bottom of the soil profile, in which the soil water content is larger than the retention capacity \(\theta_{ret}\). The water table depth is deduced by assuming that the soil is saturated below the water table, and that the water content is at the retention capacity above it. Thus, if \(\theta_l + \theta_i = \theta_p\) for the soil layer, the water table is located at the top of the soil layer; if \(\theta_l + \theta_i = \theta_{ret}\), it is located at the permeable bottom of the soil layer; and if \(\theta_l + \theta_i\) is between these two values, its location is given by:

\(z_{wt} = z_{b, w} - \Delta z_w [(\theta_l + \theta_i - \theta_{ret})/(\theta_p - \theta_{ret})]\)

where \(\Delta z_w\) is the permeable thickness of the soil layer.

The thermal conductivities of organic and mineral soils are calculated following the analysis of Côté and Konrad (2005). They model the soil thermal conductivity \(\lambda\) using the concept of a relative thermal conductivity \(\lambda_r\) which has a value of 0 for dry soils and 1 at saturation:

\(\lambda = [ \lambda_{sat} – \lambda_{dry} ] \lambda_r + \lambda_{dry}\)

The relative thermal conductivity is obtained from the degree of saturation (the water content divided by the pore volume) \(S_r\), using the following generalized relationship:

\(\lambda_r = \kappa S_r/[1 + (\kappa-1) S_r ]\)

The empirical coefficient kappa takes the following values:

Unfrozen coarse mineral soils: \(\kappa = 4.0\)
Frozen coarse mineral soils: \(\kappa = 1.2\)
Unfrozen fine mineral soils: \(\kappa = 1.9\)
Frozen fine mineral soils: \(\kappa = 0.85\)
Unfrozen organic soils: \(\kappa = 0.6\)
Frozen organic soils: \(\kappa = 0.25\)

The dry thermal conductivity \(lambda_{dry}\) is calculated using an empirical relationship based on the pore volume \(\theta_p\), with different coefficients for mineral and organic soils:

\(\lambda_{dry} = 0.75 exp(-2.76 \theta_p) (mineral)\) \(\lambda_{dry} = 0.30 exp(-2.0 \theta_p) (organic)\)

The saturated thermal conductivity \(\lambda_{sat}\) is calculated by Cote and Konrad (2005) [30] as a geometric mean of the conductivities of the soil components. However, other researchers (e.g. Zhang et al., 2008) have found the linear averaging used by de Vries (1963) [94] to be more generally accurate:

\(lambda_{sat} = lambda_w \theta_p + \lambda_s (1 - \theta_p) (unfrozen)\) \(lambda_{sat} = lambda_i \theta_p + \lambda_s (1 - \theta_p) (frozen)\)

where \(\lambda_w\) is the thermal conductivity of water, \(\lambda_i\) is that of ice and \(\lambda_s\) is that of the soil particles.

In the 500 loop, thermal conductivities are calculated for the top and bottom of each soil layer. The degree of saturation SATRAT is calculated as the sum of liquid and frozen water contents, \(\theta_w\) and \(\theta_i\), divided by the pore volume. In organic soils, if the liquid water content of the soil layer is above the retention capacity \(\theta_{ret}\), \(\theta_w\) at the top of the soil layer is assumed to be equal to \(\theta_{re}\) and \(S_r\) at the bottom of the layer is assumed to be 1. The relative liquid and frozen water contents, THLSAT and THISAT, are calculated from \(\theta_w\) and \(\theta_i\) normalized by \(\theta_w + \theta_i\). The dry thermal conductivity, and the saturated thermal conductivity for unfrozen and frozen conditions, are evaluated using the equations above. The unfrozen and frozen relative thermal conductivity, TCRATU and TCRATF, are obtained from SATRAT and the appropriate values of the empirical coefficient \(\kappa\). For mineral soils, \(\kappa\) is obtained as a weighted average over the percent sand content (ISAND converted to a real :: value) and the percentage of fine material (assumed to be 100-ISAND). The unfrozen and frozen soil thermal conductivities, TCSOLU and TCSOLF, are then calculated from TCRATU, TCRATF, and the dry and saturated thermal conductivities; and the actual thermal conductivity of the soil, TCSOIL, is determined as the average of TCSOLU and TCSOLF, weighted according to the relative liquid and frozen water contents THLSAT and THISAT. If the permeable thickness of the layer, DELZW, is greater than zero, the thermal conductivity at the top of the layer is set to TCSOIL; otherwise it is set to the rock value, TCSAND. If DELZW is less than the thermal thickness of the layer DELZ, the thermal conductivity at the bottom of the layer is set to TCSAND; otherwise it is set to TCSOIL. (In the case of organic soils in the latter branch, if \(\theta_w\) was greater than \(\theta_{ret}\), the thermal conductivity at the bottom of the layer is set to the average of the saturated unfrozen and frozen values, weighted by THLSAT and THISAT.) Finally, if there is ponded water present on the soil surface, the thermal conductivity at the top of the first soil layer is treated as varying linearly from the calculated soil thermal conductivity if the pond depth ZPOND is zero, to the thermal conductivity of water if ZPOND \(\geq 10^{-2} m\).

Finally, in loop 600, a variable TBAR1P is evaluated, representing the weighted average value of the first layer soil temperature and the ponded water, if any. (The heat capacity of the soil is determined as the weighted average of HCPG over the permeable thickness DELZW, and the heat capacity of rock, HCPSND, over the impermeable thickness, DELZ-DELZW.)