C-----------------------------------------------------------------------
	program PlotModel
C-----------------------------------------------------------------------
C OCCAM 2.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Program Revision 2.00, 1 Feb 93
c revisions 31 July to include poly file format on unit io2
c revisions 26 Feb 97 to take it out again
c revisions 28th July 98 (added depth limit)
c 
c Reads in standard model and data files and outputs a model in a 
c format suitable for plotting or printing.  Outmod() is the routine
c that does all the work, and so this is where all the problem-dependent
c stuff goes.   FindResponse does the same thing, but since it also
c does the forward computation, it is slower and more memory-hungry.
c
	include 'occamdim.inc'
c _______________________start occam includes___________________________	
c Commons Revision 2.00, 13 May 1992
	real d, sd, dp
	integer nd, ndp
	common /data/ nd, d(ndd), sd(ndd), dp(ndd,4), ndp
	real pm, dm
	integer np
	common /model/ np, pm(npp), dm(ndd)
	real ptp, wjtwj, wjtwd, pwk1, pwk2, pwk3, frac, premod, prewts
	integer idebug, nfor, iout
	common /result/ ptp(npp,npp), wjtwj(npp,npp), wjtwd(npp),
     *  pwk1(npp), pwk2(npp), pwk3(npp), frac, nfor, iout, idebug,
     *  prewts(npp),premod(npp)
C ___________________end occam includes_________________________________
	character*80 itform,descr,modfil,datfil,iterfl
	real dwk1(ndd)

c calls
c   itrin, inputd, inputm, outmod
c
C-----------------------------------------------------------------------
c read the iteration file
	write(*,*) 'Enter iteration or startup file name'
	read (*,'(a)') iterfl
	open(12,file=iterfl)
      call itrin (itform,descr,modfil,datfil,maxitr,tolreq,iruf,
     *    nit,pmu,rlast,tobt,ifftol,12)
c input the data
c dont input the data till we need station locations
	call inputd(datfil)
c input the fixed model characteristics
	call inputm(modfil)
c output the model
	open(12, file='model.ps')
	call outmod(12)
	close(12)
c
	stop
c
198	call filerr (' Iteration file ended prematurely', iof2)
199	call filerr (' Error reading iteration file', iof2)
c	
	end
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:
      implicit none
      include 'occamdim.inc'
c   uses npp
c _______________________start occam includes___________________________	
c Commons Revision 2.00, 13 May 1992
	real d, sd, dp
	integer nd, ndp
	common /data/ nd, d(ndd), sd(ndd), dp(ndd,4), ndp
	real pm, dm
	integer np
	common /model/ np, pm(npp), dm(ndd)
	real ptp, wjtwj, wjtwd, pwk1, pwk2, pwk3, frac, premod, prewts
	integer idebug, nfor, iout
	common /result/ ptp(npp,npp), wjtwj(npp,npp), wjtwd(npp),
     *  pwk1(npp), pwk2(npp), pwk3(npp), frac, nfor, iout, idebug,
     *  prewts(npp),premod(npp)
C ___________________end occam includes_________________________________
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:
      implicit none
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:
      implicit none
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 + 0.001)
	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:
      implicit none
c input arguments:
      character*(*) mssg
      integer io1
c
C-----------------------------------------------------------------------
      write(*,*) mssg
      if (io1 .gt. 0) close (io1)
      stop
      end
C________________________________________________________________________
	SUBROUTINE OUTMOD(IO)
c
c OCCAM2DMT 1.0: Steven Constable IGPP/SIO La Jolla CA 92093-0225
c Subroutine Revision 1.01, 14 July 1992
c 
C OUTMOD WRITES A 2D MT MODEL TO DISK as a colour postscript file
c (this does not work with fixed structure, but could be made to
c do so with a little bit of cunning.)
c
c
c _______________________start MT includes______________________________	
c Commons Revision 1.04, 13 Jan 93
      include 'occamdim.inc'
      include 'mt2Ddim.inc'
      integer nlay, nrcol, nrfix, icorn, nbrick, nexep
      integer ijexep, iprewt
      real thick, width, y, expen
      common /rgrid/ nlay, nrcol(nl), thick(npp), width(npp), 
     *  nrfix, y(npp+36), icorn(npp,4), 
     *  nbrick, nexep, ijexep(maxexp,2), expen(maxexp), iprewt(npp)
      integer nrc, nfre, nrl, ndata
      real sitloc, boffs, frek
      common /datinfo/ nrc, nrl(nr), ndata, nfre, boffs,  
     *  frek(nf), sitloc(nr)
      character*30 sitnam
      common /sitcom/ sitnam(nr)
      integer isx1, isx3, iscon, nfss, isfx, nsb
      real shift
      common /statics/ nsb, iscon, nfss, isx1(2*nr), isx3(2*nr), 
     *  isfx(2*nr), shift(2*nr)
      REAL YY(IPY),ZZ(IPZ),DLY(IPY),DLZ(IPZ)
      COMMON /BLKD/ YY,ZZ,DLY,DLZ
      INTEGER NPT(IPZ,4,IPY)
      REAL SIG(IPZ,4,IPY)
      COMPLEX SIGH(IPZ)
      COMMON /BLKE/ NPT,SIG,SIGH
      COMPLEX FPA(IPZ),FPB(IPZ),AFN(IPY),AFT(IPY)
      REAL PI,WU,WE
      COMMON /BLKF/ FPA,FPB,AFN,AFT,PI,WU,WE
      INTEGER IDX,IMD,NODEY,NODEZ,NTQ,NDZ,NDZ1,NTAO(IPY),NTA(IPY)
      INTEGER NZZB(2),NOLA(2),NPRM,LAPRM(IPI),ISEA
      COMMON /BLKI/ IDX,IMD,NODEY,NODEZ
      COMMON /BLKI/ NTQ,NDZ,NDZ1,NTAO,NTA
      COMMON /BLKI/ NZZB,NOLA,NPRM,LAPRM,ISEA
      COMPLEX KZY(NR,NF),AFND(IPI),AFTD(IPI)
      REAL RHTE(NR,NF),PHTE(NR,NF),RHTM(NR,NF),PHTM(NR,NF)
      REAL KZRD(IPI,NR,NF),KZID(IPI,NR,NF)
      REAL RTED(IPI,NR,NF),PTED(IPI,NR,NF)
      REAL RTMD(IPI,NR,NF),PTMD(IPI,NR,NF)
      COMMON /BLKJ/ RHTE,PHTE,RHTM,PHTM
      COMMON /BLKJ/ KZY,AFND,AFTD,KZRD,KZID
      COMMON /BLKJ/ RTED,PTED,RTMD,PTMD
      REAL DLR(NL1,2),RLR(NL,2),DDPP(NL1),RL(NL),RA,W,WM,WP
      INTEGER NLYR
      COMMON /BLKL/ DLR,RLR,DDPP,RL,RA,NLYR,W,WM,WP
      COMPLEX ARG(NL1),AEX(NL1),K10,K1(NL),RZ(NL),COA(NL),Z1,ZA,T1
      REAL H(NL1)
      COMMON /BLKP/ ARG,AEX,K10,K1,RZ,COA,Z1,ZA,T1,H
C ___________________end MT includes____________________________________
c _______________________start occam includes___________________________	
c Commons Revision 2.00, 13 May 1992
	real d, sd, dp
	integer nd, ndp
	common /data/ nd, d(ndd), sd(ndd), dp(ndd,4), ndp
	real pm, dm
	integer np
	common /model/ np, pm(npp), dm(ndd)
	real ptp, wjtwj, wjtwd, pwk1, pwk2, pwk3, frac, premod, prewts
	integer idebug, nfor, iout
	common /result/ ptp(npp,npp), wjtwj(npp,npp), wjtwd(npp),
     *  pwk1(npp), pwk2(npp), pwk3(npp), frac, nfor, iout, idebug,
     *  prewts(npp),premod(npp)
C ___________________end occam includes_________________________________
C
c
	real colors(3,0:255)
	real xsite(200), ysite, ddd(300), ddx(300), ddy(300)

      write(*,*) ' Enter color table;'
	write(*,*) ' 1: 256 linear greyscale'
	write(*,*) ' 2: 256 red to blue'
	write(*,*) ' 3: 256 rainbow'
	write(*,*) ' 4: 256 truncated greyscale'
	write(*,*) ' -ve number reverses any of these'
	read (*,*) icol
	if (abs(icol) .eq. 1) call palet1(colors) 
	if (abs(icol) .eq. 2) call palet2(colors) 
	if (abs(icol) .eq. 3) call palet3(colors) 
	if (abs(icol) .eq. 4) call palet4(colors) 
	if (icol .lt. 0) call revct(colors)  

	write(*,*) ' Draw mesh (1) or not (0)'
	read (*,*) ifmesh

	write(*,*) ' Enter log resistivity limits (low, high) or default (0,0)'
	read (*,*) fmin, fmax
c set plot limits to go from 1.0 to 10**4 Ohm-m
	if (fmin .ge. fmax) then
	  fmax = 4.0
	  fmin = 0.0
	end if
	range = fmax-fmin
	
	write(*,*) ' Enter vertical axis scaling; 0=log, 1=linear'
	read (*,*) ivax
	if (ivax .eq. 1) then
	  write(*,*) ' Enter vertical exaggeration'
	  read (*,*) exag
	end if
	write(*,*) ' Enter maximum depth to plot'
	read (*,*) dmax

c set up locations of the nodes:
c first get the boundaries of the actual model (less the side and bottom layers)
      ileft = icorn(1,4) 
      iright = icorn(np,2) 
	ibottom = icorn(np,1)
	itop = icorn(1,3)
c then accumualate the horizontal positions and depths
	yy(1) = 0.0
	zz(1) = 0.0
	do 1 i = 2, nodey
c make the left and right layers a reasonable size by re-setting the fe widths
	  if (i .le. ileft) then
	    yy(i) = yy(i-1) + dly(ileft)
	  else if (i .gt. iright) then
	    yy(i) = yy(i-1) + dly(iright-1)
	  else
 	    yy(i) = yy(i-1) + dly(i-1)
	  end if
1     continue
c set the origin using boffs
      y0 = yy(ileft)
	do 2 i = 1, nodey
2	  yy(i) = yy(i) - y0 + boffs
c set up the vertical dimension using the same trick to get the bottom layer
      maxz = nodez
	do 3 i = 2, nodez
	  if (i .gt. ibottom) then
	    zz(i) = zz(i-1) + dlz(ibottom-1)
	  else
	    zz(i) = zz(i-1) + dlz(i-1)
	  end if
	  if (zz(i) .gt. dmax) then
	    maxz = i
	    goto 5
	  end if
3     continue
5      if (ivax .eq. 0) then
	  do 4 i = 2, nodez
4           zz(i) = alog10(zz(i))
      end if
	    
c set scale up to give 6in by 6in plot (i.e. total nodes = 6in)
	xscale = 432./(yy(nodey) - yy(1))
	x0 = yy(1)
      z0 = zz(maxz)
c      z0 = zz(nodez)
	top = 0.7*thick(1)
c	if (ivax .eq. 0) top = 0.0
	yscale = 432./(zz(maxz))
c	yscale = 432./(zz(nodez))
      if (ivax .eq. 0) then
	  yscale = 1.4*yscale
      end if
	write(*,*) 'yscale, xscale :',yscale, xscale
	if (ivax .eq. 1) then
	  yscale = exag*xscale
	end if

c write postscript header	
      write (io,'(a4)') '%!PS'
      write (io,*) '0 setlinecap'
      write (io,*) '0 setlinejoin'
      write (io,*) '30 0  translate'
	write (io,*) '/Times-Roman findfont 12 scalefont setfont'

      write (io,*) 
     * '/box {gsave setrgbcolor moveto currentpoint translate'
      write (io,*) '      0 0      moveto'
      write (io,*) '      dup 0 exch lineto'
      write (io,*) '      dup 3 -1 roll dup 4 1 roll exch lineto'
      write (io,*) '      exch dup 0 lineto'
      write (io,*) '      0 0 lineto pop pop'
      if (ifmesh) then
      write (io,*) 
     * '      closepath gsave fill grestore 0 setgray stroke'
      else
      write (io,*) 
     * '      closepath gsave fill grestore 0 setgray'
      end if
      write (io,*) '      grestore} def'
      write (io,*) 
     * '/site {moveto -3 4 rlineto 6 0 rlineto -3 -4 rlineto stroke} def'
     
      depth = 0.0
      i = 0
	
c main loop: for every layer
      do 110 il = 1,nlay
c find the depth to the bottom of the layer:
c if this is the bottom layer reduce the thickness from infinity
	  if (il .eq. nlay) then
	    depth = depth + (nodez - ibottom)*dlz(ibottom-1)
c	    depth = depth + (maxz - ibottom)*dlz(ibottom-1)
	  else
c else just take thickness as given
	    depth = depth + thick(i+1)
	  end if
c store the depth for plotting depth labels
	  ddd(il) = depth
c internal loop: for every column
        do 100 j = 1, nrcol(il)
c increment the brick counter
          i = i + 1
c find color
	    indx = 2 + int(254.*(pm(i) - fmin)/range)
	    if (indx .lt. 2) indx = 2
	    if (indx .gt. 255) indx = 255
c find box origin coords in pixels
	    xxx = xscale*(yy(icorn(i,2))-x0) + 77.
	    yyy = yscale*(z0-zz(icorn(i,3))) + 150.
c if this is the leftmost col make the width less than infinity
          if (j .eq. 1) then
	      width(i) = (ileft - 1)*dly(ileft)
c also record coords for plotting depth labels
	      ddx(il) = xxx - 22.
		ddy(il) = yyy - 4.
	    end if
c if this is the rightmost col ditto
          if (j .eq. nrcol(il)) then
	      width(i) = (nodey - iright)*dly(iright-1)
	    end if
c if this is the bottom layer make the thickness less than infinity
          if (il .eq. nlay) then
	      thick(i) = (nodez - ibottom)*dlz(ibottom-1)
	    end if
c find width in pixels	
	    dx = xscale*width(i)
c find depth in  pixels, adjusting for a log scale if used
	    if (ivax .eq. 0) then
	      if (il .eq. 1) then 
	        dy = yscale*alog10(depth/(depth-thick(i)+top))
		  yref = yyy + dy
		else
	        dy = yscale*alog10(depth/(depth-thick(i)))
		end if
	    else
	      dy = yscale*thick(i)
	    end if
c if this is the uppermost layer and leftmost col find the site 
c coords for plotting
	    if (il = 1 .and. j .eq. 1) then
	      ysite = yyy + dy 
	      do 88 is = 1 , nrc
88		  xsite(is) = xscale*(sitloc(is)-x0) + 77.
          end if

c write out box description to postscript file
     	  write(io,'(2f6.1,2f7.1,3f6.3,'' box'')')  
     *    dx, dy, xxx, yyy, colors(1,indx),
     *    colors(2,indx), colors(3,indx)
     
100     continue
      if (depth .gt. dmax) goto 150
110   continue
  
c color bar:
150	write(io, *) ' % Color bar:'
      dx = 20
	dy = 20
	xxx = 70
	yyy = 80
	do 200 i = 0, 20
	  xxx = xxx + dx
	  val = fmin + float(i)*range/20.
	  indx = 2 + int(254.*(val - fmin)/range)
	  if (indx .lt. 2) indx = 2
	  if (indx .gt. 255) indx = 255
200	  write(io,'(2f6.1,2f7.1,3f6.3,'' box'')')  
     *    dx, dy, xxx, yyy, colors(1,indx),
     *    colors(2,indx), colors(3,indx)

       write(io,*) '225 65 moveto (Scale, log resistivity) show'
	 yyy = 106.
	 scale = 400./range
	 do 202 i = int(fmin), int(fmax)
	   xxx = 90. + (float(i) - fmin)*scale
	   write(io,'(2f7.1, a, i2, a)') xxx, yyy, ' moveto (', i, ') show'
202 	 continue
       if (int(fmin) .ne. fmin) then
	   xxx = 90. 
	   write(io,'(2f7.1, a, f4.1, a)') xxx, yyy, ' moveto (', fmin, ') show'
       end if
       if (int(fmax) .ne. fmax) then
	   xxx = 90. + (fmax - fmin)*scale
	   write(io,'(2f7.1, a, f4.1, a)') xxx, yyy, ' moveto (', fmax, ') show'
       end if

c site locations
	do 188 is = 1, nrc
188     write(io,*)  xsite(is), ysite,  'site'

c depths
	do 199 i = 1, nlay
	  if (ddd(i) .lt. 1000.) then
          write(io,'(2f7.1, a, i3, a)') 
     *      ddx(i)-12., ddy(i), ' moveto (', int(ddd(i)+0.5), ' m) show'
        else if (ddd(i) .lt. 10000.) then
          write(io,'(2f7.1, a, f3.1, a)') 
     *      ddx(i)-18., ddy(i), ' moveto (', ddd(i)/1000., ' km) show'
        else 
          write(io,'(2f7.1, a, i3, a)') 
     *      ddx(i)-18., ddy(i), ' moveto (', int(ddd(i)/1000. + .5), ' km) show'
        end if
	if (ddd(i) .gt. dmax) goto 250
199	continue	  

c write out the left and right coordinates of the mesh
250   write(io,'(2f7.1, a, f6.0, a)') xscale*(yy(ileft)-x0) + 57., ddy(i)-12, ' moveto (', yy(ileft), ') show'
      write(io,*)  xscale*(yy(ileft)-x0) + 77., ddy(i),  'site'
	write(io,'(2f7.1, a, f6.0, a)') xscale*(yy(iright)-x0) + 57., ddy(i)-12, ' moveto (', yy(iright), ') show'
      write(io,*)  xscale*(yy(iright)-x0) + 77., ddy(i),  'site'

300	write(io,*) ' showpage'

      return
	end
		
c ***************************************************************output
	subroutine revct(ctable)
c reverses color table
	
	dimension ctable(3,0:255), ct2(3,0:255)
	
	do 10 i = 2,255
	  ct2(1,i) = ctable(1, 257-i)
	  ct2(2,i) = ctable(2, 257-i)
10	  ct2(3,i) = ctable(3, 257-i)
	  
	do 20 i = 2,255
	  ctable(1,i) = ct2(1, i)
	  ctable(2,i) = ct2(2, i)
20	  ctable(3,i) = ct2(3, i)
	  
	return
	end
c ***************************************************************revct

	subroutine palet2(ctable)
c red to blue with no cyan (red-yellow, yellow-green, green-blue)
	
	dimension ctable(3,0:255)
		  
	do 11 i = 2,86
	  ctable(1,i) = 1.0
	  ctable(2,i) = (float(i)-2.)/85.
	  ctable(3,i) = 0.0
11	continue

	do 12 i = 87,170
	  ctable(1,i) = 1.0 - (float(i)-86.)/84.
	  ctable(2,i) = 1.0
	  ctable(3,i) = 0.0
12	continue
	  
	do 13 i = 171,255
	  ctable(1,i) = 0.0
	  ctable(2,i) = 1.0 - (float(i)-170.)/85.
	  ctable(3,i) = (float(i)-170.)/85.
13	continue
	
	return
	end	
c ***************************************************************palet1
	subroutine palet3(ctable)
c rainbow with cyan  (red-yellow, yellow-green, green-cyan, cyan-blue)
c but with green shortened
	
	dimension ctable(3,0:255)
		
	do 11 i = 2,68
	  ctable(1,i) = 1.0
	  ctable(2,i) = (float(i)-2.)/67.
	  ctable(3,i) = 0.0
11	continue
	  
	do 12 i = 69,128
	  ctable(1,i) = 1.0 - (float(i)-68.)/60.
	  ctable(2,i) = 1.0
	  ctable(3,i) = 0.0
12	continue
	  
	do 13 i = 129,188
	  ctable(1,i) = 0.0
	  ctable(2,i) = 1.0
	  ctable(3,i) = (float(i)-128.)/60.
13	continue
	  
	do 14 i = 189,255
	  ctable(1,i) = 0.0
	  ctable(2,i) = 1.0 - (float(i)-188.)/67.
	  ctable(3,i) =  1.0
14	continue
	  
	return
	end
c ***************************************************************palet2
	subroutine palet4(ctable)
c truncated 255 grey-scale
	
	dimension ctable(3,0:255)

	const = 0.7/255.
	do 16 i = 2,255
	  grey = 0.3 + float(i)*const
	  ctable(1,i) = grey
	  ctable(2,i) = grey
	  ctable(3,i) = grey
16	continue
	
	return
	end
c ***************************************************************palet3
	subroutine palet1(ctable)
c linear 255 grey-scale
	
	dimension ctable(3,0:255)

	const = 1.0/255.
	do 16 i = 2,255
	  grey = float(i)*const
	  ctable(1,i) = grey
	  ctable(2,i) = grey
	  ctable(3,i) = grey
16	continue
	
	return
	end
c ***************************************************************palet0
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:
      implicit none
c _______________________start MT includes______________________________	
c Commons Revision 1.04, 13 Jan 93
      include 'occamdim.inc'
      include 'mt2Ddim.inc'
      integer nlay, nrcol, nrfix, icorn, nbrick, nexep
      integer ijexep, iprewt
      real thick, width, y, expen
      common /rgrid/ nlay, nrcol(nl), thick(npp), width(npp), 
     *  nrfix, y(npp+36), icorn(npp,4), 
     *  nbrick, nexep, ijexep(maxexp,2), expen(maxexp), iprewt(npp)
      integer nrc, nfre, nrl, ndata
      real sitloc, boffs, frek
      common /datinfo/ nrc, nrl(nr), ndata, nfre, boffs,  
     *  frek(nf), sitloc(nr)
      character*30 sitnam
      common /sitcom/ sitnam(nr)
      integer isx1, isx3, iscon, nfss, isfx, nsb
      real shift
      common /statics/ nsb, iscon, nfss, isx1(2*nr), isx3(2*nr), 
     *  isfx(2*nr), shift(2*nr)
      REAL YY(IPY),ZZ(IPZ),DLY(IPY),DLZ(IPZ)
      COMMON /BLKD/ YY,ZZ,DLY,DLZ
      INTEGER NPT(IPZ,4,IPY)
      REAL SIG(IPZ,4,IPY)
      COMPLEX SIGH(IPZ)
      COMMON /BLKE/ NPT,SIG,SIGH
      COMPLEX FPA(IPZ),FPB(IPZ),AFN(IPY),AFT(IPY)
      REAL PI,WU,WE
      COMMON /BLKF/ FPA,FPB,AFN,AFT,PI,WU,WE
      INTEGER IDX,IMD,NODEY,NODEZ,NTQ,NDZ,NDZ1,NTAO(IPY),NTA(IPY)
      INTEGER NZZB(2),NOLA(2),NPRM,LAPRM(IPI),ISEA
      COMMON /BLKI/ IDX,IMD,NODEY,NODEZ
      COMMON /BLKI/ NTQ,NDZ,NDZ1,NTAO,NTA
      COMMON /BLKI/ NZZB,NOLA,NPRM,LAPRM,ISEA
      COMPLEX KZY(NR,NF),AFND(IPI),AFTD(IPI)
      REAL RHTE(NR,NF),PHTE(NR,NF),RHTM(NR,NF),PHTM(NR,NF)
      REAL KZRD(IPI,NR,NF),KZID(IPI,NR,NF)
      REAL RTED(IPI,NR,NF),PTED(IPI,NR,NF)
      REAL RTMD(IPI,NR,NF),PTMD(IPI,NR,NF)
      COMMON /BLKJ/ RHTE,PHTE,RHTM,PHTM
      COMMON /BLKJ/ KZY,AFND,AFTD,KZRD,KZID
      COMMON /BLKJ/ RTED,PTED,RTMD,PTMD
      REAL DLR(NL1,2),RLR(NL,2),DDPP(NL1),RL(NL),RA,W,WM,WP
      INTEGER NLYR
      COMMON /BLKL/ DLR,RLR,DDPP,RL,RA,NLYR,W,WM,WP
      COMPLEX ARG(NL1),AEX(NL1),K10,K1(NL),RZ(NL),COA(NL),Z1,ZA,T1
      REAL H(NL1)
      COMMON /BLKP/ ARG,AEX,K10,K1,RZ,COA,Z1,ZA,T1,H
C ___________________end MT includes____________________________________
c uses nrc, sitnam(), sitloc(), nfre, frek()
c uses idx
c _______________________start occam includes___________________________	
c Commons Revision 2.00, 13 May 1992
	real d, sd, dp
	integer nd, ndp
	common /data/ nd, d(ndd), sd(ndd), dp(ndd,4), ndp
	real pm, dm
	integer np
	common /model/ np, pm(npp), dm(ndd)
	real ptp, wjtwj, wjtwd, pwk1, pwk2, pwk3, frac, premod, prewts
	integer idebug, nfor, iout
	common /result/ ptp(npp,npp), wjtwj(npp,npp), wjtwd(npp),
     *  pwk1(npp), pwk2(npp), pwk3(npp), frac, nfor, iout, idebug,
     *  prewts(npp),premod(npp)
C ___________________end occam includes_________________________________
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      
      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
      implicit none
c _______________________start MT includes______________________________	
c Commons Revision 1.04, 13 Jan 93
      include 'occamdim.inc'
      include 'mt2Ddim.inc'
      integer nlay, nrcol, nrfix, icorn, nbrick, nexep
      integer ijexep, iprewt
      real thick, width, y, expen
      common /rgrid/ nlay, nrcol(nl), thick(npp), width(npp), 
     *  nrfix, y(npp+36), icorn(npp,4), 
     *  nbrick, nexep, ijexep(maxexp,2), expen(maxexp), iprewt(npp)
      integer nrc, nfre, nrl, ndata
      real sitloc, boffs, frek
      common /datinfo/ nrc, nrl(nr), ndata, nfre, boffs,  
     *  frek(nf), sitloc(nr)
      character*30 sitnam
      common /sitcom/ sitnam(nr)
      integer isx1, isx3, iscon, nfss, isfx, nsb
      real shift
      common /statics/ nsb, iscon, nfss, isx1(2*nr), isx3(2*nr), 
     *  isfx(2*nr), shift(2*nr)
      REAL YY(IPY),ZZ(IPZ),DLY(IPY),DLZ(IPZ)
      COMMON /BLKD/ YY,ZZ,DLY,DLZ
      INTEGER NPT(IPZ,4,IPY)
      REAL SIG(IPZ,4,IPY)
      COMPLEX SIGH(IPZ)
      COMMON /BLKE/ NPT,SIG,SIGH
      COMPLEX FPA(IPZ),FPB(IPZ),AFN(IPY),AFT(IPY)
      REAL PI,WU,WE
      COMMON /BLKF/ FPA,FPB,AFN,AFT,PI,WU,WE
      INTEGER IDX,IMD,NODEY,NODEZ,NTQ,NDZ,NDZ1,NTAO(IPY),NTA(IPY)
      INTEGER NZZB(2),NOLA(2),NPRM,LAPRM(IPI),ISEA
      COMMON /BLKI/ IDX,IMD,NODEY,NODEZ
      COMMON /BLKI/ NTQ,NDZ,NDZ1,NTAO,NTA
      COMMON /BLKI/ NZZB,NOLA,NPRM,LAPRM,ISEA
      COMPLEX KZY(NR,NF),AFND(IPI),AFTD(IPI)
      REAL RHTE(NR,NF),PHTE(NR,NF),RHTM(NR,NF),PHTM(NR,NF)
      REAL KZRD(IPI,NR,NF),KZID(IPI,NR,NF)
      REAL RTED(IPI,NR,NF),PTED(IPI,NR,NF)
      REAL RTMD(IPI,NR,NF),PTMD(IPI,NR,NF)
      COMMON /BLKJ/ RHTE,PHTE,RHTM,PHTM
      COMMON /BLKJ/ KZY,AFND,AFTD,KZRD,KZID
      COMMON /BLKJ/ RTED,PTED,RTMD,PTMD
      REAL DLR(NL1,2),RLR(NL,2),DDPP(NL1),RL(NL),RA,W,WM,WP
      INTEGER NLYR
      COMMON /BLKL/ DLR,RLR,DDPP,RL,RA,NLYR,W,WM,WP
      COMPLEX ARG(NL1),AEX(NL1),K10,K1(NL),RZ(NL),COA(NL),Z1,ZA,T1
      REAL H(NL1)
      COMMON /BLKP/ ARG,AEX,K10,K1,RZ,COA,Z1,ZA,T1,H
C ___________________end MT includes____________________________________
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()
c _______________________start occam includes___________________________	
c Commons Revision 2.00, 13 May 1992
	real d, sd, dp
	integer nd, ndp
	common /data/ nd, d(ndd), sd(ndd), dp(ndd,4), ndp
	real pm, dm
	integer np
	common /model/ np, pm(npp), dm(ndd)
	real ptp, wjtwj, wjtwd, pwk1, pwk2, pwk3, frac, premod, prewts
	integer idebug, nfor, iout
	common /result/ ptp(npp,npp), wjtwj(npp,npp), wjtwd(npp),
     *  pwk1(npp), pwk2(npp), pwk3(npp), frac, nfor, iout, idebug,
     *  prewts(npp),premod(npp)
C ___________________end occam includes_________________________________
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 + abs(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
      
      

