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

Solves surface energy balance for non-vegetated subareas. More...

Functions/Subroutines

subroutine energbalnovegsolve (ISNOW, FI, QSWNET, QLWOUT, QTRANS, QSENS, QEVAP, EVAP, TZERO, QZERO, GZERO, QMELT, CDH, CDM, RIB, CFLUX, FTEMP, FVAP, ILMO, UE, H, QLWIN, TPOTA, QA, VA, PADRY, RHOAIR, ALVISG, ALNIRG, CRIB, CPHCH, CEVAP, TVIRTA, ZOSCLH, ZOSCLM, ZRSLFH, ZRSLFM, ZOH, ZOM, FCOR, GCONST, GCOEFF, TSTART, PCPR, TRSNOWG, FSSB, ALSNO, THLIQ, THLMIN, DELZW, RHOSNO, ZSNOW, ZPOND, IWATER, IEVAP, ITERCT, ISAND, ISLFD, ITG, ILG, IG, IL1, IL2, JL, NBS, ISNOALB, TSTEP, TVIRTS, EVBETA, Q0SAT, RESID, DCFLXM, CFLUXM, WZERO, TRTOP, A, B, LZZ0, LZZ0T, FM, FH, ITER, NITER, JEVAP, KF, ipeatland, co2conc, pressg, coszs, Cmossmas, dmoss, anmoss, rmlmoss, iday, daylength, pdd)
 

Detailed Description

Solves surface energy balance for non-vegetated subareas.

Author
D. Verseghy, M. Lazare, A. Wu, P. Bartlett, Y. Delage, J. Cole, R. Brown, J. Melton, Y. Wu

Function/Subroutine Documentation

◆ energbalnovegsolve()

subroutine energbalnovegsolve ( integer, intent(in)  ISNOW,
real, dimension (ilg), intent(in)  FI,
real, dimension(ilg), intent(inout)  QSWNET,
real, dimension(ilg), intent(inout)  QLWOUT,
real, dimension(ilg), intent(inout)  QTRANS,
real, dimension (ilg), intent(inout)  QSENS,
real, dimension (ilg), intent(inout)  QEVAP,
real, dimension (ilg), intent(inout)  EVAP,
real, dimension (ilg), intent(inout)  TZERO,
real, dimension (ilg), intent(inout)  QZERO,
real, dimension (ilg), intent(inout)  GZERO,
real, dimension (ilg), intent(inout)  QMELT,
real, dimension (ilg), intent(in)  CDH,
real, dimension (ilg), intent(in)  CDM,
real, dimension (ilg), intent(in)  RIB,
real, dimension (ilg), intent(inout)  CFLUX,
real, dimension (ilg), intent(in)  FTEMP,
real, dimension (ilg), intent(in)  FVAP,
real, dimension (ilg), intent(in)  ILMO,
real, dimension (ilg), intent(in)  UE,
real, dimension (ilg), intent(in)  H,
real, dimension (ilg), intent(in)  QLWIN,
real, dimension (ilg), intent(in)  TPOTA,
real, dimension (ilg), intent(in)  QA,
real, dimension (ilg), intent(in)  VA,
real, dimension (ilg), intent(in)  PADRY,
real, dimension(ilg), intent(in)  RHOAIR,
real, dimension(ilg), intent(in)  ALVISG,
real, dimension(ilg), intent(in)  ALNIRG,
real, dimension (ilg), intent(in)  CRIB,
real, dimension (ilg), intent(in)  CPHCH,
real, dimension (ilg), intent(in)  CEVAP,
real, dimension(ilg), intent(in)  TVIRTA,
real, dimension(ilg), intent(in)  ZOSCLH,
real, dimension(ilg), intent(in)  ZOSCLM,
real, dimension(ilg), intent(in)  ZRSLFH,
real, dimension(ilg), intent(in)  ZRSLFM,
real, dimension (ilg), intent(in)  ZOH,
real, dimension (ilg), intent(in)  ZOM,
real, dimension (ilg), intent(in)  FCOR,
real, dimension(ilg), intent(in)  GCONST,
real, dimension(ilg), intent(in)  GCOEFF,
real, dimension(ilg), intent(in)  TSTART,
real, dimension (ilg), intent(in)  PCPR,
real, dimension(ilg,nbs), intent(in)  TRSNOWG,
real, dimension(ilg,nbs), intent(in)  FSSB,
real, dimension(ilg,nbs), intent(in)  ALSNO,
real, dimension(ilg,ig), intent(in)  THLIQ,
real, dimension(ilg,ig), intent(in)  THLMIN,
real, dimension(ilg,ig), intent(in)  DELZW,
real, dimension(ilg), intent(in)  RHOSNO,
real, dimension(ilg), intent(in)  ZSNOW,
real, dimension (ilg), intent(in)  ZPOND,
integer, dimension(ilg), intent(in)  IWATER,
integer, dimension (ilg), intent(inout)  IEVAP,
integer, dimension(ilg,6,50), intent(inout)  ITERCT,
integer, dimension(ilg,ig), intent(in)  ISAND,
integer, intent(in)  ISLFD,
integer, intent(in)  ITG,
integer, intent(in)  ILG,
integer, intent(in)  IG,
integer, intent(in)  IL1,
integer, intent(in)  IL2,
integer, intent(in)  JL,
integer, intent(in)  NBS,
integer, intent(in)  ISNOALB,
real, dimension (ilg), intent(inout)  TSTEP,
real, dimension(ilg), intent(inout)  TVIRTS,
real, dimension(ilg), intent(inout)  EVBETA,
real, dimension (ilg), intent(inout)  Q0SAT,
real, dimension (ilg), intent(inout)  RESID,
real, dimension(ilg), intent(inout)  DCFLXM,
real, dimension(ilg), intent(inout)  CFLUXM,
real, dimension (ilg), intent(inout)  WZERO,
real, dimension(ilg,nbs), intent(inout)  TRTOP,
real, dimension (ilg), intent(inout)  A,
real, dimension (ilg), intent(inout)  B,
real, dimension (ilg), intent(inout)  LZZ0,
real, dimension (ilg), intent(inout)  LZZ0T,
real, dimension (ilg), intent(inout)  FM,
real, dimension (ilg), intent(inout)  FH,
integer, dimension(ilg), intent(inout)  ITER,
integer, dimension(ilg), intent(inout)  NITER,
integer, dimension (ilg), intent(inout)  JEVAP,
integer, dimension(ilg), intent(inout)  KF,
integer, dimension(ilg), intent(in)  ipeatland,
real, dimension(ilg), intent(in)  co2conc,
real, dimension(ilg), intent(in)  pressg,
real, dimension(ilg), intent(in)  coszs,
real, dimension(ilg), intent(in)  Cmossmas,
real, dimension(ilg), intent(in)  dmoss,
real, dimension(ilg), intent(inout)  anmoss,
real, dimension(ilg), intent(inout)  rmlmoss,
integer, intent(in)  iday,
real, dimension(ilg), intent(in)  daylength,
real, dimension(ilg), intent(inout)  pdd 
)
Parameters
[in]isnowFlag indicating presence or absence of snow
[in]isnoalbSwitch to model snow albedo in two or more wavelength bands
[in,out]qswnetNet shortwave radiation at surface \([W m^{-2}]\)
[in,out]qlwoutUpwelling longwave radiation at surface \([W m^{-2}]\) (L)
[in,out]qtransShortwave radiation transmitted into surface \([W m^{-2}]\)
[in,out]qsensSensible heat flux from surface \([W m^{-2}] (Q_H)\)
[in,out]qevapLatent heat flux from surface \([W m^{-2}] (Q_E)\)
[in,out]evapEvaporation rate at surface \([kg m^{-2} s^{-1}] (E(0))\)
[in,out]tzeroTemperature at surface \([K] (T(0))\)
[in,out]qzeroSpecific humidity at surface \([kg kg^{-1}] (q(0))\)
[in,out]gzeroHeat flux into surface \([W m^{-2}] (G(0))\)
[in,out]qmeltHeat available for melting snow or freezing water at the surface \([W m^{-2}]\)
[in]cdhSurface drag coefficient for heat \([ ] (C_{DH}) \)
[in]cdmSurface drag coefficient for momentum [ ]
[in]ribBulk Richardson number at surface [ ]
[in,out]cfluxProduct of surface drag coefficient and wind speed \([m s^{-1}]\)
[in]ftempProduct of surface-air temperature gradient, drag coefficient and wind speed \([K m s^{-1}]\)
[in]fvapProduct of surface-air humidity gradient, drag coefficient and wind speed \([kg kg^{-1} m s^{-1}]\)
[in]ilmoInverse of Monin-Obukhov roughness length \((m-1]\)
[in]ueFriction velocity of air \([m s^{-1}]\)
[in]hHeight of the atmospheric boundary layer [m]
[in]zpondDepth of ponded water on surface [m]
[in]fiFractional coverage of subarea in question on modelled area [ ]
[in]qlwinDownwelling longwave radiation at bottom of atmosphere \([W m^{-2}]\)
[in]tpotaPotential temperature of air at reference height \([K] (T_{a,pot})\)
[in]qaSpecific humidity at reference height \([kg kg^{-1}] (q_a)\)
[in]vaWind speed at reference height \([m s^{-1}] (v_a)\)
[in]padryPartial pressure of dry air \([Pa] (p_{dry})\)
[in]rhoairDensity of air \([kg m^{-3}] (\rho_a)\)
[in]alvisgVisible albedo of ground surface [ ]
[in]alnirgNear-IR albedo of ground surface [ ]
[in]cribRichardson number coefficient \([K^{-1}]\)
[in]cphchLatent heat of vaporization at surface \([J kg^{-1}]\)
[in]cevapSoil evaporation efficiency coefficient \([ ] (\beta)\)
[in]tvirtaVirtual potential temperature of air at reference height [K]
[in]zosclhRatio of roughness length for heat to reference height for temperature and humidity [ ]
[in]zosclmRatio of roughness length for momentum to reference height for wind speed [ ]
[in]zrslfhDifference between reference height for temperature and humidity and height at which extrapolated wind speed goes to zero [m]
[in]zrslfmDifference between reference height for wind speed and height at which extrapolated wind speed goes to zero [m]
[in]zohSurface roughness length for heat [m]
[in]zomSurface roughness length for momentum [m]
[in]gconstIntercept used in equation relating surface heat flux to surface temperature \([W m^{-2}]\)
[in]gcoeffMultiplier used in equation relating surface heat flux to surface temperature \([W m^{-2} K^{-1}]\)
[in]tstartStarting point for surface temperature iteration [K]
[in]fcorCoriolis parameter \([s^{-1}]\)
[in]pcprSurface precipitation rate \([kg m^{-2} s^{-1}]\)
[in]rhosnoDensity of snow \([kg m^{-3}]\)
[in]zsnowDepth of snow pack [m] \((z_s)\)
[in]thliqVolumetric liquid water content of soil layers \([m^{3} m^{-3}]\)
[in]thlminResidual soil liquid water content remaining after freezing or evaporation \([m^{3} m^{-3}]\)
[in]delzwPermeable thickness of soil layer [m]
[in]trsnowgShort-wave transmissivity of snow pack over bare ground [ ]
[in]alsnoAlbedo of snow in each modelled wavelength band [ ]
[in]fssbTotal solar radiation in each modelled wavelength band \([W m^{-2}]\)
[in]iwaterFlag indicating condition of surface (dry, water-covered or snow-covered)
[in,out]ievapFlag indicating whether surface evaporation is occurring or not
[in,out]iterctCounter of number of iterations required to solve energy balance for four subareas
[in]isandSand content flag

For the surface temperature iteration, two alternative schemes are offered: the bisection method (selected if the flag ITG = 1) and the Newton-Raphson method (selected if ITG = 2). In the first case, the maximum number of iterations ITERMX is set to 50, and in the second case it is set to 5. An optional windless transfer coefficient EZERO is available, which can be used, following the recommendations of Brown et al. (2006) [19], to prevent the sensible heat flux over snow packs from becoming vanishingly small under highly stable conditions. Currently, EZERO is set to zero.

In the next section, calculations of surface net and transmitted shortwave radiation are done depending on whether a snow pack is present (ISNOW=1) and whether the incoming shortwave radiation is being supplied in two bands (one visible and one near-IR, ISNOALB=0) or four bands (one visible and three near-IR, ISNOALB=1). For each band the net shortwave radiation is calculated as the incoming radiation multiplied by one minus the appropriate albedo. If there is no snow present the shortwave transmissivity at the surface is set to zero and the absorbed visible radiation is calculated from the first radiation band. If ISNOALB=0 the absorbed near-IR radiation is calculated from the second band and if ISNOALB=1 it is calculated from the sum of the second to fourth bands. The total net shortwave radiation QSWNET is then evaluated. If a snow pack is present, then if ISNOALB=0 the calculations are done separately for the visible and near-IR bands; if ISNOALB=1, they are done separately over the four wavelength bands. The transmissivity is calculated as an average value for ISNOALB=0, and separately for the four wavelength ranges for ISNOALB=1. The net shortwave radiation at the surface is obtained from the sum of the net values in each wavelength band, corrected for the respective loss by transmission into the surface.

In the 50 loop, the initial value of the surface temperature TZERO is set to TSTART, which contains the value of TZERO from the previous time step, and the first step in the iteration sequence, TSTEP, is set to 1.0 K. The flag ITER is set to 1 for each element of the set of modelled areas, indicating that its surface temperature has not yet been found. The iteration counter NITER is initialized to 1 for each element. Initial values are assigned to several other variables. In particular, the maximum evaporation from the ground surface that can be sustained for the current time step is specified as the total snow mass if snow is present, otherwise as the water ponded on the surface plus the total available water in the first soil layer.

The 100 continuation line marks the beginning of the surface temperature iteration sequence. First the flags NIT (indicating that there are still locations at the beginning of the current iteration step for which the surface temperature has not yet been found) and NUMIT (indicating that there are still locations at the end of the current iteration step for which the surface temperature has not yet been found) are set to zero. Loop 150 is then performed over the set of modelled areas. If ITER=1, NIT is incremented by one, and the initial value of the surface transfer coefficient CFLUXM for this iteration pass is set to its value from the previous pass. The virtual temperature at the surface, \(T(0)_v\), is obtained using the standard expression (see documentation for subroutine energyBudgetDriver):

\(T(0)_v = T(0) [1 + 0.61 q(0)]\)

where T(0) is the surface temperature and q(0) is the specific humidity at the surface. The surface humidity can be obtained from the saturated specific humidity \(q(0)_{sat}\) by making use of the definition of the surface evaporation efficiency \(\beta\):

\(\beta = [q(0) – q_a]/[q(0)_{sat} - q_a]\)

where \(q_a\) is the specific humidity of the air at the reference height. This expression is inverted to obtain an expression for q(0). The saturated specific humidity \(q(0)_{sat}\) is determined from the mixing ratio at saturation, \(w(0)_{sat}\):

\(q(0)_{sat} = w(0)_{sat}/[1 + w(0)_{sat}]\)

The saturation mixing ratio is a function of the saturation vapour pressure \(e(0)_{sat}\) at the surface:

\(w(0)_{sat} = 0.622 e(0)_{sat}/(p_{dry})\)

where \(p_{dry}\) is the partial pressure of dry air. For the saturated vapour pressure, following Emanuel (1994) [34] \(e_{sat}\) is from the temperature \(T_a\) and the freezing point \(T_f\):

\(e_{sat} = exp[53.67957 - 6743.769 / T - 4.8451 * ln(T)] T \geq T_f\)

\(e_{sat} = exp[23.33086 - 6111.72784 / T + 0.15215 * log(T)] T < T_f\)

where \(T_f\) is the freezing point. If there is a snow cover or ponded water present on the surface (IWATER > 0), the surface evaporation efficiency EVBETA is set to 1 and q(0) is set to \(q(0)_{sat}\). Otherwise EVBETA is set to CEVAP, the value obtained in subroutine energyBudgetPrep on the basis of ambient conditions, and q(0) is calculated as above. If \(q(0) > q_a\) and the evaporation flag IEVAP has been set to zero, EVBETA is reset to zero and q(0) is reset to \(q_a\). Finally, \(T(0)_v\) is determined using the equation above.

If NIT > 0, a subroutine is called to evaluate the stability- corrected surface drag coefficients for heat and momentum. The subroutine selection is made on the basis of the flag ISLFD. If ISLFD=1, indicating that the calculations are to be consistent with CCCma conventions, subroutine DRCOEF is called; if ISLFD=2, indicating that the calculations are to be consistent with RPN conventions, subroutine FLXSURFZ is called.

In loop 175, the terms of the surface energy balance are evaluated. The energy balance equation is written as:

\(K_* + L_* - Q_H – Q_E – G(0) = 0\)

where \(K_*\) is the net shortwave radiation, \(L_*\) is the net longwave radiation, \(Q_H\) is the sensible heat flux, \(Q_E\) is the latent heat flux, and G(0) is the conduction into the surface. \(K_*\) was evaluated earlier in loop 50. \(L_*\) is obtained as the difference between the downwelling radiation \(L \downarrow\) and the upwelling radiation \(L \uparrow\), which in turn is determined using the Stefan- Boltzmann equation:

\(L \uparrow = \sigma T(0)^4\)

where \(\sigma\) is the Stefan-Boltzmann constant. (It is assumed that natural surfaces, because of their radiative complexity, act as effective black bodies, so that their emissivity can be taken to be 1.) The sensible heat flux is given by

\(Q_H = [\rho_a c_p C_{DH} v_a + \epsilon_0] [T(0) – T_{a, pot}]\)

where \(\rho_a\) is the density of the air, \(c_p\) is its specific heat, \(C_{DH}\) is the surface drag coefficient for heat, and \(v_a\) and \(T_{a, pot}\) are the wind speed and potential temperature respectively at the reference height. (Note in the code that the variable CFLUX, evaluated in subroutine DRCOEF or FLXSURFZ, represents the product of \(C_{DH}\) and \(v_a\).) The windless transfer coefficient \(\epsilon_0\), evaluated at the beginning of the subroutine, is used only under stable conditions, i.e. if \(T(0) < T_{a, pot}\). The evaporation rate at the surface, E(0), is calculated as

\(E(0) = \rho_a C_{DH} v_a [Q(0) – q_a]\)

The evaporation rate is constrained to be less than or equal to the maximum rate evaluated earlier in the subroutine. \(Q_E\) is obtained by multiplying E(0) by the latent heat of vaporization at the surface. The ground heat flux G(0) is determined as a linear function of T(0) (see documentation for subroutines soilHeatFluxPrep and snowHeatCond). It can be seen that each of the terms of the surface energy balance is a function of a single unknown, T(0) or TZERO. The residual RESID of the energy balance is now evaluated on the basis of the current estimation for TZERO. If the absolute value of RESID is less than \(5.0 W m^{-2}\), or if the absolute value of the iteration step TSTEP most recently used is less than 0.01 K, the surface temperature is deemed to have been found and ITER is set to 0. If the iteration counter NITER is equal to the maximum number and ITER is still 1, ITER is set to -1.

In the following section, the iteration sequence is moved ahead a step. If ITG = 1, then if NIT > 0 and if ITER for the array element in question is 1, the calculations for the bisection method of solution are performed. If NITER = 1 (indicating that this is the first step in the iteration), then if RESID > 0 (indicating that the current value for TZERO had undershot the correct value), TZERO is incremented by 1 K; otherwise it is decremented by 1 K. If this is not the first step in the iteration, then if RESID >0 and TSTEP < 0 (indicating that TZERO has undershot the correct value and the last temperature increment had been a negative one) or if RESID < 0 and TSTEP > 0 (indicating that TZERO has overshot the correct value and the last temperature increment had been a positive one), TSTEP is divided in half and its sign changed. TSTEP is then added to TZERO. The iteration counter NITER and the flag NUMIT are each incremented by one. The next loop contains an optional set of print statements that are executed if ITER=-1, that is if a solution for the surface temperature has not been found within the prescribed maximum number of iteration steps. Finally, if NUMIT > 0, the iteration cycle is repeated from line 100 on.

If ITG = 2, then if NIT > 0 and if ITER for the array element in question is 1, the calculations for the Newton-Raphson method of iteration are performed. In this approach, the value \(x_{n+1}\) used at each iteration step is obtained from the value \(x_n\) at the previous step as follows:

\(x_{n+1} = x_n – f(x_n)/f'(x_n)\)

Identifying \(x_n\) with TZERO and \(f(x_n)\) with the surface energy balance equation, it can be seen that the second term on the right-hand side corresponds to TSTEP; the numerator is equal to RESID and the denominator to the first derivative of the energy balance equation evaluated at TZERO, which in turn is equal to the sum of the derivatives of the individual terms:

\(d( L \uparrow)/dT = -4 \sigma T(0)^3\)

\(d(Q_H)/dT = \rho_a c_p {C_{DH} v_a + [T(0) – T_{a, pot}] d(C_{DH} v_a)/dT}\)

\(d(Q_E)/dT = L_v \rho_a{C_{DH} v_a dq(0)/dT+ [q(0) – q_a] d(C_{DH} v_a)/dT}\)

and dG(0)/dT is equal to the coefficient multiplying TZERO in the equation for G(0). ( \(L_v\) is the latent heat of vaporization at the surface.) At the end of the calculations the iteration counter NITER and the flag NUMIT are each incremented by one, and upon exiting the loop, if NUMIT > 0, the iteration cycle is repeated from line 100 on.

After the iteration has been completed, NUMIT is reset to zero and a check is carried out to ascertain whether convergence has not been reached (i.e. whether ITER = -1) for any location. In such cases it is assumed that conditions of near-neutral stability at the surface are the cause of the difficulty in finding a solution. A trial value of TZERO is calculated using the virtual potential temperature of the air. If RESID > \(50 W m^{-2}\), TZERO is set to this trial value. The values of q(0) and the components of the surface energy balance are recalculated as above, except that \(Q_H\) and \(Q_E\) are assumed as a first approximation to be zero. RESID is determined on this basis. If RESID is positive, it is assigned to \(Q_E\); otherwise RESID is divided equally between \(Q_H\) and \(Q_E\), except in the case of an absolutely dry surface, in which case \(Q_E\) is set to zero and \(Q_H\) to RESID. RESID is reset to zero, E(0) is obtained from \(Q_E\), and \(T(0)_v\) is recalculated. The flag JEVAP for the location is set to 1, and NUMIT is incremented by 1. If NUMIT > 0 at the end of this loop, DRCOEF or FLXSURF are called again, and their calculations are performed for any location where JEVAP is 1, to ensure consistency with the new surface temperature and humidity.

At this point a check is performed for unphysical values of the surface temperature, i.e. for values greater than 100 C or less than -250 C. If such values are encountered, an error message is printed and a call to abort is carried out.

Finally, a check is performed to ensure that TZERO is not less than 0 C if ponded water is present on the surface (IWATER = 1) or greater than 0 C if snow is present on the surface (IWATER = 2), or greater than zero if the surface is an ice sheet (ISAND = -4). If any of these cases is true, TZERO is reset to the freezing point, and q(0) and \(T(0)_v\) are recalculated. DRCOEF or FLXSURFZ are called, and their calculations are performed for all locations meeting these criteria. The components of the surface energy balance are recalculated; the residual amount is assigned to the energy associated with phase change of water at the surface, QMELT, and RESID is set to zero. If QMELT < 0, i.e. if there is freezing taking place, the evaporation flux QEVAP is added to it and then reset to zero (since if ponded water is freezing, it is unavailable for evaporation).

In the last half of the loop, some final adjustments are made to a few variables. If the evaporation flux is vanishingly small, it is added to RESID and reset to zero. If an anomalous case has arisen in which QMELT < 0 over a snow-covered surface or QMELT > 0 over a snow-free surface, QMELT is added to the heat flux into the surface and then reset to zero. Any remaining residual flux is added to \(Q_H\). The shortwave radiation transmitted into the surface is added back to the net shortwave radiation for diagnostic purposes. The surface vapour flux is converted into units of \(m s^{-1}\). Lastly, the iteration counter ITERCT is updated for the level corresponding to the subarea type and the value of NITER.