#include #include #include #include #include #include #define PI 3.14159 // pi #define PI180 0.017453292 // pi divided by 180, radians per degree #define PI9 2.864788976 #define PI2 6.283185307 // 2 time pi /* ---------------------------------------------------- 7-23-2002, Leaf_energy_ps.C DENNIS BALDOCCHI Ecosystem Science Division Department of Environmental Science, Policy and Management & Berkeley Atmospheric Science Center 151 Hilgard Hall University of California, Berkeley Berkeley, CA 94720-3110 baldocchi@nature.berkeley.edu 510-642-2874 ------------------------------------------------- The mechanistic and biochemical photosynthesis model of Farquhar is used to scale leaf CO2 exchange to the canopy. Stomatal conductance is computed using the Ball-Berry routine, that is dependent on leaf photosynthesis, relative humidity and the CO2 concentration at the leaf's surface. Photosynthesis and stomatal conductance are computed using a cubic analytical solution I developed from the coupled models during my stay in Viterbo, 92/93 (Baldocchi, 1994). Photosynthetic parameters (Vcmax and Jmax) are scaled with height according to specific leaf weight. SLW is a surrogate for the effect of leaf nitrogen on these parameters. Kinetic photosynthetic coefficients are derived from Peter Harley's 1992 field measurements at WBW. For the seasonal runs Vcmax and Jmax are scaled to seasonal changes in leaf area index. Analytical (quadratic) solutions are used for the leaf and soil surface energy balances Information on the photosynthetic and stomatal conductance model are reported in: Baldocchi, D.D. 1994. An analytical solution for coupled leaf photosynthesis and stomatal conductance models. Tree Physiology 14: 1069-1079. Information on the leaf photosynthetic parameters can be found in: Harley, P.C. and Baldocchi, 1995.Scaling carbon dioxide and water vapor exchange from leaf to canopy in a deciduous forest:leaf level parameterization. Plant, Cell and Environment. 18: 1146-1156. Wilson, K.B., D.D. Baldocchi and P.J. Hanson. 2000. Spatial and seasonal variability of photosynthesis parameters and their relationship to leaf nitrogen in a deciduous forest. Tree Physiology. 20, 565-587. Tests of the model are reported in: Baldocchi, D.D. 1997. Measuring and modeling carbon dioxide and water vapor exchange over a temperate broad-leaved forest during the 1995 summer drought. Plant, Cell and Environment. 20: 1108-1122 Baldocchi, D.D. and P.C. Harley. 1995. Scaling carbon dioxide and water vapor exchange from leaf to canopy in a deciduous forest: model testing and application. Plant, Cell and Environment. 18: 1157-1173. Baldocchi, D.D and T.P. Meyers. 1998. On using eco-physiological, micrometeorological and biogeochemical theory to evaluate carbon dioxide, water vapor and gaseous deposition fluxes over vegetation. Agricultural and Forest Meteorology 90: 1-26. Baldocchi, D.D. Fuentes, J.D., Bowling, D.R, Turnipseed, A.A. Monson, R.K. 1999. Scaling isoprene fluxes from leaves to canopies: test cases over a boreal aspen and a mixed species temperate forest. J. Applied Meteorology. 38, 885-898. Baldocchi, D.D. and K.B.Wilson. 2001. Modeling CO2 and water vapor exchange of a temperate broadleaved forest across hourly to decadal time scales. Ecological Modeling 142: 155-184 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ DEBUGGING NOTES // 11-25-2002 fixed root3 Note done debugging. */ //***************************************************************** // Declare Subroutines // computes direct and diffuse radiation from radiation inputs void DIFFUSE_DIRECT_RADIATION(); // computes the incoming longwave radiation flux density double SKY_IR(double a); // computes longwave radiation flux density void IRFLUX(double a); // computes the leaf orientation function for diffuse radiation void G_FUNC_DIFFUSE(); // computes the leaf orientation angle for the given sun angle and // a spherical leaf inclination void GFUNC(); // sun angles void ANGLE(); // updates leaf reflectance and transmissivity void LEAF_OPTICS(); // photosynthesis, stomatal conductance and respiration // computes leaf photosynthesis with Farquhar model void PHOTOSYNTHESIS(double a,double *b,double c, double d,double *f, double *g, double *h, double *i, double *j, double *k); // Arrenhius function for temperature double TEMP_FUNC(double a,double b,double c,double d, double e); // Boltzmann function for temperature double TBOLTZ(double a,double b, double c, double d); // leaf boundary layer resistances for heat, water, CO2 void BOUNDARY_RESISTANCE(); // leaf energy balance and photosynthesis subroutines and functions void ENERGY_AND_CARBON_FLUXES(); // computes leaf energy balance void ENERGY_BALANCE(double b, double *apc,double d,double e, double f, double g, double *hpt, double *ipt, double *jpt); // humidity at leaf surface double SFC_VPD(double a, double *c); // latent heat of vaporization double LAMBDA(double a); // saturation vapor pressure function double ES(double a); // first derivative of saturation vapor pressure function double DESDT(double a); // second derivative of saturation vapor pressure function double DES2DT(double a); // declare subroutine for inputting data void INPUT_DATA(); // house keeping routines void OUTPUT_DATA(); // defines file names and opens files void FILENAMES(); // Declare Structures // meteorolotical inputs // &dayy,&hhrr,&ta,&rglobal,&parin,&pardif,&ea,&wnd,&ppt,&co2air,&press_mb,&flag struct input_variables { int dayy; // day int hhrr; // hour float ta; // air temperature, C float rglobal; // global radiation, W m-2 float parin; // photosynthetically active radiation, micromole m-2 s-1 float ea; // vapor pressure, kPa float wnd; // wind speed, m s-1 float ppt; // precipitation, mm per hour float co2air; // CO2 concentration, ppm float press_mb; // air pressure, mb long int flag; // input coding } input; // structure for time variables struct time_variables { double local_time; long int daytime; // day+hour coding, e.g. 1230830 int year; // year int days; // day int jdold; // previous day int fileyear; // year for filenames int filestart; // time for filestart int leafout; // day of leaf out int fulleaf; // date of full leaf int count; // number of iterations double lai; // lai as a function of time } time_var; // structure for meteorological variables struct meteorology { double ustar; // friction velocity, m s-1 double ustarnew; // updated friction velocity with new H, m s-1 double rhova_g; // absolute humidity, g m-3 double rhova_kg; // absolute humidity, kg m-3 double sensible_heat_flux; // sensible heat flux, W M-2 double H_old; // old sensible heat flux, W m-2 double air_density; // air density, kg m-3 double T_Kelvin; // absolute air temperature, K double press_kpa; // station pressure, kPa double press_bars; // station pressure, bars double press_Pa; // pressure, Pa double pstat273; // gas constant computations double air_density_mole; // air density, mole m-3 double relative_humidity; // relative humidity, ea/es(T) double vpd; // vapor pressure deficit double ir_in; // infrared flux density } met; // structure for surface resistances and conductances struct surface_resistances { double gcut; // cuticle conductances, mol m-2 s-1 double kballstr; // drought factor to Ball Berry Coefficient double rcuticle; // cuticle resistance, m2 s mol-1 }sfc_res; // structure for plant and physical factors struct factors { double latent; // latent heat of vaporization, J kg-1 double latent18; // latent heat of vaporization times molecular mass of vapor, 18 g mol-1 double heatcoef; // factor for sensible heat flux density double a_filt; // filter coefficients double b_filt; // filter coefficients double co2; // CO2 factor, ma/mc * rhoa (mole m-3) } fact; // structure for leaf fluxes struct leaf_fluxes{ double latent_heat_flux; double sensible_heat_flux; double IR_out; double photosynthesis_mole; double stomatal_res; double stomatal_cond; double respiration; double Wc; double Wj; double internal_CO2; double T_sfc_K; double T_sfc_C; }leaf; // structure for non dimensional variables struct non_dimensional_variables { // Prandtl Number double pr; double pr33; // Schmidt number for vapor double sc; double sc33; // Schmidt number for CO2 double scc; double scc33; // Schmidt number for ozone double sco3; double sco333; // Grasshof number double grasshof; // multiplication factors with leaf length and diffusivity double lfddv; double lfddh; } non_dim; // boundary layer resistances struct boundary_layer_resistances{ double vapor; // resistance for water vapor, s/m double heat; // resistance for heat, s/m double co2; // resistance for CO2, s/m } bound_layer_res; // radiation variables, visible, near infrared and infrared struct solar_radiation_variables { // inputs of near infrared radiation and components, W m-2 double nir_beam; // incoming beam component near infrared radiation, W m-2 double nir_diffuse; // incoming diffuse component near infrared radiation, W m-2 double nir_total; // incoming total near infrared radiaion, W m-2 // computed profiles of near infrared radiation, W m-2 double nir_dn; // downward scattered near infrared radiation double nir_up; // upward scattered near infrared radiation double nir_sun; // near infrared radiation on sunlit fraction of layer double nir_sh; // near infrared radiation on shaded fraction of layer double beam_flux_nir; // flux density of direct near infrared radiation // leaf and soil optical properities of near infrared radiation double nir_reflect; // leaf reflectance in the near infrared double nir_trans; // leaf transmittance in the near infrared double nir_soil_refl; // soil reflectance in the near infrared double nir_absorbed; // leaf absorptance in the near infrared // inputs of visible light, PAR, W m-2 double par_diffuse; // diffuse component of incoming PAR, parin double par_beam; // beam component of incoming PAR, parin // computed profiles of visible radiation, PAR, W m-2 double par_shade; // PAR on shaded fraction of layer double par_sun; // PAR on sunlit fraction of layer, beam and diffuse double beam_flux_par; // PAR in the beam component double par_down; // upward scattered PAR double par_up; // downward scattered PAR // flux densities of visible quanta on sun and shade leaves for photosynthesis // calculations, micromoles m-2 s-1 double quantum; // par on sunlit leaves // optical properties of leaves and soil for PAR double par_absorbed; // PAR leaf absorptance double par_reflect; // PAR leaf reflectance double par_trans; // PAR leaf transmittance double par_soil_refl; // PAR soil reflectance // Net radiation profiles, W m-2 double rnet; double exxpdir; // exponential transmittance of diffuse radiation through a layer double beta_rad; // solar elevation angle, radians double sine_beta; // sine of beta double beta_deg; // solar elevation angle, degrees double ratrad; // radiation ratio to detect cloud amount double ratradnoon; // radiation ratio at noon for guestimating cloud amount at night double diff_tot_par; // ratio par diffuse to total } solar; // ------------------------------------------------------ // DECLARE PARAMETER VALUES // ------------------------------------------------------ const double lai = 6.0; // original lai=6.0 Leaf area index const double vcopt = 73.0; // carboxylation rate at optimal temperature, umol m-2 s-1 const double jmopt = 170.0; // electron transport rate at optimal temperature, umol m-2 s-1 const double rd25 = .34; // dark respiration at 25 C, rd25= 0.34 umol m-2 s-1 const double pi4=12.5663706; // Universal gas constant const double rugc = 8.314; // J mole-1 K-1 const double rgc1000 = 8314; // gas constant times 1000. // Consts for Photosynthesis model and kinetic equations. // for Vcmax and Jmax. Taken from Harley and Baldocchi (1995, PCE) const double hkin = 200000.0; // enthalpy term, J mol-1 const double skin = 710.0; // entropy term, J K-1 mol-1 const double ejm = 55000.0; // activation energy for electron transport, J mol-1 const double evc = 55000.0; // activation energy for carboxylation, J mol-1 // Enzyme constants & partial pressure of O2 and CO2 // Michaelis-Menten K values. From survey of literature. const double kc25 = 274.6; // kinetic coef for CO2 at 25 C, microbars const double ko25 = 419.8; // kinetic coef for O2 at 25C, millibars const double o2 = 210.0; // oxygen concentration mmol mol-1 // tau is computed on the basis of the Specificity factor (102.33) // times Kco2/Kh2o (28.38) to convert for value in solution // to that based in air/ // The old value was 2321.1. // New value for Quercus robor from Balaguer et al. 1996 // Similar number from Dreyer et al. 2001, Tree Physiol, tau= 2710 const double tau25 = 2904.12; // tau coefficient // Arrhenius constants // Eact for Michaelis-Menten const. for KC, KO and dark respiration // These values are from Harley const double ekc = 80500.0; // Activation energy for K of CO2; J mol-1 const double eko = 14500.0; // Activation energy for K of O2, J mol-1 const double erd = 38000.0; // activation energy for dark respiration, eg Q10=2 const double ektau = -29000.0; // J mol-1 (Jordan and Ogren, 1984) const double tk_25 = 298.16; // absolute temperature at 25 C const double toptvc = 311.0; // optimum temperature for maximum carboxylation const double toptjm = 311.0; // optimum temperature for maximum electron transport const double eabole=45162; // activation energy for bole respiration for Q10 = 2.02 // Constants for leaf energy balance const double sigma = 5.67e-08; // Stefan-Boltzmann constant W M-2 K-4 const double cp = 1005.; // Specific heat of air, J KG-1 K-1 const double mass_air = 29.; // Molecular weight of air, g mole-1 const double mass_CO2=44.; // molecular weight of CO2, g mole-1 const double dldt = -2370.; // Derivative of the latent heat of vaporization const double ep = .98; // emissivity of leaves const double epm1=0.02; // 1- ep const double epsoil = .98; // Emissivity of soil const double epsigma=5.5566e-8; // ep*sigma const double epsigma2 = 11.1132e-8; // 2*ep*sigma const double epsigma4 = 22.2264e-8; // 4.0 * ep * sigma const double epsigma6 = 33.3396e-8; // 6.0 * ep * sigma const double epsigma8 = 44.448e-8; // 8.0 * ep * sigma const double epsigma12= 66.6792e-8; // 12.0 * ep * sigma const double betfact=1.5; // multiplication factor for aerodynamic // sheltering, based on work by Grace and Wilson // constants for the polynomial equation for saturation vapor pressure-T function, es=f(t) const double a1en=617.4; const double a2en=42.22; const double a3en=1.675; const double a4en=0.01408; const double a5en=0.0005818; // Ball-Berry stomatal coefficient for stomatal conductance const double kball = 9.5; // intercept of Ball-Berry model, mol m-2 s-1 const double bprime = .0175; // intercept for H2O const double bprime16 = 0.0109375; // intercept for CO2, bprime16 = bprime / 1.6; // Minimum stomatal resistance, s m-1. const double rsm = 145.0; const double brs=60.0; // curvature coeffient for light response // leaf quantum yield, electrons const double qalpha = .22; const double qalpha2 = 0.0484; // qalpha squared, qalpha2 = pow(qalpha, 2.0); // Leaf dimension. geometric mean of length and width (m) const double lleaf = .1; // leaf length, m // Diffusivity values for 273 K and 1013 mb (STP) using values from Massman (1998) Atmos Environment // These values are for diffusion in air. When used these values must be adjusted for // temperature and pressure // nu, Molecular viscosity const double nuvisc = 13.27; // mm2 s-1 const double nnu = 0.00001327; // m2 s-1 // Diffusivity of CO2 const double dc = 13.81; // mm2 s-1 const double ddc = 0.00001381; // m2 s-1 // Diffusivity of heat const double dh = 18.69; // mm2 s-1 const double ddh = 0.00001869; // m2 s-1 // Diffusivity of water vapor const double dv = 21.78; // mm2 s-1 const double ddv = 0.00002178; // m2 s-1 // Declare file pointers for I/O FILE *fptr1; main () { // Define Filenames FILENAMES(); // Constants for leaf boundary layers non_dim.lfddh=lleaf/ddh; // PRANDTL NUMBER // Prandtl Number non_dim.pr = nuvisc / dh; non_dim.pr33 = pow(non_dim.pr,.33); // DIFFUSIVITY OF WATER VAPOR, m2 s-1 non_dim.lfddv=lleaf/ddv; // SCHMIDT NUMBER FOR VAPOR non_dim.sc = nuvisc / dv; non_dim.sc33 = pow(non_dim.sc,.33); // SCHMIDT NUMBER FOR CO2 non_dim.scc = nuvisc / dc; non_dim.scc33 = pow(non_dim.scc,.33); // Grasshof Number non_dim.grasshof=9.8*pow(lleaf,3)/pow(nnu,2); /************************************************************************* MAIN PROGRAM ************************************************************************** Describe canopy attributes Flow chart of the Main Program: 1) input met values 2) Compute solar elevation angle; 3) Compute PAR and NIR 4) Compute first estimate of stomatal conductance 5) Compute first estimate of IRFLUX assuming leaf temperature equals air temperature. 6) Compute first estimate of leaf energy balance and leaf temperature 7) Compute photosynthesis, transpiration and stomatal conductance 8) Iterate until convergence 9) Thats all Folks!!!! ************************************************************************/ // set Ball-Berry stomatal factor sfc_res.kballstr = kball; LEAF_OPTICS(); // define leaf optical properties INPUT_DATA(); ANGLE(); // make sure PAR is zero at night. Some data have negative offset, which causes numerical // problems if(solar.sine_beta <= 0.05) input.parin=0; // Compute the fractions of beam and diffuse radiation from incoming measurements // Set the radiation factor to the day before for night calculations. This way if the day // was cloudy, so will IR calculations for the night reflect this. if (solar.sine_beta > 0.05) DIFFUSE_DIRECT_RADIATION(); else solar.ratrad = solar.ratradnoon; // computes leaf inclination angle distribution function, the mean direction cosine // between the sun zenith angle and the angle normal to the mean leaf // for CANOAK we use the leaf inclination angle data of Hutchison et al. 1983, J Ecology // for other canopies we assume the leaf angle distribution is spherical // Loop out if night time met.rhova_kg = met.rhova_g / 1000.; /* absolute humidity, kg m-3 */ fact.heatcoef = met.air_density * cp; // any recursive looping would occur here, below TL initialization IRFLUX(met.T_Kelvin); // iteration looping for energy fluxes and scalar fields // iterate until energy balance closure occurs or 75 iterations are // reached // Compute leaf energy balance, leaf temperature, photosynthesis and stomatal conductance. ENERGY_AND_CARBON_FLUXES(); // update long wave radiation fluxes with new leaf and air temperatures OUTPUT_DATA(); fclose(fptr1); } // end of main function // ====================================================================== // Listing of Subroutines void INPUT_DATA() { double est; // input data and check for bad data // note that the data were produced in single precision (float) // so I had to read them as single precision, otherwise I ingested // garbage time_var.year=2002; time_var.days=180; time_var.local_time=12; input.ta=25; input.rglobal=500; input.parin= (input.rglobal/2.0)*4.6; input.wnd=3.; input.ppt=0; input.co2air=370.; input.press_mb=1013; input.ea=1.0; // estimate friction velocity met.ustar=0.055*input.wnd; met.T_Kelvin=input.ta+273.15; // compute absolute air temperature met.rhova_g = input.ea * 2165/met.T_Kelvin; // compute absolute humidity, g m-3 est=ES(met.T_Kelvin); met.relative_humidity=input.ea*10./est; // relative humidity met.vpd=est-input.ea*10.; // vapor pressure deficit, mb met.press_kpa=input.press_mb/10.; // air pressure, kPa met.press_bars=input.press_mb/1000.; // air pressure, bars met.press_Pa=met.press_kpa*1000.; // pressure, Pa // combining gas law constants met.pstat273 = .022624 / (273.16 * met.press_bars); // cuticular conductance adjusted for pressure and T, mol m-2 s-1 sfc_res.gcut = bprime * met.T_Kelvin * met.pstat273; // cuticular resistance sfc_res.rcuticle = 1.0 / sfc_res.gcut; // air density, kg m-3 met.air_density = met.press_kpa * mass_air / (rugc * met.T_Kelvin); // air density, mole m-3 met.air_density_mole=met.press_kpa/ (rugc * met.T_Kelvin) * 1000.; return; } void FILENAMES() { fptr1=fopen("c:\\stomata\\energy_balance.dat","w"); fprintf(fptr1,"%s, %s, %s, %s, %s, %s\n","latent_heat_flux", "sensible_heat_flux", "photosynthesis_mole", "respiration", "stomatal_cond", "T_sfc_K"); return; } void OUTPUT_DATA() { printf("%s, %s, %s, %s, %s, %s\n","latent_heat_flux", "sensible_heat_flux", "photosynthesis_mole", "respiration", "stomatal_cond", "T_sfc_K"); printf("%f, %f, %f, %f, %f, %f\n",leaf.latent_heat_flux, leaf.sensible_heat_flux, leaf.photosynthesis_mole, leaf.respiration, leaf.stomatal_cond, leaf.T_sfc_K); fprintf(fptr1,"%f, %f, %f, %f, %f, %f\n",leaf.latent_heat_flux, leaf.sensible_heat_flux, leaf.photosynthesis_mole, leaf.respiration, leaf.stomatal_cond, leaf.T_sfc_K); return; } void ENERGY_AND_CARBON_FLUXES() { /* The ENERGY_AND_CARBON_FLUXES routine to computes coupled fluxes of energy, water and CO2 exchange, as well as leaf temperature. Computataions are performed for each layer in the canopy and on the sunlit and shaded fractions. Analytical solution for leaf energy balance and leaf temperature is used. The program is derived from work by Paw U (1986) and was corrected for errors with a re-derivation of the equations. The quadratic form of the solution is used, rather than the quartic version that Paw U prefers. Estimates of leaf temperature are used to adjust kinetic coefficients for enzyme kinetics, respiration and photosynthesis, as well as the saturation vapor pressure at the leaf surface. The Analytical solution of the coupled set of equations for photosynthesis and stomatal conductance by Baldocchi (1994, Tree Physiology) is used. This equation is a solution to a cubic form of the photosynthesis equation. The photosynthesis algorithms are from the model of Farquhar. Stomatal conductance is based on the algorithms of Ball- Berry and Collatz et al., which couple gs to A */ double T_sfc_K; // surface temperatures in Kelvin and Centigrade // initial guess of stomatal conductance leaf.stomatal_res= 50 +5000./input.parin; /* First compute energy balance. The stomatal resistances on the sunlit and shaded leaves are pre-estimated as a function of PAR using STOMATA Remember layer is two-sided so include upper and lower streams are considered. KC is the convective transfer coeff. (W M-2 K-1). A factor of 2 is applied since convective heat transfer occurs on both sides of leaves. To calculate the total leaf resistance we must combine stomatal and the leaf boundary layer resistance. Many crops are amphistomatous so KE must be multiplied by 2. Deciduous forest, on the other hand has stomata on one side of the leaf. */ // Initialize surface temperature with air temperature T_sfc_K = met.T_Kelvin; // Energy balance on sunlit leaves // update latent heat with new temperature during each call of this routine fact.latent = LAMBDA(T_sfc_K); fact.latent18=fact.latent*18.; // initial stomatal resistance as a function of PAR. Be careful and use // the light dependent version only for first iteration // Compute the resistances for heat and vapor transfer, rh and rv, // for each layer, s/m BOUNDARY_RESISTANCE(); // compute energy balance of sunlit leaves // define absorbed solar radiation, PAR + NIR + IR, work in reflectivities and emissivities // need IR on top and bottom of the leaf solar.rnet=input.rglobal+met.ir_in+met.ir_in; // + Lout ENERGY_BALANCE(solar.rnet, &leaf.T_sfc_K, met.T_Kelvin, met.rhova_kg, bound_layer_res.heat, leaf.stomatal_res, &leaf.latent_heat_flux, &leaf.sensible_heat_flux, &leaf.IR_out); // compute photosynthesis PHOTOSYNTHESIS(input.parin, &leaf.stomatal_res, input.co2air, leaf.T_sfc_K, &leaf.latent_heat_flux, &leaf.photosynthesis_mole, &leaf.respiration, &leaf.internal_CO2,&leaf.Wj,&leaf.Wc); leaf.stomatal_cond=1./leaf.stomatal_res; return; } double SKY_IR (double T) { // Infrared radiation from sky, W m-2, using algorithm from Norman double y; y = sigma * pow(T,4.) * ((1. - .261 * exp(-.000777 * pow((273.16 - T), 2.))) * solar.ratrad + 1 - solar.ratrad); return y; } void ENERGY_BALANCE (double qrad, double *tsfckpt, double taa, double rhovva, double rvsfc, double stomsfc, double *lept, double *H_leafpt, double *lout_leafpt) { /* ENERGY BALANCE COMPUTATION A revised version of the quadratic solution to the leaf energy balance relationship is used. Paw U, KT. 1987. J. Thermal Biology. 3: 227-233 H is sensible heat flux density on the basis of both sides of a leaf J m-2 s-1 (W m-2). Note KC includes a factor of 2 here for heat flux because it occurs from both sides of a leaf. */ double est, ea, tkta, le2; double tk2, tk3, tk4; double dest, d2est; double lecoef, hcoef, hcoef2, repeat, acoeff, acoef; double bcoef, ccoef, product; double atlf, btlf, ctlf,vpd_leaf,llout; double ke; // tkta=taa+273.16; tkta=taa; est = ES(tkta) * 100; // converts es(T) from mb to Pa // ea = RHOA * TAA * 1000 / 2.165 ea = rhovva * taa * 461.89; // vapor pressure above leaf // Vapor pressure deficit, Pa vpd_leaf = est - ea; if (vpd_leaf < 0.) vpd_leaf = 0; // Slope of the vapor pressure-temperature curve, Pa/C // evaluate as function of Tk dest = DESDT(tkta); // Second derivative of the vapor pressure-temperature curve, Pa/C // Evaluate as function of Tk d2est = DES2DT(tkta); // Compute products of air temperature, K tk2 = tkta * tkta; tk3 = tk2 * tkta; tk4 = tk3 * tkta; // Longwave emission at air temperature, W m-2 llout = epsigma * tk4; /* Coefficient for latent heat flux Oaks evaporate from only one side. They are hypostomatous. Cuticle resistance is included in STOM. */ ke = 1./ (rvsfc + stomsfc); lecoef = met.air_density * .622 * fact.latent * ke / met.press_Pa; // Coefficients for sensible heat flux hcoef = met.air_density*cp/bound_layer_res.heat; hcoef2 = 2 * hcoef; // The quadratic coefficients for the a LE^2 + b LE +c =0 repeat = hcoef + epsigma4 * tk3; acoeff = lecoef * d2est / (2. * repeat); acoef = acoeff / 4.; bcoef = -(repeat) - lecoef * dest / 2. + acoeff * (-qrad / 2. + llout); ccoef = repeat * lecoef * vpd_leaf + lecoef * dest * (qrad / 2. - llout) + acoeff * ((qrad * qrad) / 4. + llout * llout - qrad * llout); // LE1 = (-BCOEF + (BCOEF ^ 2 - 4 * ACOEF * CCOEF) ^ .5) / (2 * ACOEF) product = bcoef * bcoef - 4. * acoef * ccoef; // LE2 = (-BCOEF - (BCOEF * BCOEF - 4 * acoef * CCOEF) ^ .5) / (2. * acoef) le2= (-bcoef - sqrt(product)) / (2. * acoef); *lept=le2; // need to pass pointer out of subroutine // solve for Ts using quadratic solution // coefficients to the quadratic solution atlf = epsigma12 * tk2 + d2est * lecoef / 2.; btlf = epsigma8 * tk3 + hcoef2 + lecoef * dest; ctlf = -qrad + 2 * llout + lecoef * vpd_leaf; // IF (BTLF * BTLF - 4 * ATLF * CTLF) >= 0 THEN product = btlf * btlf - 4 * atlf * ctlf; // T_sfc_K = TAA + (-BTLF + SQR(BTLF * BTLF - 4 * ATLF * CTLF)) / (2 * ATLF) if (product >= 0) *tsfckpt = tkta + (-btlf + sqrt(product)) / (2 * atlf); else *tsfckpt=tkta; if(*tsfckpt < -230 || *tsfckpt > 325) *tsfckpt=tkta; // long wave emission of energy *lout_leafpt =epsigma2*pow(*tsfckpt,4); // H is sensible heat flux *H_leafpt = hcoef2 * (*tsfckpt- tkta); return; } double LAMBDA (double tak) { // Latent heat of Vaporiation, J kg-1 double y; y = 3149000. - 2370. * tak; // add heat of fusion for melting ice if(tak < 273.) y +=333; return y; } double ES(double t) { // saturation vapor pressure function (mb) // T is temperature in Kelvin double y, y1; if(t > 0) { y1 = (54.8781919 - 6790.4985 / t - 5.02808 * log(t)); y = exp(y1); } else printf("bad es calc"); return y; } double DESDT (double t) { // first derivative of es with respect to tk // routine needs to convert es(t) (mb) to Pa double y; y = ES(t)* 100* fact.latent18 / (rgc1000 * t * t); return y; } double DES2DT(double t) { // The second derivative of the saturation vapor pressure // temperature curve, using the polynomial equation of Paw U // a3en=1.675; // a4en=0.01408; // a5en=0.0005818; double y,tcel; tcel=t-273.16; y=2*a3en+6*a4en*tcel+12*a5en*tcel*tcel; return y; } void PHOTOSYNTHESIS(double Iphoton,double *rstompt, double cca,double tlk, double *leleaf, double *A_molept, double *resppt, double *cipnt, double *wjpnt, double *wcpnt) { /* This program solves a cubic equation to calculate leaf photosynthesis. This cubic expression is derived from solving five simultaneous equations for A, PG, cs, CI and GS. Stomatal conductance is computed with the Ball-Berry model. The cubic derivation assumes that b', the intercept of the Ball-Berry stomatal conductance model, is non-zero. Gs = k A rh/cs + b' We also found that the solution for A can be obtained by a quadratic equation when Gs is constant or b' is zero. The derivation is published in: Baldocchi, D.D. 1994. An analytical solution for coupled leaf photosynthesis and stomatal conductance models. Tree Physiology 14: 1069-1079. ----------------------------------------------------------------------- A Biochemical Model of C3 Photosynthesis After Farquhar, von Caemmerer and Berry (1980) Planta. 149: 78-90. The original program was modified to incorporate functions and parameters derived from gas exchange experiments of Harley, who paramertized Vc and J in terms of optimal temperature, rather than some reference temperature, eg 25C. Program calculates leaf photosynthesis from biochemical parameters rd25 - Dark respiration at 25 degrees C (umol m-2 s-1) tlk - leaf temperature, Kelvin jmax - optimal rate of electron transport vcopt - maximum rate of RuBP Carboxylase/oxygenase iphoton - incident photosynthetically active photon flux (mmols m-2 s-1) note: Harley parameterized the model on the basis of incident PAR gs - stomatal conductance (mols m-2 s-1), typically 0.01-0.20 pstat-station pressure, bars aphoto - net photosynthesis (umol m-2 s-1) ps - gross photosynthesis (umol m-2 s-1) aps - net photosynthesis (mg m-2 s-1) aphoto (umol m-2 s-1) -------------------------------------------------- iphoton is radiation incident on leaves The temperature dependency of the kinetic properties of RUBISCO are compensated for using the Arrhenius and Boltzmann equations. From biochemistry, one observes that at moderate temperatures enzyme kinetic rates increase with temperature. At extreme temperatures enzyme denaturization occurs and rates must decrease. Arrhenius Eq. f(T)=f(tk_25) exp(tk -298)eact/(298 R tk)), where eact is the activation energy. Boltzmann distribution F(T)=tboltz) Define terms for calculation of gross photosynthesis, PG PG is a function of the minimum of RuBP saturated rate of carboxylation, Wc, and the RuBP limited rate of carboxylation, Wj. Wj is limiting when light is low and electron transport, which re-generates RuBP, is limiting. Wc is limiting when plenty of RuBP is available compared to the CO2 that is needed for carboxylation. Both equations take the form: PG-photorespiration= (a CI-a d)/(e CI + b) PG-photorespiration=min[Wj,Wc] (1-gamma/Ci) Wc=Vcmax Ci/(Ci + Kc(1+O2/Ko)) Wj=J Ci/(4 Ci + 8 gamma) Ps kinetic coefficients from Harley at WBW. Gamma is the CO2 compensation point Jan 14, 1999 Updated the cubic solutions for photosynthesis. There are times when the restriction that R^2 < Q^3 is violated. I therefore need alternative algorithms to solve for the correct root. =============================================================== */ double tprime25, bc, ttemp, gammac; double jmax, vcmax, cs, ci; double kct, ko, tau; double rd, rdz; double rb_mole,gb_mole,dd,b8_dd; double rh_leaf, k_rh, gb_k_rh,ci_guess; double j_photon,alpha_ps,bbeta,gamma; double denom,Pcube,Qcube,Rcube; double P2, P3, Q, R; double root1, root2; double root3, arg_U, ang_L; double aphoto, j_sucrose, wj; double gs_leaf_mole, gs_co2,gs_m_s; double ps_1,delta_1,Aquad1,Bquad1,Cquad1; double theta_ps, wc, B_ps, a_ps, E_ps, psguess; double sqrprod, product; double rt; double rr,qqq, minroot, maxroot, midroot; rt = rugc * tlk; // product of universal gas constant and abs temperature tprime25 = tlk - tk_25; // temperature difference ttemp = exp((skin * tlk - hkin) / rt) + 1.0; // denominator term // KC and KO are solely a function of the Arrhenius Eq. kct = TEMP_FUNC(kc25, ekc, tprime25, tk_25, tlk); ko = TEMP_FUNC(ko25, eko, tprime25, tk_25,tlk); tau = TEMP_FUNC(tau25, ektau, tprime25, tk_25,tlk); bc = kct * (1.0 + o2 / ko); if(Iphoton < 1) Iphoton = 0; /* gammac is the CO2 compensation point due to photorespiration, umol mol-1 Recalculate gammac with the new temperature dependent KO and KC coefficients gammac = .5 * O2*1000/TAU */ gammac = 500.0 * o2 / tau; /* temperature corrections for Jmax and Vcmax Scale jmopt and VCOPT with a surrogate for leaf nitrogen specific leaf weight (Gutschick and Weigel). normalized leaf wt is 1 at top of canopy and is 0.35 at forest floor. Leaf weight scales linearly with height and so does jmopt and vcmax zoverh=0.65/HT=zh65 */ // spring, increase Ps capacity with leaf expansion as a function of leaf area changes // growing season, full Ps capacity (note newer data by Wilson et al shows more // dynamics // Scale rd with height via vcmax and apply temperature // correction for dark respiration rdz=vcopt * 0.004657; // reduce respiration by 40% in light according to Amthor if(Iphoton > 10) rdz *= 0.4; rd = TEMP_FUNC(rdz, erd, tprime25, tk_25, tlk); // Apply temperature correction to JMAX and vcmax jmax = TBOLTZ(jmopt, ejm, toptjm, tlk); vcmax = TBOLTZ(vcopt, evc, toptvc, tlk); /* Compute the leaf boundary layer resistance gb_mole leaf boundary layer conductance for CO2 exchange, mol m-2 s-1 RB has units of s/m, convert to mol-1 m2 s1 to be consistant with R. rb_mole = RBCO2 * .0224 * 1.01 * tlk / (met.pstat * 273.16) */ rb_mole = bound_layer_res.co2 * tlk * (met.pstat273); gb_mole = 1. / rb_mole; dd = gammac; b8_dd = 8 * dd; /*************************************** APHOTO = PG - rd, net photosynthesis is the difference between gross photosynthesis and dark respiration. Note photorespiration is already factored into PG. **************************************** coefficients for Ball-Berry stomatal conductance model Gs = k A rh/cs + b' rh is relative humidity, which comes from a coupled leaf energy balance model */ rh_leaf = SFC_VPD(tlk, leleaf); k_rh = rh_leaf * sfc_res.kballstr; // combine product of rh and K ball-berry /* Gs from Ball-Berry is for water vapor. It must be divided by the ratio of the molecular diffusivities to be valid for A */ k_rh = k_rh / 1.6; // adjust the coefficient for the diffusion of CO2 rather than H2O gb_k_rh = gb_mole * k_rh; ci_guess = cca * .7; // initial guess of internal CO2 to estimate Wc and Wj // cubic coefficients that are only dependent on CO2 levels alpha_ps = 1.0 + (bprime16 / gb_mole) - k_rh; bbeta = cca * (gb_k_rh - 2.0 * bprime16 - gb_mole); gamma = cca * cca * gb_mole * bprime16; theta_ps = gb_k_rh - bprime16; /* Test for the minimum of Wc and Wj. Both have the form: W = (a ci - ad)/(e ci + b) after the minimum is chosen set a, b, e and d for the cubic solution. estimate of J according to Farquhar and von Cammerer (1981) J photon from Harley */ if (jmax > 0) j_photon = qalpha * Iphoton / sqrt(1. +(qalpha2 * Iphoton * Iphoton / (jmax * jmax))); else j_photon = 0; wj = j_photon * (ci_guess - dd) / (4. * ci_guess + b8_dd); wc = vcmax * (ci_guess - dd) / (ci_guess + bc); // frost and end of leaf photosynthesis and respiration if(time_var.days > 300) { wj = 0; j_photon = 0; wc=0; rd = 0; } if(wj < wc) { // for Harley and Farquhar type model for Wj psguess=wj; B_ps = b8_dd; a_ps = j_photon; E_ps = 4.0; } else { psguess=wc; B_ps = bc; a_ps = vcmax; E_ps = 1.0; } /* if wj or wc are less than rd then A would probably be less than zero. This would yield a negative stomatal conductance. In this case, assume gs equals the cuticular value. This assumptions yields a quadratic rather than cubic solution for A */ if (wj <= rd) goto quad; if (wc <= rd) goto quad; /* cubic solution: A^3 + p A^2 + q A + r = 0 */ denom = E_ps * alpha_ps; Pcube = (E_ps * bbeta + B_ps * theta_ps - a_ps * alpha_ps + E_ps * rd * alpha_ps); Pcube /= denom; Qcube = (E_ps * gamma + (B_ps * gamma / cca) - a_ps * bbeta + a_ps * dd * theta_ps + E_ps * rd * bbeta + rd * B_ps * theta_ps); Qcube /= denom; Rcube = (-a_ps * gamma + a_ps * dd * (gamma / cca) + E_ps * rd * gamma + rd * B_ps * gamma / cca); Rcube /= denom; // Use solution from Numerical Recipes from Press P2 = Pcube * Pcube; P3 = P2 * Pcube; Q = (P2 - 3.0 * Qcube) / 9.0; R = (2.0 * P3 - 9.0 * Pcube * Qcube + 27.0 * Rcube) / 54.0; /* Test = Q ^ 3 - R ^ 2 if test >= O then all roots are real */ rr=R*R; qqq=Q*Q*Q; // real roots arg_U = R / sqrt(qqq); ang_L = acos(arg_U); root1 = -2.0 * sqrt(Q) * cos(ang_L / 3.0) - Pcube / 3.0; root2 = -2.0 * sqrt(Q) * cos((ang_L + PI2) / 3.0) - Pcube / 3.0; root3 = -2.0 * sqrt(Q) * cos((ang_L - PI2) / 3.0) - Pcube / 3.0; // rank roots #1,#2 and #3 according to the minimum, intermediate and maximum // value if(root1 <= root2 && root1 <= root3) {minroot=root1; if (root2 <= root3) {midroot=root2; maxroot=root3;} else {midroot=root3; maxroot=root2;} } if(root2 <= root1 && root2 <= root3) {minroot=root2; if (root1 <= root3) {midroot=root1; maxroot=root3;} else {midroot=root3; maxroot=root1;} } if(root3 <= root1 && root3 <= root2) {minroot=root3; if (root1 < root2) {midroot=root1; maxroot=root2;} else {midroot=root2; maxroot=root1;} } // end of the loop for real roots // find out where roots plop down relative to the x-y axis if (minroot > 0 && midroot > 0 && maxroot > 0) aphoto=minroot; if (minroot < 0 && midroot < 0 && maxroot > 0) aphoto=maxroot; if (minroot < 0 && midroot > 0 && maxroot > 0) aphoto=midroot; /* Here A = x - p / 3, allowing the cubic expression to be expressed as: x^3 + ax + b = 0 */ // aphoto=root3; // back to original assumption /* also test for sucrose limitation of photosynthesis, as suggested by Collatz. Js=Vmax/2 */ j_sucrose = vcmax / 2. - rd; if(j_sucrose < aphoto) aphoto = j_sucrose; cs = cca - aphoto / gb_mole; if(cs > 1000) cs=input.co2air; /* Stomatal conductance for water vapor forest are hypostomatous. Hence we don't divide the total resistance by 2 since transfer is going on only one side of a leaf. */ gs_leaf_mole = (sfc_res.kballstr * rh_leaf * aphoto / cs) + bprime; // convert Gs from vapor to CO2 diffusion coefficient gs_co2 = gs_leaf_mole / 1.6; /* stomatal conductance is mol m-2 s-1 convert back to resistance (s/m) for energy balance routine */ gs_m_s = gs_leaf_mole * tlk * met.pstat273; // need point to pass rstom out of subroutine *rstompt = 1.0 / gs_m_s; // to compute ci, Gs must be in terms for CO2 transfer ci = cs - aphoto / gs_co2; /* if A < 0 then gs should go to cuticular value and recalculate A using quadratic solution */ if(aphoto <= 0.0) goto quad; goto OUTDAT; // if aphoto < 0 set stomatal conductance to cuticle value quad: gs_leaf_mole = bprime; gs_co2 = gs_leaf_mole / 1.6; /* stomatal conductance is mol m-2 s-1 convert back to resistance (s/m) for energy balance routine */ gs_m_s = gs_leaf_mole * tlk * (met.pstat273); // need point to pass rstom out of subroutine as a pointer *rstompt = 1.0 / gs_m_s; /* a quadratic solution of A is derived if gs=ax, but a cubic form occurs if gs =ax + b. Use quadratic case when A is less than zero because gs will be negative, which is nonsense */ ps_1 = cca * gb_mole * gs_co2; delta_1 = gs_co2 + gb_mole; denom = gb_mole * gs_co2; Aquad1 = delta_1 * E_ps; Bquad1 = -ps_1 * E_ps - a_ps * delta_1 + E_ps * rd * delta_1 - B_ps * denom; Cquad1 = a_ps * ps_1 - a_ps * dd * denom - E_ps * rd * ps_1 - rd * B_ps * denom; product=Bquad1 * Bquad1 - 4.0 * Aquad1 * Cquad1; if (product >= 0) sqrprod= sqrt(product); aphoto = (-Bquad1 - sqrprod) / (2.0 * Aquad1); /* Tests suggest that APHOTO2 is the correct photosynthetic root when light is zero because root 2, not root 1 yields the dark respiration value rd. */ cs = cca - aphoto / gb_mole; ci = cs - aphoto / gs_co2; OUTDAT: /* compute photosynthesis with units of mg m-2 s-1 and pass out as pointers A_mg = APHOTO * 44 / 1000 */ *A_molept = aphoto; *resppt=rd; *cipnt=ci; *wcpnt=wc; *wjpnt=wj; /* printf(" cs ci gs_leaf_mole CA ci/CA APS root1 root2 root3\n"); printf(" %5.1f %5.1f %6.3f %5.1f %6.3f %6.3f %6.3f %6.3f %6.3f\n", cs, ci, gs_leaf_mole, cca, ci / cca,aphoto,root1, root2, root3 ); */ return; } double GAMMAF(double x) { // gamma function double y,gam; gam= (1.0 / (12.0 * x)) + (1.0 / (288.0 * x*x)) - (139.0 / (51840.0 * pow(x,3.0))); gam = gam + 1.0; if (x > 0) y = sqrt(2.0 * PI / x) * pow(x,x) * exp(-x) * gam; else printf("gamma/n"); return y; } double TBOLTZ (double rate, double eakin, double topt, double tl) { // Boltzmann temperature distribution for photosynthesis double y, dtlopt,prodt,numm,denom; dtlopt = tl - topt; prodt = rugc * topt * tl; numm = rate * hkin * exp(eakin * (dtlopt) / (prodt)); denom = hkin - eakin * (1.0 - exp(hkin * (dtlopt) / (prodt))); y = numm / denom; return y; } double TEMP_FUNC(double rate,double eact,double tprime,double tref, double t_lk) { // Arhennius temperature function double y; y = rate * exp(tprime * eact / (tref * rugc*t_lk)); return y; } void ANGLE() { // ANGLE computes solar elevation angles, // This subroutine is based on algorithms in Walraven. 1978. Solar Energy. 20: 393-397 double theta_angle,G,EL,EPS,sin_el,A1,A2,RA; double delyr,leap_yr,T_local,time_1980,leaf_yr_4, delyr4; double day_savings_time, day_local; double S,HS,phi_lat_radians,value,declination_ang,ST,SSAS; double E_ang, zenith, elev_ang_deg, cos_zenith; double radd = .017453293; double twopi = 6.28318; // Oak Ridge, TN double latitude = 35.95; // latitude double longitude= 84.28; // longitude // Eastern Standard TIME double zone = 5.0; // Five hour delay from GMT delyr = time_var.year - 1980.0; delyr4=delyr/4.0; leap_yr=fmod(delyr4,4.0); day_savings_time=0.0; // Daylight Savings Time, Dasvtm =1 // Standard time, Dasvtm= 0 T_local = time_var.local_time; day_local = time_var.days; time_1980 = delyr * 365.0 + leap_yr + (day_local - 1) + T_local / 24.0; leaf_yr_4=leap_yr*4.0; if(delyr == leaf_yr_4) time_1980 -= 1.0; if(delyr < 0.0) { if (delyr < leaf_yr_4 || delyr > leaf_yr_4) time_1980 -= 1.0; } theta_angle = (360.0 * time_1980 / 365.25) * radd; G = -.031272 - 4.53963E-7 * time_1980 + theta_angle; EL = 4.900968 + 3.6747E-7 * time_1980 + (.033434 - 2.3E-9 * time_1980) * sin(G) + .000349 * sin(2. * G) + theta_angle; EPS = .40914 - 6.2149E-9 * time_1980; sin_el = sin(EL); A1 = sin_el * cos(EPS); A2 = cos(EL); RA = atan(A1/A2); /* for ATAN2 RA = atan2(A1,A2); if(RA < 0) RA=RA+twopi; */ /* The original program was in FORTRAN. It used the ATAN2 function. In C we must find the correct quadrant and evaluate the arctan correctly. Note ATAN2 is -PI TO PI, while ATN is from PI/2 TO -PI/2 */ // QUAD II, TAN theta_angle < 0 if(A1 > 0) { if(A2 <= 0) RA += 3.1415; } // QUAD III, TAN theta_angle > 0 / if(A1 <= 0) { if(A2 <= 0) RA += 3.14159; } value = sin_el * sin(EPS); if (1.-value * value >= 0) declination_ang = atan(value/ sqrt(1. - value * value)); else printf(" bad declination_ang\n"); // declination_ang=asin(value) ST = 1.759335 + twopi * (time_1980 / 365.25 - delyr) + 3.694E-7 * time_1980; if(ST >= twopi) ST = ST - twopi; S = ST - longitude * radd + 1.0027379 * (zone - day_savings_time + T_local) * 15. * radd; if(S >= twopi) S = S - twopi; HS = RA - S; phi_lat_radians = latitude * radd; // DIRECTION COSINE SSAS = (sin(phi_lat_radians) * sin(declination_ang) + cos(phi_lat_radians) * cos(declination_ang) * cos(HS)); if(1. - SSAS * SSAS >=0) E_ang = atan(sqrt(1. - (SSAS * SSAS))/ SSAS); else printf(" bad SSAS \n"); if(SSAS < 0) E_ang=E_ang+3.1415; // E=asin(SSAS); if(E_ang < 0) E_ang=3.1415/2.; zenith = E_ang / radd; elev_ang_deg = 90. - zenith; solar.beta_rad = elev_ang_deg * radd; solar.sine_beta = sin(solar.beta_rad); cos_zenith = cos(zenith); solar.beta_deg = solar.beta_rad / PI180; return; } void IRFLUX(double x) { met.ir_in = SKY_IR(x); return; } void DIFFUSE_DIRECT_RADIATION() { double fand, fir,fv; double rdir, rdvis,rsvis,wa; double ru, rsdir, rvt,rit, nirx; double xvalue,fvsb,fvd,fansb; /* This subroutine uses the Weiss-Norman ( 1985, Agric. forest Meteorol. 34: 205-213) routine tocompute direct and diffuse PAR from total par fractions of NIR and PAR (visible) */ fir = .54; fv = .46; ru = 98.5 / (101.3 * solar.sine_beta); // visible direct PAR rdvis = 600.0 * exp(-.185 * ru) * solar.sine_beta; // potential diffuse PAR rsvis = .4 * (600.0 - rdvis) * solar.sine_beta; /* solar constant: 1320 W m-2 water absorption in NIR for 10 mm precip water */ wa = 1320.0 * .077 * pow((2. * ru),.3); /* direct beam NIR */ rdir = (720.0 * exp(-.06 * ru) - wa) * solar.sine_beta; if(rdir < 0) rdir=0; /* potential diffuse NIR */ rsdir = .6 * (720 - wa - rdir) * solar.sine_beta; if (rsdir < 0) rsdir=0; rvt = rdvis + rsvis; rit = rdir + rsdir; if(rit <= 0) rit = .1; if(rvt <= 0) rvt = .1; solar.ratrad = input.rglobal / (rvt + rit); if (time_var.local_time >= 12.00 && time_var.local_time <=13.00) solar.ratradnoon=solar.ratrad; /* ratio is the ratio between observed and potential radiation NIR flux density as a function of PAR since NIR is used in energy balance calculations convert it to W m-2: divide PAR by 4.6 */ nirx = input.rglobal - (input.parin / 4.6); // ratio = (PARIN / 4.6 + NIRX) / (rvt + rit) if (solar.ratrad >= .9) solar.ratrad = .89; if (solar.ratrad <= 0) solar.ratrad=0.22; // fraction PAR direct and diffuse xvalue=(0.9-solar.ratrad)/.70; fvsb = rdvis / rvt * (1. - pow(xvalue,.67)); if(fvsb < 0) fvsb = 0.; if (fvsb > 1) fvsb=1.0; fvd = 1. - fvsb; // note PAR has been entered in units of uE m-2 s-1 solar.par_beam = fvsb * input.parin; solar.par_diffuse = fvd * input.parin; if(solar.par_beam <= 0) { solar.par_beam = 0; solar.par_diffuse = input.parin; } if(input.parin == 0) { solar.par_beam=0.001; solar.par_diffuse=0.001; } xvalue=(0.9-solar.ratrad)/.68; fansb = rdir / rit * (1. - pow(xvalue,.67)); if(fansb < 0) fansb = 0; if(fansb > 1) fansb=1.0; fand = 1. - fansb; // NIR beam and diffuse flux densities solar.nir_beam = fansb * nirx; solar.nir_diffuse = fand * nirx; if(solar.nir_beam <= 0) { solar.nir_beam = 0; solar.nir_diffuse = nirx; } if (nirx == 0) { solar.nir_beam=0.1; solar.nir_diffuse=0.1; } return; } double SFC_VPD (double tlk, double *leleafpt) { // this function computes the relative humidity at the leaf surface for // application in the Ball Berry Equation // latent heat flux, LE, is passed through the function, mol m-2 s-1 // and it solves for the humidity at leaf surface double y, rhov_sfc,e_sfc,vpd_sfc,rhum_leaf; double es_leaf; es_leaf = ES(tlk); // saturation vapor pressure at leaf temperature rhov_sfc = (*leleafpt / (fact.latent)) * bound_layer_res.vapor + met.rhova_kg; /* kg m-3 */ e_sfc = rhov_sfc * tlk / .2165; // mb vpd_sfc = es_leaf - e_sfc; // mb rhum_leaf = 1. - vpd_sfc / es_leaf; // 0 to 1.0 y = rhum_leaf; return y; } void LEAF_OPTICS() { solar.par_reflect = .133; solar.par_trans = .072; solar.par_soil_refl = .1; solar.par_absorbed = (1. - solar.par_reflect - solar.par_trans); // NIR wave band: after Norman (1979) and NASA report solar.nir_reflect = .43; solar.nir_trans = .26; solar.nir_soil_refl = .1; // Absorbed NIR solar.nir_absorbed = (1. - solar.nir_reflect - solar.nir_trans); return; } void BOUNDARY_RESISTANCE() { /* ************************************** BOUNDARY_RESISTANCE This subroutine computes the leaf boundary layer resistances for heat, vapor and CO2 (s/m). Flat plate theory is used, as discussed in Schuepp (1993) and Grace and Wilson (1981). We consider the effects of turbulent boundary layers and sheltering. Schuepp's review shows a beta factor multiplier is necessary for SH in flows with high turbulence. The concepts and theories used have been validated with our work on HNO3 transfer to the forest. Schuepp. 1993 New Phytologist 125: 477-507 Diffusivities have been corrected using the temperature/Pressure algorithm in Massman (1998) */ double Re,Re5,Re8; // Reynolds numbers double Sh_heat,Sh_vapor,Sh_CO2; // Sherwood numbers double graf, GR25; // Grasshof numbers double deltlf; double Res_factor; double nnu_T_P, ddh_T_P, ddv_T_P, ddc_T_P, T_kelvin; // Difference between leaf and air temperature deltlf = leaf.T_sfc_C - input.ta; T_kelvin=input.ta + 273.16; if(deltlf > 0) graf = non_dim.grasshof * deltlf / T_kelvin; else graf=0; nnu_T_P=nnu*(1013./input.press_mb)*pow((T_kelvin/273.16),1.81); Re = lleaf * input.wnd / nnu_T_P; if (Re > 0.) Re5 = sqrt(Re); else { printf("bad RE in RESHEAT\n"); Re5=100.; } Re8 = pow(Re,.8); if( Re > 14000.) { Res_factor=0.036*Re8*betfact; /* turbulent boundary layer SH = .036 * Re8 * pr33*betfact; SHV = .036 * Re8 * sc33*betfact; SHCO2 = .036 * Re8 * scc33*betfact; */ Sh_heat = Res_factor * non_dim.pr33; Sh_vapor = Res_factor * non_dim.sc33; Sh_CO2 = Res_factor * non_dim.scc33; } else { Res_factor=0.66*Re5*betfact; /* laminar sublayer SH = .66 * Re5 * pr33*betfact; SHV = .66 * Re5 * sc33*betfact; SHCO2 = .66 * Re5 * scc33*betfact; */ Sh_heat = Res_factor * non_dim.pr33; Sh_vapor = Res_factor * non_dim.sc33; Sh_CO2 = Res_factor * non_dim.scc33; } // If there is free convection if(graf / (Re * Re) > 1.) { // Compute Grashof number for free convection if (graf < 100000.) GR25 = .5 * pow(graf,.25); else GR25 = .13 * pow(graf,.33); Sh_heat = non_dim.pr33 * GR25; Sh_vapor = non_dim.sc33 * GR25; Sh_CO2 = non_dim.scc33 * GR25; } // lfddx=lleaf/ddx // Correct diffusivities for temperature and pressure ddh_T_P=ddh*(1013./input.press_mb)*pow((T_kelvin/273.16),1.81); ddv_T_P=ddv*(1013./input.press_mb)*pow((T_kelvin/273.16),1.81); ddc_T_P=ddc*(1013./input.press_mb)*pow((T_kelvin/273.16),1.81); bound_layer_res.heat = lleaf/(ddh_T_P * Sh_heat); bound_layer_res.vapor = lleaf/(ddv_T_P * Sh_vapor); bound_layer_res.co2 = lleaf / (ddc_T_P * Sh_CO2); return; }