The MAESTRA Model

This website maintained by Belinda Medlyn: b.medlyn@unsw.edu.au
Last updated Sept. 01
Hosted by Biological Science at University of NSW
Mirror sites:
http://www.ed.ac.uk/~bmedlyn/maestra/
http://www.gwdg.de/~aibrom/maestra/

Some history:

MAESTRO/MAESTRA is a model of forest canopy radiation absorption and photosynthesis. The model has a long history, going back to the work of Norman & Jarvis in the 1970's and 80's. Ying Ping Wang improved and tested the model for his PhD thesis, and it was published in Wang & Jarvis (1990) Agric For Meteorol 51:257-280. A lot of other people have worked on the model over the years, and as of 1997 there were several versions of the model in existence, most of which were complicated and difficult to understand or modify. For the purposes of the ECOCRAFT project, I have revised a version of the MAESTRO model which I obtained in early 1997 from Ying-Ping Wang. The revised model was renamed to MAESTRA. The major revisions included: (i) removing redundant and ‘spaghetti’ code; (ii) modularising the code to make the program easier to understand and modify; (iii) incorporating standard formulations of physiological sub-modules (Farquhar - von Caemmerer, Ball-Berry).

Model description:

The model is implemented in DOS-based Fortran, with text input and output files. At present there are no plans to develop a Windows version.

Would you like to be kept informed of changes to the MAESTRA program? If you write me a note I will add your name to my mailing list. This mailing list will be advised whenever an updated version of MAESTRA becomes available from these pages. Please note that all ECOCRAFT members are automatically on my mailing list!

Met data file. The start of the actual met data in the file met.dat must be preceded by the line
DATA STARTS
(note: caps essential; no tabs or extra blank spaces). This ensures the file is read from the correct place, regardless of system.

Soil water content and its effect on stomatal conductance. Soil water content may be specified in the met.dat file by adding a column 'SW'. This will be converted into a soil moisture deficit using the two additional parameters SWMAX and SWMIN which must be given in the namelist ENVIRON in met.dat. The soil moisture deficit is given by SMD = (SWMAX - SW) / (SWMAX - SWMIN). Thus SW should be in the same units as SWMAX and SWMIN.
The effect of soil moisture on stomatal condutance is calculated as f(SMD) = 1 - SMD1*exp(SMD2*SMD) where SMD1, SMD2 are parameters. This function is implemented for all three stomatal conductance models. The parameters SMD1, SMD2 should be added to the relevant namelists in phy.dat.
In the code you will note that work is also under way to add back in the water balance routines used by Barton & Massheder - however this has not yet been completed.

Quantum yield of electron transport. Previously this was specified as AJQ in namelist JMAXPARS in phy.dat. Alternatively, you can now specify different values according to date, layer, and age class, as you can for JMAX, VCMAX, etc. This is done by adding namelists
&AJQCON    NODATES, NOLAYERS, NOAGES \
&AJQ            DATES, VALUES \
to phy.dat.

Alternative formulations of the temperature dependences of Km and Gamma*.  Use the temperature dependences found by Bernacchi et al (2001) Plant Cell Environ. 24:253-260 by adding the parameter IECO = 0 in namelist JMAXPARS in file phy.dat. The default IECO = 1 is to use the temperature dependences opted for by the ECOCRAFT group.

Estimating incident radiation. The program can use the formula of Bristow & Campbell (1984) Ag For Met 31: 159-166 to calculate incident PAR from air temperature. To use this option it is necessary to put the namelist
& BRISTO    DELTAT    \
in the met.dat file. DELTAT should have 12 values representing the mean monthly daily temperature amplitude for each month of the year.

Stem maintenance respiration on an area basis. To specify stem maintenance respiration on an area basis rather than a mass basis, change the namelist WRESP in phy.dat to have the following parameters:
RMA - maintenance respiration per unit stem surface area (instead of RM, the rate per unit mass)
STEMFORM - used to calculate total stem surface area. Stem SA = STEMFORM * PI * HT * DBH^2
Parameters RTEMP, Q10W, EFFY stay the same.

Outputs. Mean canopy temperature (TCAN) is output on an hourly basis. There is also a column for water evaporation but this is non-operational at the moment.

Bugs. In the version of 14/3/99, a cosine was added to the calculation of WNS1 in TRANSD. This change should not have been made and has been corrected in the current version.

Thermal radiation. The formula used to estimate incident thermal radiation was changed to that of Monteith & Unsworth (1990) Principles of Environmental Physics V2 p53 - this formula takes into account higher radiation from cloudy skies.

Understorey. A program has been written to calculate radiation absorption and photosynthesis by a grass understorey. This program is not available on this web site but may be obtained by writing to Belinda Medlyn.

CHANGES TO 12/1/00
Fixed millennium bug. Dates should still be specified as ‘DD/MM/YY’. YY = 50 - 99, years 1950-1999 assumed; YY = 00 - 49, years 2000 – 2049 assumed. Someone else can change the program come 2050 …

Powerstation Version: Changes made so that program compiles in Microsoft Fortran PowerStation. IMPORTANT:  in input files, arrays involving strings (e.g. dates) must come at the end of the namelist, after all numbers – this means a change to the input file.

Added incident total short-wave radiation to hourly input met data options. Column ‘RAD’ in W m-2.

Added new VPD function for Jarvis stomatal conductance.
    F(VPD) = min(1,1/(VK1*VPD^VK2))
with parameters VK1 and VK2 in namelist JARGS in phy.dat. To use, set j = 3 in namelist MODELGS in confile.dat.

Bugs Fixed: (in order of importance)
(1) The photosynthesis subroutine PHOTOSYN was returning gross photosynthesis, not net, for the Ball-Berry-Leuning case.
(2) A correction to the subroutine DIST in radn.for to handle the case where the ray is completely outwith the crown.
(3) Handle correctly the cases where Jmax or Vcmax <= 0, where photosynthesis is below the light compensation point, and to calculate Gamma* when T < 0.
(4) If physiology is only specified for one age class, but LAD distributions are specified for several, the physiology data must be copied into arrays for the other age classes.
(5) In the calculation of the beam fraction of incident radiation, a more accurate estimate of the zenith angle is now used.

New Parameters Added:
(1) The initial slope of the light-response curve of electron transport is now the parameter AJQ (namelist: JMAXPARS in phy.dat). Default value: ALPHAQ = 0.425 mol mol-1 (defined in maestcom).
(2) Added parameters to force Jmax and Vcmax to go to zero at low temperatures. They begin to decline linearly to zero at T = TVJUP °C and reach zero at T = TVJDN °C. TVJUP, TVJDN are in namelist VCMAXPARS in phy.dat. Default values: -100.

BUG in TRANSD fixed. There was a missing cosine in the calculation of WNS1 in TRANSD, which affected the calculation of absorbed diffuse radiation.

The MAESTEST program. I added a new program to the family, which prints out radiation incident at specified grid points. It is useful for testing measurements made with PAR sensors. The input file, points.dat, contains the co-ordinates of the sensors. The output file, testflx.dat, gives beam, diffuse, and total transmission of PAR to the sensors for each half-hour.

Direct beam fraction at low sun angles. The program uses Spitter's routine to differentiate beam & diffuse fractions. This was only developed for sun angles > 80 deg (from vertical). The problem was how to estimate the diffuse fraction for very low sun angles. I have now assumed that all radiation is diffuse > 80 deg.

Respiration is now set to zero if the air temperature is < TBELOW degrees.TBELOW is specified in the RDPARS namelist in phy.dat.

Daily shortwave irradiance. In daily met files, you can specifiy incident radiation as global shortwave radiation (in MJ m-2 d-1) - the column heading should be 'SI'. Either 'SI' or 'PAR' must be specified.

The calculation of photosynthesis on sunlit foliage can now be done separately for each angle class, by setting the switch MODELSS = 2 (in MODEL namelist, confile.dat).

Incident radiation. There were a couple of problems with the calculation of PAR in W m-2 and the estimation of incident NIR and longwave radiation. The following changes were made:
  1. the conversion factor from umol to J (UMOLPERJ) in maestcom was changed from 2.1 to 4.57 (see Jones 1992 Table 2.2).
  2. the incident NIR is estimated from PAR using the fraction FPAR (in maestcom, currently 0.5), which is the fraction of total incident shortwave which is PAR.
  3. the equation used to calculate incident thermal radiation was changed to: (see Leuning et al 1995 PC&E 18:1183-1200) thermal = 0.642*(EA/TK(TAIR(I)))**(1./7.)*SIGMA*(TK(TAIR(I))**4) where EA is the atmospheric water vapour pressure. These changes will only affect calculation of transpiration when leaf boundary layers conductances are very low.
Effects of OTC. If met data has been measured outside the chamber, the effects of the chamber on climate can be incorporated using the following namelist (in confile.dat): If temperature is affected, the absolute humidity is assumed constant and VPD and RH adjusted accordingly.

Leaf phenology. The model of leaf expansion used by Wang, Rey & Jarvis (1998, GCB) has been added. This can be used instead of specifying leaf areas. Note, however, that it does assume that all trees have the same leaf area. Instead of the namelist INDIVLAREA or ALLLAREA, add the following to trees.dat:

Respiration outputs. Daily respiration fluxes are all output to a separate file, resp.dat. To get this output you must add IORESP = 1 to the CONTROL namelist in confile.dat. There is a also a separate program (resp.exe) which only calculates the respiration fluxes, allowing you to do these calculations independently of radiation interception & photosynthesis.

Respiration. The program can now calculate growth and maintenance respiration for foliage, stem, branch, fine root and coarse root separately. Foliage maintenance respiration, and stem growth and maintenance are calculated as before. The following namelists give the parameters required to calculate branch and root maintenance respiration rates (in phy.dat):

Note that root respiration rates depend on the soil temperature. Soil T is now calculated as the mean daily air temperature. The biomass of roots and branches is estimated from stem height and diameter according to the allometic relation: Fine roots are assumed to be a constant fraction of the total root biomass. The following namelists are used to specify the parameters of this relation (in str.dat): Note that an intercept has also been added to this relation for stem biomass; the parameter WINTERC can be added to the namelist ALLOM in str.dat. The biomass increment of each component is also calculated and used to estimate growth respiration according to where EFFY is in g g-1 C. Only two efficiencies are used: one for foliage and fine roots (EFFYRF) and one for stem, branch and coarse roots (EFFY). EFFYRF is added to the RDPARS namelist and EFFY is added to the WRESP namelist, both in phy.dat. Foliage biomass is calculated from leaf area and SLA.

SLA can be specified by layer, age and date, in a similar fashion to JMAX, VCMAX & LEAFN. Use the following namelists (in phy.dat):

This value of SLA is only used for the calculation of foliage growth respiration, nothing else.

PAR Histogram. The program can print out a frequency histogram of PAR within each tree crown during the simulation. To get the histogram, put IOHIST = 1 in the CONTROL namelist in confile.dat. Then add a namelist which gives the size of the histogram intervals (in umol m-2 s-1) to confile.dat. An example:

Climate change scenarios. The program can add a fixed amount to CO2 concentration and temperature from the met file. The following namelist should be put in confile.dat: The temperature increase is assumed to affect both air temperature and soil temperature. Absolute air humidity (g m-3) is assumed to remain constant; RH & VPD are adjusted accordingly.

Stomatal conductance (Jarvis model). I added some alternative functions for responses in the Jarvis model.

To indicate which model is wished: MODELGS (in confile.dat) can be up to three digits long, e.g. ijk. If: Pressure, vapour pressure mole fraction deficit, dewpoint temp. All of these can now be entered as hourly values in met data files. The column headings are 'PRESS' (Pa), 'MFD' (mmol mol-1), 'TDEW' (deg C). If RH and VPD are both missing, MFD or TDEW will be used to calculate them.

Sumtrees program. This program takes normal output files, in which fluxes are specified for each tree (dayflx.dat, hrflux.dat,histo.dat) and computes totals on an area basis. You can specify which trees you want summed, and weightings for them (e.g. for diameter classes) by using an input file sumcon.dat with the following namelist:

where ITARGETS are the numbers of the trees you want summed, and weights are the weights needed. If this file is not available, all trees are summed and given equal weights.   Bug fix: Foliage respiration was being subtracted twice from the gross photosynthetic rate. Correction to subroutine PHOTOSYN in physiol.for.   In the met data file, hourly relative humidity can be entered as % if wished. The column should be called 'RH%'.

  The crown can have a box shape. Put CSHAPE = 'BOX' in the CANOPY namelist in str.dat.

  You can specify a list of target trees to use, rather than just one. In confile.dat, add ITARGETS = list of target tree numbers to the TREESCON namelist (instead of NOTARGET or NORANDOM). e.g.

&TREESCON
NOTREES = 50
ITARGETS = 5 10 15 20

/

  The wind speed declines exponentially with depth in the canopy. The exponential coefficient is given in STR.DAT:

&AERO
EXTWIND

/

where the windspeed at canopy depth D, W(D) is given by W(D) = W(0)*exp(-EXTWIND*D). Default: EXTWIND = 0.

  Transpiration rate is also calculated from the conductance of the whole canopy. This is printed out in hrflux.dat to compare with the transpiration rate calculated by applying the Penman-Monteith equation to each grid point. The canopy stomatal conductance is calculated as the sum of the individual leaf conductances. The canopy boundary layer conductance is calculated from

U is the wind speed, K is von Karman's constant (added to MAESTCOM), and new parameters required are (in trees.dat):
&AERODYN
ZHT (measurement height, m)
ZPD (zero-plane displacement, m)
Z0HT (roughness length, m)

/

If these parameters aren't specified, this calculation is not done. This part is not fully tested yet.

  Entered the equations for stem respiration rate from Collelongo (Valentini et al 1996, Global Change Biol 2: 199 - 207). The stem respiration rate is given by

The parameters required (in phy.dat) are
&COLLWRESP
CA (constant: paper gives 0.217)
CK (constant: paper gives 0.065)
STEMSDW (the stem surface area: dry weight ratio in m2 kg-1)

/

To indicate that you want to use this model, add MODELRW = 1 to the MODEL namelist in confile.dat. The parameters of the temperature dependence, ie RTEMPW and Q10W should be entered using the WRESP namelist as usual.

  You can now enter the mean leaf inclination angle and the program will calculate the fractions of foliage area in each angle class assuming an elliptical distribution (instead of having to enter the ELP parameter). Just add the parameter AVGANG (in degrees) to the LIA namelist in str.dat. (NB previous methods of entering the LIA distribution are still supported).

  Leaf temperature is now printed out in PTFLX.DAT.

The program will run for more than one tree at a time. Using controls in confile.dat, you can make the program do the calculations for a specified number of randomly-chosen trees, or for all trees in the plot. If you do not want to include trees around the edge of the plot, you can specify the edge width, and trees within this distance from the edge will not be used as target trees.

A minimum (cuticular) conductance was added to the Jarvis model, to avoid the stomatal conductance going to zero.

The units of the output files were changed to more sensible units (mainly umol m-2 s-1). The new units are shown in the output files.

The photosynthesis calculations can be done for sunlit and shaded foliage separately, or done assuming that the absorbed radiation at a grid point is spread evenly over the foliage at that grid point. Use the flag MODELSS in confile.dat to control this.

The latitude and longitude are now specified in the met.dat file.

You can read the hourly or daily CO2 concentrations from the met file. Add a column with the title 'CO2' to MET.DAT. Values should be in umol mol-1.

Added Spitter's formula to calculate the diffuse fraction of radiation. For hourly data, FBEAM is no longer necessary. For daily data, the incident radiation should be entered in columns 'PAR' and 'FBEAM' (rather than RBM and RDF), where FBEAM is not necessary. (Spitters et al 1986, Agric For Meteorol 38: 217-229)

Bug fix - Fixed a bug in the subroutine SUN - the date used to calculate the daylength was fixed! Added a routine to calculate the julian date (DAYJUL). Also added a routine to calculate the Julian day (JDATE).

Night-time foliar respiration is calculated and subtracted from hourly and daily CO2 flux.

Woody maintenance respiration is calculated from the woody biomass, which in turn is calculated using an allometric relation with height and diameter. The diameter should be input in the file trees.dat, in a similar format to the height and radius, i.e.

NAMELIST /ALLDIAM/ NODATES, DATES, VALUES

or

NAMELIST /INDIVDIAM/ NODATES, DATES, VALUES

where the diameter is in m. The allometric relationship is given by:

WBIOM = COEFFT * HEIGHT * DIAM^EXPONT (1)

where: WBIOM is the woody biomass (in kg DW tree-1), HEIGHT is the total height of the tree in m (green height + trunk length), DIAM is the diameter in m, and COEFFT and EXPONT are parameters. The diameters should be specified at the height for which they were used to develop this allometric relationship. The allometric parameters should be added to the file STR.DAT using the namelist:

NAMELIST /ALLOM/ COEFFT, EXPONT

The woody maintenance respiration rate is given by a Q10 relationship:

RESPW = RMW * exp(Q10W * (TAIR - RTEMPW)) * WBIOM (2)

where RMW is the maintenance respiration rate (in umol kg-1 DW s-1 or nmol g-1 DW s-1) at temperature RTEMPW and Q10W is the temperature response factor. The maintenance respiration is calculated on an hourly basis and output in the hourly and daily files. The daily woody growth respiration is calculated from the increment in woody biomass on that day according to

RESPWG = WBINC * EFFY

where EFFY is the efficiency of woody biomass growth (mol CO2 g-1 DW) and WBINC is the daily increment in woody biomass (g DW d-1) which is calculated from the changes in height and diameter, using the allometric relation (1). The parameters required for the calculation of woody respiration rates should be added to PHY.DAT using

NAMELIST /WRESP/ EFFY,RMW,RTEMPW,Q10W

Note that the photosynthesis values reported in the hourly and daily output files are net of foliar respiration but not of woody respiration. The woody maintenance respiration is reported in both the hourly and daily output files. The woody growth respiration is calculated on a daily basis & so is reported only in the daily output file. Note that Franz says woody growth respiration will probably lag the actual growth by approx. one week. We could add a lag to the calculation of the growth rate. Also we could add some kind of fix to make a relationship with temperature. However we don't have the time course of the growth in any fine detail so it's probably not worth making these corrections to the growth respiration calculation.

The temperature at which respiration rates (foliage, woody) are specified must be given - RTEMP (deg C) in the physiology file.

The program now integrates the physiology parameters to give the correct value for each layer (you used to have to have the radiation layers a multiple of the physiology layers). This is for the reflectance, transmittance, Jmax, Vcmax , leaf N and Rd. In the input files:

Output in layers and for individual points. If you set IOHRLY (in confile.dat) to 2, the values of absorbed PAR, photosynthesis and transpiration for each layer for each hour will be printed to the file layflx.dat. There is also another main program which only runs one hour and prints out values for individual points. When compiling, replace maestra.for by maeshr.for. When running maeshr, it will ask which hour you want to do the simulation for. It will then output results for individual points (by layer) to the file ptflx.dat.

Added new shapes for the canopy. The new possible shapes are:

Made some changes to the way BETA is calculated. Abandoned the beta interpolation table - don’t think it helps much anyway? Also, note that the integral of BETA in the vertical direction must be 1, and in the horizontal direction 2 PI. - this determines what a and d must be in the beta formula.

Re-structured the input files to make them more self-contained. I did not want to change the format of the input files but the division of parameters between files was not very logical. I think it is now very logical and should not have to be changed again. Now:

Worked through the scattered radiation routines, adding comments, re-structuring code - and fixed one bug (subtracted sky diffuse radiation from downward diffuse flux). The scattered radiation routines still need more work though! Added parameter DIFSKY to namelist ENVIRON in met.dat. This parameter controls the distribution of diffuse radiation incident from the sky. The distribution function is given by

G(zenith) = (1 + DIFSKY*cos(zenith)) / (1 + 2/PI * DIFSKY)

For a uniform overcast sky, DIFSKY = 0.0 (default). See Wang & Jarvis (1990) Silva Carelica 15:167-180 and Steven & Unsworth (1980) Quart. J. Roy. Meteorol. Soc. 105: 593-602 for more information about DIFSKY.

Included a simulation of shadecloth covering a greenhouse around the trees. I assume that the greenhouse is rectangular in shape and the corners are at co-ordinates (0,0), (0,YMAX), (XMAX,0), (XMAX,YMAX). I added a parameter SHADEHT to namelist PLOT in trees.dat, which is the height of shadecloth on the greenhouse (in m). It is assumed that no rays pass through the shadecloth.