!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine leafanetmodel(ilimittype,aPPFDlf,templeaf,pco2i,po2i, &chlflphips2,gascon,resistwp25,resistch25,ha_gmeso,hd_gmeso, &sv_gmeso,vcmax25,ha_vcmax,hd_vcmax,sv_vcmax,fjmax25,ha_jmax, &hd_jmax,sv_jmax,tpu25,ha_tpu,hd_tpu,sv_tpu,alpha,rdlight25, &ha_darkresp,stargamma25,ha_stargamma,fkc25,ha_kc,fko25,ha_ko, &phifactor,thetafactor,betaPSII,assim_net,rdlight,iphotolimit,atp, &resistwp,resistch,stargamma,pco2c,realizedfjelect,assim_net_flu, &pco2c_flu) implicit none ! !------------ Inputs ------------------- !ilimittype=1: Rubisco,RuBp and TPU limitations !ilimittype=2: Rubisco and RuBp limitations only !ilimittype=3: Rubisco and TPU limitations only !ilimittype=4: RuBp and TPU limitations only !ilimittype=5: Rubisco limitation only !ilimittype=6: RuBp limitation only !ilimittype=7: TPU limitation only ! aPPFDlf: absorbed photosynthetic photon flux density by leaf (umol m-2 s-1) ! templeaf: leaf temperature [K] ! resistwp25: resistance to CO2 via cell walls and plasmalemma [umol-1m2sPa] ! resistch25: resistance to CO2 via chloroplastic envelope [umol-1m2sPa] ! vcmax25: maximum RuBP saturated rate of carboxylation at 25oC ! of leaf temperature [umol m-2 s-1] ! fjmax25: maximum rate of electron transport at 25oC ! of leaf temperature [umol m-2 s-1] ! rdlight25: Mitochondrial respiration rate in the light at 25oC ! tpu25: tpu at 25oC, [umol m-2 s-1] ! alpha: The fraction of glycolate carbon not returned to the chloroplat ! pco2i: intercellular CO2 partial pressure (Pa). ! po2i: intercellular O2 partial pressure (Pa, often taking the ambient value). ! chlflphips2: Photochemical efficiency of photosystem II measured with fluorescence, if provided (NA) ! gascon: universal gas constant [JK-1 mol-1]. ! ha_vcmax: Ha in Vcmax temperature response function ~73637.0d0/1000.0d0 ! hd_vcmax: Hd in Vcmax temperature response function ~149252.0d0/1000.0d0 ! sv_vcmax: Sv in Vcmax temperature response function ~486.0d0/1000.0d0 ! ha_jmax: Ha in Jmax temperature response function ~50300.0d0/1000.0d0 ! hd_jmax: Hd in Jmax temperature response function ~152044.0d0/1000.0d0 ! sv_jmax: Sv in Jmax temperature response function ~495.0d0/1000.0d0 ! ha_tpu: ha in tpu temperature response function ~53100/1000 ! hd_tpu: hd in tpu temperature response function ~201800/1000 ! sv_tpu: Sv in tpu temperature response function ~650/1000 ! ha_darkresp: Mitochondrial respiation parameter deltaha = 46.39d0 ! stargamma25: the compensation point at 25oC [Pa] ! ha_stargamma: parameter in the temperature response function, typically around ! 37.83 kJmol-1 ! fkc25: the Michaelis constant for CO2 at 25oC ! ha_kc: a CO2 Michaelis temp coefficient (~79.43 KJmol-1) ! fko25: the Michaelis constant for O2 at 25oC ! ha_ko: a O2 Michaelis temp coefficient (~36.38 KJmol-1) ! ha_gmeso: The activation energy (Ha) in gmeso temperature response function ~49.6d0[kJmol-1] ! hd_gmeso: The deactivation energy (Hd) in gmeso temperature response function ~437.4d0[kJmol-1] ! sv_gmeso: The entropy term (Sv) in gmeso temperature response function ~1.4d0 [kJmol-1K-1] ! phifactor: modifies Bernacchi phiPSIImax ! thetafactor: modifies Bernacchi thetaPSII ! betaPSII: fraction of absorbed light that reaches photosystem II (~0.5) !!3/29/2014. Bernacchi method cannot fit A/PAR curves well. phifactor and thetafactor are !used to modify his method when A/PAR curves are present. integer ilimittype double precision aPPFDlf,templeaf,vcmax25,fjmax25,rdlight25, &tpu25,alpha,pco2i,po2i,chlflphips2,gascon,ha_vcmax,hd_vcmax, &sv_vcmax,ha_jmax,hd_jmax,sv_jmax,ha_tpu,hd_tpu,sv_tpu, &ha_darkresp,stargamma25,ha_stargamma,fkc25,ha_kc,fko25,ha_ko, &resistwp25,resistch25,ha_gmeso,hd_gmeso,sv_gmeso,phifactor, &thetafactor,betaPSII !------------ Outputs ------------------- ! assim_net: net rate of CO2 uptake per unit leaf area ! [umol m-2 s-1] ! rdlight: Mitochondrial respiration rate in the light ! [umol m-2 s-1] ! iphotolimit: index of photosynthesis limiting process ! 1 = Rubisco-limited photosynthesis ! 2 = RuBP-Limited photosynthesis ! 3 = TPU-limited photosynthesis ! stargamma: CO2 compensation point (Pa) ! atp: the concentration of ATP (mmol ATP m-2) ! resistwp: resistance to CO2 via cell walls and plasmalemma [umol-1m2sPa] ! resistch: resistance to CO2 via chloroplastic envelope [umol-1m2sPa] ! realizedfjelect: the realized electron transport rate, <=fjelect from light equation (=when ! RuBP regeneration limits photosynthesis) (umolm-2s-1) ! assim_net_flu: assim_net calculated with chlflphips2 (if provided), umolm-2s-1 ! pco2c_flu: pco2c calculated from chlflphips2 (if provided), Pa. double precision assim_net,rdlight,atp,pco2c,realizedfjelect, &assim_net_flu,pco2c_flu integer iphotolimit !------------ Local variables ----------- ! wc: Rubisco limited rate of carboxylation [umol m-2 s-1] ! wj: RuBP limited rate of carboxylation [umol m-2 s-1] ! tpu: rate of triose phosphate utilization [umol m-2 s-1] ! atpcon: the concentration of ATP (mmol ATP m-2) ! fkco ! fkc: Michaelis constant for CO2 (Pa) ! fko: Michaelis constant for O2 (Pa) ! vcmax: maximum RuBP saturated rate of carboxylation [umol m-2 s-1] ! fjelect: rate of electron transport [umol m-2 s-1] double precision wc,wj,tpu,atpcon,fkco,gmeso,resistwp,resistch, &fkc,fko,vcmax,fjelect,stargamma,term integer i call vcmaxontemp(templeaf,vcmax25,gascon,ha_vcmax,hd_vcmax, &sv_vcmax,vcmax) call jontemp(aPPFDlf,templeaf,fjelect,fjmax25,ha_jmax,hd_jmax, &sv_jmax,gascon,phifactor,thetafactor,betaPSII) call tpuontemp(templeaf,gascon,tpu25,ha_tpu,hd_tpu,sv_tpu,tpu) if(resistwp25.gt.0.0d0)then call gmesoontemp(templeaf,1.0d0,gascon,ha_gmeso,hd_gmeso, &sv_gmeso,gmeso) resistwp=resistwp25/gmeso else resistwp=0.0d0 endif if(resistch25.gt.0.0d0)then call gmesoontemp(templeaf,1.0d0,gascon,ha_gmeso,hd_gmeso, &sv_gmeso,gmeso) resistch=resistch25/gmeso else resistch=0.0d0 endif call resp_mitocho(templeaf,rdlight25,ha_darkresp,gascon,rdlight) call co2compens(templeaf,stargamma25,ha_stargamma,gascon, &stargamma) call MichaelisCO2(templeaf,fkc25,ha_kc,gascon,fkc) call MichaelisO2(templeaf,fko25,ha_ko,gascon,fko) fkco=fkc*(1.0d0+po2i/fko) call Anet_Final(vcmax,fjelect,tpu,resistwp,resistch,stargamma, &fkco,pco2i,alpha,rdlight,ilimittype,iphotolimit,assim_net, &pco2c,realizedfjelect) wc=vcmax*pco2c/(pco2c+fkco) wj=fjelect*pco2c/(4.0d0*pco2c+8.0d0*stargamma) atp=atpcon(wc,wj,vcmax) assim_net_flu=-9999.0d0 pco2c_flu=-9999.0d0 if(chlflphips2.gt.0.0d0)then fjelect=betaPSII*chlflphips2*aPPFDlf call Anet_Final(vcmax,fjelect,tpu,resistwp,resistch,stargamma, &fkco,pco2i,alpha,rdlight,6,i,assim_net_flu,pco2c_flu,term) endif return end subroutine leafanetmodel !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% double precision function threeexp(fact1,fact2,fact3) implicit none ! ! This function computes ! threeexp=(dexp(fact3)+dexp(fact2+fact3))/(1.0d0+dexp(fact1)) ! in a way that prevents overflow ! double precision fact1,fact2,fact3,term1,term2 if(fact3.gt.fact1)then if(fact3.gt.-10.0d0)then if((dexp(-fact3)+dexp(fact1-fact3)).lt.1.0d-20)then term1=1.0d+30 else term1=1.0d0/(dexp(-fact3)+dexp(fact1-fact3)) endif else term1=dexp(fact3)/(1.0d0+dexp(fact1)) endif else if(fact1.gt.-10.0d0)then term1=dexp(fact3-fact1)/(dexp(-fact1)+1.0d0) else term1=dexp(fact3)/(1.0d0+dexp(fact1)) endif endif if((fact2+fact3).gt.fact1)then if((fact2+fact3).gt.-10.0d0)then if((dexp(-fact2-fact3)+dexp(fact1-fact2-fact3)) & .lt.1.0d-20)then term2=1.0d+30 else term2=1.0d0/ & (dexp(-fact2-fact3)+dexp(fact1-fact2-fact3)) endif else term2=dexp(fact2+fact3)/(1.0d0+dexp(fact1)) endif else if(fact1.gt.-10.0d0)then term2=dexp(fact2+fact3-fact1)/(dexp(-fact1)+1.0d0) else term2=dexp(fact2+fact3)/(1.0d0+dexp(fact1)) endif endif threeexp=term1+term2 return end !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine co2compens(templeaf,gammas25,ha_gammas, & gascon,gammas) implicit none ! purpose: calculates CO2 compensation point. !------- Methods ------------------- ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 double precision templeaf,gascon, & gammas,gammas25,ha_gammas !------- inputs -------------------- ! templeaf: leaf temperature [K] ! gascon: universal gas constant [JK-1 mol-1]. ! gammas25: the compensation point at 25oC [Pa] ! ha_gammas: parameter in the temperature response function, typically around ! 37.83 kJmol-1 !------- outputs ------------------- ! gammas: CO2 compensation point (Pa) gammas=gammas25*dexp((ha_gammas/gascon)* & (1000.0d0/298.15d0-1000.0d0/templeaf)) return end ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine resp_mitocho(templeaf,rdlight25, & ha_darkresp,gascon,rdlight) implicit none ! purpose: calculates mitochondrial respiration rate and its ! derivative wrt. leaf temperature. !------- Methods ------------------- ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 double precision templeaf,rdlight25,gascon, & rdlight,ha_darkresp,term !------- inputs -------------------- ! templeaf: leaf temperature [K] ! rdlight25: Mitochondrial respiration rate in the light at 25oC ! of leaf temperature [umol m-2 s-1] ! gascon: universal gas constant [JK-1 mol-1]. ! ha_darkresp: Mitochondrial respiation parameter deltaha = 46.39d0 unit in ! Jmol-1/1000 !------- outputs ------------------- ! rdlight: Mitochondrial respiration rate in the light [umol m-2 s-1] if(ha_darkresp.lt.0.0d0)then !The calling program uses the Q10 form. ha_darkresp is really -Q10 term=-ha_darkresp rdlight=rdlight25* & (term**((templeaf-298.15d0)/10.0d0)) else rdlight=rdlight25*dexp((ha_darkresp/gascon)* & (1000.0d0/298.15d0-1000.0d0/templeaf)) endif return end ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine MichaelisCO2(templeaf,fkc25,ha_kc,gascon, &fkc) implicit none ! purpose: calculates Michaelis constant for CO2 at a given leaf temperature. !------- Methods ------------------- ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 double precision templeaf,fkc25,ha_kc, & gascon,fkc !------- inputs -------------------- ! templeaf: leaf temperature [K] ! gascon: universal gas constant [JK-1 mol-1]. ! fkc25: the Michaelis constant for CO2 at 25oC ! ha_kc: a coefficient (~79.43 KJmol-1) !------- outputs ------------------- ! fkc: Michaelis constant for CO2 (Pa) fkc=fkc25*dexp((ha_kc/gascon)* & (1000.0d0/298.15d0-1000.0d0/templeaf)) return end ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine MichaelisO2(templeaf,fko25,ha_ko,gascon, &fko) implicit none ! purpose: calculates Michaelis constant for O2 at a given leaf temperature. !------- Methods ------------------- ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 double precision templeaf,fko25,ha_ko, & gascon,fko !------- inputs -------------------- ! templeaf: leaf temperature [K] ! gascon: universal gas constant [JK-1 mol-1]. ! fko25: the Michaelis constant for O2 at 25oC ! ha_ko: a coefficient in the temp function (~36.38 KJmol-1) !------- outputs ------------------- ! fko: Michaelis constant for CO2 (Pa) fko=fko25*dexp((ha_ko/gascon)* & (1000.0d0/298.15d0-1000.0d0/templeaf)) return end ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! subroutine tpuontemp(templeaf,gascon,tpu25, & dha,dhd,ds,tpu) implicit none ! !(in) templeaf: leaf temperature [K] !(in) gascon: universal gas constant [JK-1 mol-1]. !(in) tpu25: tpu at 25 oC [umol co2 m-2 s-1] !(in) dha ~ 53.1d0 (53100/1000) !(in) dhd ~ 201.8d0 (201800/1000) !(in) ds ~ 0.65d0 (650/1000) !(out) tpu: tpu at the given leaf temperature double precision templeaf,gascon,tpu25, & dha,dhd,ds,tpu ! !-------------- Locals -------------------- double precision t25k,univR,fact1,fact2,fact3,threeexp t25k=273.15d0+25.0d0 univR=gascon/1000.0d0 fact1=(ds*templeaf-dhd)/(univR*templeaf) fact2=(ds*t25k-dhd)/(univR*t25k) fact3=dha/(univR*t25k)-dha/(univR*templeaf) tpu=tpu25*threeexp(fact1,fact2,fact3) return end !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine vcmaxontemp(templeaf,vcmax25,gascon, & ha_vcmax,hd_vcmax,sv_vcmax,vcmax) implicit none ! purpose: calculates maximum RuBP saturated rate of carboxylation. !------- Methods ------------------- ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 double precision templeaf,vcmax25,gascon,ha_vcmax,hd_vcmax, & sv_vcmax,vcmax !------- inputs -------------------- ! templeaf: leaf temperature [K] ! vcmax25: maximum RuBP saturated rate of carboxylation at 25oC ! of leaf temperature [umol m-2 s-1] ! gascon: universal gas constant [JK-1 mol-1]. ! ha_vcmax: ! hd_vcmax: ! sv_vcmax: !------- outputs ------------------- ! vcmax: maximum RuBP saturated rate of carboxylation [umol m-2 s-1] !------- Local variables ----------- double precision threeexp double precision T0,univR,fact1,fact2,fact3 univR=gascon/1000.0d0 T0=25.0d0+273.15d0 fact1=(sv_vcmax*templeaf-hd_vcmax)/(univR*templeaf) fact2=(sv_vcmax*T0-hd_vcmax)/(univR*T0) fact3=ha_vcmax*(1.0d0-T0/templeaf)/(univR*T0) ! threeexp=(dexp(fact3)+dexp(fact2+fact3))/(1.0d0+dexp(fact1)) vcmax=vcmax25*threeexp(fact1,fact2,fact3) return end ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine gmesoontemp(templeaf,gmeso25,gascon, & ha_gmeso,hd_gmeso,sv_gmeso,gmeso) implicit none ! purpose: calculates mesophyll CO2 conductance at a given temperature !------- Methods ------------------- ! Based on functions in Bernacchi, et al. (2002), ! Plant Physiology, 130, 1992-1998, or Scafaro et al. (2011), PCE 34: 1999-2008 double precision templeaf,gmeso25,gascon,ha_gmeso,hd_gmeso, & sv_gmeso,gmeso !------- inputs -------------------- ! templeaf: leaf temperature [K] ! ha_gmeso: The activation energy (Ha) in gmeso temperature response function [kJmol-1] ! hd_gmeso: The deactivation energy (Hd) in gmeso temperature response function [kJmol-1] ! sv_gmeso: The entropy term (Sv) in gmeso temperature response function [kJmol-1K-1] ! gascon: The universal gas constant (JK-1mol-1) !------- outputs ------------------- ! gmeso: Mesophyll CO2 transfer conductance [umol m-2 s-1 Pa-1] !------- Local variables ----------- double precision T0,univR,fact1,fact2, & fact3,threeexp univR=gascon/1000.0d0 T0=25.0d0+273.15d0 ! if(dabs(sv_gmeso+9999.0d0).lt.1.0d-5.or. &dabs(hd_gmeso+9999.0d0).lt.1.0d-5)then gmeso=gmeso25*dexp(ha_gmeso*(1.0d0-T0/templeaf)/(univR*T0)) return endif fact1=(sv_gmeso*templeaf-hd_gmeso)/(univR*templeaf) fact2=(sv_gmeso*T0-hd_gmeso)/(univR*T0) fact3=ha_gmeso*(1.0d0-T0/templeaf)/(univR*T0) ! threeexp=(dexp(fact3)+dexp(fact2+fact3))/(1.0d0+dexp(fact1)) gmeso=gmeso25*threeexp(fact1,fact2,fact3) return end ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine jmaxontemp(templeaf,fjmax25,gascon, & ha_jmax,hd_jmax,sv_jmax,fjmax) implicit none ! purpose: calculates maximum rate of electron transport ! and its derivative wrt. leaf temperature. !------- Methods ------------------- ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 ! Sept. 16, 2005. Bernacchi et al. functions are replaced by ! June et al. function, June T., J. R. Evans, G. D. Farquhar (2004) ! A simple new equation for the reversible temperature dependence of ! photosynthetic electron transport: a study on soybean leaf, ! Functional Plant Biology 31, 275-283. double precision templeaf,fjmax25,gascon, & ha_jmax,hd_jmax,sv_jmax,fjmax !------- inputs -------------------- ! templeaf: leaf temperature [K] ! fjmax25: maximum rate of electron transport at 25oC ! of leaf temperature [umol m-2 s-1] ! gascon: universal gas constant [JK-1 mol-1]. ! ha_jmax ~ temperature response function parameter ! hd_jmax ~ temperature response function parameter ! sv_jmax ~ temperature response function parameter !------- outputs ------------------- ! fjmax: maximum rate of electron transport [umol m-2 s-1] !------- Local variables ----------- double precision T0,omega,threeexp double precision univR,fact1,fact2,fact3 univR=gascon/1000.0d0 T0=25.0d0+273.15d0 fact1=(sv_jmax*templeaf-hd_jmax)/(univR*templeaf) fact2=(sv_jmax*T0-hd_jmax)/(univR*T0) fact3=ha_jmax*(1.0d0-T0/templeaf)/(univR*T0) ! threeexp=(dexp(fact3)+dexp(fact2+fact3))/(1.0d0+dexp(fact1)) fjmax=fjmax25*threeexp(fact1,fact2,fact3) ! June et al. function ! ha_jmax = 30.0d0+273.15d0 ! hd_jmax = 11.6d0 ! sv_jmax = 0.18d0 ! T0=ha_jmax-273.15d0 ! omega=hd_jmax+sv_jmax*T0 ! fjmax=fjmax25*dexp( ! & (templeaf-273.15d0+25.0d0-2.0d0*T0)* ! & (25.0d0-templeaf+273.15d0)/(omega*omega)) return end ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine jontemp(aPPFDlf,templeaf,fjelect,fjmax25,ha_jmax, &hd_jmax,sv_jmax,gascon,phifactor,thetafactor,betaPSII) implicit none ! purpose: calculates rate of electron transport ! and its derivative wrt. leaf temperature. !------- Methods ------------------- ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 !3/29/2014. Bernacchi method cannot fit A/PAR curves well. phifactor and thetafactor are !used to modify his method when A/PAR curves are present. double precision aPPFDlf,templeaf,fjelect,fjmax25,gascon, &ha_jmax,hd_jmax,sv_jmax,phifactor,thetafactor,betaPSII !------- inputs -------------------- ! aPPFDlf: absorbed photosynthetic photon flux density by leaf (umol m-2 s-1) ! templeaf: leaf temperature [K] ! fjmax25: maximum rate of electron transport at 25oC ! of leaf temperature [umol m-2 s-1] ! gascon: universal gas constant [JK-1 mol-1]. ! ha_jmax ~ 50300.0d0/1000.0d0 ! hd_jmax ~ 152044.0d0/1000.0d0 ! sv_jmax ~ 495.0d0/1000.0d0 ! phifactor: modifier for phiPSIImax ! thetafactor: modifer for thetaPSII ! betaPSII: fraction of absorbed light that reaches photosystem II (~0.5) !------- outputs ------------------- ! fjelect: rate of electron transport [umol m-2 s-1] !------- Local variables ----------- double precision Q2,phiPSIImax, &thetaPSII,zeroT,b2_4ac,b2_4ac_sqrt,fjmax parameter (zeroT=273.15d0) ! ! beta: fraction of absorbed qunata reaching pSII [--] ! Q2: maximum incident quanta that can be utilized in ! electron transport (umol m-2 s-1) ! phiPSIImax: Maximum dark-adapted quantum yield of PSII [--] ! thetaPSII: convexity term for electron transport rates [--] ! zeroT: zero temperature [K] ! b2-4ac: temporay variable ! b2_4ac_sqrt: sqrt of b2_4ac ! fjmax: maximum rate of electron transport, computed in the ! subroutine of jmaxontime [umol m-2 s-1] call thetaphipsii(templeaf,phiPSIImax,thetaPSII) phiPSIImax=phiPSIImax*phifactor thetaPSII=thetaPSII*thetafactor Q2=aPPFDlf*phiPSIImax*betaPSII call jmaxontemp(templeaf,fjmax25,gascon, & ha_jmax,hd_jmax,sv_jmax,fjmax) b2_4ac=(Q2+fjmax)*(Q2+fjmax)-4.0d0*thetaPSII*Q2*fjmax if(b2_4ac.lt.0.0d0)then b2_4ac=0.0d0 endif b2_4ac_sqrt=dsqrt(b2_4ac) fjelect=(Q2+fjmax-b2_4ac_sqrt)/(2.0d0*thetaPSII) return end !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine jmaxfromj(aPPFDlf,templeaf,fjelect,fjmax, &phifactor,thetafactor,betaPSII) implicit none ! purpose: calculates jmax from the rate of electron transport, ! leaf temperature and absorbed PAR !------- Methods ------------------- ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 !3/29/2014. Bernacchi method cannot fit A/PAR curves well. phifactor and thetafactor are !used to modify his method when A/PAR curves are present. ! double precision aPPFDlf,templeaf,fjelect,fjmax, &phifactor,thetafactor,betaPSII !------- inputs -------------------- ! aPPFDlf: absorbed photosynthetic photon flux density by leaf (umol m-2 s-1) ! templeaf: leaf temperature [K] ! fjelect: rate of electron transport [umol m-2 s-1] ! phifactor: modifier for phiPSIImax ! thetafactor: modifer for thetaPSII ! betaPSII: the fracction of absorbed light that reaches photosystem II (~0.5) !------- outputs ------------------- ! fjmax: maximum rate of electron transport ! of leaf temperature [umol m-2 s-1] !------- Local variables ----------- double precision Q2,phiPSIImax,thetaPSII,zeroT parameter (zeroT=273.15d0) ! ! beta: fraction of absorbed qunata reaching pSII [--] ! Q2: maximum incident quanta that can be utilized in ! electron transport (umol m-2 s-1) ! phiPSIImax: Maximum dark-adapted quantum yield of PSII [--] ! thetaPSII: convexity term for electron transport rates [--] ! zeroT: zero temperature [K] ! call thetaphipsii(templeaf,phiPSIImax,thetaPSII) phiPSIImax=phiPSIImax*phifactor thetaPSII=thetaPSII*thetafactor Q2=aPPFDlf*phiPSIImax*betaPSII fjmax=(thetaPSII*fjelect-Q2)*fjelect/(fjelect-Q2) return end subroutine thetaphipsii(templeaf0,phiPSIImax,thetaPSII) implicit none ! Based on functions in Bernacchi, Pimentel and Long (2003), ! Plant, Cell and Environment, 26, 1419-1430 ! templeaf0: leaf temperature [K] !------- outputs ------------------- ! phiPSIImax: Maximum dark-adapted quantum yield of PSII [--] ! thetaPSII: convexity term for electron transport rates [--] double precision phiPSIImax,thetaPSII,templeaf0 double precision zeroT,templeaf parameter (zeroT=273.15d0) ! zeroT: zero temperature [K] templeaf=dmax1(zeroT+5.0d0,templeaf0) templeaf=dmin1(zeroT+45.0d0,templeaf) phiPSIImax=0.352d0+0.022d0*(templeaf-zeroT)-3.4d-4* & (templeaf-zeroT)*(templeaf-zeroT) thetaPSII=0.76d0+0.018d0*(templeaf-zeroT)-3.7d-4* & (templeaf-zeroT)*(templeaf-zeroT) return end !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! double precision function atpcon(wc,wj,vcmax) implicit none ! ! Calculate the concentration of ATP provided by photophosphorylation ! based on Buckley et al. (2003), Plant, Cell and Environment, 26, 1767-1785. ! ! ------ Inputs --------------------------------------------------- ! wc: Rubisco limited rate of carboxylation [umol m-2 s-1] ! wj: RuBP limited rate of carboxylation [umol m-2 s-1] ! vcmax: maximum RuBP saturated rate of carboxylation [umol m-2 s-1] double precision wc,wj,vcmax ! ! ------ Outputs --------------------------------------------------- ! atpcon: the concentration of ATP (mmol ATP m-2) ! ------ Local variables -------------------------------------------- ! at: at in Buckley et al (2003), total concentration of adenylates ! [mmol AxP m-2] ! p: p in Buckley et al (2003), concentration of photophosphorylation ! sites [mmol sites m-2] ! Communication with Buckley reveals error in Table 1 in Buckley paper ! at and p should be divided by 100. ! vr: Vr in Buckley et al (2003), carboxylation rate limited by potential ! RuBP pool size only [umol m-2 s-1] ! tau0: basal ATP level provided by other processes [mmol ATP m-2] ! tauc: atpcon when wc limiting mmol m-2 ! tauj: atpcon when wj limiting mmol m-2 double precision at,p,vr,tau0,tauc,tauj parameter(tau0=1.6d0) if(wj.le.0.0d0)then atpcon=tau0 else at=12.6d0*vcmax/100.0d0 p=2.5d0*vcmax/100.0d0 vr=2.27d0*vcmax tauc=at-p*wc/wj tauj=(at-p)*(vr/vcmax-1.0d0)/(wc*vr/(wj*vcmax)-1.0d0) if(wc.lt.wj)then atpcon=tau0+tauc else atpcon=tau0+tauj endif endif return end ! !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine QSat (T, p, es, esdT, qs, & qsdT) implicit none !----------------------------------------------------------------------- ! ! CLMCLMCLMCLMCLMCLMCLMCLMCLMCL A community developed and sponsored, freely ! L M available land surface process model. ! M --COMMUNITY LAND MODEL-- C ! C L ! LMCLMCLMCLMCLMCLMCLMCLMCLMCLM ! !----------------------------------------------------------------------- ! Purpose: ! Computes saturation mixing ratio and the change in saturation ! mixing ratio with respect to temperature. ! Method: ! Reference: Polynomial approximations from: ! Piotr J. Flatau, et al.,1992: Polynomial fits to saturation ! vapor pressure. Journal of Applied Meteorology, 31, 1507-1513. ! ! Author: ! 15 September 1999: Yongjiu Dai; Initial code ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision ! April 2002: Vertenstein/Oleson/Levis; Final form ! !----------------------------------------------------------------------- ! $Id: QSat.F90,v 1.1.8.4 2002/04/10 19:25:07 mvertens Exp $ !----------------------------------------------------------------------- ! use precision ! use shr_const_mod, ONLY: SHR_CONST_TKFRZ ! implicit none !----Arguments---------------------------------------------------------- DOUBLE PRECISION T, p DOUBLE PRECISION es, esdT, qs, qsdT DOUBLE PRECISION SHR_CONST_TKFRZ PARAMETER (SHR_CONST_TKFRZ=273.16d0) ! real(r8), intent(in) :: T ! temperature (K) ! real(r8), intent(in) :: p ! surface atmospheric pressure (pa) ! real(r8), intent(out) :: es ! vapor pressure (pa) ! real(r8), intent(out) :: esdT ! d(es)/d(T) ! real(r8), intent(out) :: qs ! humidity (kg vapor/kg moist air) ! real(r8), intent(out) :: qsdT ! d(qs)/d(T) !----Local Variables---------------------------------------------------- DOUBLE PRECISION T_limit DOUBLE PRECISION td,vp,vp1,vp2 DOUBLE PRECISION a0,a1,a2,a3,a4,a5,a6,a7,a8 DOUBLE PRECISION b0,b1,b2,b3,b4,b5,b6,b7,b8 DOUBLE PRECISION c0,c1,c2,c3,c4,c5,c6,c7,c8 DOUBLE PRECISION d0,d1,d2,d3,d4,d5,d6,d7,d8 ! real(r8) T_limit ! limitation on valid temperatures [C] ! real(r8) td,vp,vp1,vp2 ! real(r8) a0,a1,a2,a3,a4,a5,a6,a7,a8 ! real(r8) b0,b1,b2,b3,b4,b5,b6,b7,b8 ! real(r8) c0,c1,c2,c3,c4,c5,c6,c7,c8 ! real(r8) d0,d1,d2,d3,d4,d5,d6,d7,d8 ! ! For water vapor (temperature range 0C-100C) ! data a0/6.11213476d0/,a1/0.444007856d0/,a2/0.143064234d-01/, & a3/0.264461437d-03/,a4/0.305903558d-05/,a5/0.196237241d-07/, & a6/0.892344772d-10/,a7/-0.373208410d-12/,a8/0.209339997d-15/ ! ! For derivative:water vapor ! data b0/0.444017302d0/,b1/0.286064092d-01/,b2/0.794683137d-03/, & b3/0.121211669d-04/,b4/0.103354611d-06/,b5/0.404125005d-09/, & b6/-0.788037859d-12/,b7/-0.114596802d-13/,b8/0.381294516d-16/ ! ! For ice (temperature range -75C-0C) ! data c0/6.11123516d0/,c1/0.503109514d0/,c2/0.188369801d-01/, & c3/0.420547422d-03/,c4/0.614396778d-05/,c5/0.602780717d-07/, & c6/0.387940929d-09/,c7/0.149436277d-11/,c8/0.262655803d-14/ ! ! For derivative:ice ! data d0/0.503277922d0/,d1/0.377289173d-01/,d2/0.126801703d-02/, & d3/0.249468427d-04/,d4/0.313703411d-06/,d5/0.257180651d-08/, & d6/0.133268878d-10/,d7/0.394116744d-13/,d8/0.498070196d-16/ !----End Variable List-------------------------------------------------- T_limit = T - SHR_CONST_TKFRZ if (T_limit .GT. 100.0d0) T_limit=100.0d0 if (T_limit .LT. -75.0d0) T_limit=-75.0d0 td = T_limit if(td .GE. 0.0d0) then es=a0+td*(a1+td*(a2+td*(a3+td*(a4+ & td*(a5 + td*(a6 + td*(a7 + td*a8))))))) esdT=b0+td*(b1+td*(b2+td*(b3+td*(b4+ & td*(b5+td*(b6+td*(b7+td*b8))))))) else es=c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4 & + td*(c5 + td*(c6 + td*(c7 + td*c8))))))) esdT=d0+td*(d1+td*(d2 + td*(d3 + td*(d4 & + td*(d5 + td*(d6 + td*(d7 + td*d8))))))) endif es = es * 100.0d0 ! pa esdT = esdT * 100.0d0 ! pa/K vp = 1.00d0 / (p - 0.378d0*es) vp1 = 0.622d0 * vp vp2 = vp1 * vp qs = es * vp1 ! kg/kg qsdT = esdT * vp2 * p ! 1 / K end double precision function esat(T, p) implicit none !----------------------------------------------------------------------- ! ! CLMCLMCLMCLMCLMCLMCLMCLMCLMCL A community developed and sponsored, freely ! L M available land surface process model. ! M --COMMUNITY LAND MODEL-- C ! C L ! LMCLMCLMCLMCLMCLMCLMCLMCLMCLM ! !----------------------------------------------------------------------- ! Purpose: ! Computes saturation mixing ratio and the change in saturation ! mixing ratio with respect to temperature. ! Method: ! Reference: Polynomial approximations from: ! Piotr J. Flatau, et al.,1992: Polynomial fits to saturation ! vapor pressure. Journal of Applied Meteorology, 31, 1507-1513. ! ! Author: ! 15 September 1999: Yongjiu Dai; Initial code ! 15 December 1999: Paul Houser and Jon Radakovich; F90 Revision ! April 2002: Vertenstein/Oleson/Levis; Final form ! !----------------------------------------------------------------------- ! $Id: QSat.F90,v 1.1.8.4 2002/04/10 19:25:07 mvertens Exp $ !----------------------------------------------------------------------- ! use precision ! use shr_const_mod, ONLY: SHR_CONST_TKFRZ ! implicit none !----Arguments---------------------------------------------------------- DOUBLE PRECISION T, p DOUBLE PRECISION es DOUBLE PRECISION SHR_CONST_TKFRZ PARAMETER (SHR_CONST_TKFRZ=273.15d0) ! real(r8), intent(in) :: T ! temperature (K) ! real(r8), intent(in) :: p ! surface atmospheric pressure (pa) ! real(r8), intent(out) :: es ! vapor pressure (pa) !----Local Variables---------------------------------------------------- DOUBLE PRECISION T_limit DOUBLE PRECISION td DOUBLE PRECISION a0,a1,a2,a3,a4,a5,a6,a7,a8 DOUBLE PRECISION b0,b1,b2,b3,b4,b5,b6,b7,b8 DOUBLE PRECISION c0,c1,c2,c3,c4,c5,c6,c7,c8 ! real(r8) T_limit ! limitation on valid temperatures [C] ! real(r8) td ! real(r8) a0,a1,a2,a3,a4,a5,a6,a7,a8 ! real(r8) c0,c1,c2,c3,c4,c5,c6,c7,c8 ! For water vapor (temperature range 0C-100C) ! data a0/6.11213476d0/,a1/0.444007856d0/,a2/0.143064234d-01/, & a3/0.264461437d-03/,a4/0.305903558d-05/,a5/0.196237241d-07/, & a6/0.892344772d-10/,a7/-0.373208410d-12/,a8/0.209339997d-15/ ! ! For ice (temperature range -75C-0C) ! data c0/6.11123516d0/,c1/0.503109514d0/,c2/0.188369801d-01/, & c3/0.420547422d-03/,c4/0.614396778d-05/,c5/0.602780717d-07/, & c6/0.387940929d-09/,c7/0.149436277d-11/,c8/0.262655803d-14/ !----End Variable List-------------------------------------------------- T_limit = T - SHR_CONST_TKFRZ if (T_limit .GT. 100.0d0) T_limit=100.0d0 if (T_limit .LT. -75.0d0) T_limit=-75.0d0 td = T_limit if(td .GE. 0.0d0) then es=a0+td*(a1+td*(a2+td*(a3+td*(a4+ & td*(a5 + td*(a6 + td*(a7 + td*a8))))))) else es=c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4 & + td*(c5 + td*(c6 + td*(c7 + td*c8))))))) endif es = es * 100.0d0 ! pa esat=es !The following is what used in Li-Cor 6400 so we use it for A/Ci analysis esat=613.65d0*dexp(17.502d0*(T-273.15d0)/ & (240.97d0+T-273.15d0)) return end ! !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine metvarconverter(airtemp,pres_moistair,vpd, & pres_h2ovapor,rh,gascon,vapormixingratio,spechumid, & cp_moistair,dens_moistair,dewpoint,wmoistair) implicit none ! ! =============Inputs=================================== ! airtemp: air temperature [K] ! pres_moistair: air pressure [Pa] ! rh or vpd or pres_h2ovapor: one of the three must be inputs ! rh: relative humidity [%] ! vpd: vapor pressure deficit [Pa] ! pres_h2ovapor: water vapor partial pressure [Pa] ! gascon: universal gas constant [J K-1 mol-1] ! double precision airtemp,pres_moistair,vpd, & pres_h2ovapor,rh,gascon !==============Outputs=================================== ! rh, vpd, and pres_h2ovapor: if not provided as inputs ! vapormixingratio: water vapor mixing ratio [g vapor / g dry air] ! spechumid: specific humidity [g vapor / g moist air] ! cp_moistair: specific heat of the moist air [J g-1 K-1] ! dens_moistair: density of the moist air [gm-3] ! dewpoint: dewpoint temperature [K] ! wmoistair: molecular weight of moist air [g mol-1] ! double precision vapormixingratio,spechumid, & cp_moistair,dens_moistair,dewpoint,wmoistair !==============Locals=================================== ! dryairmw: molecular weight of the dry air [g mol-1] ! h2omw: molecular weight of water [g mol-1] ! cp_dryair: specific heat of dry air [J g-1 K-1] ! cp_watervapor: specific heat of water vapor [J g-1 K-1] ! es: saturation vapor pressure [Pa] ! esdT: not used ! qs: specific humidity at saturation [g vapor / g moist air] ! qsdT: not used ! pres_dryair: partial pressure of the dry air [Pa] ! dens_dryair: density of the dry air [gm-3] ! dens_vapor: density of water vapor [gm-3] ! ws: mixing ratio at saturation [ density vapor / density dry air] ! abszero: absolute zero temperature [K] ! latent: latent heat of vaporization [J g-1] double precision dryairmw,h2omw,cp_dryair,cp_watervapor, & es,esdT,qs,qsdT,pres_dryair,dens_dryair,dens_vapor,esat, & ws,abszero,latent parameter(dryairmw=28.9644d0,h2omw=18.016d0) parameter(cp_dryair=1.00467d0,cp_watervapor=1.84d0) parameter(abszero=273.15d0) ! if((rh.le.0.0d0.or.rh.gt.100.02d0).and.(vpd.lt.0.0d0.or. & vpd.gt.9000.0d0).and.(pres_h2ovapor.le.0.0d0.or. & pres_h2ovapor.gt.9000.0d0))then write(*,*)'At least of the one variables VPD, RH, and & water vapor partial pressure must be valid. & check your data. Program stops' stop endif call QSat(airtemp,pres_moistair,es,esdT,qs,qsdT) ws=qs/(1.0d0-qs) if((rh.le.0.0d0.or.rh.gt.100.0d0).or.(vpd.le.0.0d0.or. & vpd.gt.9000.0d0).or.(pres_h2ovapor.le.0.0d0.or. & pres_h2ovapor.gt.9000.0d0))then if(pres_h2ovapor.gt.0.0d0.and.pres_h2ovapor.lt.9000.0d0)then pres_dryair=pres_moistair-pres_h2ovapor dens_dryair=pres_dryair*dryairmw/(gascon*airtemp) dens_vapor=pres_h2ovapor*h2omw/(gascon*airtemp) vapormixingratio=dens_vapor/dens_dryair vpd=es-pres_h2ovapor vpd=dmin1(es,vpd) vpd=dmax1(0.0d0,vpd) rh=vapormixingratio*100.0d0/ws rh=dmin1(100.0d0,rh) rh=dmax1(0.0d0,rh) else if(rh.gt.0.0d0.and.rh.le.100.01d0)then vapormixingratio=ws*dmin1(rh,100.0d0)/100.0d0 pres_h2ovapor=vapormixingratio*pres_moistair/( & vapormixingratio+h2omw/dryairmw) vpd=es-pres_h2ovapor else if(vpd.ge.0.0d0.and.vpd.lt.9000.0d0)then pres_h2ovapor=es-vpd pres_h2ovapor=dmin1(es,pres_h2ovapor) pres_h2ovapor=dmax1(0.0d0,pres_h2ovapor) pres_dryair=pres_moistair-pres_h2ovapor dens_dryair=pres_dryair*dryairmw/(gascon*airtemp) dens_vapor=pres_h2ovapor*h2omw/(gascon*airtemp) vapormixingratio=dens_vapor/dens_dryair rh=vapormixingratio*100.0d0/ws endif endif endif endif pres_dryair=pres_moistair-pres_h2ovapor dens_dryair=pres_dryair*dryairmw/(gascon*airtemp) dens_vapor=pres_h2ovapor*h2omw/(gascon*airtemp) dens_moistair=dens_dryair+dens_vapor vapormixingratio=dens_vapor/dens_dryair spechumid=dens_vapor/dens_moistair cp_moistair=(cp_dryair*dens_dryair+ & cp_watervapor*dens_vapor)/dens_moistair latent=1.91846d+3*(airtemp/(airtemp-33.91d0))* & (airtemp/(airtemp-33.91d0)) call QSat(abszero,pres_moistair,es,esdT,qs,qsdT) dewpoint=1.0d0/abszero-(gascon/(h2omw*latent))* & dlog(pres_h2ovapor/es) dewpoint=1.0d0/dewpoint wmoistair=dens_moistair*gascon*airtemp/pres_moistair return end ! c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& double precision function latent(temp) c c This function calculates latent heat from temperature. The formulae c is taken from Henderson-sellers, 1984, A new formulae for latent heat c of vaporation of water as a function of temperature, Quart. J. R. Met. c Soc. 1984, 110 pp. 1186-1190. c c temp: temperature in K. c latent: latent heat of water in J g-1. implicit double precision (a-h,l,o-z) parameter(abszero=273.15d0,fusionheat=333.7d0) ! latent=1.91846d+3*(temp/(temp-33.91d0))* & (temp/(temp-33.91d0)) ! ! assume freezing or sublimation occurs at 273.15d0 K ! latent heat of fusion is 333.7 J/g ! if(temp.le.abszero)then latent=latent+fusionheat endif return end c !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$