! canopy_BBstd_19aug11_D.f90 ~/job/pecan_modelling
! Added, relative to 8aug11 model: writing extra info (gcan, etc.) to a unit 13
! Aim of the model:
! * Predict photosynthesis (A) and transpiration (E) of a tree in an orchard
! replete with other trees that affect its light interception and
! energy balance...
! * While accounting for the distribution of realistic microenvironments of leaves
! in diverse positions and orientations
! * Using process-based models of physiological and biophysical processes,
! rather than statistical models.
! * The model is intended for:
! * Simulating field measurements of managed tree crops and wild lands
! (forests), to test our understanding of the roles of physiological
! controls and canopy structure in determining A and E
! * Predictive use (e.g., irrigation management), but more often when it
! is used to parametrize more empirical models
! * Exploring the roles of biophysical and physiological controls over A
! and E, in general
!
! Features:
! * Use of real weather data, with estimation of skylight conditions
! * Light interception by leaves of varied canopy positions and orientations with:
! * Orchard geometry - in a flexible specification, using ellipsoids of
! arbitrary orientation
! * Leaf angle distribution
! * Direct and diffuse light, with accounting of solar angles
! * Resolution of PAR, NIR, and TIR fluxes - all are important for energy balance
! * Enzyme-kinetic formulation of PS,including co-limitation by Rubisco and by
! RuBP regeneration (but not resolving triose-P shuttle limitation)
! * Hence, accounting for effects of temperature (short of deactivation)
! * Also accounting to first order for the role of leaf N content, which can
! be modelled as conferring PS capacity with
! * Stomatal control with refined empirical model of conductance (gs)
! * Ball-Berry form, modified acc. to Dewar (acronym = BBD) to replace Cs in
! denominator of index with Ci, and to use (Agross+Rd) instead of Anet in
! numerator...but I'm keeping the factor hs, not using 1/(1+Ds/D0)
! ** 25 Apr. 11: I reverted, in this program version, to using Cs, because
! the fit to David Johnson's data with the BBD form was very poor
! (there may have been issue with incomplete equilibration of gs in
! his data, but, in any case, we can't use it)
! * Response to root and leaf water potential
! * Computation of water potentials from soil WP and transport resistance
! (soil, stem)
! * In turn, computation of soil WP from soil water content and soil texture
! class. The soil water content is tracked with a soil water balance model
! * Resolution of leaf temperatures using energy balance
! * Accounting for classes of temperatures acc. to PAR, NIR, TIR interceptions,
! thus, for variation according to positio
! * Accounting for air temperature, transpirational cooling
! * Simultaneous solution of PS, gs, and energy balance equations with
! efficient methods
! * Canopy aerodynamic resistances that cause offset of in-canopy air
! temperature and humidity as sensible heat and water vapor are transported out
! * Reporting of very many diagnostics of canopy performance - fluxes, energy
! balance, flux control coefficients for sto
! * Leaves' photosynthetic capacity is represented as acclimated
! to local mean irradiance
!
! Limitations:
! * Soil energy balance and evaporation not accounted currently (this can be
! added, but stability of iterations for canopy eair and Tair is reduced;
! also, a separate model of soil transport of heat and water is needed; I
! have developed some, but they're a bit tedious to implement; use Ted
! Sammis' models, in the medium term)
! * Canopy structure is static (prescribed); growth is not tracked
! * Leaves are all the same physiologically and optically, rather than acclimated
! to local mean irradiance (I could use U:lo Niinemets' model - calculate
! the mean irradiance over a day before starting iterations, and then
! scale Vcmax25 to a [linear] function of this Qsol. As it stands now, the
! leaves with higher average shade contribute too much to leaf respiration).
! * Similarly, leaf capacities do not age.
! * Transient irradiance effects on photosynthetic rates not accounted (for
! overstory plants such as we are modelling here, these are probably of
! small impact, integrated over the tree).
! * Light interception does not account for leaf clumping. Empirical corrections
! have been developed (incl. Gutschick, 1991, book chapter). Clumping alters
! the proportion of sunlit to shaded leaves and its effect on PS rates
! declines in canopies of high leaf area index.
! * Ice, snow, dew, and canopy interception of precipitation are not accounted
! (This model is meant for our semi-arid climate and the growing season, only)
! * Currently, accounting for respiration is limited to leaf respiration.
! Empirical corrections can approximate the autotrophic respiration in other
! tissues and heterotrophic respiration in soil (including decomposition of
! dead roots and litter)
! * The model is not intended to track allocation of photosynthate to growth
! and maintenance (we can patch this into Ted's model, or vice versa).
! A model that relieved all these limitations would be humongous and not necessarily
! much more accurate for the present purposes.
!
! An outline of what the main program and the modules do, and how they interact,
! follows the description of the process models (it only makes sense after
! seeing this description)
!
! A lengthy description of the model follows, in comment lines. The program
! statements begin on line 392, with type declarations for variables
!
! The model simulates the CO2 assimilation (A), transpiration (E), and
! sensible heat flux (H) from each leaf in a tree, and then is able to
! compose A, E, and from the whole tree and of the whole stand per
! unit ground area.
! To do this, we need to estimate the complete distributions of PAR, NIR,
! TIR, u, eair, and Tair around all leaves in canopy.
! We sample leaf performance at nradius x ntheta x nphi locations in a
! given central tree crown, with points centered in equal volumes of
! the whole canopy
!
! Micromet boundary conditions
! Weather data
! Basic weather data (hourly averages of Tair, RHair, u, solar shortwave
! energy flux density) are obtained from one of 2 weather stations
! on the Bosque) (there are problems in translating conditions from
! these locations to the top of canopy, above canopy b.l.; we will
! deal with these later. Basically, we have to extrapolate eair and
! Tair to a reference height above the tree canopy, using fluxes
! Estand and Hstand' from the surface beneath the weather station;
! these fluxes are poorly known)
! Data are converted to eair = partial pressure of w.v. and to SI units
! uniformly
! Data are processed further to estimate direct beam flux density in the PAR
! (I0), diffuse skylight PAR flux density (D0), and fraction of the
! time that sky is clear vs. uniformly overcast (Pclear). This
! involved method is described separately in the file XXXX.
! All micromet quantities are interpolated linearly between hourly
! records linearly, which are presumed to represent averages centered
! 0.5 h earlier than the time recorded
! Estimate of canopy boundary-layer conductance, gbcan, is
! gbcan = (rho/C1)*u, in molar units (mol m-2 s-1),
! where u = windspeed above canopy at ref. height. This relation is
! based on K-theory (with lots of adjustments to simplest theory as
! developed by P. Sellers et al. [1996] in their SiB2 model). Value of
! constant C1 depends upon mean canopy height and ref. height; we did not
! recalculate it for our geometry (esp. since the ref. height is not
! well-defined), but used values in range used in SiB2
! Solar angles (thetas, phis) are calculated for each time of each day for
! each day, using the Smithsonian formulae (we have not yet added the
! equation of time correction; it's simple)
!
! Transport processes
! Light interception: predicting (1) probability of direct light reaching
! a given canopy location; (2) same, for diffuse skylight; (3)
! distribution of total direct + diffuse irradiance in PAR, including
! accounting for varied leaf orientations
! To get penetration probabilities: ray-tracing of direct beam for
! first interceptions (diffuse skylight is modelled as composed of
! collimated beams from 5 x 5 angles theta, phi distributed
! uniformly around the sky hemisphere)
! Generalization of method of Welles and Norman (1983) for uniform
! array of spherical (or upright ellipsoidal) tree crowns
! Efficient geometric calculation for deciding if point is in a given
! crown inside an ellipsoid of arbitrary dimensions and angular
! orientation
! For lpen = penetration probability from outside entire stand to any
! location in central tree, for any solar angles: trace from loc. out to
! edge or top of canopy in 0.1 m steps, testing for being in
! canopy of any tree at each step
! Complex geometry of crowns is allowed (tilted ellipsoids; arbitrary
! locations of their centers; different sizes and angular orientations
! are allowed for each tree)
! Parametrized from field data at site
! Irradiance on single leaves is sum of diffuse flux density (D0) and
! direct flux density (I0).
! Diffuse flux arrives deterministically and is essentially same for
! all leaves of all angular orientations at a given canopy
! location (shown to be usably accurate by Gutschick, 1984)
! Direct beam arrives probabilistically, with probability
! lpen = exp(-k*pathlength)
! Thus, total flux density at a canopy location on a leaf oriented
! to give a cosine xls between it and the direct beam direction is
! D0, with probability 1 - lpen
! D0 + I0*xls, with probability lpen
! Distribution of leaf orientation angles gives rise to
! a further distribution in IL; currently, use uniform distribution
! in cos(theta) and phi, so that IL is distributed uniformly
! over the range D0 to D0+I0
! Foliage density (fd) is uniformly distributed in volume of ellipsoids, which
! have individually distinct values of the vertical axis atree,
! horiz.axis btree, and tilt angles thetac, phic
! We measured fd with LAI-2000 and corroborated stand-mean values with
! average LAI measured by Kustas and Norman
! We have to account for light (PAR, NIR) that is scattered within the crowns
! and then is (partly) absorbed by leaves within the crown. This scattered
! light is a moderate contribution (about 12% in a dense canopy) to the
! PAR flux that drives photosynthesis, and bigger part (typically nearly 30%)
! of the energy load on leaves, thus, in the determination of
! leaf temperatures, hence, photosynthetic and transpiration rates and WUE.
! The spatial pattern of scattered light and its absorption by leaves is
! in general, very difficult to calculate accurately. We use instead an
! approximation, that the crown (of the central tree) has scattering similar
! to that in a uniform layered canopy (ULC), with the same cross-sectional
! area as the projected area of the tree crown, and a depth (as leaf area
! index; call it Lavg, for "average") equal to the average in the
! ellipsoidal crown (= total leaf area / projected area).
! The solution for scattered fluxes in the ULC is readily done (I have a
! separate writeup for this; the algebra and differential equations are
! a bit too lengthy to present here). It uses similar assumptions (random
! leaf angles, unclumped leaves) as are used in the calculation of
! directly-intercepted fluxes. An additional (good) assumption is that
! the fluxes, which are angle-dependent, can be represented well with just
! upwelling and downwelling parts.
! We don't bother calculating scattered TIR fluxes, since only 4% is scattered.
! The method of estimating scattered fluxes and their absorption is developed
! thus:
! * From the consideration of interception of the direct solar beam, and
! diffuse skylight, we have the penetration probabilities to the 125 or
! so locations we resolve in the ellipsoidal canopy.
! * The sum of these, multiplied by the flux density in the direct (or
! diffuse) beam, is the total quanta (PAR) or total energy (NIR) absorbed
! by the crown directly. In a ULC,the equivalent absorption is an
! effective (top-incident) flux density, I00eff, multiplied by the fraction
! intercepted, which is 1 - exp(-Keff*Lavg)
! * Keff depends on the angle of incidence at the top of the canopy and the
! effects of neighbor canopies that may partly shade the central tree.
! We let it be set by the pattern of penetration probabilities, Ppen
! (Ppendiffuse, for diffuse skylight):
! * We rank-order the values of Ppen and fit them, as a cumulative
! distribution, to the distribution of Ppen in the equivalent ULC.
! In this ULC, Ppen is exp(-Keff*L)and we express L as Lavg*CF, with
! CF = cumulative fraction. Then, the regression of Ppen vs. CF has
! the slope Keff*Lavg. We have already set Lavg, so we now get
! Keff. There are some interesting subtleties in Keff for different
! crown shapes and spacings, all of which confirm that this regression
! is the best was to get Keff.
! * We use only the top half of the Ppen values, which represent a
! well-exposed crown and thus the real Keff.
! * We run the ULC calculation of scattered fluxes in a subroutine,
! providing it with Keff, Lavg, the number of steps in leaf area index
! (= the number of canopy locations in our ellipsoidal crown), and the
! optical properties of the leaves (fractions absorbed, scattered
! forward, and scattered backward). This gives us the amount of diffuse
! scattered flux absorbed at each depth, L, in the ULC, at an even spacing
! of depths. Call this ASR = absorbed scattered radiation.
! * Now we associate each location in our ellipsoidal crown with an equivalent
! depth in the the ULC. We use Ppen at each location and say that it is
! for L such that Ppen = exp(-Keff*L)
! * One modest complication: there are canopy locations that are at a
! depth greater than Lavg, because there are chords through the diameters
! longer than those corresponding to Lavg. We have no ULC results to
! match to these. However, we can extrapolate ASR to depths greater
! than Lavg, because the analytic formulas for the diffuse fluxes, hence,
! for ASR can be evaluated at L values beyond Lavg. They are
! well-behaved up to at least 1.4*Lavg.
! * One small note about the regression of ln(Ppen) on CF: the plot of
! ln(Ppen) vs. CF looks like a nice nearly straight line, except at very
! low Ppen, where we have cut off the search for optical path through
! the crowns when Ppen would be less than about 0.002 (we set
! a max. path length, stmax, for which Ppen achieves the value
! exp(-0.5*stmax). So, for the lowest Ppen we can get a flat line
! vs. CF. Therefore, we take care only to use ln(Ppen) > -0.5*stmax.
! * We run the estimates of ASR separately for PAR and NIR (which have quite
! different scattering amounts; NIR is rather weakly absorbed at each
! encounter). We can run direct and diffuse radiation in a single
! simulation, because the analytical solution for fluxes in the ULC is
! defined for any combination of direct and diffuse fluxes.!
! * Finally, we also distinguish an upwelling diffuse flux, that flux
! originating from radiation that hit the ground and got scattered.
! We make a few approximations here: the flux, even though it is in patches,
! is distributed as if uniform across the ground; it is uniform in angle
! (Lambertian), like diffuse skylight.
! * Clearly, we have to assign ASR from the ground-reflected radiation as
! a mirror image of that from the skylight. All we have to do is assign
! the ASR to a crown location that is at the same radius in the canopy
! and same azimuth, but the negative of the zenith angle.
! * Finally, all the diffuse flux is treated as if all leaves at a location,
! regardless of orientation, get the same amount of this flux. I proved
! that this is a reasonable approximation for diffuse skylight in a 1984
! paper in Agric. Meteorol 30: 327-341. Thus, we calculate all the ASR
! contributions and add them as a uniform amount at a given canopy location.
! That is, we we add them to the diffuse skylight contribution.
!
! Wind penetration to leaves and turbulent transport of w.v., heat
! Canopy resolved simply as one layer with uniform u, eair, Tair
! Canopy boundary layer presents a resistance 1/gbcan to transport of
! w.v., heat from canopy to free air at ref. height above canopy
! Because eair and Tair affect leaf sources of w.v. and heat, a
! self-consistent solution is iterated until errors,
! eair - (eair(0) + Estand*Pair/gbcan)
! Tair - (Tair(0) + Hstand/(Cp*gbcan)
! are acceptably small (80 Pa in eair error, 1oC in Tair error).
! In order of occurrence, we try simple iteration (Estand, Hstand averaged
! over most recent 2 intervals); 2nd iteration of same type but using
! Estand and Hstand from 1st iteration; 3rd iteration that
! extrapolates errors eair - (eair(0) + Estand*Pair/gbcan) and analog
! for Tair to zero, linearly extrapolating in eair and Tair; 4th and
! 5th iterations use numerical derivatives dE/d(eair), etc. to
! generate simultaneous linear eqns. for eair, Tair
! Transport of heat and water vapor in the canopy:
! As with well-used models such as SiB (SiB2, Sib3), the model effectively uses
! K-theory. This is accurate for sources modelled in fwd. direction to fluxes.
! We resolve only one layer in canopy; not bad, but could do better
! In the iterations, corrective steps in eair and Tair are constrained
! to maximal magnitudes (less than 2000 Pa in eair, less than 6oC in
! Tair), so that solutions will not go into unphysical regions from
! using linear approximation on nonlinear phenomena
! Thermal IR (TIR) fluxes
! Sky TIR flux is taken as a hemispherically uniform source with flux density
! QTIR=eps,sky*sigma*Tabs,sky**0.25. Here, eps,sky is the weighted average of
! the clear-sky value, epsskyclear,and the cloudy-sky value, 0.97 (really, this
! should be adjusted for the offset of cloud-base T from Tair at screen height),
! as Pclear*eps,sky,clear + 0.97*(1.-Pclear).
! In turn, the clear-sky emissivity is estimated (Norman and Campbell,
! Intro. to Envir'l Biophysics, pp. 163-4) as 1.72*(eair(kPa)/Tair,abs)**(1/7)
! I also allow the sky to be made 'hotter', to simulate the effect of
! poorer models that don't resolve the offset in effective radiative T between
! sky and surrounding veg.+soil. I read in a parameter fepsmove, between
! 0 (no adjustment) and 1 (move to full cloudiness, high sky T), and apply it
! to eps,sky,clear --> eps,sky,clear,use = (1-fepsmove)*eps,sky,clear+fepsmove
! NIR energy fluxes
! These are taken as equal to PAR energy fluxes. This approximation
! only works well near the edges of the canopy, because NIR is notably
! less attenuated than is PAR (higher-order scattering is pronounced)
!
! Leaf physiological and biophysical processes
! Photosynthesis (or A = CO2 assimilation) - important to get correctly, because
! stomatal conductance gs scales to it, and, thus, so does E (not
! linearly, but significantly). (Note: gs models that don't respond to
! A or to humidity only change E about 5-10% when they give the same
! average A over a day; evolutionary selection pressure must be responding
! to modest differences!)
! There are 2 regimes, light-limited,
! A,LL = Q0*IL,
! and light-saturated,
! A,sat -> Vcmax*(Ci-gamma)/(Ci+KCO).
! Here, Q0 is the initial quantum yield, and it responds to leaf-
! internal CO2 concentration and to gamma, the compensation point
! in CO2 (at which A -> 0), as
! Q0 = Q00*(Ci - gamma)/(Ci + 2*gamma).
! Gamma is a function only of temperature, as is KCO, and both
! are essentially universal (same magntiude for all plants). The
! T-dependences of gamma and KCO are given by de Pury and Farquhar (1997).
! In the light-saturated rate, Vcmax is the photosynthetic capacity.
! It is scaled from its value at a ref. T (say, 25oC) by another
! universal function of T. Vcmax25 is specific to a leaf, being
! higher in some species than others, higher in sun-exposed leaves,
! etc. It, and Q0, gamma, and KCO, are essentially unaffected by
! water stress until stress is severe. Q0, and later, Vcmax25,
! can be decreased by excess-light stress (photoinhibition).
! Generally, photoinhibition reduces Q0 so that A,LL matches A,sat.
! We can calculate A,ll and A,sat without much reference to light
! stress, except for transient changes in light level.
! There is a moderately sharp transition from A,LL to A,sat,
! described by the Johnson-Thornley equation,
! theta*A**2 - A*(A,LL + A,sat) + A,LL*A,sat = 0
! When theta = 1, the transition is sharp: A = min(A,LL, A,sat).
! When theta = 0, we obtain the rectangular hyperbola,
! A = A,sat*A,LL/(A,sat + A,LL)
! For most plants, theta is about 0.8.
! Aside from leaf temperature Tleaf, and IL, we must provide the value
! of Ci in order to solve for A. Ci is determined by the CO2 partial
! pressure in ambient air (Ca), by A itself, and by
! the conductance for CO2, gtotCO2, via the transport relation,
! Ci = Ca - A*Pair/gtotCO2,
! with Pair = total air pressure. To finish the calculation, we
! need the stomatal and leaf boundary-layer conductances, gs and
! gbL, respectively, to compose
! gtotCO2 = 1/ (1.6/gs + 1.37/gbL).
! The factors 1.6 and 1.37 convert from conductances for w.v. to
! conductances for CO2. Also, gbL is a function of leaf physical
! dimension, dleaf, and local windspeed, u, approx. as
! gbL = 0.264(mol m-2 s-0.5)*sqrt(u/dleaf)
! Leaf respiration
! A finding of ours in field work, and also a finding of others, is that
! dark respiration of leaves is about 8-10% of max. PS rate, when the
! leaf is at the mean daytime temperature for the past 2 weeks or so.
! We calculate this RdTmean.
! T varies during the day, so we scale Rd as RdTmean*exp(0.07*(TL-Tmean)
! for each leaf class
! Currently, we do not correct for effects of water stress on respiration.
! These appear to be weaker than effects on stomatal conductance.
! As of 8 Feb. 11 ff.: leaf photosynthetic capacity (and dark respiration)
! are scaled to mean fractional penetration of direct sun to each
! location in the canopy, over a representative day. PS capacity is
! represented by Vcmax25, the value of Vcmax at reference T = 25oC.
! The value read in as a physiological datum is that for a fully
! sun-exposed leaf; we might call it Vcmax25_0. At any given canopy
! location, at which mean fraction of direct-beam penetration over a day
! is avgPpen, I scale Vcmax25 = Vcmax25_0*[sqrt(avgPpen)*(1-bscale) + bscale],
! following modelling results of Gutschick and Wiegel (1988), with
! a typical value of bscale = 0.3. The same scaling is applied to
! the dark respiration parameter, RdTmean, at each canopy location,
! on the premise that both Vcmax25 and respiration rates are linearly
! related to leaf N content per unit area.
! Stem and root autotrophic respiration and soil heterotrophic respiration
! is not modelled. This can be estimated (as for comparing simulations to
! eddy-covariance measurements) empirically, taking the "semi-gross"
! PS rates (corrected only for leaf respiration) and debiting them by an
! empirical fraction of leaf photosynthate allocated to these functions.
!
! Now we must proceed to modelling gs.
! Stomatal control
! We have modified the Ball-Berry model, which, in its original form, is
! gs = m A hs / Cs + b
! Here, m and b are empirical constants, A is the CO2 assimilation rate,
! and hs and Cs are, respectively, the relative humidity and the CO2
! partial pressure at the leaf surface, beneath the boundary layer:
! Cs = Ca - 1.37 A/ gbL
! hs solved from two equivalent expressions for leaf E that are
! formulated in terms of w.v. partial pressures inside the leaf
! (eL = sat'd pressure at Tleaf) and at leaf surface (es):
! E = gs*(eL - es)/Pair = gb*(es - eair)/Pair
! -> hs = es/eL = (eair/eL + gs/gb)/(1 + gs/gb)
! This model fits our data very well at every site we have studied
! (BOsque, crops, desert shrubs, sunflowers, ...). It is an implicit
! equation for gs,in one sense, since A responds to gs, all else equal.
! Nonetheless, there are efficient ways to solve for gs, A, and Tleaf
! altogether, as will be shown later.
! We obtain the slope and intercept, m and b, for individual leaves from
! gas-exchange measts. We use mean values of m,b over all leaves on
! a species; later, we discuss why the variations in m, b over leaves
! has little effect.
! We modify this equation, to allow hs to appear to a power other than 1:
! gs = mBB (Agross+Rd) hs**hBB /Cs + bBB
! Several researchers, including Gutschick
! and Simonneau (2002), argued for a higher power than 1 for the
! response to humidity - thus, we add the exponent hBB, typically about
! 1-2.
! We also tried modifying the equation to a mixed form, changing the CO2
! response of stomata, from a response to mixing ratio, Cs, at the leaf
! surface (beneath the boundary layer) to a response to leaf-internal
! partial pressure. R. C. Dewar (Plant Cell Environ. 25[2002]: 1383-1398)
! argued that K. A. Mott showedthat gs responded to Ci, not Cs. However,
! in the data obtained by David Johnson (which may have equilibration
! problems), the fit to the data is much worse using Ci.
! Dewar also argued that gs responds to
! a total demand for photoreductant that is closer to being proportional
! to Agross + Rd than to Anet = Agross - Rd.
! A further modification is to represent the effect of water stress on gs,
! hence, also on PS and transpiration.
! The entire BB relation is multiplied by a factor
! facgs=exp(beta*psiroot*exp(-delta*psileaf))
! where psi = water potential. This is a form used by many reseachers,
! incl. R. C. Dewar. In earlier 2002 work, I used the beta factor only,
! with conc. of ABA in xylem sap, `a la Tardieu et al. This is
! appropriate for herbaceous plants such as maize; the delta factor is
! a better representation for trees. I, as others, applied the factor
! only to the slope, mBB, but this generates unrealistically modest
! responses to high stress, so now I apply facgs to the whole expression
! for gs.
!
! Water balance and water potentials:
! Soil water content is calculated in a "bucket" model Uniform soil of
! depth Sdepth has a volumetric water content theta. Total water
! content in the ground area allocated to one tree, Aground, is
! Aground*Sdepth*theta. The time derivative of this is the whole-tree
! transpiration rate, and the decrement in theta is d(theta)/dt
! multiplied by the 1-h simulation interval
! Soil water potential, psisoil, is computed from theta using the
! empirical equations of van Genuchten. The program allows the user
! to set the soil (texture) type and it selects the appropriate
! parameter values.
! The catenary of water flow is soil-->root-->leaf-->air
! The flow soil-->root is through a resistance, Rsoil, that depends on
! root length density, soil hydraulic conductivity, etc. The form
! derived elsewhere, is
! Rsoil = C*rhor*rr**2*ln(d/rr) / (2*mr*khPm), where
! C = a root-clumping adjustment (Tardieu's F, set as high as 5)
! rhor = dry matter density of live roots, about 1/4 that of water,
! or about 250 kg m-3
! rr = mean fine-root radius, about 1 mm
! d = mean spacing between fine roots = 1/sqrt(root length density),
! where RLD = (total root system length)/(total soil vol. occupied)
! and (total root system length) =
! (total live-root volume)/(pi*rr**2)
! and (total live-root volume) = (total root dry mass)/rhor
! and (total soil vol. occupied)
! = pi*(canopy radius)**2 * (soil depth)
! mr = total root dry mass
! khPm = soil hydraulic conductivity, in units of SI pressure
! gradient (Pa m-1). It equals khv/g, where khv is the
! conductivity in head units (m per m) for the gradient and
! velocity units (m s-1) for the flow; g is the acceleration of
! gravity.
! Given Rsoil and the whole-tree transpiration rate in L/h as EtreeLh,
! the program calculates psiroot = psisoil - Rsoil*EtreeLh, The units
! of Rsoil are chosen to give the decrement in MPa.
! The drop in water potential from roots to leaves is taken as a root-to-leaf
! transport resistance, Rstem (stem+branches+petiole, really), multiplied
! by EtreeLh.
!
! Leaf energy balance: we use a standard formulation for steady-state Tleaf:
! 0 = QSW + QTIRin - QTIRout -QE - Qcc
! Here,
! QSW = shortwave radiant energy gain = aPAR*EPAR + aNIR*ENIR,
! with absorptivities aPAR, aNIR and corresp. energy flux densities.
! We scale EPAR as 2.2x10^5 (J mol-1)*IL and ENIR = EPAR
! QTIRin = TIR energy flux density on leaf (factor of epsilon = 0.96
! is left out, roughly accounting for second scatterings); the
! calc. of TIR flux densities is discussed earlier
! QTIRout = (epsilon) sigma Tleaf,abs^4, with sigma = Stefan-Boltzmann
! const. and Tleaf,abs = Kelvin T
! QE = latent heat-loss rate, lambda*E
! Qcc = sensible heat-loss rate = gbL Cp (Tleaf - Tair)
!
! Simultaneous solution of the 3 nonlinear/transcendental/implicit eqns.
! to yield gs, A, E, H: nontrivial, since they all involve each other,
! as well as being nonlinear, etc. Here is a method I have perfected:
! Pick gs, from a range in which gs should lie (estimates can be made,
! and the upper and lower limits of gs search can be reset according
! to search results)
! Solve energy-balance for Tleaf (all other quantities are known, such as
! IL, u, ...)
! We also can solve for hs from this.
! Use Tleaf to calculate gamma, KCO, Vcmax in assimilation equation
! Express A as the carboxylation rate and equate it to the transport rate.
! For the simple example of A = A,sat
! A = (Ca -Ci)*gtotCO2 = Vcmax*(Ci - gamma)/(Ci + KCO),
! we can rewrite this as a quadratic in Ci, then solve for Ci and
! back-substitute to get A. For the general Johnson-Thornley eqn.,
! we need A^2, which generates a quartic in Ci that we solve by a
! binary search.
! We get Cs also.
! Compose IBBtest = gs - [ mBB A hs**hBB / Cs + bBB]. This function is zero
! when the solution is correct. It is negative when gs is too low
! and positive when gs is too high; we use this behavior to run a
! binary search over gs, to get gs within a tolerance 'tol.'
!
! Accounting
! Model reports many rates and diagnostics at each time interval:
! Etree (l h-1), Atree (mol s-1), control coefficients of gs over
! E and A (see below), eair and Tair both above canopy and in
! canopy interior, residual errors in consistency of eair and of
! Tair, mean leaf T (weighted by leaf transpiration), number of
! eair - Tair iteration in which convergence was obtained
! Many other intermediate data are available, including mean leaf
! conductances, etc.
! Model reports total water use and C gain, as integrals of E*dt and A*dt,
! and the water use as mm/d on a daily basis for quick checking.
! It also calculates, for each time interval's E and A and for int(E*dt) and
! int(A*dt), the stomatrol coefficient. This is the relative change in
! E (or A, or integrated values) divided by a relative change in stomatal
! conductance gs. Having solved for gs (and E, A, H) of all leaves at
! a given time, the model then resets gs to 1.1*gs and re-solves for
! E, A< and H, to get the control coefficients; for control over E, for
! example, this is
! CE = [dE/E] / [d(gs)/gs]
! = approx. [DELTA E/ E] / 0.1 = 10 * [DELTA E ]/E
!
! Complex geometry of light (and QTIRsky) interception vs. simplifications:
! The model with ray-tracing in arrays of ellipsoidal crowns can be
! compared with a model of layered canopy (laterally uniform). The latter
! model does all the other computations and process representations the
! same way, other than for light interception.
!
! Program prompts user for data, in specified physical units
! Simulation is for every hour from midnight of the initial day, for
! selected Julian days.
! In analyses of sensitivity of results to, say, stomatal control parameters,
! the model is used with same geometry and light-interception regimes;
! to save computation time, the light-interception probabilities are
! writeable to disk the first time and reused in subsequent runs
!
! Summary of what the MAIN program and the subroutines do, and how they communicate
! MAIN:
! Reads in input data: physiological, meteorological, and edaphic/hydrologic
! Allows iinteractive resetting of (almost) any parameter (not very useful
! in batch mode, which is the mode mandated by the need to enter the
! locations and geometries of many trees - 81, in our current data file)
! Retrieves (or, later, computes) the probability of light penetration by
! the direct beam to any of the 125 or so locations that sample the canopy.
! This computation is specific to a given day and hour.
! Also retrieves or computes the (nearly deterministic) penetration of diffuse
! skylight to each location. (Effectively indep. of hour and day)
! Computes mean leaf respiration for use in all later calculations.
! Loops over days and hours of a chosen simulation interval in Julian days
! Calculates water status and stress-response factor of stomatal conductance.
! Calculates microenvironmental conditions, from processed weather data file.
! Calls subroutine Ecalc that computes full radiation and weather
! microenvironments at each of the (125 or so) canopy locations.
! Ecalc subsequently calls for computation of PS and E and other statuses
! for various classes of environments.
! Composes whole-canopy PS and E (transpiration).
! Iterates the whole-canopy calculation for consistent values of air
! temperature and w.v. partial pressure in the canopy.
! Re-solves for all fluxes at 10% higher values of stomatal conductance, in
! order to estimate control coefficients.
! Computes results of interest and a number of diagnostics (some left
! in but commented out)
! Computes soil water status for next hour of simulation.
! Ecalc: computes whole-tree PS and E, iteratively (iterating for consistency
! between A, E, and the microenvironment in the canopy affected by A, E, and
! sensible heat flux H)
! Loops over the (125 or so) discrete locations that sample the canopy of the
! central tree of interest
! If necessary (not done before and retrieved from a file), it calculates
! the diffuse and direct sunlight penetration factors to each location,
! for this hour and date
! Accounting for both varied leaf orientations and estimated fraction of the
! simulation interval (hour) as overcast, it computes the fraction of leaf
! area at this location that falls into each of 'nbin' different ranges of
! PAR irradiance.
! For each irradiance bin, it calls BBsolve, which then computes leaf energy
! balance (hence, leaf T that affects PS enzymes and dark respir.),
! stomatal conductance, PS rate, and transpir. rate per leaf area.
! Contributions are summed up and returned.
!
! BBsolve: finds the value of stomatal conductance (gs) that satisfies the
! simultaneous equations for leaf temperature (energy balance), PS rate, and
! Ball-Berry-Dewar equation for stomatal control, given a specified
! microenvironment (and leaf and root water status)
! Performs a binary search in gs values.
! At each value of gs:
! Computes leaf T (it knows all the other conditions), thus, the thermally
! activated values of PS enzyme-kinetic parameters and of dark respir.
! Computes the value of PS rate, Aleaf, that is consistent with the
! enz.-kinetic formula for PS and the transport relation for Aleaf,
! by finding the appropriate value of leaf-internal CO2 partial pressure,
! Ci (it calls subr. CiACssolve to do this).
! Computes surface relative humidity, hs, also.
! Tests if gs meets the BB relation
! Stops iterating after a fixed number of iterations, set by user-chosen tolerance
!
! CiACssolve: finds the value of Ci that satisfies simultaneously the enzyme-kinetic
! transport equations for leaf PS rate.
! The equation includes both light-limited and light-saturated rates, combined via
! the Johnson-Thornley transition.
! Computes coefficients of Ci in resulting quartic equation.
! Solves for Ci by a hobbled Newton-Raphson search.
!
! Tleaf: solves for leaf temperature.
! It is provided with the radiant fluxes (SW and TIR), boundary-layer conductance,
! etc. The value of stomatal conductance, gs, is the only variable.
! It iterates the eqn. for static energy balance, by a Newton-Raphson method.
!
! esat: computes the partial pressure of water vapor over pure water at
! temperature T, using the empirical relation in the Smithsonian tables
!
! esatpr: computes the derivative of esat with respect to T, using the same relation
!
! julian: computes the Julian date from a string specifying the date as "mmdd"
!
! iterate: interpolates the canopy air T and w.v. partial pressure (eair) from the
! 3 recent estimates that gave the least error in T and eair (error taken as
! the shift in both, between values used to start an iteration and values
! computed from resulting E and H fluxes and canopy aero. resistance).
!
! control: hobbles the shift in canopy eair to be used in next iteration for
! canopy fluxes, so that extreme oscillations are damped
!
! Tcontrol: does same for shift in canopy air T
!
! psisoilcalc: calculates soil water potential from volumetric soil water content
! and soil texture class, using van Genuchten's parametrization.
!
! Rsoilcalc: calculates root-to-soil hydraulic resistance, from soil water potential
! and root geometry, using a steady-state model of flow from soil to effectively
! isolated cylindrical fine roots.
! Uses van Genuchten's empirical equations to relate soil hydraulic conductivity
! to soil water content, for a given soil texture type
!
! Who calls whom: the figure below is readable in a fixed-width font, only. It
! does not represent the calls added recently for calculating diffuse scattered
! fluxes of radiation, which are graphed in a subsequent figure.
!
! MAIN
! / | \
! / | \ (for processing weather data)
! / | |\
! / | | \
! / | | \|/ \
! / | \|/ julian \
! / \|/ Ecalc \
! | (for | \ \
! | RdTmean) | \ \
! | | | \(for \
! | \ | +10% gs) |
! | _\| | \ |
! | BBsolve \ |
! | | \ |
! | |\ \ \|/
! | | \ \
! | | \ \ ----------
! | | _\| \|/ |
! | | Tleaf ------->esat |
! | | \ |
! |\ | \-------->esatpr |
! | \ | /
! | \ \ /
! | \ \ /
! | \ \ |/_
! | \ \------------->CiACssolve
! \|/ \
! iterate \
! | \
! | \
! |\ |
! | \---> control <-
! | |
! | |
! \ |
! \--->Tcontrol <-
!
!
! Who calls whom for the calculation of diffuse scattered fluxes of radiation
! XXXXXXXXXXXXXXXXXXXXX TO BE DONE SHORTLY XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!
! Variable and array declarations below include units (e.g., m s-1; a "-" means unitless) and a description
parameter (numwea=10000, maxtimes=10000, maxtrees=100)
parameter (maxsamples=20)
! Type declarations for simple variables
character*2 ampm ! (-)AM vs. PM designator in hourly weather records
character*3 test ! (-) 1st 3 characters of keyboard input, for manual changing
! of any input parameters
character*8 date0 ! (mm-dd-yy) date for beginning of simulation
character*8 date0wea ! (mm-dd-yy) date at which the weather record begins
character*8 datef ! (mm-dd-yy) date for end of simulation
character*8 dt0 ! (text) literal entry from keyboard to indicate you want to change date0
character*8 dtf ! (text) literal entry from keyboard to indicate you want to change datef
character*20 chaval ! (text) combination of test+sval, for changing an input parameter
character*96 fmt ! (text) final format for results, for chosen degree of completeness
character*132 fmtt ! (text) final format for writing out title (and units) of results
character*31 form1, form2*21, form3*15, form4*28
character*48 infile ! (text) name of data-input file, if used
character*4 infix ! (text) infix to be inserted into output file names, making each run's
! filenames unique so that they don't overwrite each other
character*48 lpenfile ! (text) name of file that has previously (and tediously) calculated
! array of light-penetration probabilities to each canopy location
character*48 newlpenfile ! (text) name of new such array, if calcs. must be done anew
character*40 outfile ! (text) name of output file for computed results
character*20 scratchname ! (text) expanded name of an output file, including the infix
character*17 sval ! (various) new numerical value for any input parameter being changed
character*47 title_1,title_2*33,title_3*18,title_4*35 ! (text) successive pieces of the title
! for the results being printed out
character*47 title2_1,title2_2*33,title2_3*18,title2_4*29 ! (text)successive pieces of the
! units printout
character*134 titlefull, title2full*128 ! (text) full title and unit specifications, for chosen
! degree of completeness of results printout
character*48 weafile ! (text) name of file with weather data
integer iac ! (-) accumulator for number of days in previous month (to get to right day in file)
integer icall ! (-) icall=1 indicates the very first call to Ecalc, at which time
! it is appropriate to calculate the factors 'scaling'
integer iconv ! (-) number of the iteration of canopy eair, Tair in which convergence occurred
integer iday ! (-) for a weather record, the number of the specific day within the month
integer iflag ! (-) not used
integer ihi ! (-) intermediate calculation - 1st digit in counter for # of char. in
! a format for a literal
integer ihr ! (-) for a weather record, the hour of the day
integer ihrcorr ! (-) ihr, corrected for 12 h addition in the afternoon
integer iinfo ! (-) selector of the degree of completeness of results to be printed
integer iinfoextra ! (-) selector for even more results to be written, but only to disk
integer iline ! (-) index of the current line number being read in the weather file
integer ilo ! (-) intermediate calculation - lowest digit in counter for # of char. in
! a format for a literal
integer imid ! (-) intermediate calculation - middle digit in counter for # of char. in
! a format for a literal, when # char. is 100 or more (<200)
integer index1 ! (-) index number of the hourly weather record to retrieve when interpolating data
integer indexx ! (-) used in debugging, to pick out one sky angle for checking penetr. calc.
integer infixlen ! (-) number of nonblank characters in the infix to be inserted into names of
! output files, so that each run in unique and outputs don't overwrite each other
integer ipensave ! (-) switch: set to 1 means 'save new calc. of light-penetration probabilities;
! set to 0 means don't save
integer ipre ! (-) switch: set to 1 means that light-penetr. array has been calculated in prior run,
! and you want to use them; set to 0 means they don't exist and need to be calculated
integer iRHpct ! (-) air relative humidity, as percent, in original weather data file
integer isoiltype ! (-) identifier of soil texture type, for hydrology and waer potential
integer iter ! (-) number of the current iteration in solving for eair, Tair in canopy
integer jd,jd0,jd0wea,jdf ! (-) Julian dates of current weather record, day of start of simulation,
! day on which weather record begins, and final day of simulation
integer JDdiff ! (-) offset between 1) the day used for scaling PS capacity according to
! mean light penetration to each location and 2) the day that starts the
! weather record (must be .ge.0)
integer lentitle ! (-) number of characters in the final title used for writing out hourly resutls
integer month !(-) month of the year, in a weather record being read
integer nbins ! (-) number of discrete leaf irradiance levels to resolve
integer ncha ! (-) ID number of an input parameter being changed manually
integer nd ! (-) number of steps of size 'stepsize' to take through canopy
! while accumulating total length of path within tree crowns
integer ndmax ! (-) max. number of steps (prevents one from wasting time
! looking outside all known trees)
integer nhr ! (-) do-loop index: hour of the current day in the simulation
integer nphi,nradius,ntheta ! (-) number of discrete azimuth angles/radii/zenith angles
! to resolve in sampling locations within the canopy for light penetration
integer ntimes ! (-) number of the current simulation interval
integer ntimetot ! (-) number of the final simulation interval
integer ntrees ! (-) number of trees to be described around central tree
logical BBswitch ! (-) switch: set to 1 means we will use Ball-Berry humidity response; 0 - will not
logical wrtest ! (-) switch: if this is a multiple of 6 hours after 1st hour, write headers
real AAtree ! (mol s-1) accumulator for CO2 assimilation rate of whole central tree
real AAtree10 ! (mol s-1) same, with gs increased 10%, for sensitivity calculation
real aax ! (m) a distance scalefor sampling inside the ellipsoid representing the central tree
real Agrand ! (mol) accumulator for CO2 assimilation of central tree, integrated over whole
! simulation interval
real Agrand10 ! (mol) same, for tree with gs increased 10%
real Aground ! (m2) ground area covered per tree
real Aleaf ! (mol m-2 s-1) photosynthetic rate per leaf area for leaf under consideration
real alph ! (m) distance scale for sampling inside the central tree
real aNIR ! (-) absorptivity of leaf for NIR radiation
real aPAR ! (-) same, for PAR
real arg ! (-) temporary variable, in several locations
real avgI ! (mol m-2 s-1) average PAR irradiance over all loc. in tree
real bax ! (m) a distance scale for sampling inside central tree's ellipsoid
real bBB ! (mol m-2 s-1) intercept in Ball-Berry model of gs
real bet ! (m) intermediate in calc. of arbitary tree's projected canopy length along sun direction
real beta ! (MPa-1) coefficient of response of stomatal conductance to root water potential
real bscale ! (-) intercept in the scaling of PS capacity to mean light penetration fraction
real bterm ! (-) intermediate term in computing equation-of-time correction for solar time
real BTIR ! (W m-2K-1) derivative w/r to T of QTIR flux density from leaf
real Ca ! (Pa) partial pressure of CO2 in ambient air (at ref. height above canopy)
real CAlocal ! (-) control coefficient for gs over A of tree, for the current simulation interval
real CElocal ! (-) same, for conrol of gs over E of tree
real Ci ! (Pa) partial pressure of CO2 inside leaf
real Clump ! (-) clumping factor for fine roots, as multiplier of soil-to-root hydr. resis.
real cosalphas ! (-) cosine of solar elevation angle
real costhetaD ! (-) cosine of zenith angle of sky direction being sampled while computing
! diffuse skylight penetration
real costhetape ! (-) dot product of direction vectors, in computation of distance that
! direct solar beam travels inside tree canopy
real Cpair ! (J mol-1 K-1) heat capacity of air
real dFR ! (m) mean distance between fine roots
real Cs ! (Pa) partial pressure of CO2 at surface of leaf, beneath leaf boundary layer
real D000 ! (micromol m-2 s-1) initial estimate of adiant flux density, D00 (not radiance)
! in diffuse skylight at top of canopy
real D0PAR ! (micromol m-2 s-1) same as D00, but interpolated to this specific time
! from weather data
real Daa ! (Pa) air-to-air vapor pressure deficit above canopy
real deair ! (Pa) increment in eair within canopy for this iteration
real delta ! (-) coefficient of leaf WP amplification of water-stress
! response of stomatal conductance
real derrordeair ! (-) derivative of error in eair of canopy with respect to increment
! in eair
real dleaf ! (m) characteristic linear dimension of leaf
real dpc ! (m) distance of a point from center of selected tree
real dpen ! (-) penetration probability of direct beam to canopy location
real drow ! (m) (mean) distance between tree rows
real Ds ! radians) solar declination N of equator
real dt ! (m) distance along path from chosen point in tree to top of canopy
real dTair ! (oC or K) increment in canopy air temperature during iterative solution
real derrorTdTair ! derivative of error in Tair of canopy w/r to increment in Tair
real dtheta ! (-) change in volumetric soil water content per h interval
real dtime ! (h) cumulative (decimal) time from beginning of year
real dtree! (m) (mean) distance between columns of trees
real dummy ! (-) dummy argument to function Ecalc (copious arguments are passed in
! common blocks, instead)
real dx,dy,dz ! (m) x/y/z-components of step along direct beam path while
! checking for interception by trees
real Edaysum ! (L) Cumulative water use over current day of simulation
real Egrand10Lh ! (l h-1) transpiration rate of central tree, integrated
! over entire simulation time, after gs is increased 10% at each time
real EgrandLh ! (l h-1) same, at correct gs
real ELaavg ! (mol m-2 s-1) average transpiration rate per unit leaf area, in whole canopy
real ELaavg10 ! (mol m-2 s-1) same, after gs is increased 10%
real enlocations ! (-) number of discrete locations sampled in central tree
real EOTcorr ! (h) equation-of-time correction to civil time
real epsleaf ! (-) TIR emissivity and absorptivity of leaf (not used yet)
real epssky ! (-) effective TIR emissivity of the sky at T=Tair
real epsskyclear ! (-) true TIR emissivity of the *clear* sky at T=Tair
real epsskyclearuse ! (-) TIR emissivity of the clear sky, possibly adjusted higher to
! mimic the effect of less complete models that don't use full radiative balance
real errorenow ! (Pa) current value of error in consistency in eair in canopy
real errorTnow ! (K) same, for error in Tair
real Estandlastavg ! (mol m-2 s-1) stand tranpsiration rate, averaged
! over previous two simulation intervals
real Estandm1 ! (mol m-2 s-1) stand transpir. rate, at immediately preceding time interval
real Estandm2 ! (mol m-2 s-1) same, two intervals ago
real Etree ! (mol s-1) transpiration rate of central tree
real Etree10 ! (mol s-1) same, when gs is artificially increased 10%
real Etree10Lh ! (L h-1) same, in units of liters/h
real EtreeLh ! (L h-1) transpiration rate of central tree, in these units
real facgs ! (-) factor by which the Ball-Berry-Dewar slopr for stomatal conductance
! is cut by water stress
real Fcrown ! (-) fraction of total solar flux intercepted by central
! tree over its ground area, integrated over whole simulation time
real fd ! (m-1)= m3 m-2) foliage density inside crowns of trees
real fepsmove ! (-) a 'slider' that allows adjustment of clear-sky TIR emissivity
! toward the max. value of 0.97 for cloudy sky, to test rad. balance approxs.
real fNIR ! (-) transmissivity of leaves in the NIR, as a fraction
real fPAR ! (-) transmissivity of leaves in the PAR, as a fraction
real fraci ! intermediate in constructing zenith angles to sample
! within central tree
real gbair ! (mol m-2 s-1) conductance of leaf boundary layer, molar units
real gbcan ! (mol m-2 s-1) same, for canopy b. l.
real gcan ! (mol m-2 s-1) total conductance for w.v. of all leaves in canopy
real gtotCO2 ! (mol m-2 s-1) conductance of stomata plus leaf b. l. for CO2
real ha ! (-) notional (default) average relative humidity at leaf surface
real hBB ! (-) exponent of surface relative humidity in the response of
! stomatal conductance to the environment
real hcan ! (m) max. height of canopy (among all trees)
real hour ! (h) hour angle at simulation time, relative to local solar noon
real hs0 ! (-) mean rel. humidity at leaf surface; neg. value means
! model should use Ball-Berry model
real Hstandlastavg ! (W m-2) = sens. heat flux of stand per ground
! area, averaged over previous two time intervals
real Hstandm1 ! (W m-2) sens. heat flux of stand at previous time interval
real Hstandm2 ! (W m-2) same, 2 intervals ago
real Htree ! (W m-2) sensible heat flux from central tree, currently
! real htree ! (m) (m) - height of individual tree, when finding max. in canopy,
! hcan - has same name as transpir. rate but use is finished
! before Htree is calc'd
real Hvap ! (J mol-1) heat of vap. of water
real I00 ! (micromol m-2 s-1) radiant flux density in PAR of direct
! solar beam, normal to its direc. of propag.
real I00start ! (micromol m-2 s-1) initial est. for I00
real I0PAR ! (micromol m-2 s-1) I00, after interpolation from weather data
real Iground ! (micromol m-2) sum over all simul. intervals of
! total PAR flux density incident per ground area of stand
real ILavg ! (micromol m-2 s-1) average PAR irradiance on all leaves in
! central tree
real Isum ! (micromol m-2 s-1) = sum of PAR irradiance over all canopy locations
real Kc25 ! (Pa) Michaelis constant for binding of CO2 to Rubisco, at ref. T
! of 25oC
real Kdiff ! (-) extinction coefficient for propagation of diffuse skylight in
! reference ULC (uniform layered canopy) - used in computing absorbed
! scattered radiation deriving from scattering of solar direct and diffuse rad.
real Kdiffstore ! (-) Kdiff value, for passing in common to Ecalc
real Ko25 ! (Pa) same, for binding of O2
real kpen ! (-) average cosine of angle between leaves and sun
real kTIR ! not used now
real lat ! (deg., then converted to radians) latitude of site, N of equator
real longcorr ! (h) correction of civil time for offset of site from central
! meridian of its time zone
real longmerid ! (deg.) longitude of the central meridian for the time zone
! at the study site
real longstation ! (deg.) longitude of station (study site), with neg. sign if west
real Leafarea ! (m2) total leaf area on central tree
real mBB ! (-) slope in Ball-Berry model
real mr ! (kg) otal dry mass of fine roots on central tree
real Mcrown ! (micromol h-1) total light intercepted by all leaves in central tree
real Mground ! micromol h-1) total light intercepted by flat ground
real overrad ! (radians-1) conversion factor from deg. to radians
real Pair ! (Pa) total pressure of air
real Pclear ! (-) = fraction of time in simulation interval that sky is
! clear (estimated)
real phican ! (radians) = currently selected azimuth within central tree
real phiD ! (radians) = currently selected azimuth within sky for calc. of
! diffuse skylight penetration
real phis ! (radians) = solar azmiuth
real pi ! (-) math. constant, 3.141592...
real Ppenmin ! (-) min. value of penetration probability; used in calculating
! stmax, the longest occluded path we'll bother tracing through crowns
real PPFD ! (mol m-2 s-1) PAR irradiance on leaf
real psiroot ! (MPa) root water potential, at fine roots
real psisoil ! (MPa) soil water potential, taken as uniform over depth
real Q00 ! (mol CO2 [mol photons]-1) initial quantum yield of CO2
! assimilation, at saturating CO2
real Qrad ! (W m-2) total radiative energy gain per area of leaf
real Qsky ! (W m-2) TIR energy flux density from sky
real Qsol ! (MJ m-2 h-1) solar energy flux density
real Qveg ! (W m-2) TIR energy flux density from surface at T of vegetation
real ra ! (s m-1) - leaf b. l. resistance in old resistance units
real radmax ! (m) max. radius from location in central tree to which to
! search for other tree crowns
real RdovrA ! (-) dark respir.of leaf as fraction of saturated A at
! Tleaf = recent mean T for last 2 weeks
real rdsol ! (m) radius of canopy along solar direction
real Rdswitch ! (-) turns on or off the inclusion of 2*Rd in the
! Ball-Berry-Dewar equation for stomatal conductance
real RdTmean ! (mol m-2 s-1) leaf respir. rate at Tleaf = Tmean
real reflsoilNIR ! (-) Lambertian reflectivity of soil in the NIR, as a fraction
real reflsoilPAR ! (-) Lambertian reflectivity of soil in the PAR, as a fraction
real remain ! (-) remainder in converting decimal time to closest
! interval in weather data
real rhoair ! (mol m-3) molar density of air
real rhoovrC1 ! (mol m-3) proportionality between gbcan and windspeed
real rhor ! (kg m-3) dry mass density of live fine roots
real RLD ! (m m-3 = m-2) mean root length density
real rr ! (m) mean radius of fine roots
real Rsoil ! (m-1 s-1) soil-to-root hydraulic resistance
real Rstem ! (MPa / [L h-1] ) hydraulic resistance of stem+branches+petioles
real Sdepth ! (m) Depth of soil in orchard
real sigma ! (W m-2 K-4) Stefan-Boltzmann constant
real sinalphas ! (-) sine of solar elevation angle
real sinthetaD ! (-) = sine of solar zenith angle currently being sampled
! in calc. of diffuse skylight penetration
real st ! (m) cumulative path of sun through crowns of all trees in stand
! along current solar direction
real stepsize ! (m) increment in path search along solar beam direction
real stmax ! (m) max. value of st, at which search quits (penetration
! probability is very low)
real sumI ! (mol m-2 s-1) sum of PAR irradiance over all leaves over all
! simulation intervals
real sumpen ! (-) sum of fractions of diffuse skylight penetrating to
! canopy location
real Tdrop ! (oC) difference, Tsky,effec - Tair, from using calculated
! TIR emissivity of the sky
real theta ! current value of volumetric soil water content (-)
real thetaPS ! (-) convexity parameter, for sharpness of transition from
! light-limited to light-sat'd CO2 assimilation
real TLavg ! (oC) leaf temperature, averaged over all leaves in central tree
real Tleafval ! (oC) leaf temperature, renamed so it can be passed in common
real TLtree ! (oC) sum of (leaf T) x (leaf conductance) over all
! locations in tree
real Tmean ! oC) mean leaf T during photoperiod over last 2 weeks
! (T to which leaf respiration has acclimated)
real ToF ! (oF) Fahrenheit air T in current weather record
real tol ! (mol m-2 s-1) tolerance for error in solving for gs
real tole ! (Pa) tolerance for error in solving for eair in canopy
real tolT ! (oC or K) tolerance for error in solving for Tair in canopy
real u ! (m s-1) windspeed at screen height
real umph ! (mi h-1) = windspeed, in barbaric English units
real Vcmax25 ! (micromol m-2 s-1, immed. converted to mol m-2 s-1) photosyn.
! capacity per leaf area (light- and CO2-saturated), at 25oC
real VPDavg ! (Pa) average vapor-pressure deficit from leaf to air, over
! whole canopy
real wt ! (-) fraction of leaf area at this location in shade only
real wt1 ! (-) weighting function for previous hourly weather record, in
! interpolation of weather variables
real wt2 ! (-) = same, for next hourly record
real x ! (radians) sine of solar azimuth; also (m), x-coordinate of
! current point in central tree
real xaux ! (m) x-coordinate of center of central tree
real xnew ! (m) x-coordinate of current position in whole canopy during
! search for interception by crowns
real xo ! (m) x-coordinate of currently-sampled point in central tree
real yaux ! (m) y-coordinate of center of central tree
real ynew ! (m) y-coordinate of current position in whole canopy during
! search for interception by crowns
real yo ! (m) y-coordinate of currently-sampled point in central tree
real zaux ! (m) z-coordinate of center of central tree
real znew ! (m) z-coordinate of current position in whole canopy during
! search for interception by crowns
real zo ! m) z-coordinate of currently-sampled point in central tree
!
! Type declaration and dimensioning for arrays
character*3 ID(13) ! (-) 3-letter keyword to recognize parameter that
! the user wants to change
integer itime(15000) ! (min) time of simulation interval, measured from
! midnight on first day of weather record
integer ndm(12) ! d) number of days in given month
real ae(0:10) ! (Pa) absolute value of error in solving for eair in canopy
real aeT(0:10) ! (oC) absolute value of error in solving for Tair in canopy
real alphas(maxtimes) ! (radians) solar elevation angle for
! each simulation interval
real alphaswea(numwea)! (radians) solar elevation angle for each weather
! record
real Asave(2000) ! (mol m-2 s-1) value of leaf PS rate at all locations in tree, for
! one specific simulation time [highly specialized info; not of general utility]
real ASRlocdiffNIR(maxsamples,maxsamples,maxsamples) ! (-) absorbed scattered
! radiation from unit amount of diffuse radiation in the NIR, at canopy locations
! specified by 3 indices, for radius, zenith angle, and azimuthal angle
real ASRlocdiffPAR(maxsamples,maxsamples,maxsamples) ! (-) absorbed scattered
! radiation from unit amount of diffuse radiation in the PAR, at canopy locations
! specified by 3 indices, for radius, zenith angle, and azimuthal angle
real ASRlocdirNIR(maxsamples,maxsamples,maxsamples) ! (-) absorbed scattered
! radiation from unit amount of direct radiation in the NIR, at canopy locations
! specified by 3 indices, for radius, zenith angle, and azimuthal angle
real ASRlocdirPAR(maxsamples,maxsamples,maxsamples) ! (-) absorbed scattered
! radiation from unit amount of direct radiation in the PAR, at canopy locations
! specified by 3 indices, for radius, zenith angle, and azimuthal angle
real atree(maxtrees) ! (m) horizontal full axis of ellipsoid for given tree
real btree(maxtrees) ! (m) vertical full axis of ellipsoid for given tree
real costheta(maxsamples) ! (-) cosine of zenith angle of point
! being sampled in central tree
real D0PARwea(numwea) ! (micromol m-2 s-1) diffuse skylight flux density in
! PAR, estimated for given weather record
real Daawea(numwea) ! (Pa) air-to-air vapor pressure deficit for given
! weather record
real E(10) ! (mol m-2 s-1) unnormalized whole-tree transpiration, during iterations
real eair(0:10) ! (Pa) estimated canopy w.v. partial pressure during given iteration
real eairwea(numwea) ! (Pa) w.v. partial pressure deficit for given
! weather record
real errore(0:10) ! (Pa) error in consistency of eair with transport eqn., for
! any given iteration
real errorT(0:10) ! (K) error in consistency of Tair with transport eqn., for
! any given iteration
real Estand(0:10) ! (mol m-2 s-1) estimated transpiration per ground area, for
! any given iteration
real Hstand(0:10) ! (W m-2) estimated sensible heat flux density per ground area, for
! any given iteration
real I0PARwea(numwea) ! (micromol m-2 s-1) direct-beam flux density in
! PAR, estimated for given weather record
real lpenshort(maxsamples) ! (-) single records of penetration probability to canopy locations
real lpenstore(maxsamples,maxsamples,maxsamples) !
! (-) = penetration probability for direct beam to all canopy locations, for current simulation interval
real PARwea(numwea) ! (micromol m-2 s-1) total PAR irradiance on flat
! ground area, for each weather record
real Pclearwea(numwea) ! (-) estimated frac. of time that sky is clear
! (direct beam is present), for given weather record
real phi(maxsamples) ! (radians) sky azimuths being sampled for diffuse skylight
real phic(maxtrees) ! (radians) = azimuthal inclination of ellipsoid of each tree
real phiswea(numwea) ! (radians) solar azimuth for each weather record
real Ppendiffuse(maxsamples,maxsamples,maxsamples) !
! (-) total penetration fraction of diffuse skylight, for each location in central tree
real r(maxsamples,maxsamples) ! (m) array of radii being sampled
! in central tree, for each choice of theta and phi
real scaling(maxsamples,maxsamples,maxsamples) ! (-) Downscaling factor for
! photosynthetic and respiratory capacity at canopy location, tracking
! mean exposure to light
real serrorboth(10) ! (-) composite error in eair and Tair of canopy, for
! this iteration
real sintheta(maxsamples) ! (-) sine of zenith angle of point being
! sampled in central tree
real Tair(0:10) ! (oC) estimated canopy air T during this iteration
real Tairwea(numwea) ! (oC) air temperature above canopy,
! for each weather record
real thetac(maxtrees) ! (radians) zenith inclination of ellipsoid of each tree
real timerel(maxtimes) ! (h) time of current simul. interval, after start
! of simulation
real uwea(numwea) ! (m s-1) windspeed for each weather record
real xtree(maxtrees),ytree(maxtrees),ztree(maxtrees)
! x/y/z-coordinates of center of each tree
! Type declarations for do-loop indices
integer i,i2,iphi,iradius,itheta,itree,j,nt
! Type declarations for functions
! real control, Ecalc, esat, esatpr, Tcontrol, Tleaf
!!!
! At beginning, declare storage arrays for histograms:
real histoILnoscatt(120), histoILall(120), histoTLraw(120)
real histoAraw(120), histoEraw(120), histogtotraw(120)
real histoTL_Awtd(120), histoTL_Ewtd(120)
real histogtot_Awtd(120), histogtot_Ewtd(120)
real histoIL_Awtd(120), histoIL_Ewtd(120)
logical histotime
! Later, I might want histograms of TIR flux density
!!!
!
! common blocks
common/main2Tleaf/sigma,Hvap,BTIR
common/main2BBsolve/mBB,Rdswitch,hBB,bBB,BBswitch,ha,hs0,tol,facgs
common/main2Cisolve/Kc25,Ko25,Vcmax25,Ca,Q00,thetaPS,&
& RdTmean,Tmean,iflag
common/main2all/Pair,Tairuse,Tleafval,eairuse,gbair,Aleaf,Ci,Cs,&
& gtotw,gtotCO2, PPFD, Qrad,CioverCa,Rd
common/main2Ecalc/nbins,ipre,ndmax,iter,ipensave,icall,JDdiff,& ! integers
& bscale,aNIR,aPAR,Qsky,Etree,Etree10,AAtree,AAtree10,gtottree,& ! real variables
& TLtree,sumI,Htree,Ppenmin,fPAR,fNIR,reflsoilPAR,reflsoilNIR,Kdiffstore,&
& enlocations,Citree,hstree,&
& lpenstore,scaling,& ! real arrays
& histotime,& ! logical variable
& infix,infixlen ! info about infix to be inserted into output file names
common/main2Rsoil/Clump,rhor,rr,dFR,mr
common/special/ntimes
common/centertreedescrip/nradius,ntheta,nphi,costheta,sintheta,&
& r,phi,xo,yo,zo
common/alltrees/ntrees,atree,btree,thetac,phic,xtree,ytree,ztree,&
& hcan,kpen,fd
common/other/sinalphas,cosalphas,phis,ra,rhoair,D0PAR,I0PAR,&
& Pclear,wt,stepsize,stmax
! Later commons among subroutines only
! common/BBsolve2Tleaf/eL
! common/BBsolve2Cisolve/gtotCO2
!
data ID/'vcm','tps','rda','tme','apa','mbb','hbb','bbb','hs0','tol',&
& 'ca_','pai','dle'/ ! Changed 11 Oct.
data ndm/31,28,31,30,31,30,31,31,30,31,30,31/
title_1='Day Hr Etree Atree SWC psisoil psiroot'
title_2=' Ctrl,E Ctrl,A Rsoil facgs'
title_3=' eairf Tairf u '
title_4=' eairi Tairi avgI erre errT iter'
title2_1=' (L/h) (kg/h) (MPa) (MPa) '
title2_2=' (m-1 s-1) '
title2_3=' (Pa) (oC) (m/s)'
title2_4=' (Pa) (oC) (umol) (Pa) (oC)'
form1='(i3,i4,f8.1,f8.3,f7.3,f8.3,f8.3'
form2=',f8.3,f8.3,e10.2,f7.2'
form3=',f7.0,f6.1,f6.1'
form4=',f7.0,f6.1,f5.0,f6.0,f6.1,i4'
!
!
write(6,'("Enter name of input data file")')
read(5,*)infile
open(9,file=infile)
write(6,'("Enter infix for labelling output files (4 char. limit)")')
write(6,'("Remember to put the same infix into the name of the results file, in the input file")')
read(5,*)infix
infixlen=len_trim(infix)
!!!
read(9,*)ntimeshisto
!!!
open(14,file='wea.txt') ! This is redundant, really - simply conversions to new units'
!
scratchname='daily_fluxes'//infix(1:infixlen)//'.txt'
open(16,file=scratchname)
! open(16,file='daily_fluxes.txt')
!
write(6,'("Enter choice of information to print out each hour:")')
write(6,'(" 1=Onl\y state variables (day, hour, E, A, SWC, psi values)")')
write(6,'(" 2=Add control information (E, A control coeffs, Rsoil, gs stress factor)")')
write(6,'(" 3=Add local met drivers (ecan, Tcan, u)")')
write(6,'(" 4=Add diagnostics (initial ecan, Tcan, avg. leaf I, e&T errors, interations)")')
read(9,*)iinfo
! Compose the titles, units notes, format, and list of variables being printed
if(iinfo.eq.1)then
titlefull=title_1
title2full=title2_1
fmt=form1//')'
lentitle=47
elseif(iinfo.eq.2)then
titlefull=title_1//title_2
title2full=title2_1//title2_2
fmt=form1//form2//')'
lentitle=80
elseif(iinfo.eq.3)then
titlefull=title_1//title_2//title_3
title2full=title2_1//title2_2//title2_3
fmt= form1//form2//form3//')'
lentitle=98
elseif(iinfo.eq.4)then
titlefull= title_1//title_2//title_3//title_4
title2full=title2_1//title2_2//title2_3//title2_4
fmt= form1//form2//form3//form4//')'
lentitle=131
else
write(6,*)"Error; 'iinfo' not in [1,4]"
stop 30 ! Or put in goto XX to re-enter, but if one is using an input file
! in batch mode, there is no re-setting that's practical
endif
! Compose a format appropriate to the length of the titles
if(lentitle.gt.99)then ! Less than 200, of course
imid=(lentitle-100)/10
ilo=lentitle-100-10*imid
fmtt='(a1'//char(imid+48)//char(ilo+48)//')'
else
ihi=lentitle/10
ilo=lentitle-10*ihi
fmtt='(a'//char(ihi+48)//char(ilo+48)//')'
endif
write(6,'("Enter 0 or 1; 0=don''t write extra diagnostics (VPDavg, TLavg,&
&gcan,hsavg,Ciavg)")')
read(9,*)iinfoextra
if(iinfoextra.ne.0)then
scratchname='infoextra'//infix(1:infixlen)//'.txt'
open(13,file=scratchname)
! open(13,file='infoextra.txt')
endif
40 write(6,'("Enter nbins(=< 100), max. D00 (umol...)",&
& ", I00 (umol...)")')
read(9,*)nbins,D000,I00start
if(nbins.gt.100)then
write(6,'("Too many bins (>100); try again")')
goto 40
endif
! Some math constants
pi=3.14159265
overrad=pi/180.
! Some physical constants
sigma=5.67e-8
epsleaf=0.96 ! (This is not used yet, for TIR balance)
Hvap=4.5e4
Cpair=29.
Kc25=40.4
Ko25=24800.
! Q00=0.0732*aPAR/0.85 ! Corrects initial quantum yield from nominal
! value of 0.05 at standard Tleaf, standard Ci, and standard aPAR=0.85
! However, this has to be done after the value of aPAR is read in
! Set tolerance (convergence interval) for E = Estand, in mol m-2 s-1
! That is, set the fraction, f, of flux E that is the acceptable error
write(6,'("Enter the fraction, f, for error control")')
read(9,*)f
! tole=80. ! 80 Pa error in self-consistency of eair is tolerated
! tolT=1. ! 1oC error is tolerable
!
! First, read in weather data from file 'wtrdata.txt' that Ernesto pulled
! from weather.nmsu.edu - convert units while doing so
write(6,'("month day hr Tair eair u PAR",/,&
& " oC Pa m/s umol...")')
write(6,'("Enter latitude, longitude of station, & longitude of time-zone central",&
& " meridian, all in degrees")')
read(9,*)lat,longstation,longmerid
! Convert latitude to radians
lat=lat*overrad ! latitude of Las Cruces in radians
! Correct for the offset in longitude from the meridian for this time zone
longcorr=(longmerid-longstation)/15.
write(6,'("Enter name of file with weather data")')
read(9,*)weafile
open(8,file=weafile)
write(6,'("Enter date of first entry in weather file, as ",&
& "mm-dd-yy")')
read(9,*)date0wea
jd0wea=julian(date0wea,ndm)
do iline=1,50000
read(8,10,err=50,end=50)month,iday,ihr,ampm,ToF,iRHpct,umph,&
& Qsol
10 format(i2,1x,i2,4x,i2,a2,1x,f5.1,1x,i3,7x,f5.1,5x,f5.2)
if(Qsol.lt.0.)Qsol=0.0 ! Correction for non-working sensor recorded
! as Qsol=-1.00
Tairuse=(ToF-32.)/1.8
Tairwea(iline)=Tairuse
umph=amax1(umph,1.) ! Don't let effective u get too low (free convection
! takes over at low speeds)
uwea(iline)=umph/2.2
PARwea(iline)=625.*Qsol
! Qsol(MJ m-2/3600s)*(10^6 J/MJ)*(0.5 in PAR, rest in NIR)*(4.5 umol s-1/W)=PAR
! or Qsol*625 = PAR
esatair=esat(Tairuse)
RHfrac=float(iRHpct)/100.
eairwea(iline)=esatair*RHfrac
Daawea(iline)=esatair*(1.-RHfrac)
! Added on eve., 4 Mar. 99: use Qsol to extract estimate of I0, D0 in PAR
! First, compute solar zenith and azimuth for this date and hour; use
! the Julian day to get 1st part
iac=0
do i=1,12 ! I should replace this with retrieval from data statement
if(month.eq.i)goto 11
iac=iac+ndm(i)
enddo
11 jd=iac+iday
Ds=23.5*overrad*cos(2.*pi*(jd-172.)/365.)
! Calculate the equation of time as part of correction from civil time to solar time
bterm=0.01726*float(jd-81)
EOTcorr=0.1645*sin(2.*bterm)-0.1255*cos(bterm)-0.025*sin(bterm)
! Time in Ted's weather records i MST (not MDT), approx. local solar time
if(ampm.eq."AM")then
if(ihr.ne.12)then
ihrcorr=ihr
else
ihrcorr=12 ! Ted uses 12 AM as noon!
endif
else
if(ihr.ne.12)then
ihrcorr=ihr+12
else
ihrcorr=24 !Should be 25,but does this cause problems
endif
endif
! Take mean time of record as 0.5 h before time assigned (that is, a
! record listed as 4 PM is really for 3-4 PM, for a mean time of 3:30 PM)
hour=(float(ihrcorr)-12.-0.5-0.5+EOTcorr)*pi/12. ! This is the hour angle
! relative to local solar noon. The -0.5 term is
! to account for the mean time in a weather record being
! 1/2 h earlier than the time of final recording
! Calculate solar angles at this time and day, using Smithsonian algorithm
! Calculate sine of solar elevation (radians)
sinalphas=sin(lat)*sin(Ds) + cos(lat)*cos(Ds)*cos(hour)
! Calculate solar ELEVATION angle (Norman and Welles use this, not zenith angle)
alphaswea(iline)=atan(sinalphas/sqrt(1.-sinalphas*sinalphas))
! Calculate sine of solar azimuth, then the azimuth itself
! x=-cos(Ds)*sin(hour)/cos(alphaswea(iline))
! phiswea(iline)=-atan(x/sqrt(1-x*x)) ! convention is meauring CW from South
phiswea(iline)=hour
! Change to std. convention (CCW from E)
if((phiswea(iline).gt.0.5*pi).and.(phiswea(iline).lt.pi))then
phiswea(iline)=1.5*pi-phiswea(iline)
else
phiswea(iline)=-0.5*pi-phiswea(iline)
endif
! phi(solar) to times of PD measts.
! Make simple estimate of I0PAR, D0PAR
if(alphaswea(iline).gt.0.)then
airmass=1./sin(alphaswea(iline))
D00=200.*airmass**(-0.3) ! I can't recall the more exact formula
I00=I00start
Pclearwea(iline)=(PARwea(iline)-D00)/(I00*sinalphas)
if(Pclearwea(iline).gt.1.)then
Pclearwea(iline)=1.
I00=(PARwea(iline)-D00)/sinalphas
! Do a simple correction for sinalphas not being averaged, as it is in
! real wea. station data - e.g., if sinalphas is near 0 at mean time of
! weather record, then mean solar SW input is modest...but if I use
! sinalphas at the mean time of the record, I figure that I00 in I=I00*sinalphas
! has to be huge. I could average sinalphas, but the effort adds little
! accuracy
I00=amin1(I00,2000.)
endif
if(Pclearwea(iline).lt.0.)then
Pclearwea(iline)=0.
D00=PARwea(iline)
I00=0.
endif
else ! sun is below horizon
I00=0.
D00=0.
endif
I0PARwea(iline)=I00
D0PARwea(iline)=D00
write(14,'(i3,f9.1,i5,f9.3,i4,6f9.3,2f9.1,f9.3)')&
& iline,PARwea(iline),jd,Ds,ihrcorr,hour,sinalphas,alphaswea(iline),&
& x,phiswea(iline),airmass,D00,I00,Pclearwea(iline) ! again, the redundant weather data
enddo
50 close(8)
close(14)
scratchname='EH.surface'//infix(1:infixlen)//'.txt'
open(14,file=scratchname)
! open(14,file='EH.surface')
write(6,'("PAR,I0PAR,D0PAR for ilines=1-24:",(3f10.1))')&
& (PARwea(i),I0PARwea(i),D0PARwea(i),i=1,24)
! Record number of weather records
iwea=iline-1
open(8,file='scratch') ! For writing and rereading text to
! convert from characters to numbers
!
!
! To sample locations in canopy, we will select points at various radii and
! zenith, azimuth angle in canopy; choose how many of each:
write(6,'("Enter number of radii, zeniths, and azimuths",&
& " to use in canopy")')
read(9,*)nradius,ntheta,nphi
! The next few sections allow one flexibility in entering parameters.
! The first round, one must type in all parameters. Subsequently, one
! can enter a 3-letter keyword and a value, and only that parameter
! is changed in value. (One has the option of specifying 'all' to
! re-enter everything.) Entering a null keyword (just a carriage return)
! asks for the computation to proceed with current values.
!
! Ask if the light penetration array for all locations in the tree have
! already been calculated in a previous run
write(6,'("Is the array lpen(ntimes,iradius,itheta,iphi) ",&
& "available as a file already?")')
write(6,'(" (and Ppendiffuse, too) Yes=1, no-0")')
read(9,*)ipre
if(ipre.eq.1)then
write(6,'(" Enter name of the file with lpen already",&
& " calculated FOR THE SAME GEOMETRY")')
read(9,*)lpenfile
open(10,file=lpenfile)
write(6,'("Enter optional data for scaling PS & Rd capacities:",/,&
& " JDdiff = offset of scaling data and start of weather data;",/,&
& " a value <0 indicates no scaling will be done;",/,&
& " bscale = intercept (about 0.3)",/,&
& "Both values must be entered to avoid reading errors",/,&
& "and scaling can only be done if lpen array has already been "&
& "calculated in a prior run of the program")')
read(9,*)JDdiff,bscale
if(JDdiff.ge.0)then
! Now calculate the PS scaling factor
! First, test if the date for averaging Ppen is after the beginning of the
! weather record; if so, do dummy reads of file 10 to get to the right date
idum=(JDdiff+1)*24*nradius*ntheta
do i=1,idum
read(10,*)(lpenshort(iphi),iphi=1,nphi)
enddo
! Zero out the accumulators for average penetration fraction (use the same
! variable 'scaling' to store it)
do iradius=1,nradius
do itheta=1,ntheta
do iphi=1,nphi
scaling(iradius,itheta,iphi)=0.
enddo
enddo
enddo
! Start the accumulation over the 24 hours (we could skip the night hours,
! but why put in complicated sunrise/sunset calculations?
! Also, we need to scale relative to the max. value
scalemax=0.
do ihr=1,24
do iradius=1,nradius
do itheta=1,ntheta
read(10,*)(lpenshort(iphi),iphi=1,nphi)
do iphi=1,nphi
scaling(iradius,itheta,iphi)=scaling(iradius,itheta,iphi)+&
& lpenshort(iphi)
scalemax=amax1(scalemax,scaling(iradius,itheta,iphi))
enddo
enddo
enddo
enddo
write(6,'("Computing scaling using JDdiff=",i3)')JDdiff
! Compute the scaling functions
scalingmin=10.
scalingmax=0.
do iradius=1,nradius
do itheta=1,ntheta
do iphi=1,nphi
value=scaling(iradius,itheta,iphi)/scalemax
atemp=value*(1.-bscale)+bscale
scalingmin=amin1(scalingmin,atemp)
scalingmax=amax1(scalingmax,atemp)
scaling(iradius,itheta,iphi)=atemp
! Print out diagnostic
! write(6,'("Scaling at iradius/itheta/iphi=",3i3," = ",f6.4)')&
! & iradius,itheta,iphi,scaling(iradius,itheta,iphi)
enddo
enddo
enddo
write(13,'("Max and min scaling factors for Vcmax25 are ",2f6.2)')scalingmax,scalingmin
rewind(10) ! Else, future use is positioned at wrong record
endif ! End of block-if for calculating scaling
else
write(6,'("Program will calculate Ppendiffuse and lpen ",&
& "arrays. Do you want to save them",/,3x,"in a file?",&
& " 1=yes, 0=no")')
read(9,*)ipensave
if(ipensave.eq.1)then
write(6,'("Enter name of file in which to save Ppendiffuse",&
& " and lpen arrays")')
read(9,*)newlpenfile
open(10,file=newlpenfile)
write(6,'("Opening new light pen. file ",a48)')newlpenfile
endif
! For now, set scaling to none
do iradius=1,nradius
do itheta=1,ntheta
do iphi=1,nphi
scaling(iradius,itheta,iphi)=1.
enddo
enddo
enddo
write(6,'("Setting all scaling to 1.0, for this round")')
! Dummy read
read(9,*)JDdiff,bscale
pause
endif
!
! Set parameters that are 'invariant' for this plant
! Read in leaf parameters
110 write(6,1)
1 format(1x,"Enter aPAR,fPAR,aNIR,fNIR,reflsoilPAR,reflsoilNIR")
read(9,*)aPAR,fPAR,aNIR,fNIR,reflsoilPAR,reflsoilNIR
write(6,'(1x,"Enter Vcmax25(umol...) thetaPS Rd/A@Tmean Tmean")')
read(9,*)Vcmax25,thetaPS,RdovrA,Tmean ! Changed 11 Oct.
write(6,'("Enter mBB Rdswitch hBB bBB hs0 tol",&
& /,5x,"(hs<0 means use true BB response;", &
& " tol=precis. of gs sol (mol m-2 s-1).")')
read(9,*)mBB,Rdswitch,hBB,bBB,hs0,tol
write(6,'("Water-stress response: enter beta, delta, Rstem")')
read(9,*)beta,delta,Rstem
write(6,3)
3 format(1x,'Enter Ca(Pa), Pair (Pa) dleaf (m) fepsmove (-)')
read(9,*)Ca,Pair,dleaf,fepsmove
85 write(6,'("To modify parameters, enter the keyword and ",&
& "the value, such as ''vcm 74.5''",/,5x,"To see the ",&
& "keywords and the other commands, type ''key''")')
87 write(6,'("Enter the keyword and value (total up to ",&
& "20 characters)")')
read(9,'(a20)',err=88)chaval
goto 89 ! If there's not an error, continue the treatment
88 write(6,'("Error; try again")')
goto 87
89 test=chaval(1:3)
ncha=1
if(test.eq.'key')then
write(6,888)
888 format(1x,"Keywords are:",/,5x,"apa -> aPAR; fpa -> fPAR;",&
& "ani -> aNIR; fni -> fNIR; rsp -> reflsoilPAR;",&
& "rsn -> reflsoilNIR",/,5x,"vcm -> Vcmax25; tps-> thetaPS",&
& ";rda - > Rd/A; ; tme -> Tmean for Rd",/,5x,"apa -> aPAR; ",&
& " mbb -> mBB; rds -> Rdswitch; hbb -> hBB; bbb -> bBB; "&
& "hs0 -> hs0;",/,5x,"tol -> tol; ",&
& "ca_ -> Ca; pai -> Pair; dle -> dleaf; fep -> fepsmove",//,5x,"sto -> STOP; ",&
& "see -> print all values; all -> change all values",/,10x,&
& "(nothing, except a carriage return) -> run with current ",&
& "values",/) ! Changed 11 Oct.
goto 85
elseif(test.eq.'all')then
ncha=0
goto 110
elseif(test.eq.'sto')then
stop
elseif(test.eq.'see')then
write(6,'("aPAR=",f6.3," fPAR=",f6.3," aNIR=",f6.3," fNIR=",f6.3,&
& " reflsoilPAr=",f6.3," reflsoilNIR=",f6.3)')aPAR,fPAR,aNIR,&
& fNIR,reflsoilPAR,reflsoilNIR
write(6,25)Vcmax25,thetaPS,RdovrA,Tmean,mBB,Rdswitch,hBB,bBB,hs0,tol,Ca,&
& Pair,dleaf,fepsmove ! Changed 11 Oct.
25 format(3x,"Vcmax25=",f7.2," umol m-2 s-1; theta(PS)=",f6.3,&
& "; Rd/A=",f6.4,"; Tmean for Rd=",f6.2," oC",/,3x,"aPAR=",f6.4,&
& "; mBB=",f6.2,"; Rdswitch=",f4.1,"; hBB=",f5.1,"; bBB=",f7.4,&
& " mol m-2 s-1;",/,3x,&
& "hs0=",f7.3,"; tol=",f6.4," mol m-2 s-1; Ca=",f7.3,&
& " Pa; Pair=",f8.1," Pa; d(leaf)=",f6.4," m; fepsmove=",f6.3,/) ! Changed 11 Oct.
goto 85
elseif(test.eq.'')then
goto 90
else ! Change a specific parameter
sval=chaval(4:)
rewind 8
write(8,*)sval
rewind 8
read(8,*,err=88)val
if(test.eq.'apa')then
aPAR=val
elseif(test.eq.'fpa')then
fPAR=val
elseif(test.eq.'ani')then
aNIR=val
elseif(test.eq.'fni')then
fNIR=val
elseif(test.eq.'rsp')then
reflsoilPAR=val
elseif(test.eq.'rsn')then
reflsoilNIR=val
elseif(test.eq.'vcm')then
Vcmax25=val
elseif(test.eq.'tps')then ! Changed 11 Oct.
thetaPS=val
elseif(test.eq.'rda')then
RdovrA=val
elseif(test.eq.'tme')then
Tmean=val ! Down to here
elseif(test.eq.'mbb')then
mBB=val
elseif(test.eq.'rds')then
Rdswitch=val
elseif(test.eq.'hbb')then
hBB=val
elseif(test.eq.'bbb')then
bBB=val
elseif(test.eq.'hs0')then
hs0=val
elseif(test.eq.'tol')then
tol=val
elseif(test.eq.'ca_')then
Ca=val
elseif(test.eq.'pai')then
Pair=val
elseif(test.eq.'dle')then
dleaf=val
elseif(test.eq.'fep')then
fepsmove=val
else
write(6,'("Keyword unrecognizable; try again")')
endif
goto 85 ! Allow further changes
endif
! Convert to SI
90 Vcmax25=1.e-6*Vcmax25
BBswitch=.false.
if(hs0.lt.0.) BBswitch=.true.
Q00=0.0732*aPAR/0.85 ! Corrects initial quantum yield from nominal
! value of 0.05 at standard Tleaf, standard Ci, and standard aPAR=0.85
! Calculate Rd at Tmean; 1st, calc. A at Tmean, standard Ca, 50% RH,
! PPFD=1000 umol...,Tveg=Tmean, Tsky=Tmean-12.4 (calculated from eps,sky
! at 30oC and 50% RH --> eair=2.123 kPa), and u=1 m s-1:
Tairuse=Tmean
ha=0.5
eairuse=esat(Tmean)*0.5
gbair=0.264*sqrt(1./dleaf)
PPFD=1.e-3
Qsky=sigma*(Tairuse+273.2-12.4)**4
Qveg=sigma*(Tairuse+273.2)**4
QSW=PPFD*2.2e5*(aPAR+aNIR)
Qrad=0.5*Qsky+1.5*Qveg+QSW ! Mean PS-effective leaf sees 0.5 hemisphere
! sky, 1.5 hemisphere veg.
! NOTE: Qsky is coupled to Tair(0), not Tairuse, but the effect is minor
facgs=1. ! Initialize this - compute Rd for case of no water stress
! Set scaling factor for Vcmax25 to unity, for this fully sunlit leaf
dummy=1.
call BBsolve(dummy) ! This returns A through common
RdTmean=RdovrA*Aleaf
! write(6,'("DIAGNOSTIC: RdTmean=",f13.8)')RdTmean
! Enter descriptors of the canopy
210 write(6,'("Enter foliage density (m2 m-3)")')
read(9,*)fd
! fd=1.3 ! m2 m-3 , as measured by Jesus Arreola; 0.82, by Ernesto, later
kpen=0.5 ! for random foliage angles (assumption)
Ppenmin=0.002
stmax=log(Ppenmin)/(-kpen*fd) ! st for a ligth penetration prob=0.002
write(6,'("Enter rho/C1 for canopy boundary-layer resistance")')
read(9,*)rhoovrC1
write(6,'("Enter number of trees ")')
read(9,*)ntrees
do itree=1,ntrees
! write(6,'("Coordinates of canopy center, x(m),y(m),z(m): ")')
read(9,*)xtree(itree),ytree(itree),ztree(itree)
! write(6,'("Vertical and horizontal axes of canopy in (m): ")')
read(9,*)atree(itree),btree(itree)
! write(6,'("Zenith(0-180) and azimuth(0-360) angles of major",&
! & " axis in (degrees): ")')
read(9,*)thetac(itree),phic(itree)
thetac(itree)=thetac(itree)*overrad
phic(itree)=phic(itree)*overrad
enddo
write(6,'("Finished reading locations and geometries of ",i3," trees")')ntrees
! Read in step size that will be used to sample the light path
write(6,'("Enter stepsize")')
read(9,*)stepsize
! Read in the maximum distance (radmax) between any point in the central
! tree crown and any point in any other tree, so that we can limit the
! search path in calculations of Ppendiffuse or lpen
write(6,'("Enter radmax = max. distance (m) between any point",&
& " in central tree",/,3x,"and any point in most distant tree")')
read(9,*)radmax
ndmax=ifix(radmax/stepsize)+1
! Compute relative coordinates with respect to the center of the main
! tree canopy
xaux=xtree(1)
yaux=ytree(1)
zaux=ztree(1)
do i=1,ntrees
xtree(i)=xtree(i)-xaux
ytree(i)=ytree(i)-yaux
ztree(i)=ztree(i)-zaux
enddo
! Calculate maximum height of canopy
hcan=0.0
do i=1,ntrees
htree=ztree(i)+0.5*cos(thetac(i))*atree(i)
if(htree.gt.hcan) hcan=htree
enddo
185 write(6,'("To modify canopy parameters, enter the keyword and ",&
& "the value, such as ''fd 0.82''",/,5x,"To see the ",&
& "keywords and the other commands, type ''key''")')
187 write(6,'("Enter the keyword and value (total up to ",&
& "20 characters)")')
read(9,'(a20)',err=188)chaval
goto 189 ! If there's not an error, continue the treatment
188 write(6,'("Error; try again")')
goto 187
189 test=chaval(1:3)
ncha=1
if(test.eq.'key')then
write(6,988)
988 format(1x,"Keywords are:",/,5x,"fd_ -> fd;",&
& /,5x,"nra -> nradius; nth -> ntheta;",&
& " nph -> nphi; ",//,5x,"sto -> STOP; see -> ",&
& "print all values; all -> change all values",/,5x,"phy ->",&
& " change physiology params.; tim -> change time params.",/,10x,&
& "(nothing, except a carriage return) -> run with current ",&
& "values",/)
goto 185
elseif(test.eq.'all')then
ncha=0
goto 210
elseif(test.eq.'sto')then
stop
elseif(test.eq.'see')then
write(6,125)fd,nradius,ntheta,nphi
125 format(3x,"fd=",f5.3," m2 m-3",/,5x,"(nradius,ntheta,nphi)=(",&
& 3i5,")")
goto 185
elseif(test.eq.'')then
goto 190
else ! Change a specific parameter
sval=chaval(4:)
rewind 8
write(8,*)sval
rewind 8
read(8,*,err=88)val
if(test.eq.'fd_')then
fd=val
elseif(test.eq.'nra')then
nradius=ifix(val)
elseif(test.eq.'nth')then
ntheta=ifix(val)
elseif(test.eq.'nph')then
nphi=ifix(val)
else
write(6,'("Keyword unrecognizable; try again")')
endif
goto 185 ! Allow further changes
endif
! Select values of theta=zenith angle equally spaced in cos(theta)
! Also select values of radii so that each 'shell' has same volume
190 do iradius=1,nradius
arg=(float(iradius)-0.5)/float(nradius)
aax=0.5*atree(1)*arg**0.3333
bax=0.5*btree(1)*arg**0.3333
alph=(aax*aax)/(bax*bax)
bet=alph-1.
do itheta=1,ntheta
fraci=(float(itheta)-0.5)/float(ntheta)
arg=1.-2.*fraci
costheta(itheta)=sqrt(alph*arg*arg/(1.+bet*arg*arg))
if(fraci.gt.0.5)costheta(itheta)=-costheta(itheta)
sintheta(itheta)=sqrt(1.-costheta(itheta)**2)
!ZZZZ if(costheta(itheta).lt.0.)sintheta(itheta)=&
!ZZZZ & -sintheta(itheta)
r(itheta,iradius)=aax/sqrt(1.-sintheta(itheta)*&
& sintheta(itheta)*(1.-alph))
! r(itheta,iradius)=sqrt(1.0/(sintheta(itheta)**2.0/bax**2.0+&
! & costheta(itheta)**2.0/aax**2.0))
enddo
enddo
! Select equal increments in azimuthal angle, for same reason
do iphi=1,nphi
phi(iphi)=pi*float(2*iphi-1)/float(nphi)
enddo
enlocations=float(nradius*ntheta*nphi)
! Read in soil data
write(6,'("Soil types 1 through 3 are: sand,loamy_sand, sandy_loam.",/,&
& " Types 4 through 6 are: loam,silt_loam, sandy_clay_loam.",/,&
& " Types 7 through 9 are clay_loam,silty_clay_loam,silty_clay.",/,&
& " Types 10 and 11 are: sandy_clay,clay.")')
write(6,'("Enter isoiltype, soil depth (m), initial vol. SWC (-)")')
read(9,*)isoiltype,Sdepth,theta
! Read in data on roots and stem hydraulic resistance
write(6,'("Enter total dry mass of fine roots on tree (kg),"/,&
& " mean DM density of fine roots (kg m-3), mean FR radius (m)",/,&
& " mean root length density (m m-3), and clumping factor (-)")')
read(9,*)mr,rhor,rr,RLD,Clump
! For the moment, skip options to reset any of the soil, root, and stem parameters;
! option is not needed in short term, and is not used well in batch mode, anyway
! Calculate mean spacing between fine roots
dFR=1./sqrt(RLD)
! Read data to set up grid of times
312 write(6,'("Enter initial(dt0) and final(dtf) dates, mm-dd-yy: ")')
read(9,'(a8,1x,a8)')date0,datef
write(6,'("Enter dtree and drow")')
read(9,*)dtree,drow
Aground=dtree*drow ! Orchard ground area covered by each tree
285 write(6,'("To modify parameters, enter the keyword and ",&
& "the value, such as ''dt0 02-15-98''",/,5x,"To see the ",&
& "keywords and the other commands, type ''key''")')
287 write(6,'("Enter the keyword and value (total up to ",&
& "20 characters)")')
read(9,'(a20)',err=288)chaval
goto 289 ! If there's not an error, continue the treatment
288 write(6,'("Error; try again")')
goto 287
289 test=chaval(1:3)
ncha=1
if(test.eq.'key')then
write(6,1088)
1088 format(1x,"Keywords are:",/,5x,"dt0 -> date0; dtf -> datef; ",&
& //,5x,"sto -> STOP; see -> ",&
& "print all values; all -> change all values",/,5x,&
& "phy -> change physiol. params.; can -> change canopy ",&
& "params.",/,10x,&
& "(nothing, except a carriage return) -> run with current ",&
& "values",/)
goto 285
elseif(test.eq.'all')then
ncha=0
goto 312
elseif(test.eq.'sto')then
stop
elseif(test.eq.'see')then
write(6,'(3x,"dt0,dtf=",2a8)')dt0,dtf
goto 285
elseif(test.eq.'')then
goto 290
else ! Change a specific parameter
sval=chaval(5:)
if(test.eq.'dt0')then
dt0=sval(1:8)
elseif(test.eq.'dtf')then
dtf=sval(1:8)
else
write(6,'("Keyword unrecognizable; try again")')
endif
goto 285 ! Allow further changes
endif
!Be sure that times are referred to same base as weather data, now
! that we have changed run over times to not have same beginning as
! weather data
!
! Set up grid of times to cover chosen dates, every 10 minutes; first, read in
! an offset for the beginning time for the HPSG's, which time may not have
! been a simple integer multiple of 10 min. from midnight
290 continue
jd0=julian(date0,ndm)
jdf=julian(datef,ndm)
ntimetot=0
do jd=jd0,jdf
do nhr=1,24 ! Loop over hours in each day
! in 10-min. increments
ntimetot=ntimetot+1
! Now compute itime(ntimetot) = minutes after beg. of weather data
itime(ntimetot)=1440*(jd-jd0wea)+60*nhr ! Cumulative minutes in simulation
! measured from start of weather data, not start of simulation
timerel(ntimetot)=24.*float(jd-jd0)+float(nhr) ! Cumulative hours in
! simulation, measured from start of simulation
enddo
enddo
!
! Loop over canopy locations, to compute diffuse skylight penetration fraction
! (independent of time of day, actual day, etc.)
do iradius=1,nradius
do itheta=1,ntheta
costhetacan=costheta(itheta)
sinthetacan=sintheta(itheta)
radius=r(itheta,iradius)
! SKIP THE CALCULATIONS, if the file lpenfile exists
if(ipre.eq.1)then
read(10,*)(Ppendiffuse(iradius,itheta,iphi),iphi=1,nphi)
else
do iphi=1,nphi
! write(6,*)'iradius,itheta,iphi=',iradius,itheta,iphi!XXXX
phican=phi(iphi)
sumpen=0.
! Calculate x,y,z for canopy location!!
! First compute coordinates of this point with respect to the center
! of the central canopy as if this were vertical
xo=radius*sinthetacan*cos(phican)
yo=radius*sinthetacan*sin(phican)
zo=radius*costhetacan
! Compute the coordinates for the real canopy by rotating the point
! write(6,*)cos(phic(1)),sin(phic(1)),cos(thetac(1)),&
! & sin(thetac(1))
x=cos(phic(1))*(xo*cos(thetac(1))+zo*sin(thetac(1)))&
& -yo*sin(phic(1))
y=sin(phic(1))*(xo*cos(thetac(1))+zo*sin(thetac(1)))&
& +yo*cos(phic(1))
z=-xo*sin(thetac(1))+zo*cos(thetac(1))
!
do i=1,5 ! loop over sky zeniths
costhetaD=0.2*float(i)-0.1 ! Even steps in costhetaD, to get
! about equal coverage of solid angles in sky
sinthetaD=sqrt(1.-costhetaD*costhetaD) ! = cos(alphaD)
! write(6,'("New sky zenith angle; costhetaD(",i2,")=",&
! & f6.4)')i,costhetaD
do j=1,5 ! loop over sky azimuths
indexx=625*(iradius-1)+125*(itheta-1)&
& +25*(iphi-1)+5*(i-1)+j
phiD=2.*pi*(2.*float(j)-5.)/10.
! write(6,'("New sky azimuth anglem phi(",i2,")=",f6.3)')&
! & j,phiD
! Calculate the light path length within the original canopy
! First calculate total distance between canopy point and top of canopy
dt = (hcan - z) / costhetaD
!
! Initializing light path length and coordinates of new point
st = 0.0
xnew=x
ynew=y
znew=z
! Compute increments in x, y, and z through light path
dx=stepsize*sinthetaD*cos(phiD)
dy=stepsize*sinthetaD*sin(phiD)
dz=stepsize*costhetaD
! Compute the number of distance increments through the light path
nd = ifix(dt/stepsize)+1
! Limit the number of steps to not exceed the number that puts us
! outside the most distant crown
nd=min0(nd,ndmax)
! Loop for distance increments
do i2=1,nd
! Compute coordinates of new point
xnew=xnew+dx
ynew=ynew+dy
znew=znew+dz
! Loop for trees in the quadrant
do nt=1,ntrees
! Compute distance from new point to the center of each tree
dpc=sqrt((xnew-xtree(nt))**2+(ynew-ytree(nt))**2+&
& (znew-ztree(nt))**2)
! Compute dot product of the direction vector of the canopy symmetry
! axis with the direction vector from the center of the canopy to
! the new point
costhetape=(sin(thetac(nt))*cos(phic(nt))*(xnew-&
& xtree(nt))+sin(thetac(nt))*sin(phic(nt))*(ynew-ytree(nt))+&
& cos(thetac(nt))*(znew-ztree(nt)))/dpc
! Compute the canopy radius for the sun path direction
bet=1.-btree(nt)*btree(nt)/(atree(nt)*atree(nt))
! It is not necessary to convert canopy axis to semiaxis in this last
! equation
rdsol=0.5*btree(nt)/sqrt(1.-bet*costhetape**2) ! 0.5 to
! convert canopy axis to canopy semiaxis
! Compare the distance from the new point to the canopy center with
! the canopy radius for the sun path direction
if(dpc.le.rdsol)then
st =st+stepsize
! write(6,'("At x,y,z =",3f8.2,", we are in tree ",&
! & "number",i3)')xnew,ynew,znew,nt
endif
! if((indexx.eq.1).and.(nt.eq.1))write(6,'("st=",f6.2)')st
enddo
! Compare st with stmax for a light penetration probability= 0.01
if (st.gt.stmax)goto 701
enddo
! We've searched all canopies that might be in the way, so we have
! the total path length for diffuse light from this direction
! (specified by thetaD and phiD) to point x,y,z in canopy centered
! at the origin. Now calculate penetration probability for
! this skylight reaching (x,y,z)...and increment the accumulator
! for penetration from all skylight directions
701 dpen = EXP(-kpen * fd * st)
! write(6,*)' i,j,hcan-z,dt,nd,st',i,j,hcan-z,dt,nd,st
sumpen = sumpen + dpen
! write (6,*)st,dpen,sumpen
enddo ! End of loop over sky azimuths
enddo ! End of loop over sky zeniths
! Divide by (#diffuse zeniths)*(# diffuse azimuths) = 25
Ppendiffuse(iradius,itheta,iphi)=sumpen/25.
enddo ! End of loop over canopy azimuths
if(ipensave.eq.1)write(10,*)&
& (Ppendiffuse(iradius,itheta,iphi),iphi=1,nphi)
endif
enddo ! End of loop over canopy zeniths
enddo ! End of loop over canopy radii
! Compute how first interceptions of diffuse skylight generate scattered flux at
! all canopy locations, using the propagation in a reference canopy (uniform layered
! canopy, ULC) as a model
atreepass=atree(1) ! Need a separate variable to pass value of a variable that is
! already in common
call calc_ASR_EC(1,Ppendiffuse,Ppenmin,ASRlocdiffPAR,aPAR,fPAR,reflsoilPAR,&
& 1.,0.,atreepass,Kdiff) ! 1st arg=1 says this is for diffuse skylight only;
! args 1., 0. are nominal intensities (flux densities) of diffuse=1, direct=0
call calc_ASR_EC(1,Ppendiffuse,Ppenmin,ASRlocdiffNIR,aNIR,fNIR,reflsoilNIR,&
& 1.,0.,atreepass,Kdiff)
Kdiffstore=Kdiff ! This is a one-shot calculation of Kdiff, which needs to be
! used in later calls when the direct-beam is also being simulated
! The penetration probabilities, Ppendiffuse, are used to set up the ULC; the values
! values ASRlocdiff(iradius,itheta,iphi) are the diffuse fluxes arriving at locations
! in the real ellipsoidal canopy (EC) per unit of incident diffuse skylight
! The first argument, 1, means use diffuse light
!
! Open file to store results
write(6,'("Enter name of output file for results (40 characters or less)")')
read(9,*)outfile
open(12,file=outfile)
! Initialize accumulators for grand totals of A, E, and Iintercepted
! Possibly return here, if we're doing a rerun with new mBB, hBB, bBB, Vcmax25
8000 Iground=0.0 !!!!
EgrandLh=0.
Agrand=0.
Egrand10Lh=0. ! Etree over whole interval, with 10% higher gs at all times
Agrand10=0.
!
Isum=0.0d0
Estandm2=0.
Hstandm2=0.
! NEW 29June 09: set up tolerance for eair iterations
tole=10. ! Pa
tolT=0.1 ! oC or K
! Once and for all, compute total leaf area on tree
Leafarea=fd*0.125*1.3333*pi*atree(1)*btree(1)**2
! Factor of 0.125 converts axes to semiaxes !YYYYYYYYY
!
!
! First, write headers on data, for the results written to disk
write(12,fmtt)titlefull
write(12,fmtt)title2full
! Start the run over the simulation interval, in hours'
! Initialize the counter of calls to Ecalc (icall=1 means very first time,
! for doing calculation of PS-scaling functions
icall=0
do jd=jd0,jdf
! Initialize the accumulator for daily total water use
Edaysum=0.
Adaysum=0.
do nhr=1,24
ntimes=ntimes+1
!!!
if(ntimes.eq.ntimeshisto)then
histotime=.true.
else
histotime=.false. ! This branch is necessary, to reset histotime
endif
!!!
! Compute the current soil water potential from the volumetric SWC
psisoil=psisoilcalc(theta,isoiltype)
! Also compute root water potential - first, calculate soil-to-root
! hydraulic resistance
Rsoil=Rsoilcalc(theta,isoiltype)
psiroot=psisoil-Rsoil*EtreeLh*1.e-6 ! units are kg m-1 s-2 *10^-6
! Note that EtreeLh is not set on first interval; accept this
! (we start at midnight, when E is very small)
! There is a more accurate form that relates psiroot to psisoil, when
! psisoil is changing "fast"; it uses derivatives d(Etree)/d(psiroot)
! and d(Rsoil)/d(psisoil); the former is very difficult to get without
! tedious iterations, and the correction should be small for 1-h intervals
! For this time of day : pick up the
! light-penetration values, if they were calculated in an
! earlier run...and then, in all cases, pick up the weather
! Calculate the factor by which stomatal conductance is cut by root and
! leaf water potentials
psileaf=psiroot-Rstem*EtreeLh ! Units of last term are MPa divided
! by whole-tree transpiration in L/h
facgs=exp(beta*psiroot+delta*psileaf)
! data and compute Daa, ra, kTIR
if(ipre.eq.1)then
do iradius=1,nradius
do itheta=1,ntheta
read(10,*)(lpenstore(iradius,itheta,iphi),&
& iphi=1,nphi) ! Note that we can drop the index
! ntimes, since the arrays lpen(iradius...) have been
! stored in order
enddo
enddo
endif
! Need pick up Tair, eair, and u, for the current hour (the simulation is for
! the center of the hour interval - e.g., for hour 6 on the first day of the
! simulation, we want to pick up the 6th record in the weather file....IF the
! weather file begins on the same day as the simulation. If not, the
! simulation can only start at a date later than in the weather file, by the
! number of days jd0-jd0wea
index1=ntimes+24*(jd0-jd0wea)
! write(6,'("Using weather record # ",i5)')index1
Tair(0)=Tairwea(index1)
eair(0)=eairwea(index1)
u=uwea(index1)
Daa=Daawea(index1)/(8.314*(Tairuse+273.2)) ! in mol m-3
! Compute canopy boundary layer conductance then recompute eair as ecan
gbcan=u*rhoovrC1 ! here 0.2 is rhoovrC1, added 22 May 01
! Interpolate and store the I0PAR and D0PAR values, for later use in
! comparing predicted vs. actual light interception
I0PAR=I0PARwea(index1)
D0PAR=D0PARwea(index1)
! write(6,'(10x,"ON DAY ",i3," AT HOUR ",i2,", I0PAR=",f8.1,"; D0PAR=",f8.1)')&
! & jd,nhr,I0PAR,D0PAR
Pclear=Pclearwea(index1)
! Do same for solar angle
alphas(ntimes)=alphaswea(index1)
sinalphas=sin(alphas(ntimes))
cosalphas=cos(alphas(ntimes))
phis=phiswea(index1)
! Compute the TIR flux density from the sky, Qsky, using its relation to Tair and eair
epsskyclear=1.72*(0.001*eair(0)/(Tair(0)+273.2))**0.143
! Allow intervention to make sky T closer to air T (to test effect on radiative balance, as in simpler models)
! We replaced read-in of Tsky offset with read-in of fepsmove
epsskyclearuse=(1.-fepsmove)*epsskyclear+fepsmove
epssky=Pclear*epsskyclearuse+0.97*(1.-Pclear)
Qsky=epssky*sigma*(Tair(0)+273.2)**4
! Also report how much the effectdive Tsky is below Tair
Tdrop=(epssky**0.25-1.)*(Tair(0)+273.2)
if(ntimes.eq.ntimeshisto)write(6,'("Drop from Tair to Tsky = ",f8.2)')Tdrop
! We then pass Qsky to Ecalc in common block main2ecalc
! write(11,*)nhr,index1,alphas,sinalphas,cosalphas,phis !DDDDDDDDDDDDDDDDDD
Iground=Iground+D0PAR+I0PAR*sinalphas*Pclear
! Compute derivatives d(eair)/dT, d(QTIR)/dT
kTIR=4.*0.96*(Tairuse+273.2)**3*5.67e-8
rhoair=88600./(8.314*(Tairuse+273.3))
ra=151.*sqrt(dleaf/u)
! write(6,'("At jd=",i4," nhr=",i3,"in MAIN: u=",f6.2," & ra=",f8.2)')&
! & jd,nhr,u,ra
! Estimate eair inside canopy, including self-humidification by E (estimated
! on first iteration from E in previous two time intervals, ntimes)
Estandm1=Etree*Leafarea/(enlocations*dtree*drow) ! Use Etree rather than
Hstandm1=Htree*Leafarea/(enlocations*dtree*drow) ! Use Etree rather than
! tryin to flag whether the converged answer is E(1), E(2), E(3), or E(4)
Estandlastavg=0.5*(Estandm2+Estandm1)
Hstandlastavg=0.5*(Hstandm2+Hstandm1)
deair=Estandlastavg*Pair/gbcan
! Call for eair, to routine that constrains step in eair and min. value of eair
eair(1)=control(deair,eair(0) )
eairuse=eair(1) ! Pass this to Ecalc routine in variable 'eair'
dTair= Hstandlastavg/(gbcan*29.) ! Cpair=29 J mol-1 K-1
Tair(1)=Tcontrol(dTair,Tair(0))
Tairuse=Tair(1)
! From "crude" estimate of eair inside canopy (eair(1)), compute the
! corresponding transpiration rate of the stand, E(1)
iter=1
icall=icall+1
E(1)=Ecalc(Ppendiffuse,ASRlocdiffPAR,ASRlocdiffNIR,ASRlocdirPAR,&
& ASRlocdirNIR)
! This is actually sum(E(locations); convert now to
! average E per m2 of stand
Estand(1)=E(1)*Leafarea/(enlocations*drow*dtree)
Hstand(1)=Htree*Leafarea/(enlocations*drow*dtree)
errore(1)=eair(1)-(eair(0)+Estand(1)*Pair/gbcan) ! inconsistency measure
errorT(1)=Tair(1)-(Tair(0)+Hstand(1)/(gbcan*29.))
iconv=1 ! Will report convergence on 1st iteration if test passed
ae(1)=abs(errore(1))
aeT(1)=abs(errorT(1))
errorenow=errore(1)
errorTnow=errorT(1)
! Test if the errors in eair, Tair are less than tolerances. We would like to
! estimate tolerances as those changes in eair, Tair that shift E by a fraction
! f, say, 1%...but to estimate the sensitivities, partial E / partial ( ), we
! would have to run at least 2 more simulations at ( eair(0), Tair(1) ) and
! ( eair(1), Tair(0) ), and this would slow the program markedly. Thus, we
! will use standard tolerances, 10 Pa in eair and 0.1oC in Tair, until we have
! enough iterations (three) to compute the partial derivatives
! Estimate the composite error from both eair and Tair
serrorboth(1)=sqrt((errorenow/tole)**2+(errorTnow/tolT)**2 )
if((ae(1).lt.tole).and.(aeT(1).lt.tolT))goto 700
! Now use this Estand(1) to compute a new value of eair = eair(2), ad
! compute the corresponding E = E(2)
deair=(Estand(1)-Estandlastavg)*Pair/gbcan
eair(2)=control(deair,eair(1) )
eairuse=eair(2) ! Pass this to all subroutines
dTair=(Hstand(1)-Hstandlastavg)/(gbcan*29.)
Tair(2)=Tcontrol(dTair,Tair(1))
Tairuse=Tair(2)
iter=2
E(2)=Ecalc(Ppendiffuse,ASRlocdiffPAR,ASRlocdiffNIR,ASRlocdirPAR,&
& ASRlocdirNIR)
Estand(2)=E(2)*Leafarea/(enlocations*drow*dtree)
Hstand(2)=Htree*Leafarea/(enlocations*drow*dtree)
errore(2)=eair(2)-(eair(0)+Estand(2)*Pair/gbcan)
errorT(2)=Tair(2)-(Tair(0)+Hstand(2)/(gbcan*29.))
iconv=2
ae(2)=abs(errore(2))
aeT(2)=abs(errorT(2))
errorenow=errore(2)
errorTnow=errorT(2)
serrorboth(2)=sqrt((errorenow/tole)**2+(errorTnow/tolT)**2 )
if((ae(2).lt.tole).and.(aeT(2).lt.tolT))goto 700
! With two data points, try extrapolating eair to give zero error in eair
! ...as if all error is only from error in eair! This should be changed!!
! but it's not obvious how to do this; if we pick arbitrary steps in eair only
! for iter. 2, Tair only for iter. 3, we are almost certain not to hit the
! acceptable solution...and it seems that iter. 2 and 3 ARE often close enough
derrordeair=(errore(2)-errore(1))/(eair(2)-eair(1))
derrordeair=amax1(derrordeair,1.)
deair=-errore(2)/derrordeair
eair(3)=control(deair,eair(2) )
eairuse=eair(3)
derrorTdTair=(errorT(2)-errorT(1))/(Tair(2)-Tair(1))
dTair=-errorT(2)/derrorTdTair
Tair(3)=Tcontrol(dTair,Tair(2))
Tairuse=Tair(3)
iter=3
E(3)=Ecalc(Ppendiffuse,ASRlocdiffPAR,ASRlocdiffNIR,ASRlocdirPAR,&
& ASRlocdirNIR)
Estand(3)=E(3)*Leafarea/(enlocations*drow*dtree)
errore(3)=eair(3)-(eair(0)+Estand(3)*Pair/gbcan)
Hstand(3)=Htree*Leafarea/(enlocations*drow*dtree)
errorT(3)=Tair(3)-(Tair(0)+Hstand(3)/(gbcan*29.))
iconv=3
ae(3)=abs(errore(3))
aeT(3)=abs(errorT(3))
errorenow=errore(3)
errorTnow=errorT(3)
serrorboth(3)=sqrt((errorenow/tole)**2+(errorTnow/tolT)**2 )
if((ae(3).lt.tole).and.(aeT(3).lt.tolT))then
goto 700
else
do iconv=4,10
! Try 4th iteration, or more, up to a max. of 10
! Pass the values of eair, Tair, E, H, errore, and errorT to subroutine
! "iterate", which will find the 3 iterations with the smallest errors and
! use them to estimate the partial derivatives, dE/d(eair), dE/d(Tair),
! dH/d(eair), and dH/d(Tair). Then, we can estimate E and H at any
! eair and Tair, thus, we can solve for the values of eair and Tair that
! make eair and Tair "exactly" consistent with the transport of w.v. and heat
! through the canopy b.l. resistance
call iterate(eair,Tair,Estand,Hstand,errore,errorT,&
& serrorboth,f,tole,tolT,iconv,Pair,gbcan)
! Use the estimates that "iterate" computed for eair and Tair to compute
! the fluxes and the errors (inconsistencies) of eair and Tair with
! transport through the boundary layer
eairuse=eair(iconv)
Tairuse=Tair(iconv)
Estand(iconv)=Ecalc(Ppendiffuse,ASRlocdiffPAR,ASRlocdiffNIR,&
& ASRlocdirPAR,ASRlocdirNIR)*Leafarea/(enlocations*drow*dtree)
Hstand(iconv)=Htree*Leafarea/(enlocations*drow*dtree)
errore(iconv)=eair(iconv)-(eair(0)+Estand(iconv)*Pair/gbcan)
errorT(iconv)=Tair(iconv)-(Tair(0)+Hstand(iconv)/&
& (29.*gbcan))
errorenow=errore(iconv) !!!!
errorTnow=errorT(iconv) !!!!
! Use the tolerances, tole and tolT, computed in "iterate" to test if the
! solution (eair(iconv), Tair(iconv)) is acceptable
scalederrore=errore(iconv)/tole
scalederrorT=errorT(iconv)/tolT
! if(ntimes.eq.121)write(14,*)'iconv=',iconv,'scaled',&
! & scalederrore,scalederrorT
serrorboth(iconv)=sqrt(scalederrore**2 + scalederrorT**2)
! if(ntimes.eq.121)then
! write(14,*)'iconv',iconv,'scaled',abs(scalederrore),&
! & abs(scalederrorT)
! endif
if((abs(scalederrore).lt.1.).and.(abs(scalederrorT).lt.1.))&
& goto 700
! if(abs(scalederrore).lt.1.)then
! write(14,*)'found scalederrore <1'
! if(abs(scalederrorT).lt.1.)then
! write(14,*)'found scalederrorT <1'
! goto 700
! endif
! endif
! if(ntimes.eq.123)stop
enddo
iconv=10 ! Reverses increment to iconv=11 by do-loop ( a WUn Fortran
! quirk? !!!!
! If we reached iteration 10 and still have not converged, we should
! accept the results but print a warning
write(6,'("Still no convergence; accepting bigger",&
& " errors=",f6.1," Pa and ",f6.2," oC")')errore(10),&
& errorT(10)
write(6,'("eair(0,7,8,9,10) were ",5f8.1,/,"Estand(7,8,",&
& "9,10) ",4f10.6,/,"Tair(0,7,8,9,10) were ",5f8.1,/,&
& "Hstand(7,8,9,10) were ",4f8.2,/,"Errors in eair were",&
& 4f10.1,/,"Errors in Tair were ",4f8.2)')eair(0),eair(71),&
& eair(8),eair(9),eair(10),(Estand(i),i=7,10),Tair(0),&
& (Tair(i),i=7,10),(Hstand(i),i=7,10),(errore(i),i=7,10),&
& (errorT(i),i=7,10)
write(12,'("Still no convergence; accepting bigger",&
& " errors=",f6.1," Pa and ",f6.2," oC")')errore(10),&
& errorT(10)
write(12,'("eair(0,7,8,9,10) were ",5f8.1,/,"Estand(7,8,",&
& "9,10) ",4f10.6,/,"Tair(0,7,8,9,10) were ",5f8.1,/,&
& "Hstand(7,8,9,10) were ",4f8.2,/,"Errors in eair were",&
& 4f10.1,/,"Errors in Tair were ",4f8.2)')eair(0),eair(7),&
& eair(8),eair(9),eair(10),(Estand(i),i=7,10),Tair(0),&
& (Tair(i),i=7,10),(Hstand(i),i=7,10),(errore(i),i=7,10),&
& (errorT(i),i=7,10)
! New diagnostic for convergence: write out the E, H, eair, Tair surfaces
do iii=0,iconv
write(14,*)eair(iii),Tair(iii),Estand(iii),Hstand(iii),&
& errore(iii),errorT(iii),tole,tolT
enddo
endif
700 continue
!
! More diagnostics, if E, H are anomalously high
! if(ntimes.eq.121)then
! do iii=0,iconv
! write(14,*)eair(iii),Tair(iii),Estand(iii),Hstand(iii),&
! & errore(iii),errorT(iii),ae(iii),aeT(iii)
! enddo
! ELaavg=Etree/enlocations
! EtreeLh=ELaavg*Leafarea*64.8
! write(14,*)'As EtreeLh, this is',EtreeLh
! close(14)
! endif
! Update Isum
Isum=Isum+dble(sumI)
! Convert sum to average ELa and then to Etree in l h-1
ELaavg=Etree/enlocations
ELaavg10=Etree10/enlocations
VPDavg=Pair*Etree/gtottree
gcan=gtottree*Leafarea/enlocations
TLavg=TLtree/gtottree
hsavg=hstree/enlocations !DDDDDDDDDDDDDDDDDDD
Ciavg=Citree/gtottree !DDDDDDDDDDDDDDDDDDDD
!e write(6,199)VPDavg,TLavg,gcan
EtreeLh=ELaavg*Leafarea*64.8
! Factor of 64.8 converts mol s-1 to L h-1 (0.018*3600)
Etree10Lh=ELaavg10*Leafarea*64.8
AAtree=AAtree*Leafarea*108./enlocations ! 108 converts from
! mol CO2 s-1 to kg PSate h-1 (0.030 kg/mol * 3600 s/h)
AAtree10=AAtree10*Leafarea*108./enlocations
EgrandLh=EgrandLh+EtreeLh
Agrand=Agrand+AAtree
Egrand10Lh=Egrand10Lh+Etree10Lh
Agrand10=Agrand10+AAtree10
! Last factor converts as 0.018 l mol-1 x 3600 s h-1
! write(6,200)ntimes,ELaavg,EtreeLh,AAtree
! write(12,200)ntimes,ELaavg,EtreeLh,AAtree
200 format("For ntimes=",i4,", ELa(avg)=",f7.5," mol m-2 s-1",&
& " and Etree=",f6.2," l h-1",/,5x,"A(tree)=",f10.6,&
& " kgPSate h-1")
! Compute the local control coefficients (local to this time of the run)
CAlocal=(AAtree10/AAtree-1.)*10.
CElocal=(Etree10Lh/EtreeLh-1.)*10.
! Compute new volumetric soil water content, theta, now that we know whole-tree tpn.
! Rate of loss of water = Etree (kg s-1) = Etree (mol s-1)*(MW of water)(kg mol-1)
! = (tree ground area)(m2)* (d/dt)(Sdepth(m)*theta) * (density of water)(10^3 kg m-3)
! --> (d/dt)theta = Etree(mol s-1)*kgMW(kg mol-1)/[Sdepth(m)*Atree(m2)*rhowater(kg m-3)]
dtheta= -0.001*EtreeLh/(Sdepth*Aground) ! (-);
! EtreeLh, computed for diagnostics, is in L h-1 = kg h-1
theta=theta+dtheta
! Note: this assumes no irrigation and no rain; we're simulating time between irrigations
! Compute the irradiance, averaged over locations in the canopy
avgI=sumI/enlocations
if(mod(ntimes,24).eq.1)then
write(6,fmtt)titlefull
write(6,fmtt)title2full
endif
if(iinfo.eq.1)then
write(6,fmt)jd,nhr,EtreeLh,AAtree,theta,psisoil,psiroot
write(12,fmt)jd,nhr,EtreeLh,AAtree,theta,psisoil,psiroot
elseif(iinfo.eq.2)then
write(6,fmt)jd,nhr,EtreeLh,AAtree,theta,psisoil,psiroot,&
&CElocal,CAlocal,Rsoil,facgs
write(12,fmt)jd,nhr,EtreeLh,AAtree,theta,psisoil,psiroot,&
&CElocal,CAlocal,Rsoil,facgs
elseif(iinfo.eq.3)then
write(6,fmt)jd,nhr,EtreeLh,AAtree,theta,psisoil,psiroot,&
&CElocal,CAlocal,Rsoil,facgs,eairuse,Tairuse,u
write(12,fmt)jd,nhr,EtreeLh,AAtree,theta,psisoil,psiroot,&
&CElocal,CAlocal,Rsoil,facgs,eairuse,Tairuse,u
elseif(iinfo.eq.4)then
write(6,fmt)jd,nhr,EtreeLh,AAtree,theta,psisoil,psiroot,&
&CElocal,CAlocal,Rsoil,facgs,eairuse,Tairuse,u,&
&eair(0),Tair(0),avgI,errorenow,errorTnow,iter
write(12,fmt)jd,nhr,EtreeLh,AAtree,theta,psisoil,psiroot,&
&CElocal,CAlocal,Rsoil,facgs,eairuse,Tairuse,u,&
&eair(0),Tair(0),avgI,errorenow,errorTnow,iter
endif
if(iinfoextra.eq.1)write(13,199)VPDavg,TLavg,gcan,hsavg,Ciavg
199 format(1x,"VPDavg= ",f10.2,3x,"TLavg (gs-wtd.)= ",f10.2,3x,&
& "gcan= ",f10.2,", hsavg=",f6.3,", Ciavg=",f6.2)
! write(6,'(" *** sumI,Isum=",2f14.1)')sumI,Isum
Estandm2=Estandm1
Hstandm2=Hstandm1
Edaysum=Edaysum+EtreeLh
Adaysum=Adaysum+AAtree
enddo ! End loop over hours of the day
write(6,'(//,"For Julian day ",i3,", total water use = ",f8.2," L --> ",&
& f6.2," mm/d; PSate/area=",f8.5," kg/m^2"//)')jd,Edaysum,&
& Edaysum/Aground,Adaysum/Aground
write(16,'(i3,2x,f8.2,2x,f6.2,2x,f8.5," JD, daily Etree (L), ET (mm),"&
& " PSate/area (kg m-2)")')jd,Edaysum,Edaysum/Aground,Adaysum/Aground
enddo ! End loop over days of the simulation
! I did not indent doubly after I changed loop over ntimes into 2 nested loops,
! outer over days and inner over hours - too much work!
! Enter some information for testing vs. spherical crown program
! write(6,'("Enter dtree and drow")')
! read(9,*)dtree,drow
Mcrown=Isum*Leafarea*3600.0*1.e-6/enlocations
! Factor of 3600 converts umol s-1 umol h-1
Mground=Iground*3600.0*Aground*1.e-6
Fcrown=Mcrown/Mground
ILavg=Isum/(enlocations*float(ntimetot))
! ILavg=Isum/(enlocations*float(945-888+1))
write(6,202)Leafarea,Aground,ILavg,Mcrown,Mground,Fcrown,dtree,&
& drow
write(12,202)Leafarea,Aground,ILavg,Mcrown,Mground,Fcrown,dtree,&
& drow
202 format("Leafarea= ",f6.1,2x,"Aground= ",f5.1,2x,"ILavg= ",f7.1,2x,&
& "Mcrown= ",f10.1,2x,"Mground= ",f10.1,2x,"Fcrown= ",f6.4,2x,&
& "dtree= ",f7.2,2x,"drow= ",f4.2)
write(6,'("Total water use=",f8.0," liters; total PSate=",&
& f6.2," kg PSate; WUE=",f7.5)')EgrandLh,Agrand,Agrand/EgrandLh
write(12,'("Total water use=",f8.0," liters; total PSate=",&
& f6.2," kg PSate; WUE=",f7.5)')EgrandLh,Agrand,Agrand/EgrandLh
CAtotal=(Agrand10/Agrand -1.)*10.
CEtotal=(Egrand10Lh/EgrandLh -1.)*10.
write(6,'("Control coeff. for gs over Atotal is ",f6.4,&
& " and over Etotal is ",f6.4)')CAtotal,CEtotal
write(12,'("Control coeff. for gs over Atotal is ",f6.4,&
& " and over Etotal is ",f6.4)')CAtotal,CEtotal
if(iinfoextra.ne.0)close(13)
! Enable a rerun with all else the same but new values of mBB, hBB, bBB, Vcmax25
read(9,*)Vcmax25,mBB,hBB,bBB
Vcmax25=1.e-6*Vcmax25
if(Vcmax25.gt.0.)then
! Get ready to reuse the data on penetration of diffuse light
rewind 10
do i=1,25
read(10,*)dummy ! 1st 25 lines are locations of r, theta, phi
enddo
goto 8000
endif
close(10)
close(12)
close(14)
close(16)
stop
end
!
!
subroutine BBsolve(scalinghere)
! Type declarations for variables (all are simple, not arrays)
integer nruns ! (-) number of binary searches to do to narrow down the
! solution for gs
logical BBswitch ! (-) switch: 1 = use Ball-Berry humidity response,
! 0 = don't
real Aleaf ! (mol m-2 s-1) CO2 assim. rate of leaf
real bBB ! (mol m-2 s-1) intercept in Ball-Berry-Dewar relation for
! stomatal conductance
real Cs ! (Pa) = partial pressure of CO2 at leaf surface, beneath
! leaf boundary layer
real eair ! (Pa) partial pressure of water vapor in ambient air around leaf
real eL ! (Pa) saturated partial pressure of w.v., at leaf temperature
real facgs ! (-) factor by which the Ball-Berry-Dewar slopr for stomatal conductance
! is cut by water stress
real gbair ! (mol m-2 s-1) conductance of leaf boundary layer, molar units
real gs0 ! (mol m-2 s-1) initial estimate of gs, stomatal conductance
real gshi,gslo,gsmid ! (mol m-2 s-1) upper, lower, midpoint in
! current search for gs
real gsmin ! (mol m-2 s-1) minimum value of gs allowed physically
real gsminus,gsplus ! (mol m-2 s-1) refined lower, upper limits in
! search for gs, during binary search
real gtotCO2 ! (mol m-2 s-1) conductance of stomata plus leaf b. l.
! for CO2
real gtotw ! (mol m-2 s-1) same, for water vapor
real ha ! (-) relative humidity of air
real hBB ! (-) exponent of surface relative humidity in the response
! of stomatal conductance to environment
real hs ! (-) relative humidity at leaf surface, beneath leaf
! boundary layer
real hs0 ! (-) notional (average) value of hs
real IBBhi,IBBlo,IBBmid ! (mol m-2 s-1) Ball-Berry index for
! gshi, gslo, gsmid
real mBB ! (mol m-2 s-1) slope in Ball-Berry model
real Pair ! Pa) total pressure of air
real PPFD ! (mol m-2 s-1) PAR irradiance on leaf
real Qrad ! (W m-2) total radiative energy gain per area of leaf
real Rdswitch ! (-) turns on or off the inclusion of 2*Rd in the
! Ball-Berry-Dewar equation for stomatal conductance
real scalinghere ! (-) downscaling of PS capacity at current canopy location
real Tair ! (oC) air temperature inside canopy
real testhi,testlo,testmid ! mol m-2 s-1) error in solution for
! gs at gshi, gslo, gsmid
real TL ! (oC) leaf temperature
real tol ! (mol m-2 s-1) tolerance for error in solving for gs
! Type declaration for do-loop index
integer n
! common blocks
common/main2BBsolve/mBB,Rdswitch,hBB,bBB,BBswitch,ha,hs0,tol,facgs
common/main2all/Pair,Tair,TL,eair,gbair,Aleaf,Ci,Cs,&
& gtotw,gtotCO2, PPFD, Qrad,CioverCa,Rd
common/BBsolve2Tleaf/eL
common/ecalc2BBsolve/gsplus,gsminus,hs
!Test
! write(6,'("In BBsolve, scalinghere=",f6.4)')scalinghere
!
!
! Method: assume a gs (at two endpoints)
! Calculate Tleaf, thus, E, (hs,) Vcmax, and Q00
! Solve for Asat(Ci),TL)=gtotCO2*(Ca-Ci)/P
! Solve for ALL(Ci)=Q00*(Ci-gamma)/(Ci+2*gamma)*IL=gtot*(Ca-Ci)/P
! Choose lower of two
! Test if gs >, <, or = m (Agross+Rd) hs/Ci + b
!
! Get a starting guess for gs: take Tleaf=Tair, solve for A; then
! compute Cs and compose A*hs0/Ci to estimate gs; take 0.7 and
! 1.5x this as endpoints
TL=Tair
! Need an estimate of gtotCO2; set it for infinite gs
gtotCO2=1./(1.+1.37/gbair)
call CiACssolve(scalinghere)
if(BBswitch)then
hs=ha
else
hs=hs0
endif
! Use air humidity as guess of hs, if using full BB response
gs0=facgs*(mBB*Aleaf*hs**hBB*Pair/Cs) ! Not using Agross+Rd here; we overcorrect, to be sure
gs0=amax1(gs0,bBB)
gslo=0.7*gs0
gshi=1.5*gs0
! Do gslo first, in a generic form
gtotCO2=1./(1.37/gbair+1.6/gslo)
TL=Tleaf(gslo)
! Solve Vcmax*(Ci-gamma)=(Ci+Kco)*gtotCO2*(Ca-Ci)/Pair for Ci, then
! substitute to get A
! Then solve Q00*(Ci-gamma)/(Ci+2*gamma) = gtotCO2*(Ca-Ci)/Pair for Ci,
! then substitute Ci to get new A; choose lower of this, first A
call CiACssolve(scalinghere)
if(BBswitch)then
hs=(gbair*eair/eL+gslo)/(gbair+gslo)
else
hs=hs0
endif
IBBlo=facgs*(mBB*(Aleaf+Rdswitch*2.*Rd)*hs**hBB*Pair/Cs+bBB)
testlo=gslo-IBBlo
! Now do gshi
95 gtotCO2=1./(1.37/gbair+1.6/gshi)
TL=Tleaf(gshi)
call CiACssolve(scalinghere)
if(BBswitch)then
hs=(gbair*eair/eL+gshi)/(gbair+gshi)
else
hs=hs0
endif
IBBhi=facgs*(mBB*(Aleaf+Rdswitch*2.*Rd)*hs**hBB*Pair/Cs+bBB)
testhi=gshi-IBBhi
! Enter decision loop
! Test if gs > mBB*(Aleaf+2*Rd)*hs**hBB/Cs+b
97 if(testlo.lt.0.)then
if(testhi.lt.0.)then
! write(6,90)
! 90 format(1x,'No root evident in interval')
! Fix this problem: remake the search interval
gslo=gshi
testlo=testhi
gshi=1.5*gslo
go to 95
elseif(testhi.eq.0.d0)then
! write(6,100)gshi
gsplus=gshi ! be sure to return answer in gsplus!
gsminus=gshi
goto 200
100 format(1x,'Exact solution at gs=',f12.6)
else
gsplus=gshi
gsminus=gslo
endif
elseif(testlo.eq.0.d0)then
! write(6,100)gslo
gsplus=gslo
gsminus=gslo
goto 200
else
if(testhi.lt.0.)then
gsminus=gshi
gsplus=gslo
elseif(testhi.eq.0.)then
! write(6,100)gshi
gsplus=gshi
gsminus=gshi
goto 200
else
! write(6,90)
! Fix this problem: remake the search interval
gshi=gslo
gsmin=amax1(bBB,0.003)
! ...but be sure that gs does not become too small!! 30 June 01
if(gshi.le.gsmin)then
gsplus=gsmin
gsminus=gsmin
gtotCO2=1./(1.37/gbair+1.6/gsmid)
TL=Tleaf(gsplus)
call CiACssolve(scalinghere)
if(BBswitch)then
hs=(gbair*eair/eL+gsmid)/(gbair+gsmid)
else
hs=hs0
endif
return
endif
testhi=testlo
gslo=0.6*gshi
! Solve for testlo now; saves redoing calcs. from stmt. 95
gtotCO2=1./(1.37/gbair+1.6/gslo)
TL=Tleaf(gslo)
call CiACssolve(scalinghere)
if(BBswitch)then
hs=(gbair*eair/eL+gslo)/(gbair+gslo)
else
hs=hs0
endif
IBBlo=facgs*(mBB*(Aleaf+Rdswitch*2.*Rd)*hs**hBB*Pair/Cs+bBB)
testlo=gslo-IBBlo
go to 97
endif
endif
! Now run through binary search; 1st, figure out how many times
nruns=alog((gshi-gslo)/tol)/0.693
do n=1,nruns
gsmid=0.5*(gsplus+gsminus)
gtotCO2=1./(1.37/gbair+1.6/gsmid)
TL=Tleaf(gsmid)
call CiACssolve(scalinghere)
if(BBswitch)then
hs=(gbair*eair/eL+gsmid)/(gbair+gsmid)
else
hs=hs0
endif
IBBmid=facgs*(mBB*(Aleaf+Rdswitch*2.*Rd)*hs**hBB*Pair/Cs+bBB)
testmid=gsmid-IBBmid
if(testmid.lt.0.)then
gsminus=gsmid
elseif(testmid.eq.0.)then
! write(6,100)gsmid
gsplus=gsmid
gsminus=gsmid
goto 200
else
gsplus=gsmid
endif
enddo
200 return
end
!
!
function Tleaf(gs)
! Type declaration for variables (all simple; no arrays)
real A ! (mol m-2 s-1) CO2 assim. rate of leaf
real Bcc ! (W m-2 K-1) heat transfer coefficient for leaf
real Bnet ! (W m-2 K-1) derivative of total energy balance w/r to Tleaf
real BTIR ! (W m-2K-1) derivative w/r to T of QTIR flux density from leaf
real Cs ! (Pa) partial pressure of CO2 at surface of leaf, beneath
! leaf boundary layer
real dTleaf ! (oC or K) increment in Tleaf
real eair ! (Pa) partial pressure of water vapor in ambient air
! around leaf
real eL ! (Pa) saturated partial pressure of w.v., at leaf temperature
real Eleaf ! (mol m-2 s-1) transpiration rate of chosen leaf
real ep ! (Pa K-1) derivative of eL w/r to Tleaf
real gbair ! (mol m-2 s-1) conductance of leaf boundary layer, molar units
real gtotCO2 ! (mol m-2 s-1) conductance of stomata plus leaf b. l.
! for CO2
real gtotw ! (mol m-2 s-1) conductance of stomata plus leaf b. l.,
! for w.v.
real Hvap ! (J mol-1) heat of vap. of water
real Pair ! (Pa) total pressure of air
real PPFD ! (mol m-2 s-1) PAR irradiance on leaf
real Qmcc ! (W m-2) rate of sensible heat loss from leaf
real QmE ! (W m-2) rate of latent heat loss from leaf
real Qmrad ! (W m-2) rate of TIR radiant heat loss from leaf
real Qnet ! (W m-2) net energy flux imbalance of leaf
real Qrad ! (W m-2) total radiative energy gain per area of leaf
real sigma ! (W m-2 K-4) Stefan-Boltzmann constant
real Tair ! (oC) air temperature inside canopy
real TL ! (oC) ! leaf temperature, as passed in common
real Tleaf ! (oC) leaf temperature
! common blocks
common/main2Tleaf/sigma,Hvap,BTIR
common/main2all/Pair,Tair,TL,eair,gbair,A,Ci,Cs,&
& gtotw,gtotCO2, PPFD, Qrad, CioverCa,Rd
common/BBsolve2Tleaf/eL
common/Tleaf2Ecalc/Bcc,Qmcc,Eleaf
!
Tleaf=Tair ! initial guess
gtotw=1./(1./gs+1./gbair)
10 Qmrad=2.*sigma*(Tleaf+273.2)**4
BTIR=4.*Qmrad/(Tleaf+273.2)
eL=esat(Tleaf)
Eleaf=gtotw*(eL-eair)/Pair
QmE=Hvap*Eleaf
Qmcc=Bcc*(Tleaf-Tair)
Qnet=Qrad-Qmrad-QmE-Qmcc
ep=esatpr(Tleaf)
Bnet=BTIR+Hvap*gtotw*ep/Pair+Bcc
dTleaf=Qnet/Bnet
Tleaf=Tleaf+dTleaf
if(abs(dTleaf).gt.0.01)go to 10
! if(Tleaf.gt.49.)then
! write(6,'("Tleaf=",f6.2,"; Qmrad,BTIR,QmE,Qmcc,Qnet,Bnet=",6f9.1)')&
! &Tleaf,Qmrad,BTIR,QmE,Qmcc,Qnet,Bnet
! pause
! endif
return
end
!
!
subroutine CiACssolve(scalinghere)
! Type declaration for simple variables
integer iflag ! not used currenlty
integer iter ! (-) iteration in solution for eair, Tair in canopy
real Aleaf ! (mol m-2 s-1) CO2 assimilation rate of leaf under
! consideration
real ALL ! (mol m-2 s-1) light-limited rate of CO2 assimilation
real Amax ! (mol m-2 s-1) light-saturated rate of CO2 assimilation
real c0,c1,c2,c3,c4 ! (monster units) coefficients in quartic solution
! for Ci
real Ca ! (Pa) partial pressure of CO2 in ambient air ! (at ref. height
! above canopy)
real Ci ! (Pa) partial pressure of CO2 inside leaf
real Cs ! (Pa) partial pressure of CO2 at surface of leaf, beneath
! leaf boundary layer
real dCi ! (Pa) increment in solution for Ci
real denom ! (J mol-1 K) intermediate in calc. of temperature activation
! functions for CO2 assimilation
real dfac ! (K2) similar
real eair ! (Pa) partial pressure of w.v. inside canopy
real F ! (monster units) function whose root is sought to solve for Ci
real G ! (Pa) shorthand for gamma
real gamma ! (Pa) CO2 compensation partial pressure
real gbair ! (mol m-2 s-1) conductance of leaf boundary layer, molar units
real gt ! (mol m-2 s-1) shorthand for gtotCO2
real gt2 ! (mol2 m-4 s-2) shorthand for square of gtotCO2
real gtotCO2 ! (mol m-2 s-1) conductance of stomata plus leaf b. l.
! for CO2
real gtotw ! (mol m-2 s-1) conductance of stomata plus leaf b. l.,
! for w.v.
real IL ! (mol m-2 s-1) PAR irradiance on leaf
real K ! (Pa) shorthand for Kco
real Kc ! (Pa) Michaelis constant for binding of CO2 to Rubisco, at
! current leaf temperature
real Kc25 ! (Pa) Michaelis constant for binding of CO2 to Rubisco,
! at ref. T of 25oC
real Kco ! (Pa) effective Michaelis constant for binding of CO2 to
! Rubisco, in presence of current O2 partial pressure
real Ko ! (Pa) Michaelis constant for binding of O2 to Rubisco, at
real Ko25 ! (Pa) Michaelis constant for binding of O2 to Rubisco,
! at ref. T of 25oC
real P ! (Pa) shorthand for Pair
real P2 ! (Pa2) shorthand for square of Pair
real Pair ! (Pa) total pressure of air
real PPFD ! (mol m-2 s-1) PAR irradiance on leaf
real Q00 ! (mol CO2 [mol photons]-1) initial quantum yield of CO2
! assimilation, at saturating CO2
real QI ! (mol CO2 s-1) shorthand for Q00*IL
real Qrad ! (W m-2) total radiative energy gain per area of leaf
real Rd ! (mol m-2 s-1) dark respiration rate of leaf
real RdTmean ! (mol m-2 s-1) leaf respir. rate at Tleaf Tmean
real scalinghere ! (-) downscaling of PS capacity at current canopy location
real T ! (-) shorthand for thetaPS
real Tair ! (oC) air temperature inside canopy
real thetaPS ! (-) convexity parameter, for sharpness of transition from
! light-limited to light-sat'd CO2 assimilation
real TL ! (oC) leaf temperature
real Tmean ! (oC) mean leaf T during photoperiod over last 2 weeks
! ! (T to which leaf respiration has acclimated)
real V ! (mol m-2 s-1) shorthand for Vcmax
real Vcmax ! (micromol m-2 s-1, immed. converted to mol m-2 s-1) photosyn.
! capacity per leaf area! (light- and CO2-saturated)
real Vcmax25 !(micromol m-2 s-1, immed. converted to mol m-2 s-1)
! photosyn. capacity per leaf area (light- and CO2-saturated), at 25oC
! common blocks
common/main2Cisolve/Kc25,Ko25,Vcmax25,Ca,Q00,&
& thetaPS,RdTmean,Tmean,iflag
common/main2all/Pair,Tair,TL,eair,gbair,Aleaf,Ci,Cs,&
& gtotw,gtotCO2, PPFD, Qrad, CioverCa,Rd
common c4,c3,c2,c1,c0
!
dfac=298.*8.314
denom=dfac*(TL+273.2)
Vcmax=scalinghere*Vcmax25*exp(64800.*(TL-25.)/denom)
if(TL.gt.60.)then
Vcmax=amax1(Vcmax*(1.-0.05*(TL-49.)-0.005*(TL-49.)**2),0.)
write(6,'("TL=",f6.2,"; Vcmax=",6Pf6.1,"; TL=",f6.2)')TL,Vcmax,TL
pause
endif
Kc=Kc25*exp(59400.*(TL-25.)/denom)
Ko=Ko25*exp(36000.*(TL-25)/denom)
Kco=Kc*(1.+21000./Ko)
gamma=3.69+(TL-25.)*(0.188+0.00036*(TL-25.))
Rd=scalinghere*RdTmean*exp(0.07*(TL-Tmean))
! Solve the Johnson-Thornley equation for A,
! thetaPS*A**2 - A*(Amax+Q0*IL) - sqrt((Amax+Q0*IL)**2-4*thetaPS*
! Amax*Q0*IL)) = 0, with
! Amax=Vcmax*(Ci-gamma)/(Ci+Kco),
! Q0=Q00*(Ci-gamma)/(Ci+2*gamma)
! Write A as gtotCO2*(Ca-Ci).
! After writing all terms explicitly in terms of Ci, we multiply
! through by terms in the denominator, obtaining a quartic in Ci.
! Solve for Ci by a binary search.
! Substitute this value of Ci in A = gtotCO2*(Ca-Ci) to get A.
! Solve trivially for Cs = Ca - 1.37*A/gbair.
!
! Now compute terms in quartic; first, convert long names to short names
! (rather than converting names in well-debugged code and possibly
! introducing errors)
gt=gtotCO2
P=Pair
IL=PPFD
K=Kco
G=gamma
V=Vcmax
T=thetaPS
b=2. !!!!
!
gt2=gt*gt
P2=P*P
QI=Q00*IL
! Start the iterations of the Newton-Raphson method
Ci=Ca
! !!!! From here to about 20 lines down is new method, using quadratic
! solution vs. quartic equation that introduces false roots for Ci
! (must be > Ca when Q0*IL is less than Rd)
nfail=0 ! number of failures to find physiologically realistic
90 iter=0
100 Amax = V*(Ci-G)/(Ci+K)
Q0IL = QI*(Ci-G)/(Ci+b*G)
AJT = (0.5/T)*(Amax+Q0IL - sqrt((Amax+Q0IL)**2-4.*T*Amax*Q0IL))
F = AJT - Rd - gt*(Ca-Ci)/P
Amaxp=V*(K+G)/(Ci+K)**2
Q0ILp = QI*(b+1.)*G/(Ci+b*G)**2
AJTp = (AJT*(Amaxp+Q0ILp) - Amaxp*Q0IL - Amax*Q0ILp)/&
& (2.*T*AJT - (Amax+Q0IL))
Fprime = AJTp +gt/P
dCi=-F/Fprime
! Hobble the Newton-Raphson step to be less than 1/2 current value of Ci
dCiabs=abs(dCi)
dCi=sign(amin1(dCiabs,0.5*Ci),dCi)
Ci=Ci+dCi
iter=iter+1
if(abs(dCi).gt.0.01)goto 100
! Converged on a math'l soln. - test if it's physiological for low-light cases
! That is, test if Q0*IL < Rd, need Ci > Ca; if Ci is not greater than Ca,
! reset Ci to 1.5*Ca and do iterations again
! This may be unnecessary when we use the quadratic equation, but it's safe
if( (QI*(Ci-G).lt.Rd*(Ci+b*G) ) .and.(Ci.lt.Ca))then
nfail=nfail+1
Ci=(1.+0.2*float(nfail))*Ca
goto 90
endif
! !!!!
! Compute A and compare to Amax(Ci) and Q0(Ci)*IL
Amax=Vcmax*(Ci-gamma)/(Ci+Kco)
ALL=Q00*(Ci-gamma)*PPFD/(Ci+gamma)
! Solve quadratic for A:
Aleaf=(0.5/T)*(Amax+ALL-sqrt((Amax+ALL)**2-4.*T*&
& Amax*ALL)) - Rd
Cs=Ca-1.37*Aleaf*Pair/gbair
CioverCa=Ci/Ca
return
end
!
!
FUNCTION ESAT(T)
real T ! (oC) temperature (of leaf)
ESAT=610.78*EXP(17.269*T/(T+237.3))
! = SAT'D VAPOR P. OF WATER (PASCALS) AT AMBIENT TEMP. T (DEG. C)
RETURN
END
!
!
FUNCTION ESATPR(T)
real T ! (oC) temperature (of leaf)
ESATPR=610.78*EXP(17.269*T/(T+237.3))*4097.9&
& /(T+237.3)**2
! = (D/DT) PSAT(T)
RETURN
END
!
!
function julian(date,ndm)
! Type declarations for variables
character*8 date ! (text) current date, in mm-dd-yy format
integer iac ! (-) = accumulator of number of days in preceding months
integer nday ! (-) = number of day within current month
integer nday1 ! (-) = 1st digit in nday
integer nday2 ! (-) = 2nd digit in nday
integer nmonth ! (-) = number of current month in year
integer nmonth1 ! (-) = 1st digit of nmonth
integer nmonth2 ! (-) = 2nd digit of nmonth
! Type declaration of one array
integer ndm(12) ! (-) = number of days in each month
!
nmonth1=ichar(date(1:1))
nmonth2=ichar(date(2:2))
nmonth=10*(nmonth1-48)+nmonth2-48
nday1=ichar(date(4:4))
nday2=ichar(date(5:5))
nday=10*(nday1-48)+nday2-48
iac=0
do i=1,12
if(nmonth.eq.i)goto 99
iac=iac+ndm(i)
enddo
99 julian=iac+nday
return
end
!
!
real function Ecalc(Ppendiffuse,ASRlocdiffPAR,ASRlocdiffNIR,ASRlocdirPAR,&
& ASRlocdirNIR)
parameter (maxtrees=100, maxsamples=20)
! Type declarations for simple variables
character*4 infix !(text) infix to be inserted into output file names, making each run's
! filenames unique so that they don't overwrite each other
character*20 scratchname ! (text) expanded name of an output file, including the infix
integer infixlen ! (-) number of nonblank characters in the infix to be inserted into names of
! output files, so that each run in unique and outputs don't overwrite each other
integer ipensave ! (-) switch: 1 = save newly calculated array of
! light-penetration probabilities, 0 = don't
integer ipre ! (-) switch: 1 light-penetration probability arrays
! exist in a file; use them; 0 they don't exist
integer iter ! (-) iteration in solution for eair, Tair in canopy
integer nbins ! (-) number of bins to use in histogram of leaf
! PAR irradiance, IL
integer nd ! (-) number of steps of size 'stepsize' to take through canopy
! while accumulating total length of path within tree crowns
integer ndmax ! (-) max. number of steps ! (prevents one from wasting time
! looking outside all known trees)
integer nphi ! (-) number of azimuth angles to use in sampling locations
! within central tree
integer nradius ! (-) number of radii from center to use in sampling locations
! within central tree
integer ntheta ! (-) number of zenith angles to use in sampling locations
! within central tree
integer ntimes ! (-) number of the current simulation interval
integer ntrees ! (-) = number of trees around the central tree that we will
! describe
real AAtree ! (mol s-1) accumulator for CO2 assim. rate of entire tree
! Converted at end of hourly simulation loop to kg(PSate) h-1
real AAtree10 ! (mol s-1) same, when gs is increased by 10%)
real Aleaf ! (mol m-2 s-1) photosynthetic rate of leaf
real aNIR ! (-) absorptivity of leaf for NIR
real aPAR ! (-) absorptivity of leaf for PAR
real Asum ! (mol m-2 s-1) sum of leaf assim. rate over different
! bins of PAR irradiance
real Asum10 ! (mol m-2 s-1) same, for gs increased by 10%
real Bcc ! (W m-2 K-1) heat transfer coefficient for leaf
real bet ! (m) intermediate in calc. of arbitary tree's projected
! canopy length along sun direction
real cosalphas ! (-) cosine of solar elevation angle
real costhetape ! (-) dot product of direction vectors, in computation
! of distance that direct solar beam travels inside tree canopy
real Cs ! (Pa) partial pressure of CO2 at surface of leaf, beneath
! leaf boundary layer
real D0PAR ! (micromol m-2 s-1) same as D00, but interpolated to
! this specific time from weather data
real dI ! (micromol m-2 s-1) increment in PAR irradiance, between bins
real dpc ! (m) distance of a point from center of selected tree
real dt ! (m) distance along path from chosen point in tree to top of canopy
real dx ! (m) x-component of step along direct beam path while checking
! for interception by trees
real dy ! (m) same, for y-component
real dz ! (m) same, for z-component
real eair ! (Pa) partial pressure of w.v. inside canopy
real fluxdirNIR ! (W m-2) energy flux density in the NIR at this time
real Eleaf ! (mol m-2 s-1) transpiration rate of chosen leaf
real ENIR ! (W m-2) energy flux density incident in the NIR on leaf in
! the shade (taken as uniform among all leaves in the shade - no
! distribution of values)
real ENIRshade ! (W m-2) energy flux density incident in the NIR
real Esum ! (mol m-2 s-1) sum of leaf transpir. rate over different
! bins of PAR irradiance
real Esum10 ! (mol m-2 s-1) same, for gs increased by 10%
real Etree ! (mol s-1) transpiration rate of central tree
real Etree10 ! (mol s-1) same, with gs increased by 10% on all leaves
real fbin ! (-) fraction of total irrad. range per bin
real fd ! (m-1) foliage density inside crowns of trees
real fNIR ! (-) forward scattering by (transmissivity of) leaves in NIR
real fPAR ! (-) forward scattering by (transmissivity of) leaves in PAR
real gbair ! (mol m-2 s-1) conductance of leaf boundary layer, molar units
real gsminus ! (mol m-2 s-1) current lower limit in search for gs,
! during binary search
real gsplus ! (mol m-2 s-1) current upper limit in search for gs,
! during binary search
! real gstest ! (mol m-2 s-1) test value of gs, 10% higher than true sol.
real gtot ! (mol m-2 s-1) conductance of stomata plus leaf b. l.,
! for w.v.
real gtotCO2 ! (mol m-2 s-1) conductance of stomata plus leaf b. l.
! l for CO2
real gtotsum ! (mol m-2 s-1) accumulator for total leaf conductance
! for this location in canopy
real gtottree ! (mol m-2 s-1) =accumulator for conductance of all leaves on tree
real gtotw ! (mol m-2 s-1) conductance of stomata plus leaf b. l.,
! for w.v. - only passed to keep common block intact
real hcan ! (m) max. height of canopy ! (among all trees)
real hs ! (-) relative humidity at leaf surface, beneath leaf boundary layer
real Hsum ! (W m-2) sensible heat flux from leaf at this location
real Htree ! (mol s-1) transpir. rate of central tree
real I0PAR ! (micromol m-2 s-1) I00, after interpolation from weather data
real Iadd ! (micromol m-2 s-1) average PAR irradiance for leaves at this
! location in canopy
real IL ! (micromol m-2 s-1) PAR irradiance on leaf, in particular bin
! real Isum ! (micromol m-2 s-1) PAR irradiance sum for leaf at this
! location
real ILshade ! (micromol m-2 s-1) PAR irradiance on leaf in the shade (taken
! as uniform among all leaves in the shade - no distribution of values)
real Kdiff ! (-) extinction coefficient for diffuse skylight
real kpen ! (-) average cosine of angle between leaves and sun
real lpen ! (-) probability of direct-beam penetration to this
! l location in canopy at this simulation interval
real Pair ! (Pa) total pressure of air
real Pclear ! (-) fraction of time in simulation interval that sky is
! clear ! (estimated)
real phican ! (radians) currently selected azimuth within central tree
real Ppenmin ! (-) cutoff (min.) value of Ppen, the penetration
! probability for solar radiation (needed to handle lower limits in
! calculation of absorbed scattered radiation)
real PPFD ! (mol m-2 s-1) PAR irradiance on leaf
real Qmcc ! (W m-2) rate of sensible heat loss from leaf
real Qrad ! (W m-2) total radiative energy gain per area of leaf
real Qsky ! (W m-2) TIR energy flux density from sky
real QSW ! (W m-2) total shortwave ! (PAR + NIR) energy flux density
! absorbed by this leaf
real QTIR ! (W m-2) total TIR energy flux density absorbed by this leaf
real Qveg ! (W m-2) TIR energy flux density from surface at T of
! vegetation
real ra ! (s m-1) leaf b.l. resistance in old units
real radius ! (m) radial distance in central tree, at which leaf
! performance is being sampled
real rdsol ! (m) radius of canopy along solar direction
real reflsoilNIR ! (-) reflectivity of soil (ground, incl. cover) in NIR
real reflsoilPAR ! (-) reflectivity of soil (ground, incl. cover) in PAR
real rhoair ! (mol m-3) molar density of air
real scalinghere ! (-) value of the PS-scaling factor at the current
! canopy location
real sigma ! (W m-2 K-4) Stefan-Boltzmann constant
real sinalphas ! (-) sine of solar elevation angle
real Skyview ! (-) fraction of diffuse skylight reaching this location
! in central tree
real st ! (m) cumulative path of sun through crowns of all trees in stand
! along current solar direction
real stepsize ! (m) increment in path search along solar beam direction
real stmax ! (m) max. value of st, at which search quits ! (penetration
! probability is very low)
real sumI ! (mol m-2 s-1) sum of PAR irradiance over all leaves over all
! simulation intervals
real Tair ! (oC) air temperature inside canopy
real Tleafval ! (oC) leaf temperature, renamed so it can be passed
! in common
real TLsum ! (W m-2) accumulator for average !(leaf T) x !(leaf conductance)
! at this location
real TLtree ! (oC) sum of ! (leaf T) x ! (leaf conductance) over all
! locations in tree
real Tsky ! (oC) effective radiative temperature of sky
real Tveg ! (oC) effective radiative temperature of vegetation
real Wmol ! (J mumol-1) ! conversion from quantum flux density in PAR
! to energy flux density in PAR
real wt ! (-) fraction of leaf area at this location in shade only
real x ! (radians) sine of solar azimuth; also ! (m), x-coordinate of
! current point in central tree
real xnew ! (m) x-coordinate of current position in whole canopy during
! search for interception by crowns
real xo ! (m) x-coordinate of currently-sampled point in central tree
real y ! (m) y-coordinate of current point in central tree
real ynew ! (m) y-coordinate of current position in whole canopy during
! search for interception by crowns
real yo ! (m) y-coordinate of currently-sampled point in central tree
real z ! (m) z-coordinate of current point in central tree
real znew ! (m) z-coordinate of current position in whole canopy during
! search for interception by crowns
real zo ! (m) z-coordinate of currently-sampled point in central tree
! Type declaration and dimensioning of arrays
real ASRlocdiffNIR(maxsamples,maxsamples,maxsamples) ! (-) absorbed scattered
! radiation at each canopy location, arising from unit input of diffuse
! sky radiation (in the NIR) at the top of the canopy - estimated with
! uniform layered canopy
real ASRlocdiffPAR(maxsamples,maxsamples,maxsamples)! (-) absorbed scattered
! radiation at each canopy location, arising from unit input of diffuse
! sky radiation (in the PAR) at the top of the canopy - estimated with
! uniform layered canopy
real ASRlocdirNIR(maxsamples,maxsamples,maxsamples)! (-) absorbed scattered
! radiation at each canopy location, arising from unit input of direct-
! beam radiation (in the NIR) at the top of the canopy, at the current
! solar geometry - estimated with uniform layered canopy
real ASRlocdirPAR(maxsamples,maxsamples,maxsamples)! (-) absorbed scattered
! radiation at each canopy location, arising from unit input of direct-
! beam radiation (in the PAR) at the top of the canopy, at the current
! solar geometry - estimated with uniform layered canopy
real atree(maxtrees) ! (m) horizontal full axis of ellipsoid for given tree
real btree(maxtrees) ! (m) vertical full axis of ellipsoid for given tree
real costheta(maxsamples) ! (-) cosine of zenith angle of point being sampled in
! central tree
real lpenshort(maxsamples) ! (-) lpenstore for current iradius,itheta, over all iphi
real lpenstore(maxsamples,maxsamples,maxsamples) !
! (-) penetration probability for direct beam
! to all canopy locations, for current simulation interval
real phi(maxsamples) ! (radians) sky azimuths being sampled for diffuse skylight
real phic(maxtrees) ! (radians) azimuthal inclination of ellipsoid of each tree
real Ppendiffuse(maxsamples,maxsamples,maxsamples)
! (-) total penetration fraction of diffuse
! skylight, for each location in central tree
real r(maxsamples,maxsamples) ! (m) array of radii being sampled in
! central tree, for each choice of theta and phi
real scaling(maxsamples,maxsamples,maxsamples) ! (-) Downscaling factor for
! photosynthetic and respiratory capacity at canopy location, tracking
! mean exposure to light
real sintheta(maxsamples) ! (-) sine of zenith angle of point being sampled in
! central tree
real thetac(maxtrees) ! (radians) zenith inclination of ellipsoid of each tree
real xtree(maxtrees) ! (m) x-coordinate of center of each tree
real ytree(maxtrees) ! (m) y-coordinate of center of each tree
real ztree(maxtrees) ! (m) z-coordinate of center of each tree
!
! Type declaration for do-loop indices
integer i,iphi,iradius,itheta,nb,nt
integer idummy !DDDDDDDDDDDDDDDDD
real*8 xx,ran !DDDDDDDDDDDDDDDDD
!!!
logical histotime
real enlocations
! Declare storage arrays for histograms:
real histoILnoscatt(120), histoILall(120), histoTLraw(120)
real histoAraw(120), histoEraw(120), histogtotraw(120)
real histoTL_Awtd(120), histoTL_Ewtd(120)
real histogtot_Awtd(120), histogtot_Ewtd(120)
real histoIL_Awtd(120), histoIL_Ewtd(120)
! Later, I might want histograms of TIR flux density
!!!
! common blocks
common/main2Ecalc/nbins,ipre,ndmax,iter,ipensave,icall,JDdiff,& ! integers
& bscale,aNIR,aPAR,Qsky,Etree,Etree10,AAtree,AAtree10,gtottree,& ! real variables
& TLtree,sumI,Htree,Ppenmin,fPAR,fNIR,reflsoilPAR,reflsoilNIR,Kdiff,&
& enlocations,Citree,hstree,&
& lpenstore,scaling,& ! real arrays
& histotime,& ! logical variable
& infix, infixlen ! info about infix to be inserted into name of output file, histo.txt
common/Tleaf2Ecalc/Bcc,Qmcc,Eleaf
common/special/ntimes
common/main2all/Pair,Tair,Tleafval,eair,gbair,Aleaf,Ci,Cs,&
& gtotw,gtotCO2, PPFD, Qrad, CioverCa,Rd
common/ecalc2BBsolve/gsplus,gsminus,hs
common/centertreedescrip/nradius,ntheta,nphi,costheta,sintheta,&
& r,phi,xo,yo,zo
common/alltrees/ntrees,atree,btree,thetac,phic,xtree,ytree,ztree,&
& hcan,kpen,fd
common/other/sinalphas,cosalphas,phis,ra,rhoair,D0PAR,I0PAR,&
& Pclear,wt,stepsize,stmax
!
sigma=5.67e-8
nloc=nradius*ntheta*nphi
! Initialize sum of transpiration...and for sensible heat flux
Etree=0. ! sum of E over all locations on tree
Etree10=0. ! Same, with 10% higher gs at all locations in canopy
Htree=0.
gtottree=0.
TLtree=0.
hstree=0.
Citree=0.
AAtree=0.
AAtree10=0.
! Initialize accumulator for irradiance on all locations
sumI=0.
! Initialize counter over canopy locations
iloc=0
! Nested loops over locations within canopy - first pass is to calculate
! (or retrieve, from previous runs) the direct-beam penetration
! probabilities to all locations in the canopy. We need them all at
! once in order to estimate the light scattered from other leaves
! and from the ground
if((ipre.eq.0).and.(ipensave*iter.le.1))then ! Calculate the pen. probs.
sumPpen=0. !DDDDDDDDDDDDDDDDDDDDDDD
do iradius=1,nradius
do itheta=1,ntheta
costhetacan=costheta(itheta)
sinthetacan=sintheta(itheta)
radius=r(itheta,iradius) ! radius of sampled point
do iphi=1,nphi
iloc=iloc+1
! pause
! If the array of light-penetration probabilities is available from
! a file and has been read in, OR if these probabilities have been
! done in a previous iteration #1, skip this calculation
phican=phi(iphi)
! First compute coordinates of this point with respect to the center
! of the central canopy as if this were vertical
xo=radius*sinthetacan*cos(phican)
yo=radius*sinthetacan*sin(phican)
zo=radius*costhetacan
! Compute the coordinates for the real canopy by rotating the point
x=cos(phic(1))*(xo*cos(thetac(1))+zo*sin(thetac(1)))&
& -yo*sin(phic(1))
y=sin(phic(1))*(xo*cos(thetac(1))+zo*sin(thetac(1)))&
& +yo*cos(phic(1))
z=-xo*sin(thetac(1))+zo*cos(thetac(1))
! Check if neighboring trees are in the light path
! First calculate total distance between canopy point and top of canopy
dt=(hcan-z)/sinalphas
!
! Initializing light path length and coordinates of new point
st = 0.0
xnew=x
ynew=y
znew=z
dx=stepsize*cosalphas*cos(phis)
dy=stepsize*cosalphas*sin(phis)
dz=stepsize*sinalphas
! Compute the number of distance increments through the light path
nd = ifix(dt/stepsize)+1
nd=min0(nd,ndmax)
! Loop for distance increments
do i=1,nd
! Compute increments in x, y, and z through light path
! Compute coordinates of new point
xnew=xnew+dx
ynew=ynew+dy
znew=znew+dz
! Loop for trees in the quadrant
do nt=1,ntrees
! Compute distance from new point to the center of each tree
dpc=sqrt((xnew-xtree(nt))**2+(ynew-ytree(nt))**2+&
& (znew-ztree(nt))**2)
! Compute dot product of the direction vector of the canopy symmetry
! axis with the direction vector from the center of the canopy to
! the new point
costhetape=(sin(thetac(nt))*cos(phic(nt))*(xnew-&
& xtree(nt))+sin(thetac(nt))*sin(phic(nt))*(ynew-ytree(nt))+&
& cos(thetac(nt))*(znew-ztree(nt)))/dpc
! Compute the canopy radius for the sun path direction
bet=1.-btree(nt)*btree(nt)/(atree(nt)*atree(nt))
! It is not necessary to convert from canopy axis to canopy
! semiaxis in this last equation
rdsol=0.5*btree(nt)/sqrt(1.-bet*costhetape**2) ! 0.5 is
! to convert to canopy semiaxis
! Compare the distance from the new point to the canopy center with
! the canopy radius for the sun path direction
if(dpc.le.rdsol)then
st =st+stepsize
! write(6,'("At x,y,z =",3f8.2,", we are in tree ",&
! & "number",i3)')xnew,ynew,znew,nt
endif
enddo ! End loop over different trees
! Compare st with stmax for a light penetration probability= 0.002
if (st.gt.stmax)goto 801
enddo ! End loop over distance increments
!
! Calculate the light penetration probability
801 lpen = EXP(-kpen * fd * st)
lpenstore(iradius,itheta,iphi)=lpen
scaling(iradius,itheta,iphi)=1.
! Above: if penetration probabilities have never been calculated, one cannot
! compute scaling functions for PS capacity, so they must default to 1.0
sumPpen=sumPpen+lpen
enddo ! End loop over azimuthal angles in canopy, phi
enddo ! End loop over zenith angles in canopy, theta
enddo ! End loop over radii in canopy
! write(6,'(" AVERAGE Ppen=",f8.3)')sumPpen/float(nloc)
!
endif ! End of block-if that tests if these calcs. have been done before
! Now handle the arrival of diffuse scattered light at this canopy location,
! arising from scattering by leaves and soil of solar inputs (diffuse
! skylight and direct beam, in both NIR and PAR)
! For scattering of the diffuse skylight, we have, from earlier calculations done
! only once in MAIN, the fractional propagation to
! the current location, as ASRlocdiff(NIR or PAR). Multiply by the appropriate
! direct or diffuse flux density to get the contribution, as flux density
! Now we need to calculate how the scattering of the direct beam adds flux density.
! This depends on the solar geometry at this hour, so we have to call for the
! calculation each hour and day of the simulation
! MAIN has already passed to Ecalc here the values of D0PAR and I0PAR, in the
! common block "other" - we can scale ASR from the direct beams readily
! Note that ASRlocdirPAR will be in quantum flux units, mol m-2 s-1
Wmol=0.22
! If the sun is above the horizon and it's not overcast (I0PAR>0), then
! calculate the absorbed scattered radiation (PAR and NIR) arising from
! the direct solar beam entering the canopy
!
! We only need to do this calculation once (for all locations) in Ecalc, at the
! first iteration at this day and hour...provided that we make the values
! nonvolatile by having them passed to MAIN, using explicit arguments
! Actually, we could store these along with Ppenddiffuse and Ppendir; that will
! take some extra coding that I'm not going to do before all the debugging
if(iter.eq.1)then
atreepass=atree(1) ! Pick up major axis of central tree
if(I0PAR.gt.0.)then
call calc_ASR_EC(2,lpenstore,Ppenmin,ASRlocdirPAR,&
& aPAR,fPAR,reflsoilPAR,0.,1.,atreepass,Kdiff)
! Results, in ASRlocdirPAR, are per unit of dir. beam PAR
! ASRlocdirNIR will be in same form, per unit of dir. beam NIR; we assume equal energy flux in
! the PAR and NIR spectra, so that energy flux in the NIR = that in the PAR =
! Wmol*I0PAR or Wmol*D0PAR, according to whether we are considering the direct
! beam or the diffuse beam (actually, diffuse flux in the NIR is much less than
! in the PAR, but the rescaling is a little tricky and makes little difference in
! actual energy balance of leaves
! EfluxdirNIR=Wmol*I0PAR - no longer needed here
! write(6,'("Calculating ASRlocdirNIR in Ecalc, with EfluxdirNIR=",&
! & f8.2)')EfluxdirNIR
call calc_ASR_EC(2,lpenstore,Ppenmin,ASRlocdirNIR,&
& aNIR,fNIR,reflsoilNIR,0.,1.,atreepass,Kdiff)
else ! There is no direct beam at this hour and day
do iradius=1,nradius
do itheta=1,ntheta
do iphi=1,nphi
ASRlocdirPAR(iradius,itheta,iphi)=0.
ASRlocdirNIR(iradius,itheta,iphi)=0.
enddo ! End loop over phi
enddo ! End loop over itheta
enddo ! End loop over iradius
endif ! End of block-if that tests for presence or absence of direct beam
endif ! End of block-if that tests if these functions have already been
! computed in iteration #1
!!!
! Note the ranges of the IL, TL, A, E, and gs to be recorded:
! Il in (0,2500), by steps of 25
! TL in (20., 45.), by steps of 0.25
! A in (-6., 34.)e-6, by steps of 0.4e-6
! E in (0., 12.5.)e-3, by steps of 0.125e-3
! gtot in (0., 0.5), by steps of 0.005
! At beginning of leaf-flux loop, zero out the histogram accumulators
if(histotime)then
do ihisto=1,100
histoILnoscatt(ihisto)=0.
histoILall(ihisto)=0.
histoTLraw(ihisto)=0.
histoAraw(ihisto)=0.
histoEraw(ihisto)=0.
histogtotraw(ihisto)=0.
histoTL_Awtd(ihisto)=0.
histoTL_Ewtd(ihisto)=0.
histoIL_Awtd(ihisto)=0.
histoIL_Ewtd(ihisto)=0.
histogtot_Awtd(ihisto)=0.
histogtot_Ewtd(ihisto)=0. !RRRRR Revise the type declarations, incl. descriptions
! No need to sum the weights A, E ï¿½ theyï¿½re simply the tree-total values
! Nor is there a need to sum the weights by area separately for histograms;
! the sum is simply the total leaf area on the tree (or, given the way
! I worked it, simply the number of leaf locations)
enddo
endif
!!!
! Now restart the loop over canopy location, armed with all the radiant flux information
do iradius=1,nradius
do itheta=1,ntheta
do iphi=1,nphi
lpen=lpenstore(iradius,itheta,iphi)
scalinghere=scaling(iradius,itheta,iphi)
! Compute some microenvironmental quantities
gbair=rhoair/ra
Bcc=29.*gbair
Tveg=Tair
! Tsky=Tair-25.
! Here, I removed recalculation of Qsky, which should be done with original
! above-canopy values of eair, Tair; this is done in MAIN now
Qveg=sigma*(Tveg+273.2)**4
Skyview=Ppendiffuse(iradius,itheta,iphi)
QTIR=Qsky*Skyview+Qveg*(2.-Skyview)
! Set accumulators for average E, A, and H at this location omly
Asum=0.
Esum=0.
Hsum=0.
gtotsum=0.
TLsum=0.
hssum=0. ! DDDDDDDDDDDDDDDDDDDDDDDDDDDD
Cisum=0. !DDDDDDDDDDDDDDDDDDDDDDDDDDDD
Asum10=0.
Esum10=0.
! Calculate A, E, etc. for cloudy time AND for contribution from
! clear time by shade leaves
! if(histotime)then
! write(6,'("D0PAR,I0PAR=",2f16.9)')D0PAR,I0PAR
! write(6,'("ASRlocdiffPAR,..dirPAR,...diffNIR,...dirNIR=",4f9.5)')&
! &ASRlocdiffPAR(iradius,itheta,iphi),ASRlocdirPAR(iradius,itheta,iphi),&
! &ASRlocdiffNIR(iradius,itheta,iphi),ASRlocdirNIR(iradius,itheta,iphi)
! pause
! endif
ILshade=D0PAR*skyview& ! PAR quantum flux density from 1st interception
! of diffuse skylight incident at top of canopy
& +D0PAR*ASRlocdiffPAR(iradius,itheta,iphi)/aPAR&! ..and from scatt. of
! diffuse skylight incident at top of canopy, in the PAR;
! no factor of aPAR, because this is included in calc. of ASR...
& +I0PAR*ASRlocdirPAR(iradius,itheta,iphi)/aPAR ! ...and from scatt. of
! direct solar beam incident at top of canopy, in the PAR; again,
! no factor of aPAR, because this is included in calc. of ASR...
! Note: dividing by aPAR puts ILshade back on the basis of interception
ENIRshade=Wmol*D0PAR*Skyview& ! NIR energy flux density; made equal to
! that in the PAR, with conversion from quanta to energy
& +Wmol*D0PAR*ASRlocdiffNIR(iradius,itheta,iphi)/aNIR&! ...and from scatt. of
! diffuse skylight incident at top of canopy, in the NIR
& +Wmol*I0PAR*ASRlocdirNIR(iradius,itheta,iphi)/aNIR ! ...and from scatt. of direct
! solar beam incident on canopy, in the NIR; no conversions,
! because ASR... was calculated with energy flux density itself
! Similarly, dividing by aNIR puts the NIR heat loading on an intercepted basis
QSW=Wmol*aPAR*ILshade+aNIR*ENIRshade
Qrad=QTIR+QSW
if(Qrad.gt.2000.)then
write(6,'("shade Qrad=",f8.1,"; Qsky,Qveg,QSW=",3f9.1,"; ntimes,iloc=",2i5)')&
&Qrad,Qsky,Qveg,QSW,ntimes,iloc
endif
Tleafval=Tair ! initial guess for leaf T
PPFD=1.e-6*ILshade ! convert to mol from umol
Aleaf=0.
Eleaf=0.
call BBsolve(scalinghere)
! Calculate weighting function
wt=(1.-Pclear) + Pclear*(1.-lpen)
Asum=Asum+Aleaf*wt
Esum=Esum+Eleaf*wt
Hsum=Hsum+Qmcc*wt
gtot=1./(1./gbair+1./gsplus)
!!! For histograms:
if(histotime)then
indexILnoscatt=ifix(aPAR*D0PAR*skyview/25.)+1
indexILall=ifix(PPFD/25.e-6)+1
indexTL=ifix(Tleafval/0.25)-79 ! Thus, the index=1 for TL=20.
! write(6,'("location ",3i3,"; Tleafval=",f6.2,"; indexTL=",i4,"wt=",f8.6)')&
! &iradius,itheta,iphi,Tleafval,indexTL,wt
! pause
indexA=ifix(Aleaf/0.4e-6)+16 ! Thus, the index=1 for A=-6.e-6
indexE=ifix(Eleaf/0.125e-3)+1
indexgtot=ifix(gtot/0.005)+1
histoILnoscatt(indexILnoscatt)=histoILnoscatt(indexILnoscatt)+wt
histoILall(indexILall)=histoILall(indexILall)+wt
histoTLraw(indexTL)=histoTLraw(indexTL)+wt
histoAraw(indexA)=histoAraw(indexA)+wt
histoEraw(indexE)=histoEraw(indexE)+wt
histogtotraw(indexgtot)=histogtotraw(indexgtot)+wt
! write(6,'("indexTL=",i3,";histoTLraw-->",f10.6)')indexTL,histoTLraw(indexTL)
! pause
histoTL_Awtd(indexTL)=histoTL_Awtd(indexTL)+wt*Aleaf
histoTL_Ewtd(indexTL)=histoTL_Ewtd(indexTL)+wt*Eleaf
histogtot_Awtd(indexgtot)=histogtot_Awtd(indexgtot)+wt*Aleaf
histogtot_Ewtd(indexgtot)=histogtot_Ewtd(indexgtot)+wt*Eleaf
histoIL_Awtd(indexILall)=histoIL_Awtd(indexILall)+wt*Aleaf
histoIL_Ewtd(indexILall)=histoIL_Ewtd(indexILall)+wt*Eleaf
endif
! Il in (0,2500), by steps of 25
! TL in (20., 45.), by steps of 0.25
! A in (-5., 35.)e-6, by steps of 0.4e-6
! E in (0., 12.5.)e-3, by steps of 0.125e-3
! gs in (0., 0.5), by steps of 0.005
!!!
gtotsum=gtotsum+gtot*wt
TLsum=TLsum+gtot*Tleafval*wt
hssum=hssum+hs*wt !DDDDDDDDDDDDDDDDDDDDDDDDDDDD
Cisum=Cisum+gtot*Ci*wt !DDDDDDDDDDDDDDDDDDDDDDDDD
Iadd=ILshade+Pclear*lpen*0.5*I0PAR ! This is intercepted, not absorbed, PAR
! The 0.5 is the average cosine (projection) onto all leaf orientations
sumI=sumI+Iadd
! Note: sumI is then accumulated into Isum, by being passed in common to MAIN
! Now redo with 10% higher gs; copying things from tam.simul.f
gsplus=gsplus*1.1
Tleafval=Tleaf(gsplus)
gtotCO2=1./(1.6/gsplus+1.37/gbair)
call CiACssolve(scalinghere)
Esum10=Esum10+Eleaf*wt
Asum10=Asum10+Aleaf*wt
! Calculate A, E, etc. for clear time
! Calculate contrib. from shade leaf area, of probability 1-lpen
fbin=1./float(nbins)
if(Pclear.gt.0.)then
dI=I0PAR/float(nbins)
do nb=1,nbins
IL=ILshade+dI*(float(nb)-0.5) ! RRRRRRRRRRRRR Add factor aPAR to latter
! Maye rename this ILabs, to be perfectly clear
ENIR=ENIRshade+Wmol*di*(float(nb)-0.5) ! RRRRRRR Add factor aNIR to latter - critical!
QSW=aPAR*IL*Wmol+aNIR*ENIR ! This will be correct heat loading when factor aPAR is added above
Qrad=QTIR+QSW
if(Qrad.gt.2000.)then
write(6,'("sunQrad=",f8.1,"; Qsky,Qveg,QSW=",3f9.1,"; ntimes,iloc,iter=",3i4)')&
&Qrad,Qsky,Qveg,QSW,ntimes,iloc,iter
endif
Tleafval=Tair ! initial guess for leaf T
PPFD=1.e-6*IL ! convert to mol from umol
Aleaf=0.
Eleaf=0.
call BBsolve(scalinghere)
! Calculate weighting function
wt=Pclear*lpen*fbin
Asum=Asum+Aleaf*wt
Esum=Esum+Eleaf*wt
Hsum=Hsum+Qmcc*wt
gtot=1./(1./gbair+1./gsplus)
!!! For histograms:
if(histotime)then
indexILnoscatt=ifix((aPAR*D0PAR*skyview+dI*(float(nb)-0.5))/25.)+1 ! RRRRRR Wrong
indexILall=ifix(PPFD/25.e-6)+1
indexTL=ifix(Tleafval/0.25)-79 ! Thus, the index=1 for TL=20.
indexA=ifix(Aleaf/0.4e-6)+16 ! Thus, the index=1 for A=-6.e-6
indexE=ifix(Eleaf/0.125e-3)+1
indexgtot=ifix(gtot/0.005)+1
histoILnoscatt(indexILnoscatt)=histoILnoscatt(indexILnoscatt)+wt
histoILall(indexILall)=histoILall(indexILall)+wt
histoTLraw(indexTL)=histoTLraw(indexTL)+wt
histoAraw(indexA)=histoAraw(indexA)+wt
histoEraw(indexE)=histoEraw(indexE)+wt
histogtotraw(indexgtot)=histogtotraw(indexgtot)+wt
histoTL_Awtd(indexTL)=histoTL_Awtd(indexTL)+wt*Aleaf
histoTL_Ewtd(indexTL)=histoTL_Ewtd(indexTL)+wt*Eleaf
histogtot_Awtd(indexgtot)=histogtot_Awtd(indexgtot)+wt*Aleaf
histogtot_Ewtd(indexgtot)=histogtot_Ewtd(indexgtot)+wt*Eleaf
histoIL_Awtd(indexILall)=histoIL_Awtd(indexILall)+wt*Aleaf
histoIL_Ewtd(indexILall)=histoIL_Ewtd(indexILall)+wt*Eleaf
endif
! Il in (0,2500), by steps of 25
! TL in (20., 45.), by steps of 0.25
! A in (-5., 35.)e-6, by steps of 0.4e-6
! E in (0., 12.5.)e-3, by steps of 0.125e-3
! gs in (0., 0.5), by steps of 0.005
!!!
gtotsum=gtotsum+gtot*wt
TLsum=TLsum+gtot*Tleafval*wt
hssum=hssum+hs*wt !DDDDDDDDDDDDDDDDDDDDDD
Cisum=Cisum+gtot*Ci*wt !DDDDDDDDDDDDDDDDDD
! !!! Write out Ci/Ca, WUE, and other diagnostics
! Now redo with 10% higher gs; copying things from tam.simul.f
gsplus=gsplus*1.1
Tleafval=Tleaf(gsplus)
gtotCO2=1./(1.6/gsplus+1.37/gbair)
call CiACssolve(scalinghere)
Esum10=Esum10+Eleaf*wt
Asum10=Asum10+Aleaf*wt
enddo
endif
Etree=Etree+Esum
Htree=Htree+Hsum
gtottree=gtottree+gtotsum
TLtree=TLtree+TLsum
hstree=hstree+hssum !DDDDDDDDDDDDDDD
Citree=Citree+Cisum !DDDDDDDDDDDDDDD
!
Etree10=Etree10+Esum10
AAtree10=AAtree10+Asum10
AAtree=AAtree+Asum
enddo !End of loop over canopy azimuths
enddo ! End of loop over canopy zeniths
enddo ! End of loop over canopy radii
Ecalc=Etree
!!! Finish calculating the histograms and write them out
if(histotime)then
scratchname='histo'//infix(1:infixlen)//'.txt'
open(18,file=scratchname)
! open(18,file='histol.txt')
write(18,*)'IL in (0,2500), by steps of 25'
write(18,*)'TL in (20., 45.), by steps of 0.25'
write(18,*)'A in (-5., 35.)e-6, by steps of 0.4e-6'
write(18,*)'E in (0., 12.5.)e-3, by steps of 0.125e-3'
write(18,*)'gs in (0., 0.5), by steps of 0.005'
do ihisto=1,120
! DUPE histoILnoscatt(ihisto)=histoILnoscatt(ihisto)/enlocations
! This is proper normalization, since, at each location (say, in
! sunlit area), wt was
! Cf. how Etree was calculated, as sum over nloc locations
! of Esum, and Esum is a sum at each location of
! Eleaf (E per leaf area) times wt=Pclear*lpen*fbin (for the
! sunlit cases); the wt sums to Pclear*lpen*(sum of fbin=1)
! = Pclear*lpen. With the shade fraction, of
! wt = (1-Pclear)+Pclear*(1-lpen), the sum of wt at any location
! is 1-Pclear + Pclear*(1-lpen) + Pclear*lpen = 1
! This is not fully normalized until values are returned to MAIN:
! Here, we divide Etree by enlocations (call the result ELaavg)
! and then multiply by Leafarea*64.8 (gives total E of tree in
! mol s-1 * 64.8, which converts to L h-1)
! So, in Ecalc, to fully normalize IL, I need to divide by enlocations;
! However, enlocations is not available within Ecalc unless I
! pass it in common, again
histoILnoscatt(ihisto)=histoILnoscatt(ihisto)/enlocations
histoILall(ihisto)=histoILall(ihisto)/enlocations
histoTLraw(ihisto)=histoTLraw(ihisto)/enlocations
histoAraw(ihisto)=histoAraw(ihisto)/enlocations
histoEraw(ihisto)=histoEraw(ihisto)/enlocations
histogtotraw(ihisto)=histogtotraw(ihisto)/enlocations
histoTL_Awtd(ihisto)=histoTL_Awtd(ihisto)/AAtree
! Note that AAtree is the sum of wt*Aleaf over enlocations
histoTL_Ewtd(ihisto)=histoTL_Ewtd(ihisto)/Etree
histogtot_Awtd(ihisto)=histogtot_Awtd(ihisto)/AAtree
histogtot_Ewtd(ihisto)=histogtot_Ewtd(ihisto)/Etree
histoIL_Awtd(ihisto)=histoIL_Awtd(ihisto)/AAtree
histoIL_Ewtd(ihisto)=histoIL_Ewtd(ihisto)/Etree
write(18,'(14f8.4)')histoILnoscatt(ihisto),histoILall(ihisto),histoTLraw(ihisto),&
&histoAraw(ihisto),histoEraw(ihisto),histogtotraw(ihisto),histoTL_Awtd(ihisto),&
&histoTL_Ewtd(ihisto),histogtot_Awtd(ihisto),histogtot_Ewtd(ihisto),&
&histoIL_Awtd(ihisto),histoIL_Ewtd(ihisto)
enddo
close(18)
endif
!!!
! If this is first time for calculating lpen array, and saving is
! requested, then save it
if((ipre.eq.0).and.(ipensave*iter.eq.1))then
! write(10,'("lpen values for ntimes=",i4)')ntimes
do iradius=1,nradius
do itheta=1,ntheta
write(10,*)(lpenstore(iradius,itheta,iphi),iphi=1,nphi)
enddo
enddo
endif
return
end
!
!
real function control(deair,eair)
real deair ! (Pa) test increment in eair in canopy
real eair ! (Pa) current value of eair in canopy
real eairtest ! (Pa) new test value of eair in canopy
! Restrict size of deair step
if(abs(deair).gt.2000.)then
deair=sign(2000.,deair)
endif
! Be sure that eair is not set too low
eairtest=eair+deair
control=amax1(eairtest,200.)
return
end
!
!
real function Tcontrol(dTair,Tair)
real dTair ! (oC) test increment in eair in canopy
real Tair ! (oC) current value of eair in canopy
real Tairtest ! (PoC) new test value of eair in canopy
! Restrict size of deair step
if(abs(dTair).gt.6.)then
dTair=sign(6.,dTair)
endif
! Be sure that Tair is not set too low
Tairtest=Tair+dTair
Tcontrol=amax1(Tairtest,0.)
return
end
!
!
subroutine iterate(eair,Tair,Estand,Hstand,errore,errorT,&
& serrorboth,f,tole,tolT,iconv,Pair,gbcan)
! type declarations of variables
integer i1,i2,i3 ! (-) rank-order of iterations, in sizes of errors
integer istor ! (-) temporary storage of rank order
integer ntimes ! number of the current simulation time (for special reporting
! of intermediate results)
!
real c ! (K m2 s J-1) coefficient in denominator of error derivatives
real deair ! (Pa) estimate of increment eair in canopy, by 2-D extrapolation
real denom ! (Pa K) denominator in error derivatives
real dTair ! (oC or K) estimate of increment in Tair in canopy, by 2-D extrapolation
real Ebar ! (mol m-2 s-1) average transpiration per ground area, over last 3 iterations
real Epe ! (mol m-2 s-1 Pa-1) derivative of transpir. rate per ground area
! with respect to changes in eair in canopy
real EpT ! (mol m-2 s-1 K-1) derivative of transpir. rate per ground area
! with respect to changes in Tair in canopy
real Hbar ! (W m-2) average sensible heat flux density, over last 3 iterations
real Hpe ! (W m-2 Pa-1) derivative of sensible heat flux density
! with respect to changes in eair in canopy
real HpT ! (W m-2 K-1) derivative of sensible heat flux density
! with respect to changes in Tair in canopy
real p ! (Pa mol-1 m2 s) coefficient in denominator of error derivatives
real tole ! (Pa) error tolerance for estimate of eair in canopy
real tolT ! (K) error tolerance for estimate of Tair in canopy
! type declarations of arrays
integer irank(10) ! rank-order of iteration 'i' in size of errors
real eair(0:10) ! (Pa) 'i-th' estimate (iteration) for eair in canopy
real errore(0:10) ! (Pa) error in iteration 'i' for eair in canopy
real errorsort(10) ! (-) normalized error of iteration 'i', while being sorted
! for rank order
real errorT(0:10) ! (oC or K) error in iteration 'i' for Tair in canopy
real Estand(0:10) !
real Hstand(0:10) !
real serrorboth(10) ! (-) normalized error of iteration 'i', accounting for
! both eair and Tair errors (each divided by their respective tolerances)
real tair(0:10) ! (oC) 'i-th' estimate (iteration) for Tair in canopy
!
! type declarations of do-loop indices
integer i,ifind
common/special/ntimes
! Sort the iterations by their combined scaled error, serrorboth, to
! find the three with the least error, closest therefore to the
! correct solution (giving the best extrapolation or interpolation,
! we hope)
do i=1,iconv-1
irank(i)=i
errorsort(i)=serrorboth(i)
enddo
do ifind=1,iconv-2
do i=ifind+1,iconv-1
if(errorsort(i).lt.errorsort(ifind))then
errorsortstor=errorsort(ifind)
errorsort(ifind)=errorsort(i)
errorsort(i)=errorsortstor
istor=irank(ifind)
irank(ifind)=irank(i)
irank(i)=istor
endif
enddo
enddo
! Later: it may be necessary to check that we have sufficient differences
! in eair and in Tair between iterations, so that we are not computing
! derivatives from very small steps that have numerical error associated with them
!
! Now, use the 3 best iterations to solve for the numerical derivatives,
! dE/d(eair), dE/d(Tair), dH/d(eair), dH/d(Tair) (all are actually partial
! derivatives)
! Below, i1, i2, i3 denote the iteration numbers with the lowest, 2nd lowest,
! and 3rd lowest combined errors. We will extrapolate from iteration i1
i1=irank(1)
i2=irank(2)
i3=irank(3)
denom=(eair(i3)-eair(i1))*(Tair(i2)-tair(i1)) &
& - (eair(i2)-eair(i1))*(Tair(i3)-Tair(i1))
EpT=( (eair(i3)-eair(i1))*(Estand(i2)-Estand(i1)) -(eair(i2)-&
& eair(i1))*(Estand(i3)-Estand(i1)) )/denom
Epe=-( (Tair(i3)-Tair(i1))*(Estand(i2)-Estand(i1)) - (Tair(i2)-&
& Tair(i1))*(Estand(i3)-Estand(i1)) )/denom
HpT=( (eair(i3)-eair(i1))*(Hstand(i2)-Hstand(i1)) -(eair(i2)-&
& eair(i1))*(Hstand(i3)-Hstand(i1)) )/denom
Hpe=-( (Tair(i3)-Tair(i1))*(Hstand(i2)-Hstand(i1)) - (Tair(i2)-&
& Tair(i1))*(Hstand(i3)-Hstand(i1)) )/denom
! Compute the new tolerance limits
Ebar=(Estand(i1)+Estand(i2)+Estand(i3))/3.
Hbar=(Hstand(i1)+Hstand(i2)+Hstand(i3))/3.
if(abs(Ebar).lt.0.001)Ebar=0.001 ! We don't want to impose too-strict
! standards on very small E's (have an absolute standard for these,
! f*1 mmol m-2 s-1)
tole=abs(f*Ebar/Epe)
tolT=abs(f*Ebar/EpT)
! if(ntimes.eq.121)write(14,*)f,Ebar,Epe,EpT
! if(ntimes.eq.121)write(14,*)'iconv',iconv,'tole,tolT',tole,tolT
do i=1,iconv-1
serrorboth(i)=sqrt( (errore(i)/tole)**2+(errorT(i)/tolT)**2 )
enddo
! Compute the new values of eair and Tair that should exactly solve the
! transport equations,
! eair = eair(0) + E*Pair/gbcan
! tair = Tair(0) + H/(Cpair*gbcan),
! assuming that E and H are linear functions of eair and Tair,
! E = Estand(i1)+Epe*(eair-eair(i1)) + EpT*(Tair-Tair(i1)), sim. for H
!
c=1./(29.*gbcan)
p=Pair/gbcan
denom=1.-Epe*p-HpT*c + p*c*(HpT*Epe-Hpe*EpT)
! denom=(1.-c*HpT)*(1.-p*Epe)-p*c*EpT*Hpe ! This is the same as above
deair=((Hstand(i1)*c-Tair(i1)+Tair(0))*EpT*p + (-Estand(i1)*p&
& +eair(i1)-eair(0))*(HpT*c -1.))/denom
dTair=((-eair(i1)+eair(0)+Estand(i1)*p)*Hpe*c + (Tair(i1)-Tair(0)-&
& Hstand(i1)*c)*(Epe*p -1.))/denom
eair(iconv)=control(deair,eair(i1) )
Tair(iconv)=Tcontrol(dTair,Tair(i1) )
return
end
!
!
real function psisoilcalc(theta,isoiltype) ! (MPa)
! van Genuchten parameters:
real alpha(11) ! (kPa-1)
real ms ! (-)
real ns ! (-) Value selected from nsdata, acc. to soil type
real nsdata(11) ! (-)
real R ! (-) Intermediate variable, for convenience of calc.
real thetar(11) ! (-) Residual soil water content
real thetas(11) ! (-) Saturated SWC
data thetas/0.375,0.390,0.387,0.399,0.439,0.384,0.442,&
& 0.482,0.481,0.439,0.545/
data thetar/0.053,0.049,0.039,0.061,0.065,0.063,0.079,&
& 0.090,0.111,0.117,0.098/
data alpha/0.35e3,0.35e3,0.27e3,0.11e3,0.05e3,0.21e3,0.16e3,0.08e3,&
& 0.16e3,0.33e3,0.15e3/
data nsdata/3.162,1.738,1.445,1.479,1.660,1.318,1.413,&
& 1.514,1.318,1.202,1.259/
ns=nsdata(isoiltype)
ms=1.-1./ns
R=(thetas(isoiltype)-thetar(isoiltype))/(theta-thetar(isoiltype))
! write(6,'("In psisoilcalc, R=",f8.3,", 1/ms=",f8.3,", 1/ns=",f8.3)')R,1./ms,1./ns
psisoilcalc=-(R**(1./ms)-1.)**(1./ns)/alpha(isoiltype) ! MPa
! pause
return
end
!
!
real function Rsoilcalc(theta,isoiltype) ! (m-1 s-1)
! The flow soil-->root is through a resistance, Rsoil, that depends on
! root length density, soil hydraulic conductivity, etc. The form
! derived elsewhere, is
! Rsoil = C*rhor*rr**2*ln(d/rr) / (2*mr*khPm), where
! C = a root-clumping adjustment (Tardieu's F, set as high as 5)
! rhor = dry matter density of live roots, about 1/4 that of water,
! or about 250 kg m-3
! rr = mean fine-root radius, about 1 mm
! d = mean spacing between fine roots = 1/sqrt(root length density),
! where RLD = (total root system length)/(total soil vol. occupied)
! and (total root system length) =
! (total live-root volume)/(pi*rr**2)
! and (total live-root volume) = (total root dry mass)/rhor
! and (total soil vol. occupied)
! = pi*(canopy radius)**2 * (soil depth)
! mr = total root dry mass
! khPm = soil hydraulic conductivity, in units of SI pressure
! gradient (Pa m-1). It equals khv/g, where khv is the
! conductivity in head units (m per m) for the gradient and
! velocity units (m s-1) for the flow; g is the acceleration of
! gravity.
! I have tabulations of khvs, s = saturated conditions, for various
! soil types
! I computed khv at actual water content as krel*khvs, with the
! van Genuchten form,
! krel=sqrt(R)*[1-(1-R**(1/ms))**ms]**2, and
! R = (theta-thetar)/(thetas-thetar)
real alpha(11) ! (kPa-1)
real Clump ! (-) clumping factor for fine roots, as multiplier of Rsoil
real dFR ! (m) mean distance between fine roots
real khPm ! (s) unsaturated soil hydraulic conductivity, in metric units
! of kg m-2 s-1 per MPa m-1, which is units of seconds
real krel ! (-) ratio of saturated to unsaturated soil hydr. cond.
real khvs(11) ! (m s-1) saturated soil hydraulic conductivity, in
! pressure-head units
real mr ! (kg) total dry mass of fine roots on central tree
real ms ! (-)
real ns ! (-) Value selected from nsdata, acc. to soil type
real nsdata(11) ! (-)
real R ! (-) Intermediate variable, for convenience of calc.
real rhor ! (kg m-3) dry mass density of live fine roots
real rr ! (m) mean radius of fine roots
real thetar(11) ! (-) Residual soil water content
real thetas(11) ! (-) Saturated SWC
common/main2Rsoil/Clump,rhor,rr,dFR,mr
data thetas/0.375,0.390,0.387,0.399,0.439,0.384,0.442,&
& 0.482,0.481,0.439,0.545/
data thetar/0.053,0.049,0.039,0.061,0.065,0.063,0.079,&
& 0.090,0.111,0.117,0.098/
data alpha/0.35e3,0.35e3,0.27e3,0.11e3,0.05e3,0.21e3,0.16e3,0.08e3,&
& 0.16e3,0.33e3,0.15e3/
data nsdata/3.162,1.738,1.445,1.479,1.660,1.318,1.413,&
& 1.514,1.318,1.202,1.259/
data khvs/8.24E-05,4.05E-05,1.23E-05,2.88E-06,1.25E-06,&
& 3.63E-06,7.22E-07,1.94E-07,5.56E-07,3.33E-07,5.56E-07/
R=(theta-thetar(isoiltype))/(thetas(isoiltype)-thetar(isoiltype))
ns=nsdata(isoiltype)
ms=1.-1./ns
krel=sqrt(R)*(1-(1-R**(1./ms))**ms)**2
khPm=khvs(isoiltype)*krel/9.8 ! denominator is g = accel. of gravity
Rsoilcalc=Clump*rhor*rr**2*alog(dFR/rr)*0.5/(mr*khPm) ! m-1 s-1
return
end
!
!
subroutine calc_ASR_EC(iradtype,Ppen,Ppenmin,ASRloc,aleaf,fleaf,reflsoil,&
& D0,I0,atree,Kdiff)
! Subroutine to calculate the absorbed scattered radiation to any location in
! the true ellipsoidal canopy from a unit input of either diffuse or
! direct-beam radiation. It uses as a near equivalent the scattering in
! a uniform layered canopy (ULC) for these calculations, and the appropriate
! optical properties of leaves and soil. See comments in MAIN for a fuller
! explanation
parameter (maxtrees=100, maxsamples=20)
!
! Declarations of simple variables (vs. arrays, later)
integer ihalfenloc ! (-) approximately half the total number of canopy locations
! being sampled - the number we use in regressions to get Keff (below)
integer indexraw ! (-) the index of the layer in the ULC (see iULClayer), before
! reversal of order and correction for bounds
integer iphi ! (-) index for do-loop over values of canopy azimuthal angle
integer iradius ! (-) index for do-loop over values of canopy radius
integer iradtype ! (-) 1 = diffuse only; 2 = direct-beam (and possibly
! diffuse simultaneously)
integer itheta ! (-) index for do-lopp over values of canopy zenith angle
integer iULClayer ! index of the layer in the ULC corresponding to apparentL,
! the apparent location in the ULC as leaf area index
integer locpack ! (-) holder for current value of ilocpack
integer nloc ! (-) total number of locations being sampled in the canopy
integer nphi ! (-) number of distinct azimuthal angles sampled in the canopy
integer nradius ! (-) number of distinct radii sampled in the canopy
integer ntheta ! (-) number of distinct zenith angles sampled in the canopy
real apparentL ! (-) apparent value of leaf area index at which current point
! would be found in a ULC
real atree ! (m) major axis of the canopy for the central tree
real D0 ! (variously -, mol m-2 s-1, or W m-2, according to place called from)
! flux density of diffuse skylight incident at top of canopy
real enloc ! (-) nloc converted to real, for faster calculations
real fd ! (m2 m-3 = m-1) foliage area density in canopy
real I0 ! (variously -, mol m-2 s-1, or W m-2, according to place called from)
! flux density of direct solar beam incident at top of canopy
real invKL ! (-) inverse of Keff*Lavg, as an intermediate variable
real Kdiff ! (-) extinction coefficient for diffuse skylight in equivalent ULC
real Keff ! (-) effective extinction coefficient for specified type of
! radiation (diffuse or direct)
real kpen ! (-) penetration coefficient per unit leaf area index (0.5 is
! appropriate for a random leaf angle distribution
real Lavg ! (m2 m-2 = -) effective optical depth of ULC, as leaf area index
real lnPpenmin ! (-) natural log of minimal Ppen we will resolve
real Ppenmin ! (-) min. value of Ppen worth tracking in rank-ordering of Ppen's
real slope ! (-) slope of the regression of ln(Ppen) on cumulative
! fraction of all ln(Ppen) values
!
! Declarations of arrays
integer ilocpack(8000) ! stores indices of canopy location (radius, 2 angles)
! in a packed form; follows the rank-reordering of Ppen, so we can reassociate
! ASR to proper locations in canopy
real ASR(11200) ! (-) absorbed scattered radiation generated in the ULC at
! (rank-ordered) depth as LAI, for a unit input of radiation
! Note that ASR will be calculated for L > Lavg used in the ULC simulation -
! thus, for 1.4*nloc points, or 175 in the current case
real ASRloc(maxsamples,maxsamples,maxsamples) ! (-) Same as ASR, but mapped to locations in the canopy
! specified by indices of radius, zenith angle, and azimuthal angle
real CF(8000) ! (-) cumulative fraction of all canopy locations at which the
! value of Ppen (or ln of same) falls when these are rank-ordered
real lnPpen(8000) ! (-) natural log of Ppen (see below)
real Ppen(maxsamples,maxsamples,maxsamples) ! (-) probability of penetration of radiation to canopy
! location specified by indices of radii and two angles
! Dimensions above are set by the (unlikely) use of up to 20 samples each in
! canopy radius, zenith angle, and azimuthal angle
!
! Declarations below are for variables passed in common but not used in
! any calculations in this subroutine
integer ntrees
real costheta(maxsamples), sintheta(maxsamples), r(maxsamples,maxsamples)
real phi(maxsamples), xo, yo, zo, atree_c(maxtrees), btree(maxtrees)
real thetac(maxtrees), phic(maxtrees), xtree(maxtrees)
real ytree(maxtrees), ztree(maxtrees), hcan
common/centertreedescrip/nradius,ntheta,nphi,costheta,sintheta,&
& r,phi,xo,yo,zo
common/alltrees/ntrees,atree_c,btree,thetac,phic,xtree,ytree,ztree,&
& hcan,kpen,fd
sleaf=1.-aleaf-fleaf ! backward scattering
! If iradtype=2, then Ppen was available inside the subroutine (Ecalc, that is)
! Get ln(Ppen) and make associations of the Ppen values with canopy locations
! write(6,'("nradius, ntheta, nphi=",3i4)')nradius,ntheta,nphi
! pause
! write(6,'("ASR:for aleaf,fleaf,reflsoil,D0,I0=",5f8.3)')aleaf,fleaf,reflsoil,D0,I0
iloc=0
do iradius=1,nradius
do itheta=1,ntheta
do iphi=1,nphi
iloc=iloc+1
ilocpack(iloc)=10000*iradius+100*itheta+iphi
! write(6,'("iloc=",i6,"ilocpack(iloc)=",i8)')iloc,ilocpack(iloc)
lnPpen(iloc)=log(Ppen(iradius,itheta,iphi))
enddo
enddo
enddo
nloc=iloc
enloc=float(nloc)
ihalfenloc=ifix(0.5*enloc)
! write(6,'(5f8.4)')(lnPpen(i),i=1,nloc)
! pause
lnPpenmin=log(Ppenmin)
! Rank-order the Ppen values; sort ilocpack at the same time to keep track of
! correspondences
call ssort(lnPpen,ilocpack,nloc)
! write(6,'(5f8.4)')(lnPpen(i),i=1,nloc)
! Before fitting the ln(Ppen) to exp(-0.5*Keff*CF*Lavg), purge truncated values;
! also, create cum. frac. values for each value
do i=1,nloc
CF(i)=(float(nloc-i)+0.5)/enloc ! starts at 0.5/enloc, steps of 1/enloc -
! e.g., for standard nloc=125, starts at 0.996, then 0.988, ..., to 0.004
if(lnPpen(i).le.lnPpenmin)goto 100
enddo
100 istop=i-1 ! Subtract 1 because do-loops increment register to 1 more than end
istop=min(istop,ihalfenloc) ! In any case, don't use the lowest half of the Ppen values,
! which can come from the chords of L > Lavg (up to 1/3 longer)
! write(6,'("istop=",i4)')istop
! Do lsq fit of ln(Ppen) to CF
call lsqfit(lnPpen,CF,istop,slope,aint)
! write(6,'("slope=",f8.4)')slope
Lavg=fd*2.*atree/3.
Keff=slope/Lavg
! write(6,'("Keff=",f7.3)')Keff
if(iradtype.eq.1)Kdiff=Keff ! Return this value for use in later calls
! Recall associations of location with Ppen --> It checks out
do ilocsorted=1,nloc
locpack=ilocpack(ilocsorted)
iradius=locpack/10000
itheta=(locpack-10000*iradius)/100
iphi=locpack-iradius*10000-itheta*100
! write(6,'("ilocsorted=",i4,", locpack=",i12,", iradius, itheta, iphi=",3i10)')&
! & ilocsorted,locpack,iradius,itheta,iphi
enddo
! Now run the ULC simulation to get the absorbed scattered radiation at end location
! First, set the optical properties of the leaves (leaf clusters of random
! orientation)
call calc_ASR_ULC(Lavg,nloc,Keff,Kdiff,aleaf,fleaf,sleaf,reflsoil,&
& D0,I0,ASR)
!
! Finally, at each canopy location, apply the ASR computed from the ULC.
! Run over all the ln(Ppen) values and convert them to equivalent depths, L, in the ULC.
invKL=1./Keff*Lavg
CFmax=1.-0.5/float(nloc) ! This should be CF for the highest point
do iloc=1,nloc
! Retrieve the geometric location for this index
locpack=ilocpack(iloc)
iradius=locpack/10000
itheta=(locpack-10000*iradius)/100
iphi=locpack-iradius*10000-itheta*100
! L = -(1/Keff)ln(Ppen) --> equivalent to CF=L/Lavg=-(1/[Keff*Lavg])ln(Ppen);
! this CF is not the original CF, but close to the latter, which hues to the
! line fitted with least squares
! Compute the index, iCF, of the closest CF segment that does not exceed the
! effective value of L
apparentL=-lnPpen(iloc)/Keff
indexraw=125-ifix(125.*apparentL/Lavg)
if(iradtype.ne.1)then
! write(6,'("ln(Ppen)=",f8.3,", apparent L=",f8.3,", indexraw=",i4)'),&
! & lnPpen(iloc),apparentL,indexraw
endif
iULClayer=125-indexraw
iULClayer=max(iULClayer,1)! Some points may lie
! above the regression line, as if L <0; use max0 to prevent this
! If icFnew is less than zero, it means that Ppen was below exp(-Keff*Lavg);
! we have to use the extrapolated ASR fluxes.
! Remember, CF=1.0 (or 1-0.5/nloc) is n=1 and L = 0 (or 0.5*Lavg/nloc)
! and CF=0.0 (really, 0.5/nloc) is n=nloc and L = Lavg (or short of it by
! fraction 0.5/nloc
! So, we want to retrieve ASR corresponding to the correct index: high effective
! L --> low Ppen --> very neg. ln(Ppen) --> high index n is ASR(n)
! We can go to n as large as 1.4*nloc, for which we're calculated ASR(n) as if
! by extrapolation, but no larger
iULClayer=min(iULClayer,175) ! When pasting this code into canopy_simXX.f90,
! we'll have to have the max index (here, 175) calculated in MAIN and
! passed to the subroutine that this program becomes
! Another special case: if the ray's transit pass through the canopy exceeded stmax
! and we stopped tracking it (Ppen is only 0.002 then), we need to use
! a value appropriate to this low Ppen. As many at 30% of the points can come
! out with this same low value of Ppen, and we can't rank-order them, really,
! We'll just use the apparent index iCFnew. No extra processing of iCFnew is
! needed, therefore
! Store the estimated ASR value according to indices that locate the point in
! the ellipsoidal canopy
! if(iradtype.ne.1)then
! write(6,'("iloc=",i4,", iradius,itheta,iphi=",3i5,",lnPpen=",f8.3,&
! & ", iULClayer=",i4)')iloc,iradius,itheta,iphi,lnPpen(iloc),&
! & iULClayer
! pause
! endif
ASRloc(iradius,itheta,iphi)=ASR(iULClayer)
! write(6,'(5x," ASRloc=",f8.3)')ASR(iULClayer)
! if(mod(iloc,1).eq.1)pause
enddo
return
end
!
!
subroutine calc_ASR_ULC(xf,nlayers,Kdir,Kdiff,&
& aleaf,fleaf,sleaf,r,D0,I0,ASR)
! Simulates diffuse (scattered) fluxes and their absorption in a
! uniform layered canopy (ULC). Modified to be a subroutine, from tested program
! canopy_abs+scatt.f90 ~/fortran/radtpt 21-22 Feb. 09
!
! Declaration of simple variables
! Variables without explanations are algebraic intermediates, whose definition
! requires access to the fairly lengthy exposition of the radiation-
! propagation model, which cannot be included here for lack of space
integer nmore ! (-) number of layers beyond average LAI for which the
! absorbed scattered radiation is calculated as an extrapolation
real a ! (-) absorption of diffuse radiation (either upwelling or downwelling)
! per unit LAI traversed
real abs ! (-, or mol m-2 s-1, or W m-2; varies, according to calling routine)
! running total, over all canopy layers, of absx, defined below
real abssoil ! (-, or mol m-2 s-1, or W m-2; varies, according to calling routine)
! flux density absorbed by ground or soil
real absx ! (-, or mol m-2 s-1, or W m-2; varies, according to calling routine)
! rate of absorption of flux density, per unit LAI, at this location in canopy
real adir ! (-) absorption of direct beam per unit LAI traversed
real aleaf ! (-) absorptivity of leaf
real Aminus0
real Aminusxf
real Aplus0
real Aplusxf
real bminus
real bplus
real D ! (-, or mol m-2 s-1, or W m-2; varies, according to calling routine)
! value of DOWNwelling diffuse flux density at current location
! (LAI value) in the canopy
real D0 ! (as with D) DOWNwelling diffuse flux density incident at
! top of canopy
real denom
real dx ! (-) increment in LAI per layer in the canopy
real Dxf
real Eminus
real enlayers ! (-) number of uniform layers (increments in LAI) being
! simulated in the canopy
real Eplus
real fac
real fdir ! (-) forward scattering of direct beam, per unit of leaf area
! index traversed
real fleaf ! (-) forward scattering (transmissivity) of direct beam by leaf
real I ! (-, or mol m-2 s-1, or W m-2; varies, according to calling routine)
! value of direct-beam flux density at current location
! (LAI value) in the canopy
real I0 ! (as with I) direct-beam flux density incident at top of canopy
real K
real Kdiff ! (-) extinction coefficient for diffuse flux, per unit LAI
real Kdir ! (-) extinction coefficient for direct-beam flux, per unit LAI
real ks
real P1, P2, P3
real part1
real ratio
real s ! (-) scattering of diffuse radiation (either upwelling or downwelling)
! per unit LAI traversed
real sleaf ! (-) (back)scattering of direct beam radiation by leaf
real sdir ! (-) (back)scattering of direct beam, per unit of leaf area index
! traversed
real sumabs_extra_D ! (-, or mol m-2 s-1, or W m-2; varies, according to
! calling routine) sum, over all canopy locations, of the extra DOWNwelling
! flux density generated by scattering (2nd and higher interceptions)
real sumabs_U ! (-, or mol m-2 s-1, or W m-2; varies, according to
! calling routine) sum, over all canopy locations, of the extra UPwelling
! flux density generated by scattering (2nd and higher interceptions)
real sum_asr ! (-) sum, over all canopy locations, of the flux density
! (irradiance) absorbed by leaves, in excess of first interceptions -
! that is, absorbed scattered radiation = ASR
real Totalin ! (-) total radiant flux density incident at top of canopy, as
! sum of direct beam and diffuse skylight
real Totalout ! (-) sum of fates of flux density, as absorbed (=irradiance,
! to be proper; absorption by both canopy and ground) and reflected from top
real U ! (-, or mol m-2 s-1, or W m-2; varies, according to calling routine)
! value of UPwelling diffuse flux density at current location
! (LAI value) in the canopy
real U0 ! (as with U) UPwelling diffuse flux density exiting
! top of canopy
real x ! (-) current position in canopy, as cumulative LAI from top
!
! Declaration of arrays
real ASR(*) ! (-, or mol m-2 s-1, or W m-2; varies, according to calling routine)
! flux density (irradiance) absorbed by leaves in n-th canopy layer, in excess
! of first interceptions - that is, absorbed scattered radiation = ASR
! real CF(8000) ! No longer needed
! write(6,'("calc_ASR is using Kdiff=",f7.3)')Kdiff
!
a=aleaf*Kdiff
s=sleaf*Kdiff
adir=aleaf*Kdir
sdir=sleaf*Kdir
fdir=fleaf*Kdir
part1=-(1.+a/s)
K=sqrt(2.*a*s+a*a)
bplus=part1+K/s
bminus=part1-K/s
ks=adir+sdir+fdir
Eplus=(exp(K*xf)-exp(-ks*xf))/(K+ks)
Eminus=(exp(-ks*xf)-exp(-K*xf))/(K-ks)
fac=exp(-K*xf)
P1=(fdir-bminus*sdir)*I0*Eplus*fac
ratio=(1.+r*bminus)/(1.+r*bplus)
P2=(fdir-bplus*sdir)*I0*Eminus*ratio*fac
P3=(bminus-bplus)*r*I0*exp(-ks*xf)*fac/(1.+r*bplus)
denom=bminus-bplus*exp(-2.*K*xf)*(1.+r*bminus)/(1.+r*bplus)
! write(6,'("K,bplus,bminus,ks,Eplus,Eminus=",/,2x,6f9.3)')&
! & K,bplus,bminus,ks,Eplus,Eminus
! write(6,'("fac,P1,ratio,P2,P3,denom=",/,2x,6f10.3)')&
! & fac,P1,ratio,P2,P3,denom
Aplus0=((bminus-bplus)*D0+bplus*(P2-P1+P3))/(bminus-bplus*fac*fac*ratio)
Aminus0=(Aplus0*bminus-(bminus-bplus)*D0)/bplus
U0=(Aminus0-Aplus0)/(bminus-bplus)
! write(6,'("Aplus0=",f8.3,", Aminus0=",f8.3,", U0=",f8.3)')&
! & Aplus0,Aminus0,U0
dx=xf/float(nlayers)
abs=0.
sum_asr=0.
sumabs_extra_D=0.
sumabs_U=0.
x=-0.5*dx
do n=1,nlayers
x=x+dx
! CF(n)=x/xf
Aplus=Aplus0*exp(-K*x)+(fdir-bplus*sdir)*I0*(exp(-ks*x)-exp(-K*x))/(K-ks)
Aminus=Aminus0*exp(K*x)+(fdir-bminus*sdir)*I0*(exp(K*x)-exp(-ks*x))/(K+ks)
D=(bminus*Aplus-bplus*Aminus)/(bminus-bplus)
U=(Aminus-Aplus)/(bminus-bplus)
I=I0*exp(-ks*x)
absx=a*(D+U)+adir*I
abs=abs+absx
abs_extra_D=a*(D-D0*exp(-Kdiff*x))
abs_U=a*U
sumabs_extra_D=sumabs_extra_D+abs_extra_D
sumabs_U=sumabs_U+abs_U
ASR(n)=a*(D-D0*exp(-Kdiff*x)+U) ! Removes the first interceptions, already accounted
! if(mod(n,10).eq.1)write(6,'(" x=",f6.3,": Aplus=",f8.3,", Aminus=",f8.3,/,&
! & 4x,"D=",f8.3,", U=",f8.3,", I=",f8.3,", absrate=",f8.3,/,&
! & 8x,"abs_extra_D=",f8.3,", abs_U=",f8.3,", total ASR=",f8.3)')&
! & x,Aplus,Aminus,D,U,I,absx,abs_extra_D,abs_U,ASR(n)
asr_sum=asr_sum+ASR(n)
enddo
! Print results for inspection
! do n=1,nlayers,10
! write(6,'("Layer #",i3,": ASR=",f7.4)')n,ASR(n)
! enddo
! write(6,'("Sum of ASR=",f7.3,", from extra D=",f7.3,", from U=",f7.3)')&
! & asr_sum*dx,sumabs_extra_D*dx,sumabs_U*dx
! Accounting
abs_can=abs*dx
Refl_toc=U0
Aplusxf=Aplus0*fac+(fdir-bplus*sdir)*I0*Eminus
Aminusxf=Aminus0*exp(K*xf)+(fdir-bminus*sdir)*I0*Eplus
Dxf=(bminus*Aplusxf-bplus*Aminusxf)/(bminus-bplus)
! write(6,'("Aplusxf=",f8.3,", Aminusxf=",f8.3,", Dxf=",f8.3)')&
! & Aplusxf,Aminusxf,Dxf
abssoil=(Dxf+I0*exp(-ks*xf))*(1.-r)
Totalout=Refl_toc+abs+abssoil
Totalin=I0+D0
! write(6,'("Input: I0+D0 = ",f8.3)')Totalin
! write(6,'("Outputs:",/,3x,"Refl_toc=",f8.3,/,3x,"abs_can=",f8.3,&
! & /,3x,"abs_soil=",f8.3,/,"Total outputs=",f8.3)')Refl_toc,abs_can,abssoil,&
! & Refl_toc+abs_can+abssoil
! Calculate ASR(n) for n > nloc, as an extrapolation scheme, simply by using the
! analytical forms of D(n), U(n)
x=xf-0.5*dx
! Do 40% more points
enlayers=float(nlayers)
nmore=ifix(0.4*enlayers)
! write(6,'("Doing layers of higher indices, from ",i3," to ",i3)')nlayers+1,&
! & nlayers+nmore
! pause
do n=nlayers+1,nlayers+nmore
x=x+dx
! CF(n)=x/xf ! We don't need to retain this marker for L > Lavg
Aplus=Aplus0*exp(-K*x)+(fdir-bplus*sdir)*I0*(exp(-ks*x)-exp(-K*x))/(K-ks)
Aminus=Aminus0*exp(K*x)+(fdir-bminus*sdir)*I0*(exp(K*x)-exp(-ks*x))/(K+ks)
D=(bminus*Aplus-bplus*Aminus)/(bminus-bplus)
U=(Aminus-Aplus)/(bminus-bplus)
I=I0*exp(-ks*x)
absx=a*(D+U)+adir*I
abs=abs+absx
abs_extra_D=a*(D-D0*exp(-Kdiff*x))
abs_U=a*U
sumabs_extra_D=sumabs_extra_D+abs_extra_D
sumabs_U=sumabs_U+abs_U
ASR(n)=a*(D-D0*exp(-Kdiff*x)+U) ! Removes the first interceptions, already accounted
! if(mod(n,10).eq.1)write(6,'(" x=",f6.3,"; abs_extra_D=",f8.3,", abs_U=",f8.3,&
! & ", total ASR=",f8.3)')x,abs_extra_D,abs_U,ASR(n)
enddo
return
end
!
!
subroutine lsqfit(y,x,istop,slope,aint)
! Does (uniformly weighted) least-squares fit of y to slope*x + aint
! Declaration of simple variables
integer istop ! number of points used
integer n ! do-loop index
real aint ! intercept in least-squares fit
real denom ! denominator in expressions for both slope and intercept
real enptsused ! number of points, as a real number
real slope ! obvious
real x2sum ! sum of values of x(i)**2
real xsum ! sum of values of x(i)
real xysum ! sum of values of x(i)*y(i)
real ysum ! sum of values of y(i)
! Declaration of arrays
real x(*) ! array of x values
real y(*) ! array of y values
xsum=0.
ysum=0.
x2sum=0.
xysum=0.
enptsused=float(istop)
do n=1,istop
xsum=xsum+x(n) ! extra fetches of x, y, but no big deal
ysum=ysum+y(n)
x2sum=x2sum+x(n)*x(n)
xysum=xysum+x(n)*y(n)
enddo
denom=enptsused*x2sum-xsum*xsum
slope=(enptsused*xysum-xsum*ysum)/denom
aint=(x2sum*ysum-xsum*xysum)/denom
return
end
!
!
SUBROUTINE SSORT (X, IY, N)
IMPLICIT NONE
! From http://www.personal.psu.edu/jhm/f90/progref.html
!
! Example of a Bubble Sort
!
!***BEGIN PROLOGUE SSORT
!***PURPOSE Sort an array and make the same interchanges in
! an auxiliary array. The array is sorted in
! decreasing order.
!***TYPE SINGLE PRECISION
!***KEYWORDS SORT, SORTING
!
! Description of Parameters
! X - array of values to be sorted (usually abscissas)
! IY - array to be carried with X (all swaps of X elements are
! matched in IY . After the sort IY(J) contains the original
! postition of the value X(J) in the unsorted X array.
! N - number of values in array X to be sorted
!
!***REVISION HISTORY (YYMMDD)
! 950310 DATE WRITTEN
! John Mahaffy
!***END PROLOGUE SSORT
! .. Scalar Arguments ..
INTEGER N
! .. Array Arguments ..
REAL X(*)
INTEGER IY(*)
! .. Local Scalars ..
REAL TEMP
INTEGER I, J, JMAX, ITEMP
! .. External Subroutines ..
! None
! .. Intrinsic Functions ..
! None
!
!***FIRST EXECUTABLE STATEMENT SSORT
!
JMAX=N-1
DO 200 I=1,N-1
TEMP=1.E38
DO 100 J=1,JMAX
IF(X(J).GT.X(J+1)) GO TO 100
TEMP=X(J)
X(J)=X(J+1)
X(J+1)=TEMP
ITEMP=IY(J)
IY(J)=IY(J+1)
IY(J+1)=ITEMP
100 CONTINUE
IF(TEMP.EQ.1.E38) GO TO 300
JMAX=JMAX-1
200 CONTINUE
300 RETURN
END