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
	
