      subroutine funkmin_UnivPhotoFit(ndim,beta,fvalue)
      implicit none
      include '../testarea/LeafGasParams.h'
      include '../testarea/LeafGasHybridFit.h'
      integer ndim
      double precision beta(1:ndim),fvalue
!(in)  ndim: the dimension of the parameter vector
!(in)  beta: the parameters       
!(out) fvalue: the value of the cost function at beta
!
!---------Local variables--------------------------------------------------
      integer i,n,ilimit0,nummismatch
      double precision pointfvalue
!----------- End of variables declaration ---------------------------------
!check to see if parameters are out of bounds.
      if(isitbounded.eq.1)then
        do i=1,ndim
          if(beta(i).lt.betamin(i).or.beta(i).gt.betamax(i))then
! parameter out of bound          
            fvalue=1.0d+100
            return
          endif
        enddo
      endif
      n=0
      do i=1,ntotunivparams
!replace the values in univparams with those optimized
        if(ifixunivparams(i).eq.1)then
          n=n+1
          univparams(i)=beta(n)
        endif
      enddo
      call UnivParamsAlloc(2)
      ilimit0=Currentilimittype
      fvalue=0.0d0
      nummismatch=0
      do i=1,ntotsamples
        if(Currentilimittype.le.4.and.Currentiknowlimit.eq.1)
     &ilimit0=Currentiphotolimit(i)+4
        call leafunivphotosyn(Currentiknowlimit,ilimit0,ifitmode,
     &aPPFDlf(i),templeaf(i),pco2i(i),po2i(i),chlflphips2(i),
     &anet_obs(i),weitresponses(i:i,1:1),weitresponses(i:i,1:1),
     &weitresponses(i:i,2:2),weitresponses(i:i,1:1),
     &pco2i_pred(i),anet_pred(i),Postiphotolimit(i),pco2c(i),
     &PhiPSII_pred(i),anet_pred_flu(i),pco2i_pred_flu(i),
     &pco2c_anet_flu(i),pco2c_pco2i_flu(i),pointfvalue)
       if(pco2c(i).lt.0.0d0.and.Currentiknowlimit.ne.-1)then
         fvalue=1.0d+101
         return
       endif
       fvalue=fvalue+pointfvalue
      if(Currentiknowlimit.eq.2.and.Currentiphotolimit(i).ne.
     &Postiphotolimit(i))nummismatch=nummismatch+1
      enddo
      if(nummismatch.ne.0)then
!penalize inadmissible fit
        fvalue=fvalue*(dble(nummismatch)*1000.0d0)**2+
     &dble(nummismatch)*1000.0d0
      endif
      return
      end subroutine funkmin_UnivPhotoFit
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$      
      double precision function f1dim_UnivPhotoFit(x)
      implicit none
      double precision x
CU    USES funkmin_UnivPhotoFit
      INTEGER j
!((((((((((((((((((((((((((((((((((((((((((((((((((((
      integer NMAX,ncom
      parameter(NMAX=1000)
      double precision pcom(NMAX),xicom(NMAX)
      COMMON /f1com/ pcom,xicom,ncom
      save /f1com/
!))))))))))))))))))))))))))))))))))))))))))))))))))))
      double precision xt(NMAX)
      do 11 j=1,ncom
        xt(j)=pcom(j)+x*xicom(j)
11    continue
      call funkmin_UnivPhotoFit(ncom,xt,f1dim_UnivPhotoFit)
      return
      END
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      SUBROUTINE FCN_UnivPhotoFit(N,M,NP,NQ,
     +               LDN,LDM,LDNP,
     +               BETA,XPLUSD,
     +               IFIXB,IFIXX,LDIFX,
     +               IDEVAL,F,FJACB,FJACD,
     +               ISTOP)
      implicit none
      include '../testarea/LeafGasParams.h'
      include '../testarea/LeafGasHybridFit.h'
C  SUBROUTINE ARGUMENTS
C      ==> N        NUMBER OF OBSERVATIONS
C      ==> M        NUMBER OF COLUMNS IN EXPLANATORY VARIABLE
C      ==> NP       NUMBER OF PARAMETERS
C      ==> NQ       NUMBER OF RESPONSES PER OBSERVATION
C      ==> LDN      LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING N
C      ==> LDM      LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING M
C      ==> LDNP     LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING NP
C      ==> BETA     CURRENT VALUES OF PARAMETERS
C      ==> XPLUSD   CURRENT VALUE OF EXPLANATORY VARIABLE, I.E., X + DELTA
C      ==> IFIXB    INDICATORS FOR "FIXING" PARAMETERS (BETA)
C      ==> IFIXX    INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X)
C      ==> LDIFX    LEADING DIMENSION OF ARRAY IFIXX
C      ==> IDEVAL   INDICATOR FOR SELECTING COMPUTATION TO BE PERFORMED
C     <==  F        PREDICTED FUNCTION VALUES
C     <==  FJACB    JACOBIAN WITH RESPECT TO BETA
C     <==  FJACD    JACOBIAN WITH RESPECT TO ERRORS DELTA
C     <==  ISTOP    STOPPING CONDITION, WHERE
C                     0 MEANS CURRENT BETA AND X+DELTA WERE
C                       ACCEPTABLE AND VALUES WERE COMPUTED SUCCESSFULLY
C                     1 MEANS CURRENT BETA AND X+DELTA ARE
C                       NOT ACCEPTABLE;  ODRPACK SHOULD SELECT VALUES 
C                       CLOSER TO MOST RECENTLY USED VALUES IF POSSIBLE
C                    -1 MEANS CURRENT BETA AND X+DELTA ARE
C                       NOT ACCEPTABLE;  ODRPACK SHOULD STOP

C  INPUT ARGUMENTS, NOT TO BE CHANGED BY THIS ROUTINE:
      INTEGER          I,IDEVAL,ISTOP,L,LDIFX,LDM,LDN,LDNP,M,N,NP,NQ
      DOUBLE PRECISION BETA(NP),XPLUSD(LDN,M)
      INTEGER          IFIXB(NP),IFIXX(LDIFX,M)
C  OUTPUT ARGUMENTS:
      DOUBLE PRECISION F(LDN,NQ),FJACB(LDN,LDNP,NQ),FJACD(LDN,LDM,NQ)
      integer k
      double precision fvalue
c     
      ISTOP=0
      do I=1,NP
        if(BETA(I).lt.betamin(I).or.
     &BETA(I).gt.betamax(I))then
          ISTOP=1
          return
        endif
      enddo
      do I=1,N
        pco2i(I)=XPLUSD(I,1)
        aPPFDlf(I)=XPLUSD(I,2)
        templeaf(I)=XPLUSD(I,3)
        po2i(I)=XPLUSD(I,4)
        if(Currentiknowlimit.eq.-1)chlflphips2(I)=XPLUSD(I,M)
      enddo
      IF (MOD(IDEVAL,10).GE.1) THEN
        call funkmin_UnivPhotoFit(NP,BETA,fvalue)
        if(fvalue.gt.1.0d+20)then
          ISTOP=1
          return
        endif
        DO 100 I = 1,N
          if(Currentiknowlimit.eq.-1)then
            F(I,1)=anet_pred_flu(I)
          else
            F(I,1)=anet_pred(I)
          endif
  100   CONTINUE
        if(NQ.eq.2)then
          DO 110 I = 1,N
            F(I,NQ)=PhiPSII_pred(I)
  110     CONTINUE
        endif
      END IF
      RETURN
      END
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      double precision function ff_pikaia(ndim,beta01)
      implicit none
      include '../testarea/LeafGasParams.h'
      integer ndim,i
      double precision beta01(ndim),beta(ndim),fvalue

      do i=1,ndim
!        beta01(i)=(beta(i)-betamin(i))/(betamax(i)-betamin(i))
        beta(i)=betamin(i)+beta01(i)*(betamax(i)-betamin(i))       
      enddo
      call funkmin_UnivPhotoFit(ndim,beta,fvalue)
      ff_pikaia=1.0d0/(fvalue+0.00001d0)
      return
      end
