Skip to content

refactor(LIU_AN_ET): canopy scaling fixes + full documentation #81

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
braghiere wants to merge 1 commit into
base: main
Choose a base branch
from
Open

refactor(LIU_AN_ET): canopy scaling fixes + full documentation #81

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
297 changes: 157 additions & 140 deletions C/projects/CARDAMOM_MODELS/DALEC/DALEC_ALL/LIU_AN_ET_REFACTOR.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,201 +64,218 @@ double leaf_refl_nir=A->IN.leaf_refl_nir;
double maxPevap=A->IN.maxPevap;
double precip=A->IN.precip;

//CONSTS
double Ephoton = 2.0e-25/500.0e-9;// Planck constant times speed of light (J.s*m.s-1) divided by light wavelength (m)
double NA = 6.02e23;// Avogadro's constant /mol
double R = 8.31e-3;//Gas constant, kJ/mol/K
double lambda0 = 2.26e6;
double gammaV = 100*1005/(lambda0*0.622);


//DATA contains all met forcings

double PAR; //incident PAR (?) review
double T_C; // air temperature in C
double Kc; //Michaelis-Menten parameter temperature dependent
double Ko; //Michaelis-Menten parameter temperature dependent
double cp; //specific heat capacity of air

double Vcmax; //the maximum rate of carboxylation of Rubisco
double Jmax; //the maximum rate of electron transport
double J;
double medlyn_term;
double ci;
/* --- Physical and thermodynamic constants --- */
double Ephoton = 2.0e-25/500.0e-9; // Energy per photon at 500 nm: h*c/lambda (J); h=6.626e-34 J.s, c=3e8 m/s
double NA = 6.02e23; // Avogadro's constant (mol-1)
double R = 8.31e-3; // Universal gas constant (kJ mol-1 K-1)
double lambda0 = 2.26e6; // Latent heat of vaporization of water (J kg-1)
double rho_air = 1.225; // Dry air density at sea level (kg m-3)
double Cp_air = 1005.; // Specific heat of dry air at constant pressure (J kg-1 K-1)
double eps_mol = 0.622; // Ratio of molar mass of water vapour to dry air: Mw/Md = 18.015/28.97 (-)
/* Psychrometric constant gamma = rho_air * Cp_air / (lambda0 * eps_mol) (Pa K-1);
here expressed as a dimensionless coefficient for use with kPa and m s-1 units */
double gammaV = 100.*Cp_air/(lambda0*eps_mol); // ~0.067 kPa K-1; Monteith & Unsworth 2008


/* --- FvCB photosynthesis state variables --- */
double PAR; // Incident PAR at canopy top (umol photons m-2 s-1)
double T_C; // Air temperature (deg C)
double Kc; // Michaelis-Menten constant for CO2 (umol mol-1), temperature-dependent
double Ko; // Michaelis-Menten constant for O2 (mmol mol-1), temperature-dependent
double cp; // CO2 compensation point in absence of Rd, Gamma* (umol mol-1)
double Vcmax; // Maximum rate of Rubisco carboxylation at T_C (umol CO2 m-2 s-1)
double Jmax; // Maximum rate of electron transport at T_C (umol e- m-2 s-1)
double J; // Actual electron transport rate (umol e- m-2 s-1)
double medlyn_term; // Medlyn optimal stomatal model factor: 1 + g1/sqrt(VPD); Medlyn et al. 2011
double ci; // Leaf intercellular CO2 concentration (umol mol-1)
//double C3_frac = 1.; //Fraction of C3 per gridcell
double a1;
double a2;

double Ag_C3;
double Ag_C4;

double Ag;
double An;

double Rd_C3;
double Rd_C4;

double Rd;

//Array of results with 3 positions for An, E, and Transpiration;
//static double r[3];


PAR = SRAD/(2*Ephoton*NA)*1e6;


double canopy_scale = (1. - exp(-VegK*LAI*clumping))/(VegK);
double a1; // Rubisco-limited (co-limited) gross assimilation term (umol CO2 m-2 s-1)
double a2; // Electron-transport-limited gross assimilation term (umol CO2 m-2 s-1)

double Ag_C3; // Gross assimilation, C3 pathway (leaf-level, umol CO2 m-2 s-1)
double Ag_C4; // Gross assimilation, C4 pathway (leaf-level, umol CO2 m-2 s-1)
double Ag; // Blended gross assimilation (leaf-level, umol CO2 m-2 s-1)
double An; // Net assimilation = Ag - Rd (leaf-level, umol CO2 m-2 s-1)

double Rd_C3; // Dark respiration, C3 (umol CO2 m-2 s-1)
double Rd_C4; // Dark respiration, C4 (umol CO2 m-2 s-1)
double Rd; // Blended dark respiration (umol CO2 m-2 s-1)

/* Convert incoming shortwave to incident PAR.
SRAD is total shortwave (W m-2); ~50% is PAR (400-700 nm); dividing by Ephoton*NA
converts W m-2 (= J s-1 m-2) to mol photons m-2 s-1, then *1e6 -> umol photons m-2 s-1.
The factor of 2 accounts for the ~50% PAR fraction of shortwave; Liu et al. myFun.py. */
PAR = SRAD/(2.*Ephoton*NA)*1e6;

/* --- Canopy radiation interception --- */
/* canopy_scale: integral of Beer's law profile int_0^LAI exp(-k*Omega*l) dl
= (1 - exp(-k*Omega*LAI)) / (k*Omega); denominator requires both k AND Omega (clumping);
de Pury & Farquhar (1997) Plant Cell Environ. 20:537; Bonan (2019) Ch. 14 Eq. 14.6 */
double canopy_scale = (1. - exp(-VegK*LAI*clumping))/(VegK*clumping); // Beer's law integral /(k*Omega), not /(k); de Pury & Farquhar 1997; Bonan Ch. 14 Eq. 14.6

//printf(" \n Initial LAI %f \n ", LAI);

//absorbed PAR assuming black canopy.
//PAR = PAR*(1. - exp(-LAI*VegK));

PAR *= (1. - leaf_refl_par)*(1. - exp(-VegK*LAI*clumping));
/* PAR_leaf: effective leaf-level absorbed PAR for the big-leaf model.
Derived as APAR_canopy / canopy_scale = I_0*(1-rho)*k*Omega; this keeps
all FvCB quantities (a1, a2) at leaf level so canopy_scale is applied only once at output.
Bonan (2019) Ch. 14; Liu et al. myFun.py f_carbon() */
double PAR_leaf = PAR * (1. - leaf_refl_par) * VegK * clumping; // leaf-level PAR = APAR/canopy_scale = I_0*(1-rho)*k*Omega; Bonan Ch. 14


T_C = TEMP - DGCM_TK0C; // Convert temperature to degrees C


Kc = 300.*exp(0.074*(T_C - 25.));
Ko = 300.*exp(0.015*(T_C - 25.));
cp = 36.9 + 1.18*(T_C - 25.) + 0.036*pow((T_C - 25.), 2.);
/* Temperature-dependent Michaelis-Menten kinetics; Q10-type fits to Bernacchi et al. (2001) data.
Kc ref = 300 umol mol-1, Ko ref = 300 mmol mol-1, both at 25 C. */
Kc = 300.*exp(0.074*(T_C - 25.)); // Michaelis constant for CO2 (umol mol-1)
Ko = 300.*exp(0.015*(T_C - 25.)); // Michaelis constant for O2 (mmol mol-1)
cp = 36.9 + 1.18*(T_C - 25.) + 0.036*pow((T_C - 25.), 2.); // CO2 compensation point Gamma* (umol mol-1); Bernacchi et al. 2001; Liu et al. myFun.py


//Vcmax = vcmax25*exp(50.*(TEMP - 298.)/(298.*R*TEMP));
/* --- Vcmax temperature response (Q10 formulation with upper/lower inhibition) ---
Vcmax(T) = Vcmax25 * Q10^(0.1*(T-25)) / [(1+exp(0.3*(T-Tupp))) * (1+exp(0.3*(Tdown-T)))]
Tupp, Tdown: upper/lower thermal inhibition thresholds (K); inverted by MCMC. */
double q_10 = A->IN.q10canopy;
double fT = pow(A->IN.q10canopyRd,(T_C-25.)/10.); // Q10 temperature scaling for Rd; reference T = 25 C

double fT = pow(A->IN.q10canopyRd,(T_C-25.)/10); // reference temperature is 25 degrees C



Vcmax = vcmax25*pow(q_10,0.1*(T_C-25.))/((1 + exp(0.3*(T_C-(Tupp-DGCM_TK0C))))*(1 +exp(0.3*((Tdown-DGCM_TK0C)-T_C))));
Vcmax = vcmax25*pow(q_10,0.1*(T_C-25.))/((1. + exp(0.3*(T_C-(Tupp-DGCM_TK0C))))*(1. + exp(0.3*((Tdown-DGCM_TK0C)-T_C))));

/* Jmax/Vcmax ~ e ≈ 2.72 at 25 C; Leuning (1997) J. Exp. Bot. 48:345
(Wullschleger 1993 109-species dataset rescaled to 25 C: slope = 2.68, r2 = 0.87)
Fixed ratio minimises MCMC parameter space; temperature response of Jmax tracks Vcmax. */
Jmax = Vcmax*exp(1.);

J = (0.3*PAR + Jmax - sqrt(pow(0.3*PAR + Jmax,2) - 4.*0.9*0.3*PAR*Jmax))/2./0.9;

/* Electron transport rate J: non-rectangular hyperbola (Farquhar & Wong 1984).
kai1 = 0.9 (curvature parameter), kai2 = 0.3 (quantum yield, mol e- per mol photon);
Bonan (2019) Ch. 11 Eq. 11.33; Liu et al. myFun.py */
J = (0.3*PAR_leaf + Jmax - sqrt(pow(0.3*PAR_leaf + Jmax,2) - 4.*0.9*0.3*PAR_leaf*Jmax))/2./0.9; // PAR_leaf keeps a2 leaf-level (consistent with a1); single canopy_scale at output; Bonan Ch. 11

/* Medlyn et al. (2011) optimal stomatal conductance model: gs/g0 ~ 1 + g1/sqrt(VPD).
ci derived analytically from gs = An/(co2-ci)*1.6, giving ci = co2*(1 - 1/medlyn_term).
Medlyn et al. (2011) Glob. Change Biol. 17:2134 */
medlyn_term = 1. + g1/sqrt(VPD);

ci = co2*(1. - 1./medlyn_term);

//CLM 4.5 does this
if (ci<cp){ci=cp;}

//Two terms for C3 photosythesis
a1 = Vcmax*(ci - cp)/(ci + Kc*(1.+209./Ko));
a2 = J*(ci-cp)/(4.*(ci + 2.*cp));



// An_C3 = fmax(0., fmin(a1*beta_factor,a2) - 0.015*Vcmax*beta_factor);
Ag_C3 = fmin(a1*beta_factor,a2);
Rd_C3 = A->IN.canopyRdsf*vcmax25*fT;



//Two terms for C4 photosynthesis
if (ci<cp){ci=cp;} // CLM 4.5: ci cannot drop below CO2 compensation point

/* --- C3 photosynthesis: FvCB co-limitation (Farquhar, von Caemmerer & Berry 1980) ---
a1: Rubisco-limited rate; 209 umol mol-1 = ambient O2 concentration (Coa);
a2: RuBP-regeneration (electron-transport) limited rate;
beta_factor applied to a1 only: soil moisture suppresses Rubisco carboxylation
(stomatal + non-stomatal limitation); light-limited a2 is not gated by stomata;
Bonan (2019) Ch. 11 Eqs. 11.28-11.29 */
a1 = Vcmax*(ci - cp)/(ci + Kc*(1.+209./Ko)); // Rubisco-limited (umol CO2 m-2 s-1); 209 = O2 conc. in mmol/mol; Bonan Ch. 11 Eq. 11.28
// beta_factor applied to a1 only; a2 (light-limited) unconstrained by soil moisture
a2 = J*(ci-cp)/(4.*(ci + 2.*cp)); // RuBP-limited (umol CO2 m-2 s-1); Bonan Ch. 11 Eq. 11.29

Ag_C3 = fmin(a1*beta_factor, a2); // C3 gross assimilation: co-limited minimum
Rd_C3 = A->IN.canopyRdsf*vcmax25*fT; // C3 dark respiration = canopyRdsf * Vcmax25 * Q10(T)

/* --- C4 photosynthesis: simplified two-term co-limitation ---
C4 plants have a CO2-concentrating mechanism (CCM), so Rubisco is not CO2-limited;
a1 ~ Vcmax (enzyme-limited), a2 ~ J (light-limited).
C3_frac is prescribed per pixel from remote-sensing products, not inverted by MCMC. */
a1 = Vcmax;
a2 = J;

// An_C4 = fmax(0., fmin(a1*beta_factor,a2) - 0.015*Vcmax*beta_factor);
Ag_C4 = fmin(a1*beta_factor,a2);
Rd_C4 =A->IN.canopyRdsf*vcmax25*fT;

//Total photosynthesis
Ag_C4 = fmin(a1*beta_factor, a2); // C4 gross assimilation
Rd_C4 = A->IN.canopyRdsf*vcmax25*fT; // C4 dark respiration

Ag = C3_frac*(Ag_C3) + (1. - C3_frac)*(Ag_C4);
Rd = C3_frac*(Rd_C3) + (1. - C3_frac)*(Rd_C4);
/* Fractional C3/C4 mix (C3_frac prescribed, not MCMC-inverted) */
Ag = C3_frac*(Ag_C3) + (1. - C3_frac)*(Ag_C4); // blended gross assimilation (leaf-level, umol CO2 m-2 s-1)
Rd = C3_frac*(Rd_C3) + (1. - C3_frac)*(Rd_C4); // blended dark respiration (leaf-level, umol CO2 m-2 s-1)

//Potential Rd
/* Potential canopy Rd (g C m-2 day-1):
leaf-level Rd * canopy_scale (leaf->canopy) * 12e-6 (g C per umol CO2) * 86400 (s day-1) */
double Rd_daily_potential = Rd*canopy_scale*(12.e-6)*(24.*60.*60.);

//Ensures NSCs available
//A->OUT.LEAF_MORTALITY_FACTOR=fmax( 1- A->IN.NSC/(Rd_daily_potential * A->IN.deltat) ,0);
/* Leaf mortality from NSC (non-structural carbohydrates) limitation.
Exponential formulation: LEAF_MORTALITY_FACTOR -> 0 when NSC >> Rd*deltat (no stress),
-> 1 when NSC -> 0 (full mortality); updated to exponential form 11/1/2023. */
if (Rd_daily_potential==0){
A->OUT.LEAF_MORTALITY_FACTOR=0;
} else {
A->OUT.LEAF_MORTALITY_FACTOR=(1/(exp(A->IN.NSC/(Rd_daily_potential * A->IN.deltat)))); // new exponential formulation 11/1/2023
A->OUT.LEAF_MORTALITY_FACTOR=(1./(exp(A->IN.NSC/(Rd_daily_potential * A->IN.deltat))));
}
// printf(" \n Rd_pot and NSC: %f %f \n ", Rd_daily_potential * A->IN.deltat, A->IN.NSC);
//Actual daily Rd
A->OUT.Rd =Rd_daily_potential*(1 - A->OUT.LEAF_MORTALITY_FACTOR);

//Scaling Rd to available NSCs
Rd = Rd *(1 - A->OUT.LEAF_MORTALITY_FACTOR);
/* Actual daily canopy Rd, scaled by leaf survival fraction */
A->OUT.Rd = Rd_daily_potential*(1. - A->OUT.LEAF_MORTALITY_FACTOR);

/* Scale leaf-level Rd by survival fraction before computing leaf-level An */
Rd = Rd*(1. - A->OUT.LEAF_MORTALITY_FACTOR);

//Net A
/* Net assimilation (leaf-level); both Ag and Rd are leaf-level here.
canopy_scale applied once below when writing to output. */
An = Ag - Rd;



//To scale from leaf to canopy, comment out the following line and uncomment the one after
//double canopy_scale = 1.;

//r[0] = An*canopy_scale*(12.e-6)*(24.*60.*60.); //from umolCO2m-2s-1 to gCm-2day-1
/* Scale leaf -> canopy and convert units: umol CO2 m-2 s-1 -> g C m-2 day-1
12e-6: molecular mass of C in g per umol CO2 (12 g mol-1 * 1e-6 mol umol-1)
86400: seconds per day */
//double canopy_scale = 1.; // uncomment to disable canopy scaling (diagnostic)
A->OUT.Ag = Ag*canopy_scale*(12.e-6)*(24.*60.*60.);
A->OUT.An = An*canopy_scale*(12.e-6)*(24.*60.*60.);

//##################Transpiration#################

double SRADg;
double sV; // sV is the rate of change of saturated vapor pressure with air temperature in degC.
double petVnum;
double petVnumB;
double gs; //stomatal conductance
double transp; // transpiration
double evap; // evaporation
double Psurf = 100.0; //Surface pressure in kPa
double VPD_kPa = VPD;//*Psurf; //100.0 kPa = 1000.0 hPa => Surface pressure
/* ===================== Transpiration: Penman-Monteith (Monteith 1965) =====================
Canopy transpiration E = [Delta*(Rn_veg) + rho_air*Cp*(VPD)*ga] / [lambda*(Delta + gamma*(1 + ga/gs))]
where Delta = d(esat)/dT (kPa K-1), gamma = psychrometric constant (kPa K-1),
ga = aerodynamic conductance (m s-1), gs = canopy stomatal conductance (m s-1).
Reference: Monteith & Unsworth (2008) "Principles of Environmental Physics", Ch. 13;
Liu et al. myFun.py PenmanMonteith(). */

double SRADg; // Shortwave radiation reaching the ground surface (W m-2)
double sV; // Slope of saturation vapour pressure curve, d(esat)/dT (kPa K-1)
double petVnum; // PM numerator for canopy transpiration (W m-2 equivalent)
double petVnumB; // PT numerator for ground evaporation (Priestley-Taylor 1972)
double gs; // Canopy stomatal conductance (m s-1)
double transp; // Canopy transpiration (mm hr-1)
double evap; // Ground evaporation (mm hr-1)
double VPD_kPa = VPD; // VPD already in kPa
double evap_scale_factpr;

sV = 0.04145*exp(0.06088*T_C);
/* Slope of saturation vapour pressure curve (kPa K-1); empirical fit to Magnus equation */
sV = 0.04145*exp(0.06088*T_C);

/* Canopy and ground radiation partitioning using Beer's law:
SRADg = below-canopy shortwave (W m-2); (1-0.5*(rho_PAR+rho_NIR)) = broadband canopy absorptance */
SRADg = (1. - 0.5*(leaf_refl_par+leaf_refl_nir))*SRAD*exp(-VegK*LAI*clumping);
SRAD = (1. - 0.5*(leaf_refl_par+leaf_refl_nir))*SRAD; // absorbed canopy shortwave (W m-2)

SRAD = (1. - 0.5*(leaf_refl_par+leaf_refl_nir))*SRAD;

petVnum = (sV*(SRAD-SRADg)+1.225*1005*VPD_kPa*ga)/lambda0*60*60;
petVnumB = 1.26*(sV*SRADg)/(sV+gammaV)/lambda0*60*60; //from mm.hr-1
if(beta_factor > 0 && SRAD >0){
/* PM numerator for canopy transpiration (mm hr-1):
numerator = [Delta*(Rn_canopy) + rho_air*Cp*VPD*ga] / lambda * 3600
rho_air = 1.225 kg m-3; Cp_air = 1005 J kg-1 K-1; lambda0 = 2.26e6 J kg-1 */
petVnum = (sV*(SRAD-SRADg) + rho_air*Cp_air*VPD_kPa*ga)/lambda0*60.*60.; // mm hr-1

//Option 1. gs = 1.6*Ag/(co2-ci)*LAI*0.02405;
//Option 1. gs = 1.6*An/(co2-ci)*LAI*0.02405;
gs = fmax(0,1.6*An/(co2-ci)*LAI*0.02405);
/* Priestley-Taylor ground evaporation (mm hr-1):
alpha_PT = 1.26 (empirical coefficient; Priestley & Taylor 1972 Mon. Wea. Rev. 100:81);
uses below-canopy radiation SRADg; assumes soil evaporation ~ equilibrium */
petVnumB = 1.26*(sV*SRADg)/(sV+gammaV)/lambda0*60.*60.; // mm hr-1; alpha_PT=1.26; Priestley & Taylor 1972

//transp = petVnum/(sV+gammaV*(1+ga*(1/ga+1/gs)));
transp = petVnum/(sV+gammaV*(ga*(1/ga+1/gs)));
if(beta_factor > 0 && SRAD > 0){

/* Canopy stomatal conductance to water vapour (m s-1):
Based on leaf-level gs_CO2 = An/(co2-ci) (mol CO2 m-2 s-1 per umol mol-1);
1.6: ratio of diffusivity of H2O to CO2 in air (dimensionless; Jones 1992 "Plants and Microclimate");
LAI: leaf-to-canopy upscaling (m2 leaf m-2 ground);
0.02405: molar volume of air at ~20 C = RT/P = 8.314*293/101325 = 0.02406 m3 mol-1,
converts mol CO2 m-2 s-1 -> m s-1. */
gs = fmax(0., 1.6*An/(co2-ci)*LAI*0.02405);

/* Penman-Monteith transpiration (mm hr-1).
Denominator: sV + gamma*(1 + ga/gs) = sV + gamma*ga*(1/ga + 1/gs); Monteith 1965.
Liu et al. myFun.py PenmanMonteith() */
transp = petVnum/(sV + gammaV*(ga*(1./ga + 1./gs)));

}
else{
transp = 0.0;
} else {
transp = 0.0;
}

//r[1] = transp*24; //from mm.hr-1 to mm.day-1

A->OUT.transp = transp*24;


A->OUT.transp = transp*24.; // mm hr-1 -> mm day-1

/* Ground evaporation: scale by precipitation availability (evap <= precip/maxPevap) */
evap_scale_factpr = fmin(precip/maxPevap, 1.);

evap = petVnumB*evap_scale_factpr;

//evap = petVnumB;

//r[2] = evap*24; //from mm.hr-1 to mm.day-1

A->OUT.evap = evap*24;



//printf("SRADg = %f, VegK = %f, LAI = %f\n", (SRADg/SRAD), VegK, LAI);
//printf("transp = %f, evap = %f\n", r[1], r[2] );
A->OUT.evap = evap*24.; // mm hr-1 -> mm day-1

return 0;
}
Loading