DECLARE FUNCTION TBOLTZ! (RATE!, EA!, TOPT!, TL!)

DECLARE FUNCTION TEMP! (RATE25!, EACT!, TPRIME, TREF)

DECLARE FUNCTION ACOS! (X!)

COMMON SHARED H, RUGC

CLS

'

'November 26, 1993

'

'PSCUBIC.BAS Dennis Baldocchi, NOAA/ATDD. Oak Ridge, TN

'

' This program solves a cubic equation to calculate

' leaf photosynthesis. This cubic expression is derived from solving

' five simultaneous equations for A, PG, CS, CI and GS.

' Stomatal conductance is computed with the Ball-Berry model.

' The cubic derivation assumes that b', the intercept of the Ball-Berry

' stomatal conductance model, is non-zero.

'

' Gs = k A rh/Cs + b'

'

' Derived in Viterbo, Italia

' Dec. 4, 1992

'

'*************************************

'

'

OPEN "d:\psoak\pscubic.jnk" FOR OUTPUT AS #1

PI = 3.141593

'

' Gas exchange parameters are for sunlit oak. Data

' are derived from leaf exchange measurements by

' Peter Harley, NCAR during the summer of 1992 at

' Walker Branch Watershed

'*

'* Enzyme constants & partial pressure of O2 and CO2

'*

' Michaelis-Menten K values.

'

KC25 = 275 'microbars

KO25 = 420 'millibars

O2 = 210 'umol mol-1

'

TAU25 = 2321.2

'

'*

'* Arrhenius constants

'*

'Eact for Michaelis-Menten const. for KC, KO and dark respiration

'

EKC = 80500! 'CO2 J mol-1

EKO = 14500! ' O2

ERD = 38000 'Dark respiration,

EKTAU = -29000

'

ZCK = 273.16

T25 = 298.16

EPS = .1

'

'Universal gas constant

'

RUGC = 8.314

RGC1000 = RUGC * 1000

'

'consts for Boltzmann distrib.: in eq 36 Farquhar et al. (1980)

' For Vcmax and Jmax

'

H = 200000!

S = 710

EJM = 55000

EVC = 55000

'

PI = 3.14159265#

PI4 = PI * 4

PO = PI / 180

BPRIME = .0175

BPRIME16 = BPRIME / 1.6

'

' Wj=J Ci/(4 Ci + 8 GAMMA)

'

vcopt = 73 'umol m-2 s-1

jmopt! = 170 'umol e- m-2 s-1

rd25 = .34

vcopt = 25

jmopt! = 29 + 1.64 * vcopt

rd25 = .015 * vcopt

QALPHA = .2 'mole CO2 per mole electron

QALPHA2 = QALPHA * QALPHA

'

KBALL = 9.5

'

STRTTIM = TIMER

'

'PSCUBIC.BAS Dennis Baldocchi, NOAA/ATDD. Oak Ridge, TN

'

' This program solves a cubic equation to calculate

' leaf photosynthesis. This cubic expression is derived from solving

' five simultaneous equations for A, PG, CS, CI and GS.

' Stomatal conductance is computed with the Ball-Berry model.

' The cubic derivation assumes that b', the intercept of the Ball-Berry

' stomatal conductance model, is non-zero.

'

' Gs = k A rh/Cs + b'

'

' Define terms for calculation of gross photosynthesis, PG

'

' PG is a function of the minimum of RuBP saturated rate of

'carboxylation, Wc, and the RuBP limited rate of carboxylation, Wj.

' Wj is limiting when light is low and electron transport, which

're-generates RuBP, is limiting. Wc is limiting when plenty of RuBP is

'available compared to the CO2 that is needed for carboxylation.

' Both equations take the form:

'

' PG-photorespiration= (a CI-a d)/(e CI + b)

'

' PG-photorespiration=min[Wj,Wc] (1-GAMMA/Ci)

'

' Wc=Vcmax Ci/(Ci + Kc(1+O2/Ko))

'

'*****************************************************************

'

' Viterbo, Italia

' Nov. 26, 1992

'

'*************************************

'

'

' BIOCHEMICAL MODEL OF C3 PHOTOSYNTHESIS

'

' After Farquhar, von Caemmerer and Berry (1980) Planta.

' 149: 78-90.

'*

'**************************************************************

'

'* Program calculates leaf photosynthesis from biochemical parameters

'*

'* RD25 - Dark respiration at 25 degrees C (umol m-2 s-1)

'* TLEAF - leaf temperature

'* JMAX - rate of electron transport at optimal temp degrees C

'* VCOPT - maximum rate of RuBP Carboxylase/oxygenase at opt T degrees C

'* IPHOTON - absorbed photosynthetically active photon flux (mmols m-2 s-1)

'* Gs - stomatal conductance (mols m-2 s-1), typically 0.01-0.20

'* PSTAT-station pressure, bars

'* APHOTO - net photosynthesis (umol m-2 s-1)

'* PS - gross photosynthesis (umol m-2 s-1)

'* AMG - net photosynthesis (mg m-2 s-1)

'*

'*************************************************************

'

'Leaf temperature in C and K

'

'FOR kk = 10 TO 20

TL = 0 + kk * 2.5

TL = 25

tlk = TL + 273.16

PSTAT = .975

TOPTJM = 38 + 273

TOPTVC = 38 + 273

'FOR K = 0 TO 10

'FOR J = 0 TO 4

ca = 250 + 100 * J

ca = 350

'

'Convert Par from W m^-2 to

'irradiance in uE m-2 s-1

'

FOR I = 0 TO 20

IPHOTON = I * 25' absorbed quanta is inputted

'IPHOTON = 2000

'

'IPHOTON is radiation incident on leaves

'

' The temperature dependency of the kinetic properties of

' RUBISCO are compensated for using the Arrhenius and

' Boltzmann equations. From biochemistry, one observes that

' at moderate temperatures enzyme kinetic rates increase

' with temperature. At extreme temperatures enzyme

' denaturization occurs AND rates must decrease.

'

' Arrhenius Eq.

'

'f(T)=f(T25) exp(Tk -298)Eact/(298 R Tk)), where Eact is the

'activation energy.

'

' Boltzmann distribution

'

' F(T)=TBOLTZ)

'*

'* Temperature corrections for enzyme kinetics

'*

RT = RUGC * tlk

TPRIME25 = tlk - T25

'

' KC and KO are solely a function of the Arrhenius Eq.

'

KCT = TEMP(KC25, EKC, TPRIME25, T25)

KO = TEMP(KO25, EKO, TPRIME25, T25)

TAU = TEMP(TAU25, EKTAU, TPRIME25, T25)

BC = KCT * (1 + O2 / KO)

'

'GAMMAC is the CO2 compensation point due to photorespiration, umol mol-1

'Recalculate GAMMAC with the new temperature dependent KO and KC

' coefficients

'

'GAMMAC = .5 * O2*1000/TAU ubar

'

GAMMAC = .5 * 1000 * O2 / TAU

'

' temperature corrections for Jmax and Vcmax

'

JMAXZ! = jmopt!

VCMAXZ = vcopt

'

'Apply temperature correction to

' JMAX and VCMAX

'

JMAX! = TBOLTZ(JMAXZ!, EJM, TOPTJM, tlk)

VCMAX = TBOLTZ(VCMAXZ, EVC, TOPTVC, tlk)

'

'Scale Rd with height via VCMAX and apply temperature

'correction for dark respiration

'

RD = TEMP(rd25, ERD, TPRIME25, T25)

'

'**********************************************

' GB leaf boundary layer conductance, mol m-2 s-1

'for CO2 exchange

'*

'* Compute the leaf boundary layer resistance

'*

' RB = 10 ^ (K / 3)

RB = 10

'

'RB has units of s/m, convert to mol-1 m2 s1 to be

'consistant with R.

'

RBB = RB * .022624 * tlk / (PSTAT * 273.16)

GB = 1 / RBB

'

DD = GAMMAC

'***************************************

'

' APHOTO = PG - RD, net photosynthesis is the difference

' between gross photosynthesis and dark respiration. Note

' photorespiration is already factored into PG.

'

'****************************************

'

' coefficients for Ball-Berry stomatal conductance model

'

' Gs = k A rh/Cs + b'

'

' rh is relative humidity, which comes from a coupled

'leaf energy balance model.

'

'remember Gs is for water vapor and must be converted to a

' CO2 conductance for the photosynthesis calculations

'

RH = .85

'

KRH = RH * KBALL

'

'convert to a CO2 conductance

'

KRH = KRH / 1.6

GBKRH = GB * KRH

'

CIGUESS = ca * .7

'

'cubic coefficients that are only dependent on CO2 levels

'

ALPHA = 1 + (BPRIME16 / GB) - (KRH)

BETA = ca * (GBKRH - 2 * BPRIME16 - GB)

GAMMA = ca * ca * GB * BPRIME16

THETA = GBKRH - BPRIME16

'

'Test for the minimum of Wc and Wj. Both have the form:

'

' W = (a Ci - ad)/(e Ci + b)

'

' after the minimum is chosen set a, b, e and d for the cubic solution.

'

'J photon from Harley

'

JPHOTON = QALPHA * IPHOTON / (1 + (QALPHA2 * IPHOTON * IPHOTON / (JMAX! * JMAX!))) ^ .5

WJ = JPHOTON * (CIGUESS - DD) / (4 * CIGUESS + 8 * DD)

'

WC = VCMAX * (CIGUESS - DD) / (CIGUESS + BC)

'

IF WJ < WC THEN

'

' for Harley and Farquhar type model for Wj

'

B = 8 * DD

A = JPHOTON

E = 4

'

ELSE

B = BC

A = VCMAX

E = 1

END IF

IF WJ <= 0 THEN GOTO quad

IF WC <= 0 THEN GOTO quad

'

'cubic solution:

'

' A^3 + p A^2 + q A + r = 0

'

'

DENOM = E * ALPHA

'

PCUBE = (E * BETA + B * THETA - A * ALPHA + E * RD * ALPHA)

PCUBE = PCUBE / DENOM

'

QCUBE = (E * GAMMA + (B * GAMMA / ca) - A * BETA + A * DD * THETA + E * RD * BETA + RD * B * THETA)

QCUBE = QCUBE / DENOM

'

RCUBE = (-A * GAMMA + A * DD * (GAMMA / ca) + E * RD * GAMMA + RD * B * GAMMA / ca)

RCUBE = RCUBE / DENOM

' Use solution from Numerical Recipes from Press

'

P2 = PCUBE * PCUBE

P3 = P2 * PCUBE

Q = (P2 - 3 * QCUBE) / 9

r = (2 * P3 - 9 * PCUBE * QCUBE + 27 * RCUBE) / 54

'

TEST = Q ^ 3 - r ^ 2

IF TEST <= O THEN

LPRINT IPHOTON, TL, ca, RB

GOTO outdat

END IF

'

'ALL ROOTS ARE REAL

'

ARGU = r / (Q * Q * Q) ^ .5

ANGL = ACOS(ARGU)

ROOT3 = -2 * Q ^ .5 * COS((ANGL + PI4) / 3) - PCUBE / 3

'END IF

'

' Here A = x - p / 3, allowing the cubic expression to be expressed

' as:

' x^3 + ax + b = 0

'

APHOTO = ROOT3

'

'also test for sucrose limitation of photosynthesis, as suggested by

' Collatz. Js=Vmax/2

'

JSUCROSE = VCMAX / 2 - RD

'

IF JSUCROSE < APHOTO THEN

APHOTO = JSUCROSE

END IF

'

CS = ca - APHOTO / GB

'

' Stomatal conductance for water vapor

'

'

'forest are hypostomatous.

' Hence we don't divide the total resistance

' by 2 since transfer is going on only one side of a leaf.

'

'

GS = (KBALL * RH * APHOTO / CS) + BPRIME

'

' convert Gs from vapor to CO2 diffusion coefficient

'

GSCO2 = GS / 1.6

'

'stomatal conductance is mol m-2 s-1

'convert back to resistance (s/m) for energy balance routine

'

GSS = GS * .022624 * tlk / (PSTAT * 273.16)

RSTOM = 1 / GSS

'

' to compute CI, Gs must be in terms for CO2 transfer

'

CI = CS - APHOTO / GSCO2

'

' if A < 0 then gs should go to cuticular value and recalculate A

' using quadratic solution

'

IF APHOTO < 0 THEN GOTO quad

GOTO outdat

'IF APHOTO < 0 set stomatal conductance to cuticle value

quad:

'

GS = BPRIME

GSCO2 = GS / 1.6

'

'stomatal conductance is mol m-2 s-1

'convert back to resistance (s/m) for energy balance routine

'

GSS = GS * .022624 * tlk / (PSTAT * 273.16)

RSTOM = 1 / GSS

'

PSI1 = ca * GB * BPRIME16

DELTA1 = BPRIME16 + GB

DENOM = GB * BPRIME16

'

AQUAD1 = DELTA1 * E

BQUAD1 = -PSI1 * E - A * DELTA1 + E * RD * DELTA1 - B * DENOM

CQUAD1 = A * PSI1 - A * DD * DENOM - E * RD * PSI1 - RD * B * DENOM

'

APHOTO = (-BQUAD1 - (BQUAD1 * BQUAD1 - 4 * AQUAD1 * CQUAD1) ^ .5) / (2 * AQUAD1)

'

' Tests suggest that APHOTO2 is the correct photosynthetic root when

' light is zero because root 2, not root 1 yields the dark respiration

' value RD.

'

CS = ca - APHOTO / GB

CI = CS - APHOTO / GSCO2

'END IF

'

outdat:

'

'compute photosynthesis with units of mg m-2 s-1

'

'AMG = APHOTO * 44 / 1000

'

AMG = APHOTO * .044

LOCATE 2, 3

PRINT " CS CI GS CA CI/CA WC WJ WSUC"

PRINT USING "####.# ####.# #.##### ##### #.#### ###.## ###.## ###.##"; CS; CI; GS; ca; CI / ca; WC; WJ; JSUCROSE

PRINT ""

PRINT " A CUBIC RESP "

PRINT APHOTO, RD

PRINT ""

'**********************************************

'iterative solution

'

CIGUESS = .7 * ca

'

' test for minimum of Wj and Wc. Carboxylation rates limited by

' regeneration of RuBP by electron tranport (Wj) or carboxylation

' limited by Ci (Wc). Both have the form:

'

' PG= (a CI-a d)/(CI + b)

'

' After minimum(Wj,Wc) is chosen set a b and d for quadratic solution

'

'J photon from Harley

'

JPHOTON = QALPHA * IPHOTON / (1 + (QALPHA2 * IPHOTON * IPHOTON / (JMAX! * JMAX!))) ^ .5

WJ = JPHOTON * (CIGUESS - DD) / (4 * CIGUESS + 8 * DD)

'

WC = VCMAX * (CIGUESS - DD) / (CIGUESS + BC)

'

IF WJ < WC THEN

B = 8 * DD

A = JPHOTON

E = 4

ELSE

B = BC

A = VCMAX

E = 1

END IF

'

' now do iterative solution of PG, A, CI, CS and GS and

'see if they are the same as the quadratic solution.

'

'first guess CI equal CA

'

CIIT = ca

'

PGIT = (A * CIIT - A * DD) / (E * CIIT + B)

AIT = PGIT - RD

AOLD = AIT

'

COUNT = 0

10 CSIT = ca - AIT / GB

'

IF AIT > 0 THEN

GSIT = (KBALL * RH * AIT / CSIT) + BPRIME

ELSE GSIT = BPRIME

END IF

CIIT = CSIT - 1.6 * AIT / GSIT

PGIT = (A * CIIT - A * DD) / (E * CIIT + B)

AIT = PGIT - RD

IF ABS(AIT - AOLD) > .001 THEN

AOLD = AIT

COUNT = COUNT + 1

IF COUNT > 100 THEN GOTO ENDLOOP

GOTO 10

END IF

ENDLOOP:

'LOCATE 10, 3

PRINT " A ITERATED photon flux RESP COUNT"

PRINT USING "####.### #####. #####.##### ###.## ####"; AIT; IPHOTON; RD; COUNT

PRINT ""

PRINT " CS CI GS CA CI/CA WC WJ WSUC"

PRINT USING "####.# ####.# #.##### ##### #.#### ###.## ###.## ###.##"; CSIT; CIIT; GSIT; ca; CIIT / ca; WC; WJ; JSUCROSE

PRINT ""

PRINT " GB TL"

PRINT GB, TL

'

ENDIT = TIMER

'

QUADTIME = ENDTIME - STRTTIM

ITERTIME = ENDIT - ENDTIME

'LOCATE 11, 3

PRINT ""

PRINT " QUADTIME ITERTIME"

PRINT QUADTIME, ITERTIME

'

'CLS

WRITE #1, RB, ca, tlk, IPHOTON, APHOTO, RD, AIT, GS, RSTOM

PRINT ca, IPHOTON, APHOTO, AIT, GS, RSTOM

NEXT I

'NEXT J

'NEXT K

'NEXT kk

FUNCTION ACOS (X)

'

IF ABS(X * X) > 1 THEN

PRINT " CAN'T CALCULATE ACOS"

STOP

END IF

'

ACOS = 1.570796 - ATN(X / SQR(1 - X * X))

END FUNCTION

DEFINT J

FUNCTION TBOLTZ (RATE, EA, TOPT, TL)

NUMM = RATE * H * EXP(EA * (TL - TOPT) / (RUGC * TOPT * TL))

DENOM = H - EA * (1 - EXP(H * (TL - TOPT) / (RUGC * TOPT * TL)))

TBOLTZ = NUMM / DENOM

END FUNCTION

DEFSNG J

FUNCTION TEMP (RATE, EACT, TPRIME, TREF) STATIC

SHARED RT

TEMP = RATE * EXP(TPRIME * EACT / (TREF * RT))

END FUNCTION