!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&      
      subroutine cicaoptimization5(npoints,cicameas0,
     &pco2ambient0,beta,ndim,imodel0,bmin0,bmax0)
      implicit none
      
      include '../testarea/cica.h'
c
      integer npoints
      double precision pco2ambient0(npoints),cicameas0(npoints),
     &    acica,bcica,ccica,dcica,ecica
      
      integer i,ndim,imodel0
      double precision beta(ndim),fatbeta,ftol,bmin0(ndim),
     &    bmax0(ndim),f1dim_cica
      parameter(ftol=1.0d-7)
      external funkmin_cica,f1dim_cica
                 
      nobs = npoints
      imodel=imodel0
      do i=1,npoints
        pco2ambient(i)=pco2ambient0(i)
        cicameas(i)=cicameas0(i)
      enddo
      
      do i=1,ndim
        bmin(i)=bmin0(i)
        bmax(i)=bmax0(i)
      enddo
!
! Initialize the cost function evaluation counter in the subroutine funkmin.
! The counter counts and memorizes points where the cost function is evaluated.
      
      call funkmin_cica(ndim,beta,fatbeta)
      call nongradopt(ndim,funkmin_cica,f1dim_cica,beta,
     &      bmin,bmax,ftol,fatbeta)
      if(imodel.eq.3)then
        call RepeatCompassSearch(ndim,beta,fatbeta,bmin,
     &bmax,funkmin_cica,f1dim_cica,ftol)
      else
        call nongradopt(ndim,funkmin_cica,f1dim_cica,beta,
     &bmin,bmax,ftol,fatbeta)
      endif
      return
      END
c
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      subroutine cica_ca5(imodel,ndim,beta,ambco2_in_Pa,cica,der_cica,
     &der_beta)
      implicit none
      integer imodel,ndim
! calculate Ci/Ca ratio for a given ambient CO2 partial pressure in Pa.
!
!Ci/Ca=a*exp(-b*Ca)+c/(1+exp(-(Ca-d)/e))
!
      double precision a,b,c,d,e,f,ambco2_in_Pa,cica,der_cica,
     &term1,term2,term0,crit,grad(10),twoexpfunc,beta(ndim),
     &der_beta(ndim)
      parameter(crit=300.0d0)

      a=beta(1)
      b=beta(2)
      c=beta(3)
      d=beta(4)
      e=beta(5)
      f=beta(ndim)
      if(ndim.le.5)f=0.0d0
      if(imodel.eq.1)then
!        cica=twoexpfunc(y0,a1,b1,c1,x01,
!      &     a2,b2,c2,x02,x)
        cica=twoexpfunc(e,a,b,1.d0,0.0d0,
     &     c,d,1.0d0,0.0d0,ambco2_in_Pa)
!      call gradtwoexp(y0,a1,b1,c1,x01,
!     &     a2,b2,c2,x02,x,grad)
        call gradtwoexp(e,a,b,1.d0,0.0d0,
     &     c,d,1.0d0,0.0d0,ambco2_in_Pa,grad)
!        a1<->grad(1)<->a
!        b1<->grad(2)<->b
!        c1<->grad(3)
!        x01<->grad(4)
!        y0<->grad(5)<->e
!        x<->grad(6)      
!        a2<->grad(7)
!        b2<->grad(8)
!        c2<->grad(9)
!        x02<->grad(10)
        der_cica=grad(6)
        der_beta(5)=grad(5)
        der_beta(1)=grad(1)
        der_beta(2)=grad(2)
        der_beta(3)=grad(7)
        der_beta(4)=grad(8)
        return
      endif
      if(imodel.eq.2)then
        term0=-(ambco2_in_Pa-d)/e
        term1=dexp(-b*ambco2_in_Pa)
        if(term0.lt.crit)then
          term2=dexp(term0)
          cica=a*term1+c/(1.0d0+term2)
          der_cica=-a*b*term1+(c*term2)/(e*(1.0d0+term2)**2)
          der_beta(1)=term1
          der_beta(2)=-a*ambco2_in_Pa*term1
          der_beta(3)=1.0d0/(1.0d0+term2)
          der_beta(4)=-c*term2/(e*(1.0d0+term2)**2)
          der_beta(5)=c*term2*(-ambco2_in_Pa+d)/((e*(1.0d0+term2))**2)
        else
          term2=dexp(-term0)
          cica=a*term1+c*term2/(1.0d0+term2)
          der_cica=-a*b*term1+(c*term2)/(e*(1.0d0+term2)**2)
          der_beta(1)=term1
          der_beta(2)=-a*ambco2_in_Pa*term1
          der_beta(3)=term2/(1.0d0+term2)
          der_beta(4)=-c*term2/(e*(1.0d0+term2)**2)
          der_beta(5)=c*term2*(-ambco2_in_Pa+d)/((e*(1.0d0+term2))**2)
        endif
        if(ndim.eq.6)then
          cica=cica+beta(ndim)
          der_beta(ndim)=1.0d0
        endif
      endif
      if(imodel.eq.3)then
        term1=dexp(-b*ambco2_in_Pa)
        cica=a*term1+c+d*dlog(ambco2_in_Pa)+e*dlog(ambco2_in_Pa)*
     &dlog(ambco2_in_Pa)+f*dlog(ambco2_in_Pa)*dlog(ambco2_in_Pa)*
     &dlog(ambco2_in_Pa)
        der_cica=-a*b*term1+d/ambco2_in_Pa+2.0d0*e*dlog(ambco2_in_Pa)/
     &ambco2_in_Pa+3.0d0*f*dlog(ambco2_in_Pa)*dlog(ambco2_in_Pa)/
     &ambco2_in_Pa
        der_beta(1)=term1
        der_beta(2)=-a*ambco2_in_Pa*term1
        der_beta(3)=1.0d0
        der_beta(4)=dlog(ambco2_in_Pa)
        der_beta(5)=dlog(ambco2_in_Pa)*dlog(ambco2_in_Pa)
        if(ndim.eq.6)der_beta(ndim)=dlog(ambco2_in_Pa)*
     &dlog(ambco2_in_Pa)*dlog(ambco2_in_Pa)
      endif
      return
      end