energBalVegSolve.f90 File Reference

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

Functions/Subroutines
subroutine	energbalvegsolve (ISNOW, FI, QSWNET, QSWNC, QSWNG, QLWOUT, QLWOC, QLWOG, QTRANS, QSENS, QSENSC, QSENSG, QEVAP, QEVAPC, QEVAPG, EVAPC, EVAPG, EVAP, TCAN, QCAN, TZERO, QZERO, GZERO, QMELTC, QMELTG, RAICAN, SNOCAN, CDH, CDM, RIB, TAC, QAC, CFLUX, FTEMP, FVAP, ILMO, UE, H, QFCF, QFCL, HTCC, QSWINV, QSWINI, QLWIN, TPOTA, TA, QA, VA, VAC, PADRY, RHOAIR, ALVISC, ALNIRC, ALVISG, ALNIRG, TRVISC, TRNIRC, FSVF, CRIB, CPHCHC, CPHCHG, CEVAP, TADP, TVIRTA, RC, RBCOEF, ZOSCLH, ZOSCLM, ZRSLFH, ZRSLFM, ZOH, ZOM, FCOR, GCONST, GCOEFF, TGND, TRSNOW, FSNOWC, FRAINC, CHCAP, CMASS, PCPR, FROOT, THLMIN, DELZW, RHOSNO, ZSNOW, IWATER, IEVAP, ITERCT, ISLFD, ITC, ITCG, ILG, IL1, IL2, JL, N, TSTEP, TVIRTC, TVIRTG, EVBETA, XEVAP, EVPWET, Q0SAT, RA, RB, RAGINV, RBINV, RBTINV, RBCINV, TVRTAC, TPOTG, RESID, TCANO, WZERO, XEVAPM, DCFLXM, WC, DRAGIN, CFLUXM, CFLX, IEVAPC, TRTOP, QSTOR, CFSENS, CFEVAP, QSGADD, A, B, LZZ0, LZZ0T, FM, FH, ITER, NITER, KF1, KF2, AILCG, FCANC, CO2CONC, RMATCTEM, THLIQ, THFC, THLW, ISAND, IG, COSZS, PRESSG, XDIFFUS, ICTEM, IC, CO2I1, CO2I2, ctem_on, SLAI, FCANCMX, L2MAX, NOL2PFTS, CFLUXV, ANVEG, RMLVEG, DAYL, DAYL_MAX, ipeatland, Cmossmas, dmoss, anmoss, rmlmoss, iday, pdd)

Detailed Description

Solves surface energy balance for vegetated subareas.

Author

D. Verseghy, M. Lazare, A. Wu, P. Bartlett, Y. Delage, V. Arora, E. Chan, J. Melton, Y. Wu

Function/Subroutine Documentation

energbalvegsolve()

subroutine energbalvegsolve	(	integer, intent(in)	ISNOW,
		real, dimension (ilg), intent(in)	FI,
		real, dimension(ilg), intent(out)	QSWNET,
		real, dimension (ilg), intent(inout)	QSWNC,
		real, dimension (ilg), intent(inout)	QSWNG,
		real, dimension(ilg), intent(out)	QLWOUT,
		real, dimension (ilg), intent(inout)	QLWOC,
		real, dimension (ilg), intent(inout)	QLWOG,
		real, dimension(ilg), intent(inout)	QTRANS,
		real, dimension (ilg), intent(out)	QSENS,
		real, dimension(ilg), intent(inout)	QSENSC,
		real, dimension(ilg), intent(inout)	QSENSG,
		real, dimension (ilg), intent(inout)	QEVAP,
		real, dimension(ilg), intent(inout)	QEVAPC,
		real, dimension(ilg), intent(inout)	QEVAPG,
		real, dimension (ilg), intent(inout)	EVAPC,
		real, dimension (ilg), intent(inout)	EVAPG,
		real, dimension (ilg), intent(inout)	EVAP,
		real, dimension (ilg), intent(inout)	TCAN,
		real, dimension (ilg), intent(inout)	QCAN,
		real, dimension (ilg), intent(inout)	TZERO,
		real, dimension (ilg), intent(inout)	QZERO,
		real, dimension (ilg), intent(inout)	GZERO,
		real, dimension(ilg), intent(inout)	QMELTC,
		real, dimension(ilg), intent(out)	QMELTG,
		real, dimension(ilg), intent(inout)	RAICAN,
		real, dimension(ilg), intent(inout)	SNOCAN,
		real, dimension (ilg), intent(in)	CDH,
		real, dimension (ilg), intent(in)	CDM,
		real, dimension (ilg), intent(in)	RIB,
		real, dimension (ilg), intent(inout)	TAC,
		real, dimension (ilg), intent(inout)	QAC,
		real, dimension (ilg), intent(in)	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(inout)	QFCF,
		real, dimension (ilg), intent(inout)	QFCL,
		real, dimension (ilg), intent(inout)	HTCC,
		real, dimension(ilg), intent(in)	QSWINV,
		real, dimension(ilg), intent(in)	QSWINI,
		real, dimension (ilg), intent(in)	QLWIN,
		real, dimension (ilg), intent(in)	TPOTA,
		real, dimension (ilg), intent(in)	TA,
		real, dimension (ilg), intent(in)	QA,
		real, dimension (ilg), intent(in)	VA,
		real, dimension (ilg), intent(in)	VAC,
		real, dimension (ilg), intent(in)	PADRY,
		real, dimension(ilg), intent(in)	RHOAIR,
		real, dimension(ilg), intent(in)	ALVISC,
		real, dimension(ilg), intent(in)	ALNIRC,
		real, dimension(ilg), intent(in)	ALVISG,
		real, dimension(ilg), intent(in)	ALNIRG,
		real, dimension(ilg), intent(in)	TRVISC,
		real, dimension(ilg), intent(in)	TRNIRC,
		real, dimension (ilg), intent(in)	FSVF,
		real, dimension (ilg), intent(in)	CRIB,
		real, dimension(ilg), intent(inout)	CPHCHC,
		real, dimension(ilg), intent(in)	CPHCHG,
		real, dimension (ilg), intent(in)	CEVAP,
		real, dimension (ilg), intent(in)	TADP,
		real, dimension(ilg), intent(in)	TVIRTA,
		real, dimension (ilg), intent(inout)	RC,
		real, dimension(ilg), intent(in)	RBCOEF,
		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)	TGND,
		real, dimension(ilg), intent(in)	TRSNOW,
		real, dimension(ilg), intent(inout)	FSNOWC,
		real, dimension(ilg), intent(inout)	FRAINC,
		real, dimension (ilg), intent(inout)	CHCAP,
		real, dimension (ilg), intent(in)	CMASS,
		real, dimension (ilg), intent(in)	PCPR,
		real, dimension (ilg,ig), intent(inout)	FROOT,
		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,
		integer, dimension(ilg), intent(in)	IWATER,
		integer, dimension (ilg), intent(inout)	IEVAP,
		integer, dimension(ilg,6,50), intent(inout)	ITERCT,
		integer, intent(in)	ISLFD,
		integer, intent(in)	ITC,
		integer, intent(in)	ITCG,
		integer, intent(in)	ILG,
		integer, intent(in)	IL1,
		integer, intent(in)	IL2,
		integer, intent(in)	JL,
		integer, intent(in)	N,
		real, dimension (ilg), intent(inout)	TSTEP,
		real, dimension(ilg), intent(inout)	TVIRTC,
		real, dimension(ilg), intent(inout)	TVIRTG,
		real, dimension(ilg), intent(inout)	EVBETA,
		real, dimension (ilg), intent(inout)	XEVAP,
		real, dimension(ilg), intent(inout)	EVPWET,
		real, dimension (ilg), intent(inout)	Q0SAT,
		real, dimension (ilg), intent(inout)	RA,
		real, dimension (ilg), intent(inout)	RB,
		real, dimension(ilg), intent(inout)	RAGINV,
		real, dimension (ilg), intent(inout)	RBINV,
		real, dimension(ilg), intent(inout)	RBTINV,
		real, dimension(ilg), intent(inout)	RBCINV,
		real, dimension(ilg), intent(inout)	TVRTAC,
		real, dimension (ilg), intent(inout)	TPOTG,
		real, dimension (ilg), intent(inout)	RESID,
		real, dimension (ilg), intent(inout)	TCANO,
		real, dimension (ilg), intent(inout)	WZERO,
		real, dimension(ilg), intent(inout)	XEVAPM,
		real, dimension(ilg), intent(inout)	DCFLXM,
		real, dimension (ilg), intent(inout)	WC,
		real, dimension(ilg), intent(inout)	DRAGIN,
		real, dimension(ilg), intent(inout)	CFLUXM,
		real, dimension (ilg), intent(inout)	CFLX,
		integer, dimension(ilg), intent(inout)	IEVAPC,
		real, dimension (ilg), intent(inout)	TRTOP,
		real, dimension (ilg), intent(inout)	QSTOR,
		real, dimension(ilg), intent(inout)	CFSENS,
		real, dimension(ilg), intent(inout)	CFEVAP,
		real, dimension(ilg), intent(inout)	QSGADD,
		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)	KF1,
		integer, dimension(ilg), intent(inout)	KF2,
		real, dimension(ilg,ictem), intent(in)	AILCG,
		real, dimension(ilg,ictem), intent(in)	FCANC,
		real, dimension(ilg), intent(in)	CO2CONC,
		real, dimension(ilg,ictem,ig), intent(in)	RMATCTEM,
		real, dimension(ilg,ig), intent(in)	THLIQ,
		real, dimension(ilg,ig), intent(in)	THFC,
		real, dimension(ilg,ig), intent(in)	THLW,
		integer, dimension(ilg,ig), intent(in)	ISAND,
		integer, intent(in)	IG,
		real, dimension(ilg), intent(in)	COSZS,
		real, dimension(ilg), intent(in)	PRESSG,
		real, dimension(ilg), intent(in)	XDIFFUS,
		integer, intent(in)	ICTEM,
		integer, intent(in)	IC,
		real, dimension(ilg,ictem), intent(in)	CO2I1,
		real, dimension(ilg,ictem), intent(in)	CO2I2,
		logical, intent(in)	ctem_on,
		real, dimension(ilg,ictem), intent(in)	SLAI,
		real, dimension(ilg,ictem), intent(in)	FCANCMX,
		integer, intent(in)	L2MAX,
		integer, dimension(ic), intent(in)	NOL2PFTS,
		real, dimension(ilg), intent(inout)	CFLUXV,
		real, dimension(ilg,ictem), intent(in)	ANVEG,
		real, dimension(ilg,ictem), intent(in)	RMLVEG,
		real, dimension(ilg), intent(in)	DAYL,
		real, dimension(ilg), intent(in)	DAYL_MAX,
		integer, dimension(ilg), intent(in)	ipeatland,
		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(inout)	pdd
	)

Parameters

[in]	isnow	Flag indicating presence or absence of snow
[out]	qswnet	Total net shortwave radiation of canopy and underlying surface \([W m^{-2} ]\)
[in,out]	qswnc	Net shortwave radiation on vegetation canopy \([W m^{-2} ] (K_{*c} )\)
[in,out]	qswng	Net shortwave radiation at underlying surface \([W m^{-2} ] (K_{*g} )\)
[out]	qlwout	Upwelling longwave radiation from canopy and underlying surface \([W m^{-2} ]\)
[in,out]	qlwoc	Upwelling longwave radiation from vegetation canopy \([W m^{-2} ] (L \uparrow_c, L \downarrow_c)\)
[in,out]	qlwog	Upwelling longwave radiation from underlying surface \([W m^{-2} ] (L \uparrow_g)\)
[in,out]	qtrans	Shortwave radiation transmitted into surface \([W m^{-2} ] \)
[out]	qsens	Sensible heat flux from canopy and underlying surface \([W m^{-2} ] (Q_H)\)
[in,out]	qsensc	Sensible heat flux from vegetation canopy \([W m^{-2} ] (Q_{H,c} )\)
[in,out]	qsensg	Sensible heat flux from underlying surface \([W m^{-2} ] (Q_{H,g} )\)
[in,out]	qevap	Latent heat flux from canopy and underlying surface \([W m^{-2} ] (Q_E)\)
[in,out]	qevapc	Latent heat flux from vegetation canopy \([W m^{-2} ] (Q_{E,c} )\)
[in,out]	qevapg	Latent heat flux from underlying surface \([W m^{-2} ] (Q_{E,g} )\)
[in,out]	evapc	Evaporation rate from vegetation \([kg m^{-2} s^{-1} ] (E_c)\)
[in,out]	evapg	Evaporation rate from underlying surface \([kg m^{-2} s^{-1} ] (E(0))\)
[in,out]	tcan	Vegetation canopy temperature \([K] (T_c)\)
[in,out]	qcan	Saturated specific humidity at canopy temperature \([kg kg^{-1} ] (q_c)\)
[in,out]	tzero	Temperature at surface \([K] (T(0))\)
[in,out]	qzero	Specific humidity at surface \([kg kg^{-1} ] (q(0))\)
[in,out]	gzero	Heat flux into surface \([W m^{-2} ] (G(0))\)
[in,out]	qmeltc	Heat available for melting snow or freezing water on the vegetation \([W m^{-2} ]\)
[out]	qmeltg	Heat available for melting snow or freezing water on the underlying surface \([W m^{-2} ]\)
[in,out]	raican	Intercepted liquid water stored on vegetation canopy \([kg m^{-2} ]\)
[in,out]	snocan	Intercepted frozen water stored on vegetation canopy \([kg m^{-2} ]\)
[in]	cdh	Surface drag coefficient for heat [ ]
[in]	cdm	Surface drag coefficient for momentum [ ]
[in]	rib	Bulk Richardson number at surface [ ]
[in,out]	tac	Temperature of air within vegetation canopy \([K] (T_{ac} )\)
[in,out]	qac	Specific humidity of air within vegetation canopy space \([kg kg^{-1} ] (q_{ac} )\)
[in]	cflux	Product of surface drag coefficient and wind speed \([m s^{-1} ]\)
[in]	ftemp	Product of surface-air temperature gradient, drag coefficient and wind speed \([K m s^{-1} ]\)
[in]	fvap	Product of surface-air humidity gradient, drag coefficient and wind speed \([kg kg^{-1} m s^{-1} ]\)
[in]	ilmo	Inverse of Monin-Obukhov roughness length \((m^{-1} ]\)
[in]	ue	Friction velocity of air \([m s^{-1} ]\)
[in]	h	Height of the atmospheric boundary layer [m]
[in,out]	qfcf	Sublimation from frozen water on vegetation \([kg m^{-2} s^{-1} ]\)
[in,out]	qfcl	Evaporation from liquid water on vegetation \([kg m^{-2} s^{-1} ]\)
[in,out]	htcc	Internal energy change of canopy due to changes in temperature and/or mass \([W m^{-2} ]\)
[in,out]	evap	Diagnosed total surface water vapour flux over modelled area \([kg m^{-2} s^{-1} ]\)
[in]	fi	Fractional coverage of subarea in question on modelled area [ ]
[in]	qswinv	Visible radiation incident on horizontal surface \([W m^{-2} ] (K \downarrow)\)
[in]	qswini	Near-infrared radiation incident on horizontal surface \([W m^{-2} ] (K \downarrow)\)
[in]	qlwin	Downwelling longwave radiation at bottom of atmosphere \([W m^{-2} ] (L \downarrow)\)
[in]	tpota	Potential temperature of air at reference height \([K] (T_{a,pot} )\)
[in]	ta	Air temperature at reference height [K]
[in]	qa	Specific humidity at reference height \([kg kg^{-1} ] (q_a)\)
[in]	va	Wind speed at reference height \([m s^{-1} ]\)
[in]	vac	Wind speed within vegetation canopy space \([m s^{-1} ] (v_{ac} )\)
[in]	padry	Partial pressure of dry air \([Pa] (p_{dry} )\)
[in]	rhoair	Density of air \([kg m^{-3} ] ( \rho_a)\)
[in]	alvisc	Visible albedo of vegetation canopy \([ ] ( \alpha_c)\)
[in]	alnirc	Near-IR albedo of vegetation canopy \([ ] ( \alpha_c)\)
[in]	alvisg	Visible albedo of underlying surface \([ ] ( \alpha_g)\)
[in]	alnirg	Near-IR albedo of underlying surface \([ ] ( \alpha_g)\)
[in]	trvisc	Visible transmissivity of vegetation canopy \([ ] ( \tau_c)\)
[in]	trnirc	Near-IR transmissivity of vegetation canopy \([ ] ( \tau_c)\)
[in]	fsvf	Sky view factor of ground underlying canopy \([ ] ( \chi)\)
[in]	crib	Richardson number coefficient \([K^{-1} ]\)
[in,out]	cphchc	Latent heat of vaporization on vegetation canopy \([J kg^{-1} ] (L_v)\)
[in]	cphchg	Latent heat of vaporization on underlying ground \([J kg^{-1} ] (L_v)\)
[in]	cevap	Soil evaporation efficiency coefficient [ ]
[in]	tadp	Dew point of air at reference height [K]
[in]	tvirta	Virtual potential temperature of air at reference height \([K] (T_{a,v} )\)
[in,out]	rc	Stomatal resistance of vegetation \([s m^{-1}] (r_c)\)
[in]	rbcoef	Parameter for calculation of leaf boundary resistance
[in]	zosclh	Ratio of roughness length for heat to reference height for temperature and humidity [ ]
[in]	zosclm	Ratio of roughness length for momentum to reference height for wind speed [ ]
[in]	zrslfh	Difference between reference height for temperature and humidity and height at which extrapolated wind speed goes to zero [m]
[in]	zrslfm	Difference between reference height for wind speed and height at which extrapolated wind speed goes to zero [m]
[in]	zoh	Surface roughness length for heat [m]
[in]	zom	Surface roughness length for momentum \([m] (z_{0,m} )\)
[in]	fcor	Coriolis parameter \([s_{-1} ]\)
[in]	gconst	Intercept used in equation relating surface heat flux to surface temperature \([W m^{-2} ]\)
[in]	gcoeff	Multiplier used in equation relating surface heat flux to surface temperature \([W m^{-2} K^{-1} ]\)
[in]	tgnd	Starting point for surface temperature iteration [K]
[in]	trsnow	Short-wave transmissivity of snow pack [ ]
[in,out]	fsnowc	Fractional coverage of canopy by frozen water \([ ] (F_s)\)
[in,out]	frainc	Fractional coverage of canopy by liquid water \([ ] (F_l)\)
[in,out]	chcap	Heat capacity of vegetation canopy \([J m^{-2} K^{-1} ] (C_c)\)
[in]	cmass	Mass of vegetation canopy \([kg m^{-2} ]\)
[in]	pcpr	Surface precipitation rate \([kg m^{-2} s^{-1} ]\)
[in]	rhosno	Density of snow \([kg m^{-3} ]\)
[in]	zsnow	Depth of snow pack [m] \((z_s)\)
[in,out]	froot	Fraction of total transpiration contributed by soil layer [ ]
[in]	thlmin	Residual soil liquid water content remaining after freezing or evaporation \([m^{3} m^{-3} ]\)
[in]	thliq	Volumetric liquid water content of soil layers \([m^{3} m^{-3} ]\)
[in]	delzw	Permeable thickness of soil layer [m]
[in]	iwater	Flag indicating condition of surface (dry, water-covered or snow-covered)
[in,out]	ievap	Flag indicating whether surface evaporation is occurring or not output variable
[in,out]	iterct	Counter of number of iterations required to solve energy balance for four subareas

For the subcanopy surface temperature iteration, two alternative schemes are offered: the bisection method (selected if the flag ITCG = 1) and the Newton-Raphson method (selected if ITCG = 2). In the first case, the maximum number of iterations ITERMX is set to 12, and in the second case it is set to 5. An optional windless transfer coefficient EZERO is made 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. If the snow cover flag ISNOW is zero (indicating bare ground), EZERO is set to zero; if ISNOW=1, EZERO is set to \(2.0 W m^{-2} K^{-1}\) . The surface transfer coefficient under conditions of free convection, RAGCO, is set to \(1.9 x 10^{-3}\) (see the calculation of RAGINV below).

In the 50 loop, some preliminary calculations are done. The shortwave transmissivity at the surface, TRTOP, is set to zero in the absence of a snow pack, and to the transmissivity of snow, TRSNOW, otherwise. The net shortwave radiation at the surface, QSWNG, is calculated as the sum of the net visible and net near-infrared shortwave radiation. Each of these is obtained as: \(K_{*g} = K \downarrow \tau_c [1 - \alpha_g ]\) where \(K_{*g} is the net radiation at the surface, \)K \( is the incoming shortwave radiation above the canopy, \)$ is the canopy transmissivity and \(\alpha_g\) is the surface albedo. This average value is corrected for the amount of radiation transmitted into the surface, QTRANS, obtained using TRTOP. The net shortwave radiation for the vegetation canopy, QSWNC, is calculated as the sum of the net visible and net near-infrared shortwave radiation. Each of these is determined as: \(K_{*c} = K \downarrow [1 - \alpha_c ] - K_{*g}\) where \(K_{*c}\) is the net radiation on the canopy and \(\alpha_c\) is the canopy albedo. If the canopy temperature is essentially 0 K, indicating that a canopy was not present in the previous time step but has now appeared, the canopy temperature is initialized to the potential temperature of the air, TPOTA. The outgoing longwave radiation emitted upward \((L \uparrow_c)\) or downward \((L \downarrow_c)\) from the canopy is calculated using the standard Stefan-Boltzmann equation: \(L \uparrow_c = L \downarrow_c = \sigma T_c^4\) where \(\sigma\) is the Stefan-Boltzmann constant and \(T_c\) is the canopy temperature.

Virtual temperature is defined as temperature adjusted for the reduction in air density due to the presence of water vapour. This is applied in order to enable the use of the equation of state for dry air. The virtual temperature can be approximated by multiplying the actual temperature by a factor [1 + 0.61 q], where q is the specific humidity of the air in question. Thus, the virtual temperature of air at the vegetation canopy temperature, \(T_{c, v}\) , is obtained as: \(T_{c, v} = T_c [1 + 0.61 q_c ]\) where \(q_c\) is the saturated specific humidity at the canopy temperature. This is determined from the saturation mixing ratio at the canopy temperature, \(w_c\) : \(q_c = w_c /[1 + w_c ]\) The saturation mixing ratio is a function of the saturation vapour pressure \(e_c\) at the canopy temperature: \(w_c = 0.622 e_c /(p_{dry} )\) where \(p_{dry}\) is the partial pressure of dry air. A standard empirical equation for the saturation vapour pressure dependence on the temperature T is used following Emanuel (1994) [34] \(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. The virtual temperature of the air in the canopy space, \(T_{ac, v}\) , is likewise calculated from the canopy air temperature \(T_{ac}\) and the specific humidity in the canopy air space, \(q_{ac}\) , as

\(T_{ac, v} = T_{ac} [1 + 0.61 q_{ac} ]\)

If the Newton-Raphson method is being used for the canopy temperature iteration (pre-selected by setting the flag ITC to 2), the temperature of the air in the canopy space is approximated as the canopy temperature, and the specific humidity in the canopy air space is approximated as the specific humidity of the air above the canopy.

If there is intercepted snow on the vegetation, the latent heat associated with water flux from the canopy, CPHCHC, is set to the latent heat of sublimation (the sum of the heat of vaporization and the heat of melting). Otherwise the latent heat is set to that of vaporization alone. The leaf boundary layer resistance RB, and its inverse RBINV, are calculated using the wind speed in the canopy air space, VAC, and a coefficient RBCOEF evaluated in subroutine calcLandSurfParams. This coefficient is formulated after Bartlett (2004), who developed an expression for the inverse of the leaf boundary resistance, \(1/r_b\) , drawing on the analysis of Bonan (1996) [18] and McNaughton and van den Hurk (1995), of the form:

\(1/r_b = v_{ac}^{1/2} \sigma f_i \gamma_i \Lambda_i^{1/2} /0.75 [1 - exp(-0.75 \Lambda_i^{1/2})]\)

where \(v_{ac}\) is the wind speed in the canopy air space, \(f_i\) is the fractional coverage of each vegetation type i over the subarea in question, \(\Lambda_i\) is its leaf area index, and \(\gamma_i\) is a vegetation-dependent parameter which incorporates the effects of leaf dimension and sheltering. The initial value of the surface temperature TZERO is set to TGND, which contains the value of TZERO from the previous time step, and the initial temperature of the canopy, TCANO, is set to the current canopy temperature. 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 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 total available water in the first soil layer.

(At this point in the code, if CLASS is being run in coupled mode with CTEM, the CTEM subroutine photosynCanopyConduct is called. These lines are commented out for uncoupled runs.)

The 100 continuation line marks the beginning of the surface temperature iteration sequence. First the flag NUMIT (indicating whether there are still locations at the end of the current iteration step for which the surface temperature has not yet been found) is set to zero. Loop 125 is then performed over the vector of modelled areas. If ITER=1, the saturated specific humidity at the ground surface temperature, \(q(0)_{sat}\) , is calculated using the same method as that outlined above for the saturated specific humidity at the canopy temperature. 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 the surface specific humidity 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 from \(q(0)_{sat}\) by making use of the definition of the surface evaporation efficiency \(\beta\):

\(\beta = [q(0) - q_{ac} ]/[q(0)_{sat} - q_{ac} ]\)

which is inverted to obtain an expression for q(0). If \(q(0) > q_{ac}\) and the evaporation flag IEVAP has been set to zero, EVBETA is reset to zero and q(0) is reset to \(q_{ac}\).

Next, the potential temperature of air at the ground surface temperature, \(T(0)_{pot}\) , is calculated relative to the sum of the displacement height d and the roughness length for momentum \(z_{0, m}\) , the height at which the canopy temperature is assumed to apply. (This is also the height corresponding to the potential temperature of the air at the reference height above the canopy). \(T(0)_{pot}\) is found using the dry adiabatic lapse rate, \(dT/dz = -g/c_p\) , where g is the acceleration due to gravity and \(c_p\) is the specific heat at constant pressure. Since the displacement height d is assumed to lie at 0.7 of the canopy height, and the roughness length for momentum is assumed to lie at 0.1 of the canopy height, their sum can be obtained as \(8.0 z_{0, m}\) . Thus, \(T(0)_{pot} = T(0) - 8.0 z_{0, m} g/c_p\)

The virtual potential temperature at the surface is obtained using the same expression as above: \(T(0)_v = T(0)_{pot} [1 + 0.61 q(0)]\)

Since wind speeds under the canopy are expected to be relatively small, it is assumed that under stable or neutral conditions, turbulent fluxes are negligible. If the virtual potential temperature of the surface is more than 1 K greater than that of the canopy air, free convection conditions are assumed. Townsend's (1964) equation for the surface-air transfer coefficient, or the inverse of the surface resistance \(r_{a, , g}\) , is used in a form derived from the analysis of Deardorff (1972) [31] : \(1/r_{a, , g} = 1.9 x 10^{-3} [T(0)_v - T_{ac, v} ]^{1/3}\)

The first derivative of the transfer coefficient with respect to the surface temperature is calculated for use with the Newton-Raphson iteration scheme: \(d(1/r_{a, , g} )/dT = 1.9 x 10^{-3} [T(0)_v - T_{ac, v} ]^{-2/3} /3\)

If the virtual potential temperature of the surface is between 0.001 and 1 K greater than that of the canopy air, a simple diffusion relation is assumed: \(1/r_{a, , g} = 1.9 x 10^{-3} [T(0)_v - T_{ac, v} ]\) Thus, \(d(1/r_{a, , g} )/dT = 1.9 x 10^{-3}\)

The remaining terms of the surface energy balance are now evaluated. The energy balance equation is expressed as: \(K_{*g} + L_{*g} - Q_{H, g} - Q_{E, g} - G(0) = 0\) where \(K_{*g}\) is the net surface shortwave radiation, \(L_{*g}\) is the net longwave radiation, \(Q_{H, g}\) is the sensible heat flux, \(Q_{E, g}\) is the latent heat flux, and G(0) is the conduction into the surface. \(K_{*g}\) was evaluated earlier in loop 50. \(L_{*g}\) is calculated as the difference between the downwelling radiation \(L \downarrow_g\) at the surface and the upwelling radiation \(L \uparrow_g\) . The downwelling radiation incident on the surface is determined from the downwelling sky radiation above the canopy \(L \downarrow\) and the downwelling radiation from the canopy itself, \(L \downarrow_c\) , weighted according to the sky view factor \(\chi\): \(L \downarrow_g = \chi L \downarrow + [1 - \chi] L \downarrow_c\) The upwelling radiation is obtained using the Stefan-Boltzmann equation: \(L \uparrow_g = \sigma T(0)^4\)

The sensible heat flux is given by \(Q_{H, g} = \rho_a c_p [T(0)_{pot} - T_{a, c} ]/r_{a, , g}\) where \(\rho_a\) is the density of the air and \(c_p\) is its specific heat. (The windless transfer coefficient defined at the beginning of the subroutine is currently not used under vegetation canopies.) The evaporation rate at the surface, E(0), is calculated as \(E(0) = \rho_a [q(0) - q_{a, c} ]/r_{a, , g}\)

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 heat flux into the surface 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 ITCG = 1, the calculations for the bisection method of solution in loop 150 are performed, over each of the modelled areas for which ITER = 1. 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. Finally, if NUMIT > 0, the iteration cycle is repeated from line 100 on.

If ITCG = 2, the calculations for the Newton-Raphson method of iteration in loop 175 are performed, over each of the modelled areas for which ITER = 1. 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_g)/dT = -4 \sigma T(0)^3 \]

\[ d(Q_{H, g} )/dT = \rho_a c_p {1/r_{a, , g} + [T(0)_{pot} - T_{ac} ] d(1/r_{a, , g} )/dT} \]

\[ d(Q_{E, g} )/dT = L_v \rho_a {1/r_{a, , g} dq(0)/dT+ [q(0) - q_{ac} ] d(1/r_{a, , g} )/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, if the Newton-Raphson method has been used, a check is carried out in loop 200 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 canopy. If the absolute value of RESID is \(> 15 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, g}\) and \(Q_{E, g}\) are assumed as a first approximation to be zero. RESID is determined on this basis, and is then partitioned between \(Q_{H, g}\) and \(Q_{E, g}\) on the basis of the Bowen ratio B, the ratio of \(Q_{H, g}\) over \(Q_{E, g}\) . Setting the residual R equal to \(Q_{H, g} + Q_{E, g}\) , and substituting for \(Q_{H, g}\) using \(B = Q_{H, g} /Q_{E, g}\) , results in: \(Q_{E, g} = R/(1 + B)\) \(Q_{H, g}\) is then obtained as \(R - Q_{E, g}\) , the residual is reset to zero, and E(0) is recalculated from \(Q_{E, g}\) .

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 -100 C. If such values are encountered, an error message is printed and a call to abort is carried out.

Finally, clean-up calculations are performed in loop 250. A check is carried out 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). If either is the case, TZERO is reset to the freezing point, and q(0), \(T(0)_{pot}\) and \(T(0)_v\) are re-evaluated. The components of the surface energy balance are recalculated using the above equations. The residual of the energy balance equation is assigned to the energy associated with phase change of water at the surface, QMELTG, and RESID is set to zero.

In the last part of the loop, some final adjustments are made to a few other variables. If the evaporation flux is vanishingly small, it is added to RESID and reset to zero. Any remaining RESID is added to \(Q_{H,g}\) . Lastly, the iteration counter ITERCT is updated for the level corresponding to the subarea type and the value of NITER.

In the 300 loop, preliminary calculations are done in preparation for the canopy temperature iteration. The sensible heat flux from the surface is treated differently if ITC = 1 (bisection method of solution) and ITC = 2 (Newton-Raphson method of solution). In the first instance the surface sensible heat flux is applied to heating the air in the canopy space; in the second, the canopy and the air space within it are treated as one aggregated mass, and the sensible heat flux from below is assumed to be added to its energy balance. Thus, if ITC = 1, the sensible heat flux that is added to the canopy, QSGADD, is set to 0. If ITC = 2, it is set to \(Q_{H, g}\) ; and \(T_{ac}\) is set to \(T_c\) , \(q_{ac}\) to \(q_c\) , and \(T_{ac, v}\) to \(T_{c, v}\) . 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. The first step in the iteration sequence, TSTEP, is set to 1.0 K. Initial values are assigned to other variables. In the following 350 loop, the water available for transpiration, WAVAIL, is evaluated for each soil layer. Lastly, after exiting the loop, the maximum number of iterations ITERMX is set to 50 if ITC = 1, and to 5 if ITC = 2.

The 400 continuation line marks the beginning of the canopy temperature iteration sequence. First the flags NIT (indicating whether 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 whether 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 450 is then performed over the set of modelled areas. If ITER=1, NIT is incremented by one; if ITC = 1, the vegetation virtual temperature \(T_{c, v}\) is recalculated.

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.

Next, canopy parameters and turbulent transfer coefficients are calculated prior to evaluating the terms of the surface energy balance equation. If ITC = 1, the analysis of Garratt (1992) [39] is followed. The sensible and latent heat fluxes from the canopy air to the overlying atmosphere, \(Q_H\) and \(Q_E\) , are obtained as the sums of the sensible and latent heat fluxes from the canopy to the canopy air, \(Q_{H, c}\) and \(Q_{E, c}\) , and from the underlying surface to the canopy air, \(Q_{H, g}\) and \(Q_{E, g}\) :

\(Q_H = Q_{H, c} + Q_{H, g}\)

\(Q_E = Q_{E, c} + Q_{E, g}\)

where

\(Q_H = \rho_a c_p [T_{a, c} - T_{a, pot, } ]/r_a\)

\(Q_{H, c} = \rho_a c_p [T_c - T_{a, c, } ]/r_b\)

and

\(Q_E = L_v \rho_a [q_{a, c} - q_{a, } ]/r_a\)

\(Q_{E, c} = L_v \rho_a [q_{, c} - q_{a, c, } ]/(r_b + r_c)\)

The equations for the sensible and latent heat fluxes from the surface were presented above. In these expressions, \(T_{a, pot}\) and \(q_a\) are the potential temperature and specific humidity of the air overlying the canopy, \(r_a\) is the aerodynamic resistance to turbulent transfer between the canopy air and the overlying air, and \(r_c\) is the stomatal resistance to transpiration. (The term CFLUX that is generated by the subroutines DRCOEF and FLXSURFZ is equivalent to \(1/r_a\) .) Thus, \(T_{a, c}\) and \(q_{a, c}\) can be evaluated as

\(T_{a, c} = [T_{a, pot} /r_a + T_c /r_b + T(0)_{pot} /r_{a, , g} ]/[1/r_a + 1/r_b + 1/r_{a, , g} ]\)

\(q_{a, c} = [q_a /r_a + q_c /(r_b + r_c) + q(0)/r_{a, , g} ]/[1/r_a + 1/(r_b + r_c) + 1/r_{a, , g} ]\)

If the water vapour flux is towards the canopy leaves, or if the canopy is snow-covered or rain-covered, \(r_c\) is zero. If the water vapour flux is away from the canopy leaves and the fractional snow or water coverage on the canopy, \(F_s\) or \(F_l\) , is less than 1, the term \(1/(r_b + r_c)\) adjusted for the presence of intercepted snow or water, \(X_E\) (XEVAP in the code) is calculated as a weighted average, as follows. If \(F_s > 0\), the canopy must be at a temperature of 0 C or less, and so it is deduced that no transpiration can be occurring. Thus, \(X_E = (F_s + F_l)/r_b\) .

Otherwise, \(X_E\) is calculated on the assumption that the resistance is equal to \(r_b\) over the water-covered area and \(r_b + r_c\) over the rest: \(X_E = F_l /r_b + [1 - F_l ]/[r_b + r_c ]\)

In the 450 loop, XEVAP is first set as a trial value to RBINV (neglecting stomatal resistance), and QAC is calculated using this value. If \(q_{a, c} < q_c\) , indicating that the vapour flux is away from the canopy leaves, XEVAP is recalculated as above and QAC is re-evaluated. Otherwise, if \(F_s > 0\), XEVAP is scaled by \(F_s\) , since it is assumed that snow on the canopy will be at a lower temperature than the canopy itself, and deposition via sublimation will occur preferentially onto it. \(T_{a, c}\) and \(T_{ac, v}\) are calculated as above, and the canopy-air turbulent transfer coefficients for sensible and latent heat flux, CFSENS and CFEVAP, are set to RBINV and XEVAP respectively.

If ITC = 2, the canopy parameters and turbulent transfer coefficients are calculated, as noted above, on the basis of the assumption that the vegetation canopy and the air space within it can be treated as one aggregated mass, with a single representative temperature. Thus, the resistances \(r_a\) and \(r_b\) are considered as acting in series upon the sensible and latent heat fluxes between the canopy and the overlying air:

\(Q_{H, c} = \rho_a c_p [T_{, c} - T_{a, pot, } ]/(r_a + r_b)\)

\(Q_{E, c} = L_v \rho_a [q_{, c} - q_{a, } ]/(r_a + r_b + r_c)\)

In loop 500 the term \(1/(r_a + r_b)\) is calculated from the variables CFLUX (the inverse of \(r_a\) ) and RBINV (the inverse of \(r_b\) ), and assigned to the temporary variable CFLX. If the incoming visible shortwave radiation QSWINV is greater than or equal to \(25 W m^{-2}\) , the calculated value of CFLX is retained; if QSWINV is zero, it is reset to CFLUX; and between these two limits it varies linearly between the two.

Thus the effect of the calculated leaf boundary resistance is suppressed during conditions of zero or low solar heating. (This is done to avoid unrealistically low calculated turbulent fluxes at night, which can lead to anomalously low canopy temperatures.) The overall aerodynamic resistance \(r_A = r_a + r_b\) is obtained as 1/CFLX and assigned to the variable RA. As with ITC = 1, if \(q_a < q_c\) , XEVAP is recalculated as above (except that \(r_A\) is substituted for \(r_b\) ). If \(F_s > 0\), XEVAP is again scaled by \(F_s\) .

In the approach used here, the specific humidity of the aggregated canopy \(q_{0, c}\) is not assumed to be equal to the saturated specific humidity at the canopy temperature, but is rather determined using \(X_E\) . If the two methods of calculating \(Q_{E, c}\) are assumed to be analogous:

\(Q_{E, c} = L_v \rho_a X_E [q_{, c} - q_{a, } ]\)

\(Q_{E, c} = L_v \rho_a [q_{, 0, c} - q_{a, } ]/r_A\)

then solving for \(q_{0, c}\) leads to the expression

\(q_{0, c} = r_A X_E q_c + [1 - r_A X_E ]q_a\)

In the second part of the 500 loop the saturated specific humidity of the canopy is calculated as before and adjusted using the above equation to obtain the specific humidity of the aggregated canopy. This is then used to calculate \(T_{c, v}\) . The canopy-air turbulent transfer coefficients for sensible and latent heat flux, CFSENS and CFEVAP, are both set to CFLX, and for calculation purposes in the following loop \(T_{a, c}\) is set to the potential temperature of the air above the canopy, and \(q_{a, c}\) to the specific humidity of the air above the canopy.

In the 525 loop, the terms of the canopy energy balance are evaluated. The energy balance is expressed as:

\(K_{*c} + L_{*c} - Q_{H, c} + Q_{H, g+} - Q_{E, c} - \Delta Q_{S, c} = 0\)

(The term \(Q_{H, g+}\) corresponds to the variable QSGADD discussed above; the term \(\Delta Q_{S, c}\) represents the change of energy storage in the canopy.) The net shortwave radiation \(K_{*c}\) was evaluated earlier in the 50 loop. The net longwave radiation is obtained as:

\(L_{*c} = (1 - \chi) [L \downarrow + L \uparrow_g - 2 L \downarrow_c ]\)

\(Q_{H, c}\) is calculated as above. If there is intercepted liquid or frozen water on the canopy, or if the vapour flux is downward, or if the stomatal resistance is less than a limiting high value of \(5000 s m^{-1}\) , the canopy vapour flux \(E_c\) (that is, \(Q_{E, c} /L_v\) ) is calculated as above and the flag IEVAPC is set to 1; otherwise \(E_c\) and IEVAPC are both set to zero and \(q_c\) is set to \(q_a\) . If the water vapour flux is towards the canopy and the canopy temperature is greater than the air dew point temperature, the flux is set to zero. In the following IF block, checks are done to ensure that the calculated canopy evapotranspiration flux will not exceed the available water. If intercepted snow is present on the canopy, the maximum evapotranspiration is set to this amount. Otherwise, if the calculated evapotranspiration rate exceeds the liquid water stored on the canopy, the rate is tested to ensure that the additional transpiration does not exhaust the available soil water weighted according to the root distribution. \(Q_{E, c}\) is calculated from \(E_c\) and \(\Delta Q_{S, c}\) is obtained as

\(\Delta Q_{S, c} = C_c [T_c - T_{c, o} ]/ \Delta t\)

where \(C_c\) is the canopy heat capacity, \(T_{c, o}\) is the canopy temperature from the previous time step and \(\Delta t\) is the length of the time step. The residual RESID of the energy balance is evaluated on the basis of the current estimation for the canopy temperature TCAN. 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 ITC = 1, the calculations for the bisection method of solution in loop 550 are performed, over each of the modelled areas for which ITER = 1. If NITER = 1 (indicating that this is the first step in the iteration), then if RESID > 0 (indicating that the current value for TCAN had undershot the correct value), TCAN 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 TCAN has undershot the correct value and the last temperature increment had been a negative one) or if RESID < 0 and TSTEP > 0 (indicating that TCAN 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 TCAN. If TCAN is vanishingly close to 0 C, it is reset to that value. The iteration counter NITER and the flag NUMIT are each incremented by one. Finally, if NUMIT > 0, the iteration cycle is repeated from line 400 on.

If ITC = 2, the calculations for the Newton-Raphson method of iteration in loop 575 are performed, over each of the modelled areas for which ITER = 1. As outlined above, 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 TCAN 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 TCAN, which in turn is equal to the sum of the derivatives of the individual terms:

\(d(L_{*c} )/dT = -8 \sigma T_c^3 (1 - \chi)\)

\(d(Q_{H, c} )/dT = \rho_a c_p {1/r_A + [T_c - T_{a, pot} ] d(1/r_A)/dT}\)

\(d(Q_{E, c} )/dT = L_v \rho_a {X_E dq_c /dT + [q_c - q_a ] dX_E /dT}\)

\(d \Delta Q_{S, c} /dT = C_c / \Delta t\)

The term \(d(1/r_A)/dT\) is represented by the variable DCFLXM, which is approximated as the difference between CFLX and its value for the previous iteration, CFLUXM, divided by TSTEP. The term \(dX_E /dT\) is represented by the variable DXEVAP, which is approximated as the difference between XEVAP and its value for the previous iteration, XEVAPM, divided by TSTEP. The calculated value of TSTEP obtained from the above calculations is constrained to be between -10 and 5 K to dampen any spurious oscillations, and is then added to TCAN. If the resulting value of TCAN is vanishingly close to 0 C, it is reset to that value. At the end of the calculations the iteration counter NITER and the flag NUMIT are each incremented by one. The values of \(T_{a, c}\) , \(q_{a, c}\) and \(T_{ac, v}\) are reset to \(T_c\) , \(q_c\) and \(T_{c, v}\) respectively. Upon exiting the loop, if NUMIT > 0, the iteration cycle is repeated from line 400 on.

After the iteration has been completed, if the Newton-Raphson method has been used, a check is carried out in loop 600 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. The flags NUMIT and IEVAPC are set to zero, and a trial value of TCAN is calculated using the virtual potential temperature of the air and the canopy specific humidity. If the absolute value of RESID is \(> 100 W m^{-2}\) , TCAN is set to this trial value. The values of \(q_{0, c}\) and the components of the surface energy balance are recalculated as above, except that \(Q_{H, c}\) and \(Q_{E, c}\) are assumed as a first approximation to be zero. RESID is determined on this basis, and is then assigned to \(Q_{H, c}\) and \(Q_{E, c}\) . If RESID > 0, \(Q_{E, c}\) is set to this value; otherwise it is equally divided between \(Q_{H, c}\) and \(Q_{E, c}\) . The residual is then reset to zero, and E(0) and \(T_{v, c}\) are recalculated. NUMIT is incremented by 1, and the flag IEVAPC for the current location is set to 1.

After loop 600, calls to DRCOEF or FLXSURFZ are performed to re-evaluate the surface turbulent transfer coefficients for any locations for which the fluxes were modified in the previous loop, i.e. for any locations for which IEVAPC was set to 1. After this a check is performed for unphysical values of the canopy temperature, i.e. for values greater than 100 C or less than -100 C. If such values are encountered, an error message is printed and a call to abort is carried out.

Next a check is carried out to determine whether freezing or melting of intercepted water has occurred over the current time step. If this is the case, adjustments are required to the canopy temperature, the intercepted liquid and frozen water amounts and fractional coverages, and to \(C_c\) and \(\Delta Q_{S, c}\) . In loop 650, if there is liquid water stored on the canopy and the canopy temperature is less than 0 C, the first half of the adjustment to \(\Delta Q_{S, c}\) is performed, and the flags ITER and NIT are set to 1. The available energy sink HFREZ is calculated from CHCAP and the difference between TCAN and 0 C, and compared with HCONV, calculated as the energy sink required to freeze all of the liquid water on the canopy. If HFREZ \(\leq\) HCONV, the amount of water that can be frozen is calculated using the latent heat of melting. The fractional coverages of frozen and liquid water FSNOWC and FRAINC and their masses SNOCAN and RAICAN are adjusted accordingly, TCAN is set to 0 C, and the amount of energy involved is stored in the diagnostic variable QMELTC. Otherwise all of the intercepted liquid water is converted to frozen water, and the energy available for cooling the canopy is calculated as HCOOL = HFREZ - HCONV. This available energy is applied to decreasing the temperature of the canopy, using the specific heat of the canopy elements, and the amount of energy that was involved in the phase change is stored in the diagnostic variable QMELTC. In both cases \(q_c\) and \(T_{c, v}\) are recalculated, and at the end of the loop \(C_c\) and \(\Delta Q_{S, c}\) are re-evaluated.

In loop 675, if there is frozen water stored on the canopy and the canopy temperature is greater than 0 C, the first half of the adjustment to \(\Delta Q_{S, c}\) is performed, and the flags ITER and NIT are set to 1. The available energy for melting, HMELT, is calculated from CHCAP and the difference between TCAN and 0 C, and compared with HCONV, calculated as the energy required to melt all of the frozen water on the canopy. If HMELT \(\leq\) HCONV, the amount of frozen water that can be melted is calculated using the latent heat of melting. The fractional coverages of frozen and liquid water FSNOWC and FRAINC and their masses SNOCAN and RAICAN are adjusted accordingly, TCAN is set to 0 C, and the amount of energy involved is stored in the diagnostic variable QMELTC. Otherwise, all of the intercepted frozen water is converted to liquid water, and the energy available for warming the canopy is calculated as HWARM = HMELT - HCONV. This available energy is applied to increasing the temperature of the canopy, using the specific heats of the canopy elements, and the amount of energy that was involved in the phase change is stored in the diagnostic variable QMELTC. In both cases \(q_c\) and \(T_{c, v}\) are recalculated, and at the end of the loop \(C_c\) and \(\Delta Q_{S, c}\) are re-evaluated.

For the locations over which melting or freezing of water on the canopy has occurred, the surface fluxes must now be recalculated to reflect the modified canopy temperature and humidity. If NIT > 0, first DRCOEF or FLXSURFZ is called to re-evaluate the surface turbulent transfer coefficients over all locations where ITER has been set to 1. Loops 700 and 750 repeat the calculations done in loops 475 and 500, to obtain the surface transfer coefficients. Loop 800 repeats the calculations of the surface fluxes done in loop 525.

At the end of the subroutine, some final diagnostic and clean-up calculations are performed. If the evaporation rate from the canopy, EVAPC, is very small, it is added to the residual of the canopy energy balance equation, RESID, and then reset to zero. The overall residual is added to the sensible heat flux from the canopy. The energy balance of the surface underlying the canopy is then re-evaluated to take into account the new value of the canopy longwave radiation. If the surface temperature is close to 0 C, the residual of the equation is assigned to QMELTG, representing the energy associated with melting or freezing of water at the surface; otherwise it is assigned to the ground heat flux. If the water vapour flux is towards the canopy, the water sublimated or condensed is assigned to interception storage: to SNOCAN if SNOCAN > 0 (for consistency with the definition of CPHCHC), and to RAICAN otherwise. The diagnostic water vapour flux variables QFCF and QFCL, and the diagnostic change of canopy heat storage HTCC, are updated accordingly, and the canopy heat capacity CHCAP is recalculated; EVAPC is added to the diagnostic variable EVAP and then reset to zero. The net shortwave radiation, outgoing longwave radiation, sensible and latent heat fluxes are calculated for the whole canopy-ground surface ensemble; the evaporation rates are converted to \(m s^{-1}\) ; and the iteration counter ITERCT is updated for the level corresponding to the subarea type and the value of NITER. Finally, the fraction of water removed by transpiration from each soil layer is updated according to the fractions determined above on the basis of water availability.