c ****** These routines have been hackup up to do the forward 
c ****** computation only and save the large amount of memory 
c ****** needed to compute derivatives.  FORDIV has been removed 
c ****** to avoid the jacobian matrix, and PW2DI has been hacked 
c ****** to remove RJ, which is huge.
c
c-----------------------------------------------------------------------
	subroutine formod(np, nd, pm, dp, dm)
c-----------------------------------------------------------------------
c
c OCCAM2DMT 1.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 1.01, 13 Jan 1993
c
c includes:
      include 'imp.inc'
      include 'mt2D.inc'
c uses nfre, frek(), nrl, nrfix, y(), nbrick
c uses rhte(), phte(), kzy(), rhtm(), phtm()
c
c arguments:
      real pm(npp), dp(ndd,4), dm(ndd)
      integer np, nd
c local variables
      integer i, j, nexec, nres, itype, isite, ifreq
c   i,j = loop counters
c   nexec = instruction to PW2DI to perform only forward calculation
c   nres = number of resistivities in y() array
c   itype, isite, ifreq = integer versions of dp() entries
       real tesum, tmsum
c   tesum = sum of logs of TE free static shifts 
c   tmsum = sum of logs of TM free static shifts 
c 
c-----------------------------------------------------------------------
c total number of resistivities = fixed resistivities + free res'ties 
c (but not air)
      nres = nrfix + nbrick
      nexec = 1
      tesum = 0.0
      tmsum = 0.0
c      
c copy across resistivities from parameter array to y()
      do 10 i = 1, nbrick
10        y(nrfix+1+i) = 10.**pm(i)
c copy across any static shifts in the free parameters
      if (nfss .gt. 0) then
        do 12 i = 1,nsb
          if (isfx(i) .ne. 0) shift(i) = pm(isfx(i))
12        continue
      end if
c
c compute the MT response
      call PW2DI(nres,nexec)
c            
c strip off data:
c  the first column of dp() is the site number:
c  the second column of dp() is the frequency number:
c  the third column of dp() is the data type:
c   1 = TE resistivity
c   2 = TE phase
c   3 = real(tipper)
c   4 = imag(tipper)
c   5 = TM resistivity
c   6 = TM phase
c
c for every datum:
      do 150 i = 1, ndata
        isite = int(dp(i,1)+0.1)
        ifreq = int(dp(i,2)+0.1)
        itype = int(dp(i,3)+0.1)
        goto (50, 60, 70, 80, 90, 100) itype
        
50        dm(i) = alog10(rhte(isite,ifreq)) 
        goto 120
60        dm(i) = phte(isite,ifreq)
        goto 120
70        dm(i) = real(kzy(isite,ifreq))
        goto 120
80        dm(i) = aimag(kzy(isite,ifreq))
        goto 120
90        dm(i) = alog10(rhtm(isite,ifreq)) 
        goto 120
100        dm(i) = phtm(isite,ifreq)
120        continue      
c
c Search the static shifts to see if one applies to this datum
        do 140 j = 1, nsb
          if ((isx1(j).eq.isite) .and. (isx3(j).eq.itype)) then
c Apply the static shift
            dm(i) = dm(i) + shift(j)
          end if
140     continue
c
150      continue
c
c Now compute the sums for the static shift constraints 
      do 160 i = 1, nsb
        if (isfx(i) .ne. 0) then
c This shift is a free parameter: 
          if (isx3(i) .eq. 1) then
            tesum = tesum + shift(i)
          else if (isx3(i) .eq. 5) then
            tmsum = tmsum + shift(i)
          end if
        end if
160   continue
c
c Use the sums to create the synthetic constraint data
      if (iscon.eq.0) then
c no constraint: do nothing
        continue
      else if (iscon .eq. 1) then
c constraint is to sum both TE and TM shifts together 
        dm(ndata+1) = tesum + tmsum
      else if (iscon .eq. 2) then
c constraint is to sum both TE and TM shifts separately 
        dm(ndata+1) = tesum
        dm(ndata+2) = tmsum
      end if
c
      return
      end
c      
c      
c-----------------------------------------------------------------------
      subroutine inputd(datfil)
c-----------------------------------------------------------------------
c
c OCCAM2DMT 1.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 1.02, 13 Jan 1993
c
c inputd opens the data file datfil and reads in the data array
c
c on input:
c   datfil = data file name 
c
c on output:
c   lots of stuff in common blocks (see include files)
c
c calls:  
c   filerr
c
c includes:
      include 'imp.inc'
      include 'mt2D.inc'
c uses nrc, sitnam(), sitloc(), nfre, frek()
c uses idx
      include 'occam.inc'
c uses nd, d(), sd(), dp(), ndp

c arguments:
      character*(*) datfil
c local variables
      logical iffte, ifftm, iffvf
c   flags to indicate which modes are required
      integer lerr, i, j
c   error flag, loop index and handy integer
      character*80  string
c   string buffer
      integer idcode
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
c
c parameters:
      integer iof 
      parameter (iof = 15)
c   iof = io unit number for file operations
c
c-----------------------------------------------------------------------
c there are 3 data parameters in 2D MT (site, freq, data type)
      ndp = 3
c open data file 
      open (iof, file=datfil, status='old', iostat=lerr)
      if (lerr .ne. 0) then
        write(*,*) ' Error opening data file'
        stop
      end if
c read stuff
      read(iof,'(18x,a)', end=198, err=199) string
      if (string(1:16) .ne. 'OCCAM2MTDATA_1.0') call filerr(
     *  ' Data file not supported',iof)
      read(iof,'(18x,a)', end=198, err=199) string
      read(iof,'(18x,a)', end=198, err=199) string
      nrc = idcode(string, iof)
      if (nrc .gt. nr)    call filerr (' Too many sites', iof)
      do 5 i = 1, nrc
5        read(iof,'(a)', end=198, err=199) sitnam(i)
      read(iof,'(a)', end=198, err=199) string
      read(iof,*, end=198, err=199) (sitloc(i), i=1,nrc)
      read(iof,'(18x,a)', end=198, err=199) string
      nfre = idcode(string, iof)
      if (nfre .gt. nf)    call filerr (' Too many frequencies', iof)
      read(iof,*, end=198, err=199) (frek(i), i=1,nfre)
      read(iof,'(18x,a)', end=198, err=199) string
      ndata = idcode(string, iof)
      if (ndata .gt. ndd)    call filerr (' Too many data', iof)
      read(iof,'(a)', end=198, err=199) string
      read(iof,*, end=198, err=199) (dp(i,1), dp(i,2), dp(i,3), d(i), 
     *  sd(i), i=1,ndata)
      
c set idx, the mode definition for PW2DI
      iffte = .false.
      ifftm = .false.
      iffvf = .false.
      do 20 i = 1, ndata
        j = int(dp(i,3)+0.01)
        if ((j.lt.1) .or. (j.gt.6)) then
          write(*,*) ' Invalid data type: ', j
          stop
        end if
        if (j .eq. 1 .or. j .eq. 2) iffte = .true.
        if (j .eq. 3 .or. j .eq. 4) iffvf = .true.
20        if (j .eq. 5 .or. j .eq. 6) ifftm = .true.
c consider vertical field to be TE+TM calculation
      if ((iffte .and. ifftm) .or. iffvf) then
        idx = 1
      else if (iffte) then
        idx = 2
      else if (ifftm) then
        idx = 3
      end if
c 
	write(*,*) ' idebug in inputd ', idebug
      if (idebug .ge. 1) then
        write(*,*) ndata, ' data read'
      end if
c      
c initially set nd for occam as number of MT data
      nd = ndata
c      
      return
c
198      call filerr (' Data file ended prematurely', iof)
199      call filerr (' Error reading data file', iof)
c
      end
      
      
c-----------------------------------------------------------------------
      subroutine inputm(modfil)
c-----------------------------------------------------------------------
c
c OCCAM2DMT 1.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 1.04, 14 Jan 1993
c
c inputm() opens and reads the model file modfil, opens and reads the 
c mesh and static shift files pointed at by modfil, and sets up the 
c environments needed by OCCAM and PW2DI.
c
c Note: inputm() assumes that inputd() has already been called
c
c on input:
c   modfil = (character) model file name 
c on output:
c   incredible amounts of stuff in common blocks (see include files)
c calls:
c   filerr
c
      include 'imp.inc'
      include 'mt2D.inc'
c uses nlay, nrcol(), thick(), width(), nrfix, y(), 
c   icorn(), maxexp, ijexep(), nexep, expen()
c uses sitloc, boffs, nrc, nrl()
c uses dly(), dlz(), yy()
c uses laprm(), nprm
c uses npt()
      include 'occam.inc'
c uses np
c
c arguments:
      character*(*) modfil
c      
c local variables:
      integer irz(nl), iry(nl,mcol), irindx(nl,mcol), lerr, i, l, j
      integer  nsites, lidx, lfre, nexec
      integer ii, jj, inz1, inz2, iny1, iny2
      integer locexp(maxexp,4), nprej, iconb(ipy)
c   irz() = no of vert. FEs in each row
c   iry() = no of horiz FEs in each col for each row
c   irindx() = layer/column to parameter pointers.  0 if fixed block
c   lerr = iostat from file operations
c   i,l,j = loop counters
c   nsites = no of sites read from file.  Ignored after input by OCCAM
c   lidx = type of modelling; should be 0 for OCCAM
c   lfre = number of frequencies read from file. Ignored by OCCAM
c   nexec = type of action required; reset by formod and fordiv
c   ii, jj = index counters for regularization bricks
c   inz1, inz2 = top and  bottom edges of bricks defined by FE mesh
c   iny1, iny2 = left and right edges     ..    ditto          ..
c   locexp() = brick layer/col indices of exeption pairs 
c   nprej = number of entries in model prejudice 
c   iconb() = flag to indicate nodes on conductivity boundaries
      real dummy, tesum, tmsum, sserr
c   dummy = dummy variable to read frequencies into
c   tesum = sum of logs of TE free static shifts 
c   tmsum = sum of logs of TM free static shifts 
c   sserr = error term to use in conjuction with static shift constraint
      character*1 apt(ipz,4,ipy)
c   apt() = character array to store 4 resistivity codes for each element
      character*80 mshfil, mshtyp, string, ssfile, prfile
c   mshfil = mesh file name
c   mshtyp = mesh type
c   string = string buffer for free format reads
c   ssfile = static shift file name
c   prfile = file containing model prejudice and weights
c
      logical iffree
c   iffree = flag to indicate free element in regularization brick
      integer idcode
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
c
c parameters:
      integer iof 
      parameter (iof = 15)
c   iof = io unit number for file operations
c
c-----------------------------------------------------------------------      
cccccccccccccccccccc model file reading cccccccccccccccccccccccccccccccc
c
      open (iof, file=modfil, status='old', iostat=lerr)
      if (lerr .ne. 0) then
        write(*,*) ' Error opening model file'
        stop
      end if
c read character descriptors      
      read (iof,'(18x,a)', end=198, err=199) string
      if (string(1:15) .ne. 'OCCAM2MTMOD_1.0') call filerr(
     *  ' Model file not supported',iof)
      read (iof,'(18x,a)', end=198, err=199) string
      read (iof,'(18x,a)', end=198, err=199) string
      read (iof,'(18x,a)', end=198, err=199) mshfil
      read (iof,'(18x,a)', end=198, err=199) mshtyp
      read (iof,'(18x,a)', end=198, err=199) ssfile
      read (iof,'(18x,a)', end=198, err=199) prfile
c read in the binding offset and number of layers
      read (iof,'(18x,a)', end=198, err=199) string
      boffs = rdcode(string, iof)
      read (iof,'(18x,a)', end=198, err=199) string
      nlay = idcode(string, iof)
      if (nlay .gt. nl) call filerr(' NLAY exceeds dimensions', iof)
c read in regularization mesh description
      do 20 i = 1,nlay
        read (iof,*, end=198, err=199) irz(i),nrcol(i)
        if (nrcol(i) .gt. mcol) 
     *     call filerr (' Too many columns in reg. mesh', iof)
20        read (iof,*, end=198, err=199) (iry(i,j), j=1,nrcol(i))
c read in exceptions to penalty matrix
      read (iof,'(18x,a)', end=198, err=199) string
      nexep = idcode(string, iof)
      if (nexep .gt. maxexp) 
     *  call filerr (' Too many penalty exeptions', iof)
      do 21 i = 1,nexep
21      read (iof,*, end=198, err=199) (locexp(i,j),j=1,4),expen(i)
c
c we are done with the model file
      close (iof)
c
ccccccccccccccccccccc mesh file reading cccccccccccccccccccccccccccccccc
c
      if (mshtyp(1:4) .ne. 'PW2D') call filerr(
     *  ' Mesh type not supported', iof)
      open (iof, file=mshfil, status='old', iostat=lerr)
      if (lerr .ne. 0) then
        write(*,*) ' Error opening mesh file'
        stop
      end if
c      
c read in stuff
      read (iof,'(a)', end=198, err=199) string
      read (iof,*, end=188, err=189) lidx,nodey,nodez,nrfix,lfre,nexec
c
c screen input for errors:      
      if (lidx .ne. 0) call filerr(' Mesh file not OCCAM file', iof)
      if (nodey .gt. ipy) call filerr(' NODEY exceeds dimensions', iof)
      if (nrfix .ge. 36) call filerr('Too many fixed resistivities',iof)
c test size of nodez
      if (nodez+10 .gt. ipz) call filerr('NODEZ exceeds dimensions',iof)      
c input values for fixed resistivities, if any
      if (nrfix .ge. 1) then 
        read (iof,*, end=188, err=189) (y(i), i = 2,nrfix+1)
      end if
c input frequencies, if any (these are ignored)
      if (lfre .ge. 1) then
        read (iof,*, end=188, err=189) (dummy, i=1,lfre)
      end if
c input the element widths and heights
      read (iof,*, end=188, err=189) (dly(i), i = 1,nodey-1)
      read (iof,*, end=188, err=189) (dlz(i), i = 1,nodez-1)
c input site information, if any (these are ignored)      
      read (iof,*, end=188, err=189) nsites, (dummy, i = 1,nsites)
c input the resistivity codes
      do 25 i = 1, nodez-1
        do 25 l = 1,4
25        read (iof,'(170a1)',end=188,err=189) (apt(i,l,j), j=1,nodey-1)
      
c We are done with the mesh file
      close (iof)
c
cccccccccccccccccccc statics file reading cccccccccccccccccccccccccccccc
c
      if (ssfile(1:4) .eq. 'none' .or. ssfile(1:4) .eq. 'NONE') then
        iscon = 0
        nsb = 0
      else
        open (iof, file=ssfile, status='old', iostat=lerr)
        if (lerr .ne. 0) then
          write(*,*) ' Error opening static shift file'
          stop
        end if
c Read the stuff in it.
        read (iof,'(18x,a)', end=208, err=209) string
        if (string(1:17) .ne. 'OCCAM2MTSHIFT_1.0') call filerr(
     *    ' Static shift file not supported', iof)
        read (iof,'(18x,a)', end=208, err=209) string
        read (iof,'(18x,a)', end=208, err=209) string
        read (iof,'(18x,a)', end=208, err=209) string
        iscon = idcode(string, iof)
        read (iof,'(18x,a)', end=208, err=209) string
        sserr = rdcode(string, iof)
        read (iof,'(18x,a)', end=208, err=209) string
        nsb = idcode(string, iof)
        read (iof,'(a)', end=208, err=209) string
        do 30 i = 1,nsb
30        read (iof,*,end=208,err=209) isx1(i),isx3(i),shift(i),isfx(i)
        close (iof)
      end if
c
cccccccccccccccccccc prejudice file reading cccccccccccccccccccccccccccc
c
      if (prfile(1:4) .eq. 'none' .or. prfile(1:4) .eq. 'NONE') then
c no prejudice: zero out the model pejudice and associated weights
        do 40 i = 1,npp
          premod(i) = 0.0
40        iprewt(i) = 0
        nprej = 0
      else
c read the prejudice file
        open (iof, file=prfile, status='old', iostat=lerr)
        if (lerr .ne. 0) then
          write(*,*) ' Error opening model prejudice file'
          stop
        end if
c Read the stuff in it.
        read (iof,'(18x,a)', end=218, err=219) string
        if (string(1:16) .ne. 'OCCAM2MTPREJ_1.0') call filerr(
     *    ' Model prejudice file not supported', iof)
        read (iof,'(18x,a)', end=218, err=219) string
        read (iof,'(18x,a)', end=218, err=219) string
        nprej = idcode(string, iof)
        if (nprej .gt. npp)      call filerr(' Too much prejudice', iof)
        read (iof,*,end=218,err=219) (premod(i), i=1,nprej)
        read (iof,*,end=218,err=219) (iprewt(i), i=1,nprej)
        close (iof)
      end if
c
cccccccccccccccc set up mesh and penalty matrix info ccccccccccccccccccc
c
c npt() runs from 0 (0) to 35 (Z), '?' is flagged as -99, to be overwritten 
c from the parameter array later.
      do 50 j = 1,nodey-1
        do 50 l = 1,4
          do 50 i = 1,nodez-1
            npt(i,l,j) = ichar(apt(i,l,j)) - 48
            if (npt(i,l,j) .eq. 15) npt(i,l,j) = -99
            if (npt(i,l,j) .gt. 9) npt(i,l,j) = npt(i,l,j) - 7
50      continue
c
c Set up the free parameter part of npt() and also laprm().
c The following block of code finds the outline of every parameter brick,
c and if any part of it is free, adds it to npt() and laprm().
c This block of code also computes the widths and thickness of the bricks,
c which will be needed later during the computation of the penalty matrix.
c
c initialize the number of parameters and the top edge of bricks
      nprm = 0
      inz1 = 1
c for every layer in the reg. mesh
      do 80 i = 1,nlay
c find the bottom edge of bricks and initialize the right edge of first brick
        inz2 = inz1 + irz(i) - 1
        iny1 = 1
c for every brick in that layer 
        do 70 j = 1,nrcol(i)
c set irindx() assuming that this brick is fixed (will be tested later)
          irindx(i,j) = 0
c find the right edge of brick
          iny2 = iny1 + abs(iry(i,j)) - 1
c test to see if the brick is free 
          iffree = .false.
            do 52 ii = inz1,inz2
              do 52 jj = iny1, iny2
                do 52 l = 1,4
52                  if (apt(ii,l,jj) .eq. '?') iffree = .true.
          if (iffree) then
c it is: increment the parameter count and set a new element of laprm
            nprm = nprm + 1
            laprm(nprm) = nprm + nrfix 
c set irindx() to point to the new parameter
            irindx(i,j) = nprm
c find the thickness of the brick
            thick(nprm) = 0.0
            do 55 ii = inz1, inz2
55              thick(nprm) = thick(nprm) + dlz(ii)
c find the width of the brick
            width(nprm) = 0.0
            do 57 jj = iny1, iny2
57              width(nprm) = width(nprm) + dly(jj)
c record the FE nodes for the four corners of the block (top, left, bottom, right)
            icorn(nprm,1) = inz1
            icorn(nprm,2) = iny1
            icorn(nprm,3) = inz2+1
            icorn(nprm,4) = iny2+1
c over the entire brick, including sub-triangles,
c test to see if triangular element if free and set npt() if it is
            do 60 ii = inz1,inz2
              do 60 jj = iny1, iny2
                do 60 l = 1,4
                  if (apt(ii,l,jj) .eq. '?') then
                    npt(ii,l,jj) = nprm + nrfix 
                  end if
60            continue
          end if
c increment the left edge of brick
          iny1 = iny2 + 1
70        continue
c increment the top edge of bricks
        inz1 = inz2 + 1
80      continue
c
c copy /blki/nprm into /rgrid/nbrick to give my routines the number of parameters
c which are actually resistivity bricks (c.f. static shifts)
      nbrick = nprm
c
c Now, the bricks on the edges of the mesh and the bottom of the mesh will have
c widths and thicknesses that are too large to use in regularization penalty.  Set
c them equal to their inside neighbours.  The indexing in this code is the same as 
c above.
      do 85 i = 1,nlay
        do 85 j = 1, nrcol(i)
          ii = irindx(i,j)
          if (ii .gt. 0) then
            if (j .eq. 1) width(ii) = width(ii+1)
            if (j .eq. nrcol(i)) width(ii) = width(ii-1)
            if (i .eq. nlay) thick(ii) = thick(ii-nrcol(i))
          end if
85    continue
c 
c search the penalty exceptions list to find the block numbers
      do 97 i = 1,nexep
        ijexep(i,1) = irindx(locexp(i,1),locexp(i,2)) 
        ijexep(i,2) = irindx(locexp(i,3),locexp(i,4)) 
        if ((ijexep(i,1).le.0) .and. (ijexep(i,2).le.0)) call filerr( 
     * ' Penalty requested for two fixed blocks; use prejudice file',0)
97    continue
c 
cccccccccccccccccccc set up receiver locations ccccccccccccccccccccccccc
c
c first find the position of the right edge of the left terminating layer: 
      dummy = 0.0
      do 95 i = 1, abs(iry(1,1))
95         dummy = dummy + dly(i)
c now find cumulative positions of all nodes (borrowing an array 
c used for similar thing in PW2DI), zero out iconb() 
      yy(1) = boffs - dummy
      iconb(1) = 0
      do 100 i = 1,nodey-1
        iconb(i+1) = 0
100        yy(i+1) = yy(i) + dly(i)
c find the node locations for surface conductivity boundaries
      l = 1
      do 103 i = 1, nrcol(1)
        l = l + iry(1,i)
        iconb(l) = 1
103      continue
c now find the nodes closest to the station locations, avoiding conductivity
c boundaries
      do 110 i = 1, nrc
        do 105 j = 2,nodey
          if (yy(j) .ge. sitloc(i)) then
            if (yy(j)-sitloc(i) .le. sitloc(i)-yy(j-1)) then
              nrl(i) = j
              if (iconb(j) .eq. 1) nrl(i) = j-1
            else
              nrl(i) = j-1
              if (iconb(j-1) .eq. 1) nrl(i) = j
            end if
            goto 110
          end if
105        continue
110      continue
      if (idebug .ge. 1) then
        write(*,*) 'station position    node number     actual position'
        do 111 i = 1, nrc
111          write(*,*) sitloc(i), nrl(i), yy(nrl(i))
      end if
c now check to see if two stations occupy the same location
      do 115 i = 1, nrc
        do 114 j = 1, nrc
          if ((nrl(i) .eq. nrl(j)) .and. (i.ne.j)) then
            call filerr(' Two stations on same node', 0)
          end if
114        continue
115      continue
c zero out the rest of the array (because they might be referenced by pw2di)
      do 120 i = nrc+1, nr
120        nrl(i) = 0
c
cccccccccccccccccccc set up static shift arrays cccccccccccccccccccccccc
c
c Examine the statics blocks, allocate new parameters and set up the
c data constriants
      nfss = 0
      tesum = 0.0
      tmsum = 0.0
      if (nsb .eq. 0) iscon = 0
      do 140 i = 1,nsb
        if (isfx(i) .ne. 0) then
c this static shift is a free parameter: add one more to the list
          nfss = nfss + 1
          if (nbrick + nfss .gt. npp) call filerr(
     *     ' No room for free statics in parameter array',0)
          isfx(i) = nbrick + nfss
c add a term to the constraint sums
          if (isx3(i) .eq. 1) then
            tesum = tesum + shift(i)
          else if (isx3(i) .eq. 5) then
            tmsum = tmsum + shift(i)
          else 
            call filerr(' Can only apply statics to resistivities',0)
          end if
        end if
140   continue
c
c Set up the statics constraint
      if (iscon.eq.0) then
c no constraint: do nothing
        continue
      else if (iscon .eq. 1) then
c constraint is to sum both TE and TM shifts together to whatever was the
c sum in the statics file
        nd = ndata + 1
        if (nd .gt. ndd) call filerr(  
     *     ' No room for statics constraint in data array', 0)
        dp(ndata+1,3) = 7
        d(ndata+1) = tesum + tmsum
        sd(ndata+1) = sserr
      else if (iscon .eq. 2) then
c constraint is to sum both TE and TM shifts separately to whatever was the
c sum in the statics file
        nd = ndata + 2
        if (nd .gt. ndd) call filerr( 
     *    ' No room for statics constraint in data array',0)
        dp(ndata+1,3) = 8
        dp(ndata+2,3) = 9
        d(ndata+1) = tesum
        d(ndata+2) = tmsum
        sd(ndata+1) = sserr
        sd(ndata+2) = sserr
      else
        call filerr (
     *   ' Statics constraints type greater than 2 not supported', 0)
      end if
c
      np = nbrick + nfss
c
ccccccccccccccccccccccc clean up and leave cccccccccccccccccccccccccccc
c
c We should have read the iteration file already and set np.  Test for
c consistency with nbrick and nfss
      if (np .ne. nbrick+nfss) 
     *  call filerr (' Model file and NP do not match', 0)
      if ((np .ne. nprej) .and. (nprej .ne. 0)) 
     *  call filerr (' Prejudice file and NP do not match', 0)
c
c that's it
      if (idebug .ge. 1) then
        write(*,*) nodey, ' x ', nodez, ' FE mesh read'
        write(*,*) np, ' model parameters defined'
      end if
c
      return
c
188      call filerr (' Mesh file ended prematurely', iof)
189      call filerr (' Error reading mesh file', iof)
198      call filerr (' Model file ended prematurely', iof)
199      call filerr (' Error reading model file', iof)
208      call filerr (' Statics file ended prematurely', iof)
209      call filerr (' Error reading statics file', iof)
218      call filerr (' Prejudice file ended prematurely', iof)
219      call filerr (' Error reading prejudice file', iof)
c
      end
      
      
C **********************************************************************
      SUBROUTINE PW2DI(NRES,NEXEC)
C **********************************************************************
      INCLUDE 'imp.inc'
C
C========================== INCLUDE FILES ==============================
C
      INCLUDE 'mt2D.inc'
C
C============================ INTERFACE ================================
C
C WRITTEN BY:  P. WANNAMAKER, UURI.
C
C DESCRIPTION: 
c     This version for forward computation only
C
C     THIS SUBROUTINE INVERTS MAGNETOTELLURIC DATA TO GIVE 2-DIMENSIONAL
C     CONDUCTIVITY SECTIONS. FINITE ELEMENTS ARE USED WITH LINEAR BASIS
C     FUNCTIONS. MKS UNITS AND EXP(+JWT) TIME DEPENDENCE ARE ASSUMED.
C     THE + COORDINATE DIRECTIONS FOR THE PROBLEM ARE X-NORTH, Y-EAST,
C     Z-DOWN. THE STRIKE DIRECTION IS ASSUMED N-S. THE ORIGIN IS AT THE
C     CENTER OF THE MESH AT THE AIR-EARTH INTERFACE.  TOPOGRAPHIC OR
C     BATHYMETRIC VARIATIONS ACROSS THE PROFILE CAN BE INCLUDED ALSO.
C     A SECONDARY FIELD FORMULATION OF THE FINITE ELEMENT EQUATIONS HAS
C     BEEN UTILIZED HERE TO IMPROVE ACCURACY (WANNAMAKER ET AL, 1987,
C     GEOP. J. ROY. ASTR. SOC., 88).  THE LAYERED HOST FOR THE INHOMO-
C     GENEITY BY CONVENTION IS TAKEN TO BE THAT AT THE LEFT EDGE OF THE
C     MESH.  TOPOGRAPHIC MODEL STUDIES USING THE FORWARD PROGRAM ARE
C     PRESENTED BY WANNAMAKER ET AL (1986, GEOPHYSICS, 51(11)).
C
C     ROUTINES NEEDED:  PW2DI, PARFLD, COAMP, ZLFLD, GAUSS, AUXFLD,
C                       BNORM, DERIV, PARJAC, GAUSSD
C
C ON INPUTS:
C         NRES - NO. OF UNIQUE RESISTIVITY MEDIA IN THE MESH
C                (INCLUDES ANY SEAWATER BUT NOT THE AIR)
C         NFRE - NO. OF FREQUENCIES TO BE MODELED
C         NEXEC - IF .EQ. 0, JUST PRINT INPUT PARAMETERS AND MESH
C                 TO CHECK CORRECTNESS
C               - IF .EQ. 1, COMPUTE RESPONSE
C               - IF .EQ. 2, INVERT DATA
C         Y - ARRAY OF RESISTIVITIES OF MEDIA MAKING UP THE MESH;
C             LAST VALUE IS THAT OF ANY SEAWATER.
C         FREK - ARRAY OF FREQUENCIES
C         NRL - SITE LOCATIONS BY FE MESH NODE NUMBER
C
C         /BLKI/ IDX - PARAMETER DEFINING MODE OF PROBLEM TO BE SOLVED
C             - IF .EQ. 1, SOLVE BOTH TE AND TM
C             - IF .EQ. 2, SOLVE TE ONLY
C             - IF .EQ. 3, SOLVE TM ONLY
C         /BLKI/ NODEY - NO. OF NODES IN THE Y-DIRECTION OF THE FINITE
C                 ELEMENT MESH
C         /BLKI/ NODEZ - NO. OF NODES IN THE Z-DIRECTION OF THE FINITE
C                 ELEMENT MESH
C         /BLKI/ NPRM - NO. OF POLYGONAL MEDIA OF INVERSION MODEL
C         /BLKI/ LAPRM - INTEGER LABELS OF EACH INVERSION MEDIUM
C         /BLKD/ DLY - ARRAY OF ELEMENT COLUMN WIDTHS
C         /BLKD/ DLZ - ARRAY OF ELEMENT ROW HEIGHTS
C         /BLKE/ NPT - INTEGER LABELS FOR EACH TRIANGULAR ELEMENT OF EACH 
C               RECTANGULAR ELEMENT OF THE MESH.  THESE LABELS ARE 
C               ASSIGNED THE VALUES IN THE RESISTIVITY ARRAY Y.  
C               LABEL 0 IS RESERVED FOR AIR AND LABEL -1 IS RESERVED
C               FOR SEAWATER.  TRIANGULAR ELEMENT LABELING ORDER IS 
C               TOP, LEFT, BOTTOM AND RIGHT.
C
C         STORAGE PARAMETERS:
C         IPY - MAX. NO. OF NODES IN THE Y-DIRECTION.  
C         IPZ - MAX. NO. OF NODES IN THE Z-DIRECTION. 
C         IPI - MAX. NO. OF POLYGONAL MEDIA IN AN INVERSION MODEL. 
C
C============== PARAMETER AND LOCAL VARIABLE DECLARATIONS ==============
C
C     ---- PASS PARAMETERS
      INTEGER NRES,NEXEC
C
C     ---- LOCAL VARIABLES
c rj converted to local variable for dynamic memory allocation
c      COMPLEX RJ(IPI,IP2)
      COMPLEX RJ(1)
      COMPLEX FT,PF,VF,HF,TP,ZIMP,FTD,ZIMPD,KZYD
      INTEGER lidx,NODEY1,MAIR,NODEZ1
      INTEGER NTR,I,J,K,L,M,III,NDX,NMODE,LF,NZZN,NCR,NTO,IA,KP,NYRB
      REAL RAD,XEO,XMUO,AIOR,conds(ipi+36)
      COMPLEX S(IP2,IP3),R(IP2)
C
C
C=======================================================================
C
C     ----ZERO OUT LOCAL S(), R() AND IF INVERTING, PARAMETER RJ()
      DO 1000 J=1,IP3
        DO 1000 I=1,IP2
          S(I,J)=CMPLX(0.,0.)
 1000 CONTINUE
      DO 1001 J=1,IP2
        R(J)=CMPLX(0.,0.)
 1001 CONTINUE
C        
      NODEY1 = NODEY - 1
C
C     ---- DEFINE LOCATION OF Y=0 AT CENTER OF MESH
      L = 2
      YY(1) = 0.0
      DO 20 I = 1,NODEY1
         YY(L) = YY(L-1) + DLY(I)
         L = L + 1
20    CONTINUE
C
      DO 30 I = 1,NODEY
        YY(I) = YY(I) - YY(NODEY)/2.
30    CONTINUE
c
c add the air layer if necessary
      lidx = idx
      call addair(lidx, mair)
      NODEZ1 = NODEZ - 1
      NTR = MAIR + 1
C
C     ---- CHANGE FROM RESISTIVITIES INPUT TO CONDUCTIVITIES AND
C          ASSOCIATE CONDUCTIVITIES WITH THE MESH INPUT CODE
      RA = 1.0E+30
      conds(1) = 1.0E-30
      DO 90 I = 2,NRES+1
        conds(i) = 1./Y(I)
90    CONTINUE
C
      ISEA=NRES+1
      conds(ISEA+1) = 1./Y(NRES+1)
      DO 100 J = 1,NODEZ1
        DO 101 L = 1,NODEY1
          DO 102 M = 1,4
            K = NPT(J,M,L)
            SIG(J,M,L) = conds(K+1)
102       CONTINUE
101     CONTINUE
100   CONTINUE
C
C     ---- POSITION THE DATUM (HIGHEST TOPOGRAPHIC POINT)
      L = NTR
      ZZ(L) = 0.
      DO 110 I = NTR,NODEZ1
        ZZ(L+1) = ZZ(L) + DLZ(I)
        L = L + 1
110   CONTINUE
C
      IF(MAIR .NE. 0)THEN
        L = NTR
        DO 130 III = 1,MAIR
          I = MAIR + 1 - III
          ZZ(L-1) = ZZ(L) - DLZ(I)
          L = L - 1
130     CONTINUE
      END IF
C
C     ---- SPECIFY NODAL ROW DEPRESSIONS OF TOPO/BATH SURFACE.
      DO 140 I = 2,NODEY1
        NTAO(I) = 0
        DO 141 J = NTR,NODEZ1
          IF(NPT(J,4,I-1) .EQ. 0  .AND. NPT(J,2,I) .EQ. 0 )GO TO 142
          IF(NPT(J,4,I-1) .EQ. ISEA .AND. NPT(J,2,I) .EQ. ISEA)GO TO 142
          GO TO 140
142       CONTINUE
          NTAO(I) = J - MAIR
141     CONTINUE
140   CONTINUE
      NTAO(1) = NTAO(2)
      NTAO(NODEY) = NTAO(NODEY1)
C
C     ---- DISCERN BOUNDARY LAYERING PARAMETERS FROM INPUT MESH
      NYRB = 1
      DO 150 I = 1,2
        IF(I .EQ. 2)NYRB = NODEY1
        NOLA(I) = 1
        NZZB(I) = NTR
        DO 160 J = NTR,NODEZ1-1
C
C         ---- PATCH TO PREVENT ARRAY OVERRUN
          IF (NOLA(I) .GE. NL)GOTO 160
          IF(NPT(J,1,NYRB) .EQ. 0)THEN
            NZZB(I) = J + 1
          ELSE
            RLR(NOLA(I),I) = 1./SIG(J,2*I,NYRB)
            IF(NPT(J,2*I,NYRB) .NE. NPT(J+1,2*I,NYRB))THEN
              DLR(NOLA(I),I) = ZZ(J+1) - ZZ(NZZB(I))
              NOLA(I) = NOLA(I) + 1
            END IF
          END IF
160     CONTINUE
        RLR(NOLA(I),I) = 1./SIG(NODEZ1,2*I,NYRB)
150   CONTINUE
C
      PI = 3.1415927
      RAD = 180./PI
      XEO = 8.85433E-12
      XMUO = PI*4.E-07
      NDX = IDX
      NMODE = 1
      IF(NDX .EQ. 1)NMODE = 2
      DO 180 LF = 1,NFRE
        W = 2.*PI*FREK(LF)
        WU = W*XMUO
        WE = W*XEO
        IF(NDX .EQ. 1)IDX = 2
        NTQ = 1
        NDZ = NODEZ
        NDZ1 = NDZ - 1
        DO 190 I = 1,NODEY
          NTA(I) = NTAO(I) + NTR + NDZ*(I-1)
190     CONTINUE
C
        DO 200 IMD = 1,NMODE
          IF(IMD .EQ. 2)IDX = 3
          IF(IDX .EQ. 3)THEN
            NTQ = NTR
            NDZ = NODEZ - MAIR
            NDZ1 = NDZ - 1
            DO 210 I = 1,NODEY
              NTA(I) = NTAO(I) + 1 + NDZ*(I-1)
210         CONTINUE
          END IF
C
C         ---- FORM AND SOLVE GLOBAL SYSTEM OF EQUATIONS
          CALL PARFLD(S,R)
C
C         ---- CALCULATE AUXILIARY FIELDS, REAL AND IMAGINARY, 
C              AND APPARENT RESISTIVITY AND PHASE OF THE 
C              IMPEDANCE FOR PRINTOUT.
C              WRITE FORWARD PROBLEM HEADER FIRST.
          NZZN = NTA(1)
          NCR = 1
          DO 260 I = 1,NODEY
            AFN(I) = (0.,0.)
            AFT(I) = (0.,0.)
            IF(I .EQ. NRL(NCR))THEN
              CALL AUXFLD(I,0,R,RJ)
              NTO = NTAO(I) + NTR
              FT = R(NTA(I)) + FPA(NTO)
              PF = FT/FPA(NZZN)
              HF = AFT(I)/FPB(NZZN)
              VF = AFN(I)/FPB(NZZN)
              TP = VF/HF
              IF(IDX .EQ. 2)THEN
                ZIMP = FT/AFT(I)
                RHTE(NCR,LF) = CABS(ZIMP)**2/WU
                PHTE(NCR,LF) = RAD*ATAN2(AIMAG(ZIMP),REAL(ZIMP))
                KZY(NCR,LF) = TP
              ELSE
                ZIMP = AFT(I)/FT
                RHTM(NCR,LF) = CABS(ZIMP)**2/WU
                PHTM(NCR,LF) = 180. + RAD*ATAN2(AIMAG(ZIMP),REAL(ZIMP))
              END IF
              NCR = NCR + 1
            END IF
260       CONTINUE
C
C         ---- CALCULATE DATA-PARAMETER JACOBIAN 
C              IF INVERSION IS SELECTED.
C              STRIKE COMPONENT COMPUTED IN PARJAC FROM 
C              ANALYTIC SOURCES WITH FORWARD-BACK SUBSTITUTIONS 
C              IN GAUSSD.
200     CONTINUE
180   CONTINUE
      IDX=NDX
C
C     ----PUT BIG VALUES INTO OVERLAYED PARAMETER RJ() TO
C         MAKE ANY OVERLAY PROBLEMS OBVIOUS.
C
c remove the air layer if necessary
      call addair(0, mair)

      RETURN
      END
C
C **********************************************************************
C **********************************************************************
      SUBROUTINE PARFLD(S,R)
      INCLUDE 'imp.inc'
C
C========================== INCLUDE FILES ==============================
C
      INCLUDE 'mt2D.inc'
C
C============================ INTERFACE ================================
C
C WRITTEN BY:  P. WANNAMAKER, UURI.
C
C DESCRIPTION: 
C
C     THIS ROUTINE COMPUTES THE GLOBAL FINITE ELEMENT SYSTEM TO GIVE
C     THE FIELD PARALLEL TO STRIKE.  TASKS PERFORMED INCLUDE THE
C     PRIMARY FIELDS FOR THE SOURCE VECTOR AND AT EACH NODAL COLUMN,
C     LOADING THE MATRIX ELEMENTS INTO THE GLOBAL SYSTEM, INCORPOR-
C     ATING THE DIRICHLET BOUNDARY CONDITIONS, AND SOLVING THE SYSTEM.
C
C============== PARAMETER AND LOCAL VARIABLE DECLARATIONS ==============
C
C     ---- PASS PARAMETERS
      COMPLEX S(IP2,IP3),R(IP2)
C
C     ---- LOCAL VARIABLES
      COMPLEX CP(4),CP12,CP23,CP34,CP14,CP1234,A11,A12,A13,A15,A22,A24,
     &A25,A33,A34,A35,A44,A45,A55,E11,G11,C11,H11,E22,H22,C22,E33,G33,
     &E44,ZLFLD,FSE,FSH
      COMPLEX*16 CK(4),CK12,CK13,CK14,CK23,CK24,CK34
      COMPLEX FP1(IPZ),FP2(IPZ),FP3(IPZ),FP4(IPZ),DJ1,DJ2,DJ3,DJ4,
     &DSS1,DSS2,DSS3,DSS4,FP5(IPZ),FP6(IPZ),FP7(IPZ),FP8(IPZ),SM1,SM2,
     &SM3,SM4,SM5,RU1,RU2,RU3,RU4,RU5,RC1,RC2,RC3,RC4
      INTEGER NODEY1,NODEZ1,NBAND,NNODE,NLYR1
      INTEGER I,J,K,L,M,L1,L2,L3,L4,IYN,IT,LZP,NBAN1,NJI,NDB,JM1
      REAL RHA,PHA,ZT,Z2,Z3,Z5,Z6,Z7,Z8,ZL,DLZY,DLYZ
      REAL A,B,C,C4,DY4
C
C=======================================================================
C
      NODEY1 = NODEY - 1
      NODEZ1 = NODEZ - 1
      NBAND = NDZ + 2
      NNODE = NODEY*NDZ
      NDB = NNODE - NBAND
      NBAN1 = NBAND - 1
C
C     THIS SECTION CALCULATES THE PRIMARY E- OR H-FIELDS FOR THE TE AND
C     TM PROBLEMS.  FIELDS IN THE LAYERED HOST ARE CALCULATED FOR AN
C     INCIDENT E-FIELD AMPLITUDE OF (1.,0.) AT THE EARTH'S SURFACE (NOT
C     AT THE TOP OF THE MESH) POLARIZED IN THE +X DIRECTION FOR THE TE
C     MODE AND IN THE +Y DIRECTION FOR THE TM MODE.  ZLFLD(Z,1) IS FOR
C     E AND ZLFLD(Z,2) IS FOR H.
C
      WM = WU
      WP = WE
      NLYR = NOLA(1)
      NLYR1 = NLYR - 1
      IF(NLYR .EQ. 1)GO TO 10
       DO 20 I = 1,NLYR1
        DDPP(I) = DLR(I,1)
        RL(I) = RLR(I,1)
20     CONTINUE
10    CONTINUE
      RL(NLYR) = RLR(NLYR,1)
C
      CALL COAMP
      IF(IMD .EQ. 1)THEN
       FSE = ZLFLD(0.,1)
       FSH = ZLFLD(0.,2)
       RHA = CABS(FSE/FSH)**2/WU
       PHA = 57.2958*ATAN2(AIMAG(FSE/FSH),REAL(FSE/FSH))
      END IF
C
      DO 30 I = NTQ,NODEZ1
       ZT = ZZ(I) - ZZ(NZZB(1))
       Z1 = ZT + DLZ(I)/6.
       Z2 = ZT + DLZ(I)/2.
       Z3 = ZT + DLZ(I)*(5./6.)
       LZP = 0
       IF(Z2 .GE. 0.)LZP = 1
       IF(NLYR .EQ. 1)GO TO 40
        DO 50 J = 2,NLYR
         JM1 = J - 1
         IF(Z2 .GT. DDPP(JM1))LZP = J
50      CONTINUE
40     CONTINUE
       IF(LZP .EQ. 0)GO TO 60
        SIGH(I) = CMPLX(1./RL(LZP),WE)
        GO TO 70
60     CONTINUE
        SIGH(I) = CMPLX(1./RA,WE)
70     CONTINUE
       IF(IDX .EQ. 3)GO TO 80
        FP1(I) = ZLFLD(Z1,1)
        FP2(I) = ZLFLD(Z2,1)
        FP3(I) = ZLFLD(Z3,1)
        FP4(I) = FP2(I)
        GO TO 31
80     CONTINUE
        FP1(I) = -ZLFLD(Z1,2)
        FP2(I) = -ZLFLD(Z2,2)
        FP3(I) = -ZLFLD(Z3,2)
        FP4(I) = FP2(I)
        Z5 = ZT
        Z6 = ZT + DLZ(I)/4.
        Z7 = ZT + DLZ(I)*(3./4.)
        Z8 = ZT + DLZ(I)
        FP5(I) = ZLFLD(Z5,1)
        FP6(I) = ZLFLD(Z6,1)
        FP7(I) = ZLFLD(Z7,1)
        FP8(I) = ZLFLD(Z8,1)
31     CONTINUE
30    CONTINUE
      DO 90 I = NTQ,NODEZ
       ZL = ZZ(I) - ZZ(NZZB(1))
       IF(IDX .EQ. 2)THEN
        FPA(I) = ZLFLD(ZL,1)
        FPB(I) = ZLFLD(ZL,2)
       ELSE
        FPA(I) = -ZLFLD(ZL,2)
        FPB(I) = ZLFLD(ZL,1)
       END IF
90    CONTINUE
C
C     THIS SECTION CALCULATES A 4 X 4 ELEMENT MATRIX FOR A RECTANGULAR
C     ELEMENT MADE UP OF 4 TRIANGULAR ELEMENTS. THE 5 X 5 MATRIX IS
C     OBTAINED FROM PROPER COMBINATION OF THE FOUR 3 X 3 TRIANGULAR
C     ELEMENT MATRICIES AND REDUCED TO A 4 X 4 MATRIX THROUGH STATIC
C     CONDENSATION. THIS ELIMINATES THE DEGREE OF FREEDOM ASSOCIATED
C     WITH THE INTERNAL NODE OF THE RECTANGLE, THUS REDUCING THE
C     SIZE OF THE GLOBAL SYSTEM MATRIX.
C
      DO 110 I = 1,NNODE
       R(I) = (0.,0.)
110   CONTINUE
      DO 120 K = 1,NBAND
       DO 130 I = 1,NNODE
        S(I,K) = (0.,0.)
130    CONTINUE
120   CONTINUE
C
      L1 = 0
      DO 140 L = 1,NODEY1
       L1 = L1 + 1
       DO 150 M = NTQ,NODEZ1
        DLZY = DLZ(M)/(2.*DLY(L))
        DLYZ = DLY(L)/(2.*DLZ(M))
        A = -(DLZY + DLYZ)/2.
        B = (DLZY - DLYZ)/2.
        C = DLY(L)*DLZ(M)/48.
        C4 = C*4.
        DY4 = DLY(L)/4.
C
C     CALCULATE THE APPROPRIATE COEFFICIENTS OF THE HELMHOLTZ
C     EQUATION FOR INSERTION INTO THE ELEMENT MATRICIES
C     ACCORDING TO WHETHER THE TE OR TM MODE IS BEING SOLVED.
C
        IF(IDX .EQ. 2)THEN
         DO 160 I = 1,4
          CK(I) = CMPLX(0.,WU)
          CK(I) = 1./CK(I)
          CP(I) = -CMPLX(SIG(M,I,L),WE)
160      CONTINUE
        ELSE
         DO 170 I = 1,4
          CK(I) = CMPLX(SIG(M,I,L),WE)
          CK(I) = 1./CK(I)
          CP(I) = -CMPLX(0.,WU)
170      CONTINUE
        END IF
C
C     THESE ARE THE ELEMENTS OF THE 5 X 5 MATRIX.
C
        CK12 = CK(1) + CK(2)
        CK13 = CK(1) + CK(3)
        CK14 = CK(1) + CK(4)
        CK23 = CK(2) + CK(3)
        CK24 = CK(2) + CK(4)
        CK34 = CK(3) + CK(4)
        CP12 = CP(1) + CP(2)
        CP23 = CP(2) + CP(3)
        CP34 = CP(3) + CP(4)
        CP14 = CP(1) + CP(4)
        CP1234 = CP12 + CP34
        A11 = A*CK12 + 2.*C*CP12
        A12 = -B*CK(2) + C*CP(2)
        A13 = B*CK(1) + C*CP(1)
        A15 = DLYZ*CK(1) + DLZY*CK(2) + C*CP12
        A22 = A*CK23 + 2.*C*CP23
        A24 = B*CK(3) + C*CP(3)
        A25 = DLYZ*CK(3) + DLZY*CK(2) + C*CP23
        A33 = A*CK14 + 2.*C*CP14
        A34 = -B*CK(4) + C*CP(4)
        A35 = DLYZ*CK(1) + DLZY*CK(4) + C*CP14
        A44 = A*CK34 + 2.*C*CP34
        A45 = DLYZ*CK(3) + DLZY*CK(4) + C*CP34
        A55 = 2.*(C*CP1234 - DLYZ*CK13 - DLZY*CK24)
C
C     THESE ARE THE ELEMENTS OF THE 5 X 1 SOURCE VECTOR
C
        DJ1 = (CMPLX(SIG(M,1,L),WE) - SIGH(M))*FP1(M)
        DJ2 = (CMPLX(SIG(M,2,L),WE) - SIGH(M))*FP2(M)
        DJ3 = (CMPLX(SIG(M,3,L),WE) - SIGH(M))*FP3(M)
        DJ4 = (CMPLX(SIG(M,4,L),WE) - SIGH(M))*FP4(M)
        SM1 = (0.,0.)
        SM2 = (0.,0.)
        SM3 = (0.,0.)
        SM4 = (0.,0.)
        SM5 = (0.,0.)
        IF(IDX .EQ. 2)GO TO 180
         DSS1 = 1. - SIGH(M)/CMPLX(SIG(M,1,L),WE)
         DSS2 = 1. - SIGH(M)/CMPLX(SIG(M,2,L),WE)
         DSS3 = 1. - SIGH(M)/CMPLX(SIG(M,3,L),WE)
         DSS4 = 1. - SIGH(M)/CMPLX(SIG(M,4,L),WE)
         DJ1 = CMPLX(0.,WU)*DSS1*FP1(M)
         DJ2 = CMPLX(0.,WU)*DSS2*FP2(M)
         DJ3 = CMPLX(0.,WU)*DSS3*FP3(M)
         DJ4 = CMPLX(0.,WU)*DSS4*FP4(M)
         SM1 = DSS1*(2.*FP5(M) - FP6(M)) + DSS2*FP6(M)
         SM2 = -DSS2*FP7(M) - DSS3*(2.*FP8(M) - FP7(M))
         SM3 = DSS1*(2.*FP5(M) - FP6(M)) + DSS4*FP6(M)
         SM4 = -DSS3*(2.*FP8(M) - FP7(M)) - DSS4*FP7(M)
         SM5 = -FP6(M)*(2.*DSS1 - DSS2 - DSS4)
     &    - FP7(M)*(DSS2 - 2.*DSS3 + DSS4)
180     CONTINUE
        RU1 = C4*(DJ1 + DJ2) + DY4*SM1
        RU2 = C4*(DJ2 + DJ3) + DY4*SM2
        RU3 = C4*(DJ1 + DJ4) + DY4*SM3
        RU4 = C4*(DJ3 + DJ4) + DY4*SM4
        RU5 = C4*(DJ1 + DJ2 + DJ3 + DJ4) + DY4*SM5
C
C     THESE ARE THE ELEMENTS OF THE 4 X 4 CONDENSED MATRIX AND SOURCE.
C
        E11 = A11 - A15*(A15/A55)
        G11 = A12 - A15*(A25/A55)
        C11 = A13 - A15*(A35/A55)
        H11 = -A15*(A45/A55)
        E22 = A22 - A25*(A25/A55)
        H22 = -A25*(A35/A55)
        C22 = A24 - A25*(A45/A55)
        E33 = A33 - A35*(A35/A55)
        G33 = A34 - A35*(A45/A55)
        E44 = A44 - A45*(A45/A55)
        RC1 = RU1 - A15*(RU5/A55)
        RC2 = RU2 - A25*(RU5/A55)
        RC3 = RU3 - A35*(RU5/A55)
        RC4 = RU4 - A45*(RU5/A55)
C
C     THIS SECTION LOADS THE ELEMENT MATRICIES INTO THE GLOBAL MATRIX.
C
        S(L1,1) = S(L1,1) + E11
        S(L1,2) = S(L1,2) + G11
        S(L1,NBAN1) = S(L1,NBAN1) + C11
        S(L1,NBAND) = S(L1,NBAND) + H11
        L2 = L1 + 1
        S(L2,1) = S(L2,1) + E22
        S(L2,NDZ) = S(L2,NDZ) + H22
        S(L2,NBAN1) = S(L2,NBAN1) + C22
        L3 = L1 + NDZ
        S(L3,1) = S(L3,1) + E33
        S(L3,2) = S(L3,2) + G33
        L4 = L3 + 1
        S(L4,1) = S(L4,1) + E44
        R(L1) = R(L1) + RC1
        R(L2) = R(L2) + RC2
        R(L3) = R(L3) + RC3
        R(L4) = R(L4) + RC4
        L1 = L1 + 1
150    CONTINUE
140   CONTINUE
C
C     THIS SECTION COMPUTES PRIMARY FIELDS AT RIGHT SIDE SURFACE.
C
      NLYR = NOLA(2)
      NLYR1 = NLYR - 1
      IF(NLYR .EQ. 1)GO TO 190
       DO 200 I = 1,NLYR1
        DDPP(I) = DLR(I,2)
        RL(I) = RLR(I,2)
200    CONTINUE
190   CONTINUE
      RL(NLYR) = RLR(NLYR,2)
C
      CALL COAMP
      IF(IMD .EQ. 1)THEN
       FSE = ZLFLD(0.,1)
       FSH = ZLFLD(0.,2)
       RHA = CABS(FSE/FSH)**2/WU
       PHA = 57.2958*ATAN2(AIMAG(FSE/FSH),REAL(FSE/FSH))
      END IF
C
C     INCORPORATE BOUNDARY CONDITIONS INTO GLOBAL SYSTEM OF EQUATIONS.
C     METHOD 1 OF HUEBNER (1983, P. 52) IS USED, REDUCING UNDERFLOWS.
C     THIS IS NEEDED FOR THE TM MODE ONLY WHERE NEUMANN CONDITIONS
C     ARE INAPPROPRIATE FOR THE AIR-EARTH INTERFACE.
C
      IF(IDX .EQ. 3)THEN
       NJI = NNODE - NDZ + 1
       IYN = 1
       DO 270 IT = 1,NJI,NDZ
        S(IT,1) = (1.,0.)
        IF(IYN .EQ. 1)GO TO 280
         S(IT-1,2) = (0.,0.)
         S(IT-NDZ+1,NDZ) = (0.,0.)
         S(IT-NDZ,NBAN1) = (0.,0.)
         IF(IYN .EQ. 2)GO TO 280
         S(IT-NBAN1,NBAND) = (0.,0.)
280     CONTINUE
        S(IT,2) = (0.,0.)
        S(IT,NDZ) = (0.,0.)
        S(IT,NBAN1) = (0.,0.)
        S(IT,NBAND) = (0.,0.)
        R(IT) = (0.,0.)
        IYN = IYN + 1
270    CONTINUE
      END IF
C      DO 1111 I = 1,NNODE
C      RJ(1,I) = R(I)
C1111  CONTINUE
C
C     SOLVE THE GLOBAL SYSTEM USING GAUSSIAN ELIMINATION
C
      CALL GAUSS(NNODE,NBAND,S,R)
C
      RETURN
      END
C
C **********************************************************************
C **********************************************************************
      SUBROUTINE COAMP()
      INCLUDE 'imp.inc'
C
C========================== INCLUDE FILES ==============================
C
      INCLUDE 'mt2D.inc'
C
C============================ INTERFACE ================================
C
C WRITTEN BY:  P. WANNAMAKER, UURI.
C
C DESCRIPTION: 
C
C     SUBROUTINE TO COMPUTE AMPLITUDE AND REFLECTION COEFFICIENTS FOR
C     ELECTRIC FIELDS PROPAGATING IN LAYERED EARTH.  COEFFICIENTS ARE
C     PASSED TO FUNCTION ZLFLD IN COMMON BLOCK BLKP.
C
C ### RPR CHANGED VARIABLE WE TO WP TO SUPPORT COMMON INCLUDE
C ### RPR CHANGED VARIABLE R TO RZ TO SUPPORT COMMON INCLUDE
C
C============== PARAMETER AND LOCAL VARIABLE DECLARATIONS ==============
C
C     ---- PASS PARAMETERS
C
C     ---- LOCAL VARIABLES
      COMPLEX AEXI,HTN(NL1),A2EX(NL1),CNUM,DEN
      INTEGER I,J,K,NLYR1,KK1
	  COMPLEX*16 DUMMY1, DUMMY2, DUMMY3
C
C=======================================================================
C
      DUMMY1 = -CMPLX(0.,WM)
	  DUMMY2 = CMPLX(1./RA,WP)
      DUMMY3 = SQRT(DUMMY1*DUMMY2)
      K10 = CMPLX(DUMMY3)
      DO 10 I = 1,NLYR
	   DUMMY2 = CMPLX(1./RL(I),WP)
       DUMMY3 = SQRT(DUMMY1*DUMMY2)
       K1(I) = CMPLX(DUMMY3)
10    CONTINUE
      IF(NLYR .GT. 1)THEN
       NLYR1 = NLYR - 1
       H(1) = DDPP(1)
       IF(NLYR1 .EQ. 1)GO TO 30
        DO 40 I = 2,NLYR1
         H(I) = DDPP(I) - DDPP(I-1)
40      CONTINUE
30     CONTINUE
       DO 50 J = 1,NLYR1
        ARG(J) = (0.,1.)*K1(J)*H(J)
        IF(REAL(ARG(J)) .GT. 100.)GO TO 60
         AEX(J) = CEXP(-ARG(J))
         GO TO 70
60      CONTINUE
         AEX(J) = (0.,0.)
70      CONTINUE
        IF(REAL(ARG(J)) .GT. 40.)GO TO 80
         AEXI = 1./AEX(J)
         HTN(J) = (AEXI - AEX(J))/(AEXI + AEX(J))
         GO TO 81
80      CONTINUE
         HTN(J) = (1.,0.)
81      CONTINUE
        A2EX(J) = AEX(J)*AEX(J)
50     CONTINUE
       ZA = WM/K1(NLYR)
       DO 90 I = 1,NLYR1
        J = NLYR - I + 1
        Z1 = WM/K1(J-1)
        RZ(J) = (ZA - Z1)/(ZA + Z1)
        ZA = Z1*(ZA + Z1*HTN(J-1))/(Z1 + ZA*HTN(J-1))
90     CONTINUE
       Z1 = WM/K10
       COA(1) = (2.*ZA/(ZA + Z1))/(1. + RZ(2)*A2EX(1))
       T1 = COA(1)*(1. + RZ(2)*A2EX(1))
       RZ(1) = T1 - 1.
       IF(NLYR1 .EQ. 1)GO TO 100
        DO 110 K = 2,NLYR1
         KK1 = K - 1
         DEN = 1. + RZ(K+1)*A2EX(K)
         CNUM = (1. + RZ(K))*AEX(KK1)
         COA(K) = COA(KK1)*CNUM/DEN
110     CONTINUE
100    CONTINUE
       COA(NLYR) = COA(NLYR1)*AEX(NLYR1)*(1. + RZ(NLYR))
       GO TO 120
      ELSE
       COA(1) = 2.*K10/(K1(1) + K10)
       T1 = COA(1)
       RZ(1) = (K10 - K1(1))/(K10 + K1(1))
      END IF
120   CONTINUE
      RETURN
      END
C
C **********************************************************************
C **********************************************************************
      COMPLEX FUNCTION ZLFLD(Z,IFL)
      INCLUDE 'imp.inc'
C
C========================== INCLUDE FILES ==============================
C
      INCLUDE 'mt2D.inc'
C
C============================ INTERFACE ================================
C
C WRITTEN BY:  P. WANNAMAKER, UURI.
C
C DESCRIPTION: 
C
C     FUNCTION COMPUTES ELECTRIC (IFL .EQ. 1)  OR MAGNETIC (IFL .EQ. 2)
C     FIELD AT A DEPTH Z DUE TO UNIT-AMPLITUDE, DOWNGOING ELECTRIC VECTOR
C     AT THE SURFACE.  REFLECTION AND TRANSMISSION COEFFICIENTS FOR THE
C     LAYERING ARE PASSED FROM COAMP IN COMMON BLOCK BLKP.
C
C============== PARAMETER AND LOCAL VARIABLE DECLARATIONS ==============
C
C     ---- PASS PARAMETERS
C
C     ---- LOCAL VARIABLES
      COMPLEX K1WC,AR,AE,AEI,AR1,AR2
      INTEGER I,IFL,LZP,IM1,LZP1
      REAL ZD1,ZD2,Z
	  COMPLEX*16 DUMMY1, DUMMY2, DUMMY3
C
C  ### RPR CHANGED VARIABLE R TO RZ TO SUPPORT INCLUDED COMMON
C
C=======================================================================
C
      IF(Z .LT. 0.)THEN
       AR = (0.,1.)*K10*Z
       IF(IFL .EQ. 1)  THEN
         DUMMY1 = CMPLX(2.*AR)
    	 DUMMY1 = EXP(DUMMY1)
    	 DUMMY1 = 1.+RZ(1)*DUMMY1
    	 DUMMY2 = CMPLX(-AR)
    	 DUMMY2 = EXP(DUMMY2)
    	 DUMMY3 = DUMMY2*DUMMY1
         ZLFLD = CMPLX(DUMMY3)
       END IF
       IF(IFL .EQ. 2) THEN
         DUMMY1 = CMPLX(2.*AR)
    	 DUMMY1 = EXP(DUMMY1)
    	 DUMMY1 = 1.-RZ(1)*DUMMY1
    	 DUMMY2 = CMPLX(-AR)
    	 DUMMY2 = EXP(DUMMY2)
		 DUMMY3 = (K10/WM)*DUMMY2*DUMMY1
	     ZLFLD = CMPLX(DUMMY3)
       END IF
       GO TO 10
      ELSE
       IF(NLYR .GT. 1)THEN
        LZP = 1
        DO 30 I = 2,NLYR
         IM1 = I - 1
         IF(Z .GT. DDPP(IM1))LZP = I
30      CONTINUE
        K1WC = (K1(LZP)/WM)*COA(LZP)
        IF(LZP .EQ. NLYR)GO TO 40
         IF(REAL(ARG(LZP)) .GT. 40.)GO TO 50
          AR = (0.,1.)*K1(LZP)*(Z-DDPP(LZP))
          AE = CEXP(-AR)
          IF(IFL .EQ. 1)ZLFLD = COA(LZP)*(AE + RZ(LZP+1)/AE)*AEX(LZP)
          IF(IFL .EQ. 2)ZLFLD = K1WC*(AE - RZ(LZP+1)/AE)*AEX(LZP)
         GO TO 10
50       CONTINUE
         IF(LZP .EQ. 1)GO TO 60
          LZP1 = LZP - 1
          ZD1 = Z - DDPP(LZP1)
          ZD2 = Z - 2.*DDPP(LZP) + DDPP(LZP1)
         GO TO 70
60       CONTINUE
          ZD1 = Z
          ZD2 = Z - 2.*DDPP(LZP)
70       CONTINUE
         AR1 = (0.,1.)*K1(LZP)*ZD1
         AE = (0.,0.)
         AEI = (0.,0.)
         IF(REAL(AR1) .GT. 100.)GO TO 80
          AE = CEXP(-AR1)
          AR2 = (0.,1.)*K1(LZP)*ZD2
         IF(REAL(AR2) .LT. -100.)GO TO 80
          AEI = CEXP(AR2)
80       CONTINUE
          IF(IFL .EQ. 1)ZLFLD = COA(LZP)*(AE + RZ(LZP+1)*AEI)
          IF(IFL .EQ. 2)ZLFLD = K1WC*(AE - RZ(LZP+1)*AEI)
         GO TO 10
40      CONTINUE
        AR = (0.,1.)*K1(LZP)*(Z-DDPP(LZP-1))
        IF(REAL(AR) .GT. 100.)GO TO 90
         IF(IFL .EQ. 1)ZLFLD = COA(LZP)*CEXP(-AR)
         IF(IFL .EQ. 2)ZLFLD = K1WC*CEXP(-AR)
        GO TO 10
90      CONTINUE
         ZLFLD = (0.,0.)
        GO TO 10
       ELSE
        AR = (0.,1.)*K1(1)*Z
        IF(REAL(AR) .GT. 100.)GO TO 100
         IF(IFL .EQ. 1)ZLFLD = COA(1)*CEXP(-AR)
         IF(IFL .EQ. 2)ZLFLD = (K1(1)/WM)*COA(1)*CEXP(-AR)
        GO TO 10
100     CONTINUE
         ZLFLD = (0.,0.)
        GO TO 10
       END IF
      END IF
10    CONTINUE
      RETURN
      END
C
C **********************************************************************
C **********************************************************************
      SUBROUTINE GAUSS(NNODE,NBAND,S,R)
      INCLUDE 'imp.inc'
C
C========================== INCLUDE FILES ==============================
C
      INCLUDE 'mt2D.inc'
C
C============================ INTERFACE ================================
C
C WRITTEN BY:  P. WANNAMAKER, UURI.
C
C DESCRIPTION: 
C
C     ROUTINE TO REDUCE FINITE ELEMENT MATRIX.  GAUSSIAN ELIMINATION
C     WITHOUT PIVOTING BUT SPECIALIZED TO BANDED SYMMETRIC MATRICES
C     IS USED.
C
C============== PARAMETER AND LOCAL VARIABLE DECLARATIONS ==============
C
C     ---- PASS PARAMETERS
      INTEGER NNODE,NBAND
      COMPLEX S(IP2,IP3),R(IP2)
C
C     ---- LOCAL VARIABLES
      COMPLEX C
      INTEGER I,J,K,L,M,N
C
C=======================================================================
C
      DO 10 N = 1,NNODE
       DO 20  L = 2,NBAND
        IF(REAL(S(N,L)) - 0.0)30,20,30
30      CONTINUE
        I = N + L - 1
        C = S(N,L)/S(N,1)
        J = 0
        DO 40 K = L,NBAND
         J = J + 1
         S(I,J) = S(I,J) - C*S(N,K)
40      CONTINUE
        S(N,L) = C
20     CONTINUE
10    CONTINUE
      DO 50 N = 1,NNODE
       DO 60 L = 2,NBAND
        IF(REAL(S(N,L)) - 0.0)70,60,70
70      CONTINUE
        I = N + L - 1
        R(I) = R(I) - S(N,L)*R(N)
60     CONTINUE
       R(N) = R(N)/S(N,1)
50    CONTINUE
      DO 80 M = 2,NNODE
       N = NNODE + 1 - M
       DO 90 L = 2,NBAND
        IF(REAL(S(N,L)) - 0.0)100,90,100
100       CONTINUE
        K = N + L - 1
        R(N) = R(N) - S(N,L)*R(K)
90     CONTINUE
80    CONTINUE
      RETURN
      END
C
C **********************************************************************
      SUBROUTINE AUXFLD(I,NVRT,R,RJ)
      INCLUDE 'imp.inc'
C
C========================== INCLUDE FILES ==============================
C
      INCLUDE 'mt2D.inc'
C
C============================ INTERFACE ================================
C
C WRITTEN BY:  P. WANNAMAKER, UURI.
C
C DESCRIPTION: 
C
C     THIS SUBROUTINE CALCULATES THE FIELD NORMAL TO THE SLOPE AND THE
C     FIELD PARALLEL TO THE SLOPE AND TRANSVERSE TO STRIKE.  THIS IS
C     DONE THROUGH THREE-POINT DIFFERENCE APPROXIMATIONS TO MAXWELL'S
C     EQUATIONS IN SUBROUTINE DERIV.  AUXILIARY FIELDS ARE RETURNED
C     TO MAIN IN ARRAYS AFN(I) (NORMAL FIELD) AND AFT(I) (TOTAL TRANS-
C     VERSE FIELD).  FOR THE TE MODE , THE TRANSVERSE H-FIELD IS CALCU-
C     LATED IN A HORIZONTAL DIRECTION.  SIMILAR PROCEDURES ARE CARRIED
C     OUT FOR THE AUXILIARY JACOBIAN QUANTITIES IN AFND(I) AND AFTD(I).
C
C============== PARAMETER AND LOCAL VARIABLE DECLARATIONS ==============
C
C     ---- PASS PARAMETERS
      INTEGER I,NVRT
      COMPLEX R(IP2),RJ(1)
C
C     ---- LOCAL VARIABLES
      COMPLEX DERIV,DV,DH,UC,DVJ(IPI),DHJ(IPI)
      INTEGER NT,NCL,NRW1,NRW,IH,IV,NPL,NPR,NCL1,NCS,KP,NTZ,NPF
      INTEGER IZ1,ISB,IV1,IV2,IV3,IZ2,NHS,MSH,IH1,IH2,IH3,IY1,IY2
      REAL AN,SIGU
C
C=======================================================================
C
C     DEFINE NODAL ROW AND COLUMN LOCATIONS OF FIELDPOINT
C
      NT = NTA(I)
      NCL = I
      NCL1 = NCL - 1
      NRW = NT - NDZ*(I-1) + NTQ - 1
      NRW1 = NRW - 1
C
C     COMPUTE SHIFT IN INDICES OF NODAL FIELDS ALONG STRIKE USED FOR
C     DETERMINING AUXILIARY FIELDS.  NODES USED ARE KEPT WITHIN THE
C     EARTH.  FOR BATHYMETRY, NODES ARE SHIFTED BY TWO INTO SEAWATER.
C
      IH = 0
      IV = -1
      NPL = 0
      NPR = 0
      IF(I .NE. 1)NPL = NPT(NRW,1,NCL1)
      IF(I .NE. NODEY)NPR = NPT(NRW,1,NCL)
      IF(I .EQ. 1 .OR. I .EQ. NODEY)GO TO 30
       CALL BNORM(NCL,NRW,NCL1,NRW1,AN)
       IF(NPL .EQ. 0)IH = 1
       IF(NPR .EQ. 0)IH = -1
       IF(NRW .NE. 1)THEN
        IF(NPL .EQ. ISEA)IH = -1
        IF(NPR .EQ. ISEA)IH = 1
        IF(NPT(NRW1,1,NCL1) .EQ. ISEA .OR. NPT(NRW1,1,NCL).EQ.ISEA)IV=1
       END IF
      GO TO 40
30    CONTINUE
       AN = 0.
40    CONTINUE
      IH1 = NT + NDZ*(IH - 1)
      IH2 = NT + NDZ*IH
      IH3 = NT + NDZ*(IH + 1)
      IV1 = NT - IV - 1
      IV2 = NT - IV
      IV3 = NT - IV + 1
      IY1 = NCL1 + IH
      IY2 = NCL + IH
      IZ1 = NRW1 - IV
      IZ2 = NRW - IV
      ISB = IZ1 + (1 + IV)/2
C
C     COMPUTE VERTICAL AND HORIZONTAL FIELD DERIVATIVES.  AVOID JUMPS
C     IN SIG OR SIGH ONE NODE AWAY WITH FIRST-ORDER DIFFERENCES.
C
      NCS = 0
      IF(I .EQ. NODEY)NCS = 1
      IF(NVRT .EQ. 0)THEN
       DV = -IV*(R(IV2) - R(IV1))/DLZ(ISB)
      END IF
      IF(IZ1 .LT. 1)GO TO 50
       IF(NPT(ISB,2+2*NCS,NCL-NCS) .NE. NPT(ISB-IV,2+2*NCS,NCL-NCS)
     & .OR. NPT(ISB,3,1) .NE. NPT(ISB-IV,1,1))GO TO 50
       IF(NVRT .EQ. 0)THEN
        DV = DERIV(R(IV1),R(IV2),R(IV3),DLZ(IZ1),DLZ(IZ2),IV)
		end if
50    CONTINUE
      DH = (0.,0.)
      DO 48 KP = 1,NPRM
       DHJ(KP) = (0.,0.)
48    CONTINUE
      IF(I .EQ. 1 .OR. I .EQ. NODEY)GO TO 52
       IF(IH .NE. 0)THEN
        NHS = 1
        MSH = 1
        IF(IH .EQ. 1)NHS = 0
        IF(IH .EQ. 1)MSH = -1
        IF(NPT(ISB,2+IV,NCL-NHS) .NE. NPT(ISB,2+IV,NCL+1-3*NHS))THEN
         NTZ = NT + NDZ*(1-2*NHS)
         IF(NVRT .EQ. 0)THEN
          DH = MSH*(R(NT) - R(NTZ))/DLY(NCL-NHS)
         END IF
         GO TO 52
        END IF
       END IF
       IF(NVRT .EQ. 0)THEN
        DH = DERIV(R(IH1),R(IH2),R(IH3),DLY(IY1),DLY(IY2),IH)
       END IF
52    CONTINUE
C
C     TE MODE AUXILIARY FIELDS AND DERIVATIVES
C
      IF(IDX .EQ. 2)THEN
       AN = 0.
       UC = CMPLX(0.,WU)
       IF(NVRT .EQ. 0)THEN
        AFN(I) = -(FPB(NRW) - DV/UC)*SIN(AN) + (DH/UC)*COS(AN)
        AFT(I) = (FPB(NRW) - DV/UC)*COS(AN) + (DH/UC)*SIN(AN)
       END IF
C
C      TM MODE AUXILIARY FIELDS AND DERIVATIVES
C
       ELSE
        SIGU = SIG(ISB,2+2*NCS,NCL-NCS)
        UC = CMPLX(SIGU,WE)
C       AN = 0.
        IF(NVRT .EQ. 0)THEN
         AFN(I) = (-(SIGH(ISB)*FPB(NRW) + DV)*SIN(AN) - DH*COS(AN))/UC
         AFT(I) = ((SIGH(ISB)*FPB(NRW) + DV)*COS(AN) - DH*SIN(AN))/UC
        END IF
       END IF
      RETURN
      END
C
C **********************************************************************
C **********************************************************************
      SUBROUTINE BNORM(NCL,NRW,NCL1,NRW1,AN)
      INCLUDE 'imp.inc'
C
C========================== INCLUDE FILES ==============================
C
      INCLUDE 'mt2D.inc'
C
C============================ INTERFACE ================================
C
C WRITTEN BY:  P. WANNAMAKER, UURI.
C
C DESCRIPTION: 
C
C     SUBROUTINE CALCULATES DOWNWARD BISECTOR FROM FIELD NODE AT BREAK
C     IN TOPOGRAPHY.  AR, AL CLOCKWISE FROM UP, AN CLOCKWISE FROM DOWN.
C
C============== PARAMETER AND LOCAL VARIABLE DECLARATIONS ==============
C
C     ---- PASS PARAMETERS
      INTEGER NCL,NRW,NCL1,NRW1
      REAL AN
C
C     ---- LOCAL VARIABLES
      COMPLEX TAL,TAR
      INTEGER IR,I
      DIMENSION IR(8)
      REAL SM,AL,AR
C
C=======================================================================
C
C     IDENTIFY RESISTIVITY CODES OF TRIANGULAR ELEMENTS SURROUNDING
C     FIELD NODE.  UPPERMOST LEFT OCTANT IS IR(1), RIGHT IS IR(8).
C
      IR(3) = NPT(NRW,1,NCL1)
      IR(4) = NPT(NRW,4,NCL1)
      IR(5) = NPT(NRW,2,NCL)
      IR(6) = NPT(NRW,1,NCL)
      IF(NRW .EQ. 1)GO TO 10
       IR(1) = NPT(NRW1,4,NCL1)
       IR(2) = NPT(NRW1,3,NCL1)
       IR(7) = NPT(NRW1,3,NCL)
       IR(8) = NPT(NRW1,2,NCL)
      GO TO 20
10    CONTINUE
       IR(1) = 0
       IR(2) = 0
       IR(7) = 0
       IR(8) = 0
20    CONTINUE
      DO 40 I = 1,8
       IF(IR(I) .EQ. ISEA)IR(I) = 0
40    CONTINUE
C
C     DETERMINE SLOPE TO LEFT OF FIELD NODE
C
      SM = 1.0E-05
      TAL = 100.*CMPLX(1.,-SM)
      IF(NRW .EQ. 1)GO TO 50
      TAL = DLZ(NRW1)*CMPLX(1.,-SM)
      IF(IR(1) .EQ. 0 .AND. IR(2) .NE. 0)
     1  TAL = CMPLX(DLZ(NRW1),-DLY(NCL1))
50    CONTINUE
      IF(IR(2) .EQ. 0 .AND. IR(3) .NE. 0)
     1  TAL = DLY(NCL1)*CMPLX(SM,-1.)
      IF(IR(3) .EQ. 0 .AND. IR(4) .NE. 0)
     1  TAL = CMPLX(-DLZ(NRW),-DLY(NCL1))
      IF(IR(4) .EQ. 0 .AND. IR(5) .NE. 0)
     1  TAL = DLZ(NRW)*CMPLX(-1.,-SM)
C
C     DETERMINE SLOPE TO RIGHT OF FIELD NODE
C
      TAR = 100.*CMPLX(1.,SM)
      IF(NRW .EQ. 1)GO TO 60
      TAR = DLZ(NRW1)*CMPLX(1.,SM)
      IF(IR(8) .EQ. 0 .AND. IR(7) .NE. 0)
     1  TAR = CMPLX(DLZ(NRW1),DLY(NCL))
60    CONTINUE
      IF(IR(7) .EQ. 0 .AND. IR(6) .NE. 0)
     1  TAR = DLY(NCL)*CMPLX(SM,1.)
      IF(IR(6) .EQ. 0 .AND. IR(5) .NE. 0)
     1  TAR = CMPLX(-DLZ(NRW),DLY(NCL))
      IF(IR(5) .EQ. 0 .AND. IR(4) .NE. 0)
     1  TAR = DLZ(NRW)*CMPLX(-1.,SM)
      AL = ATAN2(AIMAG(TAL),REAL(TAL))
      AR = ATAN2(AIMAG(TAR),REAL(TAR))
      AN = (AL + AR)/2.
      RETURN
      END
C
C **********************************************************************
C **********************************************************************
      COMPLEX FUNCTION DERIV(R1,R2,R3,D1,D2,IF)
      INCLUDE 'imp.inc'
C
C============================ INTERFACE ================================
C
C WRITTEN BY:  P. WANNAMAKER, UURI.
C
C DESCRIPTION: 
C
C     FUNCTION ESTIMATES DERIVATIVE AT LEFT END, MIDDLE OR RIGHT END
C     OF A COMPLEX THREE-POINT SEQUENCE BY FITTING THE PARABOLA
C     R = A + B*D + C*D**2.  PARABOLA COEFFICIENTS ARE CALCULATED
C     FROM AMPLITUDES AND SEPARATIONS OF THE THREE POINTS.
C
C============== PARAMETER AND LOCAL VARIABLE DECLARATIONS ==============
C
C     ---- PASS PARAMETERS
      COMPLEX R1,R2,R3
      REAL D1,D2
      INTEGER IF
C
C     ---- LOCAL VARIABLES
      COMPLEX B,C
      REAL DEN
C
C=======================================================================
C
      DEN = D1*D2*D2 + D1*D1*D2
C     A = R2
      B = ((R3-R2)*D1*D1 - (R1-R2)*D2*D2)/DEN
      C = ((R3-R2)*D1 + (R1-R2)*D2)/DEN
      IF(IF)10,20,30
C
C     DERIVATIVE AT LEFT POINT
10    CONTINUE
       DERIV = B - 2.*C*D1
      GO TO 40
C
C     DERIVATIVE AT MIDDLE POINT
20    CONTINUE
       DERIV = B
      GO TO 40
C
C     DERIVATIVE AT RIGHT POINT
30    CONTINUE
       DERIV = B + 2.*C*D2
40    CONTINUE
      RETURN
      END

c-----------------------------------------------------------------------
      subroutine addair(lidx, mair)
c-----------------------------------------------------------------------
c
c OCCAM2DMT 1.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 1.01, 13 Jan 1993
c 
c addair adds or removes an air layer of 10 nodes 
c from the FE mesh (nodez, dlz(), npt() in common blocks)
c if a TE mode calculation is being done.  Much of this code has been
c abstracted from Phil Wannamaker's PW2DI
c
c on input:
c   lidx = flag for direction and action (local version of one in PW2DI): 
c           = 1 or 2 causes addition of air layer
c           = 3 sets mair = 0 and does nothing else
c           = otherwise removes air layer  
c   mair = unchanged from any previous call to addair
c
c on ouput:
c   mair = number of air layers, 0 for TM, 10 for TE
c   /blki/nodez, /blkd/dlz() and /blke/npt() are modified for TE calculation
c
c calls:
c   no other routines
c
c-----------------------------------------------------------------------
c includes:
      include 'imp.inc'
      include 'mt2D.inc'
c  uses dlz(), yy()
c  uses npt()
c  uses nodez, nodey
c
c arguments 
      integer lidx, mair
c local variables
      REAL DZF,XA,XB,FXA,FXM,XM
      integer i,j,l,m

c-----------------------------------------------------------------------
      if (lidx .eq. 3) then
        MAIR = 0
        return
      else if (lidx.eq.1  .or. lidx.eq.2) then
C
        MAIR = 10
        NODEZ = NODEZ + MAIR
c
c make room in dlz() and npt() for air layers (push things down mair elements)
        do 10 j = nodez-mair-1, 1, -1
          dlz(j+mair) = dlz(j)
          do 10 l = 1, nodey-1
            do 10 m = 1,4
10            npt(j+mair,m,l) = npt(j,m,l)
C
c compute the thickness of the elements in the air layer
C TOTAL AIR THICKNESS IS ONE-HALF MESH WIDTH, I.E., 
C YY(NODEY) ABOVE.  ELEMENT ROW JUST ABOVE SURFACE 
C EQUALS THAT JUST BELOW AND ROW THICKNESSES INCREASE
C EXPONENTIALLY UPWARD.  IN ALL, THERE ARE MAIR ELEMENT 
C ROWS OF AIR.
        DZF = YY(NODEY)/DLZ(MAIR + 1)
        XA = 1.
        XB = DZF
        DO 40 I = 1,50
          XM = (XA + XB)/2.
          IF((XM-XA)/XM .LT. 1.0E-04)GO TO 41
          FXA = 1.
          FXM = 1.
          DO 42 J = 1,MAIR-1
            FXA = 1. + XA*FXA
            FXM = 1. + XM*FXM
 42       CONTINUE
          FXA = FXA - DZF
          FXM = FXM - DZF
C
          IF(SIGN(1.,FXA)*SIGN(1.,FXM) .LE. 0)THEN
            XB = XM
          ELSE
            XA = XM
          END IF
 41       CONTINUE
 40     CONTINUE
C
        DLZ(MAIR) = DLZ(MAIR + 1)
        DO 43 I = 1,MAIR-1
          J = MAIR - I
          DLZ(J) = XM*DLZ(J+1)
43      CONTINUE
C
c fill the air layer lookup table with 0 (air)
        DO 50 j = 1,MAIR
          DO 50 l = 1,NODEY-1
            DO 50 m = 1,4
50            NPT(j,m,l) = 0
c
      else
c remove air layers in dlz() and npt() (push things up mair elements)
        if (mair .ne. 0) then
          do 60 j = 1, nodez-mair-1
            dlz(j) = dlz(j+mair)
            do 60 l = 1, nodey-1
              do 60 m = 1,4
60              npt(j,m,l) = npt(j+mair,m,l)
          nodez = nodez - mair
        end if
      end if
      return
      end
      