function [ps]=LeafPs(Iphoton,cca, tlk, rb_co2, P_kPa, eair_Pa,vcopt)


%         12-21-2011, Leaf_photosynthesis model   


%         DENNIS BALDOCCHI
%         Ecosystem Science Division
%         Department of Environmental Science, Policy and Management
%         & Berkeley Atmospheric Science Center
%         345 Hilgard Hall
%         University of California, Berkeley
%         Berkeley, CA 94720-3110
%         baldocchi@berkeley.edu
%         510-642-2874
       
%         
%          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'
% 
%           We also found that the solution for A can be obtained by a quadratic equation
%           when Gs is constant or b' is zero.
% 
% 
%             The derivation is published in:
% 
%             Baldocchi, D.D. 1994. An analytical solution for coupled leaf photosynthesis
%             and stomatal conductance models. Tree Physiology 14: 1069-1079.
% 
% 
% -----------------------------------------------------------------------
% 
%           A Biochemical Model of C3 Photosynthesis
% 
%             After Farquhar, von Caemmerer and Berry (1980) Planta.
%             149: 78-90.
% 
%         The original program was modified to incorporate functions and parameters
%         derived from gas exchange experiments of Harley, who paramertized Vc and J in
%         terms of optimal temperature, rather than some reference temperature, eg 25C.
% 
%         Program calculates leaf photosynthesis from biochemical parameters
% 
%         rd25 - Dark respiration at 25 degrees C (umol m-2 s-1)
%         tlk - leaf temperature, Kelvin
%         jmax - optimal rate of electron transport
%         vcopt - maximum rate of RuBP Carboxylase/oxygenase
%         iphoton - incident photosynthetically active photon flux (mmols m-2 s-1)
%    
%             note: Harley parameterized the model on the basis of incident PAR
% 
%         gs - stomatal conductance (mol 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)
%         aps - net photosynthesis (mg m-2 s-1)
%         aphoto (umol m-2 s-1)
% 
% --------------------------------------------------
% 
%         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(tk_25) exp(tk -298)eact/(298 R tk)), where eact is the
%         activation energy.
% 
%         Boltzmann distribution
% 
%         F(T)=tboltz)
% 
% 
%         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))
% 
%         Wj=J Ci/(4 Ci + 8 gamma)
% 
%         Ps kinetic coefficients from Harley at WBW.
% 
%         Gamma is the CO2 compensation point
% 
%      Information on the leaf photosynthetic parameters can be found in:
% 
%      Harley, P.C. and Baldocchi, 1995.Scaling carbon dioxide and water vapor exchange
%      from leaf to canopy in a deciduous forest:leaf level parameterization.
%      Plant, Cell and Environment. 18: 1146-1156.
% 
%      Wilson, K.B., D.D. Baldocchi and P.J. Hanson. 2000. Spatial and seasonal variability of
%      photosynthesis parameters and their relationship to leaf nitrogen in a deciduous forest.
%      Tree Physiology. 20, 565-587.
% 
% 
%      Tests of the model are reported in:
% 
%      Baldocchi, D.D. 1997. Measuring and modeling carbon dioxide and water vapor
%      exchange over a temperate broad-leaved forest during the 1995 summer drought.
%      Plant, Cell and Environment. 20: 1108-1122
% 
%      Baldocchi, D.D. and P.C. Harley. 1995. Scaling carbon dioxide and water vapor
%      exchange from leaf to canopy in a deciduous forest: model testing and application.
%      Plant, Cell and Environment. 18: 1157-1173.
% 
%      Baldocchi, D.D and T.P. Meyers. 1998. On using eco-physiological, micrometeorological
%      and biogeochemical theory to evaluate carbon dioxide, water vapor and gaseous deposition 
%      fluxes over vegetation. Agricultural and Forest Meteorology 90: 1-26.
%      
%      Baldocchi, D.D. Fuentes, J.D., Bowling, D.R, Turnipseed, A.A. Monson, R.K. 1999. Scaling
%      isoprene fluxes from leaves to canopies: test cases over a boreal aspen and a mixed species temperate
%      forest. J. Applied Meteorology. 38, 885-898.
% 
%      Baldocchi, D.D. and K.B.Wilson. 2001. Modeling CO2 and water vapor exchange of a
%      temperate broadleaved forest across hourly to decadal time scales. Ecological Modeling 
%           142: 155-184



% Debug...applied ideas of Collatz and Berry to compute a smooth transition
% between Wj and Wc...Seeing jump in the plots

%  May have units problems..Collatz states that many of the MM equations
%  are in units of partial pressure, pi, not mixing ratio, c.  Note that
%  some of the units of Ko and Kc are in units of Pa yet I am dividing by
%  units of mixing ratio..Think I have problems


      %  vcopt = 127;   % carboxylation rate at optimal temperature, umol m-2 s-1 
      %  jmopt = 277;   % electron transport rate at optimal temperature, umol m-2 s-1 
      %  rd25 = .5;     % dark respiration at 25 C, rd25= 0.34 umol m-2 s-1 
      
      jmopt=2.2*vcopt; % scale Jmax with vcmax like Wullschleger
        
        % Universal gas constant  

        param.rugc = 8.314;              % J mole-1 K-1 
        rgc1000 = 8314;            % gas constant times 1000.


%          Consts for Photosynthesis model and kinetic equations.
%          for Vcmax and Jmax.  Taken from Harley and Baldocchi (1995, PCE)
 

        param.hkin = 200000.0;    % enthalpy term, J mol-1
        param.skin = 710.0;       % entropy term, J K-1 mol-1
        ejm = 55000.0;      % activation energy for electron transport, J mol-1
        evc = 55000.0;      % activation energy for carboxylation, J mol-1
    

%         Enzyme constants & partial pressure of O2 and CO2
%         Michaelis-Menten K values. From survey of literature.


       kc25 = 274.6;   % kinetic coef for CO2 at 25 C, microbars  
        
         ko25 = 419.8;   % kinetic coef for O2 at 25C,  millibars 
         o2 = 210.0;     % oxygen concentration  mmol mol-1  
        
        
        

%         tau is computed on the basis of the Specificity factor (102.33) 
%         times Kco2/Kh2o (28.38) to convert for value of C in solution
%         to that based in air/
% 
%         New value for Quercus robor from Balaguer et al. 1996
%         Similar number from Dreyer et al. 2001, Tree Physiol, tau= 2710

        tau25 = 2904.12;    %  tau coefficient
           

%          Arrhenius constants
%          Eact for Michaelis-Menten const. for KC, KO and dark respiration
%          These values are from Harley


        ekc = 80500.0;     % Activation energy for K of CO2; J mol-1  
        eko = 14500.0;     % Activation energy for K of O2, J mol-1
        erd = 38000.0;     % activation energy for dark respiration, eg Q10=2 
        ektau = -29000.0;  % J mol-1 (Jordan and Ogren, 1984)
        tk_25 = 298.16;    % absolute temperature at 25 C
        toptvc = 311.0;    % optimum temperature for maximum carboxylation
        toptjm = 311.0;    % optimum temperature for maximum electron transport
            
        kball = 9.5;		%  Ball-Berry stomatal coefficient for stomatal conductance
        bprime = .0175;	    %  intercept of Ball-Berry model, mol m-2 s-1 
        bprime16 = 0.0109375;  % intercept for CO2, bprime16 = bprime / 1.6;
 
        rsm = 145.0;  % Minimum stomatal resistance, s m-1.
        brs=60.0;     % curvature coeffient for light response

        qalpha = .22; %   leaf quantum yield, electrons 
        qalpha2 = 0.0484;   % qalpha squared, qalpha2 = pow(qalpha, 2.0);

          
        
        rt = param.rugc * tlk;               %  product of universal gas constant and abs temperature

        tprime25 = tlk - tk_25;       % temperature difference

        ttemp = exp((param.skin * tlk - param.hkin) / rt) + 1.0; %  denominator term


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


        [kct] = TEMP_FUNC(kc25, ekc, tprime25, tk_25, tlk);
        [ko] = TEMP_FUNC(ko25, eko, tprime25, tk_25,tlk);
        [tau] = TEMP_FUNC(tau25, ektau, tprime25, tk_25,tlk);

        % bc = kct * (1.0 + o2 / ko);
        
        % fix the Ko and O2 values with same units
        
        ko25_Pa= ko25* 100;  % Pa
        o2_Pa= o2 * 101.3 ;  % Pa
        
        bc= kct * (1. + o2_Pa/ko25_Pa);

       

%         gammac is the CO2 compensation point due to photorespiration, umol mol-1
%         Recalculate gammac with the new temperature dependent KO and KC
%         coefficients..C at Vc = 0.5 Vo
%         gammac = O2/(2 tau)
%         O2 has units of kPa so multiplying 0.5 by 1000 gives units of Pa


          gammac = 500.0 * o2 / tau;   % umol/mol
  
     
%         // Scale rd with vcmax and apply temperature
%         // correction for dark respiration

        rdz= vcopt * 0.004657;


%         // reduce respiration by 40% in light according to Amthor 


        if(Iphoton > 10)
        rdz = rdz*0.4;   
        end
        
        % temperature correction to dark respiration

        [rd] = TEMP_FUNC(rdz, erd, tprime25, tk_25, tlk);

    
%        Apply temperature correction to JMAX and vcmax
    

        jmax = TBOLTZ(jmopt, ejm, toptjm, tlk);
        vcmax = TBOLTZ(vcopt, evc, toptvc, tlk);


%         T leaf boundary layer resistance
        
%         gb_mole leaf boundary layer conductance for CO2 exchange,
%         mol m-2 s-1         
%  
%         RB has units of s/m, convert to mol-1 m2 s1 to be
%         consistant with R.

      
        rb_mole = rb_co2 * tlk * 101.3* 0.022624/(273.15 * P_kPa);

        gb_mole = 1. / rb_mole;

        dd = gammac;
        b8_dd = 8 * dd;


% ***************************************
% 
%         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
% */

        % [rh_leaf]  = SFC_VPD(tlk, leleaf);
        
        rh_leaf=eair_Pa/ES(tlk);

        k_rh = rh_leaf * kball;  %// combine product of rh and K ball-berry


%         Gs from Ball-Berry is for water vapor.  It must be divided
%         by the ratio of the molecular diffusivities to be valid
%         for A

        k_rh = k_rh / 1.6;      % adjust the coefficient for the diffusion of CO2 rather than H2O

        gb_k_rh = gb_mole * k_rh;

        ci_guess = cca * .7;   % // initial guess of internal CO2 to estimate Wc and Wj

      %        // cubic coefficients that are only dependent on CO2 levels


        alpha_ps = 1.0 + (bprime16 / gb_mole) - k_rh;
        bbeta = cca * (gb_k_rh - 2.0 * bprime16 - gb_mole);
        gamma = cca * cca * gb_mole * bprime16;
        theta_ps = gb_k_rh - 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.
% 
%         estimate of J according to Farquhar and von Cammerer (1981)
% 
% 
%         J photon from Harley

%        looking at Collatz these W functions should be in terms of partial
%        pressure, units of Pa
% */

              
        if (jmax > 0)
        j_photon = qalpha * Iphoton / sqrt(1. +(qalpha2 * Iphoton * Iphoton / (jmax * jmax)));
        else
        j_photon = 0;
        end


        wj = j_photon * (ci_guess - dd) / (4. * ci_guess + b8_dd);


        wc = vcmax * (ci_guess - dd) / (ci_guess + bc);

       
        
        if(wj < wc)
        

%         // for Harley and Farquhar type model for Wj

        psguess=wj;
        
        B_ps = b8_dd;
        a_ps = j_photon;
        E_ps = 4.0;
        
        else
        
        psguess=wc;
          
        B_ps = bc;
        a_ps = vcmax;
        E_ps = 1.0;
        
        end


%         if wj or wc are less than rd then A would probably be less than zero.  This would yield a
%         negative stomatal conductance.  In this case, assume gs equals the cuticular value. This
%         assumptions yields a quadratic rather than cubic solution for A


        
%         /*
%         cubic solution:
% 
%          A^3 + p A^2 + q A + r = 0
%         */

        denom = E_ps * alpha_ps;

        Pcube = (E_ps * bbeta + B_ps * theta_ps - a_ps * alpha_ps + E_ps * rd * alpha_ps);
        Pcube = Pcube/denom;

        Qcube = (E_ps * gamma + (B_ps * gamma / cca) - a_ps * bbeta + a_ps * dd * theta_ps + E_ps * rd * bbeta + rd * B_ps * theta_ps);
        Qcube = Qcube/denom;

        Rcube = (-a_ps * gamma + a_ps * dd * (gamma / cca) + E_ps * rd * gamma + rd * B_ps * gamma / cca);
        Rcube = Rcube/denom;


%         // Use solution from Numerical Recipes from Press


        P2 = Pcube * Pcube;
        P3 = P2 * Pcube;
        Q = (P2 - 3.0 * Qcube) / 9.0;
        R = (2.0 * P3 - 9.0 * Pcube * Qcube + 27.0 * Rcube) / 54.0;


% /*
%         Test = Q ^ 3 - R ^ 2
%         if test >= O then all roots are real
% */

        rr=R*R;
        qqq=Q*Q*Q;
        
%                 // real roots 

        
                        arg_U = R / sqrt(qqq);

                        ang_L = acos(arg_U);
        
                        root1 = -2.0 * sqrt(Q) * cos(ang_L / 3.0) - Pcube / 3.0;  
                        root2 = -2.0 * sqrt(Q) * cos((ang_L + pi*2) / 3.0) - Pcube / 3.0;
                        root3 = -2.0 * sqrt(Q) * cos((ang_L - pi*2) / 3.0) - Pcube / 3.0;

%                 // rank roots #1,#2 and #3 according to the minimum, intermediate and maximum
%                 // value


                        if(root1 <= root2 && root1 <= root3)
                        minroot=root1;
                                if (root2 <= root3)
                                midroot=root2;
                                maxroot=root3;
                                 else
                                 midroot=root3;
                                 maxroot=root2;
                                end
                        end


                        if(root2 <= root1 && root2 <= root3)
                        minroot=root2;
                                if (root1 <= root3)
                                midroot=root1;
                                maxroot=root3;
                                 else
                                 midroot=root3;
                                 maxroot=root1;
                                 end
                        end


                        if(root3 <= root1 && root3 <= root2)
                                minroot=root3;
                                if (root1 < root2)
                                midroot=root1;
                                maxroot=root2;
                                 else
                                midroot=root2;
                                maxroot=root1;
                                end

                        end  %// end of the loop for real roots


       % // find out where roots plop down relative to the x-y axis
            
        
        if (minroot > 0 && midroot > 0 && maxroot > 0)
        aphoto=minroot;
        end
        
        
        if (minroot < 0 && midroot < 0 && maxroot > 0)
        aphoto=maxroot;
        end

   
        if (minroot < 0 && midroot > 0 && maxroot > 0)
        aphoto=midroot;
        end

% /*
%          Here A = x - p / 3, allowing the cubic expression to be expressed
%          as: x^3 + ax + b = 0
% */

%                 // aphoto=root3;  // back to original assumption      
% 
% /*
%         also test for sucrose limitation of photosynthesis, as suggested by
%         Collatz.  Js=Vmax/2
% */
        j_sucrose = vcmax / 2.;
        
       

        
                
       
       
        
        
        if (wj <=rd) || (wc <= rd)
            
            gs_leaf_mole = bprime;
        gs_co2 = gs_leaf_mole / 1.6;

%  /*
%        stomatal conductance is mol m-2 s-1
%         convert back to resistance (s/m) for energy balance routine
% */

        gs_m_s = gs_leaf_mole * tlk * 101.3*0.022642/(P_kPa*273.15) ;
        
        



%         a quadratic solution of A is derived if gs=ax, but a cubic form occurs
%         if gs =ax + b.  Use quadratic case when A is less than zero because gs will be
%         negative, which is nonsense

            
         ps_1 = cca * gb_mole * gs_co2;
         delta_1 = gs_co2 + gb_mole;
         denom = gb_mole * gs_co2;

         Aquad1 = delta_1 * E_ps;
         Bquad1 = -ps_1 * E_ps - a_ps * delta_1 + E_ps * rd * delta_1 - B_ps * denom;
         Cquad1 = a_ps * ps_1 - a_ps * dd * denom - E_ps * rd * ps_1 - rd * B_ps * denom;

         product=Bquad1 * Bquad1 - 4.0 * Aquad1 * Cquad1;

        if (product >= 0)
        sqrprod= sqrt(product);
        end
         
         aphoto = (-Bquad1 - sqrprod) / (2.0 * 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 = cca - aphoto / gb_mole;
        ci = cs - aphoto / gs_co2;
       
            
        end
        
       
        

        cs = cca - aphoto / gb_mole;

        if(cs > 1000)
        cs=input.co2air;
        end

% /*
%         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_leaf_mole = (kball * rh_leaf * aphoto / cs) + bprime;


%         // convert Gs from vapor to CO2 diffusion coefficient


        gs_co2 = gs_leaf_mole / 1.6;

% /*
%         stomatal conductance is mol m-2 s-1
%         convert back to resistance (s/m) for energy balance routine
% */

        gs_m_s = gs_leaf_mole * tlk * 101.3* 0.022624/(273.15 * P_kPa);
        
        ps.rstom=1/gs_m_s;
        
     
        
       

%         // to compute ci, Gs must be in terms for CO2 transfer


        ci = cs - aphoto / gs_co2;
        
         % recompute wj and wc with ci
        
        wj = j_photon * (ci - dd) / (4. * ci_guess + b8_dd);

        wc = vcmax * (ci - dd) / (ci_guess + bc);
        
         % Collatz uses a quadratic model to compute a dummy variable wp to allow
  % for the transition between wj and wc, when there is colimitation.  this
  % is important because if one looks at the light response curves of the
  % current code one see jumps in A at certain Par values
  
  % theta wp^2 - wp (wj + wc) + wj wc = 0
  
  % a x^2 + b x + c = 0
  
  % x = [-b +/- sqrt(b^2 - 4 a c)]/2a

       a=0.98;
       b= -(wj +wc);
       c=wj*wc;
       
       wp1=(-b + sqrt(b*b - 4*a*c))/(2*a);
       wp2=(-b - sqrt(b*b - 4*a*c))/(2*a);
       
       wp=min(wp1,wp2);

% beta A^2 - A (Jp+Js) + JpJs = 0

       aa = 0.95;
       bb= -(wp+ j_sucrose);
       cc = wp* j_sucrose;
       
       
       Aps1=(-bb + sqrt(bb*bb - 4*aa*cc))/(2*aa);
       Aps2=(-bb - sqrt(bb*bb - 4*aa*cc))/(2*aa);
       
       Aps=min(Aps1,Aps2);
       
       if(Aps < aphoto && Aps > 0)
           aphoto=Aps - rd;
       end
       

       
         ps.aphoto=aphoto;
         ps.ci=ci;
         ps.gs_co2=gs_co2;
         ps.wj=wj;
         ps.wc=wc;
         ps.wp=wp;
         ps.jsucrose=j_sucrose;
         

        

function [fT]=DESDT(T)

% Slope of the saturation vapor pressure-temperature curve
% 	
% Thermodynamic relationship from HESS    (Pa K-1)
param.Rstar = 8.3144; 

fT = ES(T) * LAMBDA(T) * 18. / (param.Rstar * T * T * 1000.);
 
%///////////////////////////////////////////////////////////////
function [fT]=ES(T) 
% 	Saturation vapor pressure,  Pa

fT = 100. * exp(54.8781919 - 6790.4985 / T - 5.02808 * log(T));


%//////////////////////////////////////////////
function [fT]=LAMBDA(T)

% Latent heat of vaporization, J kg-1

fT = 3149000 - 2370 * T;
        
%//////////////////////////////////////////////

    function [fT] =TBOLTZ (rate, eakin, topt, tl)
        
        % Boltzmann temperature distribution for photosynthesis 
        
        param.Rstar = 8.3144;
        param.hkin = 200000.0;    % enthalpy term, J mol-1
        param.skin = 710.0;       % entropy term, J K-1 mol-1

        dtlopt = tl - topt;
        prodt = param.Rstar * topt * tl;
        numm = rate * param.hkin * exp(eakin * (dtlopt) / (prodt));
        denom = param.hkin - eakin * (1.0 - exp(param.hkin * (dtlopt) / (prodt)));
        fT = numm / denom;
        
  %/////////////////////////////////////////////////////////////      
        function [fT] =TEMP_FUNC(rate,eact,tprime,tref,t_lk)

        %Arhennius temperature function  
        
        param.Rstar = 8.3144;
            
        fT = rate * exp(tprime * eact / (tref * param.Rstar*t_lk));
        
    %//////////////////////////////////////////////////////////
%            function[fx] =SFC_VPD (tlk,leleaf)

%              this function computes the relative humidity at the leaf surface for
%             application in the Ball Berry Equation
% 
%             latent heat flux, LE, is passed through the function, mol m-2 s-1
%             and it solves for the humidity at leaf surface

%            [es_leaf] = ES(tlk);    % saturation vapor pressure at leaf temperature
%            
%            rhov_sfc = (leleaf / LAMBDA(tlk)) * bound_layer_res.vapor + met.rhova_kg; % kg m-3 */
% 
%         e_sfc = rhov_sfc * tlk / .2165;        % mb 
%         vpd_sfc = es_leaf - e_sfc;             % mb 
%         rhum_leaf = 1. - vpd_sfc / es_leaf;    % 0 to 1.0 
%         fx = rhum_leaf;
         
       
        
		
		
		
               
       
  
       




      

                
          
