【文件属性】:
文件名称:occam一维反演
文件大小:55KB
文件格式:F
更新时间:2016-05-11 03:34:42
fortran代码 occam
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-----------------------------------------------------------------------
网友评论
- 下载了,然后自己重新写了,有参考价值
- 关于雅克比矩阵的部分在哪里啊
- 这个程序真的很不错哦,我真的受益非浅