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