*BBbinary.JT.options.f   rhino   ~/nigec/sap    23 Feb. 02 - adding
*   full J-T form, etc. from nigec/bosqueE6...f
* Program to find the value of gs satisfying BB, carbox., and
*  energy bal. equations, by binary search.
*
* 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)
*
      real mBB,Kc25,Ko25
      logical BBswitch
      character*20 chaval
      character*3 test, ID(20)
      character*17 sval
      dimension nread(20)
      common/com0/mBB,bBB,BBswitch,ha,tol
      common/com1/Tair,sigma,Hvap,gtotw,eair,Pair,Bcc,Qrad,BTIR,E,eL,
     + gbair
      common/com2/TL,Vcmax25,Kc25,Kc,Ko25,Ko,gamma,Ca,A,Q00,PPFD,
     + Cs,gtotCO2,Ci,thetaPS,RdTmean,Tmean
      common/com3/gsplus,gsminus,hs0,hs
      data ID/'vcm','the','rdo','tme','apa','mbb','bbb','hs0',
     + 'tol','tai','eai','tve','tsk','ppf','gba','cai','pai',3*''/
      data nread/1,1,1,1,1,1,2,2,2,2,3,3,3,3,6*0/
      ncha=0
*       means that this is not a re-run
*  Read in leaf parameters
  110 write(6,1)
    1 format(1x,'Enter Vcmax25(umol...)    aPAR    mBB    bBB    hs0',
     + '    tol',/,5x,'(hs<0 means use true BB response;', 
     + ' tol=precis. of gs sol.')
      read(5,*)Vcmax25u,aPAR,mBB,bBB,hs0,tol !Vcmax25u is in umol...
* Read in more leaf parameters
   11 write(6,'("Enter thetaPS (0-1), RdovrA (umol...), Tmean")')
      read(5,*)thetaPS,RdovrA,Tmean
* Read in micro-e
      write(6,'("Enter Tair(oC), eair(Pa, or as -RHfraction),",
     + "Tveg(oC), Tsky(oC)")')
      read(5,*)Tair,eairi,Tveg,Tsky !eairi is not yet assured as Pascals
      write(6,'(1x,"Enter PPFD(umol...)    gbair(mol...)    Ca(Pa)",
     +  "    Pair (Pa)")')
      read(5,*)PPFDu,gbair,Ca,Pair  ! PPFDu is PPFD in umol...
   85 write(6,'("To modify parameters, enter the keyword and ",
     + "the value, such as ''vcm 74.5''",/,5x,"To see the ",
     + "keywords and the other commands, type ''key''")')
   87 write(6,'("Enter the keyword and value (total up to ",
     + "20 characters)")')
      read(5,'(a20)',err=88)chaval
      goto 89  ! If there's not an error, continue the treatment
   88 write(6,'("Error; try again")')
      goto 87
   89 test=chaval(1:3)
      ncha=1
      if(test.eq.'key')then
        write(6,888)
  888   format(1x,"Keywords are:",/,5x,"vcm -> Vcmax25; tps-> thetaPS",
     +   ";rda - > Rd/A; ; tme -> Tmean for Rd",/,5x,"apa -> aPAR; ",
     +   " mbb -> mBB; bbb -> bBB; hs0 -> hs0;",/,5x,"tol -> tol; ",
     +   "tai->Tair; eai->eair; tve->Tveg; tsk->Tsky, ppf->PPFD",/,5x,
     +   "gba->gbair; ca_ -> Ca; pai -> Pair",//,5x,"sto -> STOP; ",
     +   "see -> print all values; all -> change all values",/,10x,
     +   "(nothing, except a carriage return) -> run with current ",
     +   "values",/)  ! Changed 11 Oct.
      elseif(test.eq.'all')then
        ncha=0
        goto 110
      elseif(test.eq.'sto')then 
        stop
      elseif(test.eq.'see')then
        write(6,25)Vcmax25u,thetaPS,RdovrA,Tmean,aPAR,mBB,bBB,hs0,tol,
     +   Tair,eairi,Tveg,Tsky,PPFDu,gbair,Ca,Pair  
   25   format(3x,"Vcmax25=",f7.2," umol m-2 s-1; theta(PS)=",f6.3,
     +   "; Rd/A=",f6.4,"; Tmean for Rd=",f6.2," oC",/,3x,"aPAR=",f6.4,
     +   "; mBB=",f6.2,"; bBB=",f7.4," mol m-2 s-1;",/,3x,
     +   "hs0=",f7.3,"; tol=",f6.4," mol m-2 s-1; Tair=",f6.2,/,3x,
     +   "eair=",f7.2,"; Tveg=",f6.2," oC; Tsky=",f6.2," oC; PPFD=",
     +   f6.1,/,3x,"gbair=",f6.3," mol m-2 s-1; Ca=",f6.2,
     +   " Pa; Pair=",f8.1," Pa",/) 
        goto 85
      elseif(test.eq.'')then
        goto 90
      else
        sval=chaval(4:)
        rewind 8
        write(8,*)sval
        rewind 8
        read(8,*,err=88)val
        if(test.eq.'vcm')then
          Vcmax25u=val
          write(6,'("***New Vcmax25=",f6.2)')val
        elseif(test.eq.'tps')then  ! Changed 11 Oct.
          thetaPS=val
          write(6,'("***New thetaPS=",f6.2)')val
        elseif(test.eq.'rda')then
          RdovrA=val
          write(6,'("***New RdovrA=",f6.4)')val
        elseif(test.eq.'tme')then
          Tmean=val 
          write(6,'("***New Tmean=",f6.2)')val
        elseif(test.eq.'apa')then
          aPAR=val
          write(6,'("***New aPAR=",f6.4)')val
        elseif(test.eq.'mbb')then
          mBB=val
          write(6,'("***New mBB=",f6.2)')val
        elseif(test.eq.'bbb')then
          bBB=val
          write(6,'("***New bBB=",f6.4)')val
        elseif(test.eq.'hs0')then
          hs0=val
          write(6,'("***New hs0=",f6.3)')val
        elseif(test.eq.'tol')then
          tol=val
          write(6,'("***New tol=",f8.5)')val
        elseif(test.eq.'tai')then
          Tair=val
          write(6,'("***New Tair=",f6.2)')val
        elseif(test.eq.'eai')then
          eairi=val
          write(6,'("***New eair=",f8.2)')val
        elseif(test.eq.'tve')then
          Tveg=val
          write(6,'("***New Tveg=",f6.2)')val
        elseif(test.eq.'tsk')then
          Tsky=val
          write(6,'("***New Tsky=",f6.2)')val
        elseif(test.eq.'ppf')then
          PPFDu=val
          write(6,'("***New PPFD=",f8.2)')val
        elseif(test.eq.'gba')then
          gbair=val
          write(6,'("***New gbair=",f7.5)')val
        elseif(test.eq.'ca_')then
          Ca=val
          write(6,'("***New Ca=",f6.2)')val
        elseif(test.eq.'pai')then
          Pair=val
          write(6,'("***New Pair=",f8.1)')val
        else
          write(6,'("Keyword unrecognizable; try again")')
        endif
        goto 85  ! Allow further changes
      endif
* Convert to SI
   90 Vcmax25=1.e-6*Vcmax25u
      PPFD=1.e-6*PPFDu
      esatair=esat(Tair)
      if(eairi.lt.0.)then
        eair=-eairi*esatair
      else
        eair=eairi
      endif
      ha=eair/esatair
* Convert other entries that have switching functions
      BBswitch=0
      if(hs0.lt.0.) BBswitch=1
      aNIR=0.35
* Set some constants
      sigma=5.67e-8
      epsleaf=0.96
*       (This is not used yet)
      Hvap=4.5e4
      Cpair=29.
      Kc25=40.4
      Ko25=24800.
      Q00=0.0732
* Calc. some derived quantities
      Bcc0=0.85*Cpair*gbair
      QSW=PPFD*2.2e5*(aPAR+aNIR)
      BTIR0=8.*sigma*(Tair+273.2)**3
* Calculate Rd at Tmean; 1st, calc. A at Tmean, standard Ca, 50% RH,
*   PPFD=1000 umol...,Tveg=Tmean, Tsky=Tmean-25, and u=1 m s-1:
* First, save original micro-e conditions
      Tair0=Tair
      ha0=ha
      eair0=eair
      gbair0=gbair
      PPFD0=PPFD
      Tair=Tmean
      ha=0.5
      eair=esat(Tmean)*0.5
      gbair=1. ! or 0.264*sqrt(1./dleaf), if we go back to using u, dleaf
      Bcc=0.85*Cpair*gbair
      PPFD=1.e-3
      Qsky=sigma*(Tair+243.2)**4
      Qveg=sigma*(Tair+273.2)**4
      QSW=PPFD*2.2e5*(aPAR+aNIR)
      BTIR=8.*sigma*(Tair+273.2)**3
      Qrad=Qsky+Qveg+QSW
*       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=Qsky+Qveg+QSW
      RdTmean=0.
      call BBsolve  ! This returns A through common
      RdTmean=RdovrA*A
      write(6,'("Rd(Tmean)=RdovrA*",f6.2)')1.e6*RdTmean
* Now reset to current micro-e and solve for A, gs, etc.
      Tair=Tair0
      ha=ha0
      eair=eair0
      gbair=gbair0
      PPFD=PPFD0
      Bcc=Bcc0
      Qsky=sigma*(Tsky+273.2)**4
      Qveg=sigma*(Tveg+273.2)**4
      QSW=PPFD*2.2e5*(aPAR+aNIR)
      Qrad=Qsky+Qveg+QSW
      call BBsolve
  200 write(6,120)gsminus,gsplus
  120 format(1x,'Exact solution lies between',2f16.6)
      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
      goto 85
      end
*
*
      subroutine BBsolve
      real mBB,IBBlo,IBBmid,IBBhi
      logical Bbswitch
      common/com0/mBB,bBB,BBswitch,ha,tol
      common/com1/Tair,sigma,Hvap,gtotw,eair,Pair,Bcc,Qrad,BTIR,E,eL,
     + gbair
      common/com2/TL,Vcmax25,Kc25,Kc,Ko25,Ko,gamma,Ca,A,Q00,PPFD,
     + Cs,gtotCO2,Ci,thetaPS,RdTmean,Tmean
      common/com3/gsplus,gsminus,hs0,hs
* Method: assume a gs (at two endpoints)
*   Calculate Tleaf, thus, E, (hs,) Vcmax, and Q00
*   Solve for Asat(Ci),TL)=gtotCO2*(Ca-Ci)/P
*   Solve for ALL(Ci)=Q00*(Ci-gamma)/(Ci+2*gamma)*IL=gtot*(Ca-Ci)/P
*   Choose lower of two
*   Test if gs >, <, or = mA/Cs
*
* Get a starting guess for gs: take Tleaf=Tair, solve for A; then
*   compute Cs and compose A*hs0/Cs to estimate gs; take 0.7 and
*   1.5x this as endpoints
      TL=Tair
*   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
      gs0=mBB*A*hs*Pair/Cs
*XXX
      gs0=max(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=mBB*A*hs*Pair/Cs+bBB
      testlo=gslo-IBBlo
* 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=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
          gsplus=gshi
          gsminus=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
        gsplus=gslo
        gsminus=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
          gsmin=amax1(bBB,0.003) 
*   Make sure gs does not become too small
          if(gshi.le.gsmin)then
            gsplus=gsmin
            gsminus=gsmin
            gtotCO2=1./(1.37/gbair+1.6/gsmid)
            TL=Tleaf(gsplus)
            call CiACssolve
            if(BBswitch)then
              hs=(gbair*eair/eL+gsmid)/(gbair+gsmid)
            else
              hs=hs0
            endif
            goto 200 ! Take this as final result
          endif
          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=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=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
          gsplus=gsmid
          gsminus=gsmid
          goto 200
        else
          gsplus=gsmid
        endif
      enddo
  200 return
      end
*
*
      function Tleaf(gs)
      common/com1/Tair,sigma,Hvap,gtotw,eair,Pair,Bcc,Qrad,BTIR,E,eL,
     + gbair
*! Do  I need TL in common, as in E6..f?
      Tleaf=Tair
      gtotw=1./(1./gs+1./gbair)
   10 Qmrad=2.*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
      if(abs(dTleaf).gt.0.01)go to 10
      return
      end
*
*
      subroutine CiACssolve
      real Kc25,Kc,Ko25,Ko,Kco,K,IL
      common/com2/TL,Vcmax25,Kc25,Kc,Ko25,Ko,gamma,Ca,A,Q00,PPFD,
     + Cs,gtotCO2,Ci,thetaPS,RdTmean,Tmean
      common/com1/Tair,sigma,Hvap,gtotw,eair,Pair,Bcc,Qrad,BTIR,E,eL,
     + gbair
      dfac=298.*8.314
      denom=dfac*(TL+273.2)
      Vcmax=Vcmax25*exp(64800.*(TL-25.)/denom)
      Kc=Kc25*exp(59400.*(TL-25.)/denom)
      Ko=Ko25*exp(36000.*(TL-25)/denom)
      Kco=Kc*(1.+0.21*Pair/Ko)  ! Changed from fixes 21kPa of O2
      gamma=3.69+(TL-25.)*(0.188+0.00036*(TL-25.))
      Rd=RdTmean*exp(0.07*(TL-Tmean)) 
*  Solve the Johnson-Thornley equation for A,
*     thetaPS*A**2 - A*(Amax+Q0*IL) - sqrt((Amax+Q0*IL)**2-4*thetaPS*
*       Amax*Q0*IL)) = 0, with
*       Amax=Vcmax*(Ci-gamma)/(Ci+Kco),
*       Q0=Q00*(Ci-gamma)/(Ci+2*gamma)
*  Write A as gtotCO2*(Ca-Ci).
*  After writing all terms explicitly in terms of Ci, we multiply
*    through by terms in the denominator, obtaining a quartic in Ci.
*  Solve for Ci by a binary search.
*  Substitute this value of Ci in A = gtotCO2*(Ca-Ci) to get A.
*  Solve trivially for Cs = Ca - 1.37*A/gbair.
*
* Now compute terms in quartic; first, convert long names to short names
*   (rather than converting names in well-debugged code and possibly
*   introducing errors)
      gt=gtotCO2
      P=Pair
      IL=PPFD
      K=Kco
      G=gamma
      V=Vcmax
      T=thetaPS
*
      gt2=gt*gt
      P2=P*P
      QI=Q00*IL
      b=2. ! Multiplier of gamma in formula for Q0 = Q00*(Ci-gamma)/(Ci*b*gamma)
      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 Newton-Raphson search YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY
      Ci=Ca
      iter=0
  100 F=c0 + Ci*(c1 + Ci*(c2 + Ci*(c3 + Ci*c4)))
      Fprime=c1 + Ci*(2.*c2 + Ci*(3.*c3 + 4.*Ci*c4))
      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/T)*(Amax+ALL-sqrt((Amax+ALL)**2-4.*T*
     + Amax*ALL)) - Rd ! ZZZZZ Lack of Rd here was corrected 15 May 01
      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