*XXX Check occurrence of g = grav. accel./ units of khs; is it why I * had to post-multiply by 0.1? * gsBBxABA.wRsol.engl.f - a program to calculate stomatal conductance * (and attendant rates of transpiration and CO2 assimilation) in * response to aerial environment of leaf plus the water-stress * response mediated by ABA. The aerial response is represented by * the Farquhar-von Caemerrer-Berry model of assimilation, the * Ball-Berry equation of stomatal control, and a standard equation * of leaf energy balance * A complete description of the model is in the file "explan.model" * (this is in "troff" format; "explan.model.pdf" is in PDF format). * These files are to be found in the same directory as this file * with Fortran code * At the end of this code is a list of all variables and arrays, * with their physical/biological meaning and their units * * Written 9 June 1999 by V. P. Gutschick while on sabbatical from * New Mexico State University, Las Cruces, NM, USA, visiting * F. C. Tardieu and T. Simonneau at the Laboratoire d'Ecophysiologie * des Plantes sur Stress Environnementaux, INRA ENSA-M, * Montpellier, France (forgive the orthographic omissions not * representable in flat ASCII) * Expanded 28 Nov. 1999 * Revised several times; most recently 20 Nov. 2001 * ___________________________________________________________ * | Vince Gutschick Internet: vince@nmsu.edu | * | Dept. of Biology, 3AF Phone 1-505-646-5661 | * | New Mexico State Univ. FAX 1-505-646-5665 | * | Las Cruces, NM 88003 U.S.A. | * | | * | Please, do not send duplicate text in HTML | * ----------------------------------------------------------- * Obtained (most likely) from Web site * http://biology-web.nmsu.edu/vince * * This version implements a very easy and general way to change * values of parameters, say, to do a run with different eair values * This version has been cleaned up, too (aPAR, aNIR set explicitly; * Pair allowed to be set; ...) * ...and at end, all info is printed out (Tleaf, E, A, hs, Ci) * * Variable type declarations (most are implicit, e.g., real unless * initial letter is i-n; full specifications are noted in the * explanation of all variables at end of code here real IBBhi,IBBlo,IBBmid real Kc25,Kcotop,Kctop,Ko25,Kotop real mBB logical BBswitch character*3 string(2) common/mainparam/hs0,u,dleaf,Tsky,Tveg,aPAR,mBB,bBB,tol,ha common/mainTleaf/epsleaf,sigma,Hvap,gtotw,Bcc,BTIR,Qrad,psir, + cABA,psiL,facABA,eL,E common/mainA/Kc25,Ko25,gtotCO2,gbair,Q00,TL,A,hs,Cs,Ci common/tous/Pair common/mainTleafparam/Tair,eair,psis,Rsol common/mainAparam/Tdeac,PPFD,Vcmax25,Ca data string/'BB ','ABA'/ open(8,file='scratch') * Set some constants aNIR=0.35 sigma=5.67e-8 epsleaf=0.96 Hvap=4.5e4 Cpair=29. Kc25=40.4 Ko25=24800. Q00=0.0732 ncha=0 istop=0 call paramread(ncha,istop) 50 BBswitch=0 if(hs0.lt.0.) BBswitch=1 * Calc. some derived quantities * Recompute gbair gbair=6.62e-3*sqrt(u/dleaf)*Pair/(8.314*(Tair+273.2)) Bcc=0.85*Cpair*gbair Qsky=sigma*(Tsky+273.2)**4 Qveg=sigma*(Tveg+273.2)**4 * These imply the leaf sees sky on one side, veg. on other, only QSW=PPFD*2.2e5*(aPAR+aNIR) * Converts from umol to mol and then mol to J; finally, it acct.s * for avg. NIR energy flux density, as equal to that in PAR Qrad=epsleaf*(Qsky+Qveg)+QSW BTIR=8.*epsleaf*sigma*(Tair+273.2)**3 * New: for the effect of ABA, calculate gsmax, with Ci=0.8*Ca, * IL=ILsat,TL=Tdeac (T of highest photosyn. capacity, just before * thermal deactivation sets it). Then, calculate the corresponding * gs from the Ball-Berry relation, m*Amax*(hs=1)/(Cs=0.8*Ca) + b * De nouveau: pour l'effet d'ABA: on calcule gsmax, a Ci=0.8*Ca, * IL=ILsat, TL=Tdeac (juste avant la deactivation; le plus haute Vcmax); * donc, on calcule m*Amax*0.8/Ca + b (on prends hs -> 0.8, pas 1): denom=298*8.314*(Tdeac+273.2) Vcmaxtop=Vcmax25*exp(64800.*(Tdeac-25.)/denom) gammatop=3.69+0.188*(Tdeac-25.)+0.00036*(Tdeac-25.)**2 Kctop=Kc25*exp(59400.*(Tdeac-25.)/denom) Kotop=Ko25*exp(36000.*(Tdeac-25.)/denom) Kcotop=Kctop*(1.+0.21*Pair/Kotop) Citop=0.8*Ca Atop=Vcmaxtop*(Citop-gammatop)/(Citop+Kcotop) gsadd=mBB*Atop*0.8*Pair/Ca gsmax=gsadd+bBB write(6,'(1x,"Calculated gsmax=",f7.4)')gsmax * Method: assume a gs (at two endpoints) * Calculate Tleaf, thus, E, (hs,) Vcmax, and Q00 * Solve for Ci in full Johnson-Thornley equation, then for A * Test if gs >, <, or = mA(*hs=1)/Cs + bBB * Added today (28 Nov. 1999): try gs - gsBB, as above, and also * gs - gs(ABA), taking the max. of the two; if it is >0, then * gs is too high; reduce it. * Ajoute' aujourd'hui: on essais gs - gs(BB), comme en haute, et * aussi gs - gs(ABA). On prends le maximum de ces quantites; * s'il est >0, le gs est trop haut, on doit le reduire - ca veut * dire, si gs est trop haut par l'un ou l'autre, on le reduit. * * Get a starting guess for gs: take Tleaf=min(Tair,Tdeac), solve for A; then * compute Cs and compose A*hs0/Cs to estimate gs; take 0.7 and * 1.5x this as endpoints * Why min(Tair, Tdeac)? Because the leaf may hit TL > Tair; an estimate * of TL=Tair makes Vcmax too low...and A, gs too low, in a feedack loop * Pourquoi min(Tair,Tdeac)? Parce que la feuille peut arriver a TL>Tair, * une estimation de TL=Tair fait Vcmax trop bas, et ca fait A, gs, ... * trop bas, dans un boucle de 'feedback' TL=amin0(Tair,Tdeac) * Need an estimate of gtotCO2; set it for infinite gs gtotCO2=0.72*gbair call CiACssolve if(BBswitch)then hs=ha else hs=hs0 endif * Use air humidity as guess of hs, if using full BB response TL=Tleaf(Tair) ! Call this so that facABA is computed !!! * !!! = changes to go from model of minimum to multiplicative model gs0=facABA*mBB*A*hs*Pair/Cs !!! gs0=amax0(gs0,bBB) gslo=0.7*gs0 gshi=1.5*gs0 * Do gslo first, in a generic form gtotCO2=1./(1.37/gbair+1.6/gslo) gtotw=1./(1./gbair+1./gslo) TL=Tleaf(gslo) * Solve Vcmax*(Ci-gamma)=(Ci+Kco)*gtotCO2*(Ca-Ci)/Pair for Ci, then * substitute to get A * Then solve Q00*(Ci-gamma)/(Ci+2*gamma) = gtotCO2*(Ca-Ci)/Pair for Ci, * then substitute Ci to get new A; choose lower of this, first A call CiACssolve if(BBswitch)then hs=(gbair*eair/eL+gslo)/(gbair+gslo) else hs=hs0 endif IBBlo=facABA*mBB*A*hs*Pair/Cs+bBB !!! testlo=gslo-IBBlo !!! (many lines deleted, for testloBB, testloABA * Now do gshi 95 gtotCO2=1./(1.37/gbair+1.6/gshi) gtotw=1./(1./gbair+1./gshi) TL=Tleaf(gshi) call CiACssolve if(BBswitch)then hs=(gbair*eair/eL+gshi)/(gbair+gshi) else hs=hs0 endif IBBhi=facABA*mBB*A*hs*Pair/Cs+bBB !!! testhi=gshi-IBBhi !!! * Enter decision loop * Test if gs > mBB*A/Cs * CUTS FROM ROOTSEARCH.F 97 if(testlo.lt.0.)then if(testhi.lt.0.)then * write(6,90) * 90 format(1x,'No root evident in interval') * Fix this problem: remake the search interval gslo=gshi testlo=testhi gshi=1.5*gslo go to 95 elseif(testhi.eq.0.d0)then write(6,100)gshi gs=gshi goto 200 100 format(1x,'Exact solution at gs=',f12.6) else gsplus=gshi gsminus=gslo endif elseif(testlo.eq.0.d0)then write(6,100)gslo gs=gslo goto 200 else if(testhi.lt.0.)then gsminus=gshi gsplus=gslo elseif(testhi.eq.0.)then write(6,100)gshi gs=gshi goto 200 else * write(6,90) * Fix this problem: remake the search interval gshi=gslo testhi=testlo gslo=0.6*gshi * Solve for testlo now; saves redoing calcs. from stmt. 95 gtotCO2=1./(1.37/gbair+1.6/gslo) gtotw=1./(1./gbair+1./gslo) TL=Tleaf(gslo) call CiACssolve if(BBswitch)then hs=(gbair*eair/eL+gslo)/(gbair+gslo) else hs=hs0 endif IBBlo=facABA*mBB*A*hs*Pair/Cs+bBB !!! testlo=gslo-IBBlo !!! go to 97 endif endif * Now run through binary search; 1st, figure out how many times nruns=alog((gshi-gslo)/tol)/0.693 do n=1,nruns gsmid=0.5*(gsplus+gsminus) gtotCO2=1./(1.37/gbair+1.6/gsmid) gtotw=1./(1./gbair+1./gsmid) TL=Tleaf(gsmid) call CiACssolve if(BBswitch)then hs=(gbair*eair/eL+gsmid)/(gbair+gsmid) else hs=hs0 endif IBBmid=facABA*mBB*A*hs*Pair/Cs+bBB !!! testmid=gsmid-IBBmid !!! if(testmid.lt.0.)then gsminus=gsmid elseif(testmid.eq.0.)then write(6,100)gsmid goto 200 else gsplus=gsmid endif enddo 200 write(6,120)gsminus,gsplus 120 format(1x,'Exact solution lies between',2f16.6) write(6,201)facABA 201 format(1x,"Limitation by ABA is f(ABA)=",f6.4) write(6,202)cABA,psir,psiL 202 format(1x,"[ABA]=",f8.2," nM; psi(root)=",f8.4," MPa; ", + " psi(leaf)=",f8.4," MPa") gs=gsplus write(6,800)TL,E,1.e6*A,hs,Ci,Cs 800 format(3x,'Tleaf=',f6.2,' E=',f8.5,' A=',f6.3,' hs=',f7.5, + ' Ci=',f6.2,' (Pa)',/,10x,'Cs=',f6.2) * Allow options of changing parameter values; pick up param. to be * changed by keyword ncha=1 call paramread(ncha,istop) if(istop.eq.1)stop goto 50 end * * function Tleaf(gs) common/mainTleaf/epsleaf,sigma,Hvap,gtotw,Bcc,BTIR,Qrad,psir, + cABA,psiL,facABA,eL,E common/mainTleafparam/Tair,eair,psis,Rsol common/Tleafparam/aABA,bABA,Rhplant,deltaABA,betaABA common/tous/Pair Tleaf=Tair 10 Qmrad=2.*epsleaf*sigma*(Tleaf+273.2)**4 eL=esat(Tleaf) E=gtotw*(eL-eair)/Pair QmE=Hvap*E Qmcc=Bcc*(Tleaf-Tair) Qnet=Qrad-Qmrad-QmE-Qmcc ep=esatpr(Tleaf) Bnet=BTIR+Hvap*gtotw*ep/Pair+Bcc dTleaf=Qnet/Bnet Tleaf=Tleaf+dTleaf * Now, calculate psiL, ... * Maintenant, on calcule psiL, [ABA], facABA Emasse=E*0.018 ! kg m-2 s-1 dpsiPa=Rsol*Emasse psir=psis-dpsiPa*1.e-6 ! en MPa cABA=-aABA*psir/(E+bABA) psiL=psir-E*Rhplant arg=exp(deltaABA*psiL) facABA=exp(betaABA*cABA*arg) if(abs(dTleaf).gt.0.01)go to 10 return end * * subroutine CiACssolve real Kc25,Kc,Ko25,Ko,Kco common/tous/Pair common/mainA/Kc25,Ko25,gtotCO2,gbair,Q00,TL,A,hs,Cs,Ci common/mainAparam/Tdeac,PPFD,Vcmax25,Ca common/Aparam/convex,Tslope dfac=298.*8.314 denom=dfac*(TL+273.2) Vcmax=Vcmax25*exp(64800.*(TL-25.)/denom) if(TL.gt.Tdeac)then amult=1.-Tslope*(TL-Tdeac) amult=amax0(amult,0.) Vcmax=Vcmax*amult endif Kc=Kc25*exp(59400.*(TL-25.)/denom) Ko=Ko25*exp(36000.*(TL-25)/denom) Kco=Kc*(1.+0.21*Pair/Ko) gamma=3.69+(TL-25.)*(0.188+0.00036*(TL-25.)) * New method: solve the quartic eq. for Ci; substitute Ci into the * eqs. for Amax, ALL, et solve the Johnson-Thornley eq. for A * (and then, for Cs) * Nouveau methode: re'soudre l'eq. quartique (mot?) pour Ci, * substituer Ci dans les eq. pour Amax, A(L-L), et resoudre * l'eq. Johnson-Thornley pour A (et puis, Cs) * Solve the quartic equation for Ci, from * convex*A**2 - A*(Amax+Q0*IL) - Q0*IL*Amax = 0, which expands as * convex*(gtotc*(Ca-Ci)/Pair)**2 -gtotc*[(Ca-Ci)/Pair]* * (Vcmax*(Ci-gamma)/(Ci+Kco) + Q00*IL*(Ci-gamma)/(Ci+gamma) ) * - Qoo*IL*Vcmax*[(Ci-gamma)/(Gi+gamma)][(Ci-gamma)/(Ci+Kco)] f1=convex*gtotCO2*gtotCO2/(Pair*Pair) f2=gtotCO2*Q00*PPFD/Pair f3=gtotCO2*Vcmax/Pair f4=Q00*PPFD*Vcmax t4=f1 t3=f1*(gamma+Kco-2.*Ca) + f2 + f3 t2=f1*(gamma*Kco+Ca*(Ca-2.*(gamma+Kco))) + f2*(Kco-Ca-gamma) + -f3*Ca + f4 t1=f1*Ca*(Ca*(gamma+Kco)-2.*gamma*Kco) + f2*(Ca*gamma- + Kco*(Ca+gamma)) -f3*gamma*gamma - f4*2.*gamma t0=f1*Ca*Ca*gamma*Kco + f2*Ca*gamma*Kco + f3*Ca*gamma*gamma + +f4*gamma*gamma * Start the iterations Ci=Ca iter=0 100 F=t0 + Ci*(t1 + Ci*(t2 + Ci*(t3 + Ci*t4))) Fprime=t1 + Ci*(2.*t2 + Ci*(3.*t3 + 4.*Ci*t4)) dCi=-F/Fprime Ci=Ci+dCi iter=iter+1 if(abs(dCi).gt.0.01)goto 100 * write(6,'("In ",i3," iterations we found Ci=",f6.3)') * + iter,Ci * Compute A and compare to Amax(Ci) and Q0(Ci)*IL Amax=Vcmax*(Ci-gamma)/(Ci+Kco) ALL=Q00*(Ci-gamma)*PPFD/(Ci+gamma) * Solve quadratic for A: A=(0.5/convex)*(Amax+ALL-sqrt((Amax+ALL)**2-4.*convex* + Amax*ALL)) * write(6,20)A,Amax,ALL * 20 format(1x,'A=',6Pf6.3,' umol m-2 s-1; Amax=',6Pf6.3, * + ', A(L-L)=',6Pf6.3) Cs=Ca-1.37*A*Pair/gbair return end * * FUNCTION ESAT(T) ESAT=610.78*EXP(17.269*T/(T+237.3)) * = SAT'D VAPOR P. OF WATER (PASCALS) AT AMBIENT TEMP. T (DEG. C) RETURN END * * FUNCTION ESATPR(T) ESATPR=610.78*EXP(17.269*T/(T+237.3))*4097.9 1 /(T+237.3)**2 * = (D/DT) PSAT(T) RETURN END * * subroutine paramread(ncha,istop) character*20 chaval, test*3, sval*17 real mr,La,Lv,mBB,khstar dimension thetasval(5),thetarval(5),alphasolval(5), + emsolval(5),ensolval(5),rhobval(5),c0val(5), + c1val(5),c2val(5),c3val(5),Fval(5) common/mainparam/hs0,u,dleaf,Tsky,Tveg,aPAR,mBB,bBB,tol,ha common/mainTleafparam/Tair,eair,psis,Rsol common/Tleafparam/aABA,bABA,Rhplant,deltaABA,betaABA common/mainAparam/Tdeac,PPFD,Vcmax25,Ca common/Aparam/convex,Tslope common/tous/Pair data pi,g/3.1415926,9.8/ ! g dans m s-2 data thetasval/0.338,4*0/ data thetarval/0.001,4*0/ data alphasolval/-172.,4*0./ data emsolval/0.1166,4*0./ data ensolval/1.132,4*0./ data rhobval/1.4,4*0./ data c0val/-15.1,4*0./ data c1val/60.4,4*0./ data c2val/-409.,4*0./ data c3val/1250.,4*0./ data Fval/5.,4*0./ * If this is the first time, ready everything * Si c'est la premiere fois, lire tous 60 if(ncha.eq.0)then write(6,'("Enter Vcmax25(umol...) Tdeac(oC) Tslope(oC-1)", + " convex aPAR")') read(5,*)Vcmax25,Tdeac,Tslope,convex,aPAR if(Vcmax25.lt.0.)stop * Convert to SI units * Convertir a SI Vcmax25=1.e-6*Vcmax25 * Enter the stomatal parameters * Entrer les params.stomatiques write(6,'(1x,"Enter mBB bBB hs0 tol",/,5x,"(hs<0 means:", + " use the true BB response; tol=precis. in gs")') * write(6,'(1x,"Entrer mBB bBB hs0 tol",/,5x,"(hs<0 veut", * + " dire: user la vraie reponse BB; tol=precis. dans gs")') read(5,*)mBB,bBB,hs0,tol write(6,'("Enter beta ((nM)-1; ca. -0.009), delta (MPa-1; ", + "ca. -1., or 0 for sunflower)")') * write(6,'("Entrer beta ((nM)-1; ca. -0.009), delta (MPa-1; ", * + "ca. -1., ou 0 chez tournesol)")') read(5,*)betaABA,deltaABA write(6,'("Enter Rhplant (MPa (mol m-2 s-1)-1; ca. 9, in sun.) + ")') * write(6,'("Entrer Rhplant (MPa (mol m-2 s-1)-1; ca. 9, chez t.) * + ")') read(5,*)Rhplant write(6,'("Enter aABA (nM mol m-2 s-1 MPa-1; ca. 2.05),", + " b (mol m-2 s-1; ca. 0.0072 in sunfl)")') * write(6,'("Entrer a (nM mol m-2 s-1 MPa-1; ca. 2.05),", * + " b (mol m-2 s-1; ca. 0.0072 chez t.)")') read(5,*)aABA,bABA write(6,'("Enter soil type (1=loam, ... the only type active", + " now)")') * write(6,'("Entrer isoltype (1=loam, ... le seul type actif", * + " maintenant)")') read(5,*)isoltype write(6,'("Enter mr (mass of roots per area, kg m-2),",/, + 5x,"rhor (dens. of r., kgDM m-3),",/,5x, + "rr (mean radius of active r., m), ",/,5x, + "Dsol (depth of soil, m)")') * write(6,'("Entrer mr (masse par aire des racine, kg m-2),",/, * + 5x,"rhor (dens. des r., kgDM m-3),",/,5x, * + "rr (moyen rayon des r. actives, m), ",/,5x, * + "Dsol (profondeur du sol, m)")') read(5,*)mr,rhor,rr,Dsol * Set the meteorological conditions * Poser les conditions micrometeorologiques * Read in micro-e, IF this is not a re-run write(6,2) 2 format(1x,"Enter Tair(oC) eair Tveg(oC) Tksy(oC)", + /,5x," (eair<0 means that it is relative humidity)") * 2 format(1x,"Entrer Tair(oC) eair Tveg(oC) Tciel(oC)", * + /,5x," (eair<0 veut dire que c''est la humidite relative") read(5,*)Tair, eair0,Tveg,Tsky write(6,3) 3 format(1x,'Enter PPFD(umol...) u (m/s) dleaf (m) ', + 'Ca(Pa) Pair (Pa)') * 3 format(1x,'Entrer PPFD(umol...) u (m/s) dfeuille (m) ', * + 'Ca(Pa) Pair (Pa)') read(5,*)PPFD,u,dleaf,Ca,Pair * Calculate gbair using u and dleaf * Calculer gbair usant u et dfeuille gbair=6.62e-3*sqrt(u/dleaf)*Pair/(8.314*(Tair+273.2)) ! From * egybal.gsBB.f * Convert to SI PPFD=1.e-6*PPFD * Set the environmental conditions of the soil * Poser les conditions environnementales du sol * psis=-1. ! MPa write(6,'("Enter psi(sol) (MPa; <0)")') * write(6,'("Entrer psi(sol) (MPa; <0)")') read(5,*)psis goto 90 else ! Laisser poser n'importe quel parametre, ou voir tous, ou * quitter, etc. * Allow setting of any parameter, or see them all, or quit, etc. 85 write(6,'("To modify parameters, enter the keyword", + " and the value, such as ''rhp 10.5''",/,5x,"To see the ", + "keywords and the other commands, enter ''cle''")') * 85 write(6,'("A modifier de parametres, entrer le mot-clef", * + " et le valeur, comme ''rhp 10.5''",/,5x,"A voir les mots-clef", * + " et les autres commandes, entrer ''cle''")') ncha=0 87 write(6,'("Enter the keyword and the value (", + "up to 20 caracteres)")') * 87 write(6,'("Entrer le mot-clef et le valeur (le total ", * + "jusqu''a 20 caracteres)")') read(5,'(a20)',err=88)chaval goto 89 ! S'il n'y a pas une erreur, continuer le traitment * If there isn't an error, continue the entries 88 write(6,'("Error; try again")') * 88 write(6,'("Erreur; essayer encore")') goto 87 ! Spaghetti-code inevitable, avec le controle des erreurs * Pick off the first 3 characters as a keyword * Enlever les premieres 3 caracteres commes un mot-clef 89 test=chaval(1:3) ncha=1 if(test.eq.'cle')then write(6,888) 888 format(1x,"The keywords (3-letter) are:",/, + 'vcm -> Vcmax25; tde -> Tdeac; tsl -> Tslope; con-> convex;', + ' apa -> aPAR;',/,10x,'mbb -> mBB; bbb -> bBB; hs0 -> hs0;', + ' tol -> tol;',/,10x,'tai -> Tair; eai -> eair tve -> Tveg;', + ' tsk -> Tsky;',/,10x,'ppf -> PPFD; u__ -> u; dle -> dleaf;', + ' cai -> Ca; pai -> Pair',/,10x,'bet -> betaABA; del -> ', + ' deltaABA', + /,5x,"rhp -> Rhplant; a__ -> aABA; b__-> bABA",/,5x, + "psi -> psis; iso -> isoltype;",/,8x,"mr_ -> mr;", + " rho -> rhor // la_ -> La; lv_ -> Lv;",/, + /,8x,"One should enter mr with rhor, or La with Lv)",/,5x, + "dso -> Dsol",/,5x,"sto -> STOP; pri -> ", + "print all values",/,5x,"all -> change all", + "(null) -> Begin run with the new values",/,5x, + "One can enter many changes, in succession") * 888 format(1x,"Les mots-clef (3-lettre) sont:",/, * + 'vcm -> Vcmax25; tde -> Tdeac; tsl -> Tslope; con-> convex;', * + ' apa -> aPAR;',/,10x,'mbb -> mBB; bbb -> bBB; hs0 -> hs0;', * + ' tol -> tol;',/,10x,'tai -> Tair; eai -> eair tve -> Tveg;', * + ' tci -> Tciel;',/,10x,'ppf -> PPFD; u__ -> u; dfe -> dfeuille;', * + ' cai -> Ca; pai -> Pair',/,10x,'bet -> betaABA; del -> ', * + ' deltaABA', * + /,5x,"rhp -> Rhplant; a__ -> aABA; b__-> bABA",/,5x, * + "psi -> psis; iso -> isoltype;",/,8x,"mr_ -> mr;", * + " rho -> rhor // la_ -> La; lv_ -> Lv;",/, * + /,8x,"On doit entrer mr avec rhor, ou La avec Lv)",/,5x, * + "dso -> Dsol",/,5x,"sto -> STOP; imp -> ", * + "iprimer tous valeurs",/,5x,"tou -> changer tous", * + "(null) -> commencer avec les nouveaux valeurs",/,5x, * + "On peut entrer beaucoup de changes, en succession") goto 85 elseif(test.eq.'all')then ! Read all values, again * elseif(test.eq.'tou')then ! Lire tous, encore ncha=0 goto 60 elseif(test.eq.'sto')then istop=1 return elseif(test.eq.'pri')then * elseif(test.eq.'imp')then write(6,25)Vcmax25,Tdeac,Tslope,convex,aPAR, + mBB,bBB,hs0,tol,Tair,eair,eair0,Tveg,Tsky,PPFD,u,dleaf, + Ca,Pair,betaABA,deltaABA,Rhplant,aABA,bABA,psis, + isoltype,mr,rhor,La,Lv,Dsol * write(6,25)Vcmax25,Tdeac,Tslope,convex,aPAR, * + mBB,bBB,hs0,tol,Tair,eair,eair0,Tveg,Tsky,PPFD,u,dfeuille, * + Ca,Pair,betaABA,deltaABA,Rhplant,aABA,bABA,psis, * + isoltype,mr,rhor,La,Lv,Dsol 25 format(3x,"Vcmax25=",6Pf7.2," umol m-2 s-1; Tdeac=", + 0Pf6.2," oC; Tslope=",f6.3," (oC)-1",/, + 5x,"convex=",f6.3,"; aPAR=",f6.3,/,3x,"mBB=",f6.2, + "; bBB=",f6.4," mol m-2 s-1; hs0=",f6.2,"; tol=", + f6.4," mol m-2 s-1",/,3x,"Tair=",f7.2,"oC; eair=",f7.1, + " Pa (eair0=",f7.1,")",/,5x,"Tveg=",f6.2,"oC; Tsky=", + f6.2,"oC; PPFD=",6Pf7.1," umol m-2 s-1",/,5x, + "u=",0Pf6.2," m s-1; dleaf=",f7.3," m; Ca=", + f7.2," Pa; Pair=",f8.1," Pa",/,3x,"betaABA=",f7.4, + " (nM)-1; deltaABA=",f7.2," MPa-1; Rhplant=",f6.1, + " MPa (mol m-2 s-1)-1",/,5x,"aABA=",f6.2," nM mol m-2", + " s-1; bABA=",f7.4," mol m-2 s-1",/,5x,"psis=",f7.2, + " MPa; isoltype=",i2,"; mr=",f6.1," kg m-2; rhor=", + f6.1," kg m-3",/,5x,"La=",f7.1," m m-2; Lv=",f7.1, + " m m-3; Dsol=",f6.3," m") * 25 format(3x,"Vcmax25=",6Pf7.2," umol m-2 s-1; Tdeac=", * + 0Pf6.2," oC; Tslope=",f6.3," (oC)-1",/, * + 5x,"convex=",f6.3,"; aPAR=",f6.3,/,3x,"mBB=",f6.2, * + "; bBB=",f6.4," mol m-2 s-1; hs0=",f6.2,"; tol=", * + f6.4," mol m-2 s-1",/,3x,"Tair=",f7.2,"oC; eair=",f7.1, * + " Pa (eair0=",f7.1,")",/,5x,"Tveg=",f6.2,"oC; Tsky=", * + f6.2,"oC; PPFD=",6Pf7.1," umol m-2 s-1",/,5x, * + "u=",0Pf6.2," m s-1; dfeuille=",f7.3," m; Ca=", * + f7.2," Pa; Pair=",f8.1," Pa",/,3x,"betaABA=",f7.4, * + " (nM)-1; deltaABA=",f7.2," MPa-1; Rhplant=",f6.1, * + " MPa (mol m-2 s-1)-1",/,5x,"aABA=",f6.2," nM mol m-2", * + " s-1; bABA=",f7.4," mol m-2 s-1",/,5x,"psis=",f7.2, * + " MPa; isoltype=",i2,"; mr=",f6.1," kg m-2; rhor=", * + f6.1," kg m-3",/,5x,"La=",f7.1," m m-2; Lv=",f7.1, * + " m m-3; Dsol=",f6.3," m") goto 85 elseif(test.eq.'')then if(ncha.eq.0)istop=1 ! Pas de changement; juste stopper goto 90 ! Calc. les param. du sol ...et retourner,a faire un * passage de nouveau, ou a stopper else * Pick up the numerical value, as a character string * Enlever le valeur numerique, comme une chaine alphabetique sval=chaval(4:) * Then, convert it by writing to disk and re-reading as a number * Puis, convertir-la par ecrire au lecteur et re-lire comme un numero rewind 8 write(8,*)sval rewind 8 read(8,*,err=88)val * Now put the numerical value with the correct parameter * Puis, attribuer le valeur au correct parametre if(test.eq.'vcm')then Vcmax25=val*1.e-6 ! Convertir a mol (de micromol) goto 85 elseif(test.eq.'tde')then Tdeac=val goto 85 elseif(test.eq.'tsl')then Tslope=val goto 85 elseif(test.eq.'con')then convex=val goto 85 elseif(test.eq.'apa')then aPAR=val goto 85 elseif(test.eq.'mbb')then mBB=val goto 85 elseif(test.eq.'bbb')then bBB=val goto 85 elseif(test.eq.'hs0')then hs0=val goto 85 elseif(test.eq.'tol')then tol=val goto 85 elseif(test.eq.'tai')then Tair=val goto 85 elseif(test.eq.'eai')then eair0=val goto 85 elseif(test.eq.'tve')then Tveg=val goto 85 elseif(test.eq.'tsk')then Tsky=val goto 85 elseif(test.eq.'ppf')then PPFD=val*1.e-6 ! Convertir a mol (de micromol) * Convert to moles (from micrmol) goto 85 elseif(test.eq.'u__')then u=val goto 85 elseif(test.eq.'dle')then dleaf=val goto 85 elseif(test.eq.'cai')then Ca=val goto 85 elseif(test.eq.'pai')then Pair=val goto 85 elseif(test.eq.'bet')then betaABA=val goto 85 elseif(test.eq.'del')then deltaABA=val goto 85 elseif(test.eq.'rhp')then Rhplant=val goto 85 elseif(test.eq.'a__')then aABA=val goto 85 elseif(test.eq.'b__')then bABA=val goto 85 elseif(test.eq.'psi')then psis=val goto 85 elseif(test.eq.'iso')then isoltype=ifix(val) goto 85 elseif(test.eq.'mr_')then mr=val *XXXXXX entrer aussi rhor... / Enter also rhor if you use this option goto 85 elseif(test.eq.'rho')then rhor=val *XXXXXX entrer aussi mr... / Enter also mr if you use this option goto 85 elseif(test.eq.'la_')then La=val write(6,'("It''s necessary to enter Lv also")') * write(6,'("Il faut entrer aussi Lv")') goto 85 elseif(test.eq.'lv_')then Lv=val write(6,'("It''s necessary to enter La also")') * write(6,'("Il faut entrer aussi La")') goto 85 elseif(test.eq.'rr_')then rr=val goto 85 elseif(test.eq.'dso')then Dsol=val goto 85 else write(6,'("Keyword unrecognizable; try again")') * write(6,'("Mot-clef meconnaisable; essayer encore")') goto 85 endif * Allow more changes, in succession * Laisser plus de changes, en succession endif endif 90 continue * Calculer eair, etc. esatair=esat(Tair) if(eair0.lt.0.)then eair=-eair0*esatair else eair=eair0 endif ha=eair/esatair * Fixer les valeurs des parametres du sol, par le type du sol thetas=thetasval(isoltype) ! teneur d'eau, sature' (-) * Soil water content, saturated (-) thetar=thetarval(isoltype) ! le meme, residuel (-) * The same, residual alphasol=alphasolval(isoltype) ! mult. de psi dan l'eq. de van G. * Multiplier of psi in the Eq. of van Genuchten ensol=ensolval(isoltype) ! expon. de psi emsol=emsolval(isoltype) ! expon. dans (1+(alpha*psi)**n)**(-m) rhob=rhobval(isoltype) c0=c0val(isoltype) c1=c1val(isoltype) c2=c2val(isoltype) c3=c3val(isoltype) theta=thetar+(thetas-thetar)*(1.+(alphasol* + psis)**ensol)**(-emsol) * Calculate La = the density of root length, per area (as m m-2 = m-1), * and d = the semi-distance among the closest roots, as if they are * spaced completely regularly * Calculer La = la densite du longueur des racines, par aire (comme * m m-2 = m-1), et d = la semi-distance entre les plus proches * racines, si ils sonts espaces completement regulairement d=0.5*rr*sqrt(rhor*pi*Dsol/mr) La=mr/(rhor*pi*rr**2) * Fix the value of the "factor of Tardieu," which multiplies Rsol because * of the clumping of roots * Fixer le valeur du "Facteur de Tardieu", F, qui multiplie Rsol * a cause du groupement des racines F=Fval(isoltype) * Calculer le constant dans R(sol,air,flux-masse) cRsol=F*g*log(d/rr)/(2.*pi*La) thetam=theta/rhob arg=c0+thetam*(c1+thetam*(c2+thetam*c3)) khstar=10.**arg ! m s-1 , dans les unites de "head" ! m s-1, in units of "head" (hydraulic) Rsol=cRsol/khstar ! la vielle R_s_a_fm [R(sol,aire,flux-masse)] * The old R_s_a_fm (in Tardieu's Excel spreadsheet) Rsol=0.1*Rsol ! XXXX peut-etre pas justifiee! * Maybe not justified! return end * * * DESCRIPTION OF VARIABLES * Collating sequence is numbers first, then letters; capital and * lower case are not distinguished unless they compete, and then * capital letter takes precedence * Units are given in parentheses; "-" means unitless * * LOGICAL * BBswitch = indicator that true Ball-Berry response is to be used * (as opposed to having responses to A and hs and Cs turned off) * * CHARACTER VARIABLES * chaval = acronym for parameter whose value is to be changed * test = holding variable for "chaval" * * CHARACTER ARRAYS * string * * REAL SCALAR VARIABLES (and specified units) * A = leaf CO2 assimilation rate (mol m-2 s-1; printed out as * micromol m-2 s-1) * aABA = const. in kinetic eq. for synthesis of ABA (nM MPa-1 mol-1 m2 s) * ALL = light-limited rate of carboxylation (mol m-2 s-1) * alphasol = for this soil, the value of parameter alpha in the * van Genuchten relation of theta to psi (MPa-1) * Amax = light-saturated rate of carboxylation (mol m-2 s-1) * amult = thermal-deactivation factor for carboxylation capacity (-) * aNIR = leaf absorptivity for PAR (-) * aPAR = leaf absorptivity for NIR (-) * arg = temporary variable for input to exponentials (-) * Atop = estimate of max. CO2 assim. rate, at T = Tdeac (mol m-2 s-1) * bABA = const. in kinetic eq. for synthesis of ABA (mol m-2 s-1) * bBB = Ball-Berry intercept (mol m-2 s-1) * Bcc = convective heat-transfer coefficient of leaf (W m-2 K-1) * betaABA = sensitivity of gs to [ABA] (nM-1) * Bnet = temperature derivative of total heat gain rate of * leaf (W m-2 K-1) * BTIR = radiation "conductance (d[QTIR]/dT) (W m-2 K-1) * c0, c1, c2, c3 = values of coefficients in empirical * formula for k,h,soil in terms of mass-based soil water * content (-; there is an implicit dimension ahead of the * exponential formula in which these are used) * Ca = partial pressure of CO2 in ambient air (Pa) * cABA = conc. of ABA in xylem stream (nM) * Citop = estimate of maximal value of Ci (Pa) * convex = convexity parameter in photosynthetic response (also * called theta in explan. of model) (-) * Cpair = molar heat capacity of air at constant pressure (J mol-1 K-1) * cRsol = intermediate quantity in calculating hydraulic resistance * from bulk soil to root surface (m) * Cs = partial pressure of CO2 at leaf surface, beneath leaf * boundary layer (Pa) * d = mean inter-root distance in soil (m) * dCi = increment in Ci, in searching for solution of quartic eq. for * Ci (Pa) * deltaABA = response param. for amplification of ABA response by low * leaf water potential (MPa-1) * denom = temporary variable in T-dependence of photosyn. params. * (J mol-1 K) * dfac = factor in T-dependence of photosyn. parameters (J mol-1) * dleaf = leaf charac. linear dimension (m) * dpsiPa = drop in water potential from bulk soil to root (Pa) * Dsol = depth of soil in rooting volume (m) * dTleaf = latest increment in solving for steady-state leaf T (K) * E = leaf transpiration rate per area (mol m-2 s-1) * NOTE: this is also set to be transpiration per ground area as a * first approx.; see the typeset (PDF) file giving the explanation * of the model * eair = partial pressure of water vapor in ambient air (Pa) * eair0 = original form of entering eair value; if <0, it means it * is to be interpreted as a relative humidity * eL = partial pressure of water vapor inside leaf (Pa) * emsol = value of exponent m in van Genuchten relation between * theta and psi, for this soil (-) * ensol = value of exponent n in van Genuchten relation between * theta and psi, for this soil (-) * ep = "e-prime" = derivative of saturated partial pressure of water * vapor with respect to temperature (Pa K-1) * epsleaf = leaf TIR emissivity (-) * esatair = saturated partial pressure of water vapor at air * temperature (Pa) * F = test for convergence of quartic eq. for Ci (-) * f1, f2, f3, f4 = convenient combinations of parameters that occur * in quartic solution for Ci (various units) * facABA = gfac in text = multiplier (<1) of gs, from ABA-induced * stomatal closure (-) * Fprime = d(F)/d(Ci) in Newton-Raphson search for solution of * quartic ef. for Ci (Pa-1) * g = acceleration of gravity (m s-1) * gamma = CO2 compensation partial pressure (Pa) * gammatop = CO2 compensation partial pressure at T = Tdeac (Pa) * gbair = boundary-layer conductance of leaf (mol m-2 s-1) * gs = final solution for stomatal conductance (mol m-2 s-1) * gs0 = estimate of max. stomatal conductance (mol m-2 s-1) * gsadd = first part of estimate of max. gs, at T= Tdeac (mol m-2 s-1) * gshi = estimated upper limit of stomatal conduc. (mol m-2 s-1) * gslo = estimated lower limit of stomatal conduc. (mol m-2 s-1) * gsmax = estimate of max. gs, at T= Tdeac (mol m-2 s-1) * gsminus = a holding variable for lower limit of gs (mol m-2 s-1) * gsplus = a holding variable for upper limit of gs (mol m-2 s-1) * gtotCO2 = total (stomatal+boundary-layer) conduc. for CO2 * (mol m-2 s-1) * gtotw = total (stomatal+boundary-layer) conduc. for water vapor * (mol m-2 s-1) * ha = relative humidity in air, as a fraction (-) * hs = relative humidity at leaf surface (-) * hs0 = reference value of rel. humidity at leaf surface, to which * the leaf has acclimated (-) * Hvap (old "lambda") = heat of vap. of water (J mol-1) * IBBhi = Ball-Berry index (A hs/Cs) for upper limit of gs (gshi) * (mol m-2 s-1) * IBBlo = similarly, for lower limit of gs (gslo) * IBBmid = similarly, for midpoint of binary search (gsmid) * Kc = Michaelis constant for CO2 binding to Rubisco (Pa) * Kc25 = Michaelis constant for CO2 binding to Rubisco, at * reference T = 25oC (Pa) * Kco = Michaelis constant for CO2 binding to Rubisco, including the * O2 competition effect (Pa) * Kcotop = effective Michaelis const. for CO2 binding (including * O2 competition effect), at T = Tdeac (Pa) * Kctop = Michaelis const. for CO2 binding, evaluated at T = Tdeac (Pa) * Kc = Michaelis constant for O2 binding to Rubisco (Pa) * khstar = hydraulic conductivity of soil (mol m-2 s-1 MPa-1) * Ko25 = Michaelis constant for O2 binding to Rubisco, at * reference T = 25oC (Pa) * Kotop = Michaelis const. for O2 binding, evaluated at T = Tdeac (Pa) * La = root length density, per ground area (m m-2, which is m-1) * Lv = XXXX * mBB = Ball-Berry slope (-) * mr = mass of roots per ground area (kg m-2) * Pair = total pressure of ambient air (Pa) * pi = 3.141592... * PPFD = photon (quantum) flux density in PAR, at leaf surface * (mol m-2 s-1; read in as micromol m-2 s-1) * psiL = leaf water potential (MPa) * psir = root water potential (MPa) * psis = soil water potential (MPa) * Q00 = CO2-saturated initial quantum yield (mol[CO2] mol[photons]-1) * Qmcc = rate of convective heat loss from leaf (W m-2) * Qme = rate of transpirational (latent) heat loss from leaf (W m-2) * Qnet = total rate of heat gain by leaf (should get to zero) (W m-2) * Qmrad = rate of TIR loss from 2 sides of leaf (W m-2) * Qrad = total (PAR+NIR+TIR) energy-flux density (W m-2) * Qsky = TIR energy-flux density from sky, onto flat surface (W m-2) * Qsky = TIR energy-flux density from surr. vegetation, onto flat * surface (W m-2) * QSW = shortwave (PAR+NIR) energy-flux density (W m-2) * rhob = specific gravity of this soil (-) * rhor = dry-matter density in fresh roots (kg m-3) * Rhplant = hydraulic resistance within the plant, from root to * leaf (MPa mol-1 m2 s) * rr = mean radius of (fine or active) roots (m) * Rsol = hydraulic resistance from bulk soil to root (Pa-1 mol-1 * m2 s) * sigma = Stefan-Boltzmann constant * t0, t1, t2, t3, t4 = coefficients in quartic equation for Ci * (units of Pa**[-n], n = 1, 1, 2, 3, 4) * Tair = temperature of ambient air (oC) * Tdeac = temperature above which carbox. capacity begins * to deactivate thermally (oC) * testhi = criterion for gs being above actual solution, evaluated * for gs = gshi (mol m-2 s-1) * testlo = criterion for gs being above actual solution, evaluated * for gs = gslo (mol m-2 s-1) * testmid = criterion for gs being above actual solution, evaluated * for gs = gsmid (mol m-2 s-1) * theta = actual volumetric water content of this soil, at specified * water potential (-) * thetam = mass-fraction water content of soil (-) * thetar = residual volumetric water content for this soil (-) * thetas = saturated volumetric water content for this soil (-) * TL = leaf temperature (oC) * tol = tolerance within which we want to solve for gs (mol m-2 s-1) * Tsky = effective radiative temperature of sky (oC) * Tslope = rate at which thermal deactivation of carboxylation capacity * occurs above T = Tdeac (K-1) * Tveg = effective radiative temperature of surrounding vegetation (oC) * u = windspeed (m s-1) * Vcmax = actual carboxylation capacity per leaf area at leaf T * (mol m-2 s-1) * Vcmax25 = max. carboxylation capacity per area of leaf * (mol m-2 s-1; converted from initial entry as micromol m-2 s-1) * Vcmaxtop = max. carboxylation capacity per area of leaf, at highest * T before thermal deactiv. starts (mol m-2 s-1) * * REAL ARRAYS * * alphasolsval = value of alpha, in van Genuchten relation between * soil water content and water potential, for each of 5 soil types (-) * c0val, c1val, c2val, c3val = values of coefficients in empirical * formula for k,h,soil in terms of mass-based soil water * content (-; there is an implicit dimension ahead of the * exponential formula in which these are used) * emsolsval = value of exponent m, in van Genuchten relation between * soil water content and water potential, for each of 5 soil types (-) * ensolsval = value of exponent n, in van Genuchten relation between * soil water content and water potential, for each of 5 soil types (-) * Fval = value of root clumping factor, F, for each of 5 soil types (-) * rhobval = value of specific gravity of bulk soil, for each of 5 * soil types (-) * thetasval = value of saturated volum. water content of soil, for * each of 5 soil types (-) * thetarval = value of residual volum. water content of soil, for * each of 5 soil types (-) * * INTEGER VARIABLES * isoltype = identifier of soil texture type * istop = indicator that run should be stopped (if it =1) * iter = counter for number of iterations of the quartic eq. for Ci * ncha = indicator of which parameter has been selected to be changed * nruns = total number of binary-search halvings in searching for * gs solution (-) * * * INTEGER DO-LOOP INDICES * n * * SUMMARY OF SUBROUTINES * paramread - reads in all parameter values, and allows any number of * changes, particularly for a new run (one need specify only the * changes) * CiACssolve - given the leaf envir., it computes Ci, A, and Cs * Tleaf - given energy-balance parameters and envir. vars., it * computes leaf temperature (and transpir. rate) * esat - computes saturated partial pressure of water vapor at * specified temperature * Fval - returns value of clumping factor, F, for this soil type