subroutine UnivPhotoFit()
      implicit none
      include '../testarea/LeafGasParams.h'
      include '../testarea/LeafGasHybridFit.h'
      integer i,j,k,Priornumrubis,Priornumrubp,Priornumtpu,
     &Currentnumrubis,Currentnumrubp,Currentnumtpu,Postnumrubis,
     &Postnumrubp,Postnumtpu,Postilimittype
      double precision term,term1,term2,term3,term4,term5,term6,
     &term7,term8,term9
!
      if(Prioriknowlimit.eq.1.or.Prioriknowlimit.eq.2)then
        ntotsamples=0
        do i=1,nFixedPoints
          ntotsamples=ntotsamples+1
          Prioriphotolimit(ntotsamples)=Fixediphotolimit(i)
        enddo
        do j=1,numACicurves
          do i=1,nACiPoints(j)
            ntotsamples=ntotsamples+1
            Prioriphotolimit(ntotsamples)=ACiiphotolimit(i,j)
          enddo
        enddo
        do j=1,numALightcurves
          do i=1,nALightPoints(j)
            ntotsamples=ntotsamples+1
            Prioriphotolimit(ntotsamples)=ALightiphotolimit(i,j)
          enddo
        enddo
        do i=1,nFreePoints
          ntotsamples=ntotsamples+1
          Prioriphotolimit(ntotsamples)=Freeiphotolimit(i)
        enddo
        call ilimittypestats(ntotsamples,Prioriphotolimit,
     &Priorilimittype,Priornumrubis,Priornumrubp,Priornumtpu)
        if(bestilimittype.gt.0.and.Priorilimittype.ne.bestilimittype)
     &return
!if bestilimittype is specified, we assume the fitting is constained to the limit type specified by bestilimittype
        if(Priorilimittype.ge.3)return
!5, 6, 7 are done with Prioriknowlimit=0. We don't consider cases with only rubp and tpu limitations but no rubisco limitations
        if(Priornumrubis.gt.0.and.Priornumrubis.lt.minimumrubis)return
        if(Priornumrubp.gt.0.and.Priornumrubp.lt.minimumfj)return
        if(Priornumtpu.gt.0.and.Priornumtpu.lt.minimumvt)return
        do i=1,ntotsamples
          Currentiphotolimit(i)=Prioriphotolimit(i)
        enddo
      endif
      Currentilimittype=Priorilimittype
      Currentiknowlimit=Prioriknowlimit
!
!-------------Test Area---------------------------------------------
!      Currentilimittype=1
!      Currentiknowlimit=1
!      do i=1,ntotsamples
!        if(i.le.6)then
!          Currentiphotolimit(i)=1
!        else
!          if(i.ge.16.and.i.le.26)then
!            Currentiphotolimit(i)=2
!          else
!            Currentiphotolimit(i)=3
!          endif
!        endif
!      enddo
!-------------------------------------------------------------------
      call DoUnivPhotoFit()
      if(Prioriknowlimit.ne.1.or.Priorilimittype.ge.5)goto 1000
!---------------------------------------------------------
!Enforce the admissibility rule.
!first get the post-fit limit type of each point. (pco2i,anet_obs) should be replaced
!by (pco2i_pred, anet_pred)
      do i=1,ntotsamples
        call leafunivphotosyn(Currentiknowlimit,Currentilimittype,
     &0,aPPFDlf(i),templeaf(i),pco2i_pred(i),po2i(i),chlflphips2(i),
     &term,weitresponses(i:i,1:1),weitresponses(i:i,1:1),
     &weitresponses(i:i,2:2),weitresponses(i:i,1:1),term1,term2,
     &Postiphotolimit(i),term3,term4,term5,term6,term7,term8,term9)
      enddo
      j=0
      do i=1,ntotsamples
        if(Postiphotolimit(i).ne.Currentiphotolimit(i))j=j+1
      enddo
!if j = 0, the fitting is admissible so go to the wrapup
      if(j.eq.0)goto 1000
      call ilimittypestats(ntotsamples,Postiphotolimit,
     &Postilimittype,Postnumrubis,Postnumrubp,Postnumtpu)
!if minimum number of points is not satisfied, go to penality fit directly.
      if(Postnumrubis.gt.0.and.Postnumrubis.lt.minimumrubis)goto 500
      if(Postnumrubp.gt.0.and.Postnumrubp.lt.minimumfj)goto 500
      if(Postnumtpu.gt.0.and.Postnumtpu.lt.minimumvt)goto 500
!check to see if the fit oscillates.
      Currentilimittype=Postilimittype
      do i=1,ntotsamples
        Currentiphotolimit(i)=Postiphotolimit(i)
      enddo
      call DoUnivPhotoFit()
      do i=1,ntotsamples
        call leafunivphotosyn(Currentiknowlimit,Currentilimittype,
     &0,aPPFDlf(i),templeaf(i),pco2i_pred(i),po2i(i),chlflphips2(i),
     &term,weitresponses(i:i,1:1),weitresponses(i:i,1:1),
     &weitresponses(i:i,2:2),weitresponses(i:i,1:1),term1,term2,
     &Postiphotolimit(i),term3,term4,term5,term6,term7,term8,term9)
      enddo
      j=0
      do i=1,ntotsamples
        if(Postiphotolimit(i).ne.Prioriphotolimit(i))j=j+1
      enddo
      if(j.eq.0)then
!Osicillation. Treat osicillating points as co-limited
        k=ntotsamples
        do i=1,ntotsamples
          if(Currentiphotolimit(i).ne.Prioriphotolimit(i))then
            k=k+1
            Currentiphotolimit(k)=Currentiphotolimit(i)
            aPPFDlf(k)=aPPFDlf(i)
            templeaf(k)=templeaf(i)
            po2i(k)=po2i(i)
            pco2i(k)=pco2i(i)
            anet_obs(k)=anet_obs(i)
            chlflphips2(k)=chlflphips2(i)
            Currentiphotolimit(i)=Prioriphotolimit(i)
          endif
        enddo
        call ilimittypestats(k,Currentiphotolimit,
     &Currentilimittype,Currentnumrubis,Currentnumrubp,Currentnumtpu)
        i=ntotsamples
        ntotsamples=k
        call DoUnivPhotoFit()
        sumsquare=sumsquare*dble(i)/dble(k)
        ntotsamples=i
        goto 1000
      else
!no osicillation
        Currentilimittype=Priorilimittype
        do i=1,ntotsamples
          Currentiphotolimit(i)=Prioriphotolimit(i)
        enddo
      endif
!-------------Penalty function fit-------------------------------------------
500   Currentiknowlimit=2
      call DoUnivPhotoFit()
!-------------Wrap up--------------------------------------------------------
1000  if(Prioriknowlimit.eq.0.and.Priorilimittype.le.4)then
        call ilimittypestats(ntotsamples,Postiphotolimit,
     &Postilimittype,Postnumrubis,Postnumrubp,Postnumtpu)
        if(Postnumrubis.gt.0.and.Postnumrubis.lt.minimumrubis-1)
     &sumsquare=1.0d+10
        if(Postnumrubp.gt.0.and.Postnumrubp.lt.minimumfj-1)
     &sumsquare=1.0d+10
        if(Postnumtpu.gt.0.and.Postnumtpu.lt.minimumvt-1)
     &sumsquare=1.0d+10
      endif
      term=1.0d0+(subbestsumsquare(Priorilimittype)-sumsquare)
      if(term.gt.1.0d0)then
        subbestsumsquare(Priorilimittype)=sumsquare
        do i=1,ntotunivparams
          subbestunivparams(i,Priorilimittype)=univparams(i)
        enddo
        do i=1,ntotsamples
          if(Prioriknowlimit.eq.0)then
            subbestiphotolimit(i,Priorilimittype)=Postiphotolimit(i)
          else
            subbestiphotolimit(i,Priorilimittype)=Prioriphotolimit(i)
          endif
        enddo
      endif
      return
      end subroutine UnivPhotoFit
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 
      subroutine DoUnivPhotoFit()
      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_UnivPhotoFit,ff_pikaia
      parameter(ftol=1.0d-7,xtol=1.0d-7)
      external funkmin_UnivPhotoFit,f1dim_UnivPhotoFit,
     &FCN_UnivPhotoFit,ff_pikaia
!find out which parameters to optimize
      call UnivParamsAlloc(1)
      ndim=0
      do i=1,ntotunivparams
        univparams(i)=subbestunivparams(i,Currentilimittype)
        if(ifixunivparams(i).eq.1)then
          ndim=ndim+1
          beta(ndim)=univparams(i)
          betamin(ndim)=univparamsmin(i)
          betamax(ndim)=univparamsmax(i)
        endif
      enddo
      isitbounded=1
      call funkmin_UnivPhotoFit(ndim,beta,sumsquare)
      if(subbestsumsquare(Currentilimittype).gt.1.0d+9)then
!global search
        do i=1,ndim
          beta0(i)=beta(i)
          history(1,i)=beta(i)
        enddo
        sumsquare0=sumsquare
        history(1,ndim+1)=sumsquare
!entrance counter
        history(1,ndim+2)=1.0d0
!failure counter
        history(1,ndim+3)=0.0d0
!Is it a competition among different initial guesses?
        icompete=0
!j the total number of calls to nongradopt; k is the number of returns to the current best and reset
!to zero if a better minumum is found; n is the number of scouting points over the landscape of the cost function.
!The first initial guess provided by the user is always part of the set of scouting points.the rest consist of outcomes
!from calls to nongradopt if they are significantly different from the current best.
        jnon=0
        isame=0
        n=1
        nave=n
        ftol_relax=ftol*1000.0d0
        discount=2.0d0
!relax the convergence criterion for scouting
30      do i=1,ndim
          betacp(i)=beta(i)
        enddo
        sumsquarecp=sumsquare
        iderivative=0
        if(ifitmode.lt.0)then
          iwrong=0
        else
          iwrong=1
        endif
        k=ifitmode
        ifitmode=-1
!ifitmode: =-2, ordinary fitting with pco2i calculated as a function of anet
!ifitmode: =-1, ordinary fitting with anet calculated as a function of pco2i
!ifitmode: =1, orthogonal fitting with anet calculated as a function of pco2i
!ifitmode: =2, orthogonal fitting with pco2i calculated as a function of anet
!for odr_leastsquare, only the predicted value of the response variable is needed, i.e., the cost function value is not needed.
!also, only anet as a function of CO2i is considered (not the other way around) because odr_leastsquare cannot handle the situation
!co2i as a function of anet for tpu limitation when alpha=0
        i=1
        if(ntotphips2.ge.1)i=2
        j=4
        if(Currentiknowlimit.eq.-1)then
!fluorescence only fit. chlflphisi2 becomes a forcing variable
          i=1
          j=5
        endif
        call odr_leastsquare(ndim,FCN_UnivPhotoFit,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)
        ifitmode=k
!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)
          po2i(i)=po2i_ori(i)
          chlflphips2(i)=chlflphips2_ori(i)
        enddo
        call funkmin_UnivPhotoFit(ndim,beta,sumsquare)
        if(isitnaninf(sumsquare).eq.1.or.sumsquare.gt.sumsquarecp)then
          do i=1,ndim
            beta(i)=betacp(i)
          enddo
          sumsquare=sumsquarecp
        else
          do i=1,ndim
            betacp(i)=beta(i)
          enddo
          sumsquarecp=sumsquare
        endif
        call nongradopt(ndim,funkmin_UnivPhotoFit,
     &f1dim_UnivPhotoFit,beta,betamin,betamax,ftol_relax,sumsquare)
        if(isitnaninf(sumsquare).eq.1.or.sumsquare.gt.sumsquarecp)then
          do i=1,ndim
            beta(i)=betacp(i)
          enddo
          sumsquare=sumsquarecp
        endif
        if(sumsquare.gt.1.0d0)then
          term1=sumsquare*ftol_relax
        else
          term1=ftol_relax*10.0d0
        endif
        if(sumsquare.gt.sumsquare0)then
!failure
          if((sumsquare-sumsquare0).gt.term1)then
            if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+1.5d0
!even though sumsquare is much worse than sumsquare0, it is an output of optimization after all so
!include it in the set if it has not already been included in the set.
            i=1
            i2=1
40          if(dabs(history(i2,i)-beta(i)).gt.ftol_relax)then
              if(dabs(history(i2,ndim+1)-sumsquare).lt.term1)then
                history(i2,ndim+3)=history(i2,ndim+3)+1.0d0
                goto 60
              endif
              if(i2.ge.n)goto 50
              i2=i2+1
              i=1
              goto 40
            else
              if(i.ge.ndim)goto 60 
              i=i+1
              goto 40
            endif
50          n=n+1
            do i=1,ndim
              history(n,i)=beta(i)
            enddo
            history(n,ndim+1)=sumsquare
            history(n,ndim+2)=0.0d0
            history(n,ndim+3)=0.0d0
!use average only when there is imporvement
            nave=n
          else
!the difference is minimal even though sumsquare is larger than sumsquare0.
!Increment the counter for arriving at the same minimum.
            if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+1.0d0
            isame=isame+1
          endif
60        do i=1,ndim
            beta(i)=beta0(i)
          enddo
          sumsquare=sumsquare0
        else
!success
          if((sumsquare0-sumsquare).lt.term1)then
!negligible improvement. Increment the counter for arriving at the same minimum.
!no increment for the set of central initial guesses
            if(icompete.eq.1)history(1,ndim+3)=history(1,ndim+3)+0.5d0
            isame=isame+1
            nave=n
          else
!reset the counter for arriving at a better minimum.
!Increment the set of central initial guesses
            if(dabs(discount-2.0d0).lt.ftol)then
              discount=dmax1(0.001d0,(sumsquare0-sumsquare)/1000.0d0)
            endif
            isame=0
            n=n+1
            do i=1,ndim
              history(n,i)=beta(i)
            enddo
            history(n,ndim+1)=sumsquare
            history(n,ndim+2)=0.0d0
            history(n,ndim+3)=0.0d0
          endif
          do i=1,ndim
            beta0(i)=beta(i)
          enddo
          sumsquare0=sumsquare
        endif
        jnon=jnon+1
        if(jnon.lt.200.and.isame.lt.3)then
!we first explore around the very first initial guess
          if(jnon.lt.10)then
            icompete=1
            term1=0.05d0+dmin1(history(1,ndim+3)*0.1d0,0.9d0)
            history(1,ndim+2)=history(1,ndim+2)+1.0d0
            do i=1,ndim
              lower=history(1,i)-term1*(history(1,i)-betamin(i))
              upper=history(1,i)+term1*(betamax(i)-history(1,i))
              beta(i)=lower+ran2()*(upper-lower)
            enddo
            goto 70
          endif
!try average if n is incremented
          if(n.gt.nave)then
            term1=1.0d0/(history(1,ndim+1)+1.0d-5)
            do i=2,n
              term1=term1+1.0d0/(history(i,ndim+1)+1.0d-5)
            enddo
            do i=1,ndim
              beta(i)=history(1,i)/(term1*(history(1,ndim+1)+1.0d-5))
              do icompete=2,n
                beta(i)=beta(i)+history(icompete,i)/
     &(term1*(history(icompete,ndim+1)+1.0d-5))
              enddo
            enddo
            nave=n
            icompete=0
            goto 70
          endif
!try different initial guesses
          if(ran2().gt.0.2d0)then
!guess around the best
            icompete=1
            term1=history(1,ndim+1)+
     &discount*history(1,ndim+2)*history(1,ndim+3)
            do i=2,n
              term2=history(i,ndim+1)+
     &discount*history(i,ndim+2)*history(i,ndim+3)
              if(term2.le.term1)then
                term1=term2
                do i2=1,ndim+3
                  history(n+1,i2)=history(i,i2)
                  history(i,i2)=history(1,i2)
                  history(1,i2)=history(n+1,i2)
                enddo
              endif
            enddo
            term1=0.05d0+dmin1(history(1,ndim+2)*history(1,ndim+3)*
     &0.015d0,0.9d0)
            history(1,ndim+2)=history(1,ndim+2)+1.0d0
            do i=1,ndim
              lower=history(1,i)-term1*(history(1,i)-betamin(i))
              upper=history(1,i)+term1*(betamax(i)-history(1,i))
              beta(i)=lower+ran2()*(upper-lower)
            enddo
          else
!completely random guess
            do i=1,ndim
              beta(i)=betamin(i)+ran2()*(betamax(i)-betamin(i))
            enddo
            icompete=0
          endif
70        call funkmin_UnivPhotoFit(ndim,beta,sumsquare)
          goto 30
        else
          if((ftol_relax-ftol).gt.ftol)then
            if(isame.le.1)then
              n=n+1
              do i=1,ndim+3
                history(n,i)=history(1,i)
              enddo
              do i=1,ndim
                history(1,i)=beta(i)
              enddo
              history(1,ndim+1)=sumsquare
              history(1,ndim+2)=0.0d0
              history(1,ndim+3)=0.0d0
              do i=1,n
                do icompete=1,ndim
                  betacp(icompete)=history(i,icompete)
                enddo
                sumsquarecp=history(i,ndim+1)
                call RepeatCompassSearch(ndim,betacp,sumsquarecp,
     &betamin,betamax,funkmin_UnivPhotoFit,f1dim_UnivPhotoFit,
     &ftol_relax)
                call funkmin_UnivPhotoFit(ndim,betacp,sumsquarecp)
                if(isitnaninf(sumsquarecp).eq.0.and.sumsquarecp.lt.
     &sumsquare)then
                  do icompete=1,ndim
                    beta(icompete)=betacp(icompete)
                  enddo
                  sumsquare=sumsquarecp
                endif
              enddo
              do i=1,ndim
                beta0(i)=beta(i)
              enddo
              sumsquare0=sumsquare
              jnon=0
              icompete=1
            else
              icompete=0
            endif
            ftol_relax=ftol
            goto 30
          endif
        endif

        goto 110

        do i=1,ndim
          betacp(i)=(beta(i)-betamin(i))/(betamax(i)-betamin(i))
        enddo
        sumsquarecp=sumsquare
        call pikaia(ff_pikaia,ndim,gacontrol,betacp,sumsquarecp,i)
        if(i.eq.0)then
          do i=1,ndim
            betacp(i)=betamin(i)+betacp(i)*(betamax(i)-betamin(i))
          enddo
          call funkmin_UnivPhotoFit(ndim,betacp,sumsquarecp)
          if(isitnaninf(sumsquarecp).eq.0.and.sumsquare.gt.sumsquarecp)
     &then
            do i=1,ndim
              beta(i)=betacp(i)
            enddo
            sumsquare=sumsquarecp
          endif
        endif
      else
!local search
        call nongradopt(ndim,funkmin_UnivPhotoFit,f1dim_UnivPhotoFit,
     &beta,betamin,betamax,ftol,sumsquare)
        call funkmin_UnivPhotoFit(ndim,beta,sumsquare)
      endif
      do i=1,ndim
        beta0(i)=beta(i)
      enddo
      sumsquare0=sumsquare
      iderivative=0
      if(ifitmode.lt.0)then
        iwrong=0
      else
        iwrong=1
      endif
      k=ifitmode
      ifitmode=-1
!ifitmode: =-2, ordinary fitting with pco2i calculated as a function of anet
!ifitmode: =-1, ordinary fitting with anet calculated as a function of pco2i
!ifitmode: =1, orthogonal fitting with anet calculated as a function of pco2i
!ifitmode: =2, orthogonal fitting with pco2i calculated as a function of anet
!for odr_leastsquare, only the predicted value of the response variable is needed, i.e., the cost function value is not needed.
!also, only anet as a function of CO2i is considered (not the other way around) because odr_leastsquare cannot handle the situation
!co2i as a function of anet for tpu limitation when alpha=0
      i=1
      if(ntotphips2.ge.1)i=2
      j=4
      if(Currentiknowlimit.eq.-1)then
!fluorescence only fit. chlflphisi2 becomes a forcing variable
        i=1
        j=5
      endif
      call odr_leastsquare(ndim,FCN_UnivPhotoFit,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)
      ifitmode=k
!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)
        po2i(i)=po2i_ori(i)
        chlflphips2(i)=chlflphips2_ori(i)
      enddo
      call funkmin_UnivPhotoFit(ndim,beta,sumsquare)
      if(isitnaninf(sumsquare).eq.1.or.sumsquare.gt.sumsquare0)then
        do i=1,ndim
          beta(i)=beta0(i)
        enddo
        sumsquare=sumsquare0
      endif
      j=0
100   jnon=0
105   sumsquare0=sumsquare
      do i=1,ndim
        beta0(i)=beta(i)
      enddo
      call nongradopt(ndim,funkmin_UnivPhotoFit,f1dim_UnivPhotoFit,
     &beta,betamin,betamax,ftol,sumsquare)
      call funkmin_UnivPhotoFit(ndim,beta,sumsquare)
      if(jnon.le.2.and.(sumsquare0-sumsquare).gt.ftol)then
        jnon=jnon+1
        goto 105
      endif
      if(sumsquare.eq.sumsquare0)goto 110
      if(isitnaninf(sumsquare).eq.1.or.sumsquare.gt.sumsquare0)then
        do i=1,ndim
          beta(i)=beta0(i)
        enddo
        sumsquare=sumsquare0
      endif
      sumsquarecp=sumsquare
      do i=1,ndim
        betacp(i)=beta(i)
      enddo
      call RepeatCompassSearch(ndim,betacp,sumsquarecp,betamin,
     &betamax,funkmin_UnivPhotoFit,f1dim_UnivPhotoFit,xtol)
      call funkmin_UnivPhotoFit(ndim,betacp,sumsquarecp)
      if(sumsquare.eq.sumsquarecp)goto 110
      if(isitnaninf(sumsquarecp).eq.0.and.sumsquarecp.lt.sumsquare)then
        do i=1,ndim
          beta(i)=betacp(i)
        enddo
        sumsquare=sumsquarecp
      endif
      j=j+1
      if(j.le.2.and.(sumsquare0-sumsquare).gt.ftol)goto 100
!
!------------------------------------------------------
110   call funkmin_UnivPhotoFit(ndim,beta,sumsquare)
      return
      END subroutine DoUnivPhotoFit