!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      subroutine leafanetmodel(ilimittype,aPPFDlf,templeaf,pco2i,po2i,
     &chlflphips2,gascon,resistwp25,resistch25,ha_gmeso,hd_gmeso,
     &sv_gmeso,vcmax25,ha_vcmax,hd_vcmax,sv_vcmax,fjmax25,ha_jmax,
     &hd_jmax,sv_jmax,tpu25,ha_tpu,hd_tpu,sv_tpu,alpha,rdlight25,
     &ha_darkresp,stargamma25,ha_stargamma,fkc25,ha_kc,fko25,ha_ko,
     &phifactor,thetafactor,betaPSII,assim_net,rdlight,iphotolimit,atp,
     &resistwp,resistch,stargamma,pco2c,realizedfjelect,assim_net_flu,
     &pco2c_flu)
      implicit none
!
!------------ Inputs -------------------
!ilimittype=1: Rubisco,RuBp and TPU limitations
!ilimittype=2: Rubisco and RuBp limitations only
!ilimittype=3: Rubisco and TPU limitations only
!ilimittype=4: RuBp and TPU limitations only
!ilimittype=5: Rubisco limitation only
!ilimittype=6: RuBp limitation only
!ilimittype=7: TPU limitation only
! aPPFDlf:       absorbed photosynthetic photon flux density by leaf (umol m-2 s-1)
! templeaf:      leaf temperature [K]
! resistwp25:    resistance to CO2 via cell walls and plasmalemma [umol-1m2sPa]
! resistch25:    resistance to CO2 via chloroplastic envelope [umol-1m2sPa]
! vcmax25:       maximum RuBP saturated rate of carboxylation at 25oC
!                    of leaf temperature [umol m-2 s-1]
! fjmax25:       maximum rate of electron transport at 25oC
!                    of leaf temperature [umol m-2 s-1]
! rdlight25:     Mitochondrial respiration rate in the light at 25oC
! tpu25:         tpu at 25oC, [umol m-2 s-1]      
! alpha:         The fraction of glycolate carbon not returned to the chloroplat 
! pco2i:         intercellular CO2 partial pressure (Pa).
! po2i:          intercellular O2 partial pressure (Pa, often taking the ambient value).
! chlflphips2: Photochemical efficiency of photosystem II measured with fluorescence, if provided (NA)
! gascon:        universal gas constant [JK-1 mol-1].
! ha_vcmax:  Ha in Vcmax temperature response function ~73637.0d0/1000.0d0
! hd_vcmax:  Hd in Vcmax temperature response function ~149252.0d0/1000.0d0
! sv_vcmax:  Sv in Vcmax temperature response function ~486.0d0/1000.0d0
! ha_jmax:   Ha in Jmax temperature response function ~50300.0d0/1000.0d0
! hd_jmax:   Hd in Jmax temperature response function ~152044.0d0/1000.0d0
! sv_jmax:   Sv in Jmax temperature response function ~495.0d0/1000.0d0
! ha_tpu:    ha in tpu temperature response function ~53100/1000
! hd_tpu:    hd in tpu temperature response function ~201800/1000
! sv_tpu:    Sv in tpu temperature response function ~650/1000
! ha_darkresp: Mitochondrial respiation parameter deltaha = 46.39d0
! stargamma25: the compensation point at 25oC [Pa]
! ha_stargamma: parameter in the temperature response function, typically around
!       37.83 kJmol-1
! fkc25:      the Michaelis constant for CO2 at 25oC
! ha_kc:      a CO2 Michaelis temp coefficient (~79.43 KJmol-1)
! fko25:      the Michaelis constant for O2 at 25oC
! ha_ko:      a O2 Michaelis temp coefficient (~36.38 KJmol-1)
! ha_gmeso:  The activation energy (Ha) in gmeso temperature response function ~49.6d0[kJmol-1]
! hd_gmeso:  The deactivation energy (Hd) in gmeso temperature response function ~437.4d0[kJmol-1]
! sv_gmeso:  The entropy term (Sv) in gmeso temperature response function ~1.4d0 [kJmol-1K-1]
! phifactor: modifies Bernacchi phiPSIImax
! thetafactor: modifies Bernacchi thetaPSII
! betaPSII:    fraction of absorbed light that reaches photosystem II (~0.5)
!!3/29/2014. Bernacchi method cannot fit A/PAR curves well. phifactor and thetafactor are
!used to modify his method when A/PAR curves are present.
      integer ilimittype
      double precision aPPFDlf,templeaf,vcmax25,fjmax25,rdlight25,
     &tpu25,alpha,pco2i,po2i,chlflphips2,gascon,ha_vcmax,hd_vcmax,
     &sv_vcmax,ha_jmax,hd_jmax,sv_jmax,ha_tpu,hd_tpu,sv_tpu,
     &ha_darkresp,stargamma25,ha_stargamma,fkc25,ha_kc,fko25,ha_ko,
     &resistwp25,resistch25,ha_gmeso,hd_gmeso,sv_gmeso,phifactor,
     &thetafactor,betaPSII
!------------ Outputs -------------------
! assim_net:     net rate of CO2 uptake per unit leaf area
!                     [umol m-2 s-1]
! rdlight:       Mitochondrial respiration rate in the light
!                     [umol m-2 s-1]
! iphotolimit:   index of photosynthesis limiting process
!                    1 = Rubisco-limited photosynthesis
!                    2 = RuBP-Limited photosynthesis
!                    3 = TPU-limited photosynthesis
! stargamma:     CO2 compensation point (Pa)      

! atp:           the concentration of ATP (mmol ATP m-2)
! resistwp:    resistance to CO2 via cell walls and plasmalemma [umol-1m2sPa]
! resistch:    resistance to CO2 via chloroplastic envelope [umol-1m2sPa]
! realizedfjelect: the realized electron transport rate, <=fjelect from light equation (=when
!                  RuBP regeneration limits photosynthesis) (umolm-2s-1)
! assim_net_flu: assim_net calculated with chlflphips2 (if provided), umolm-2s-1
! pco2c_flu: pco2c calculated from chlflphips2 (if provided), Pa.
      double precision assim_net,rdlight,atp,pco2c,realizedfjelect,
     &assim_net_flu,pco2c_flu
      integer iphotolimit
!------------ Local variables -----------
! wc:           Rubisco limited rate of carboxylation [umol m-2 s-1] 
! wj:           RuBP limited rate of carboxylation [umol m-2 s-1]
! tpu:          rate of triose phosphate utilization [umol m-2 s-1]
! atpcon:       the concentration of ATP (mmol ATP m-2)
! fkco
! fkc:           Michaelis constant for CO2 (Pa)  
! fko:           Michaelis constant for O2 (Pa)
! vcmax:         maximum RuBP saturated rate of carboxylation [umol m-2 s-1] 
! fjelect:       rate of electron transport [umol m-2 s-1] 
      double precision wc,wj,tpu,atpcon,fkco,gmeso,resistwp,resistch,
     &fkc,fko,vcmax,fjelect,stargamma,term
      integer i
      call vcmaxontemp(templeaf,vcmax25,gascon,ha_vcmax,hd_vcmax,
     &sv_vcmax,vcmax)
      call jontemp(aPPFDlf,templeaf,fjelect,fjmax25,ha_jmax,hd_jmax,
     &sv_jmax,gascon,phifactor,thetafactor,betaPSII)
      call tpuontemp(templeaf,gascon,tpu25,ha_tpu,hd_tpu,sv_tpu,tpu)
      if(resistwp25.gt.0.0d0)then
        call gmesoontemp(templeaf,1.0d0,gascon,ha_gmeso,hd_gmeso,
     &sv_gmeso,gmeso)
        resistwp=resistwp25/gmeso
      else
        resistwp=0.0d0
      endif
      if(resistch25.gt.0.0d0)then
        call gmesoontemp(templeaf,1.0d0,gascon,ha_gmeso,hd_gmeso,
     &sv_gmeso,gmeso)
        resistch=resistch25/gmeso
      else
        resistch=0.0d0
      endif
      call resp_mitocho(templeaf,rdlight25,ha_darkresp,gascon,rdlight)
      call co2compens(templeaf,stargamma25,ha_stargamma,gascon,
     &stargamma)
      call MichaelisCO2(templeaf,fkc25,ha_kc,gascon,fkc)
      call MichaelisO2(templeaf,fko25,ha_ko,gascon,fko)
      fkco=fkc*(1.0d0+po2i/fko)
      call Anet_Final(vcmax,fjelect,tpu,resistwp,resistch,stargamma,
     &fkco,pco2i,alpha,rdlight,ilimittype,iphotolimit,assim_net,
     &pco2c,realizedfjelect)
      wc=vcmax*pco2c/(pco2c+fkco)
      wj=fjelect*pco2c/(4.0d0*pco2c+8.0d0*stargamma)
      atp=atpcon(wc,wj,vcmax)
      assim_net_flu=-9999.0d0
      pco2c_flu=-9999.0d0
      if(chlflphips2.gt.0.0d0)then
        fjelect=betaPSII*chlflphips2*aPPFDlf
        call Anet_Final(vcmax,fjelect,tpu,resistwp,resistch,stargamma,
     &fkco,pco2i,alpha,rdlight,6,i,assim_net_flu,pco2c_flu,term)
      endif
      return
      end subroutine leafanetmodel
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      double precision function threeexp(fact1,fact2,fact3)
      implicit none
!
! This function computes
!     threeexp=(dexp(fact3)+dexp(fact2+fact3))/(1.0d0+dexp(fact1))
! in a way that prevents overflow
!
      double precision fact1,fact2,fact3,term1,term2
 
      if(fact3.gt.fact1)then
        if(fact3.gt.-10.0d0)then
          if((dexp(-fact3)+dexp(fact1-fact3)).lt.1.0d-20)then
            term1=1.0d+30
          else
            term1=1.0d0/(dexp(-fact3)+dexp(fact1-fact3))
          endif
        else
          term1=dexp(fact3)/(1.0d0+dexp(fact1))
        endif
      else
        if(fact1.gt.-10.0d0)then
          term1=dexp(fact3-fact1)/(dexp(-fact1)+1.0d0)
        else
          term1=dexp(fact3)/(1.0d0+dexp(fact1))
        endif
      endif
      
      if((fact2+fact3).gt.fact1)then
        if((fact2+fact3).gt.-10.0d0)then
          if((dexp(-fact2-fact3)+dexp(fact1-fact2-fact3))
     &        .lt.1.0d-20)then
            term2=1.0d+30
          else
            term2=1.0d0/
     &            (dexp(-fact2-fact3)+dexp(fact1-fact2-fact3))
          endif
        else
          term2=dexp(fact2+fact3)/(1.0d0+dexp(fact1))
        endif
      else
        if(fact1.gt.-10.0d0)then
          term2=dexp(fact2+fact3-fact1)/(dexp(-fact1)+1.0d0)
        else
          term2=dexp(fact2+fact3)/(1.0d0+dexp(fact1))
        endif
      endif
      threeexp=term1+term2
      return
      end      
     
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      subroutine co2compens(templeaf,gammas25,ha_gammas,
     &       gascon,gammas)
      implicit none
      
! purpose: calculates CO2 compensation point.      
      
!------- Methods -------------------
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430

       double precision templeaf,gascon,
     &      gammas,gammas25,ha_gammas
     
!------- inputs -------------------- 
! templeaf:   leaf temperature [K]
! gascon:     universal gas constant [JK-1 mol-1].
! gammas25: the compensation point at 25oC [Pa]
! ha_gammas: parameter in the temperature response function, typically around
!       37.83 kJmol-1
!------- outputs -------------------
! gammas:    CO2 compensation point (Pa)

      gammas=gammas25*dexp((ha_gammas/gascon)*
     &    (1000.0d0/298.15d0-1000.0d0/templeaf))
      return
      end
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           
      subroutine resp_mitocho(templeaf,rdlight25,
     &      ha_darkresp,gascon,rdlight)
      implicit none
      
! purpose: calculates mitochondrial respiration rate and its 
!          derivative wrt. leaf temperature.      
      
!------- Methods -------------------
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430

       double precision templeaf,rdlight25,gascon,
     &                 rdlight,ha_darkresp,term
     
!------- inputs -------------------- 
! templeaf:   leaf temperature [K]
! rdlight25:  Mitochondrial respiration rate in the light at 25oC
!             of leaf temperature [umol m-2 s-1]
! gascon:     universal gas constant [JK-1 mol-1].
! ha_darkresp: Mitochondrial respiation parameter deltaha = 46.39d0 unit in
!              Jmol-1/1000

!------- outputs -------------------
! rdlight:    Mitochondrial respiration rate in the light [umol m-2 s-1]
      if(ha_darkresp.lt.0.0d0)then
!The calling program uses the Q10 form. ha_darkresp is really -Q10
        term=-ha_darkresp
        rdlight=rdlight25*
     &  (term**((templeaf-298.15d0)/10.0d0))
      else
        rdlight=rdlight25*dexp((ha_darkresp/gascon)*
     &    (1000.0d0/298.15d0-1000.0d0/templeaf))
      endif
      return
      end
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      subroutine MichaelisCO2(templeaf,fkc25,ha_kc,gascon,
     &fkc)
      implicit none
      
! purpose: calculates Michaelis constant for CO2 at a given leaf temperature.      
      
!------- Methods -------------------
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430

      double precision templeaf,fkc25,ha_kc,
     &   gascon,fkc
     
!------- inputs -------------------- 
! templeaf:   leaf temperature [K]
! gascon:     universal gas constant [JK-1 mol-1].
! fkc25:      the Michaelis constant for CO2 at 25oC
! ha_kc:      a coefficient (~79.43 KJmol-1)
!------- outputs -------------------
! fkc:      Michaelis constant for CO2 (Pa)      
      
      fkc=fkc25*dexp((ha_kc/gascon)*
     &    (1000.0d0/298.15d0-1000.0d0/templeaf))      
      return
      end 
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      subroutine MichaelisO2(templeaf,fko25,ha_ko,gascon,
     &fko)
      implicit none
      
! purpose: calculates Michaelis constant for O2 at a given leaf temperature.      
      
!------- Methods -------------------
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430

      double precision templeaf,fko25,ha_ko,
     &   gascon,fko
     
!------- inputs -------------------- 
! templeaf:   leaf temperature [K]
! gascon:     universal gas constant [JK-1 mol-1].
! fko25:      the Michaelis constant for O2 at 25oC
! ha_ko:      a coefficient in the temp function (~36.38 KJmol-1)
!------- outputs -------------------
! fko:      Michaelis constant for CO2 (Pa)      
      
      fko=fko25*dexp((ha_ko/gascon)*
     &    (1000.0d0/298.15d0-1000.0d0/templeaf))      
      return
      end
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!    
      subroutine tpuontemp(templeaf,gascon,tpu25,
     &                     dha,dhd,ds,tpu)
      implicit none
!
!(in) templeaf: leaf temperature [K]
!(in) gascon:   universal gas constant [JK-1 mol-1].
!(in) tpu25:    tpu at 25 oC [umol co2 m-2 s-1]
!(in) dha ~   53.1d0  (53100/1000)
!(in) dhd ~   201.8d0 (201800/1000)
!(in) ds  ~   0.65d0  (650/1000)
!(out) tpu:      tpu at the given leaf temperature     
     
      double precision templeaf,gascon,tpu25,
     &                 dha,dhd,ds,tpu
!
!-------------- Locals --------------------
      double precision t25k,univR,fact1,fact2,fact3,threeexp
      
      t25k=273.15d0+25.0d0
      univR=gascon/1000.0d0      
      fact1=(ds*templeaf-dhd)/(univR*templeaf)
      fact2=(ds*t25k-dhd)/(univR*t25k)
      fact3=dha/(univR*t25k)-dha/(univR*templeaf)
         
      tpu=tpu25*threeexp(fact1,fact2,fact3)    
      return
      end

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           
      subroutine vcmaxontemp(templeaf,vcmax25,gascon,
     &  ha_vcmax,hd_vcmax,sv_vcmax,vcmax)
      implicit none
      
! purpose: calculates maximum RuBP saturated rate of carboxylation.      
      
!------- Methods -------------------
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430

       double precision templeaf,vcmax25,gascon,ha_vcmax,hd_vcmax,
     &                  sv_vcmax,vcmax
!------- inputs -------------------- 
! templeaf:   leaf temperature [K]
! vcmax25:    maximum RuBP saturated rate of carboxylation at 25oC
!             of leaf temperature [umol m-2 s-1]
! gascon:     universal gas constant [JK-1 mol-1].
! ha_vcmax:   
! hd_vcmax:   
! sv_vcmax:

!------- outputs -------------------
! vcmax:      maximum RuBP saturated rate of carboxylation [umol m-2 s-1]    

!------- Local variables -----------      
      double precision threeexp
      double precision T0,univR,fact1,fact2,fact3
      
      univR=gascon/1000.0d0

      T0=25.0d0+273.15d0
      fact1=(sv_vcmax*templeaf-hd_vcmax)/(univR*templeaf)
      fact2=(sv_vcmax*T0-hd_vcmax)/(univR*T0)
      fact3=ha_vcmax*(1.0d0-T0/templeaf)/(univR*T0)

!     threeexp=(dexp(fact3)+dexp(fact2+fact3))/(1.0d0+dexp(fact1))

      vcmax=vcmax25*threeexp(fact1,fact2,fact3)
      
      return
      end
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           
      subroutine gmesoontemp(templeaf,gmeso25,gascon,
     &  ha_gmeso,hd_gmeso,sv_gmeso,gmeso)
      implicit none
      
! purpose: calculates mesophyll CO2 conductance at a given temperature     
      
!------- Methods -------------------
! Based on functions in Bernacchi, et al. (2002),
! Plant Physiology, 130, 1992-1998, or Scafaro et al. (2011), PCE 34: 1999-2008

       double precision templeaf,gmeso25,gascon,ha_gmeso,hd_gmeso,
     &                  sv_gmeso,gmeso
!------- inputs -------------------- 
! templeaf:   leaf temperature [K]
! ha_gmeso:  The activation energy (Ha) in gmeso temperature response function [kJmol-1]
! hd_gmeso:  The deactivation energy (Hd) in gmeso temperature response function [kJmol-1]
! sv_gmeso:  The entropy term (Sv) in gmeso temperature response function [kJmol-1K-1]
! gascon: The universal gas constant (JK-1mol-1)

!------- outputs -------------------
! gmeso: Mesophyll CO2 transfer conductance [umol m-2 s-1 Pa-1]

!------- Local variables -----------      
      double precision T0,univR,fact1,fact2,
     &   fact3,threeexp
      
      univR=gascon/1000.0d0
      T0=25.0d0+273.15d0
!
      if(dabs(sv_gmeso+9999.0d0).lt.1.0d-5.or.
     &dabs(hd_gmeso+9999.0d0).lt.1.0d-5)then
        gmeso=gmeso25*dexp(ha_gmeso*(1.0d0-T0/templeaf)/(univR*T0))
        return
      endif
      fact1=(sv_gmeso*templeaf-hd_gmeso)/(univR*templeaf)
      fact2=(sv_gmeso*T0-hd_gmeso)/(univR*T0)
      fact3=ha_gmeso*(1.0d0-T0/templeaf)/(univR*T0)
!     threeexp=(dexp(fact3)+dexp(fact2+fact3))/(1.0d0+dexp(fact1))
      gmeso=gmeso25*threeexp(fact1,fact2,fact3)
      return
      end
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           
      subroutine jmaxontemp(templeaf,fjmax25,gascon,
     &   ha_jmax,hd_jmax,sv_jmax,fjmax)
      implicit none
      
! purpose: calculates maximum rate of electron transport
!          and its derivative wrt. leaf temperature.      
      
!------- Methods -------------------
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430

! Sept. 16, 2005. Bernacchi et al. functions are replaced by
! June et al. function, June T., J. R. Evans, G. D. Farquhar (2004)
! A simple new equation for the reversible temperature dependence of
! photosynthetic electron transport: a study on soybean leaf,
! Functional Plant Biology 31, 275-283.


       double precision templeaf,fjmax25,gascon,
     &         ha_jmax,hd_jmax,sv_jmax,fjmax
     
!------- inputs -------------------- 
! templeaf:   leaf temperature [K]
! fjmax25:    maximum rate of electron transport at 25oC
!             of leaf temperature [umol m-2 s-1]
! gascon:     universal gas constant [JK-1 mol-1].
! ha_jmax ~ temperature response function parameter
! hd_jmax ~ temperature response function parameter
! sv_jmax ~ temperature response function parameter

!------- outputs -------------------
! fjmax:      maximum rate of electron transport [umol m-2 s-1]    
            
!------- Local variables -----------
      double precision T0,omega,threeexp
      double precision univR,fact1,fact2,fact3
 
      univR=gascon/1000.0d0
      
      T0=25.0d0+273.15d0
      fact1=(sv_jmax*templeaf-hd_jmax)/(univR*templeaf)
      fact2=(sv_jmax*T0-hd_jmax)/(univR*T0)
      fact3=ha_jmax*(1.0d0-T0/templeaf)/(univR*T0)

!     threeexp=(dexp(fact3)+dexp(fact2+fact3))/(1.0d0+dexp(fact1))

      fjmax=fjmax25*threeexp(fact1,fact2,fact3)

! June et al. function
!      ha_jmax = 30.0d0+273.15d0
!      hd_jmax = 11.6d0
!      sv_jmax = 0.18d0
!      T0=ha_jmax-273.15d0
!      omega=hd_jmax+sv_jmax*T0
!      fjmax=fjmax25*dexp(
!     &  (templeaf-273.15d0+25.0d0-2.0d0*T0)*
!     &  (25.0d0-templeaf+273.15d0)/(omega*omega))

      return
      end
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           
      subroutine jontemp(aPPFDlf,templeaf,fjelect,fjmax25,ha_jmax,
     &hd_jmax,sv_jmax,gascon,phifactor,thetafactor,betaPSII)
      implicit none
      
! purpose: calculates rate of electron transport
!          and its derivative wrt. leaf temperature.      
      
!------- Methods -------------------
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430
!3/29/2014. Bernacchi method cannot fit A/PAR curves well. phifactor and thetafactor are
!used to modify his method when A/PAR curves are present.
       double precision aPPFDlf,templeaf,fjelect,fjmax25,gascon,
     &ha_jmax,hd_jmax,sv_jmax,phifactor,thetafactor,betaPSII
!------- inputs -------------------- 
! aPPFDlf:       absorbed photosynthetic photon flux density by leaf (umol m-2 s-1)
! templeaf:      leaf temperature [K]
! fjmax25:       maximum rate of electron transport at 25oC
!                  of leaf temperature [umol m-2 s-1]
! gascon:        universal gas constant [JK-1 mol-1].
! ha_jmax ~ 50300.0d0/1000.0d0
! hd_jmax ~ 152044.0d0/1000.0d0
! sv_jmax ~ 495.0d0/1000.0d0
! phifactor:     modifier for phiPSIImax
! thetafactor:   modifer for thetaPSII
! betaPSII:      fraction of absorbed light that reaches photosystem II (~0.5)
!------- outputs -------------------
! fjelect:      rate of electron transport [umol m-2 s-1]    

!------- Local variables -----------    
      double precision Q2,phiPSIImax,
     &thetaPSII,zeroT,b2_4ac,b2_4ac_sqrt,fjmax
      parameter (zeroT=273.15d0)
!
! beta:        fraction of absorbed qunata reaching pSII [--]
! Q2:          maximum incident quanta that can be utilized in
!                   electron transport (umol m-2 s-1)
! phiPSIImax:  Maximum dark-adapted quantum yield of PSII [--]
! thetaPSII:   convexity term for electron transport rates [--]
! zeroT:       zero temperature [K]
! b2-4ac:      temporay variable
! b2_4ac_sqrt:      sqrt of b2_4ac
! fjmax:       maximum rate of electron transport, computed in the
!                subroutine of jmaxontime [umol m-2 s-1]    
      call thetaphipsii(templeaf,phiPSIImax,thetaPSII)
      phiPSIImax=phiPSIImax*phifactor
      thetaPSII=thetaPSII*thetafactor
      Q2=aPPFDlf*phiPSIImax*betaPSII
      call jmaxontemp(templeaf,fjmax25,gascon,
     &  ha_jmax,hd_jmax,sv_jmax,fjmax)
      b2_4ac=(Q2+fjmax)*(Q2+fjmax)-4.0d0*thetaPSII*Q2*fjmax
      if(b2_4ac.lt.0.0d0)then
         b2_4ac=0.0d0
      endif
      b2_4ac_sqrt=dsqrt(b2_4ac)
      fjelect=(Q2+fjmax-b2_4ac_sqrt)/(2.0d0*thetaPSII)
      return
      end
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      subroutine jmaxfromj(aPPFDlf,templeaf,fjelect,fjmax,
     &phifactor,thetafactor,betaPSII)
      implicit none      
! purpose: calculates jmax from the rate of electron transport,
! leaf temperature and absorbed PAR
      
!------- Methods -------------------
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430
!3/29/2014. Bernacchi method cannot fit A/PAR curves well. phifactor and thetafactor are
!used to modify his method when A/PAR curves are present.
!
       double precision aPPFDlf,templeaf,fjelect,fjmax,
     &phifactor,thetafactor,betaPSII
!------- inputs -------------------- 
! aPPFDlf:       absorbed photosynthetic photon flux density by leaf (umol m-2 s-1)
! templeaf:      leaf temperature [K]
! fjelect:       rate of electron transport [umol m-2 s-1]    
! phifactor:     modifier for phiPSIImax
! thetafactor:   modifer for thetaPSII
! betaPSII:      the fracction of absorbed light that reaches photosystem II (~0.5)
!------- outputs -------------------
! fjmax:       maximum rate of electron transport
!                  of leaf temperature [umol m-2 s-1]

!------- Local variables -----------    
      double precision Q2,phiPSIImax,thetaPSII,zeroT
      parameter (zeroT=273.15d0)
!
! beta:        fraction of absorbed qunata reaching pSII [--]
! Q2:          maximum incident quanta that can be utilized in
!                   electron transport (umol m-2 s-1)
! phiPSIImax:  Maximum dark-adapted quantum yield of PSII [--]
! thetaPSII:   convexity term for electron transport rates [--]
! zeroT:       zero temperature [K]
!
      call thetaphipsii(templeaf,phiPSIImax,thetaPSII) 
      phiPSIImax=phiPSIImax*phifactor
      thetaPSII=thetaPSII*thetafactor
      Q2=aPPFDlf*phiPSIImax*betaPSII
      fjmax=(thetaPSII*fjelect-Q2)*fjelect/(fjelect-Q2)
      return
      end

      subroutine thetaphipsii(templeaf0,phiPSIImax,thetaPSII)
      implicit none
! Based on functions in Bernacchi, Pimentel and Long (2003),
! Plant, Cell and Environment, 26, 1419-1430

! templeaf0:      leaf temperature [K] 
!------- outputs -------------------
! phiPSIImax:  Maximum dark-adapted quantum yield of PSII [--]
! thetaPSII:   convexity term for electron transport rates [--]
      double precision phiPSIImax,thetaPSII,templeaf0   
      double precision zeroT,templeaf
      parameter (zeroT=273.15d0)
! zeroT:       zero temperature [K]

      templeaf=dmax1(zeroT+5.0d0,templeaf0)
      templeaf=dmin1(zeroT+45.0d0,templeaf)
      phiPSIImax=0.352d0+0.022d0*(templeaf-zeroT)-3.4d-4*
     &  (templeaf-zeroT)*(templeaf-zeroT)
      thetaPSII=0.76d0+0.018d0*(templeaf-zeroT)-3.7d-4*
     &  (templeaf-zeroT)*(templeaf-zeroT)
      return
      end
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
      double precision function atpcon(wc,wj,vcmax)
      implicit none
!
!  Calculate the concentration of ATP provided by photophosphorylation
!  based on Buckley et al. (2003), Plant, Cell and Environment, 26, 1767-1785.
!
! ------ Inputs ---------------------------------------------------
! wc:           Rubisco limited rate of carboxylation [umol m-2 s-1] 
! wj:           RuBP limited rate of carboxylation [umol m-2 s-1]
! vcmax:         maximum RuBP saturated rate of carboxylation [umol m-2 s-1]
      double precision wc,wj,vcmax
!
! ------ Outputs ---------------------------------------------------
! atpcon: the concentration of ATP (mmol ATP m-2)

! ------ Local variables --------------------------------------------
! at:   at in Buckley et al (2003), total concentration of adenylates 
!                  [mmol AxP m-2]
! p:    p in Buckley et al (2003), concentration of photophosphorylation 
!                  sites [mmol sites m-2]
!          Communication with Buckley reveals error in Table 1 in Buckley paper
!          at and p should be divided by 100.
! vr:   Vr in Buckley et al (2003), carboxylation rate limited by potential
!                  RuBP pool size only [umol m-2 s-1]
! tau0: basal ATP level provided by other processes [mmol ATP m-2]
! tauc: atpcon when wc limiting mmol m-2
! tauj: atpcon when wj limiting mmol m-2

      double precision at,p,vr,tau0,tauc,tauj
      parameter(tau0=1.6d0)
      
      if(wj.le.0.0d0)then
        atpcon=tau0
      else
        at=12.6d0*vcmax/100.0d0
        p=2.5d0*vcmax/100.0d0
        vr=2.27d0*vcmax
        tauc=at-p*wc/wj
        tauj=(at-p)*(vr/vcmax-1.0d0)/(wc*vr/(wj*vcmax)-1.0d0)
        if(wc.lt.wj)then
          atpcon=tau0+tauc
        else
          atpcon=tau0+tauj
        endif
      endif      
      return
      end
!
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      subroutine QSat (T, p, es, esdT, qs,
     &                 qsdT)
      implicit none
!-----------------------------------------------------------------------
!
!  CLMCLMCLMCLMCLMCLMCLMCLMCLMCL  A community developed and sponsored, freely
!  L                           M  available land surface process model.
!  M --COMMUNITY LAND MODEL--  C
!  C                           L
!  LMCLMCLMCLMCLMCLMCLMCLMCLMCLM
!
!-----------------------------------------------------------------------
! Purpose:
! Computes saturation mixing ratio and the change in saturation
! mixing ratio with respect to temperature.

! Method:
! Reference:  Polynomial approximations from:
!             Piotr J. Flatau, et al.,1992:  Polynomial fits to saturation 
!             vapor pressure.  Journal of Applied Meteorology, 31, 1507-1513.
!
! Author:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision 
! April 2002: Vertenstein/Oleson/Levis; Final form
!
!-----------------------------------------------------------------------
! $Id: QSat.F90,v 1.1.8.4 2002/04/10 19:25:07 mvertens Exp $
!-----------------------------------------------------------------------

!  use precision
!  use shr_const_mod, ONLY: SHR_CONST_TKFRZ
!  implicit none

!----Arguments----------------------------------------------------------

      DOUBLE PRECISION T, p
      DOUBLE PRECISION es, esdT, qs, qsdT
      DOUBLE PRECISION SHR_CONST_TKFRZ
      PARAMETER (SHR_CONST_TKFRZ=273.16d0)
      
!  real(r8), intent(in)  :: T        ! temperature (K)
!  real(r8), intent(in)  :: p        ! surface atmospheric pressure (pa)

!  real(r8), intent(out) :: es       ! vapor pressure (pa)
!  real(r8), intent(out) :: esdT     ! d(es)/d(T)
!  real(r8), intent(out) :: qs       ! humidity (kg vapor/kg moist air)
!  real(r8), intent(out) :: qsdT     ! d(qs)/d(T)

!----Local Variables----------------------------------------------------

      DOUBLE PRECISION T_limit
      DOUBLE PRECISION td,vp,vp1,vp2
      DOUBLE PRECISION a0,a1,a2,a3,a4,a5,a6,a7,a8
      DOUBLE PRECISION b0,b1,b2,b3,b4,b5,b6,b7,b8
      DOUBLE PRECISION c0,c1,c2,c3,c4,c5,c6,c7,c8
      DOUBLE PRECISION d0,d1,d2,d3,d4,d5,d6,d7,d8

!  real(r8) T_limit                  ! limitation on valid temperatures [C]

!  real(r8) td,vp,vp1,vp2
!  real(r8) a0,a1,a2,a3,a4,a5,a6,a7,a8
!  real(r8) b0,b1,b2,b3,b4,b5,b6,b7,b8

!  real(r8) c0,c1,c2,c3,c4,c5,c6,c7,c8
!  real(r8) d0,d1,d2,d3,d4,d5,d6,d7,d8

!
! For water vapor (temperature range 0C-100C)
!
      data a0/6.11213476d0/,a1/0.444007856d0/,a2/0.143064234d-01/,
     & a3/0.264461437d-03/,a4/0.305903558d-05/,a5/0.196237241d-07/,
     & a6/0.892344772d-10/,a7/-0.373208410d-12/,a8/0.209339997d-15/

!
! For derivative:water vapor
!
      data b0/0.444017302d0/,b1/0.286064092d-01/,b2/0.794683137d-03/,
     & b3/0.121211669d-04/,b4/0.103354611d-06/,b5/0.404125005d-09/,
     & b6/-0.788037859d-12/,b7/-0.114596802d-13/,b8/0.381294516d-16/

!
! For ice (temperature range -75C-0C) 
!
      data c0/6.11123516d0/,c1/0.503109514d0/,c2/0.188369801d-01/,
     & c3/0.420547422d-03/,c4/0.614396778d-05/,c5/0.602780717d-07/,
     & c6/0.387940929d-09/,c7/0.149436277d-11/,c8/0.262655803d-14/

!
! For derivative:ice  
!
      data d0/0.503277922d0/,d1/0.377289173d-01/,d2/0.126801703d-02/,
     & d3/0.249468427d-04/,d4/0.313703411d-06/,d5/0.257180651d-08/,
     & d6/0.133268878d-10/,d7/0.394116744d-13/,d8/0.498070196d-16/

!----End Variable List--------------------------------------------------

      T_limit = T - SHR_CONST_TKFRZ
      if (T_limit .GT. 100.0d0) T_limit=100.0d0
      if (T_limit .LT. -75.0d0) T_limit=-75.0d0

      td   = T_limit
      if(td .GE. 0.0d0) then
         es=a0+td*(a1+td*(a2+td*(a3+td*(a4+
     &      td*(a5 + td*(a6 + td*(a7 + td*a8)))))))
         esdT=b0+td*(b1+td*(b2+td*(b3+td*(b4+
     &        td*(b5+td*(b6+td*(b7+td*b8)))))))
      else
         es=c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4
     &     + td*(c5 + td*(c6 + td*(c7 + td*c8)))))))
         esdT=d0+td*(d1+td*(d2 + td*(d3 + td*(d4
     &     + td*(d5 + td*(d6 + td*(d7 + td*d8)))))))
      endif

      es    = es    * 100.0d0
!                            pa
      esdT  = esdT  * 100.0d0
!                            pa/K
      vp    = 1.00d0   / (p - 0.378d0*es)
      vp1   = 0.622d0 * vp
      vp2   = vp1   * vp

      qs    = es    * vp1
!                        kg/kg
      qsdT  = esdT  * vp2 * p
!                            1 / K
      end

      double precision function esat(T, p)
      implicit none
!-----------------------------------------------------------------------
!
!  CLMCLMCLMCLMCLMCLMCLMCLMCLMCL  A community developed and sponsored, freely
!  L                           M  available land surface process model.
!  M --COMMUNITY LAND MODEL--  C
!  C                           L
!  LMCLMCLMCLMCLMCLMCLMCLMCLMCLM
!
!-----------------------------------------------------------------------
! Purpose:
! Computes saturation mixing ratio and the change in saturation
! mixing ratio with respect to temperature.

! Method:
! Reference:  Polynomial approximations from:
!             Piotr J. Flatau, et al.,1992:  Polynomial fits to saturation 
!             vapor pressure.  Journal of Applied Meteorology, 31, 1507-1513.
!
! Author:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision 
! April 2002: Vertenstein/Oleson/Levis; Final form
!
!-----------------------------------------------------------------------
! $Id: QSat.F90,v 1.1.8.4 2002/04/10 19:25:07 mvertens Exp $
!-----------------------------------------------------------------------

!  use precision
!  use shr_const_mod, ONLY: SHR_CONST_TKFRZ
!  implicit none

!----Arguments----------------------------------------------------------

      DOUBLE PRECISION T, p
      DOUBLE PRECISION es
      DOUBLE PRECISION SHR_CONST_TKFRZ
      PARAMETER (SHR_CONST_TKFRZ=273.15d0)
      
!  real(r8), intent(in)  :: T        ! temperature (K)
!  real(r8), intent(in)  :: p        ! surface atmospheric pressure (pa)

!  real(r8), intent(out) :: es       ! vapor pressure (pa)

!----Local Variables----------------------------------------------------

      DOUBLE PRECISION T_limit
      DOUBLE PRECISION td
      DOUBLE PRECISION a0,a1,a2,a3,a4,a5,a6,a7,a8
      DOUBLE PRECISION b0,b1,b2,b3,b4,b5,b6,b7,b8
      DOUBLE PRECISION c0,c1,c2,c3,c4,c5,c6,c7,c8


!  real(r8) T_limit                  ! limitation on valid temperatures [C]

!  real(r8) td
!  real(r8) a0,a1,a2,a3,a4,a5,a6,a7,a8

!  real(r8) c0,c1,c2,c3,c4,c5,c6,c7,c8

! For water vapor (temperature range 0C-100C)
!
      data a0/6.11213476d0/,a1/0.444007856d0/,a2/0.143064234d-01/,
     & a3/0.264461437d-03/,a4/0.305903558d-05/,a5/0.196237241d-07/,
     & a6/0.892344772d-10/,a7/-0.373208410d-12/,a8/0.209339997d-15/

!
! For ice (temperature range -75C-0C) 
!
      data c0/6.11123516d0/,c1/0.503109514d0/,c2/0.188369801d-01/,
     & c3/0.420547422d-03/,c4/0.614396778d-05/,c5/0.602780717d-07/,
     & c6/0.387940929d-09/,c7/0.149436277d-11/,c8/0.262655803d-14/


!----End Variable List--------------------------------------------------

      T_limit = T - SHR_CONST_TKFRZ
      if (T_limit .GT. 100.0d0) T_limit=100.0d0
      if (T_limit .LT. -75.0d0) T_limit=-75.0d0

      td   = T_limit
      if(td .GE. 0.0d0) then
         es=a0+td*(a1+td*(a2+td*(a3+td*(a4+
     &      td*(a5 + td*(a6 + td*(a7 + td*a8)))))))
      else
         es=c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4
     &     + td*(c5 + td*(c6 + td*(c7 + td*c8)))))))
      endif

      es    = es    * 100.0d0
!                            pa
      esat=es

!The following is what used in Li-Cor 6400 so we use it for A/Ci analysis
      esat=613.65d0*dexp(17.502d0*(T-273.15d0)/
     &  (240.97d0+T-273.15d0))
      return
      end
!
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      subroutine metvarconverter(airtemp,pres_moistair,vpd,
     &    pres_h2ovapor,rh,gascon,vapormixingratio,spechumid,
     &    cp_moistair,dens_moistair,dewpoint,wmoistair)
      implicit none     
!     
! =============Inputs===================================
! airtemp:       air temperature [K]
! pres_moistair: air pressure [Pa]
! rh or vpd or pres_h2ovapor: one of the three must be inputs
!    rh:            relative humidity [%]
!    vpd:           vapor pressure deficit [Pa]
!    pres_h2ovapor: water vapor partial pressure [Pa]
! gascon: universal gas constant [J K-1 mol-1]
!
      double precision airtemp,pres_moistair,vpd,
     &    pres_h2ovapor,rh,gascon
     
!==============Outputs===================================
! rh, vpd, and pres_h2ovapor: if not provided as inputs
! vapormixingratio: water vapor mixing ratio [g vapor / g dry air]
! spechumid:        specific humidity [g vapor / g moist air]
! cp_moistair:      specific heat of the moist air [J g-1 K-1]
! dens_moistair:    density of the moist air [gm-3] 
! dewpoint:         dewpoint temperature [K]
! wmoistair:        molecular weight of moist air [g mol-1]
!        
      double precision vapormixingratio,spechumid,
     &    cp_moistair,dens_moistair,dewpoint,wmoistair
     
!==============Locals===================================
! dryairmw:      molecular weight of the dry air [g mol-1]
! h2omw:         molecular weight of water [g mol-1]
! cp_dryair:     specific heat of dry air [J g-1 K-1]
! cp_watervapor: specific heat of water vapor [J g-1 K-1]
! es:            saturation vapor pressure [Pa]
! esdT:          not used
! qs:            specific humidity at saturation [g vapor / g moist air]
! qsdT:          not used
! pres_dryair:   partial pressure of the dry air [Pa]
! dens_dryair:   density of the dry air [gm-3]
! dens_vapor:    density of water vapor [gm-3]
! ws:            mixing ratio at saturation [ density vapor / density dry air]
! abszero:       absolute zero temperature [K]
! latent:        latent heat of vaporization [J g-1]

      double precision dryairmw,h2omw,cp_dryair,cp_watervapor,
     & es,esdT,qs,qsdT,pres_dryair,dens_dryair,dens_vapor,esat,
     & ws,abszero,latent
      parameter(dryairmw=28.9644d0,h2omw=18.016d0)
      parameter(cp_dryair=1.00467d0,cp_watervapor=1.84d0)
      parameter(abszero=273.15d0)
      
!           
      if((rh.le.0.0d0.or.rh.gt.100.02d0).and.(vpd.lt.0.0d0.or.
     &  vpd.gt.9000.0d0).and.(pres_h2ovapor.le.0.0d0.or.
     &  pres_h2ovapor.gt.9000.0d0))then
        write(*,*)'At least of the one variables VPD, RH, and
     & water vapor partial pressure must be valid.
     & check your data. Program stops'
        stop
      endif
      call QSat(airtemp,pres_moistair,es,esdT,qs,qsdT)
      ws=qs/(1.0d0-qs)
      
      if((rh.le.0.0d0.or.rh.gt.100.0d0).or.(vpd.le.0.0d0.or.
     &  vpd.gt.9000.0d0).or.(pres_h2ovapor.le.0.0d0.or.
     &  pres_h2ovapor.gt.9000.0d0))then
        if(pres_h2ovapor.gt.0.0d0.and.pres_h2ovapor.lt.9000.0d0)then
          pres_dryair=pres_moistair-pres_h2ovapor
          dens_dryair=pres_dryair*dryairmw/(gascon*airtemp)
          dens_vapor=pres_h2ovapor*h2omw/(gascon*airtemp)
          vapormixingratio=dens_vapor/dens_dryair
          vpd=es-pres_h2ovapor
          vpd=dmin1(es,vpd)
          vpd=dmax1(0.0d0,vpd)
          rh=vapormixingratio*100.0d0/ws
          rh=dmin1(100.0d0,rh)
          rh=dmax1(0.0d0,rh)
        else
          if(rh.gt.0.0d0.and.rh.le.100.01d0)then
            vapormixingratio=ws*dmin1(rh,100.0d0)/100.0d0
            pres_h2ovapor=vapormixingratio*pres_moistair/(
     &           vapormixingratio+h2omw/dryairmw)
            vpd=es-pres_h2ovapor
          else
            if(vpd.ge.0.0d0.and.vpd.lt.9000.0d0)then
              pres_h2ovapor=es-vpd
              pres_h2ovapor=dmin1(es,pres_h2ovapor)
              pres_h2ovapor=dmax1(0.0d0,pres_h2ovapor)
              pres_dryair=pres_moistair-pres_h2ovapor
              dens_dryair=pres_dryair*dryairmw/(gascon*airtemp)
              dens_vapor=pres_h2ovapor*h2omw/(gascon*airtemp)
              vapormixingratio=dens_vapor/dens_dryair        
              rh=vapormixingratio*100.0d0/ws  
            endif
          endif
        endif
      endif      
      pres_dryair=pres_moistair-pres_h2ovapor
      dens_dryair=pres_dryair*dryairmw/(gascon*airtemp)
      dens_vapor=pres_h2ovapor*h2omw/(gascon*airtemp)
      dens_moistair=dens_dryair+dens_vapor
      vapormixingratio=dens_vapor/dens_dryair
      spechumid=dens_vapor/dens_moistair
      cp_moistair=(cp_dryair*dens_dryair+
     &          cp_watervapor*dens_vapor)/dens_moistair
      
      latent=1.91846d+3*(airtemp/(airtemp-33.91d0))*
     &       (airtemp/(airtemp-33.91d0))
      call QSat(abszero,pres_moistair,es,esdT,qs,qsdT)
      dewpoint=1.0d0/abszero-(gascon/(h2omw*latent))*
     &         dlog(pres_h2ovapor/es)
      dewpoint=1.0d0/dewpoint
      wmoistair=dens_moistair*gascon*airtemp/pres_moistair      
      return
      end
!
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      double precision function latent(temp)
c
c This function calculates latent heat from temperature. The formulae
c is taken from Henderson-sellers, 1984, A new formulae for latent heat
c of vaporation of water as a function of temperature, Quart. J. R. Met.
c Soc. 1984, 110 pp. 1186-1190.
c
c temp:  temperature in K.
c latent: latent heat of water in J g-1.

      implicit double precision (a-h,l,o-z)
      parameter(abszero=273.15d0,fusionheat=333.7d0)
!
      latent=1.91846d+3*(temp/(temp-33.91d0))*
     &       (temp/(temp-33.91d0))
!
! assume freezing or sublimation occurs at 273.15d0 K
! latent heat of fusion is 333.7 J/g
!
      if(temp.le.abszero)then
        latent=latent+fusionheat
      endif
      return
      end
c   

!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
