CODE OUTLINE

Obtain Fixed Input Parameters

• Open file c6.filist [360]
  • Read # of requested runs and associated input and output files

• Read the input data
  • fixed data (infix)
  • soil data (insoil)
  • sprinkler data (inspkl)
  • plant data (inplnt)
  • initial soil water content profile (initc)

Begin daily loop [580-2850]:

• Calculate RLEAF, TLEAF and ALEAF based on current vegetation status
  • takes into account live and dead leaf fraction

• Calculate R/TLFHEM(K,ANGLE) - leaf hemispherical reflectance and transmittance (Starks et al, 1991) [630]
  • Here, ANGLE is the angle between the leaf normal and the sun (???)

• Calculate R/TLFDIF(K), which is R/TLFHEM integrated over the hemisphere (over all ANGLEs). [650]

• Read in the hourly weather data (hour index = IHR): [680]
  • Input variables:
    • IYEAR Calendar year of data
    • ICUMDY Cummulative day (Jan 1 = 0)
    • TIMLOC(IHR) Local time
    • WIND(IHR) Wind speed (m/s)
    • RADTOP(K,IHR) Radiation at top of canopy (W/m2)
    • FBEAM1(K,IHR) Fraction of direct beam radiation/total radiation at wavelengths k, measured above canopy
    • TEMAIR(IHR) Air temperature (C)
    • VPAIR(IHR) Air vapor pressure (mbar) at reference height above the canopy
    • PRECIP(IHR) Precipitation (mm) for hourly increment
    • IRRCHK(IHR) Flag to specify type of precipitation.
  • Opportunity to correct TEMAIR, VPAIR and WIND, if necessary [720]
    • E.g., FIFE measurements were made at non-standard reference height (?) so we can calculate an input wind profile correction factor [710]
  • Calculate vapor pressure deficit (VPDIN) and relative humidity (RELHUM)
  • Continue reading until current day does not equal previous day

• Evaluate day, month and year quantities, if they have not been set (date)

• Find the declination of the sun for each hour (declin)

• Determine sun zenith angle (zenith)
  • also determines first (STRT) and last (END) daylight hours of day

• Set NIR and VIS radiation prior to sunrise and after sunset to 0 [790]

• For each hour, calculate the distribution of leaf-normals to sun angle (DISTLS(IANGL,IHR)). (dstlit)
  • This also returns the weighted mean sunlit leaf reflectances (RLFDIR) and transmittances (TLFDIR) (weighted by sunlit LAI).

• If the fractions direct/diffuse for PAR and NIR haven't been specified in the input, calculate them (from the total solar radiation, in this case). (radin4) [870]

• Calculate the thermal sky emission, if it wasn't specified (skyir)

• Calculate a factor (PATH(IHR)) used to determine the pathlength of a direct sunbeam through the canopy for each hour (calcpath):
  • PATH * H = pathlength through canopy

• If the canopy is shaded by a hill at a given hour (based on PATH), remove the direct beam radiation from the PAR and NIR (leave the thermal alone) [910]

• Recalculate RADTOP and FBEAM1 values, taking into account the fraction of the hemisphere that is visible for a given slope and aspect. [960]

• Given the specified minimum value for LAI contained in one layer (DFMIN), calculate the total number of layers (JTOT).
  • This is constrained to lie between JMIN and JMAX.

• Calculate interception factors for diffuse radiation in a canopy layer (difint).

• Determine the function z/h(LAI) - the fractional canopy height as a function of LAI increment. (height)
  • Something unusual can happen at CUMLAI = TOTLAI/2 if non-uniform distribution was specified in the input file.

• Calculate wind factors (FWIND(J)) for determining the wind speeds within the canopy. (wndfac)

• Fill in the height arrays with physical heights (hite2):
  • Z(JZ)
  • ZMID(JZ) = [Z(JZ) + Z(JZ+1)]/2
  • DELTAZ(JZ) = [Z(JZ+1) - Z(JZ-1)]/2

• Calculate the fraction of the root system which extracts water (FROOT) and the net hydraulic root resistance (RESROT) at each soil node. [1050]
  • ROOTSM = sum-over-soil-layers(1/RESROT)

• If JTOT has changed, JZ indexing system is screwed up. Adjust indices for last hours soil parameters accordingly so comparison is proper. [1070]

Begin hourly loop [1213-2550]:

• Set up first guess for this hour for soil water content, leaf temperature, air temp and air vapor pressure:
  • Departures of TAIR-TEMAIR and TEMPLF-TEMAIR and EAIR-VPAIR set to same magnitude as was found last hour.

• Execute the canopy radiation model (radiat)

• Calculate the leaf boundary layer resistance to heat flow (see Goudriaan 1977, Eq 3.5) (rbound).

• Execute the appropriate photosynthesis model:
  • Farquhar model (photks)
  • Berry C4 plants (c4phot)
  • Berry C3 plants (c3phot)
  • (The stomatal conductance determined here is needed to calculate the leaf energy budget (transpiration rate)).

• Calculate the angle of incidence for precipitation droplets at the top of the canopy (preang)
  • This will depend on size of the drop (ie thru its terminal velocity) and on the wind speed.

• Calculate canopy layer precipitation interception weighting function (cpywt)
  • A layer is more likely to intercept drops if they are coming in at a more horizontal angle - the pathlength through the layer is longer.
  • The amount of precip decreases as you move down through canopy due to interception by higher layers.

• Determine the amount of precipitation intercepted by each canopy layer (interc).

• Run the leaf energy balance routine (lfebal). [1410]

• Based on the ET returned by lfebal, calculate the amount of intercepted precipitation that runs down the stems, drips, or stays on the leaves (drips).

• If the leaves were found by drips to be drying, call lfebal again to include this effect in the energy budget (using new FRWET from drips).

Begin Main Convergence Loop [1420-1770]:

"Solving the leaf energy budget, which involves radiative, sensible and latent heat exchanges, requires a knowledge of the radiative, temperature, humidity and wind environments adjacent to the leaf. Environmental conditions adjacent to leaves are obtained from one-dimensional radiation (deWit 1965) and turbulent transfer equations (Goudriaan 1977) which in turn need the collective results of the leaf energy budgets for all the leaves. Therefore leaf energy budgets and vertical profile equations must be solved simultaneously." (Norman 1989)

Begin PSI Convergence Loop [1420-1570]:

Begin Leaf Temperature Convergence Loop [1430-1480]:

• Move TEMPLF (from lfebal) 10% closer to TEMLF1 (input to radiat). [1433]
• Set TEMLF1 = TEMPLF.
• Run the canopy radiation model (radiat) (input temp is stored in TEMLF1).
  • Temp loop input: TEMPLF (lfebal)
• Run the appropriate photosynthesis model.
  • Temp loop input: TEMPLF (lfebal), DSTRAD (radiat)
• Get the intercepted precipitation (interc).
• Run the leaf energy budget routine (lfebal) and get a new TEMPLF.
  • Temp loop input: RSLEAF (photosyn model), FRWET & PINT1 (interc), DSTNET, DSTRAD, FRAREA (radiat)
• Determine the drippage (drips).
  • Temp loop input: EVINT (lfebal), PINT (interc)
• If leaves are drying, call lfebal again.
  • Temp loop input: Same as above + new FRWET from drips and DSTNET and DSTRAD thermal components have been updated by last call to lfebal
• If, for any leaf angle class or layer, abs(TEMPLF - TEMLF1) > 0.2, recycle the leaf temp conv. loop.

End Leaf Temperature Convergence Loop

• Calculate the water extraction rate by root system at each soil level and the water potential at the top (???) of the canopy (rootex)
  • ROOTUP(JZ) = (PN(JZ) - PSIXY)/RESROT(JZ) where PN/PSIXY are the soil/xylem water potentials (Ohm's Law for water).
  • PSITOP = PSIXY - [PSIXY->TOP] = PSIXY - CPYTR*RROOT/6 (assuming plant resistance is RROOT/6).
  • PSI Loop input: TRAN (lfebal), FRAREA (radiat)

• Calculate a bunch of fluxes per unit leaf area (sensible heat, ET, etc) (sumflx). [1490]
  • PSI Loop input: GEVAP, GHEAT, GHWATR (lfebal)

• If the number of iterations through the PSI loop exceeds 50 (ITERH>50), escape this loop. [1500]

• If this is our first iteration through the subloop for this hour (ITERH=0):
  • Do some things here including setting PSTIP1=PSITOP.
  • We've just been setting some preliminary estimates. Go back to beginning of the PSI loop.

• If at this point abs(PSITOP-PSITP1) < 0.1, we can escape the PSI loop. [1525]

• Increment PSITOP (NR technique?) [1545]

• After 10th iteration of subloop (ITERH > 10), begin averaging current PSITOP with value from last iteration. [1550]

• Store PSI and TR values:
  • PSITOP->PSITP1->PSISV1
  • CPYTR->TRSAV1

• PSI hasn't stabilized yet, so repeat PSI convergence loop .

End PSI Convergence Loop

• Calculate soil and canopy profiles (profl2). [1580]

• Evaluate the stability of the air vapor pressure and temperature estimates:
  • if [[abs(TAIR1-TAIR) < 0.02] OR
  • if [ ITER2>40 && abs(TAIR1-TAIR) < 0.2]]
  • AND
  • if [[abs(EAIR1-EAIR) < 0.02] OR
  • if [ ITER2>40 && abs(EAIR1-EAIR) < 0.2]]
  • then
    • We've converged! Escape this convergence loop!
  • else
    • Average various quantities with previous estimates to hasten convergence. [1700]
    • Start main convergence loop all over again.

End Main Convergence Loop

• If condensation has been found in a given layer for this hour, increment that layer's leaf wetness duration. [1780]

• A few steps here [1800-2010].

• Calculate the increment of plant dry matter for the given photosynthetic rate (drywt). [2010]

• Do some conversions and write some hourly output to the output file [2025].

• Scale hourly leaf quantities to canopy measures (???) (scaleh).

• Calculate bi-directional thermal reflectance function for the canopy (bdrtm). [2450]

• Calculate the hourly evapotranspiration rate (penmon). [2515]

• Verify that the temperature and VP estimates do not exceed their maximum limits in any layer. [2526]

End Hourly Loop (go to next hour)

• Check for energy balance over the whole system. [2610]

• Scale daily parameters (scaled). [2820]

End Daily Loop (go to next day)