1042 lines
40 KiB
FortranFixed
1042 lines
40 KiB
FortranFixed
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
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
|
|
|
|
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|