CLASSIC
Canadian Land Surface Scheme including Biogeochemical Cycles
iceSheetBalance.f90 File Reference

>Performs temperature stepping and surface runoff calculations over ice sheets. More...

Functions/Subroutines

subroutine icesheetbalance (TBAR, TPOND, ZPOND, TSNOW, RHOSNO, ZSNOW, HCPSNO, ALBSNO, HMFG, HTCS, HTC, WTRS, WTRG, GFLUX, RUNOFF, TRUNOF, OVRFLW, TOVRFL, ZPLIM, GGEO, FI, EVAP, R, TR, GZERO, G12, G23, HCP, QMELT, WSNOW, ZMAT, TMOVE, WMOVE, ZRMDR, TADD, ZMOVE, TBOT, DELZ, ISAND, ICONT, IWF, IG, IGP1, IGP2, ILG, IL1, IL2, JL, N)
 

Detailed Description

>Performs temperature stepping and surface runoff calculations over ice sheets.

Author
D. Verseghy, M. Lazare

Function/Subroutine Documentation

◆ icesheetbalance()

subroutine icesheetbalance ( real, dimension (ilg,ig), intent(inout)  TBAR,
real, dimension (ilg), intent(inout)  TPOND,
real, dimension (ilg), intent(inout)  ZPOND,
real, dimension (ilg), intent(inout)  TSNOW,
real, dimension(ilg), intent(inout)  RHOSNO,
real, dimension (ilg), intent(inout)  ZSNOW,
real, dimension(ilg), intent(inout)  HCPSNO,
real, dimension(ilg), intent(inout)  ALBSNO,
real, dimension (ilg,ig), intent(inout)  HMFG,
real, dimension (ilg), intent(inout)  HTCS,
real, dimension (ilg,ig), intent(inout)  HTC,
real, dimension (ilg), intent(inout)  WTRS,
real, dimension (ilg), intent(inout)  WTRG,
real, dimension (ilg,ig), intent(inout)  GFLUX,
real, dimension(ilg), intent(inout)  RUNOFF,
real, dimension(ilg), intent(inout)  TRUNOF,
real, dimension(ilg), intent(inout)  OVRFLW,
real, dimension(ilg), intent(inout)  TOVRFL,
real, dimension (ilg), intent(in)  ZPLIM,
real, dimension (ilg), intent(in)  GGEO,
real, dimension (ilg), intent(in)  FI,
real, dimension (ilg), intent(in)  EVAP,
real, dimension (ilg), intent(in)  R,
real, dimension (ilg), intent(in)  TR,
real, dimension (ilg), intent(in)  GZERO,
real, dimension (ilg), intent(in)  G12,
real, dimension (ilg), intent(in)  G23,
real, dimension (ilg,ig), intent(in)  HCP,
real, dimension (ilg), intent(inout)  QMELT,
real, dimension (ilg), intent(inout)  WSNOW,
real, dimension (ilg,igp2,igp1), intent(inout)  ZMAT,
real, dimension (ilg,igp2), intent(inout)  TMOVE,
real, dimension (ilg,igp2), intent(inout)  WMOVE,
real, dimension (ilg,igp1), intent(inout)  ZRMDR,
real, dimension (ilg), intent(inout)  TADD,
real, dimension (ilg), intent(inout)  ZMOVE,
real, dimension (ilg), intent(inout)  TBOT,
real, dimension (ig), intent(in)  DELZ,
integer, dimension (ilg,ig), intent(in)  ISAND,
integer, dimension (ilg), intent(inout)  ICONT,
integer, intent(in)  IWF,
integer, intent(in)  IG,
integer, intent(in)  IGP1,
integer, intent(in)  IGP2,
integer, intent(in)  ILG,
integer, intent(in)  IL1,
integer, intent(in)  IL2,
integer, intent(in)  JL,
integer, intent(in)  N 
)
Parameters
[in,out]tbarTemperature of ice layer \([C] (T_{av}( \Delta z))\)
[in,out]hmfgEnergy associated with freezing or thawing of water in ice layer \([W m^{-2}]\)
[in,out]htcInternal energy change of ice layer due to conduction and/or change in mass \([W m^{-2}]\) (Ij)
[in,out]gfluxHeat flow between ice layers \([W m^{-2}] (G(\Delta z))\)
[in,out]tpondTemperature of ponded water [C]
[in,out]zpondDepth of ponded water [m]
[in,out]tsnowTemperature of the snow pack [C]
[in,out]rhosnoDensity of snow pack \([kg m^{-3}]\)
[in,out]zsnowDepth of snow pack [m]
[in,out]hcpsnoHeat capacity of snow pack \([J m^{-3} K^{-1}]\)
[in,out]albsnoAlbedo of snow [ ]
[in,out]htcsInternal energy change of ice layer due to conduction and/or change in mass \([W m^{-2}] (I_j) \)
[in,out]wtrsWater transferred into or out of the snow pack \([kg m^{-2} s^{-1}]\)
[in,out]wtrgWater transferred into or out of the ice \([kg m^{-2} s^{-1}]\)
[in,out]runoffTotal runoff from ice column [m]
[in,out]trunofTemperature of total runoff from ice column [K]
[in,out]ovrflwOverland flow from top of ice column [m]
[in,out]tovrflTemperature of overland flow from top of ice column [K]
[in,out]qmeltEnergy available for melting of ice \([W m^{-2}]\)
[in,out]wsnowLiquid water content of snow pack \([kg m^{-2}]\)
[in]fiFractional coverage of subarea in question on modelled area \([ ] (X)\)
[in]evapEvaporation rate from ice surface \([m s^{-1}]\)
[in]rRainfall rate at ice surface \([m s^{-1}]\)
[in]trTemperature of rainfall [C]
[in]gzeroHeat flow into ice surface \([W m^{-2}] (G(0))\)
[in]g12Heat flow between first and second ice layers \([W m^{-2}] (G(\Delta z1))\)
[in]g23Heat flow between second and third ice layers \([W m^{-2}] (G(\Delta z2))\)
[in]zplimLimiting depth of ponded water [m]
[in]ggeoGeothermal heat flux at bottom of modelled ice profile \([W m^{-2}]\)
[in]hcpHeat capacity of ice layer \([J m^{-3} K^{-1}]\)
[in]delzOverall thickness of ice layer \([m] (\Delta z)\)
[in]isandSand content flag

In the 100 loop, any rainfall or snowmelt R reaching the ice surface is added to the ponded water on the surface. The ponded water temperature is calculated as the weighted average of the existing pond and the rainfall or snowmelt added, and the change in internal energy HTC of the first ice layer is updated using the temperature of the added water.

If a full-scale hydrological modelling application is not being run, that is, if only vertical fluxes of energy and moisture are being modelled, the flag IWF will have been pre-set to zero. In this case, overland flow of water is treated using a simple approach: if the ponded depth of water on the soil surface ZPOND exceeds a pre-determined limiting value ZPLIM, the excess is assigned to overland flow. The total runoff from the ice sheet, RUNOFF, is incremented by the excess of the ponded water, and the overland flow for the whole grid cell OVRFLW is incremented by the product of the excess ponded water and the fractional area of the grid cell. The temperature of the overall runoff from the modelled area TRUNOF, and the temperature of the overland flow for the grid cell TOVRFL, are calculated as weighted averages over their previous values and the ponded water temperature TPOND. The internal energy change HTC of the first soil layer is adjusted for the amount of water lost, and ZPOND is set to ZPLIM.

If the temperature of the remaining ponded water is greater than 0 C, the sink of energy required to cool it to 0 C, HCOOL, is calculated and compared with the amount of energy required to warm the first ice layer to 0 C, HWARM. If HWARM > HCOOL, the energy sink of the first layer is used to cool the ponded water to 0 C, and the layer temperature is updated accordingly. Otherwise, the ponded water temperature and the temperature of the first ice layer are both set to 0 C, and the excess energy source given by HCOOL-HWARM is added to the heat available for melting ice, QMELT.

In loop 125, if the temperature of the first ice layer is less than -2 C after the above operations (i.e. if it is not very close to 0 C), and if the ponded water depth is not vanishingly small, freezing of the ponded water can take place. The energy sink required to freeze all of the ponded water, HFREZ, is calculated and compared with HWARM, the amount of energy required to raise the temperature of the first ice layer to 0 C. If HWARM > HFREZ, then HFREZ is converted into an equivalent temperature change using the heat capacity of ice, and added to the temperature of the first ice layer. HFREZ is also used to update HMFG, the diagnosed energy used for phase changes of water in the first ice layer. The internal energy of the first soil layer, HTC, is adjusted to account for the loss of the ponded water, which is assumed to be added to the snow pack. The ponded water is converted into a frozen depth ZFREZ, which is used to update the internal energy of the snow pack, HTCS. If there is not a pre-existing snow pack, the snow albedo is set to the limiting low value of 0.50. The temperature and density of the snow pack are recalculated as weighted averages over the original values and the frozen amount that has been added. ZFREZ is added to the snow depth ZSNOW. The snow heat capacity is recalculated using the new value of the snow density. Finally, the diagnostic amounts of water transferred to the snow pack, WTRS, and from the ice, WTRG, are updated using ZFREZ.

In the 150 loop the heat fluxes between the ice layers are calculated for the optional multiple-layer configuration, in which the standard third layer, normally with a thickness of 3.75 m, can be subdivided into smaller layers and extended to a greater depth if desired. (The heat fluxes at the ice surface, and between the first and second and the second and third layers, were already calculated in the energyBudgetDriver subroutines.) The remaining fluxes are calculated by using a simple linearization of the soil temperature profile. The expression for the ground heat flux at a depth z, G(z), which depends on the thermal conductivity \(\lambda(z)\) and the temperature gradient, is written as:

\(G(z) = \lambda(z) dT(z)/dz\)

The linearized form is thus:

\(G_j = 2 \lambda_i (T_{j-1} – T_j) / (\Delta z_{j-1} + \Delta z_j)\)

where \(G_j\) is the heat flux at the top of layer j, \(T_j\) and \(\Delta z_j\) refer to the temperature and thickness of the layer,and \(\lambda_i\) is the thermal conductivity of ice.

In the 200 loop a value is assigned to the temperature at the bottom of the ice profile, TBOT, and the temperatures for the first three ice layers are stepped ahead. If the standard three-layer configuration is being used, TBOT is obtained by making use of the assumption (see documentation for subroutine soilHeatFluxPrep) that the variation of temperature T with depth z within each soil layer can be modelled using a quadratic equation:

\(T(z) = 1/2 a z^2 + b z +c\)

It can be shown that the temperature at the bottom of a given soil layer, \(T(\Delta z)\), is related to the temperature at the top of the layer T(0) and the heat fluxes at the top and bottom of the layer, G(0) and \(G(\Delta z)\), as follows:

\(T(\Delta z) = T(0) - (\Delta z/ 2 \lambda_i)[G(0) + G(\Delta z)]\)

Making use of the continuity requirement that the heat flux and temperature at the bottom of a given layer must be equal to the heat flux and temperature at the top of the layer beneath it, an expression for the temperature at the bottom of the third ice layer can be obtained as a function of the temperature at the surface and the fluxes between the ice layers:

\(T(\Delta z_3) = T(0) – {G(\Delta z_2) [\Delta z_3 + \Delta z_2] + G(\Delta z_1) [\Delta z_2 + \Delta z_1] + G(0) \Delta z_{21}} / 2 \lambda_i\)

The surface temperature T(0) is obtained by integrating the equation for T(z) to obtain an expression for the average layer temperature \(T_{av}(\Delta z)\), and then inverting this to solve for T(0):

\(T(0) = T_{av}(\Delta z_1) + (\Delta z_1/ 3\lambda_i) [G(0) + 1/2 G(\Delta z_1)]\)

The third layer temperature is then updated using the geothermal flux GGEO. If the optional multiple- layer configuration is used, TBOT is simply set to the temperature of the lowest layer. In either the standard or the multiple-layer case, the first three layer temperatures are updated using the heat fluxes at the surface, between the first and second layers and between the second and third layers, which were determined in the energyBudgetDriver subroutines. Finally, the latter fluxes are assigned to the appropriate levels in the diagnostic GFLUX vector.

In the 250 loop, the GFLUX values determined in the 150 loop are used to update the temperatures in the third and lower ice layers for the multiple-layer configuration. The calculations are bracketed by a determination of the change of internal energy \(I_j\) of the ice layers as a result of the heat fluxes, obtained as the difference in \(I_j\) between the beginning and end of the calculations:

\(\Delta I_j = X_i \Delta[C_i \Delta z_j T_{av}(\Delta z_j)]/ \Delta t\)

where \(C_i\) is the heat capacity of ice, \(Delta t\) is the length of the time step, and \(X_i\) is the fractional coverage of the subarea under consideration relative to the modelled area.

In the 300 loop,checks are carried out to determine whether any of the ice layer temperatures has overshot 0 C as a result of the calculations in the previous loop. If so,the excess energy is assigned to a temporary variable QADD and is also added to the total heat available for melting of the ice,QMELT. QADD is subtracted from the internal energy of the layer in questions and is added to the internal energy of the first layer,since melting is assumed to proceed from the top downwards. The temperature of the layer is reset to 0 C. Finally,the first half of a calculation of the change of internal energy of each ice layer is performed,to be completed at the end of the subroutine.

In the next three loops, the ice layers are adjusted downward to account for the removal of mass at the surface by melting or sublimation. First, the temperature TMOVE of the ice into which the layer is moving is set to the temperature of the layer below it. TMOVE for the bottom layer is set to TBOT. A depth of ice ZMELT is calculated as the amount of ice for which QMELT is sufficient to both raise its temperature to 0 C (if necessary) and melt it. The temperature of the overall runoff and the overland flow are updated as averages of the original temperature and the meltwater temperature (assumed to be at the freezing point), weighted according to their respective amounts, and the overall runoff and overland flow are incremented by ZMELT (with ZMELT converted to an equivalent water depth). The energy used for the melting of ice is calculated from ZMELT and added to HMFG for the first layer, and the amount of energy used to raise the layer temperature to 0 C is added to HTC for the first layer. The total depth of downward adjustment of the ice layers, ZMOVE, is obtained as the sum of ZMELT and the sublimation depth, calculated from EVAP. This amount is added to the diagnostic variable WTRG. Finally, the new temperature of each ice layer is calculated over the layer thickness DELZ as the average of the original temperature weighted by DELZ-ZMOVE and TMOVE weighted by ZMOVE.

In the next loops, the ice layers are adjusted upward to account for addition of mass at the surface by conversion of snow to ice. Snow is converted to ice if the depth of the snow pack exceeds 10 m, or if the density exceeds \(900 kg m^{-3}\) (approaching that of ice). In the first case the excess over and above 10 m is converted; in the second, the whole snow pack is converted. These calculations are performed in the 500 loop, bracketed by a calculation of the change in internal energy of the snow pack, HTCS. In both cases the first level of the ice level movement matrix WMOVE is set to the amount of snow that is converted, expressed as a depth of ice, and the first level of the temperature matrix TMOVE is set to the snow temperature. The amount of converted snow is added to the diagnostic variables WTRS and WTRG. The depth, density and heat capacity of the snow are recalculated. If the entire snow pack is being converted and the water content WSNOW was non-zero, WSNOW is subtracted from WTRS and added to WTRG, and is also added to the total runoff and the overland flow. The runoff and overland flow temperatures are updated accordingly. The snow temperature and water content are reset to zero. The amount of ice that is lost to the bottom of the profile is added to the total runoff, and the runoff temperature is updated accordingly.

In the remaining parts of the code, the actual adjustment of ice layer positions is performed. First the levels of the available depth matrix ZRMDR are set to the ice layer thicknesses; the matrix ZMAT is initialized to zero; and each level J of WMOVE and TMOVE from 2 to the bottom of the ice profile is set to the value of DELZ and TBAR respectively of the J-1 ice level. The ZMAT matrix represents the depth of each ice layer J that is occupied by ice from level K in the layer movement matrix WMOVE after the layer adjustments are complete. In the 700 loop, starting at the top of the ice profile, an attempt is made to assign each layer of WMOVE in turn to the K, J level of ZMAT. If the calculated value of ZMAT is greater than the available depth ZRMDR of the layer, ZMAT is set to ZRMDR, WMOVE is decremented by ZRMDR, and ZRMDR is set to zero. Otherwise the calculated value of ZMAT is accepted, ZRMDR is decremented by ZMAT, and WMOVE is set to zero.

Finally, the 900 loop is performed over each ice layer J to determine the new layer temperature. The temperature adjustment variable TADD is initialized to zero, and then incremented by the temperature TMOVE of each level K weighted by the corresponding K, J level of ZMAT, and finally by the original layer temperature weighted by the corresponding level of ZRMDR. The layer temperature is reset to TADD normalized by DELZ, and the updating of HTC, begun at the end of the 300 loop, is completed.