!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 
      subroutine C4PhotoFit()
      implicit none
      include '../testarea/LeafGasParams.h'
      include '../testarea/LeafGasHybridFit.h'
      integer i,ndim,k,j,iderivative,iwrong,jnon,n,icompete,isame,i2,
     &isitnaninf,nave
      double precision beta(20),sumsquare0,beta0(20),sumsquarecp,
     &betacp(20),ftol,xtol,shortx(maxobs,10),shorty(maxobs,5),
     &ftol_relax,term1,term2,ran2,history(2000,25),discount,upper,lower,
     &f1dim_C4Fit,ff_pikaia
      parameter(ftol=1.0d-7,xtol=1.0d-7)
      external funkmin_C4Fit,f1dim_C4Fit,FCN_C4Fit,ff_pikaia
      ndim=3
      beta(1)=vcmax25_ori
      beta(2)=c4aparslope_ori
      beta(3)=c4kp25_ori
      betamin(1)=0.0d0
      betamax(1)=500.0d0
      betamin(2)=0.0d0
      betamax(2)=10.0d0
      betamin(3)=0.0d0
      betamax(3)=200000.0d0*betamax(1)
      if(idord.eq.1)then
        ndim=4
        beta(ndim)=rdlight25_ori
        betamin(ndim)=0.0d0
        betamax(ndim)=15.0d0
      endif
      isitbounded=1
      call funkmin_C4Fit(ndim,beta,sumsquare)
      do i=1,ndim
        beta0(i)=beta(i)
      enddo
      sumsquare0=sumsquare
      do i=1,ntotsamples
        responses(i,1)=anet_obs(i)
        forcings(i,1)=pco2i(i)
        forcings(i,2)=aPPFDlf(i)
        forcings(i,3)=templeaf(i)
        forcings(i,4)=pres_air(i)
        do j=1,4
          weitforcings(i,j)=1.0d0
        enddo
        weitresponses(i,1)=1.0d0
      enddo
      j=4
      i=1
      iderivative=0
      iwrong=0
      call odr_leastsquare(ndim,FCN_C4Fit,beta,ntotsamples,
     &forcings(1:ntotsamples,1:j),j,responses(1:ntotsamples,1:i),i,
     &weitforcings(1:ntotsamples,1:j),weitresponses(1:ntotsamples,1:i),
     &iderivative,shortx(1:ntotsamples,1:j),shorty(1:ntotsamples,1:i),
     &sumsquare,iwrong)
!after odr_leastsquare, forcing variables are destroyed. restore to the origninals
      do i=1,ntotsamples
        pco2i(i)=pco2i_ori(i)
        aPPFDlf(i)=aPPFDlf_ori(i)
        templeaf(i)=templeaf_ori(i)
        pres_air(i)=pres_air_ori(i)
      enddo
      call funkmin_C4Fit(ndim,beta,sumsquare)
      do i=1,10
        call nongradopt(ndim,funkmin_C4Fit,f1dim_C4Fit,beta,betamin,
     &betamax,ftol,sumsquare)
        call funkmin_C4Fit(ndim,beta,sumsquare)
        call RepeatCompassSearch(ndim,beta,sumsquare,
     &betamin,betamax,funkmin_C4Fit,f1dim_C4Fit,ftol)
        call funkmin_C4Fit(ndim,beta,sumsquare)
      enddo
      call ilimittypestats(ntotsamples,Postiphotolimit,
     &bestilimittype,bestnumrubis,bestnumrubp,bestnumtpu)
      return
      END subroutine C4PhotoFit