1811 lines
76 KiB
FortranFixed
1811 lines
76 KiB
FortranFixed
subroutine SetUpLeafGasFit(icurveno_usr,curvename,ntotsamples0,
|
|
&CurveTypeID,anet_obs0,pco2i0,templeaf0,PARi0,pres_air0,po2i0,
|
|
&chlflphips20,pco2ambient0,trmmol0,gswmeas0,vpdl0,tempair0,
|
|
&eambient0,fo_pam0,fm_pam0,fs_pam0,pam_measlight0,stargamma25_usr,
|
|
&fkc25_usr,fko25_usr,rdlight25_usr,alpha25_usr,resistwp25_usr,
|
|
&resistch25_usr,isitmassbased,indexunit,
|
|
&siteID,Latitude,Longitude,Elevation,yearsampled,sampledoy,
|
|
&GrowingSeasonStart,GrowingSeasonEnd,standage,CanopyHeight,
|
|
&LeafAreaIndex,species,avetimeresolution,avetimesampled,
|
|
&SampleHeight,Needleage,specificLAI,nitrogencontent,carboncontent,
|
|
&phoscontent,woodporosity,sapwooddensity,leafratio)
|
|
implicit none
|
|
include '../testarea/LeafGasParams.h'
|
|
include '../testarea/LeafGasHybridFit.h'
|
|
!--------------------Inputs--------------------------------------------------------------
|
|
!None of the inputs is changed by this subroutine
|
|
!icurveno_usr(int): the curve number
|
|
!curvename(char): the curve name
|
|
!ntotsamples0: the total number of data points
|
|
!CurveTypeID =1-3: Any measurements where limitation states are known:
|
|
! =1 limited by Rubisco
|
|
! =2 limited by RuBp regeneration
|
|
! =3 limited by TPU
|
|
! =11-25: ACi Curves. Each different CurveTypeID number represents a different A/Ci curve (i.e., different PAR levels).
|
|
! For example, five different PAR levels are used to measure five A/Ci curves with PAR = 200, 400, 600, 800,
|
|
! 1000.Use 11, 12, 13, 14, 15 to indentify points of each curve. Maxumum 15 A/Ci curves.
|
|
! The curves must be numbered consecutively.
|
|
! =31-45: ALight Curves. Each different CurveTypeID number represents a different A/Light curve (i.e., different ambient CO2 levels).
|
|
! For example, five different ambient levels are used to measure five A/Light curves with CO2a= 100, 200, 300, 400, 500.
|
|
! Use 31, 32, 33, 34, 35 to indentify points of each ALight curve. The curves must be ordered consecutively.
|
|
! =-9999: all other types of measurements.
|
|
!anet_obs0: Net photosynthetic rate (umol m-2 s-1)
|
|
!pco2i0: Intercellular CO2 concentration (Pa)
|
|
!templeaf0: leaf temperature (K)
|
|
!PARi0: The PAR level inside the chamber to which photosynthesis responds (umolm-2s-1)
|
|
!pres_air0: Atmospheric pressure (Pa)
|
|
!po2i0: Oxygen partial presssure (Pa)
|
|
!chlflphips20: Chlorophyll fluorescence (NA), that is, DeltaF/Fm, the fraction of
|
|
! absorbed PSII photons that are used in photochemistry
|
|
!pco2ambient0: Ambient CO2 partial pressure (Pa)
|
|
!trmmol0: Transpiration rate (mmolm-2s-1)
|
|
!gswmeas0: Stomatal conductance for water vapor (molm-2s-1)
|
|
!vpdl0: Water vapor pressure difference between the leaf and chamber air (Pa)
|
|
!tempair0: Air temperature inside the chamber (K)
|
|
!eambient0: Water vapor pressure inside the chamber (Pa).
|
|
!fo_pam0: fo (dark adapated) or fo' (actinic light turned off, far red light on to drain electrons from PSII) from pulse amplitude modulation (arbitrary unit).
|
|
!fm_pam0: fm (dark adapated with saturation pulse) or fm' (actinic light with saturation pulse) from pulse amplitude modulation (arbitrary unit).
|
|
!fs_pam0: steady state fluorescence from pulse amplitude modulation (arbitrary unit).
|
|
!pam_measlight0: the measuring light level (umolm-2s-1)
|
|
!stargamma25_usr: Chloroplastic CO2 compenstation point at 25oC provided by the user (Pa), set to -9999 if not available
|
|
!fkc25_usr: the Michaelis constant for CO2 at 25oC provided by the user (Pa), set to -9999 if not available
|
|
!fko25_usr: the Michaelis constant for O2 at 25oC provided by the user (Pa), set to -9999 if not available
|
|
!rdlight25_usr: Leaf dark respiration at 25oC provided by user (Pa), set to -9999 if not available
|
|
!alpha25_usr: The fraction of glycolate carbon not returned to the chloroplast at 25oC provided by user (NA), set to -9999 if not available
|
|
!resistwp25_usr: resistance to CO2 via cell walls and plasmalemma provided by user [umol-1msPa], set to -9999 if not available
|
|
!resistch25_usr: resistance to CO2 via chloroplastic envelope provided by user[umol-1msPa], set to -9999 if not available
|
|
!isitmassbased: = 0, area-based (typical)
|
|
! = 1. mass-based (atypical)
|
|
!paramunit: file unit number to write ouputs
|
|
!compareunit: file unit number to write ouputs
|
|
!stomwuecicaoutunit: file unit number to write ouputs
|
|
!wuecicacompunit: file unit number to write ouputs
|
|
!stomcompunit: file unit number to write ouputs
|
|
!fluorescenceunit: file unit number to write outputs for comparison from fluorescence fit
|
|
!fluoresparamunit: file unit number to write parameters from fluorescence fit
|
|
!General information,not used but recorded in the output files
|
|
! & siteID,Latitude,Longitude,Elevation,yearsampled,
|
|
! & sampledoy,GrowingSeasonStart,GrowingSeasonEnd,
|
|
! & standage,CanopyHeight,LeafAreaIndex,species,
|
|
! & avetimeresolution,avetimesampled,SampleHeight,
|
|
! & Needleage,specificLAI,nitrogencontent,carboncontent,
|
|
! & phoscontent,woodporosity,sapwooddensity,leafratio)
|
|
integer icurveno_usr,ntotsamples0,isitmassbased,indexunit(20)
|
|
character*100 curvename
|
|
character siteID*(*),species*(*),woodporosity*(*)
|
|
double precision CurveTypeID(ntotsamples0),
|
|
&anet_obs0(ntotsamples0),pco2i0(ntotsamples0),
|
|
&templeaf0(ntotsamples0),PARi0(ntotsamples0),
|
|
&pres_air0(ntotsamples0),po2i0(ntotsamples0),
|
|
&chlflphips20(ntotsamples0),pco2ambient0(ntotsamples0),
|
|
&trmmol0(ntotsamples0),gswmeas0(ntotsamples0),vpdl0(ntotsamples0),
|
|
&tempair0(ntotsamples0),eambient0(ntotsamples0),
|
|
&fo_pam0(ntotsamples0),fm_pam0(ntotsamples0),
|
|
&fs_pam0(ntotsamples0),pam_measlight0(ntotsamples0),
|
|
&stargamma25_usr,fkc25_usr,fko25_usr,
|
|
&rdlight25_usr,alpha25_usr,resistwp25_usr,resistch25_usr,
|
|
!General information,not used but recorded in the output files
|
|
&Latitude,Longitude,Elevation,yearsampled,sampledoy,
|
|
&GrowingSeasonStart,GrowingSeasonEnd,standage,CanopyHeight,
|
|
&LeafAreaIndex,avetimeresolution,avetimesampled,SampleHeight,
|
|
&Needleage,specificLAI,nitrogencontent,carboncontent,
|
|
&phoscontent,sapwooddensity,leafratio
|
|
!------------------------------------------------------------------------------------------
|
|
character*30 modeltype,fourchars(20)
|
|
dimension modeltype(0:10)
|
|
|
|
integer i,j,k,m,n,idorwp0,idorch0,irchoption1,irchoption2,i2ndary,
|
|
&numrubis,numrubp,numtpu,INFO,iderivative,idoalpha0,
|
|
&ioriorder(3*ntotsamples0),ibelong(3*ntotsamples0),ACiID(15),
|
|
&ALightID(15),paramunit,compareunit,stomwuecicaoutunit,
|
|
&stomcompunit,wuecicacompunit,fluorescenceunit,
|
|
&fluoresparamunit,aciempfitunit,alightempfitunit,idotempcoeff,
|
|
&idomeso,idohavjt
|
|
|
|
double precision vcmax25_ini,fjmax25_ini,tpu25_ini,rdlight25_ini,
|
|
&stargamma25_ini,fkc25_ini,fko25_ini,alpha25_ini,resistwp25_ini,
|
|
&resistch25_ini,resiststomco20(ntotsamples0),term,term1,term2,
|
|
&aPPFDlf0(ntotsamples0),weitx(ntotsamples0),xmin(ntotsamples0),
|
|
&xmax(ntotsamples0),weity(ntotsamples0),beta(20),starco2i(15),
|
|
&der_starco2i(15),Amax_ACi(15),ACiinter(15),der_ACiinter(15),
|
|
&der_ACiend(15),starPAR(15),der_starPAR(15),Asat_ALight(15),
|
|
&ALightinter(15),der_ALightinter(15),der_ALightend(15),
|
|
&PhiPSIIzero_ACi(15),der_PhiPSIIzero_ACi(15),PhiPSIImax_ACi(15),
|
|
&PhiPSIIinter_ACi(15),der_PhiPSIIinter_ACi(15),
|
|
&der_PhiPSIIend_ACi(15),ExcessLightFactor(15),
|
|
&der_PhiPSII1000_ALight(15),PhiPSIIinter_ALight(15),
|
|
&der_PhiPSIIinter_ALight(15),amaxave,recycleratio(6,ntotsamples0),
|
|
&stargamma25fit(6),ACiavetempleaf(15),ACiaveaPPFDlf(15),
|
|
&ACiavepo2i(15),ALightavetempleaf(15),ALightaveCO2ambient(15),
|
|
&ALightavepo2i(15),co2c_Pa(4,ntotsamples0),co2imany(500),
|
|
&critdelPAR,critdelCi_Pa,rdlight,atp,resistwp,resistch,stargamma,
|
|
&ccc,ccj,cct,ac,aj,at,phifactor_ini,thetafactor_ini,betaPSII_ini,
|
|
&realizedfjelect,xvector(ntotsamples0),yvector(ntotsamples0),
|
|
&fvector(ntotsamples0),gvector(ntotsamples0),hvector(ntotsamples0),
|
|
&zvector(ntotsamples0),wvector(ntotsamples0),uvector(ntotsamples0),
|
|
&fo_dark,fm_dark,resp_dark,tempK_dark,ACimaxcurvature(15),
|
|
&ACimaxcurvpco2i(15),PhiPSIImaxcurvature_ACi(15),
|
|
&PhiPSIImaxcurv_ACi(15),ALightmaxcurvature(15),
|
|
&ALightmaxcurvPAR(15),PhiPSIImaxcurvature_ALight(15),
|
|
&PhiPSIImaxcurv_ALight(15),co2iRubismax25,co2iRuBpmax25,
|
|
&anetRubismax25,anetRuBpmax25,starco2a(15),der_starco2a(15),
|
|
&Amax_ACa(15),ACainter(15),der_ACainter(15),der_ACa400ppm(15),
|
|
&anet_ACa400ppm(15),PhiPSIImax_ACa(15),PhiPSIIinter_ACa(15),
|
|
&der_PhiPSIIinter_ACa(15),der_PhiPSIIend_ACa(15),
|
|
&ACamaxcurvature(15),ACamaxcurvpco2a(15),
|
|
&PhiPSIImaxcurvature_ACa(15),PhiPSIImaxcurv_ACa(15),
|
|
&PhiPSIIzero_ACa(15),der_PhiPSIIzero_ACa(15),ha_darkresp_ini,
|
|
&ha_stargamma_ini,ha_vcmax_ini,ha_jmax_ini,ha_tpu_ini,ha_gmeso_ini
|
|
parameter(critdelPAR=-2.0d0,critdelCi_Pa=-2.0d0)
|
|
!use positive critdelCi_Pa and critdelPAR to indicate absolute distance
|
|
!use negative critdelCi_Pa and critdelPAR to indicate relative distance (percentage value)
|
|
!End of declaration=======================================================================
|
|
paramunit=indexunit(1)
|
|
compareunit=indexunit(2)
|
|
stomwuecicaoutunit=indexunit(3)
|
|
stomcompunit=indexunit(4)
|
|
wuecicacompunit=indexunit(5)
|
|
fluorescenceunit=indexunit(6)
|
|
fluoresparamunit=indexunit(7)
|
|
aciempfitunit=indexunit(8)
|
|
alightempfitunit=indexunit(9)
|
|
!-----------------------------------------------------------------------------------------
|
|
call commonparameters(stargamma25_ini,fkc25_ini,fko25_ini,
|
|
&alpha25_ini,ha_vcmax_ini,hd_vcmax,sv_vcmax,ha_jmax_ini,hd_jmax,
|
|
&sv_jmax,ha_tpu_ini,hd_tpu,sv_tpu,ha_gmeso_ini,hd_gmeso,sv_gmeso,
|
|
&ha_darkresp_ini,ha_stargamma_ini,ha_kc,ha_ko,abspt_lf_par,
|
|
&gascon,phifactor_ini,thetafactor_ini,betaPSII_ini)
|
|
ha_darkresp=ha_darkresp_ini
|
|
ha_stargamma=ha_stargamma_ini
|
|
ha_vcmax=ha_vcmax_ini
|
|
ha_jmax=ha_jmax_ini
|
|
ha_tpu=ha_tpu_ini
|
|
ha_gmeso=ha_gmeso_ini
|
|
call pam_parameters(ntotsamples0,fo_pam0,fm_pam0,fs_pam0,
|
|
&pam_measlight0,anet_obs0,PARi0,templeaf0,yield_ps2,yield_npq,
|
|
&qlake,qpuddle,kps2_norm,knpq_norm,fo_dark,fm_dark,resp_dark,
|
|
&tempK_dark)
|
|
j=0
|
|
do i=1,ntotsamples0
|
|
!this is needed because the calling routine passes any data that have valid PAM measurements.
|
|
k=0
|
|
if(dabs(anet_obs0(i)+9999.0d0).lt.0.01d0)k=1
|
|
if(dabs(pco2i0(i)+9999.0d0).lt.0.01d0)k=1
|
|
if(dabs(templeaf0(i)+9999.0d0).lt.0.01d0)k=1
|
|
if(k.eq.0)then
|
|
j=j+1
|
|
anet_obs0(j)=anet_obs0(i)
|
|
pco2i0(j)=pco2i0(i)
|
|
templeaf0(j)=templeaf0(i)
|
|
PARi0(j)=PARi0(i)
|
|
pres_air0(j)=pres_air0(i)
|
|
po2i0(j)=po2i0(i)
|
|
chlflphips20(j)=chlflphips20(i)
|
|
pco2ambient0(j)=pco2ambient0(i)
|
|
trmmol0(j)=trmmol0(i)
|
|
gswmeas0(j)=gswmeas0(i)
|
|
vpdl0(j)=vpdl0(i)
|
|
tempair0(j)=tempair0(i)
|
|
eambient0(j)=eambient0(i)
|
|
!
|
|
fo_pam0(j)=fo_pam0(i)
|
|
fm_pam0(j)=fm_pam0(i)
|
|
fs_pam0(j)=fs_pam0(i)
|
|
pam_measlight0(j)=pam_measlight0(i)
|
|
yield_ps2(j)=yield_ps2(i)
|
|
yield_npq(j)=yield_npq(i)
|
|
qlake(j)=qlake(i)
|
|
qpuddle(j)=qpuddle(i)
|
|
kps2_norm(j)=kps2_norm(i)
|
|
knpq_norm(j)=knpq_norm(i)
|
|
endif
|
|
enddo
|
|
ntotsamples0=j
|
|
!
|
|
vcmax25_ini=50.0d0
|
|
fjmax25_ini=1.1d0*vcmax25_ini
|
|
tpu25_ini=0.07d0*fjmax25_ini
|
|
rdlight25_ini=0.015d0*vcmax25_ini
|
|
if(resp_dark.gt.0.0d0)then
|
|
!data contain dark-adapted rd
|
|
call resp_mitocho(tempK_dark,1.0d0,ha_darkresp,gascon,term)
|
|
rdlight25_ini=resp_dark/term
|
|
if(rdlight25_usr.le.0.0d0)rdlight25_usr=rdlight25_ini
|
|
endif
|
|
resistwp25_ini=0.1d0
|
|
resistch25_ini=0.1d0
|
|
resistwp25max=100.0d0
|
|
resistwp25min=0.0d0
|
|
resistch25max=100.0d0
|
|
resistch25min=0.0d0
|
|
rdlight25max=10.d0
|
|
rdlight25min=1.0d-7
|
|
stargamma25max=10.0d0
|
|
stargamma25min=1.0d-7
|
|
vcmax25max=700.0d0
|
|
vcmax25min=0.0d0
|
|
fkc25max=100.0d0
|
|
fkc25min=5.0d0
|
|
fko25max=20000.0d0
|
|
fko25min=10000.0d0
|
|
fjmax25max=800.0d0
|
|
fjmax25min=0.0d0
|
|
tpu25max=20.0d0
|
|
tpu25min=0.0d0
|
|
alpha25max=10.0d0
|
|
alpha25min=0.0d0
|
|
alpha25_ini=0.001d0
|
|
phifactormin=1.0d-5
|
|
phifactormax=2.0d0
|
|
thetafactormin=1.0d-5
|
|
thetafactormax=1.2d0
|
|
betaPSIImin=0.0d0
|
|
betaPSIImax=1.0d0
|
|
if(ha_darkresp.gt.0.0d0)then
|
|
ha_darkrespmin=5.0d0
|
|
ha_darkrespmax=200.0d0
|
|
else
|
|
!-Q10
|
|
ha_darkrespmin=-200.0d0
|
|
ha_darkrespmax=0.0d0
|
|
endif
|
|
ha_stargammamin=5.0d0
|
|
ha_stargammamax=200.0d0
|
|
ha_vcmaxmin=40.0d0
|
|
ha_vcmaxmax=100.0d0
|
|
ha_jmaxmin=20.0d0
|
|
ha_jmaxmax=100.0d0
|
|
ha_tpumin=20.0d0
|
|
ha_tpumax=100.0d0
|
|
ha_gmesomin=20.0d0
|
|
ha_gmesomax=100.0d0
|
|
if(isitmassbased.eq.1)then
|
|
vcmax25max=2000.0d0
|
|
fjmax25max=2000.0d0
|
|
tpu25max=100.0d0
|
|
rdlight25max=30.d0
|
|
endif
|
|
nFixedPoints=0
|
|
numACicurves=0
|
|
numALightcurves=0
|
|
nFreePoints=0
|
|
do i=1,ntotsamples0
|
|
aPPFDlf0(i)=PARi0(i)*abspt_lf_par
|
|
if(gswmeas0(i).gt.0.0d0)then
|
|
resiststomco20(i)=1.6d0/gswmeas0(i)
|
|
!unit is 1/(mol/m2/s). Now we need to change it to 1.0d0/(umol/m2/s/Pa)
|
|
resiststomco20(i)=resiststomco20(i)*pres_air0(i)*1.0d-6
|
|
else
|
|
resiststomco20(i)=-9999.0d0
|
|
endif
|
|
j=idnint(CurveTypeID(i)+0.1d0)
|
|
if(j.eq.1.or.j.eq.2.or.j.eq.3)then
|
|
!points whose limitation states are known.
|
|
nFixedPoints=nFixedPoints+1
|
|
Fixedanet_obs(nFixedPoints)=anet_obs0(i)
|
|
Fixedpco2i(nFixedPoints)=pco2i0(i)
|
|
Fixedtempleaf(nFixedPoints)=templeaf0(i)
|
|
FixedaPPFDlf(nFixedPoints)=aPPFDlf0(i)
|
|
Fixedpres_air(nFixedPoints)=pres_air0(i)
|
|
Fixedpo2i(nFixedPoints)=po2i0(i)
|
|
Fixedchlflphips2(nFixedPoints)=chlflphips20(i)
|
|
Fixedpco2ambient(nFixedPoints)=pco2ambient0(i)
|
|
Fixedtrmmol(nFixedPoints)=trmmol0(i)
|
|
Fixedgswmeas(nFixedPoints)=gswmeas0(i)
|
|
Fixedvpdl(nFixedPoints)=vpdl0(i)
|
|
Fixedtempair(nFixedPoints)=tempair0(i)
|
|
Fixedeambient(nFixedPoints)=eambient0(i)
|
|
!
|
|
Fixedfo_pam(nFixedPoints)=fo_pam0(i)
|
|
Fixedfm_pam(nFixedPoints)=fm_pam0(i)
|
|
Fixedfs_pam(nFixedPoints)=fs_pam0(i)
|
|
Fixedpam_measlight(nFixedPoints)=pam_measlight0(i)
|
|
Fixedyield_ps2(nFixedPoints)=yield_ps2(i)
|
|
Fixedyield_npq(nFixedPoints)=yield_npq(i)
|
|
Fixedqlake(nFixedPoints)=qlake(i)
|
|
Fixedqpuddle(nFixedPoints)=qpuddle(i)
|
|
Fixedkps2_norm(nFixedPoints)=kps2_norm(i)
|
|
Fixedknpq_norm(nFixedPoints)=knpq_norm(i)
|
|
!
|
|
Fixedresiststomco2(nFixedPoints)=resiststomco20(i)
|
|
Prioriphotolimit(nFixedPoints)=j
|
|
else
|
|
if(j.ge.11.and.j.le.25)then
|
|
!A/Ci curves without knowing limitation states of points.
|
|
m=0
|
|
do k=1,numACicurves
|
|
if(j.eq.ACiID(k))then
|
|
nACiPoints(k)=nACiPoints(k)+1
|
|
ACianet_obs0(nACiPoints(k),k)=anet_obs0(i)
|
|
ACipco2i0(nACiPoints(k),k)=pco2i0(i)
|
|
ACitempleaf0(nACiPoints(k),k)=templeaf0(i)
|
|
ACiaPPFDlf0(nACiPoints(k),k)=aPPFDlf0(i)
|
|
ACipres_air0(nACiPoints(k),k)=pres_air0(i)
|
|
ACipo2i0(nACiPoints(k),k)=po2i0(i)
|
|
ACichlflphips20(nACiPoints(k),k)=chlflphips20(i)
|
|
ACipco2ambient0(nACiPoints(k),k)=pco2ambient0(i)
|
|
ACitrmmol0(nACiPoints(k),k)=trmmol0(i)
|
|
ACigswmeas0(nACiPoints(k),k)=gswmeas0(i)
|
|
ACivpdl0(nACiPoints(k),k)=vpdl0(i)
|
|
ACitempair0(nACiPoints(k),k)=tempair0(i)
|
|
ACieambient0(nACiPoints(k),k)=eambient0(i)
|
|
!
|
|
ACifo_pam0(nACiPoints(k),k)=fo_pam0(i)
|
|
ACifm_pam0(nACiPoints(k),k)=fm_pam0(i)
|
|
ACifs_pam0(nACiPoints(k),k)=fs_pam0(i)
|
|
ACipam_measlight0(nACiPoints(k),k)=pam_measlight0(i)
|
|
ACiyield_ps20(nACiPoints(k),k)=yield_ps2(i)
|
|
ACiyield_npq0(nACiPoints(k),k)=yield_npq(i)
|
|
ACiqlake0(nACiPoints(k),k)=qlake(i)
|
|
ACiqpuddle0(nACiPoints(k),k)=qpuddle(i)
|
|
ACikps2_norm0(nACiPoints(k),k)=kps2_norm(i)
|
|
ACiknpq_norm0(nACiPoints(k),k)=knpq_norm(i)
|
|
!
|
|
ACiresiststomco20(nACiPoints(k),k)=resiststomco20(i)
|
|
m=1
|
|
endif
|
|
enddo
|
|
if(m.eq.0)then
|
|
!A new ACi curve
|
|
numACicurves=numACicurves+1
|
|
nACiPoints(numACicurves)=1
|
|
ACiID(numACicurves)=j
|
|
ACianet_obs0(1,numACicurves)=anet_obs0(i)
|
|
ACipco2i0(1,numACicurves)=pco2i0(i)
|
|
ACitempleaf0(1,numACicurves)=templeaf0(i)
|
|
ACiaPPFDlf0(1,numACicurves)=aPPFDlf0(i)
|
|
ACipres_air0(1,numACicurves)=pres_air0(i)
|
|
ACipo2i0(1,numACicurves)=po2i0(i)
|
|
ACichlflphips20(1,numACicurves)=chlflphips20(i)
|
|
ACipco2ambient0(1,numACicurves)=pco2ambient0(i)
|
|
ACitrmmol0(1,numACicurves)=trmmol0(i)
|
|
ACigswmeas0(1,numACicurves)=gswmeas0(i)
|
|
ACivpdl0(1,numACicurves)=vpdl0(i)
|
|
ACitempair0(1,numACicurves)=tempair0(i)
|
|
ACieambient0(1,numACicurves)=eambient0(i)
|
|
!
|
|
ACifo_pam0(1,numACicurves)=fo_pam0(i)
|
|
ACifm_pam0(1,numACicurves)=fm_pam0(i)
|
|
ACifs_pam0(1,numACicurves)=fs_pam0(i)
|
|
ACipam_measlight0(1,numACicurves)=pam_measlight0(i)
|
|
ACiyield_ps20(1,numACicurves)=yield_ps2(i)
|
|
ACiyield_npq0(1,numACicurves)=yield_npq(i)
|
|
ACiqlake0(1,numACicurves)=qlake(i)
|
|
ACiqpuddle0(1,numACicurves)=qpuddle(i)
|
|
ACikps2_norm0(1,numACicurves)=kps2_norm(i)
|
|
ACiknpq_norm0(1,numACicurves)=knpq_norm(i)
|
|
!
|
|
ACiresiststomco20(1,numACicurves)=resiststomco20(i)
|
|
endif
|
|
else
|
|
if(j.ge.31.and.j.le.45)then
|
|
!A/Light curves without knowing limitation states of points.
|
|
m=0
|
|
do k=1,numALightcurves
|
|
if(j.eq.ALightID(k))then
|
|
nALightPoints(k)=nALightPoints(k)+1
|
|
ALightanet_obs0(nALightPoints(k),k)=anet_obs0(i)
|
|
ALightpco2i0(nALightPoints(k),k)=pco2i0(i)
|
|
ALighttempleaf0(nALightPoints(k),k)=templeaf0(i)
|
|
ALightaPPFDlf0(nALightPoints(k),k)=aPPFDlf0(i)
|
|
ALightpres_air0(nALightPoints(k),k)=pres_air0(i)
|
|
ALightpo2i0(nALightPoints(k),k)=po2i0(i)
|
|
ALightchlflphips20(nALightPoints(k),k)=chlflphips20(i)
|
|
ALightpco2ambient0(nALightPoints(k),k)=pco2ambient0(i)
|
|
ALighttrmmol0(nALightPoints(k),k)=trmmol0(i)
|
|
ALightgswmeas0(nALightPoints(k),k)=gswmeas0(i)
|
|
ALightvpdl0(nALightPoints(k),k)=vpdl0(i)
|
|
ALighttempair0(nALightPoints(k),k)=tempair0(i)
|
|
ALighteambient0(nALightPoints(k),k)=eambient0(i)
|
|
!
|
|
ALightfo_pam0(nALightPoints(k),k)=fo_pam0(i)
|
|
ALightfm_pam0(nALightPoints(k),k)=fm_pam0(i)
|
|
ALightfs_pam0(nALightPoints(k),k)=fs_pam0(i)
|
|
ALightpam_measlight0(nALightPoints(k),k)=
|
|
&pam_measlight0(i)
|
|
ALightyield_ps20(nALightPoints(k),k)=yield_ps2(i)
|
|
ALightyield_npq0(nALightPoints(k),k)=yield_npq(i)
|
|
ALightqlake0(nALightPoints(k),k)=qlake(i)
|
|
ALightqpuddle0(nALightPoints(k),k)=qpuddle(i)
|
|
ALightkps2_norm0(nALightPoints(k),k)=kps2_norm(i)
|
|
ALightknpq_norm0(nALightPoints(k),k)=knpq_norm(i)
|
|
!
|
|
ALightresiststomco20(nALightPoints(k),k)=
|
|
&resiststomco20(i)
|
|
m=1
|
|
endif
|
|
enddo
|
|
if(m.eq.0)then
|
|
!A new A/Light curve
|
|
numALightcurves=numALightcurves+1
|
|
nALightPoints(numALightcurves)=1
|
|
ALightID(numALightcurves)=j
|
|
ALightanet_obs0(1,numALightcurves)=anet_obs0(i)
|
|
ALightpco2i0(1,numALightcurves)=pco2i0(i)
|
|
ALighttempleaf0(1,numALightcurves)=templeaf0(i)
|
|
ALightaPPFDlf0(1,numALightcurves)=aPPFDlf0(i)
|
|
ALightpres_air0(1,numALightcurves)=pres_air0(i)
|
|
ALightpo2i0(1,numALightcurves)=po2i0(i)
|
|
ALightchlflphips20(1,numALightcurves)=chlflphips20(i)
|
|
ALightpco2ambient0(1,numALightcurves)=pco2ambient0(i)
|
|
ALighttrmmol0(1,numALightcurves)=trmmol0(i)
|
|
ALightgswmeas0(1,numALightcurves)=gswmeas0(i)
|
|
ALightvpdl0(1,numALightcurves)=vpdl0(i)
|
|
ALighttempair0(1,numALightcurves)=tempair0(i)
|
|
ALighteambient0(1,numALightcurves)=eambient0(i)
|
|
!
|
|
ALightfo_pam0(1,numALightcurves)=fo_pam0(i)
|
|
ALightfm_pam0(1,numALightcurves)=fm_pam0(i)
|
|
ALightfs_pam0(1,numALightcurves)=fs_pam0(i)
|
|
ALightpam_measlight0(1,numALightcurves)=
|
|
&pam_measlight0(i)
|
|
ALightyield_ps20(1,numALightcurves)=yield_ps2(i)
|
|
ALightyield_npq0(1,numALightcurves)=yield_npq(i)
|
|
ALightqlake0(1,numALightcurves)=qlake(i)
|
|
ALightqpuddle0(1,numALightcurves)=qpuddle(i)
|
|
ALightkps2_norm0(1,numALightcurves)=kps2_norm(i)
|
|
ALightknpq_norm0(1,numALightcurves)=knpq_norm(i)
|
|
!
|
|
ALightresiststomco20(1,numALightcurves)=
|
|
&resiststomco20(i)
|
|
endif
|
|
else
|
|
nFreePoints=nFreePoints+1
|
|
Freeanet_obs(nFreePoints)=anet_obs0(i)
|
|
Freepco2i(nFreePoints)=pco2i0(i)
|
|
Freetempleaf(nFreePoints)=templeaf0(i)
|
|
FreeaPPFDlf(nFreePoints)=aPPFDlf0(i)
|
|
Freepres_air(nFreePoints)=pres_air0(i)
|
|
Freepo2i(nFreePoints)=po2i0(i)
|
|
Freechlflphips2(nFreePoints)=chlflphips20(i)
|
|
Freepco2ambient(nFreePoints)=pco2ambient0(i)
|
|
Freetrmmol(nFreePoints)=trmmol0(i)
|
|
Freegswmeas(nFreePoints)=gswmeas0(i)
|
|
Freevpdl(nFreePoints)=vpdl0(i)
|
|
Freetempair(nFreePoints)=tempair0(i)
|
|
Freeeambient(nFreePoints)=eambient0(i)
|
|
!
|
|
Freefo_pam(nFreePoints)=fo_pam0(i)
|
|
Freefm_pam(nFreePoints)=fm_pam0(i)
|
|
Freefs_pam(nFreePoints)=fs_pam0(i)
|
|
Freepam_measlight(nFreePoints)=pam_measlight0(i)
|
|
Freeyield_ps2(nFreePoints)=yield_ps2(i)
|
|
Freeyield_npq(nFreePoints)=yield_npq(i)
|
|
Freeqlake(nFreePoints)=qlake(i)
|
|
Freeqpuddle(nFreePoints)=qpuddle(i)
|
|
Freekps2_norm(nFreePoints)=kps2_norm(i)
|
|
Freeknpq_norm(nFreePoints)=knpq_norm(i)
|
|
!
|
|
Freeresiststomco2(nFreePoints)=resiststomco20(i)
|
|
endif
|
|
endif
|
|
endif
|
|
enddo
|
|
!-----------------------------------------------------------------------
|
|
!Average clusters and then sort of ACi and ALight points. No need to cluster or sort fixed and free points
|
|
do i=1,numACicurves
|
|
call clustering(nACiPoints(i),1,ACipco2i0(1:nACiPoints(i),i:i),
|
|
&critdelCi_Pa,k,ibelong)
|
|
if(k.lt.nACipoints(i))then
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACipco2i0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACipco2ambient0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACiaPPFDlf0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACitempleaf0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACipres_air0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACianet_obs0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACipo2i0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACitrmmol0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACigswmeas0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACivpdl0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACitempair0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACieambient0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACichlflphips20(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
!
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACifo_pam0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACifm_pam0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACifs_pam0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACipam_measlight0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACiyield_ps20(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACiyield_npq0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACiqlake0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACiqpuddle0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACikps2_norm0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACiknpq_norm0(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
!
|
|
call aftercluster(nACiPoints(i),1,
|
|
&ACiresiststomco20(1:nACiPoints(i),i:i),k,ibelong,fvector)
|
|
nACiPoints(i)=k
|
|
endif
|
|
!sort CO2i from low to high
|
|
do j=1,nACiPoints(i)
|
|
ACipco2i(j,i)=ACipco2i0(j,i)
|
|
enddo
|
|
call sort_shell(nACiPoints(i),ACipco2i(1:nACiPoints(i),i:i),
|
|
&ioriorder)
|
|
do j=1,nACiPoints(i)
|
|
ACianet_obs(j,i)=ACianet_obs0(ioriorder(j),i)
|
|
ACitempleaf(j,i)=ACitempleaf0(ioriorder(j),i)
|
|
ACiaPPFDlf(j,i)=ACiaPPFDlf0(ioriorder(j),i)
|
|
ACipo2i(j,i)=ACipo2i0(ioriorder(j),i)
|
|
ACipres_air(j,i)=ACipres_air0(ioriorder(j),i)
|
|
ACipco2ambient(j,i)=ACipco2ambient0(ioriorder(j),i)
|
|
ACitrmmol(j,i)=ACitrmmol0(ioriorder(j),i)
|
|
ACigswmeas(j,i)=ACigswmeas0(ioriorder(j),i)
|
|
ACivpdl(j,i)=ACivpdl0(ioriorder(j),i)
|
|
ACitempair(j,i)=ACitempair0(ioriorder(j),i)
|
|
ACieambient(j,i)=ACieambient0(ioriorder(j),i)
|
|
ACichlflphips2(j,i)=ACichlflphips20(ioriorder(j),i)
|
|
!
|
|
ACifo_pam(j,i)=ACifo_pam0(ioriorder(j),i)
|
|
ACifm_pam(j,i)=ACifm_pam0(ioriorder(j),i)
|
|
ACifs_pam(j,i)=ACifs_pam0(ioriorder(j),i)
|
|
ACipam_measlight(j,i)=ACipam_measlight0(ioriorder(j),i)
|
|
ACiyield_ps2(j,i)=ACiyield_ps20(ioriorder(j),i)
|
|
ACiyield_npq(j,i)=ACiyield_npq0(ioriorder(j),i)
|
|
ACiqlake(j,i)=ACiqlake0(ioriorder(j),i)
|
|
ACiqpuddle(j,i)=ACiqpuddle0(ioriorder(j),i)
|
|
ACikps2_norm(j,i)=ACikps2_norm0(ioriorder(j),i)
|
|
ACiknpq_norm(j,i)=ACiknpq_norm0(ioriorder(j),i)
|
|
!
|
|
ACiresiststomco2(j,i)=ACiresiststomco20(ioriorder(j),i)
|
|
enddo
|
|
enddo
|
|
do i=1,numALightcurves
|
|
call clustering(nALightPoints(i),1,
|
|
&ALightaPPFDlf0(1:nALightPoints(i),i:i),critdelPAR,k,ibelong)
|
|
if(k.lt.nALightpoints(i))then
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightpco2i0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightpco2ambient0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightaPPFDlf0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALighttempleaf0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightpres_air0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightanet_obs0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightpo2i0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALighttrmmol0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightgswmeas0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightvpdl0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALighttempair0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALighteambient0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightchlflphips20(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
!
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightfo_pam0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightfm_pam0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightfs_pam0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightpam_measlight0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightyield_ps20(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightyield_npq0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightqlake0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightqpuddle0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightkps2_norm0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightknpq_norm0(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
!
|
|
call aftercluster(nALightPoints(i),1,
|
|
&ALightresiststomco20(1:nALightPoints(i),i:i),k,ibelong,fvector)
|
|
nALightPoints(i)=k
|
|
endif
|
|
!sort PAR from low to high
|
|
do j=1,nALightPoints(i)
|
|
ALightaPPFDlf(j,i)=ALightaPPFDlf0(j,i)
|
|
enddo
|
|
call sort_shell(nALightPoints(i),
|
|
&ALightaPPFDlf(1:nALightPoints(i),i:i),ioriorder)
|
|
do j=1,nALightPoints(i)
|
|
ALightanet_obs(j,i)=ALightanet_obs0(ioriorder(j),i)
|
|
ALighttempleaf(j,i)=ALighttempleaf0(ioriorder(j),i)
|
|
ALightpco2i(j,i)=ALightpco2i0(ioriorder(j),i)
|
|
ALightpo2i(j,i)=ALightpo2i0(ioriorder(j),i)
|
|
ALightpres_air(j,i)=ALightpres_air0(ioriorder(j),i)
|
|
ALightpco2ambient(j,i)=ALightpco2ambient0(ioriorder(j),i)
|
|
ALighttrmmol(j,i)=ALighttrmmol0(ioriorder(j),i)
|
|
ALightgswmeas(j,i)=ALightgswmeas0(ioriorder(j),i)
|
|
ALightvpdl(j,i)=ALightvpdl0(ioriorder(j),i)
|
|
ALighttempair(j,i)=ALighttempair0(ioriorder(j),i)
|
|
ALighteambient(j,i)=ALighteambient0(ioriorder(j),i)
|
|
ALightchlflphips2(j,i)=ALightchlflphips20(ioriorder(j),i)
|
|
!
|
|
ALightfo_pam(j,i)=ALightfo_pam0(ioriorder(j),i)
|
|
ALightfm_pam(j,i)=ALightfm_pam0(ioriorder(j),i)
|
|
ALightfs_pam(j,i)=ALightfs_pam0(ioriorder(j),i)
|
|
ALightpam_measlight(j,i)=ALightpam_measlight0(ioriorder(j),i)
|
|
ALightyield_ps2(j,i)=ALightyield_ps20(ioriorder(j),i)
|
|
ALightyield_npq(j,i)=ALightyield_npq0(ioriorder(j),i)
|
|
ALightqlake(j,i)=ALightqlake0(ioriorder(j),i)
|
|
ALightqpuddle(j,i)=ALightqpuddle0(ioriorder(j),i)
|
|
ALightkps2_norm(j,i)=ALightkps2_norm0(ioriorder(j),i)
|
|
ALightknpq_norm(j,i)=ALightknpq_norm0(ioriorder(j),i)
|
|
!
|
|
ALightresiststomco2(j,i)=ALightresiststomco20(ioriorder(j),i)
|
|
enddo
|
|
enddo
|
|
!-----------------------------------------------------------------------
|
|
idoalpha0=1
|
|
do i=1,numACicurves
|
|
amaxave=0.0d0
|
|
n=3
|
|
do j=nACiPoints(i)-n+1,nACiPoints(i)
|
|
amaxave=amaxave+ACianet_obs(j,i)
|
|
enddo
|
|
amaxave=amaxave/dble(n)
|
|
!the sigmoidal function has better asymptotic behaviour so
|
|
!it is used for estimating anetmaxs.
|
|
iderivative=1
|
|
INFO=0
|
|
!INFO =0, ordinary distance regression
|
|
!INFO =1, explicit orthogonal distance regression with shortest distance within iteration
|
|
!INFO =2, explicit orthogonal distance regression with x positions as parameters
|
|
beta(1)=dabs(amaxave)
|
|
if(amaxave.lt.0.0d0)then
|
|
betamin(1)=amaxave
|
|
else
|
|
betamin(1)=0.5d0*amaxave
|
|
endif
|
|
betamax(1)=1000.0d0
|
|
beta(2)=1.5d0
|
|
betamin(2)=1.0d-5
|
|
betamax(2)=1000.0d0
|
|
beta(3)=0.1d0
|
|
betamin(3)=0.0d0
|
|
betamax(3)=100.0d0
|
|
beta(4)=30.0d0
|
|
betamin(4)=0.0d0
|
|
betamax(4)=5000.0d0
|
|
beta(5)=-10.0d0
|
|
betamin(5)=-1000.0d0
|
|
betamax(5)=1000.0d0
|
|
k=0
|
|
n=0
|
|
do j=1,nACiPoints(i)
|
|
weitx(j)=1.0d0
|
|
xmin(j)=dmax1(0.0d0,ACipco2i(j,i)-20.0d0)
|
|
xmax(j)=ACipco2i(j,i)+20.0d0
|
|
weity(j)=1.0d0
|
|
if(ACichlflphips2(j,i).gt.0.0d0)then
|
|
k=k+1
|
|
yvector(k)=ACichlflphips2(j,i)
|
|
xvector(k)=ACipco2i(j,i)
|
|
uvector(k)=ACipco2ambient(j,i)
|
|
endif
|
|
if(ACipco2ambient(j,i).gt.0.0d0)then
|
|
n=n+1
|
|
zvector(n)=ACianet_obs(j,i)
|
|
wvector(n)=ACipco2ambient(j,i)
|
|
endif
|
|
enddo
|
|
!
|
|
call GenericRegres(nACiPoints(i),1,
|
|
&ACianet_obs(1:nACiPoints(i),i:i),1,ACipco2i(1:nACiPoints(i),i:i),
|
|
&weity,weitx,5,beta,betamin,betamax,xmin,xmax,iderivative,INFO,
|
|
&fvector,gvector,sumsquare)
|
|
call properties_surffunc(5,beta,starco2i(i),der_starco2i(i),
|
|
&Amax_ACi(i),ACiinter(i),der_ACiinter(i),
|
|
&ACipco2i(nACiPoints(i),i),der_ACiend(i),term,
|
|
! &ACipco2i(nACiPoints(i):nACiPoints(i),i:i),der_ACiend(i),term,
|
|
&ACipco2i(1:1,i:i),ACipco2i(nACiPoints(i):nACiPoints(i),i:i),
|
|
&ACimaxcurvature(i),ACimaxcurvpco2i(i))
|
|
if(n.ge.5)then
|
|
call GenericRegres(n,1,zvector,1,wvector,weity,weitx,5,beta,
|
|
&betamin,betamax,xmin,xmax,iderivative,INFO,fvector,gvector,
|
|
&sumsquare)
|
|
call properties_surffunc(5,beta,starco2a(i),der_starco2a(i),
|
|
&Amax_ACa(i),ACainter(i),der_ACainter(i),40.0d0,der_ACa400ppm(i),
|
|
&anet_ACa400ppm(i),wvector(1),wvector(n),ACamaxcurvature(i),
|
|
&ACamaxcurvpco2a(i))
|
|
else
|
|
der_starco2a(i)=-9999.0d0
|
|
Amax_ACa(i)=-9999.0d0
|
|
ACainter(i)=-9999.0d0
|
|
der_ACainter(i)=-9999.0d0
|
|
der_ACa400ppm(i)=-9999.0d0
|
|
anet_ACa400ppm(i)=-9999.0d0
|
|
ACamaxcurvature(i)=-9999.0d0
|
|
ACamaxcurvpco2a(i)=-9999.0d0
|
|
endif
|
|
if(Amax_ACi(i).lt.50.0d0)amaxave=Amax_ACi(i)
|
|
j=min0(5,nACiPoints(i))
|
|
call y_aPLUSbx(j,ACipco2i(1:j,i:i),ACianet_obs(1:j,i:i),ac,at)
|
|
!fit for y=ac+at*x
|
|
if(ac.lt.0.0d0.and.dabs(ac).lt.rdlight25max)then
|
|
rdlight25_ini=dabs(ac)
|
|
if((-ac/at).lt.stargamma25max.and.
|
|
&(-ac/at).gt.stargamma25_ini)stargamma25max=-ac/at
|
|
endif
|
|
if(amaxave.gt.0.0d0)then
|
|
fjmax25_ini=(amaxave+rdlight25_ini)*4.0d0+10.0d0
|
|
vcmax25_ini=fjmax25_ini/1.1d0
|
|
tpu25_ini=(amaxave+rdlight25_ini)/3.0d0
|
|
else
|
|
fjmax25_ini=10.0d0
|
|
vcmax25_ini=fjmax25_ini/1.1d0
|
|
tpu25_ini=1.0d0
|
|
endif
|
|
if(k.ge.5)then
|
|
! beta(1)=0.50d0
|
|
! betamin(1)=0.0d0
|
|
! betamax(1)=1000.0d0
|
|
! beta(2)=5.50d0
|
|
! betamin(2)=0.0d0
|
|
! betamax(2)=1000.0d0
|
|
! beta(3)=1.50d0
|
|
! betamin(3)=-10.0d0
|
|
! betamax(3)=10.0d0
|
|
|
|
beta(1)=0.4d0
|
|
betamin(1)=0.0d0
|
|
betamax(1)=2.0d0
|
|
beta(2)=1.5d0
|
|
betamin(2)=1.0d-5
|
|
betamax(2)=1000.0d0
|
|
beta(3)=0.1d0
|
|
betamin(3)=0.0d0
|
|
betamax(3)=100.0d0
|
|
beta(4)=30.0d0
|
|
betamin(4)=0.0d0
|
|
betamax(4)=5000.0d0
|
|
beta(5)=0.1d0
|
|
betamin(5)=-5.0d0
|
|
betamax(5)=5.0d0
|
|
do j=1,k
|
|
xmin(j)=dmax1(0.0d0,xvector(j)-20.0d0)
|
|
xmax(j)=xvector(j)+20.0d0
|
|
enddo
|
|
call GenericRegres(k,1,yvector,1,xvector,weity,weitx,5,
|
|
&beta,betamin,betamax,xmin,xmax,iderivative,INFO,fvector,gvector,
|
|
&sumsquare)
|
|
call properties_surffunc(5,beta,PhiPSIIzero_ACi(i),
|
|
&der_PhiPSIIzero_ACi(i),PhiPSIImax_ACi(i),
|
|
&PhiPSIIinter_ACi(i),der_PhiPSIIinter_ACi(i),xvector(k),
|
|
&der_PhiPSIIend_ACi(i),term,xvector(1),xvector(k),
|
|
&PhiPSIImaxcurvature_ACi(i),PhiPSIImaxcurv_ACi(i))
|
|
call GenericRegres(k,1,yvector,1,uvector,weity,weitx,5,
|
|
&beta,betamin,betamax,xmin,xmax,iderivative,INFO,fvector,gvector,
|
|
&sumsquare)
|
|
call properties_surffunc(5,beta,PhiPSIIzero_ACa(i),
|
|
&der_PhiPSIIzero_ACa(i),PhiPSIImax_ACa(i),
|
|
&PhiPSIIinter_ACa(i),der_PhiPSIIinter_ACa(i),uvector(k),
|
|
&der_PhiPSIIend_ACa(i),term,uvector(1),uvector(k),
|
|
&PhiPSIImaxcurvature_ACa(i),PhiPSIImaxcurv_ACa(i))
|
|
else
|
|
PhiPSIIinter_ACi(i)=-9999.0d0
|
|
der_PhiPSIIinter_ACi(i)=-9999.0d0
|
|
PhiPSIIzero_ACi(i)=-9999.0d0
|
|
der_PhiPSIIzero_ACi(i)=-9999.0d0
|
|
PhiPSIImax_ACi(i)=-9999.0d0
|
|
der_PhiPSIIend_ACi(i)=-9999.0d0
|
|
PhiPSIImaxcurvature_ACi(i)=-9999.0d0
|
|
PhiPSIImaxcurv_ACi(i)=-9999.0d0
|
|
PhiPSIIinter_ACa(i)=-9999.0d0
|
|
der_PhiPSIIinter_ACa(i)=-9999.0d0
|
|
PhiPSIIzero_ACa(i)=-9999.0d0
|
|
der_PhiPSIIzero_ACa(i)=-9999.0d0
|
|
PhiPSIImax_ACa(i)=-9999.0d0
|
|
der_PhiPSIIend_ACa(i)=-9999.0d0
|
|
PhiPSIImaxcurvature_ACa(i)=-9999.0d0
|
|
PhiPSIImaxcurv_ACa(i)=-9999.0d0
|
|
endif
|
|
!
|
|
n=nACiPoints(i)
|
|
call y_aPLUSbxrsq(n,ACipco2i(1:n,i:i),ACianet_obs(1:n,i:i),
|
|
&ac,at,term)
|
|
resistwp25_ini=3.0d0*term**6
|
|
resistch25_ini=term**6
|
|
|
|
if(term.lt.0.9d0)then
|
|
if(Amax_ACi(i).gt.0.0d0.and.Amax_ACi(i).lt.100.0d0)then
|
|
resistwp25_ini=
|
|
&resistwp25_ini*dmin1(20.0d0/Amax_ACi(i),3.0d0)
|
|
resistch25_ini=
|
|
&resistch25_ini*dmin1(20.0d0/Amax_ACi(i),2.0d0)
|
|
else
|
|
if(Amax_ACi(i).le.0.0d0)then
|
|
resistwp25_ini=6.0d0
|
|
resistch25_ini=4.0d0
|
|
endif
|
|
endif
|
|
endif
|
|
!almost a straightline
|
|
!determine the absolute last point of rubisco or rubp for an A/Ci curve
|
|
k=4
|
|
10 if(n.le.k)goto 20
|
|
if(ACianet_obs(n,i).gt.ACianet_obs(n-1,i).and.
|
|
&ACianet_obs(n-1,i).gt.ACianet_obs(n-2,i))goto 20
|
|
do j=1,k
|
|
gvector(j)=ACipco2i(n-j+1,i)
|
|
fvector(j)=ACianet_obs(n-j+1,i)
|
|
enddo
|
|
call y_aPLUSbx(k,gvector,fvector,ac,at)
|
|
!fit for y=a+bx
|
|
if(at.gt.0.0d0)goto 20
|
|
n=n-1
|
|
goto 10
|
|
20 nendaci(i)=n
|
|
if(ACianet_obs(n,i).le.ACianet_obs(n-1,i).and.
|
|
&ACianet_obs(n-1,i).le.ACianet_obs(n-2,i))nendaci(i)=nendaci(i)-1
|
|
n=nACiPoints(i)-nendaci(i)
|
|
if(n.ge.3)then
|
|
do j=1,n
|
|
gvector(j)=ACipco2i(nendaci(i)+j,i)
|
|
fvector(j)=ACianet_obs(nendaci(i)+j,i)
|
|
enddo
|
|
call y_aPLUSbx(n,gvector,fvector,ac,at)
|
|
!fit for y=a+bx
|
|
if(dabs(at).le.1.0d-5)idoalpha0=0
|
|
endif
|
|
!Beyond nendaci, the points can only be limited by TPU
|
|
!
|
|
!Determine the point before which all points are limited by Rubisco and/or RuBP regeneration and after which some points might be
|
|
!limited by Rubisco and/or RuBP regeneration and/or TPU until nendaci after which all points are limited by TPU.
|
|
n=1
|
|
aj=-1.0d+20
|
|
22 if(n.ge.(nendaci(i)-3))goto 24
|
|
do j=1,k
|
|
gvector(j)=ACipco2i(n+j-1,i)
|
|
fvector(j)=ACianet_obs(n+j-1,i)
|
|
enddo
|
|
call y_aPLUSbx(k,gvector,fvector,ac,at)
|
|
!fit for y=a+bx
|
|
if(at.le.0.0d0)goto 24
|
|
if(at.gt.aj)then
|
|
aj=at
|
|
else
|
|
if(at.lt.aj/5.0d0)goto 24
|
|
endif
|
|
n=n+1
|
|
goto 22
|
|
24 nstartaci(i)=n-1
|
|
!
|
|
n=nACiPoints(i)
|
|
if(n.ge.4)then
|
|
if(ACianet_obs(n,i).gt.ACianet_obs(n-1,i).and.
|
|
&ACianet_obs(n-1,i).gt.ACianet_obs(n-2,i).and.
|
|
&ACianet_obs(n-2,i).gt.ACianet_obs(n-3,i))then
|
|
nstartaci(i)=n-1
|
|
nendaci(i)=n
|
|
!only the last point can be possibly tpu
|
|
endif
|
|
endif
|
|
if((nendaci(i)-nstartaci(i)).le.2)goto 29
|
|
25 n=nstartaci(i)
|
|
if(ACianet_obs(n+1,i).gt.ACianet_obs(n,i))then
|
|
!if anet continues to increase, the point is not tpu limited
|
|
if((nendaci(i)-nstartaci(i)).gt.2)then
|
|
nstartaci(i)=n+1
|
|
goto 25
|
|
endif
|
|
else
|
|
nstartaci(i)=nstartaci(i)-1
|
|
nstartaci(i)=max0(nstartaci(i),0)
|
|
endif
|
|
29 continue
|
|
!before nstartaci, no TPU points can occur
|
|
enddo
|
|
!
|
|
do i=1,numALightcurves
|
|
amaxave=0.0d0
|
|
n=3
|
|
do j=nALightPoints(i)-n+1,nALightPoints(i)
|
|
amaxave=amaxave+ALightanet_obs(j,i)
|
|
enddo
|
|
amaxave=amaxave/dble(n)
|
|
!the sigmoidal function has better asymptotic behaviour so
|
|
!it is used for estimating anetmaxs.
|
|
iderivative=1
|
|
INFO=0
|
|
beta(1)=dabs(amaxave)
|
|
if(amaxave.lt.0.0d0)then
|
|
betamin(1)=amaxave
|
|
else
|
|
betamin(1)=0.5d0*amaxave
|
|
endif
|
|
betamax(1)=200.0d0
|
|
beta(2)=1.5d0
|
|
betamin(2)=1.0d-5
|
|
betamax(2)=1.0d+5
|
|
beta(3)=0.1d0
|
|
betamin(3)=0.0d0
|
|
betamax(3)=5000.0d0
|
|
beta(4)=30.0d0
|
|
betamin(4)=-1000.0d0
|
|
betamax(4)=1000.0d0
|
|
beta(5)=-10.0d0
|
|
betamin(5)=-100.0d0
|
|
betamax(5)=100.0d0
|
|
k=0
|
|
do j=1,nALightPoints(i)
|
|
hvector(j)=ALightaPPFDlf(j,i)/abspt_lf_par
|
|
weitx(j)=1.0d0
|
|
xmin(j)=dmax1(0.0d0,hvector(j)-20.0d0)
|
|
xmax(j)=hvector(j)+20.0d0
|
|
weity(j)=1.0d0
|
|
if(ALightchlflphips2(j,i).gt.0.0d0)then
|
|
k=k+1
|
|
yvector(k)=ALightchlflphips2(j,i)
|
|
xvector(k)=hvector(j)
|
|
endif
|
|
enddo
|
|
call GenericRegres(nALightPoints(i),1,
|
|
&ALightanet_obs(1:nALightPoints(i),i:i),1,hvector(i),weity,weitx,5,
|
|
&beta,betamin,betamax,xmin,xmax,iderivative,INFO,fvector,gvector,
|
|
&sumsquare)
|
|
call properties_surffunc(5,beta,starPAR(i),der_starPAR(i),
|
|
&Asat_ALight(i),ALightinter(i),der_ALightinter(i),
|
|
&hvector(nALightPoints(i)),der_ALightend(i),term,
|
|
&hvector(1),hvector(nALightPoints(i)),ALightmaxcurvature(i),
|
|
&ALightmaxcurvPAR(i))
|
|
if(Asat_ALight(i).lt.50.0d0)amaxave=Asat_ALight(i)
|
|
j=min0(5,nALightPoints(i))
|
|
call y_aPLUSbx(j,hvector(1:j),ALightanet_obs(1:j,i:i),ac,at)
|
|
!fit for y=ac+at*x
|
|
if(ac.lt.0.0d0.and.dabs(ac).lt.rdlight25max)
|
|
&rdlight25_ini=dabs(ac)
|
|
if(amaxave.gt.0.0d0)then
|
|
fjmax25_ini=(amaxave+rdlight25_ini)*4.0d0+10.0d0
|
|
vcmax25_ini=fjmax25_ini/1.1d0
|
|
tpu25_ini=(amaxave+rdlight25_ini)/3.0d0
|
|
else
|
|
fjmax25_ini=10.0d0
|
|
vcmax25_ini=fjmax25_ini/1.1d0
|
|
tpu25_ini=1.0d0
|
|
endif
|
|
if(k.ge.5)then
|
|
beta(1)=0.50d0
|
|
betamin(1)=-1000.0d0
|
|
betamax(1)=0.0d0
|
|
beta(2)=5.50d0
|
|
betamin(2)=0.0d0
|
|
betamax(2)=1000.0d0
|
|
beta(3)=1.50d0
|
|
betamin(3)=-10.0d0
|
|
betamax(3)=10.0d0
|
|
do j=1,k
|
|
xmin(j)=dmax1(0.0d0,xvector(j)-20.0d0)
|
|
xmax(j)=xvector(j)+20.0d0
|
|
enddo
|
|
call GenericRegres(k,1,yvector,1,xvector,weity,weitx,3,
|
|
&beta,betamin,betamax,xmin,xmax,iderivative,INFO,fvector,gvector,
|
|
&sumsquare)
|
|
call properties_surffunc(3,beta,term,term1,term2,
|
|
&PhiPSIIinter_ALight(i),der_PhiPSIIinter_ALight(i),
|
|
&1000.0d0,der_PhiPSII1000_ALight(i),ExcessLightFactor(i),
|
|
&xvector(1),xvector(k),PhiPSIImaxcurvature_ALight(i),
|
|
&PhiPSIImaxcurv_ALight(i))
|
|
der_PhiPSIIinter_ALight(i)=der_PhiPSIIinter_ALight(i)*1000.0d0
|
|
der_PhiPSII1000_ALight(i)=der_PhiPSII1000_ALight(i)*1000.0d0
|
|
ExcessLightFactor(i)=1.0d0-ExcessLightFactor(i)/0.83d0
|
|
else
|
|
PhiPSIIinter_ALight(i)=-9999.0d0
|
|
der_PhiPSIIinter_ALight(i)=-9999.0d0
|
|
der_PhiPSII1000_ALight(i)=-9999.0d0
|
|
ExcessLightFactor(i)=-9999.0d0
|
|
PhiPSIImaxcurvature_ALight(i)=-9999.0d0
|
|
PhiPSIImaxcurv_ALight(i)=-9999.0d0
|
|
endif
|
|
!determine the absolute last point of rubp for an A/Light curve
|
|
k=4
|
|
n=nALightPoints(i)
|
|
30 if(n.le.k)goto 40
|
|
if(ALightanet_obs(n,i).gt.ALightanet_obs(n-1,i).and.
|
|
&ALightanet_obs(n-1,i).gt.ALightanet_obs(n-2,i))goto 40
|
|
do j=1,k
|
|
gvector(j)=ALightaPPFDlf(n-j+1,i)
|
|
fvector(j)=ALightanet_obs(n-j+1,i)
|
|
enddo
|
|
call y_aPLUSbx(k,gvector,fvector,ac,at)
|
|
!fit for y=ac+at*x
|
|
if(at.gt.0.0d0)goto 40
|
|
n=n-1
|
|
goto 30
|
|
40 nendalight(i)=n
|
|
if(ALightanet_obs(n,i).le.ALightanet_obs(n-1,i).and.
|
|
&ALightanet_obs(n-1,i).le.ALightanet_obs(n-2,i))
|
|
&nendalight(i)=nendalight(i)-1
|
|
!Beyond nendalight, the points can only be limited by Rubisco or TPU because they have constant or decreasing anet with inceased light
|
|
!
|
|
!Determine the point before which all points are limited by RuBP regeneration and after which some points might be limited by RuBP until
|
|
!nendalight.
|
|
n=1
|
|
aj=-1.0d+20
|
|
50 if(n.ge.(nendalight(i)-3))goto 55
|
|
do j=1,k
|
|
gvector(j)=ALightaPPFDlf(n+j-1,i)
|
|
fvector(j)=ALightanet_obs(n+j-1,i)
|
|
enddo
|
|
call y_aPLUSbx(k,gvector,fvector,ac,at)
|
|
!fit for y=ac+at*x
|
|
if(at.lt.1.0d-4)goto 55
|
|
if(at.gt.aj)then
|
|
aj=at
|
|
else
|
|
if(at.lt.aj/5.0d0)goto 55
|
|
endif
|
|
n=n+1
|
|
goto 50
|
|
55 if(n.ge.(nendalight(i)-1))then
|
|
n=nendalight(i)-1
|
|
goto 56
|
|
endif
|
|
if(ALightanet_obs(n,i).lt.ALightanet_obs(n+1,i))then
|
|
n=n+1
|
|
goto 55
|
|
endif
|
|
56 nstartalight(i)=n-1
|
|
!before nstartalight, no rubisco or tpu points can occur because anet increases with increased light, indicating RuBP regeneration
|
|
!limitation
|
|
!
|
|
n=nALightPoints(i)
|
|
if(n.ge.4)then
|
|
if(ALightanet_obs(n,i).gt.ALightanet_obs(n-1,i).and.
|
|
&ALightanet_obs(n-1,i).gt.ALightanet_obs(n-2,i).and.
|
|
&ALightanet_obs(n-2,i).gt.ALightanet_obs(n-3,i))then
|
|
if(ALightpco2i(n,i).le.ALightpco2i(n-1,i).and.
|
|
&ALightpco2i(n-1,i).le.ALightpco2i(n-2,i).and.
|
|
&ALightpco2i(n-2,i).le.ALightpco2i(n-3,i))then
|
|
nstartalight(i)=n-1
|
|
nendalight(i)=n
|
|
!only the last point can be possibly Rubico or TPU because anet continues to rise while Ci is constant or decreasing
|
|
endif
|
|
endif
|
|
endif
|
|
if((nendalight(i)-nstartalight(i)).le.2)goto 64
|
|
62 n=nstartalight(i)
|
|
if(ALightanet_obs(n+1,i).gt.ALightanet_obs(n,i).and.
|
|
&ALightpco2i(n+1,i).le.ALightpco2i(n,i))then
|
|
!continue until we reach the point when anet does not increase while pco2i does not decrease, i,e, if anet continues
|
|
!to increase while pco2i continues to decrease, we assumue this point is still limited by rubp regeneration.
|
|
if((nendalight(i)-nstartalight(i)).gt.2)then
|
|
nstartalight(i)=n+1
|
|
goto 62
|
|
endif
|
|
else
|
|
nstartalight(i)=nstartalight(i)-1
|
|
nstartalight(i)=max0(nstartalight(i),0)
|
|
endif
|
|
64 continue
|
|
!we generally assume the points of an A/Light curve are lined up in a sequence of (RuBP, TPU and Rubisco), which is indicated by
|
|
!ialightorder=0. However, if Ci increases steadily from nstartalight to the end, we assume a sequence of (RuBP, Rubisco and TPU),
|
|
!which is indicated by ialightorder=2.
|
|
ialightorder(i)=2
|
|
do j=nstartalight(i)+1,nALightPoints(i)
|
|
if(ALightpco2i(j,i).lt.ALightpco2i(j-1,i))ialightorder(i)=0
|
|
enddo
|
|
enddo
|
|
!------------------------------------------------------------------------------------
|
|
!Merge Fixed points, ACi points, ALight points, and Free points into single arrays. Do not change this order.
|
|
ntotsamples=0
|
|
do i=1,nFixedPoints
|
|
ntotsamples=ntotsamples+1
|
|
anet_obs(ntotsamples)=Fixedanet_obs(i)
|
|
pco2i(ntotsamples)=Fixedpco2i(i)
|
|
templeaf(ntotsamples)=Fixedtempleaf(i)
|
|
aPPFDlf(ntotsamples)=FixedaPPFDlf(i)
|
|
pres_air(ntotsamples)=Fixedpres_air(i)
|
|
po2i(ntotsamples)=Fixedpo2i(i)
|
|
chlflphips2(ntotsamples)=Fixedchlflphips2(i)
|
|
pco2ambient(ntotsamples)=Fixedpco2ambient(i)
|
|
trmmol(ntotsamples)=Fixedtrmmol(i)
|
|
gswmeas(ntotsamples)=Fixedgswmeas(i)
|
|
vpdl(ntotsamples)=Fixedvpdl(i)
|
|
tempair(ntotsamples)=Fixedtempair(i)
|
|
eambient(ntotsamples)=Fixedeambient(i)
|
|
!
|
|
fo_pam(ntotsamples)=Fixedfo_pam(i)
|
|
fm_pam(ntotsamples)=Fixedfm_pam(i)
|
|
fs_pam(ntotsamples)=Fixedfs_pam(i)
|
|
pam_measlight(ntotsamples)=Fixedpam_measlight(i)
|
|
yield_ps2(ntotsamples)=Fixedyield_ps2(i)
|
|
yield_npq(ntotsamples)=Fixedyield_npq(i)
|
|
qlake(ntotsamples)=Fixedqlake(i)
|
|
qpuddle(ntotsamples)=Fixedqpuddle(i)
|
|
kps2_norm(ntotsamples)=Fixedkps2_norm(i)
|
|
knpq_norm(ntotsamples)=Fixedknpq_norm(i)
|
|
!
|
|
resiststomco2(ntotsamples)=Fixedresiststomco2(i)
|
|
enddo
|
|
do i=1,numACicurves
|
|
ACiavetempleaf(i)=0.0d0
|
|
ACiaveaPPFDlf(i)=0.0d0
|
|
ACiavepo2i(i)=0.0d0
|
|
do j=1,nACiPoints(i)
|
|
ntotsamples=ntotsamples+1
|
|
anet_obs(ntotsamples)=ACianet_obs(j,i)
|
|
pco2i(ntotsamples)=ACipco2i(j,i)
|
|
templeaf(ntotsamples)=ACitempleaf(j,i)
|
|
aPPFDlf(ntotsamples)=ACiaPPFDlf(j,i)
|
|
pres_air(ntotsamples)=ACipres_air(j,i)
|
|
po2i(ntotsamples)=ACipo2i(j,i)
|
|
chlflphips2(ntotsamples)=ACichlflphips2(j,i)
|
|
pco2ambient(ntotsamples)=ACipco2ambient(j,i)
|
|
trmmol(ntotsamples)=ACitrmmol(j,i)
|
|
gswmeas(ntotsamples)=ACigswmeas(j,i)
|
|
vpdl(ntotsamples)=ACivpdl(j,i)
|
|
tempair(ntotsamples)=ACitempair(j,i)
|
|
eambient(ntotsamples)=ACieambient(j,i)
|
|
!
|
|
fo_pam(ntotsamples)=ACifo_pam(j,i)
|
|
fm_pam(ntotsamples)=ACifm_pam(j,i)
|
|
fs_pam(ntotsamples)=ACifs_pam(j,i)
|
|
pam_measlight(ntotsamples)=ACipam_measlight(j,i)
|
|
yield_ps2(ntotsamples)=ACiyield_ps2(j,i)
|
|
yield_npq(ntotsamples)=ACiyield_npq(j,i)
|
|
qlake(ntotsamples)=ACiqlake(j,i)
|
|
qpuddle(ntotsamples)=ACiqpuddle(j,i)
|
|
kps2_norm(ntotsamples)=ACikps2_norm(j,i)
|
|
knpq_norm(ntotsamples)=ACiknpq_norm(j,i)
|
|
!
|
|
resiststomco2(ntotsamples)=ACiresiststomco2(j,i)
|
|
ACiavetempleaf(i)=ACiavetempleaf(i)+ACitempleaf(j,i)
|
|
ACiaveaPPFDlf(i)=ACiaveaPPFDlf(i)+ACiaPPFDlf(j,i)
|
|
ACiavepo2i(i)=ACiavepo2i(i)+ACipo2i(j,i)
|
|
enddo
|
|
ACiavetempleaf(i)=ACiavetempleaf(i)/dble(nACiPoints(i))
|
|
ACiaveaPPFDlf(i)=ACiaveaPPFDlf(i)/dble(nACiPoints(i))
|
|
ACiavepo2i(i)=ACiavepo2i(i)/dble(nACiPoints(i))
|
|
enddo
|
|
do i=1,numALightcurves
|
|
ALightavetempleaf(i)=0.0d0
|
|
ALightaveCO2ambient(i)=0.0d0
|
|
ALightavepo2i(i)=0.0d0
|
|
do j=1,nALightPoints(i)
|
|
ntotsamples=ntotsamples+1
|
|
anet_obs(ntotsamples)=ALightanet_obs(j,i)
|
|
pco2i(ntotsamples)=ALightpco2i(j,i)
|
|
templeaf(ntotsamples)=ALighttempleaf(j,i)
|
|
aPPFDlf(ntotsamples)=ALightaPPFDlf(j,i)
|
|
pres_air(ntotsamples)=ALightpres_air(j,i)
|
|
po2i(ntotsamples)=ALightpo2i(j,i)
|
|
chlflphips2(ntotsamples)=ALightchlflphips2(j,i)
|
|
pco2ambient(ntotsamples)=ALightpco2ambient(j,i)
|
|
trmmol(ntotsamples)=ALighttrmmol(j,i)
|
|
gswmeas(ntotsamples)=ALightgswmeas(j,i)
|
|
vpdl(ntotsamples)=ALightvpdl(j,i)
|
|
tempair(ntotsamples)=ALighttempair(j,i)
|
|
eambient(ntotsamples)=ALighteambient(j,i)
|
|
!
|
|
fo_pam(ntotsamples)=ALightfo_pam(j,i)
|
|
fm_pam(ntotsamples)=ALightfm_pam(j,i)
|
|
fs_pam(ntotsamples)=ALightfs_pam(j,i)
|
|
pam_measlight(ntotsamples)=ALightpam_measlight(j,i)
|
|
yield_ps2(ntotsamples)=ALightyield_ps2(j,i)
|
|
yield_npq(ntotsamples)=ALightyield_npq(j,i)
|
|
qlake(ntotsamples)=ALightqlake(j,i)
|
|
qpuddle(ntotsamples)=ALightqpuddle(j,i)
|
|
kps2_norm(ntotsamples)=ALightkps2_norm(j,i)
|
|
knpq_norm(ntotsamples)=ALightknpq_norm(j,i)
|
|
!
|
|
resiststomco2(ntotsamples)=ALightresiststomco2(j,i)
|
|
ALightavetempleaf(i)=ALightavetempleaf(i)+ALighttempleaf(j,i)
|
|
ALightaveCO2ambient(i)=ALightaveCO2ambient(i)+
|
|
&ALightpco2ambient(j,i)
|
|
ALightavepo2i(i)=ALightavepo2i(i)+ALightpo2i(j,i)
|
|
enddo
|
|
ALightavetempleaf(i)=ALightavetempleaf(i)/dble(nALightPoints(i))
|
|
ALightaveCO2ambient(i)=ALightaveCO2ambient(i)/
|
|
&dble(nALightPoints(i))
|
|
ALightavepo2i(i)=ALightavepo2i(i)/dble(nALightPoints(i))
|
|
enddo
|
|
do i=1,nFreePoints
|
|
ntotsamples=ntotsamples+1
|
|
anet_obs(ntotsamples)=Freeanet_obs(i)
|
|
pco2i(ntotsamples)=Freepco2i(i)
|
|
templeaf(ntotsamples)=Freetempleaf(i)
|
|
aPPFDlf(ntotsamples)=FreeaPPFDlf(i)
|
|
pres_air(ntotsamples)=Freepres_air(i)
|
|
po2i(ntotsamples)=Freepo2i(i)
|
|
chlflphips2(ntotsamples)=Freechlflphips2(i)
|
|
pco2ambient(ntotsamples)=Freepco2ambient(i)
|
|
trmmol(ntotsamples)=Freetrmmol(i)
|
|
gswmeas(ntotsamples)=Freegswmeas(i)
|
|
vpdl(ntotsamples)=Freevpdl(i)
|
|
tempair(ntotsamples)=Freetempair(i)
|
|
eambient(ntotsamples)=Freeeambient(i)
|
|
!
|
|
fo_pam(ntotsamples)=Freefo_pam(i)
|
|
fm_pam(ntotsamples)=Freefm_pam(i)
|
|
fs_pam(ntotsamples)=Freefs_pam(i)
|
|
pam_measlight(ntotsamples)=Freepam_measlight(i)
|
|
yield_ps2(ntotsamples)=Freeyield_ps2(i)
|
|
yield_npq(ntotsamples)=Freeyield_npq(i)
|
|
qlake(ntotsamples)=Freeqlake(i)
|
|
qpuddle(ntotsamples)=Freeqpuddle(i)
|
|
kps2_norm(ntotsamples)=Freekps2_norm(i)
|
|
knpq_norm(ntotsamples)=Freeknpq_norm(i)
|
|
!
|
|
resiststomco2(ntotsamples)=Freeresiststomco2(i)
|
|
enddo
|
|
ntotphips2=0
|
|
term1=1.0d+99
|
|
term2=-1.0d+99
|
|
do i=1,ntotsamples
|
|
pco2i_ori(i)=pco2i(i)
|
|
templeaf_ori(i)=templeaf(i)
|
|
if(templeaf(i).lt.term1)term1=templeaf(i)
|
|
if(templeaf(i).gt.term2)term2=templeaf(i)
|
|
aPPFDlf_ori(i)=aPPFDlf(i)
|
|
pres_air_ori(i)=pres_air(i)
|
|
po2i_ori(i)=po2i(i)
|
|
chlflphips2_ori(i)=chlflphips2(i)
|
|
pco2ambient_ori(i)=pco2ambient(i)
|
|
trmmol_ori(i)=trmmol(i)
|
|
gswmeas_ori(i)=gswmeas(i)
|
|
vpdl_ori(i)=vpdl(i)
|
|
tempair_ori(i)=tempair(i)
|
|
eambient_ori(i)=eambient(i)
|
|
resiststomco2_ori(i)=resiststomco2(i)
|
|
if(chlflphips2_ori(i).gt.0.0d0)then
|
|
ntotphips2=ntotphips2+1
|
|
endif
|
|
enddo
|
|
idotempcoeff=0
|
|
if((term2-term1).gt.5.0d0)idotempcoeff=1
|
|
!If temperature variation in the dataset is larger enough, try to estimate parameters in temperature response functions
|
|
!All variables are now in the right order. All ACi curves are ordered and All ALight curves are ordered.
|
|
!-------------------------------------------------------------------------------------------------------
|
|
!
|
|
do i=1,ntotsamples
|
|
anet_pred(i)=-9999.0d0
|
|
pco2i_pred(i)=-9999.0d0
|
|
pco2c(i)=-9999.0d0
|
|
anet_pred_flu(i)=-9999.0d0
|
|
pco2i_pred_flu(i)=-9999.0d0
|
|
pco2c_anet_flu(i)=-9999.0d0
|
|
pco2c_pco2i_flu(i)=-9999.0d0
|
|
enddo
|
|
if(ntotphips2.gt.5)then
|
|
do idorch=0,0
|
|
!we do a fluorescence only fit
|
|
Prioriknowlimit=-1
|
|
ifitmode=1
|
|
!ifitmode: =-2, ordinary fitting with pco2i calculated as a function of anet
|
|
!ifitmode: =-1, ordinary fitting with anet calculated as a function of pco2i
|
|
!ifitmode: =1, orthogonal fitting with anet calculated as a function of pco2i
|
|
!ifitmode: =2, orthogonal fitting with pco2i calculated as a function of anet
|
|
idorwp=0
|
|
resistwp25_ori=resistwp25_ini
|
|
if(idorch.eq.1)then
|
|
resistch25_ori=resistch25_ini
|
|
else
|
|
resistch25_ori=0.0d0
|
|
endif
|
|
if(rdlight25_usr.le.0.0d0)then
|
|
idord=1
|
|
rdlight25_ori=rdlight25_ini
|
|
else
|
|
idord=0
|
|
rdlight25_ori=rdlight25_usr
|
|
endif
|
|
idostargamma=1
|
|
idobetaPSII=1
|
|
idoha_darkresp=idotempcoeff
|
|
idoha_stargamma=idotempcoeff
|
|
idoha_gmeso=idotempcoeff
|
|
stargamma25_ori=stargamma25_ini
|
|
betaPSII_ori=betaPSII_ini
|
|
fjmax25_ori=fjmax25_ini
|
|
phifactor_ori=phifactor_ini
|
|
thetafactor_ori=thetafactor_ini
|
|
ha_darkresp_ori=ha_darkresp_ini
|
|
ha_stargamma_ori=ha_stargamma_ini
|
|
ha_gmeso_ori=ha_gmeso_ini
|
|
ha_jmax_ori=ha_jmax_ini
|
|
call HybridCombinatorial()
|
|
do j=1,ntotsamples
|
|
call gmesoontemp(templeaf(j),1.0d0,gascon,ha_gmeso,
|
|
&hd_gmeso,sv_gmeso,term)
|
|
resistwp=resistwp25/term
|
|
resistch=resistch25/term
|
|
call resp_mitocho(templeaf(j),rdlight25,ha_darkresp,
|
|
&gascon,rdlight)
|
|
call co2compens(templeaf(j),stargamma25,
|
|
&ha_stargamma,gascon,stargamma)
|
|
write(fluorescenceunit,370)trim(curvename),pco2i_ori(j),
|
|
&pco2i_pred(j),pco2c(j),anet_obs(j),anet_pred(j),
|
|
&pco2ambient_ori(j),po2i_ori(j)/1000.0d0,eambient_ori(j)/1000.0d0,
|
|
&pres_air_ori(j)/1000.0d0,vpdl_ori(j)/1000.0d0,
|
|
&aPPFDlf(j)/abspt_lf_par,templeaf_ori(j)-273.15d0,
|
|
&tempair_ori(j)-273.15d0,trmmol_ori(j),gswmeas_ori(j),
|
|
&chlflphips2_ori(j),rdlight25,resistwp25,resistch25,
|
|
&stargamma25,betaPSII,sumsquare,ha_darkresp,resistwp,resistch,
|
|
&ha_stargamma,fo_pam(j),fm_pam(j),fs_pam(j),pam_measlight(j),
|
|
&yield_ps2(j),yield_npq(j),qlake(j),qpuddle(j),kps2_norm(j),
|
|
&knpq_norm(j)
|
|
enddo
|
|
if(idorch.eq.0)then
|
|
fvector(1)=rdlight25
|
|
fvector(2)=resistwp25
|
|
fvector(3)=stargamma25
|
|
fvector(4)=betaPSII
|
|
endif
|
|
enddo
|
|
if(ntotlights.gt.0)then
|
|
!Jmax estimation with fluorescence data.
|
|
!Only points before nstartalight are used because these points are apparently limited by RuBP regeneration and therefore
|
|
!the electron transport equation applies. ntotlights is the number of points that are clearly limited by RuBP regeneration.
|
|
modeltype(0)='PARi'
|
|
modeltype(1)='TempLeaf'
|
|
modeltype(2)='PhiPSII_obs'
|
|
modeltype(3)='PhiPSII_pred'
|
|
modeltype(4)='Jmax25'
|
|
modeltype(5)='phifactor'
|
|
modeltype(6)='thetafactor'
|
|
modeltype(7)='ha_jmax'
|
|
modeltype(8)='SumSquare'
|
|
write(fluorescenceunit,305)(trim(modeltype(j)),j=0,8)
|
|
do j=1,ntotlights
|
|
write(fluorescenceunit,306)aparlights(j)/abspt_lf_par,
|
|
&templflights(j)-273.15d0,flphips2lights(j),PhiPSIIlights_pred(j),
|
|
&fjmax25,phifactor,thetafactor,ha_jmax,flujmaxfval
|
|
enddo
|
|
else
|
|
fjmax25=-9999.0d0
|
|
phifactor=-9999.0d0
|
|
thetafactor=-9999.0d0
|
|
flujmaxfval=-9999.0d0
|
|
endif
|
|
term=tempK_dark-273.15d0
|
|
if(term.lt.-10000.0d0)term=-9999.0d0
|
|
write(fluoresparamunit,380)trim(curvename),fjmax25,rdlight25,
|
|
&fvector(1),resistwp25,fvector(2),resistch25,stargamma25,
|
|
&fvector(3),phifactor,thetafactor,betaPSII,fvector(4),fo_dark,
|
|
&fm_dark,resp_dark,term,sumsquare,flujmaxfval
|
|
endif
|
|
!----------------------------------------------------------------
|
|
idophifactor=0
|
|
idothetafactor=0
|
|
idobetaPSII=0
|
|
ifitmode=-1
|
|
!ifitmode: =-2, ordinary fitting with pco2i calculated as a function of anet
|
|
!ifitmode: =-1, ordinary fitting with anet calculated as a function of pco2i
|
|
!ifitmode: =1, orthogonal fitting with anet calculated as a function of pco2i
|
|
!ifitmode: =2, orthogonal fitting with pco2i calculated as a function of anet
|
|
if(numALightcurves.ge.1)then
|
|
idophifactor=1
|
|
idothetafactor=1
|
|
endif
|
|
if(ntotphips2.ge.2)idobetaPSII=1
|
|
modeltype(0)='NoSuitModel'
|
|
modeltype(1)='RubiscoRuBpTpu'
|
|
modeltype(2)='RubiscoRuBp'
|
|
modeltype(3)='RubiscoTpu'
|
|
modeltype(4)='RuBpTpu'
|
|
modeltype(5)='Rubisco'
|
|
modeltype(6)='RuBp'
|
|
modeltype(7)='Tpu'
|
|
fourchars(1)='CO2i'
|
|
fourchars(2)='CO2cc'
|
|
fourchars(3)='Ac'
|
|
fourchars(4)='CO2cj'
|
|
fourchars(5)='Aj'
|
|
fourchars(6)='CO2ct'
|
|
fourchars(7)='At'
|
|
do k=1,4
|
|
do j=1,ntotsamples
|
|
co2c_Pa(k,j)=-9999.0d0
|
|
recycleratio(k,j)=-9999.0d0
|
|
recycleratio(5,j)=-9999.0d0
|
|
recycleratio(6,j)=-9999.0d0
|
|
enddo
|
|
enddo
|
|
do idorwp0=0,1
|
|
!When resistwp is estimated, we either estimate or not estimate resistch. But when resistwp is not estimated, neither is resistch.
|
|
if(idorwp0.eq.1)then
|
|
irchoption1=0
|
|
irchoption2=1
|
|
else
|
|
irchoption1=0
|
|
irchoption2=0
|
|
endif
|
|
do idorch0=irchoption1,irchoption2
|
|
do i2ndary=0,1
|
|
if(i2ndary.eq.0)then
|
|
idostargamma=0
|
|
idokc=0
|
|
idoko=0
|
|
idoalpha=0
|
|
idoha_stargamma=0
|
|
idoha_vcmax=0
|
|
idoha_jmax=0
|
|
idoha_tpu=0
|
|
else
|
|
idostargamma=1
|
|
idokc=1
|
|
!when data are sufficient (e.g. when multiple levels of oxygen are used in the measurements, set idoko=1
|
|
idoko=0
|
|
idoalpha=idoalpha0
|
|
idoha_stargamma=idotempcoeff
|
|
idoha_vcmax=idotempcoeff
|
|
idoha_jmax=idotempcoeff
|
|
idoha_tpu=idotempcoeff
|
|
endif
|
|
if(idorwp0.eq.0)then
|
|
idorwp=0
|
|
resistwp25_ori=0.0d0
|
|
idoha_gmeso=0
|
|
else
|
|
idoha_gmeso=idotempcoeff
|
|
if(resistwp25_usr.ge.0.0d0)then
|
|
!User provides a valid estimate of resistwp25 so don't estimate it even when idorwp0=1
|
|
idorwp=0
|
|
resistwp25_ori=resistwp25_usr
|
|
else
|
|
idorwp=1
|
|
resistwp25_ori=resistwp25_ini
|
|
endif
|
|
endif
|
|
if(idorch0.eq.0)then
|
|
idorch=0
|
|
resistch25_ori=0.0d0
|
|
idoha_gmeso=0
|
|
else
|
|
idoha_gmeso=idotempcoeff
|
|
if(resistch25_usr.ge.0.0d0)then
|
|
!User provides a valid estimate of resistch25 so don't estimate it even when idorch0=1
|
|
idorch=0
|
|
resistch25_ori=resistch25_usr
|
|
else
|
|
idorch=1
|
|
resistch25_ori=resistch25_ini
|
|
endif
|
|
endif
|
|
!rd has to be provided by user or to be estimated under any circumstances
|
|
idoha_darkresp=idotempcoeff
|
|
if(rdlight25_usr.gt.0.0d0)then
|
|
idord=0
|
|
rdlight25_ori=rdlight25_usr
|
|
else
|
|
idord=1
|
|
rdlight25_ori=rdlight25_ini
|
|
endif
|
|
!if kc25 is provided by the user, it is not estimated under any circumstances. Otherwise it is estimated or not estimated depending on i2ndary
|
|
if(fkc25_usr.gt.0.0d0)then
|
|
idokc=0
|
|
fkc25_ori=fkc25_usr
|
|
else
|
|
fkc25_ori=fkc25_ini
|
|
endif
|
|
!as long as oxygen is constant, it is hard to estimate ko25
|
|
if(fko25_usr.gt.0.0d0)then
|
|
idoko=0
|
|
fko25_ori=fko25_usr
|
|
else
|
|
fko25_ori=fko25_ini
|
|
endif
|
|
!if stargamma25 is provided by the user, it is not changed under any circumstances. Otherwise it is estimated or not estimated depending on i2ndary
|
|
if(stargamma25_usr.gt.0.0d0)then
|
|
idostargamma=0
|
|
stargamma25_ori=stargamma25_usr
|
|
else
|
|
stargamma25_ori=stargamma25_ini
|
|
endif
|
|
!if alpha25 is provided by the user, it is not changed under any circumstances. Otherwise it is to be estimated or not to be
|
|
!estimated depending on i2ndary; when it is not to be estimated, it is set to zero; when it is be to be estimated, it is initialized to a positive
|
|
!value because it is not good to initialize an unknown to be zero.
|
|
if(alpha25_usr.ge.0.0d0)then
|
|
idoalpha=0
|
|
alpha25_ori=alpha25_usr
|
|
else
|
|
if(idoalpha.eq.0)then
|
|
alpha25_ori=0.0d0
|
|
else
|
|
alpha25_ori=alpha25_ini
|
|
endif
|
|
endif
|
|
vcmax25_ori=vcmax25_ini
|
|
fjmax25_ori=fjmax25_ini
|
|
tpu25_ori=tpu25_ini
|
|
phifactor_ori=phifactor_ini
|
|
thetafactor_ori=thetafactor_ini
|
|
betaPSII_ori=betaPSII_ini
|
|
ha_darkresp_ori=ha_darkresp_ini
|
|
ha_stargamma_ori=ha_stargamma_ini
|
|
ha_vcmax_ori=ha_vcmax_ini
|
|
ha_jmax_ori=ha_jmax_ini
|
|
ha_tpu_ori=ha_tpu_ini
|
|
ha_gmeso_ori=ha_gmeso_ini
|
|
Prioriknowlimit=-9999
|
|
if((nFixedPoints+numACicurves+nFreePoints).eq.0)then
|
|
!If only light response curves are available, we don't estimate Kc and Ko because variations in Ci in light response curves are limited.
|
|
idokc=0
|
|
idoko=0
|
|
endif
|
|
!-------------------------------------------------------------------
|
|
call HybridCombinatorial()
|
|
!-------------------------------------------------------------------
|
|
if(bestnumtpu.eq.0)then
|
|
!use the asymptote of the A/Ci curve whose PAR is the highest to estimate tpu25
|
|
j=1
|
|
do i=2,numACicurves
|
|
if(ACiaveaPPFDlf(i).gt.ACiaveaPPFDlf(j))j=i
|
|
enddo
|
|
call resp_mitocho(ACiavetempleaf(j),rdlight25,
|
|
&ha_darkresp,gascon,rdlight)
|
|
tpu25=(Amax_ACi(j)+rdlight)/3.0d0
|
|
call tpuontemp(ACiavetempleaf(j),gascon,1.0d0,ha_tpu,
|
|
&hd_tpu,sv_tpu,term)
|
|
tpu25=tpu25/term
|
|
alpha25=0.0d0
|
|
idoalpha=0
|
|
endif
|
|
if(bestnumrubis.eq.0)then
|
|
idokc=0
|
|
idoko=0
|
|
endif
|
|
!Calculation of CO2 recycling ratio
|
|
if(idorwp0.eq.1)then
|
|
if(idorch0.eq.0.and.i2ndary.eq.0)k=1
|
|
if(idorch0.eq.0.and.i2ndary.eq.1)k=2
|
|
if(idorch0.eq.1.and.i2ndary.eq.0)k=3
|
|
if(idorch0.eq.1.and.i2ndary.eq.1)k=4
|
|
else
|
|
if(i2ndary.eq.0)k=5
|
|
if(i2ndary.eq.1)k=6
|
|
endif
|
|
stargamma25fit(k)=stargamma25
|
|
idomeso=(idorwp+1)*100+(idorch+1)*10+idoha_gmeso+1
|
|
j=0
|
|
if(bestilimittype.le.3.or.bestilimittype.eq.5)j=idoha_vcmax
|
|
idohavjt=(j+1)*100
|
|
j=0
|
|
if(bestilimittype.le.2.or.bestilimittype.eq.4.or.
|
|
&bestilimittype.eq.6)j=idoha_jmax
|
|
idohavjt=idohavjt+(j+1)*10
|
|
j=0
|
|
if(bestilimittype.eq.1.or.bestilimittype.eq.3.or.
|
|
&bestilimittype.eq.4.or.bestilimittype.eq.7)j=idoha_tpu
|
|
idohavjt=idohavjt+j+1
|
|
! = 1, Rubisco+RuBP+TPU
|
|
! = 2, Rubisco+RuBP
|
|
! = 3, Rubisco+TPU
|
|
! = 4, RuBP+TPU
|
|
! = 5, Rubisco Only
|
|
! = 6, RuBP Only
|
|
! = 7, TPU Only
|
|
idostargamma=(idostargamma+1)*10+idoha_stargamma+1
|
|
idokc=(idokc+1)*10+1
|
|
idoko=(idoko+1)*10+1
|
|
idord=(idord+1)*10+(idoha_darkresp+1)
|
|
idoalpha=(idoalpha+1)
|
|
do j=1,ntotsamples
|
|
if(k.le.4)co2c_Pa(k,j)=pco2c(j)
|
|
if(resiststomco2(j).ge.0.0d0)then
|
|
call gmesoontemp(templeaf(j),1.0d0,gascon,ha_gmeso,
|
|
&hd_gmeso,sv_gmeso,term)
|
|
resistwp=resistwp25/term
|
|
resistch=resistch25/term
|
|
call resp_mitocho(templeaf(j),rdlight25,ha_darkresp,
|
|
&gascon,rdlight)
|
|
call co2compens(templeaf(j),stargamma25,
|
|
&ha_stargamma,gascon,stargamma)
|
|
call co2recyclingratio(resistwp,resistch,
|
|
! &resiststomco2(j),stargamma,rdlight,pco2c(j),anet_obs(j),
|
|
&resiststomco2(j),stargamma,rdlight,pco2c(j),anet_pred(j),
|
|
&recycleratio(k:k,j:j))
|
|
else
|
|
recycleratio(k,j)=-9999.0d0
|
|
endif
|
|
write(compareunit,300)trim(curvename),idomeso,idohavjt,
|
|
&idostargamma,idokc,idoko,idord,idoalpha,idobetaPSII+1,
|
|
&pco2i_ori(j),pco2i_pred(j),pco2c(j),anet_obs(j),anet_pred(j),
|
|
&bestiphotolimit(j),recycleratio(k,j)*100.0d0,pco2ambient_ori(j),
|
|
&po2i_ori(j)/1000.0d0,eambient_ori(j)/1000.0d0,
|
|
&pres_air_ori(j)/1000.0d0,vpdl_ori(j)/1000.0d0,
|
|
&aPPFDlf(j)/abspt_lf_par,templeaf_ori(j)-273.15d0,
|
|
&tempair_ori(j)-273.15d0,trmmol_ori(j),gswmeas_ori(j),
|
|
&chlflphips2_ori(j),PhiPSII_pred(j),pco2i_pred_flu(j),
|
|
&anet_pred_flu(j),pco2c_pco2i_flu(j),pco2c_anet_flu(j)
|
|
enddo
|
|
!Generate mono-limiting curves
|
|
k=nFixedPoints
|
|
do i=1,numACicurves
|
|
n=k+nACiPoints(i)
|
|
j=n-k
|
|
call ilimittypestats(j,bestiphotolimit(k+1:n),
|
|
&Currentilimittype,numrubis,numrubp,numtpu)
|
|
write(compareunit,310)(trim(fourchars(j)),j=1,7)
|
|
co2imany(1)=1.0d0
|
|
co2imany(2)=2.0d0
|
|
co2imany(3)=3.0d0
|
|
co2imany(4)=4.0d0
|
|
co2imany(5)=5.0d0
|
|
co2imany(6)=6.0d0
|
|
m=6
|
|
term=ACipco2i(nACiPoints(i),i)+10.0d0
|
|
do while (co2imany(m).le.term)
|
|
m=m+1
|
|
co2imany(m)=co2imany(m-1)+2.5d0
|
|
enddo
|
|
do j=1,m
|
|
ccc=co2imany(j)
|
|
ccj=co2imany(j)
|
|
cct=co2imany(j)
|
|
if(numrubis.gt.0)then
|
|
Currentilimittype=5
|
|
call leafanetmodel(Currentilimittype,ACiaveaPPFDlf(i),
|
|
&ACiavetempleaf(i),co2imany(j),ACiavepo2i(i),-9999.0d0,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,alpha25,
|
|
&rdlight25,ha_darkresp,stargamma25,ha_stargamma,fkc25,
|
|
&ha_kc,fko25,ha_ko,phifactor,thetafactor,betaPSII,
|
|
&ac,rdlight,m,atp,resistwp,resistch,stargamma,ccc,realizedfjelect,
|
|
&term1,term2)
|
|
else
|
|
ac=-9999.0d0
|
|
ccc=-9999.0d0
|
|
endif
|
|
if(numrubp.gt.0)then
|
|
Currentilimittype=6
|
|
call leafanetmodel(Currentilimittype,ACiaveaPPFDlf(i),
|
|
&ACiavetempleaf(i),co2imany(j),ACiavepo2i(i),-9999.0d0,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,alpha25,
|
|
&rdlight25,ha_darkresp,stargamma25,ha_stargamma,fkc25,
|
|
&ha_kc,fko25,ha_ko,phifactor,thetafactor,betaPSII,
|
|
&aj,rdlight,m,atp,resistwp,resistch,stargamma,ccj,realizedfjelect,
|
|
&term1,term2)
|
|
else
|
|
aj=-9999.0d0
|
|
ccj=-9999.0d0
|
|
endif
|
|
Currentilimittype=7
|
|
call leafanetmodel(Currentilimittype,ACiaveaPPFDlf(i),
|
|
&ACiavetempleaf(i),co2imany(j),ACiavepo2i(i),-9999.0d0,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,alpha25,
|
|
&rdlight25,ha_darkresp,stargamma25,ha_stargamma,
|
|
&fkc25,ha_kc,fko25,ha_ko,phifactor,thetafactor,
|
|
&betaPSII,at,rdlight,m,atp,resistwp,resistch,stargamma,cct,
|
|
&realizedfjelect,term1,term2)
|
|
if(cct.le.((1.0d0+3.0d0*alpha25)*stargamma))then
|
|
cct=-9999.0d0
|
|
at=-9999.0d0
|
|
endif
|
|
write(compareunit,320)co2imany(j),ccc,ac,ccj,aj,cct,at
|
|
k=n
|
|
enddo
|
|
enddo
|
|
write(compareunit,*)
|
|
!-------------------------------------------------------------------------------
|
|
!Compute Rubisco-RuBP-TPU transitional points at 25oC and PAR = 1000 umolm-2s-1
|
|
j=bestilimittype
|
|
!tpu and alpha are always available
|
|
if(bestilimittype.eq.2)j=1
|
|
if(bestilimittype.eq.5)j=3
|
|
if(bestilimittype.eq.6)j=4
|
|
if(fjmax25.gt.0.0d0)then
|
|
call jontemp(1000.0d0*abspt_lf_par,298.15d0,term,fjmax25,
|
|
&ha_jmax,hd_jmax,sv_jmax,gascon,phifactor,thetafactor,betaPSII)
|
|
endif
|
|
call EqualPoints(vcmax25,term,tpu25,resistwp25,resistch25,
|
|
&stargamma25,fkc25,fko25,21000.0d0,alpha25,rdlight25,j,
|
|
&co2iRubismax25,co2iRuBpmax25,anetRubismax25,anetRuBpmax25)
|
|
!------------------------------------------------------------------------------
|
|
write(paramunit,330)trim(curvename),idomeso,idohavjt,
|
|
&idostargamma,idokc,idoko,idord,idoalpha,idobetaPSII+1,
|
|
&trim(modeltype(bestilimittype)),vcmax25,fjmax25,
|
|
&rdlight25,resistwp25,resistch25,tpu25,
|
|
&stargamma25,fkc25,fko25,alpha25,ha_vcmax,
|
|
&hd_vcmax,sv_vcmax,ha_jmax,hd_jmax,sv_jmax,ha_tpu,hd_tpu,
|
|
&sv_tpu,ha_gmeso,hd_gmeso,sv_gmeso,ha_darkresp,ha_stargamma,
|
|
&ha_kc,ha_ko,phifactor,thetafactor,betaPSII,
|
|
&bestnumrubis,bestnumrubp,bestnumtpu,ntotsamples,bestsumsquare,
|
|
&co2iRubismax25,co2iRuBpmax25,anetRubismax25,anetRuBpmax25,
|
|
&trim(siteID),Latitude,Longitude,Elevation,yearsampled,sampledoy,
|
|
&GrowingSeasonStart,GrowingSeasonEnd,standage,CanopyHeight,
|
|
&LeafAreaIndex,trim(species),avetimeresolution,avetimesampled,
|
|
&SampleHeight,Needleage,specificLAI,nitrogencontent,carboncontent,
|
|
&phoscontent,trim(woodporosity),sapwooddensity,leafratio
|
|
enddo
|
|
enddo
|
|
enddo
|
|
if(numACicurves.gt.0)then
|
|
do i=1,numACicurves
|
|
write(aciempfitunit,390)trim(curvename),i,starco2i(i),
|
|
&der_starco2i(i),Amax_ACi(i),ACiinter(i),der_ACiinter(i),
|
|
&der_ACiend(i),PhiPSIImax_ACi(i),PhiPSIIinter_ACi(i),
|
|
&der_PhiPSIIinter_ACi(i),der_PhiPSIIend_ACi(i),
|
|
&ACimaxcurvature(i),ACimaxcurvpco2i(i),
|
|
&PhiPSIImaxcurvature_ACi(i),PhiPSIImaxcurv_ACi(i),
|
|
&starco2a(i),der_starco2a(i),Amax_ACa(i),ACainter(i),
|
|
&der_ACainter(i),der_ACa400ppm(i),anet_ACa400ppm(i),
|
|
&PhiPSIImax_ACa(i),PhiPSIIinter_ACa(i),der_PhiPSIIinter_ACa(i),
|
|
&der_PhiPSIIend_ACa(i),ACamaxcurvature(i),ACamaxcurvpco2a(i),
|
|
&PhiPSIImaxcurvature_ACa(i),PhiPSIImaxcurv_ACa(i),ACiavetempleaf(i)
|
|
&-273.15d0,ACiaveaPPFDlf(i)/abspt_lf_par,ACiavepo2i(i),
|
|
&trim(siteID),Latitude,Longitude,Elevation,yearsampled,sampledoy,
|
|
&GrowingSeasonStart,GrowingSeasonEnd,standage,CanopyHeight,
|
|
&LeafAreaIndex,trim(species),avetimeresolution,avetimesampled,
|
|
&SampleHeight,Needleage,specificLAI,nitrogencontent,carboncontent,
|
|
&phoscontent,trim(woodporosity),sapwooddensity,leafratio
|
|
enddo
|
|
endif
|
|
if(numALightcurves.gt.0)then
|
|
do i=1,numALightcurves
|
|
write(alightempfitunit,360)trim(curvename),i,starPAR(i),
|
|
&der_starPAR(i),Asat_ALight(i),ALightinter(i),der_ALightinter(i),
|
|
&der_ALightend(i),PhiPSIIinter_ALight(i),
|
|
&der_PhiPSIIinter_ALight(i),ExcessLightFactor(i),
|
|
&der_PhiPSII1000_ALight(i),ALightmaxcurvature(i),
|
|
&ALightmaxcurvPAR(i),PhiPSIImaxcurvature_ALight(i),
|
|
&PhiPSIImaxcurv_ALight(i),ALightavetempleaf(i)-273.15d0,
|
|
&ALightaveCO2ambient(i),ALightavepo2i(i),
|
|
&trim(siteID),Latitude,Longitude,Elevation,yearsampled,sampledoy,
|
|
&GrowingSeasonStart,GrowingSeasonEnd,standage,CanopyHeight,
|
|
&LeafAreaIndex,trim(species),avetimeresolution,avetimesampled,
|
|
&SampleHeight,Needleage,specificLAI,nitrogencontent,carboncontent,
|
|
&phoscontent,trim(woodporosity),sapwooddensity,leafratio
|
|
enddo
|
|
endif
|
|
call LeafGasFit_Stom(stomwuecicaoutunit,wuecicacompunit,
|
|
&stomcompunit,icurveno_usr,curvename,ntotsamples,aPPFDlf,templeaf,
|
|
&tempair,pco2i,pco2ambient,pres_air,anet_obs,gswmeas,vpdl,trmmol,
|
|
&abspt_lf_par,co2c_Pa(1:4,1:ntotsamples),
|
|
&recycleratio(1:6,1:ntotsamples),stargamma25fit,ha_stargamma,
|
|
&siteID,Latitude,Longitude,Elevation,yearsampled,sampledoy,
|
|
&GrowingSeasonStart,GrowingSeasonEnd,standage,CanopyHeight,
|
|
&LeafAreaIndex,species,avetimeresolution,avetimesampled,
|
|
&SampleHeight,Needleage,specificLAI,nitrogencontent,carboncontent,
|
|
&phoscontent,woodporosity,sapwooddensity,leafratio)
|
|
return
|
|
300 format(a,',',8(i0,','),5(f0.6,','),i0,',',16(f0.6,','),f0.6)
|
|
305 format(8(a,','),a)
|
|
306 format(8(f0.6,','),f0.6)
|
|
310 format(6(a,','),a)
|
|
320 format(6(f0.6,','),f0.6)
|
|
330 format(a,',',8(i0,','),a,',',29(f0.6,','),4(i0,','),
|
|
&5(f0.6,','),a,',',10(f0.6,','),a,',',8(f0.6,','),a,',',f0.6,
|
|
&',',f0.6)
|
|
360 format(a,',',i0,',',17(f0.6,','),a,',',10(f0.6,','),a,',',
|
|
&8(f0.6,','),a,',',f0.6,',',f0.6)
|
|
370 format(a,',',35(f0.6,','),f0.6)
|
|
380 format(a,',',17(f0.6,','),f0.6)
|
|
390 format(a,',',i0,',',32(f0.6,','),a,',',
|
|
&10(f0.6,','),a,',',8(f0.6,','),a,',',f0.6,',',f0.6)
|
|
end subroutine SetUpLeafGasFit
|