* microe.quartic.f    rhino    ~/gx   14 April 2000
*
* Program to compute Ci and A in any microenvir. (especially the
*  original, pre-cuvette conditions), given those conditions, plus
*  gs, Tleaf, Vcmax25, Rd25, theta (convexity), and Q00
* - by solving Johnson-Thornley equation, theta*A**2 + (Amax+Q0*IL)*A
*  +Amax*Q0*IL = 0, with A=Agross in the above = (Ca-Ci)*gt/P + Rd
*  and Amax = Vcmax*(Ci - gamma)/(Ci + Kco) and Q0 = Q00*(Ci-gamma)/
*  (Ci + b*gamma).  On multiplying through by the common denominator,
*  (Ci + b*gamma)*(Ci + Kco), we get a quartic equation in Ci to solve
*  and then we substitute this Ci to get A.
*
* Shorthand notation: T = theta, K = Kco, G = gamma, gt = g(tot,CO2),
*  V = Vcmax(Tleaf), Rd=Rd(Tleaf), P = Pair 
*
      real IL, ILumol, Kc, Ko, K, IBB
      common c4,c3,c2,c1,c0
      b=2.
      Rgas=8.314
      write(6,'("Enter leaf physiological parameters and temperature:",
     + /,3x,"gs (mol m-2 s-1), Tleaf (oC), Vcmax25 (umol m-2 s-1), ",
     + "Rd25 (same),",/,3x,"theta (unitless), Q00 (mol CO2/mol ",
     + "photons), dleaf (m)")')
      read(5,*)gs,Tleaf,Vcmax25umol,Rd25umol,T,Q00,dleaf
      write(6,'("Enter environmental conditions:",/,5x,"Pair (Pa), ",
     + "Ca (Pa), IL (umol m-2 s-1), u (m s-1), eair(Pa)")')
      read(5,*)P, Ca, ILumol, u, eair
* Convert to standard SI units (umol -> mol, ...)
   50 Vcmax25=1.e-6*Vcmax25umol
      IL=1.e-6*ILumol
      Rd25=1.e-6*Rd25umol
* Use Tleaf to compute Kco=K, gamma=G
      dT=Tleaf-25.
      Tabs=Tleaf+273.2
      G=3.69+dT*(0.188 + 0.0036*dT)
      fac=dT/(Rgas*298.2*Tabs)
      Kc=40.4*exp(59400.*fac)
      Ko=24800.*exp(36000.*fac)
      K=Kc*(1.+P*0.21/Ko)
      V=Vcmax25*exp(64000.*fac)
      Rd=Rd25*exp(0.07*dT)
* Compute gb from u, dleaf and then compute gt
      gb=0.264*sqrt(u/dleaf)
      gt=1./(1.6/gs+1.37/gb)
* Now compute terms in quartic
      gt2=gt*gt
      P2=P*P
      QI=Q00*IL
      c4=T*gt2/P2
      c3=T*gt2*(K+b*G-2.*Ca)/P2 -2.*T*gt*Rd/P +V*gt/P + QI*gt/P
      c2=T*gt2*(K*b*G -2.*Ca*(K+b*G)+Ca*Ca)/P2
     + +2.*T*gt*Rd*(Ca-K-b*G)/P +T*Rd*Rd -V*gt*(Ca-(b-1.)*G)/P
     + -V*Rd -QI*gt*(Ca-K+G)/P -QI*Rd +V*QI
      c1=T*gt2*(Ca*Ca*(K+b*G) - 2.*Ca*K*b*G) /P2
     + +2.*T*gt*Rd*(Ca*(K+b*G)-K*b*G)/P +T*Rd*Rd*(K+b*G)
     + -V*gt*(b*G*G+Ca*(b-1.)*G)/P -V*Rd*(b-1.)*G 
     + -QI*gt*(G*K+Ca*(K-G))/P -QI*Rd*(K-G) -V*QI*2.*G
      c0=T*gt2*Ca*Ca*K*b*G/P2 + 2.*T*gt*Rd*Ca*K*b*G/P
     + +T*Rd*Rd*K*b*G +V*gt*Ca*b*G*G/P +V*Rd*b*G*G +QI*gt*G*K*Ca/P
     + +QI*Rd*G*K +V*QI*G*G
* Solve by binary search
      nmax=20
      n=0
      Cilo=G
      Cihi=1.5*Ca
      Flo=F(Cilo)
      Fhi=F(Cihi)
  100 if(n.eq.0)write(6,'("F(Cilo)=",e14.6,", F(Cihi)=",e14.6)')Flo,Fhi
      if(Flo*Fhi.gt.0.)then
        write(6,'("No root, or 2 roots, in interval; stop")')
        write(6,'("Enter new guesses for Cilo, Cihi (Pa)")')
        read(5,*)Cilo,Cihi
        Flo=F(Cilo)
        Fhi=F(Cihi)
        goto 100
      elseif(n.le.nmax)then
        Cimid=0.5*(Cilo+Cihi)
        n=n+1
        Fmid=F(Cimid)
        if(Fmid*Flo.lt.0.)then
          Cihi=Cimid
          Fhi=Fmid
        else
          Cilo=Cimid
          Flo=Fmid
        endif
        goto 100
      endif
* Finished iterations
      Anet=(Ca-Cimid)*gt/P
        Amax=V*(Cimid-G)/(Cimid+K)
        ALL=Q00*IL*(Cimid-G)/(Cimid+b*G)
        Ag=0.5*(Amax+ALL-sqrt((Amax+ALL)**2-4*T*Amax*ALL))/T
        write(6,'("Agross,direct calc., is ",f7.2)')1.e6*Ag
      Agross=Anet+Rd
      write(6,'("Ci=",f7.2," Pa; A(net)=",6Pf7.2," umol m-2 s-1; ",
     + "A(gross)=",6Pf7.2)')Cimid,Anet,Agross
      hs=(eair/esat(Tleaf)+gs/gb)/(1.+gs/gb)
      Cs=Ca-Anet*P*1.37/gb
      IBB=Anet*hs*P/Cs
      write(6,'(3x,"hs=",0Pf7.4,", Cs=",f6.2," Pa, IBB=",f7.4)')
     + hs,Cs,IBB
      write(6,'(3x,"A(gross)=",6Pf7.2,", Vcmax=",f7.2,", Cifac=",
     + 0Pf7.4)')Agross,V,(Cimid-G)/(Ci+K)
      write(6,'("Want to rerun with changed IL? Enter IL>0 (<0 means",
     + " stop)")')
      read(5,*)ILumol
      if(ILumol.gt.0.)goto 50
      stop
      end      
*
*
      real function F(Ci)
      common c4,c3,c2,c1,c0
      F=c0+Ci*(c1+Ci*(c2+Ci*(c3+Ci*c4)))
      return
      end
*
*
      real function esat(T)
      esat=610.8*exp(17.269*T/(237.2+T))
      return
      end