      PROGRAM XYCORR2

C  Same as XYCORR, except certain fit parameters are held fixed.

C **************Compile with option -lodr ************************

C  SET PARAMETERS FOR MAXIMUM PROBLEM SIZE HANDLED BY THIS DRIVER
C  WHERE  MAXN IS THE MAXIMUM NUMBER OF OBSERVATIONS ALLOWED 
C         MAXM IS THE MAXIMUM NUMBER OF COLUMNS IN THE
C              INDEPENDENT VARIABLE ALLOWED
C        MAXNP IS THE MAXIMUM NUMBER OF FUNCTION PARAMETERS
C              ALLOWED 
C          LDX IS THE LEADING DIMENSION OF ARRAY X
C       LDSCLD IS THE LEADING DIMENSION OF ARRAY SCLD
C         LDWD IS THE LEADING DIMENSION OF ARRAY WD
C        LDIFX IS THE LEADING DIMENSION OF ARRAY IFIXX
C        LWORK IS THE DIMENSION OF VECTOR WORK
C       LIWORK IS THE DIMENSION OF VECTOR IWORK

C...PARAMETERS
      INTEGER
     +   MAXN,MAXM,MAXNP,LDSCLD,LDIFX,LDWD,LWORK,LIWORK,PARFIX
      PARAMETER
     +   (MAXN=2000,
     +   MAXM=1,
     +   MAXNP=3,
C     +   LDX=MAXN,
     +   LDSCLD=1,
     +   LDWD=1,
     +   LDIFX=1,
     +   LWORK=17 + 7*MAXN + 10*MAXN*MAXM + 2*MAXN*MAXNP + 8*MAXNP,
     +   LIWORK=19 + 2*MAXNP + MAXM)

C  DECLARE USER-SUPPLIED SUBROUTINES AND
C  ALL OTHER NECESSARY VARIABLES AND ARRAYS

C...LOCAL SCALARS
      DOUBLE PRECISION
     +   PARTOL,SSTOL,TAUFAC
      INTEGER
     +   I,INFO,IPRINT,J,JOB,LDX,LUNERR,LUNRPT,M,MAXIT,N,NDIGIT,NP

C...LOCAL ARRAYS
      DOUBLE PRECISION
     +   BETA(MAXNP),WD(LDWD,MAXM),SCLB(MAXNP),
     +   SCLD(LDSCLD,MAXM),W(MAXN),WORK(LWORK),X(MAXN,MAXM),Y(MAXN)
      INTEGER
     +   IFIXB(MAXNP),IFIXX(LDIFX,MAXM),IWORK(LIWORK)

C   The data in 'xycorr.dat' are in the format: time (or number), xdata, ydata
      DOUBLE PRECISION TIME(MAXN)

C...EXTERNAL SUBROUTINES
      EXTERNAL
     +   DODRC,FUN,JAC
 
C  Number of columns of data in the independent variable matrix X
C  Should be <= MAXM
      DATA M/1/
C  Invoke default option for scale values of errors in independent variable.
      DATA SCLD(1,1)/-1./
C  Invoke default option for weights W.
      DATA W(1)/-1./

      OPEN(UNIT=11,FILE='/users/pieter/odr/xycorr.dat')
C  Read observations from file. N is number of points, NP is order of function
C  to be fitted (1 for linear, 2 for quadratic, etc.).
C  WD(1,1)=DELTA weights (relative size of errors in Y and X coordinates)
C  Units of CO and CH4 are in ppbv, of CO2 in ppmv.

      READ (11,*) N,NP,ERRX,ERRY,PARFIX
      IF(PARFIX .EQ. 0) STOP 'XYCORR2 expects PARFIX=1'
      NP=NP+1
      WD(1,1)=ERRY/ERRX
      IF(N .GT. LDX) STOP 'Increase array sizes for observations'
      IF(NP .GT. MAXNP) STOP 'Increase parameter MAXNP'
      READ (11,*) (BETA(J),J=1,NP)
      READ (11,*) (SCLB(J),J=1,NP)
      READ (11,*) (IFIXB(J),J=1,NP)
      READ (11,*) (TIME(I),X(I,1),Y(I),I=1,N)
      CLOSE(UNIT=11)
c      WRITE(6,*) 'INPUT DATA:'
c      DO 100 I=1,N
c        WRITE(6,'(F8.4,2F8.3)') TIME(I),X(I,1),Y(I)
c100   CONTINUE

      IFIXX(1,1) = -1

C  SET CONTROL VALUES AND STOPPING CRITERIA
      JOB = -1
      NDIGIT = -1
      TAUFAC = -1.0D0
      SSTOL = -1.0D0
      PARTOL = -1.0D0
      MAXIT = -1
      IPRINT = -1
      LUNERR = -1
      LUNRPT = 12
      OPEN(UNIT=LUNRPT,FILE='/users/pieter/odr/report')

C  COMPUTE ODR SOLUTION USING USER-SUPPLIED ANALYTIC DERIVATIVES

      CALL DODRC
     +   (FUN,JAC,
     +   N,M,NP,
     +   X,LDX,IFIXX,LDIFX,SCLD,LDSCLD,
     +   Y,
     +   BETA,IFIXB,SCLB,
     +   WD,LDWD,W,
     +   JOB,NDIGIT,TAUFAC,
     +   SSTOL,PARTOL,MAXIT,
     +   IPRINT,LUNERR,LUNRPT,
     +   WORK,LWORK,IWORK,LIWORK,
     +   INFO)

      END
      SUBROUTINE FUN(N,NP,M,BETA,XPLUSD,LDXPD,F,ISTOPF)

C  Computes predicted values, F, of the dependent variable, given the current
C  values of the independent variable, XPLUSD=X+DELTA, and the function
C  parameters, BETA.  The subroutine argument list and declaration statements
C  must be exactly as shown below.

C  A polynomial of degree NP-1 is calculated

C  Input arguments
C  (which must not be changed by this routine)
C     INTEGER N,NP,M,LDXPD
C     DOUBLE PRECISION BETA(NP),XPLUSD(LDXPD,M)
C  Output arguments
C     DOUBLE PRECISION F(N)
C     INTEGER ISTOPF

      INTEGER
     +   ISTOPF,LDXPD,M,N,NP

      DOUBLE PRECISION
     +   BETA(NP),F(N),XPLUSD(LDXPD,M)

      INTEGER
     +   I

 
      DO 20 I = 1, N
        F(I)=BETA(1)
        DO 10 J = 2, NP
          F(I) = F(I)+BETA(J)*XPLUSD(I,1)**(NP-1)
   10   CONTINUE
   20 CONTINUE
      ISTOPF = 0

      RETURN
      END
