C----------------------------------------------------------------------- SUBROUTINE FORMOD(NP,ND,PM,DP,DM) C----------------------------------------------------------------------- c OCCAM1DMT: Steven Constable 10 March 1993 C FORMOD IS AN M.T. FORWARD MODELLING MODULE FOR USE BY OCCAM. C COMPUTES MT RESPONSES IN LOG DOMAIN FOR EITHER C RESPONSES C OR APPARENT RESISTIVITY AND PHASE. C C CALLS: MTLAY C C ON INPUT: C NP = NUMBER OF PARAMETERS IN THE PARAMETER ARRAY PM() C ND = NUMBER OF DATA ELEMENTS (I.E. TWICE THE NUMBER OF FREQUENCIES) C PM() = ARRAY OF PARAMETER ELEMENTS, LOG10(LAYER RESISTIVITIES) C STARTING AT THE SURFACE LAYER C DP() = ARRAY ALLOWING UP TO FOUR PARAMETERS ASSOCIATED WITH EACH c DATUM. Two are used for 1D MT: C data type in DP(I,1): C 1 = LOG10 RESISTIVITY C 2 = PHASE C 3 = LINEAR REAL C C 4 = LINEAR -(IMAG. C) C 5 = LOG10 REAL C C 6 = LOG10 -(IMAG. C) C 7 = LINEAR RESISTIVITY C and PERIOD IN SECONDS STORED IN DP(I,2), I=1,..ND. C /FMODEL/ in mt1D.inc IS USED TO CARRY AROUND THE ACTUAL RESISTIVITIES C AND THICKNESSES OF A LAYERED MODEL, SINCE THE PARAMETERISATION IN C PM() IS IN LOG10(RESISTIVITIES) WITH NO THICKNESSES BEING CHANGED C IN THE INVERSION. IT IS THE CALLING PROGRAM'S RESPONSIBILITY TO C SET: C NLAYER : NUMBER OF LAYERS IN THE MODEL C THICK(I), I=1,NLAYER-1 : LAYER THICKNESSES IN METRES c iff(NP), a lookup table for the resistivities which are allowed to c vary. C C ON OUTPUT: C DM() = ARRAY OF RESPONSE DATA ELEMENTS: C include 'occamdim.inc' include 'mt1D.inc' COMPLEX C REAL PM(NPP),DP(NDD,4),DM(NDD) C RMAX IS USED TO RESTRICT THE MAXIMUM RESISTIVITY TO 10**+-20 PARAMETER (RMAX = 20., RMIN = -20.) C C MODEL IS PARAMETERISED AS LOGS IN PM(). COPY ACTUAL C RESISTIVITIES INTO RHO() FOR FORWARD COMPUTATION. DO 5 I = 1,np 5 RHO(iff(i)) = 10.**MIN(RMAX,MAX(RMIN,PM(I))) C DO 100 J = 1, ND ITYPE = INT(DP(J,2) + 0.1) FREQ = 1./DP(J,1) CALL MTLAY(FREQ,C,NLAYER,RHO,THICK) GOTO (10, 20, 30, 40, 50, 60, 70) ITYPE 10 DM(J) = LOG10(7.8957E-06*(ABS(C)**2)/DP(J,1)) GOTO 100 20 DM(J) = 57.29578*ATAN2(REAL(C),-AIMAG(C)) GOTO 100 30 DM(J) = REAL(C) GOTO 100 40 DM(J) = -AIMAG(C) GOTO 100 50 DM(J) = LOG10(REAL(C)) GOTO 100 60 DM(J) = LOG10(-AIMAG(C)) GOTO 100 70 DM(J) = 7.8957E-06*(ABS(C)**2)/DP(J,1) 100 CONTINUE RETURN END C C C----------------------------------------------------------------------- SUBROUTINE FORDIV(NP,ND,PM,DP,DM,G) C----------------------------------------------------------------------- c OCCAM1DMT: Steven Constable 10 March 1993 C FORDIV IS AN M.T. FORWARD MODELLING MODULE FOR USE BY OCCAM AND OTHER C SUBROUTINES. RETURNS THE RESPONSE OF A MODEL AS WELL AS THE PARTIALS C MATRIX. COMPUTES MT RESPONSES IN LOG DOMAIN FOR EITHER C RESPONSES C OR APPARENT RESISTIVITY AND PHASE. C C SEE FORMOD FOR PARAMETER LISTS, EXCEPT ADD: C CALLS: MTDIV C ON OUTPUT: C G(ND,NP) = PARTIALS MATRIX: C G(I,J) = D(DATUM(I))/D(LOG10(RESISTIVITY(J))) C include 'occamdim.inc' include 'mt1D.inc' REAL PM(NPP),DP(NDD,4),DM(NDD),G(NDD,NPP) COMPLEX C,DCDP(NPP) C RMAX IS USED TO RESTRICT THE MAXIMUM RESISTIVITY TO 10**10 PARAMETER (RMAX = 20., RMIN = -20.) C C MODEL IS PARAMETERISED AS LOGS IN PM(). COPY ACTUAL C RESISTIVITIES INTO RHO() FOR FORWARD COMPUTATION. DO 5 I = 1,np 5 RHO(iff(i)) = 10.**MIN(RMAX,MAX(RMIN,PM(I))) C DO 200 J = 1, ND ITYPE = INT(DP(J,2) + 0.1) FREQ = 1./DP(J,1) CALL MTDIV(FREQ,C,DCDP,NLAYER,RHO,THICK) AC2 = ABS(C)**2 REALC = REAL(C) AIMGC = AIMAG(C) GOTO (10, 20, 30, 40, 50, 60, 70) ITYPE 10 DM(J) = LOG10(7.8957E-06*(ABS(C)**2)*FREQ) GOTO 100 20 DM(J) = 57.29578*ATAN2(REAL(C),-AIMAG(C)) GOTO 100 30 DM(J) = REALC GOTO 100 40 DM(J) = -AIMGC GOTO 100 50 DM(J) = LOG10(REALC) GOTO 100 60 DM(J) = LOG10(-AIMGC) GOTO 100 70 DM(J) = 7.8957E-06*(ABS(C)**2)*FREQ 100 CONTINUE DO 190 I = 1, NP REALD = REAL(DCDP(I)) AIMGD = AIMAG(DCDP(I)) GOTO (110, 120, 130, 140, 150, 160, 170) ITYPE 110 G(J,I) = 2.*RHO(I)*(REALC*REALD + AIMGC*AIMGD)/AC2 GOTO 190 120 G(J,I) = 131.9284*RHO(I)*(REALC*AIMGD - AIMGC*REALD)/AC2 GOTO 190 130 G(J,I) = RHO(I)*REALD/0.43429 GOTO 190 140 G(J,I) = -RHO(I)*AIMGD/0.43429 GOTO 190 150 G(J,I) = RHO(I)*REALD/REALC GOTO 190 160 G(J,I) = RHO(I)*AIMGD/AIMGC GOTO 190 170 G(J,I) = 3.63610E-05*RHO(I)*(REALC*REALD+AIMGC*AIMGD)*FREQ 190 CONTINUE 200 CONTINUE RETURN END C C C----------------------------------------------------------------------- SUBROUTINE MTLAY(F,Carg,NLAYER,RHO,THICK) C----------------------------------------------------------------------- c OCCAM1DMT: Steven Constable 10 April 1992 C C MTLAY COMPUTES THE COMPLEX C VALUE OVER A LAYERED EARTH C S.C. CONSTABLE, S.I.O. LA JOLLA CA92093, AUGUST 1986 C C CALLS: NO OTHER ROUTINES C ON INPUT: C F = FREQUENCY IN HERTZ C NLAYER = NUMBER OF LAYERS IN MODEL, INCLUDING TERMINATING HALF SPACE C RHO = ARRAY OF LAYER RESISTIVITIES, OHM-M (LINEAR DOMAIN) C THICK = ARRAY OF LAYER THICKNESSES, M (LINEAR DOMAIN) C ON OUTPUT: C C = COMPLEX C VALUE, M (LINEAR DOMAIN) C DIMENSION THICK(NLAYER-1),RHO(NLAYER) complex carg c Local variables are double precision. This is necessary. COMPLEX*16 CK,OMI,C OMEGA = 6.283185308*F OMI = dCMPLX(OMEGA)*1.256637062E-06*(0.,1.) C = 1./SQRT(OMI/dCMPLX(RHO(NLAYER))) DO 100 I = NLAYER-1,1,-1 CK = SQRT(OMI/dcmplx(RHO(I))) C = EXP(-2.*CK*dcmplx(THICK(I)))*(CK*C-1.)/(CK*C+1.) C = (1.+C)/(1.-C)/CK 100 CONTINUE carg = c RETURN END C C C C----------------------------------------------------------------------- SUBROUTINE MTDIV(FREQ,Carg,Darg,NLAYER,RHO,THICK) C----------------------------------------------------------------------- c OCCAM1DMT: Steven Constable 10 April 1992 C C MTDIV COMPUTES THE COMPLEX C VALUE OVER A LAYERED EARTH AND DC/DRHO FOR C EACH LAYER (I.E. A ROW OF THE JACOBIAN MATRIX) C S.C. CONSTABLE, S.I.O. LA JOLLA CA92093, AUGUST 1986 C C CALLS: NO OTHER ROUTINES C ON INPUT: C FREQ = FREQUENCY IN HERTZ C NLAYER = NUMBER OF LAYERS IN MODEL, INCLUDING TERMINATING HALF SPACE C RHO = ARRAY OF LAYER RESISTIVITIES, OHM-M (LINEAR DOMAIN) C THICK = ARRAY OF LAYER THICKNESSES, M (LINEAR DOMAIN) C ON OUTPUT: C C = COMPLEX C VALUE, M (LINEAR DOMAIN) C D = ARRAY OF COMPLEX DERIVATIVES OF C W.R.T. LAYER RESISTIVITIES, C M/OHM-M C include 'occamdim.inc' DIMENSION THICK(NLAYER-1),RHO(NLAYER) COMPLEX Darg(NLAYER),Carg c Local variables are double precision. This is necessary. COMPLEX*16 D(npp),C COMPLEX*16 OMI,CKM1,CKP1,OPEXP,OMEXP,CK,DCIIP1,THEEXP,RI C OMEGA = 6.283185308*FREQ OMI = dcmplx(OMEGA)*1.256637062E-06*(0.,1.) C FIRST OF ALL DC/DRHO FOR I=1,2,..NLAYER WILL BE STORED IN D() C THEN PROPAGATED TO THE SURFACE USING DC(I)/DC(I-1) AT EACH LAYER. C C ESTABLISH THE VALUES FOR THE TOP OF THE TERMINATING HALF-SPACE CK = SQRT(OMI/dcmplx(RHO(NLAYER))) C = 1./CK D(NLAYER) = 1./(CK*dcmplx(RHO(NLAYER))*2.) C NOW SWEEP UP THROUGH THE LAYERS DO 100 I = NLAYER-1,1,-1 CK = SQRT(OMI/dcmplx(RHO(I))) CKP1 = C*CK + 1. CKM1 = C*CK - 1. RI = CKM1/CKP1 THEEXP = EXP(-2.*CK*dcmplx(THICK(I))) OPEXP = THEEXP*RI OMEXP = 1. - OPEXP OPEXP = 1. + OPEXP C COMPUTE DC/DRHO FOR THE (I)TH LAYER D(I) = 2.*THEEXP*(C/CKP1**2 - dcmplx(THICK(I))*RI)/OMEXP**2 D(I) = OPEXP/OMEXP/CK/2. - D(I) D(I) = D(I)/dcmplx(RHO(I)) C COMPUTE DC(I)/DC(I+1) FOR THE (I)TH LAYER DCIIP1 = 4.*THEEXP/(CKP1*OMEXP)**2 C NOW COMPUTE THE VALUE OF C AT THE TOP OF THIS LAYER C = OPEXP/OMEXP/CK C MULTIPLY THE PRODUCT SUMS FOR ALL THE LAYERS BELOW THIS ONE, WHICH HAVE C ALREADY BEEN PARTIALLY ACCUMULATED. DO 20 J = I+1,NLAYER 20 D(J) = D(J)*DCIIP1 100 CONTINUE c copy double complex results into single complex arguments. carg = c do 200 i=1,nlayer 200 darg(i) = d(i) C AND THAT'S IT RETURN END C C C----------------------------------------------------------------------- SUBROUTINE OBJMAT(IRUF, DEL, NPEN) C----------------------------------------------------------------------- c OCCAM1DMT: Steven Constable 10 March 1993 C CONSTRUCTS PENALTY MATRIX (DEL). C C THIS VERSION IS FOR 1D CASE C C ON INPUT: C IRUF = 0 FOR NO ROUGHNENING (I.E. MIN. NORM MODELS) C = 1 FOR 1ST DERIVATIVE ROUGHNESS C 2 FOR 2ND DERIVATIVE ROUGHNESS C /MODEL/ NP = NUMBER OF MODEL PARAMETERS C ON OUTPUT: C DEL() = PENALTY MATRIX C NPEN = NUMBER OF ROWS IN PENALTY MATRIX C CALLS: C MULT, INITM C INCLUDE 'imp.inc' INCLUDE 'occamdim.inc' INCLUDE 'occam.inc' include 'mt1D.inc' C C ARGUMENTS: REAL DEL(NNP,NPP) INTEGER IRUF, NPEN C LOCAL VARIABLES INTEGER I, J REAL DEL2(NNP,NPP) C C FOR THE 1D CASE THE ROW COUNT OF PENALTY MATRIX IS JUST NP NPEN = NP CALL INITM(NNP,NPEN,NP,DEL) C MIN. NORM MODELS JUST PICK OFF THE SIZE OF THE MODEL PARAMETERS IF (IRUF .EQ. 0) THEN DO 5 I = 1,NP 5 DEL(I,I) = 1.*FLOAT(ipen1(iff(i))) ELSE C MIN. FIRST DERIVATIVE TAKE DIFFERENCES OF ADJACENT LAYERS DO 10 I = 1,NP-1 DEL(I,I) = -1.*FLOAT(ipen1(iff(i))) 10 DEL(I,I+1) = 1.*FLOAT(ipen1(iff(i))) END IF C SQUARE MATRIX AND COPY BACK INTO DEL IF SECOND DERIV. SMOOTHING REQUIRED IF (IRUF .EQ. 2) THEN CALL MULT(NNP,NPEN,NNP,NPEN,NPEN,DEL,DEL,DEL2) DO 30 I = 1,NPEN DO 30 J = 1,NPEN 30 DEL(J,I) = DEL2(J,I) END IF C SET THE WEIGHTS FOR THE PREJUDICAL MODEL DO 40 I = 1,NP 40 PREWTS(I) = FLOAT(IPEN2(iff(i))) RETURN END C----------------------------------------------------------------------- SUBROUTINE INPUTD(datfil) C----------------------------------------------------------------------- c OCCAM1DMT: Steven Constable 10 March 1993 C C INPUTD READS IN A STANDARD DATA SET FROM FILE. C INPUTD DOES NO SCALING OR TRANSFORMATION (I.E. IF LOG(DATA) C ARE REQUIRED, LOG(DATA) MUST BE SUPPLIED). C EACH LINE IS ASSUMED TO CONTAIN 1 PERIOD, A DATUM, AND DATA C ERROR (ABSOLUTE) IN THAT ORDER. THE READ IS FREE FORMAT. C C ON OUTPUT: C data in occam's common block C include 'occamdim.inc' include 'occam.inc' c input arguments: character*(*) datfil c datfil = data file name c local variables character*80 dform, ddescr, string c dform = data format c ddescr = data description c string = string buffer for free format reads c calls: c filerr C if (idebug .ge. 1) write(*,*) ' reading data file...' c open data file open (15, file=datfil, status='old', iostat=lerr) if (lerr .ne. 0) * call filerr (' Error opening data file', 0) ndp = 2 C c read character descriptors read (15,'(18x,a)', end=198, err=199) dform if (dform(1:15) .ne. 'OCCAM1MTDAT_2.0') * call filerr (' Wrong data file type', 15) read (15,'(18x,a)', end=198, err=199) ddescr c read in the number of layers read (15,'(18x,a)', end=198, err=199) string read (string,*, end=198, err=199) nd if (nd .gt. ndd) call filerr (' Max no. of data exceeded', 15) read (15,'(a)', end=198, err=199) string do 20 i = 1, nd 20 read (15,*,end=198,err=199) dp(i,1), dp(i,2), d(i), sd(i) close(15) if (idebug .ge. 1) write(*,*) nd, ' data read' return C 198 call filerr (' Data file ended prematurely', 15) 199 call filerr (' Error reading data file', 15) end C----------------------------------------------------------------------- subroutine inputm(modfil) C----------------------------------------------------------------------- c OCCAM1DMT: Steven Constable 10 March 1993 c inputm opens the model file modfil c and reads in the model description. c c arguments: character*(*) modfil c modfil = model file name c includes: include 'occamdim.inc' include 'occam.inc include 'mt1D.inc' c uses np c local variables: character*80 mform, mtitle, mdescr, string c mform = model format c mtitle = model title c mdescr = model description c string = string buffer for free format reads real ppmod(npp) c ppmod is a tempory resting place for the model prejudice c calls: c filerr c if (idebug .ge. 1) write(*,*) ' reading model file...' c open model file open (15, file=modfil, status='old', iostat=lerr) if (lerr .ne. 0) * call filerr (' Error opening model file', 0) c read character descriptors read (15,'(18x,a)', end=198, err=199) mform if (mform(1:15) .ne. 'OCCAM1MTMOD_2.0') * call filerr (' Wrong model file type', 15) read (15,'(18x,a)', end=198, err=199) mtitle read (15,'(18x,a)', end=198, err=199) mdescr c read in the number of layers read (15,'(18x,a)', end=198, err=199) string read (string,*, end=198, err=199) nlayer if (nlayer .gt. npp) * call filerr (' Max no. of model params exceeded', 15) read (15,'(a)', end=198, err=199) string c read in the thicknesses etc do 10 i = 1, nlayer 10 read (15, *, end=198, err=199) thick(i),rho(i), ipen1(i), * ppmod(i),ipen2(i) c we are done with the model file close (15) c c now sort out fixed versus free resistivities: nfree = 0 do 20 i = 1, nlayer if (rho(i) .le. 0.0) then nfree = nfree + 1 iff(nfree) = i end if 20 continue if (nfree .ne. np) * call filerr (' Inconsistent number of free parameters', 0) c set the prejudical model DO 40 I = 1,NP 40 premod(I) = ppmod(iff(i)) c c return c 198 call filerr (' Model file ended prematurely', 15) 199 call filerr (' Error reading model file', 15) end