C-----------------------------------------------------------------------
      program runocc
C-----------------------------------------------------------------------
c
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Program Revision 2.01, 13 Jan 1993
c
c runocc is the calling program for OCCAM v2.0. 
C
C REFERENCES: CONSTABLE, PARKER & CONSTABLE, 1987: GEOPHYSICS 52, 289-300.
C             DEGROOT-HEDLIN & CONSTABLE, 1990: GEOPHYSICS 55, 1613-1624.
C             CONSTABLE, 1991: GEOPHYS. J. INT. 106, 387-388.
C             CONSTABLE, 1992: OCCAM DISTRIBUTION NOTES.
C
C REVISION HISTORY:
C     MARCH 1986,                                         (VERSION 1.2)
C     OCTOBER 1987,                                       (VERSION 1.3)
C     AUGUST 1988,                                        (VERSION 1.4)
C     JANUARY 1989 (WITH THANKS TO C.DEGROOT-HEDLIN),     (VERSION 1.5)
C     SEPTEMBER 1989,                                     (VERSION 1.5.2)
C     APRIL 1992 (MORE THANKS TO CATHERINE AND ALSO WSE)  (VERSION 2.0)
C
C IF YOU OBTAIN THIS CODE FROM A THIRD PERSON, PLEASE SEND YOUR NAME AND ADDRESS
C TO S. CONSTABLE.  YOU WILL THEN RECEIVE UPDATES, NEWS ON BUGS, ETC.
C
C VERSION 2.0 IS A RE-WRITE TO MAKE OCCAM INDEPENDENT OF THE DIMENSIONALITY
C OF THE FORWARD PROBLEM, AND TO MAKE MAXIMUM USE OF DYNAMIC MEMORY ALLOCATION, 
C AS THE ARRAY SIZES FOR THE 2D PROBLEM ARE NOW GETTING TOO LARGE FOR COMFORT.
C INCLUDE FILES ARE NOW USED TO STREAMLINE CHANGES OF DIMENSION.
C IMPLICIT DOUBLE PRECISION HAS BEEN REMOVED, AGAIN FOR STORAGE REASONS.
C CALLING PROGRAM IS NOW AN INTEGRAL PART OF THE PACKAGE AS ALL THE MODEL-
C DEPENDENT STUFF HAS BEEN SPUN OFF.
C
C IF YOU OBTAIN THIS CODE FROM A THIRD PERSON, PLEASE SEND YOUR NAME AND ADDRESS
C TO S. CONSTABLE.  YOU WILL THEN RECEIVE UPDATES, NEWS ON BUGS, ETC.
C
C SUBROUTINES WHICH MUST BE SUPPLIED BY THE USER:
C
C  FORMOD(NP,ND,PM,DP,DM), COMPUTES THE FORWARD FUNCTION FOR MODEL PM() AT
C    THE DATA PARAMETERS DP() AND RETURNS IT IN DM().
C  FORDIV(NP,ND,PM,DP,DM,WJ), COMPUTES THE FORWARD FUNCTION AS DOES FORMOD()
C    BUT ALSO RETURNS THE MATRIX OF PARTIALS, WJ(ND,NP).
C  OBJMAT(IRUF, DEL, NPEN), ASSEMBLES THE PENALTY MATRIX (THIS WILL BE DEPENDENT ON MODEL 
C    TYPE AND DIMENSION).  OBJMAT IS ALSO RESPONSIBLE FOR SETTING ANY
C    WEIGHTS, /RESULT/PREWTS(), ASSOCIATED WITH THE MODEL PREJUDICE, ALTHOUGH
C    THESE DEFAULT TO ZERO IF NOT SET.
c  DATIME(DATETM), RETURNS AN 80 CHARACTER STRING WITH CURRENT DATE AND TIME.  CAN 
C    BE A DUMMY ROUTINE IF SYSTEM PROGRAMMING IS TOO MUCH.
C  INPUTD(FILENAME), TAKES THE NAMED FILE AND USES IT TO GET THE DATA IMAGE INTO 
C    THE COMMON BLOCK /DATA/.  SEE OCCAM HEADER.
C  INPUTM(FILENAME), TAKES THE NAMED FILE AND READS IN ANY FIXED ATTRIBUTES
C    OF THE MODEL, SUCH AS LAYER THICKNESSES, FINITE ELEMENT MESH, ETC.  COMMUNICATION
C    BETWEEN INPUTM() AND FORMOD() AND FORDIV() MUST BE VIA THE USER'S OWN COMMON
C    BLOCKS.  WITH STATIC COMMONS (CHECK YOU COMPILER!) THESE NEED NOT BE DECLARED
C    IN THE CALLING PROGRAMS.  INPUTM IS ALSO RESPONSIBLE FOR SETTING ANY
C    MODEL PREJUDICE IN /RESULT/ PREMOD(), ALTHOUGH THIS DEFAULTS TO ZERO
C    IF NOT SET.
c
c Files to be supplied by the user:
c
c  STARTUP must be supplied, in the correct format, to provide inversion parameters.
c  data and model files pointed at by STARTUP
c
c Files generated by OCCAM:
c
c  LOGFIL provides a blow by blow account of the inversion
c  ITERxx provides the inversion parameters, including the model, at every iteration.
c     ITER00 is a copy of the startup file at iteration 0.  Any ITER file can be re-
c     named startup to continue an inversion.
c
c As the response of the model in dm() is not necessarily the response of the final
c model found, but rather the response of the last model tested in the univariate
c optimization routines, we leave it to the user to find the response of
c the final model, using low cunning or the program 'FindResponse.f'.
c
c includes:
      include 'imp.inc'
      include 'occamdim.inc'
c   uses npp
      include 'occam.inc'
c   uses np, pm(), idebug, iout
c
c local variables:
      integer lerr, maxitr, iruf, nit, i, itmax, konv, ifftol
c   lerr = error flag for file operations
c   maxitr = maximum number of subsequent iterations by this run
c   iruf = roughness measure, 0 (size penalty), 1 (1st deriv.), 2 (2nd deriv.)
c   nit = current iteration number (0=start)
c   i = loop index
c   itmax = maximum iteration number
c   konv = convergence flag from OCCAM
c   ifftol = convergence info from OCCAM (has the desired misfit been achieved?)
      real tolreq, tobt, pmu, stepsz, rlast
c   tolreq = desired misfit, in RMS standard errors
c   tobt = current misfit
c   pmu = current Lagrange multiplier
c   stepsz = current RMS step though model space
c   rlast = last RMS measure of roughness
      character*80 itform, descr, modfil, datfil
c   itform = format number of iteration files
c   descr = description line in iteration files
c   modfil = model file name
c   datfil = data file name
c
c parameters
      integer iof1, iof2
      parameter (iof1 = 13, iof2 = 15)
c   iof1 = unit number for the LOGFIL (passed to /result/iout)
c   iof2 = unit number for reading and writing ITER files
c
c calls:
c     occam, filerr, itrout, itrin, inputm, inputd
c      
C-----------------------------------------------------------------------
c Open startup file and read stuff.
c (Startup and iteration files should be machine written, and so we do
c not necessarily test the validity of the input data unless arrays are 
c involved (np and pm)).
c
      open (iof2, file='startup', status='old', iostat=lerr)
      if (lerr .ne. 0) then
        write(*,*) ' Error opening startup file'
        stop
      end if
      call itrin (itform,descr,modfil,datfil,maxitr,tolreq,iruf,
     *    nit,pmu,rlast,tobt,ifftol,iof2)      
      close(iof2)
c            
c if this is zeroth iteration copy the startup file into ITER00
      if (nit .eq. 0) then
c initialize the inversion parameters
        pmu = 5.0
        rlast = 1.0e+10
        tobt = 1000.0
        ifftol = 0
        call itrout (itform,descr,modfil,datfil,maxitr,tolreq,
     *    iruf,nit,pmu,rlast,tobt,ifftol,iof2)      
      end if
c
c zero out the model prejudice and weights in case inputm() and objmat()
c do not take care of setting these.
      do 2 i = 1, np
        prewts(i) = 0.0
2        premod(i) = 0.0
c
      itmax = maxitr + nit
c
c input the data
      if (abs(idebug) .eq. 1) write(*,*) 'reading data..'
      call inputd(datfil)
c input the fixed model characteristics
      if (abs(idebug) .eq. 1) write(*,*) 'reading model..'
      call inputm(modfil)
c
      iout = iof1
      open (iout, file='LOGFIL', iostat=lerr)
      if (lerr .ne. 0) then
        write(*,*) ' Error opening log file'
        stop
      end if
      write(*,*) 'Major I/O accomplished:'
      write(*,*) '  sit back and relax.......'
c
c do it
5      call occam(iruf,tolreq,tobt,ifftol,nit,pmu,rlast,stepsz,konv)
c
c write an iteration file
      call itrout (itform,descr,modfil,datfil,maxitr,tolreq,iruf,
     *    nit,pmu,rlast,tobt,ifftol,iof2)      
c
c let the console know roughly what is going on
      write(*,*) ' **ITERATION ', nit
      write(*,*) ' RMS MISFIT = ', tobt
      write(*,*) ' STEPSIZE = ', stepsz
      write(*,*) ' '
c
c test for convergence: continue inversion if 
c   1) iteration number does not exceed maximum, and
c   2) stepsize is still significant, and
c   3) rms misfit larger than required, and 
c   4) there are no irrecoverable convergence problems
c
      if (konv .ne. 0) then
        write (iof1,*) ' Convergence problems of type ', konv
        close(iof1)
        stop
      else if (nit .ge. itmax) then
        write (iof1,*) ' Stop on maximum iterations'
        close(iof1)
        stop
      else if ((tobt .ge. 1.01*tolreq) .or. (stepsz .ge. 0.0001)) then
        goto 5
      else
        write (iof1,*) ' Stop on normal convergence'
        close(iof1)
        stop
      end if
c      
c      
      end
c
C-----------------------------------------------------------------------
      SUBROUTINE OCCAM(IRUF,TOLREQ,TOBT,IFFTOL,NIT,PMU,RLAST,STEPSZ,
     * KONV)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.01, 20 Jan 1993
C
C OCCAM EXECUTES ONE ITERATION OF A SMOOTH MODEL FINDER
C
C ON INPUT:
C
C  INCLUDE FILE 'OCCAM.INC':
C   /DATA/ CONTAINS
C     ND = NUMBER OF DATA
C     D(NDD) = VECTOR OF DATA VALUES
C     SD(NDD) = VECTOR OF ABSOLUTE STANDARD ERRORS IN THE DATA
C     DP(NDD,4) = MATRIX CONTAINING UP TO FOUR PARAMETERS ASSOCIATED WITH
C        EACH DATUM (RANGE, FREQUENCY, DATA TYPE ETC)
C     NPD = NUMBER OF PARAMETERS ASSOCIATED WITH EACH DATUM (E.G. FOR FREQ. ONLY
C              NPD = 1, FOR FREQ. AND RANGE NPD = 2)
C
C   /MODEL/ CONTAINS
C     NP = NUMBER OF MODEL PARAMETERS, USUALLY THE NUMBER OF LAYERS
C     PM(NPP) = VECTOR OF INITIAL MODEL PARAMETERS, USUALLY LOG10(LAYER RES'IES)
C
C   /RESULT/ CONTAINS
C     IDEBUG = 0 NORMAL, FASTEST AND QUIETEST EXECUTION
C            = 1 A LOW LEVEL OF VERBOSITY LETS YOU KNOW HOW THINGS ARE GOING
C            = 2 A TABLE OF RMS MISFIT VERSUS LAGRANGE MULTIPLIER IS PRODUCED
C        (USEFUL FOR TRACKING DOWN CONVERGENCE PROBLEMS).  THIS TABLE IS USED 
C        TO BRACKET THE MINIMUM
C            = 3 IF ONE ALSO WANTS A DUMP OF THE MODEL AT EACH LAGRANGE MULT IN
C        THE TABLE (NOT LIKELY TO BE REQUIRED BY END USERS)
c        if idebug is negative, the above absolute values are operational but
c        the step-reduction algorithm is disabled.
C     IOUT = FORTRAN UNIT NUMBER FOR OUTPUT, CAN EITHER BE TERMINAL OR A FILE
C
C  ARGUMENTS:
C    IRUF = 0 FOR MODEL MINIMIZATION
C         = 1 FOR FIRST DERIVATIVE MINIMIZATION
C         = 2 FOR SECOND DERIVATIVE MINIMIZATION
C     TOLREQ = REQUIRED RMS MISFIT
C     TOBT = RMS MISFIT FROM LAST ITERATION (LARGE INITIALLY)
C     IFFTOL = RECORD OF WHETHER FEASIBLE MODEL HAS BEEN ATTAINED (0=NO,
C         1=YES).  INITIALLY SET 0 AND THEN LEFT ALONE.
C     NIT = THE NUMBER OF PREVIOUS ITERATIONS IN THIS INVERSION. SET TO 0 ON FIRST 
C         CALL, IT WILL BE UPDATED BY OCCAM. 
C     PMU = LOG10(LAGRANGE) MULTIPLIER TO START SEARCHES WITH.  INITIALLY 
C         SET LARGE AND THEN LEFT FROM LAST CALL FOR EFFICIENCY
C     RLAST = ROUGHNESS OF LAST MODEL.  INITIALLY SET LARGE AND THEN LEFT.        
C
C ON OUTPUT:
C
C  /MODEL/ CONTAINS
C     PM(NPP) = VECTOR OF UPDATED MODEL PARAMETERS
C     DM(NDD) = VECTOR CONTAINING *APPROXIMATE* RESPONSE OF NEW MODEL
C  ARGUMENTS:
C     TOBT = RMS MISFIT OF NEW MODEL RESPONSE.
C     NIT = NUMBER OF PREVIOUS CALLS TO OCCAM1 DURING THIS INVERSION
C     STEPSZ = RMS MEASURE OF THE CHANGES IN THE MODEL PARAMETERS 
C     KONV = STATUS FLAG:
C        0 = NORMAL EXIT FOR A SUCCESSFUL ITERATION,
C        1 = A PERFECTLY SMOOTH MODEL HAS BEEN FOUND FOR THE REQUIRED TOLERANCE
C        2 = CONVERGENCE PROBLEMS.  UNABLE TO FIND A MODEL WITH AN R.M.S. MISFIT
C          LESS THAN OR EQUAL TO THAT OF THE LAST ITERATION'S
C        3 = CONVERGENCE PROBLEMS.  UNABLE TO FIND A MODEL WITH AN R.M.S. MISFIT
C          OF TOLREQ AND SMOOTHER THAN THE LAST ITERATION'S
C
C SUBROUTINES REQUIRED AND SUPPLIED:
C
C   MAKPTP, A SUBROUTINE WHICH CREATES THE ROUGHENING MATRIX TIMES ITS TRANS.
C   TOFMU(AMU), A FUNCTION WHICH RETURNS THE RMS MISFIT OF THE RESPONSE
C      OF THE MODEL CONSTRUCTED USING THE LAGRANGE MULTIPLIER AMU.
C      (CALLS CHOLIN,CHOLSL,FORMOD,ANORM)
C   INITM, TRMULT, MULT, ATAMULT AND ANORM, SUBROUTINES FOR MATRIX OPERATIONS
C   FMIN, A SUBROUTINE WHICH MINIMISES A UNIVARIATE FUNCTION
C   FROOT, A SUBROUTINE WHICH FINDS THE ROOT OF A UNIVARIATE FUNCTION
C   MINBRK, A SUBROUTINE WHICH BRACKETS A UNIVARIATE MINIMUM
C   SCANMU, A SUBROUTINE WHICH GENERATES THE MISFIT TABLE WHEN IBUG .GT. 0
C
C SUBROUTINES WHICH MUST BE SUPPLIED BY THE USER:
C
C  FORMOD(NP,ND,PM,DP,DM), COMPUTES THE FORWARD FUNCTION FOR MODEL PM() AT
C    THE DATA PARAMETERS DP() AND RETURNS IT IN DM().
C  FORDIV(NP,ND,PM,DP,DM,WJ), COMPUTES THE FORWARD FUNCTION AS DOES FORMOD()
C    BUT ALSO RETURNS THE MATRIX OF PARTIALS, WJ(ND,NP).
C  OBJMAT(IRUF, DEL, NPEN), ASSEMBLES THE PENALTY MATRIX (THIS WILL BE DEPENDENT ON MODEL 
C    TYPE AND DIMENSION)
C
C
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      INTEGER IRUF,IFFTOL,NIT,KONV
      REAL TOLREQ,TOBT,PMU,RLAST,STEPSZ
C LOCAL VARIABLES:
      INTEGER I
      REAL TOL0, TMIN, TINT
      REAL AMU1, AMU2, AMU3, TAMU1, TAMU2, TAMU3
      REAL RUF, PMU2
C FUNCTIONS
      REAL TOFMU, FMIN, FROOT, ANORM, FNDRUF
      EXTERNAL TOFMU
C
C LOCAL PARAMETERS
      REAL TOLM, TOLI
C TOLM AND TOLI ARE THE TOLERANCES FOR FMIN AND FROOT RESPECTIVELY
C DECREASING THEM WILL IMPROVE ACCURACY AT THE COST OF MORE FORWARD 
C CALCULATIONS
      PARAMETER (TOLM= 0.1, TOLI = 0.001)
C
C-----------------------------------------------------------------------
      IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ENTERING OCCAM..'
C SET THINGS UP FOR THIS ITERATION:
      KONV = 0
C FRAC CONTROLS THE STEP SIZE; NORMALLY WILL REMAIN AT 1.0 UNLESS WE HAVE
C CONVERGENCE PROBLEMS
      FRAC = 1.0
C CONSTRUCT WJTWJ, WJTWD, AND FIND INITIAL MISFIT
      IF (abs(idebug) .EQ. 1) WRITE(*,*) 
     *  ' CONSTRUCTING DERIVATIVE DEPENDENT MATRICES..'
      CALL MAKJTJ(TOL0)
C
      IF (NIT .EQ. 0) THEN
        WRITE(IOUT,*) 'STARTING R.M.S. = ',TOL0
	  if ((TOL0 .le. TOLREQ) .and. (IFFTOL .eq. 0)) then
	    write (*,*) ' This model is feasible'
	    IFFTOL = 1
	  end if
      END IF
      IF (abs(idebug) .EQ. 1) WRITE(*,*) 'STARTING R.M.S. = ',TOL0
      WRITE(IOUT,*) ' '
      WRITE(IOUT,*) '** ITERATION ',NIT+1,' **'
C
C CONSTRUCT THE PENALTY MATRIX MULTIPLIED INTO ITS TRANSPOSE
      IF (abs(idebug) .EQ. 1) WRITE(*,*)' CONSTRUCTING PENALTY MATRIX..'
      CALL MAKPTP(IRUF)
C
C THE NEXT BLOCK OF CODE CONTROLS THE SELECTION OF THE LAGRANGE MULTIPLIER,
C PMU:
      IF (abs(idebug) .EQ. 1) 
     * WRITE(*,*) ' SEARCHING LAGRANGE MULTIPLIER:'
120   NFOR = 0
C PRODUCE THE MISFIT FUNCTION IF REQUIRED (THIS WILL BE USED TO BRACKET MIN)
      IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ..BRACKETING MINIMUM..'
      IF (abs(idebug) .GE. 2) THEN
        CALL SCANMU(AMU1,AMU2,AMU3,TAMU2,TOLREQ)
      ELSE
C BRACKET THE MINIMUM USING MINBRK AND TWO GUESSES
        AMU1 = PMU - 1.0
        AMU2 = PMU
        CALL MINBRK(AMU1,AMU2,AMU3,TAMU1,TAMU2,TAMU3,TOFMU)
      END IF
C FIND THE MINIMUM
      IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ..FINDING MINIMUM..'
      IF (TAMU2 .LT. TOLREQ) THEN
C WE'VE BEEN LUCKY AND FOUND AN ACCEPTABLE MINIMUM USING MINBRK
        TMIN = TAMU2
        PMU = AMU2
        WRITE(IOUT,*) 'MINIMUM TOL FROM MINBRK IS AT MU =',PMU
      ELSE
        TMIN = FMIN(AMU1,AMU2,AMU3,TAMU2,TOFMU,TOLM,PMU)
        WRITE(IOUT,*) 'MINIMUM TOL FROM FMIN IS AT MU =',PMU
      END IF
      WRITE(IOUT,*) 'AND IS =',TMIN
      WRITE(IOUT,*) 'USING ',NFOR,' EVALUATIONS OF FUNCTION'
C IF THE NEW MINIMUM TOLERANCE IS GREATER THAN THE TOLERANCE FROM THE
C PREVIOUS MODEL, WE ARE HAVING CONVERGENCE PROBLEMS. CUT THE STEP SIZE
      IF (((TMIN .GT. TOL0) .AND. (IFFTOL .EQ. 0) .OR.
     *    (TMIN .GT. TOLREQ) .AND. (IFFTOL .EQ. 1)) .and. 
     *    idebug .ge. 0) THEN
        WRITE(IOUT,*) 'DIVERGENCE PROBLEMS, CUTTING STEP SIZE'
        IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ..DIVERGENCE PROBLEMS..'
        FRAC = FRAC*2.0
        IF (FRAC .GT. 1.0E+05) THEN
C WE HAVE CUT STEP SIZE A LOT TO NO EFFECT: GIVE UP
          KONV = 2
          RETURN
        ELSE
          GOTO 120
        END IF
      END IF
C
      IF (TMIN .LT. TOLREQ) THEN
        IF (abs(idebug) .EQ. 1) WRITE(*,*) ' ..FINDING INTERCEPT..'
C TOLERANCE IS BELOW THAT REQUIRED; FIND INTERCEPT.
C THE LOWER VALUE OF MU BRACKETING THE INTERCEPT IS EASILY FOUND: IT IS
C JUST THE MINIMUM
        AMU1 = PMU
        TAMU1 = TMIN
C THE UPPER BOUND IS FOUND BY TESTING SUCCESSIVELY GREATER MU
        AMU2 = PMU
        TAMU2 = TMIN
        NFOR = 0
140     IF (TAMU2 .LT. TOLREQ) THEN
C SUCCESSIVELY DOUBLE MU:
          AMU2 = AMU2 + 0.30103
          TAMU2 = TOFMU(AMU2)
          GOTO 140
        END IF
        PMU2 = FROOT(TOFMU,AMU1,AMU2,TAMU1,TAMU2,TOLREQ,TOLI)
        TINT = TAMU2 + TOLREQ
        WRITE(IOUT,*) 'INTERCEPT IS AT MU = ',PMU2
        WRITE(IOUT,*) 'AND IS = ',TINT
        WRITE(IOUT,*) 'USING ',NFOR,' FUNCTION EVALUATIONS'
      ELSE
C TOLERANCE IS NO YET SMALL ENOUGH.  WE WILL KEEP THE MINIMUM.
C SINCE FMIN RETURNS WITH PWK1() INSTEAD OF PWK2() WE NEED THIS COPY
        DO 150 I = 1,NP
150       PWK2(I) = PWK1(I)
      END IF
C***END LAGRANGE MULTIPLIER SELECTION
C
C COMPUTE ROUGHNESS.  WE DO THIS BY A FUNCTION CALL TO AVOID HAVING THE PENALTY
C MATRIX HANG AROUND TAKING UP MEMORY.
      IF (abs(idebug) .EQ. 1) WRITE(*,*) 'COMPUTING ROUGHNESS..'
      RUF = FNDRUF(IRUF)
C
C IF WE ATTAINED THE INTERCEPT LAST ITERATION BUT THE MODEL IS GETTING 
C ROUGHER WE HAVE PROBLEMS WITH CONERGENCE.
      IF ((RUF.GT.1.01*RLAST) .AND. (IFFTOL.EQ.1) .and. (idebug.ge.0)) 
     * THEN
        WRITE(IOUT,*) 'ROUGHNESS PROBLEMS, CUTTING STEP SIZE'
        FRAC = FRAC*2.0
C CHECK TO SEE IF ALL IS HOPELESS
        IF (FRAC .GT. 1.0E+05) THEN
          KONV = 3
          RETURN
        END IF
C OTHERWISE PLOUGH ON
        GOTO 120
      END IF
C
C SAVE NEW MODEL AND COMPUTE STEP SIZE
      IF (abs(idebug) .EQ. 1) WRITE(*,*) 'COMPUTING STEPSIZE..'
      DO 180 I = 1,NP
        PWK3(I) = PWK2(I) - PM(I)
180     PM(I) = PWK2(I)
C THE STEPSIZE IS THE ACTUAL CHANGE IN THE MODEL, NORMALISED BY NP:
      STEPSZ = SQRT(ANORM(NP,PWK3)/NP)
      WRITE(IOUT,*) 'STEPSIZE IS = ',STEPSZ
      WRITE(IOUT,*) 'ROUGHNESS IS = ',RUF
C TIDY UP
      NIT = NIT + 1
      RLAST = RUF
      IF (TMIN .LT. TOLREQ) THEN
        IFFTOL = 1
        TOBT = TINT
      ELSE
        TOBT = TMIN
      END IF
C SEE IF WE HAVE A PERFECLY SMOOTH MODEL THAT FITS DATA:
      IF (RUF .LT. 1.0E-5 .AND. TOBT .LE. 1.01*TOLREQ) KONV = 1
      IF (abs(idebug) .EQ. 1) WRITE(*,*) 'LEAVING OCCAM'
      RETURN
      END
C
C-----------------------------------------------------------------------
       SUBROUTINE SCANMU(AM1,AM2,AM3,T2,T0)
C-----------------------------------------------------------------------
C 
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
C SCANMU PRODUCES A REPORT OF MISFIT VERSUS LAGRANGE MULTIPLIER.
C IT CALLS THE FORWARD ROUTINE A LOT SO IS REALLY ONLY USED FOR DEBUGGING
C AND PARTICULARLY NASTY CASES
C IT LOOKS FOR THE GLOBAL MINIMUM AS FOLLOWS:
C   STARTING AT THE LARGEST MU USED FOR THE DEBUG REPORT, WE LOOK FOR THE 
C   FIRST MINIMUM.  WE ALSO LOOK FOR THE GLOLBAL MIN.  WE WANT THE GLOBAL 
C   MIN UNLESS THE LOCAL MIN AT LARGER MU IS LOWER THAN THE REQUIRED 
C   TOLERANCE. 
C 
C ON INPUT:
C   T0 IS THE REQUIRED TOLERANCE
C   /RESULT/ IDEBUG IS THE PRINT LEVEL FROM OCCAM1
C   /RESULT/ IOUT IS A UNIT NUMBER FOR OUTPUT
C ON OUTPUT:
C   AM1, AM2, AM3 ARE THREE LAGRANGE MULTIPLIERS WHICH BRACKET A MINIMUM
C   T2 IS THE TOLERANCE AT AM2.
C   /RESULT/ ..PMW1().. HAS THE MODEL ASSOCIATED WITH AM2,T2.
C SUBROUTINES CALLED: 
C   TOFMU
C
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      REAL AM1,AM2,AM3,T2,T0
C LOCAL VARIABLES
      LOGICAL SEARCH
      INTEGER I,J,K
      REAL T1,TG2,T3,TL2,AML3,AML2,AML1,AMG3,AMG2,AMG1,AMU,T
C FUNCTIONS
      REAL TOFMU
C
C-----------------------------------------------------------------------
      WRITE(IOUT,*) 'MISFIT AS A FUNCTION OF MU:'
      T1 = 1.0E+09
      T2 = 1.1E+09
      TG2 = 0.9E+09
      AM2 = 16.
      AM1 = 16.
      SEARCH = .TRUE.
C
C START AT THE TOP (SMOOTHEST AND PRESUMABLY WORST FIT), SWEEPING FROM
C LOG10(MU) = 16 TO LOG10(MU) = -4 AT TWO VALUES PER DECADE
      DO 100 K = 16,-4,-1
C KEEP THE LAST TWO VALUES FOR MINIMUM SEARCH
        T3 = T2
        T2 = T1
        AM3 = AM2
        AM2 = AM1
C KEEP THE MODEL ASSOCIATED WITH THE MIDDLE VALUE IF WE ARE SEARCHING
C FOR THE MINIMUM
        IF (SEARCH) THEN
          DO 10 I = 1,NP
10          PWK1(I) = PWK2(I)
        END IF
C KEEP THE MODEL ASSOCIATED WITH THE MIDDLE VALUE ANYWAY IN CASE
C GLOBAL AND LOCAL MINIMA ARE FOUND
        DO 20 I = 1,NP
20        PWK3(I) = PWK2(I)
C GET THE NEXT MU AND ASSOCIATED TOLERANCE
        AM1 = FLOAT(K)/2.
        T1 = TOFMU(AM1)
C
        IF ((T2 .LE. T1) .AND. (T2 .LE. T3) .AND. SEARCH) THEN
C WE HAVE FOUND THE FIRST MINIMUM. KEEP IT AND TURN OFF SEARCH
C (MODEL IS ALREADY STASHED IN PM1())
          SEARCH = .FALSE.
          TL2 = T2
          AML3 = AM3
          AML2 = AM2
          AML1 = AM1
        END IF
C KEEP AN EYE OUT FOR THE GLOBAL MIN IN CASE IT IS DIFFERENT
        IF (T2 .LE. TG2) THEN
          TG2 = T2
          AMG3 = AM3
          AMG2 = AM2
          AMG1 = AM1
C COPY ACCROSS MODEL ASSOCIATED WITH GLOBAL MIN.
          DO 30 I = 1,NP
30          PWK1(I) = PWK3(I)
        END IF
C RECORD IN THE LOG FILE
        WRITE(IOUT,*) AM1,T1
C
C OUTPUT THE MODEL AS WELL IF REQUIRED
        IF (abs(idebug) .GE. 3) THEN
          DO 60 J = 1,NP
60          WRITE(IOUT,*) '   ',PWK2(J)
        END IF
100   CONTINUE
C WRITE OUT THE GAUSS STEP AT HIGHEST DEBUG LEVEL
      IF (abs(idebug) .GE. 3) THEN
        AMU = -99.0
        T = TOFMU(AMU)
        WRITE(IOUT,*) AM1,T1
        DO 110 J = 1,NP
110       WRITE(IOUT,*) '   ',PWK2(J)
      END IF
C
C NOW CHECK OUT THE MINIMA
C  IF THE FIRST MINIMUM ACHIEVES THE TOLERANCE REQUIRED WE WANT IT
      IF (TL2 .LE. T0) THEN
        T2 = TL2
        AM1 = AML1
        AM2 = AML2
        AM3 = AML3
      ELSE 
C OTHERWISE THE GLOBAL MIN IS JUST FINE
        T2 = TG2
        AM1 = AMG1
        AM2 = AMG2
        AM3 = AMG3
      END IF
      RETURN
      END
C
C-----------------------------------------------------------------------
      REAL FUNCTION TOFMU(ALOGMU)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.01, 13 Jan 1993
C
C FUNCTION TOFMU RETURND THE RMS MISFIT OF THE RESPONSE OF A MODEL CONSTRUCTED
C USING THE GIVEN VALUE OF LAGRANGE MULTIPLIER
C
C ON INPUT:
C   ALOGMU = LOG10 OF THE LAGRANGE MULTIPLIER 
C            (-99. IF A ZERO LAGRANGE MULTIPLIER NEEDED)
C   /RESULT/ WJTWJ(), WJTWD(), PTP(), PREWTS(), PREMOD(), NFOR
C   /DATA/ DP(), D(), SD(), ND
C   /MODEL/ NP
C
C ON OUTPUT
C   TOFMU = R.M.S. MISFIT OR 1.0E+10 IF CHOLESKY DECOMPOSITION FAILED
C   /RESULT/ PWK2() CONTAINS THE MODEL WHICH PRODUCES THE MISFIT TOFMU
C   /MODEL/ DM() CONTIANS THE RESPONSE OF THE MODEL 
C
C SUBROUTINES CALLED:
C   CHOLIN, CHOLSL ARE R.L. PARKER'S SUBROUTINES TO PERFORM CHOLESKY DECOMP.
C       AND THEN SOLVE A LINEAR SYSTEM BY BACK SUBSTITUTION.
C   ANORM, FORMOD ARE EXPLAINED IN OCCAM1
C
C  AT VERSION 1.4 LOG10(MU) USED TO IMPROVE PERFORMANCE OF 1D OPTIMISATION
c
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      REAL ALOGMU
C LOCAL VARIABLES
      REAL A(NPP,NPP+1), DWK1(NDD), AMU, PART, PWK4(NPP)
      INTEGER I, J, ISTAT
C FUNCTIONS:
      REAL ANORM
C
C-----------------------------------------------------------------------
      IF (ALOGMU. EQ. -99.) THEN
C THE GAUSS STEP IS REQUIRED
        AMU = 0.0
      ELSE
        AMU = 10.**(ALOGMU)
      END IF
      NFOR = NFOR + 1
C ADD AMU.TRANS(PARTIAL).PARTIAL TO TRANS(W.J).W.J
C ALSO ADD MODEL PREJUDICE TO WJTWD()
      DO 20 I = 1,NP
        DO 20 J = 1,NP
          PWK4(J) = WJTWD(J) + AMU*PREWTS(J)*PREMOD(J)
20        A(I,J) = AMU*PTP(I,J) + WJTWJ(I,J)
C
C SOLVE THE LINEAR SYSTEM
      CALL CHOLIN(NPP,NP,A,1,ISTAT)
      IF (ISTAT .NE. 0) THEN
C IF DECOMP FAILED WE GIVE WARNING BUT MIGHT STILL GET BY IF TOFMU
C IS SET TO A LARGE VALUE
        WRITE(*,*) 
     *    'Warning:- Cholesky decomp. failed at mu = ',AMU
        TOFMU = 1.0E+10
        RETURN
      END IF
      CALL CHOLSL(NPP,NP,A,PWK2,PWK4,0)
C CUT STEP SIZE IF NECESSARY
      IF (FRAC .GT. 1.001) THEN
        PART = 1./FRAC
        DO 30 I = 1,NP
30        PWK2(I) = (1.-PART)*PM(I) + PART*PWK2(I)
      END IF
C CALCULATE THE MODEL RESPONSE AND MISFIT
      CALL FORMOD(NP,ND,PWK2,DP,DM)
      DO 40 I = 1,ND
40      DWK1(I) = (D(I) - DM(I))/SD(I)
      TOFMU = SQRT(ANORM(ND,DWK1)/ND)
      IF (abs(idebug) .EQ. 1) 
     *  WRITE(IOUT,*) 'TOFMU =', TOFMU,' AMU =', ALOGMU
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      REAL FUNCTION ANORM(N,D)
C-----------------------------------------------------------------------
C RETURNS THE SQUARE OF THE EUCLIDEAN NORM OF A VECTOR
      INCLUDE 'imp.inc'
C
      INTEGER I, N
      REAL D(N)
C
      ANORM = 0.0
      DO 10 I = 1,N
        ANORM = ANORM + D(I)*D(I)
10    CONTINUE
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      SUBROUTINE INITM(MD,M,N,A)
C-----------------------------------------------------------------------
C ZEROS A 2D MATRIX
      INCLUDE 'imp.inc'
C
      INTEGER MD, M, N, I, J
      REAL A(MD,N)
C
      DO 10 I = 1,M
        DO 10 J = 1,N
10        A(I,J) = 0.0
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      SUBROUTINE MULT(MAD,MA,NAD,NA,NB,A,B,C)
C-----------------------------------------------------------------------
C MULTIPLIES TWO 2D MATRICES
      INCLUDE 'imp.inc'
C
      INTEGER MAD, MA, NAD, NA, NB, I, J, K
      REAL A(MAD,NA),B(NAD,NB),C(MAD,NB), CIJ
C
      DO 12 I = 1,MA
        DO 12 J = 1,NB
          CIJ = 0.0
          DO 10 K = 1,NA
10          CIJ = A(I,K)*B(K,J) + CIJ
12        C(I,J) = CIJ
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      SUBROUTINE TRMULT(MAD,MA,NAD,NA,NB,A,B,C)
C-----------------------------------------------------------------------
C MULTIPLIES THE TRANSPOSE OF A 2D MATRIX BY ANOTHER MATRIX
      INCLUDE 'imp.inc'
C
      INTEGER MAD, MA, NAD, NA, NB, I, J, K
      REAL A(MAD,NA),B(MAD,NB),C(NAD,NB), CIJ
      DO 12 I = 1,NA
        DO 12 J = 1,NB
          CIJ = 0.0
          DO 10 K = 1,MA
10          CIJ = A(K,I)*B(K,J) + CIJ
12        C(I,J) = CIJ
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      SUBROUTINE MINBRK(AX,BX,CX,FA,FB,FC,FUNC)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
C MINBRK BRACKETS A UNIVARIATE MINIMUM OF A FUNCTION. TO BE USED PRIOR TO
C A UNIVARIATE MINIMISATION ROUTINE.
C BASED ON A ROUTINE IN NUMERICAL RECIPES BY PRESS ET AL
C MODIFIED SO THAT THE MODEL ASSOCIATED WITH THE MISFITS IS CARRIED AROUND
C FOR USE IN THE MINIMISATION ROUTINES AND POSSIBLY ULTIMATELY KEPT AS THE
C RESULT OF THIS ITERATION;
C
C ON INPUT:
C   AX, BX = TWO DISTINCT ESTIMATES OF THE MINIMUM'S WHEREABOUTS
C   FUNC = THE FUNCTION IN QUESTION
C   /MODEL/ NP = NUMBER OF MODEL PARAMETERS
C ON OUTPUT:
C   AX,BX,CX = THREE NEW POINTS WHICH BRACKET THE MINIMUM
C   FA,FB,FC = FUNC(AX), FUNC(BX) ETC.
C   /RESULT/ PWK1() = THE MODEL ASSOCIATED WITH BX AND FB
C
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      REAL AX,BX,CX,FA,FB,FC
C LOCAL VARIABLES:
      INTEGER I
      REAL DUM,R,Q,ULIM,U,FU
C LOCAL PARAMETERS
      REAL GOLD,GLIMIT,TINY
      PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-32)
C FUNCTIONS:
      REAL FUNC
C
C-----------------------------------------------------------------------
      FB = FUNC(BX)
C MAKE PWK1 ASSOCIATED WITH B
      DO 4 I = 1,NP
4       PWK1(I) = PWK2(I)
C PM2 WILL BE ASSOCIATED WITH A
      FA = FUNC(AX)
      IF (FB.GT.FA) THEN
C SWITCH A AND B
        DUM = AX
        AX = BX
        BX = DUM
        DUM = FB
        FB = FA
        FA = DUM
C KEEP PWK1 ASSOCIATED WITH B (A WILL BE LOST SO DOESN'T MATTER)
        DO 5 I = 1,NP
5         PWK1(I) = PWK2(I)
      ENDIF
      CX = BX + GOLD*(BX - AX)
      FC = FUNC(CX)
C MAIN TEST. IF IT FAILS WE EXIT WITH B
10    IF (FB.GE.FC) THEN
C KEEP PWK1 ASSOCIATED WITH C 
        DO 11 I = 1,NP
11        PWK1(I) = PWK2(I)
        R = (BX - AX)*(FB - FC)
        Q = (BX - CX)*(FB - FA)
        U = BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R))
        ULIM = BX + GLIMIT*(CX - BX)
        IF ((BX - U)*(U - CX).GT.0.) THEN
          FU = FUNC(U)
          IF (FU.LT.FC) THEN
            AX = BX
            FA = FB
            BX = U
            FB = FU
C KEEP PWK1 ASSOCIATED WITH B 
            DO 12 I = 1,NP
12            PWK1(I) = PWK2(I)
            GO TO 10
          ELSE IF (FU.GT.FB) THEN
            CX = U
            FC = FU
            GO TO 10
          ENDIF
          U = CX + GOLD*(CX - BX)
          FU = FUNC(U)
        ELSE IF ((CX - U)*(U - ULIM).GT.0.) THEN
          FU = FUNC(U)
          IF (FU.LT.FC) THEN
            BX = CX
            CX = U
            U = CX + GOLD*(CX - BX)
            FB = FC
            FC = FU
C KEEP PWK1 ASSOCIATED WITH B 
            DO 13 I = 1,NP
13            PWK1(I) = PWK2(I)
            FU = FUNC(U)
          ENDIF
        ELSE IF ((U - ULIM)*(ULIM - CX).GE.0.) THEN
          U = ULIM
          FU = FUNC(U)
        ELSE
          U = CX + GOLD*(CX - BX)
          FU = FUNC(U)
        ENDIF
        AX = BX
        BX = CX
        CX = U
        FA = FB
        FB = FC
        FC = FU
        GO TO 10
      ENDIF
      RETURN
      END
C
C
C
C-----------------------------------------------------------------------
      REAL FUNCTION FMIN(AX,BX,CX,FBX,F,TOL,XMIN)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
C FMIN RETURNS THE MINIMUM VALUE OF A FUNCTION
C BASED ON A ROUTINE IN NUMERICAL RECIPES BY PRESS ET AL 
C MODIFIED SO THAT THE MODEL ASSOCIATED WITH THE MISFITS IS CARRIED AROUND
C FOR USE POSSIBLE ULTIMATE USE AS THE
C RESULT OF THIS ITERATION;
C
C ON INPUT:
C  AX,BX,CX = INDEPENDENT VARIABLE WHICH BRACKET THE MINIMUM
C  FBX = F(BX) (USUALLY AVAILABLE FROM THE BRACKETING PROCEDURE)
C  F = FUNCTION IN QUESTION
C  TOL = FRACTIONAL TOLERANCE REQUIRED IN THE INDEPENDENT VARIABLE
C  /RESULT/ PWK1 = MODEL ASSOCIATED WITH BX AND FBX
C  /MODEL/ NP = NUMBER OF MODEL PARAMETERS
C ON OUTPUT:
C  XMIN = ABSCISSA OF MINIMUM
C  FMIN = F(XMIN)
C  /RESULT/ PWK1 = MODEL ASSOCIATED WITH XMIN AND FMIN
C
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      REAL AX,BX,CX,FBX,TOL,XMIN
C LOCAL VARIABLES
      REAL AA,B,V,W,X,E,FX,FV,FW,TOL1,Q,R,P,ETEMP,DD,FU,U,XM,TOL2
      INTEGER I,ITER
C LOCAL PARAMETERS
      INTEGER ITMAX
      REAL CGOLD,ZEPS
      PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10)
C FUNCTION:
      REAL F
C
C-----------------------------------------------------------------------
      AA = MIN(AX,CX)
      B = MAX(AX,CX)
      V = BX
      W = V
      X = V
      E = 0.
      FX = FBX
      FV = FX
      FW = FX
      DO 11 ITER = 1,ITMAX
        XM = 0.5*(AA + B)
        TOL1 = TOL*ABS(X) + ZEPS
        TOL2 = 2.*TOL1
        IF (ABS(X - XM).LE.(TOL2 - .5*(B - AA))) GOTO 30
        IF (ABS(E).GT.TOL1) THEN
          R = (X - W)*(FX - FV)
          Q = (X - V)*(FX - FW)
          P = (X - V)*Q - (X - W)*R
          Q = 2.*(Q - R)
          IF (Q.GT.0.) P =  - P
          Q = ABS(Q)
          ETEMP = E
          E = DD
          IF (ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(AA - X) .OR.
     *        P.GE.Q*(B - X)) GOTO 1
          DD = P/Q
          U = X + DD
          IF (U - AA.LT.TOL2 .OR. B - U.LT.TOL2) DD = SIGN(TOL1,XM - X)
          GOTO 2
        ENDIF
1       IF (X.GE.XM) THEN
          E = AA - X
        ELSE
          E = B - X
        ENDIF
        DD = CGOLD*E
2       IF (ABS(DD).GE.TOL1) THEN
          U = X + DD
        ELSE
          U = X + SIGN(TOL1,DD)
        ENDIF
        FU = F(U)
        IF (FU.LE.FX) THEN
          IF (U.GE.X) THEN
            AA = X
          ELSE
            B = X
          ENDIF
          V = W
          FV = FW
          W = X
          FW = FX
          X = U
          FX = FU
          DO 5 I = 1,NP
5           PWK1(I) = PWK2(I)
        ELSE
          IF(U.LT.X) THEN
            AA = U
          ELSE
            B = U
          ENDIF
          IF (FU.LE.FW .OR. W.EQ.X) THEN
            V = W
            FV = FW
            W = U
            FW = FU
          ELSE IF (FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN
            V = U
            FV = FU
          ENDIF
        ENDIF
11    CONTINUE
      WRITE(*,*) 'MAXIMUM ITERATIONS EXCEEDED IN FMIN'
30    XMIN = X
      FMIN = FX
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      REAL FUNCTION FROOT(FUNC,X1,X2,FA,FB,OFF,TOL)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
C FROOT FINDS THE POINT AT WHICH A UNIVARIATE FUNCTION ATTAINS A GIVEN
C  VALUE
C BASED ON A ROUTINE IN NUMERICAL RECIPES BY PRESS ET AL 
C MODIFIED SO THAT THE MODEL ASSOCIATED WITH THE MISFITS IS CARRIED AROUND
C FOR USE POSSIBLE ULTIMATE USE AS THE
C RESULT OF THIS ITERATION;
C
C ON INPUT:
C  FUNC = THE FUNCTION IN QUESTION
C  X1,X2 = INDEPENDENT VARIABLES BRACKETING THE ROOT
C  FA,FB = FUNC(X1),FUNC(X2) (USUALLY AVAILABLE FROM THE BRACKETING)
C  OFF = VALUE OF THE FUNCTIONAL AT THE REQUIRED ABSCISSA
C  TOL = FRACTIONAL TOLERANCE REQUIRED FOR ABSCISSA
C  /RESULT/ PWK2 = MODEL ASSOCIATED WITH X2 AND FB
C  /RESULT/ PWK1 = MODEL ASSOCIATED WITH X1 AND FA
C  /MODEL/ NP = NUMBER OF MODEL PARAMTERS
C
C ON OUTPUT:
C  FB = FUNC(FROOT) - OFF
C  FROOT = ABSCISSA REQUIRED
C  /RESULT/ PWK2 = MODEL ASSOCIATED WITH FB, FROOT
C
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      REAL X1,X2,FA,FB,OFF,TOL
C LOCAL VARIABLES
      REAL AA,B,FC,C,E,DD,TOL1,P,Q,R,S,XM
      INTEGER I,ITER
C LOCAL PARAMETERS
      INTEGER ITMAX
      REAL EPS
      PARAMETER (ITMAX = 100,EPS = 3.E-8)
C FUNCTIONS:
      REAL FUNC
C
C-----------------------------------------------------------------------
      AA = X1
      B = X2
      FA = FA - OFF
      FB = FB - OFF
      IF (FB*FA.GT.0.) WRITE(*,*) 'ROOT NOT BRACKETED IN FROOT'
      FC = FB
C MODELS WILL BE CARRIED AROUND AS:
C   AA: PWK1()
C   B: PWK2()
C   C: PWK3()
C SO COPY PM2 (WITH B) INTO WK (WITH C)
      DO 1 I = 1,NP
1       PWK3(I) = PWK2(I)
      DO 11 ITER = 1,ITMAX
C IF B AND C ARE NOT ON DIFFERENT SIDES OF ROOT THEN MAKE THEM  SO
        IF (FB*FC.GT.0.) THEN
          C = AA
          FC = FA
          DD = B - AA
          E = DD
C COPY PMF (WITH AA) INTO WK (WITH C)
          DO 2 I = 1,NP
2           PWK3(I) = PWK1(I)
        ENDIF
        IF (ABS(FC).LT.ABS(FB)) THEN
C IF B IS NOT CLOSER THAN C MAKE IT SO
          AA = B
          B = C
          C = AA
          FA = FB
          FB = FC
          FC = FA
C SWITCH PM2 (WITH B) AND WK (WITH C) (AA WILL BE LOST IF WE DON'T EXIT)
          DO 3 I = 1,NP
            PWK1(I) = PWK2(I)
            PWK2(I) = PWK3(I)
3           PWK3(I) = PWK1(I)
        ENDIF
        TOL1 = 2.*EPS*ABS(B) + 0.5*TOL
        XM = .5*(C - B)
        IF (ABS(XM).LE.TOL1 .OR. FB.EQ.0.)  THEN
          FROOT = B
          RETURN
        ENDIF
        IF (ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
          S = FB/FA
          IF (AA.EQ.C) THEN
            P = 2.*XM*S
            Q = 1. - S
          ELSE
            Q = FA/FC
            R = FB/FC
            P = S*(2.*XM*Q*(Q - R) - (B - AA)*(R - 1.))
            Q = (Q - 1.)*(R - 1.)*(S - 1.)
          ENDIF
          IF (P.GT.0.) Q =  - Q
          P = ABS(P)
          IF (2.*P .LT. MIN(3.*XM*Q - ABS(TOL1*Q),ABS(E*Q))) THEN
            E = DD
            DD = P/Q
          ELSE
            DD = XM
            E = DD
          ENDIF
        ELSE
          DD = XM
          E = DD
        ENDIF
        AA = B
        FA = FB
C COPY PM2 (WITH B) INTO PMF (WITH AA)
          DO 4 I = 1,NP
4           PWK1(I) = PWK2(I)
        IF (ABS(DD) .GT. TOL1) THEN
          B = B + DD
        ELSE
          B = B + SIGN(TOL1,XM)
        ENDIF
C THE NEW FUNCION EVALUATION AUTOMATICALLY USES PM2 FOR B
        FB = FUNC(B) - OFF
11    CONTINUE
      WRITE(*,*) 'MAXIMUM ITERATIONS EXCEEDED IN FROOT'
      FROOT = B
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      SUBROUTINE CHOLSL(LA,N, A, X, B, ISTAGE)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
C WRITTEN BY:   R.L. PARKER
C
C$$$$$ CALLS NO OTHER ROUTINES
C  ROUTINE FOR BACK-SUBSTITUTION SOLUTION OF A LINEAR SYSTEM THAT HAS
C  ALREADY BEEN FACTORED INTO CHOLESKY FACTORS BY 'CHOLIN' (Q.V.)
C  LA      LEADING DIMENSION OF A
C  N       THE ORDER OF MATRIX
C  A       ARRAY CONTAINING ON AND ABOVE ITS DIAGONAL THE TRANSPOSE OF
C          THE CHOLESKY FACTOR.  THE NORMAL WAY OF FINDING THE FACTORS
C          IS TO CALL  CHOLIN  AS FOLLOWS
C     CALL CHOLIN(N, A, 1, IFBAD)
C     CALL CHOLSL(N, A, X, B, 0)
C  X       THE UNKNOWN VECTOR.  SOLUTION TO EQUATIONS IS RETURNED HERE
C  B       THE RIGHT SIDE OF SYSTEM
C  ISTAGE  THE MODE OF THE SOLUTION.  IF=1 DO NOT COMPLETE BACK-
C     SUBSTITUTION.  INSTEAD RETURN INVERSE(L-TRANS)*B .  OTHER VALUES
C     FORM THE COMPLETE SOLUTION.
C R.L. PARKER
C
C INCLUDES:
      INCLUDE 'imp.inc'
C ARGUMENTS
      INTEGER LA,N,ISTAGE
      REAL A(LA,N+1),X(N),B(N)
C LOCAL VARIABLES
      INTEGER I,K,II,KK,I1
      REAL Z
C
C-----------------------------------------------------------------------
      DO 2500 I=1,N
      Z=B(I)
      IF (I.EQ.1) GO TO 2500
      DO 2200 KK=2,I
      K=I-KK+1
 2200 Z=Z - A(K,I+1)*X(K)
 2500 X(I)=Z/A(I,I+1)
      IF (ISTAGE.EQ.1) RETURN
C  COMPLETE BACK-SUBSTITUTION
      DO 2700 II=1,N
      I=N-II+1
      Z=X(I)
      I1=I+1
      IF (I.EQ.N) GO TO 2700
      DO 2600 K=I1,N
 2600 Z= Z - A(I,K+1)*X(K)
 2700 X(I)=Z/A(I,I+1)
      RETURN
      END                                                              
C
C
C
C-----------------------------------------------------------------------
      SUBROUTINE CHOLIN(LA,N, A, JSTAGE, IFBAD)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
C WRITTEN BY:   R.L. PARKER
C
C$$$$$ CALLS NO OTHER ROUTINES
C  INVERSION AND FACTORIZATION OF POSITIVE DEFINITE SYMMETRIC MATRICES
C  FROM ALGOL ROUTINE OF WILKINSON & REINSCH 'LINEAR ALGEBRA' PP16-17.
C  EXPLANATION
C  CHOLIN  CALCULATES THE FACTOR  L  IN THE CHOLESKY FACTORIZATION
C  A = L L-TRANS, OR THE INVERSE OF L , OR THE INVERSE OF  A.
C  TRANSPOSED STORAGE IS FOR THE CONVENIENCE OF FORTRAN USERS.  HOWEVER,
C  IT MEANS THAT L-TRANS AND INVERSE OF L-TRANS ARE COMPUTED.
C  ARGUMENTS
C    A     AN ARRAY DIMENSION N,N+1.  ON ENTRY THE INPUT MATRIX SHOULD
C          SHOULD OCCUPY THE 1ST  N COLS, OR IF DESIRED, JUST THE LOWER
C          LOWER LEFT TRIANGLE.  ON EXIT THE UPPER TRIANGLE OF THE ARRAY
C          DEFINED BY COLS 2,3...N+1 CONTAINS THE DESIRED RESULT.
C          THIS CAN BE CONVENIENTLY ACCESSED AS A(1,2), OR A(N+1) IF
C          SINGLY DIMENSIONED IN CALLING PROGRAM.
C    N     THE ORDER OF THE MATRIX
C    LA    LEADING DIMENSION OF A
C  JSTAGE  INDICATES DESIRED CALCULATION - IF=1 UPPER TRIANGLE CONTAINS
C          L-TRANS.  IF=2 UPPER TRIANGLE CONTAINS INVERSE OF L-TRANS.
C          IF=0 OR 3 CONTAINS THE UPPER HALF OF THE INVERSE OF A.  NOTE
C          LOWER TRIANGLE IS NOT DESTROYED. IF=-2,-3, THE CHOLESKY
C          FACTORS ARE ASSUMED TO BE IN PLACE ON ENTRY.
C  IFBAD   ERROR FLAG, SET TO 1 IF MATRIX NOT POSITIVE DEFINITE. SET TO
C          0 IF OK.
C R.L. PARKER
C
C INCLUDES:
      INCLUDE 'imp.inc'
C PASS PARAMETERS
      INTEGER LA,N,JSTAGE,IFBAD
      REAL A(LA,N+1)
C LOCAL VARIABLES
      INTEGER I,J,K,KBACK,ISTAGE,I1,J1,J2,N1
      REAL X,Y,Z
C
C-----------------------------------------------------------------------
C  STAGE 1 - FORM FACTORS (HERE L-TRANS IS STOTED ABOVE DIAGONAL)
      ISTAGE=IABS(JSTAGE)
      IF (JSTAGE.LT.0) GO TO 1550
      IFBAD=0
      DO 1500 I=1,N
      I1=I + 1
      DO 1400 J=I,N
      J1=J + 1
      X=A(J,I)
      IF (I.EQ.1) GO TO 1200
      DO 1100 KBACK=2,I
      K=I-KBACK+1
 1100 X=X - A(K,J1)*A(K,I1)
 1200 IF (I.NE.J) GO TO 1300
      IF (X.LE.0.0) GO TO 4000
      Y=1.0/SQRT(X)
      A(I,I1)=Y
      GO TO 1400
 1300 A(I,J1)=X*Y
 1400 CONTINUE
 1500 CONTINUE
      IF (ISTAGE.NE.1) GO TO 2000
 1550 DO 1600 I=1,N
 1600 A(I,I+1)=1.0/A(I,I+1)
      IF (ISTAGE.EQ.1) RETURN
C  STAGE 2 - INVERSION OF L-TRANS
 2000 DO 2600 I=1,N
      I1=I+1
      IF (I1.GT.N) GO TO 2600
      DO 2500 J=I1,N
      Z=0.0
      J1=J + 1
      J2=J - 1
      DO 2200 KBACK=I,J2
      K=J2-KBACK+I
 2200 Z=Z - A(K,J1)*A(I,K+1)
      A(I,J1)=Z*A(J,J1)
 2500 CONTINUE
 2600 CONTINUE
      IF (ISTAGE.EQ.2) RETURN
C  STAGE 3  - CONSTRUCTION OF IVERSE OF  A  ABOVE DIAGONAL.
      DO 3500 I=1,N
      DO 3500 J=I,N
      Z=0.0
      N1=N + 1
      J1=J + 1
      DO 3200 K=J1,N1
 3200 Z=Z + A(J,K)*A(I,K)
      A(I,J+1)=Z
 3500 CONTINUE
      RETURN
C  ERROR EXIT
 4000 IFBAD=1
C ERROR MESSAGE SUPRESSED - S.CONSTABLE
C      WRITE(*,400)
C 400  FORMAT(44H0CHOLIN FAILS - MATRIX NOT POSITIVE DEFINITE           )
      RETURN
      END  
C
C
C-----------------------------------------------------------------------
      SUBROUTINE MAKJTJ(TOL0)
C-----------------------------------------------------------------------
C 
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
C CONSTRUCTS MATRICES THAT DEPEND ON THE JACOBIAN (WJTWJ AND WJTWD).
C THE JACOBIAN ONLY NEEDS TO BE AROUND LONG ENOUGH TO DO
C THIS, SO BY STUFFING THIS INTO A SUBROUTINE DYNAMIC MEMORY ALLOCATION
C WILL MAKE WJ DISSAPEAR AFTER USE.
C 
C ON INPUT:
C   CURRENT MODEL AND DATA STORED IN COMMON /MODEL/ AND /DATA/ 
C ON OUTPUT:
C   /RESULT/ WJTWJ = W.J.TRANS(W.J)
C   /RESULT/ WJTWD = WEIGHTED, TRANSLATED DATA PREMULTIPLIED BY TRANS(W.J)
C   TOL0 = RMS MISFIT ASSOCIATED WITH CURRENT MODEL
C CALLS:
C   ATAMUL, MULT, TRMULT, FORDIV, ANORM
C
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      REAL TOL0
C LOCAL VARIABLES
      REAL WJ(NDD,NPP), DWK1(NDD), DHAT(NDD)
      INTEGER I, J
C FUNCTIONS 
      REAL ANORM
C
C-----------------------------------------------------------------------
C CALCULATE THE MATRIX OF PARTIALS AND MODEL RESPONSE
      CALL FORDIV(NP,ND,PM,DP,DM,WJ)
C CALC MISFIT VECTOR AND MISFIT
      DO 60 I = 1,ND
60      DWK1(I) = (D(I) - DM(I))/SD(I)
      TOL0 = SQRT(ANORM(ND,DWK1)/ND)
C WEIGHT THE JACOBIAN MATRIX
      DO 70 I = 1,ND
        DO 70 J = 1,NP
70        WJ(I,J) = WJ(I,J)/SD(I)
C FORM W.J.TRANS(W.J)
      CALL ATAMUL(NDD,ND,NPP,NP,WJ,WJTWJ)
C FORM THE WEIGHTED, TRANSLATED DATA AND PREMULTIPLY BY TRANS(W.J)
      CALL MULT(NDD,ND,NPP,NP,1,WJ,PM,DHAT)
      DO 80 I = 1,ND
80      DHAT(I) = DWK1(I) + DHAT(I)
      CALL TRMULT(NDD,ND,NPP,NP,1,WJ,DHAT,WJTWD)
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      SUBROUTINE MAKPTP(IRUF)      
C-----------------------------------------------------------------------
C 
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
C CONSTRUCTS PENALTY MATRIX (DEL) MULTIPLIED BY ITS TRANSPOSE.
C THE PENALTY MATRIX ONLY NEEDS TO BE AROUND LONG ENOUGH TO DO
C THIS, SO BY STUFFING THIS INTO A SUBROUTINE DYNAMIC MEMORY ALLOCATION
C CAN GET RID OF IT AFTERWARDS
C 
C ON INPUT:
C   IRUF = 0 FOR NO ROUGHNENING (I.E. MIN. NORM MODELS) (ADDED IN V1.5)
C        = 1 FOR 1ST DERIVATIVE ROUGHNESS
C          2 FOR 2ND DERIVATIVE ROUGHNESS
C   /MODEL/ NP = NUMBER OF MODEL PARAMETERS
C ON OUTPUT:
C   /RESULT/ PTP() = DEL(TRANS)*DEL
C CALLS:
C   OBJMAT, TRMULT
C
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      INTEGER IRUF
C LOCAL VARIABLES
      INTEGER NPEN, I
      REAL DEL(NNP,NPP)
C
C-----------------------------------------------------------------------
C MAKE PARTIALS MATRIX:
      CALL OBJMAT(IRUF, DEL, NPEN)
      IF (abs(idebug) .EQ. 1) 
     *   WRITE(*,*) ' NUMBER OF PENALTY TERMS =', NPEN
C FORM TRANS(DEL).DEL
      CALL TRMULT(NNP,NPEN,NPP,NP,NP,DEL,DEL,PTP)
C ADD THE WEIGHTS FOR THE MODEL PREJUDICE
      DO 10 I = 1,NP
10      PTP(I,I) = PTP(I,I) + PREWTS(I)*PREWTS(I)
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      REAL FUNCTION FNDRUF(IRUF)      
C-----------------------------------------------------------------------
C 
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
C COMPUTES THE ROUGHNESS OF THE MODEL. FOR MEMORY ALLOCATION CONSIDERATIONS,
C THE PENALTY MATRIX HAS NOT BEEN KEPT AROUND, SO WE NEED TO REFORM IT.
C 
C ON INPUT:
C   IRUF = 0 FOR NO ROUGHNENING (I.E. MIN. NORM MODELS) (ADDED IN V1.5)
C        = 1 FOR 1ST DERIVATIVE ROUGHNESS
C          2 FOR 2ND DERIVATIVE ROUGHNESS
C   /RESULT/ PWK2() = LOCATION OF THE MODEL WE WANT TO MEASURE
C   /MODEL/ NP = NUMBER OF MODEL PARAMETERS
C ON OUTPUT:
C   FNDRUF = ROUGHNESS OF MODEL
C CALLS:
C   OBJMAT, TRMULT
C
C INCLUDES:
      INCLUDE 'imp.inc'
      INCLUDE 'occamdim.inc'
      INCLUDE 'occam.inc'
C
C ARGUMENTS:
      INTEGER IRUF
C LOCAL VARIABLES
      INTEGER I, K, NPEN
      REAL DEL(NNP,NPP), SUM, VECT(NNP)
C FUNCTIONS
      REAL ANORM
C
C-----------------------------------------------------------------------
C MAKE PARTIALS MATRIX:
      CALL OBJMAT(IRUF, DEL, NPEN)
      DO 170 I = 1,NPEN
        SUM = 0.0
        DO 165 K = 1,NP
165       SUM = SUM + DEL(I,K)*PWK2(K)
170     VECT(I) = SUM
      FNDRUF = ANORM(NPEN,VECT)
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      SUBROUTINE ATAMUL(MAD,MA,NAD,NA,A,C)
C-----------------------------------------------------------------------
      INCLUDE 'imp.inc'
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
C WRITTEN BY:  C. DEGROOT-HEDLIN   
C
C DESCRIPTION: 
C
C MULTIPLIES THE TRANSPOSE OF A 2D MATRIX BY ITSELF
C TAKES ADVANTAGE OF SYMMETRY ABOUT DIAGONAL
C
C
C ARGUMENTS
      INTEGER MAD,MA,NAD,NA
      REAL A(MAD,NA),C(NAD,NA)
C
C LOCAL VARIABLES
      INTEGER I,J,K
      REAL CIJ
C
C-----------------------------------------------------------------------
      DO 12 I = 1,NA
        DO 12 J = I,NA
          CIJ = 0.0
          DO 10 K = 1,MA
10          CIJ = A(K,I)*A(K,J) + CIJ
12        C(I,J) = CIJ
      DO 15 I=2,NA
        DO 15 J=1,I-1
15        C(I,J) = C(J,I)
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      subroutine itrout (itform,descr,modfil,datfil,maxitr,tolreq,
     *    iruf,nit,pmu,rlast,tobt,ifftol,iof2)      
C-----------------------------------------------------------------------
c
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
c writes an iteration file for OCCAM
c
c includes:
      include 'imp.inc'
      include 'occamdim.inc'
      include 'occam.inc'
c   uses np, pm(), idebug
c arguments:
      integer maxitr, iruf, nit, ifftol,iof2
      real tolreq, pmu, rlast, tobt
      character*80 modfil, datfil, itform, descr, datetm
c local variables
      integer lerr, i
c   lerr = error flag on file operation
c   i = loop index
      character*2 num
c
C-----------------------------------------------------------------------
c create file name suffix
      write (num, '(i2)') nit
      if (nit .le. 9) num(1:1) = '0'
c open iteration file
      open (unit=iof2, file='ITER'//num, iostat=lerr)
      if (lerr .ne. 0) then
c we could be in the middle of an iteration, so don't die on error, but
c send a message and forge on regardless
        write(*,*) ' Error opening iteration file'
        return
      end if
      call datime(datetm)
c
c write stuff      
      write(unit=iof2,fmt='(a)') 'FORMAT:           '//itform
      write(unit=iof2,fmt='(a)') 'DESCRIPTION:      '//descr
      write(unit=iof2,fmt='(a)') 'MODEL FILE:       '//modfil
      write(unit=iof2,fmt='(a)') 'DATA FILE:        '//datfil
      write(unit=iof2,fmt='(a)') 'DATE/TIME:        '//datetm
      write(unit=iof2,fmt='(a,i3)')    'MAX ITER:         ',maxitr
      write(unit=iof2,fmt='(a,f6.3)')  'REQ TOL:          ',tolreq
      write(unit=iof2,fmt='(a,i3)')    'IRUF:             ',iruf
      write(unit=iof2,fmt='(a,i3)')    'DEBUG LEVEL:      ',idebug
      write(unit=iof2,fmt='(a,i3)')    'ITERATION:        ',nit
      write(unit=iof2,fmt='(a,g15.7)') 'PMU:              ',pmu
      write(unit=iof2,fmt='(a,g15.7)') 'RLAST:            ',rlast
      write(unit=iof2,fmt='(a,g15.7)') 'TLAST:            ',tobt
      write(unit=iof2,fmt='(a,i3)')    'IFFTOL:           ',ifftol
      write(unit=iof2,fmt='(a,i4)')    'NO. PARMS:        ',np
      write(unit=iof2,fmt='(4g15.7)') (pm(i), i=1,np)
C
      close(unit=iof2)
C      
      return
      end
c
C
C-----------------------------------------------------------------------
      subroutine itrin (itform,descr,modfil,datfil,maxitr,tolreq,
     *    iruf,nit,pmu,rlast,tobt,ifftol,iof2)      
C-----------------------------------------------------------------------
c
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.01, 13 Jan 1993
C
c reads a STARTUP or ITER file for OCCAM
c
c includes:
      include 'imp.inc'
      include 'occamdim.inc'
c   uses npp
      include 'occam.inc'
c   uses np, pm(), idebug
c arguments:
      integer maxitr, iruf, nit, ifftol,iof2
      real tolreq, pmu, rlast, tobt
      character*80 modfil, datfil, itform, descr
c local variables
      integer i, idcode
c   i = loop index
c   idcode = function to perform the equivalent of a free-format internal
c            read for an integer
      real rdcode
c   rdcode = function to perform the equivalent of a free-format internal
c            read for a real
      character*80 string
C
C-----------------------------------------------------------------------
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) itform
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) descr
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) modfil
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) datfil
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
c
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      maxitr = idcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      tolreq = rdcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      iruf = idcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      idebug = idcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      nit = idcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      pmu = rdcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      rlast = rdcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      tobt = rdcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      ifftol = idcode(string,iof2)
      read(unit=iof2,fmt='(18x,a)', end=198, err=199) string
      np = idcode(string,iof2)
      if (np .gt. npp) call filerr (' Too many parameters', iof2)
      read (unit=iof2,fmt= *, end=198, err=199) (pm(i), i=1,np)
c
      return
c
198      call filerr (' Startup file ended prematurely', iof2)
199      call filerr (' Error reading startup file', iof2)
      end
c
C
C-----------------------------------------------------------------------
      real function rdcode(string, iof2)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 Jan 1993
C
c rdcode provides a patch to get around the lack of an internal
c free-format (list-directed) read on some compilers, notably the
c cray.  This version for reals.  Butchered from decode by Bob Parker.
c 
c on input:
c    string = character string containing number to be read
c    iof2 = unit number of open file, to be closed on error 
c on output:
c    rdcode contains a real number from string
c calls:
c    filerr
c
c includes:
      include 'imp.inc'
c input arguments:
      character*(*) string
	integer iof2
c local variables
      character*40 local
	integer k, k1, kn, l
c
c  Terminate line at % sign
      kn=index(string//'%', '%') - 1
      k1=1
      do 1100 k=k1, kn
        if (string(k:k) .ne. ' ') goto 1200
 1100 continue
      call filerr (' No number found by rdcode', iof2)
 1200 do 1300 l=k, kn
        if (string(l:l).eq. ',' .or. string(l:l) .eq. ' ') goto 1500
 1300 continue
 1500 local=string(k:l-1)
      read (local, '(f40.0)', err=1900) rdcode
	return
 1900 call filerr (' Error reading string in rdcode', iof2)
      end
C-----------------------------------------------------------------------
      integer function idcode(string, iof2)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 Jan 1993
C
c idcode provides a patch to get around the lack of an internal
c free-format (list-directed) read on some compilers, notably the
c cray.  This version for integers.  Butchered from decode by Bob Parker.
c 
c on input:
c    string = character string containing number to be read
c    iof2 = unit number of open file, to be closed on error 
c on output:
c    idcode contains an integer number from string
c calls:
c    filerr
c
c includes:
      include 'imp.inc'
c input arguments:
      character*(*) string
	integer iof2
c local variables
      character*40 local
	integer k, k1, kn, l
	real realno
c
c  Terminate line at % sign
      kn=index(string//'%', '%') - 1
      k1=1
      do 1100 k=k1, kn
        if (string(k:k) .ne. ' ') goto 1200
 1100 continue
      call filerr (' No number found by idcode', iof2)
 1200 do 1300 l=k, kn
        if (string(l:l).eq. ',' .or. string(l:l) .eq. ' ') goto 1500
 1300 continue
 1500 local=string(k:l-1)
      read (local, '(f40.0)', err=1900) realno
	idcode = int(realno)
	return
 1900 call filerr (' Error reading string in idcode', iof2)
      end
C-----------------------------------------------------------------------
      subroutine filerr(mssg, io1)
C-----------------------------------------------------------------------
C
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 2.00, 13 May 1992
C
c filerr prints an error message, closes file io1, and stops
c
c on input:
c    mssg = character string containing error message
c    i01 = unit number of file to close (0 = no file open)
c on output:
c    outputs nothing
c calls:
c    no other routines
c
c
c includes:
      include 'imp.inc'
c input arguments:
      character*(*) mssg
      integer io1
c
C-----------------------------------------------------------------------
      write(*,*) mssg
      if (io1 .gt. 0) close (io1)
      stop
      end
C-----------------------------------------------------------------------
      subroutine datime(datetm)
c-----------------------------------------------------------------------
c OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 1.00, 18 May 1992
c
c Dummy date and time utility for OCCAM 2.0.  
c
c on input:
c   nothing
c on output:
c   datetm = 80 character string containing message saying that date and time
c           are not available
c      
      character*80 datetm
c      
c-----------------------------------------------------------------------
      datetm = 'No date and time available on this system'      
      return
      end
c-----------------------------------------------------------------------
      
