c***************
c exact solution for talwani 3d magnetics. plouff 7-77
c modified by r.godson for denver multics 3/81
c ability to process larger grids by row added by j.phillips 6/94
c
      character*4 unix,uniz
      character*8 name
      character*18 title
      character*40 idplot
      character*50 ifile,ifile2
      character*50 cafile,rfile,pfile
      character*80 ifmt
      character*32 id1(2),id2(2)
      common/st/cafile,rfile,pfile,ifile,ifile2,name,ifmt,
     1 idplot,unix,uniz
      common/maginp/iplotr,ibody,lsqs,dec,fincl,
     1 datum,xscale,dc,height,
     2 sus,rem,d,zi,iobs,icalc,ires
      common/field/v(1000,5),r(1000),q(1000,3)
      common/contg/nrow,ncol,dy,dx,yo,xo,xma,yma
      common/fpts/xp(1001),yp(1001),zp(1000)
      common /side/xb(900),yb(900), ll(50)
      common /bodk/ mm,nlim,jbody,short,nend
      common /tjcal/tjx(50),tjy(50),tjz(50),mtest
      common /vec/cde,cie,sde,sie,ex,ey,ez,earth,cdci,sdci
      common /trans/ inpt,iout,jin,iprint
      common /bod/ dstnc(900),ifprnt(50),zee(50),zm(50)
      dimension t(5),flsq(7),a(1000),b(1000)
      data flsq/'y','x','z','h','t','t','t'/
c dimension limits for 3 parameters.  corners, bodies, fieldpoints.
      data ncorn,nbod,ier,tx,ty,tz
     1 / 900,50,0,0.0,0.0,0.0/
      data in1/10/,in2/11/,iout1/12/,iout2/13/,nz/1/,dum/0./,
     1 id1/'magpoly calculated values from ',
     2 'model                           '/,
     3 id2/'magpoly residual values(observe',
     4 'd-calculated                    '/
      data ddval/1.e30/
c
c     get input parameters
c
      ierr=0
      call minp(*390)
      write(iout,10) idplot
   10 format (' exact solution for 3d magnetic bodies bounded by ',
     1 'horizontal and vertical planes',/,5x,a40)
      if(lsqs.ne.0) go to 20
      go to 40
   20 write (iout,30) flsq(lsqs)
   30 format (5x,'expect determination of standard deviation using ',
     1'inputted ',a1,'-component of magnetic field.')
      if (lsqs .lt. 6) go to 40
      if (earth  .eq. 0.0) earth=-50000.0
      if (fincl  .ne. 0.0) go to 40
      if (dec    .ne. 0.0) go to 40
      fincl=89.999
   40 facx=0.0
      facz=0.0
      call unit(unix,uniz,facx,facz)
      if (facx .eq. 0.0) go to 370
      zfac=1.0/facz
   50 call compon(1,0,dec,fincl,0.0,0.0,0.0,0.0,0.0)
      mtest=0
      short=1.0e20
      zdel =1.0e20
c
c     read body coordinates
c
      ierr=0
      call bodm(facz,nbod,ncorn,zdel,tx,ty,tz,ierr)
      zdel=abs(facz*zdel)
c
c mm (common) is total number of layers
      if (mm .eq. 0) go to 370
      if (mtest .eq. 0) go to 70
      if (mtest .eq.mm) go to 70
      write (iprint,60)
   60 format (5x,'magnetization not specified for all bodies')
      go to 370
   70 kk=0
c
c     read fieldpoint locations and maybe values
c
      call ptsm(kk,*390)
      title='body corners'
      if(ibody.ne.0) call mag3d(kk,*390)
      if(name.ne.'byrow') then
        print*,'Begin calculations...'
        go to 33
      endif
      open(iout1,file=cafile,status='unknown',form='unformatted')
      write(iout1) id1,ncol,nrow,nz,xo,dx,yo,dy
      if(lsqs.eq.1) then
        open(iout2,file=rfile,status='unknown',form='unformatted')
        write(iout2) id2,ncol,nrow,nz,xo,dx,yo,dy
      endif
      print*,'Calculating magnetic field by rows'
      print*
      mcol=ncol+1
      xo1=xo-(2*dx)
      yo1=yo-dy
      i=0
   31 i=i+1
      if(i.gt.nrow) go to 390
      write(6,420) i,nrow
      kk=0
      call rowin(xp,mcol,in1)
      if(ifile2.ne.' ') call rowin(yp,mcol,in2)
      do 32 j=2,mcol
      if(lsqs.ne.0) then
           if(xp(j).ge.ddval.or.yp(j).ge.ddval) stop 'dvals found'
      endif
      kk=kk+1
      r(kk)=xp(j)
      if(ifile2.ne.' ') height=yp(j)
      zp(kk)=height
      yp(kk)=xo1+(j*dx)
      xp(kk)=yo1+(i*dy)
c      write(iout,30) yf(k),xf(k),zf(k),r(k)
   32 continue
   33 if (kk .lt. 1) go to 370
c
c     check if printer plot of bodies & fieldpoints
c     desired if fieldpoints are in grid form
c
c      title='body corners'
c     if(ibody.ne.0.and.iplotr.eq.9.and.name.eq.'gridded')
c    1 call cplot(dc,a,1,0,title,*365)
c
c     check if location plot of bodies & fieldpoints on
c     plotting device wanted
c
c      if(ibody.ne.0) call mag3d(kk,*390)
c
c     check if contour plot of observed values desired
c
      if(lsqs.eq.0) go to 80
      title='observed values'
c     if(iobs.ne.0.and.iplotr.eq.9.and.name.eq.'gridded')
c    1 call cplot(dc,r,kk,1,title,*365)
c
   80 nk1=kk+1
      do 90 k=1,kk
   90 zp(k)=zp(k)*facz
      if(lsqs.eq.0) go to 110
      do 100 k=1,kk
      r(k)=r(k)-datum
  100 continue
  110 write (iout,120) kk
  120 format (//,' total of',i5,' fieldpoints.')
      call calc(kk,mm,short,zdel)
      if(lsqs.lt.6) go to 220
      go to 160
  130 if(ierr.ne.0) go to 360
      read (jin,*,end=360) icalc,earth,sus,rem,d,zi
      write(iout,150)
  150 format('1')
      call compon(1,0,dec,fincl,0.0,0.0,0.0,0.0,0.0)
      go to 210
  160 d=0.0
      zi=0.0
      rem=0.0
      if (nk1 .gt. 3) go to 180
      write (iprint,170)
  170 format ('0stop. insufficient fieldpoints for least-squares.')
      go to 370
  180 call tlsqs(zi,d,rem,nk1,lsqs,datum)
      if (nk1 .gt. 0) go to 190
      go to 370
  190 if (lsqs .eq. 6) go to 200
      sus=rem
      rem=0.0
      go to 210
  200 sus=0.0
  210 call compon(2,ier,d,zi,sus,rem,tx,ty,tz)
      if (ier .eq. -1) go to 370
      mtest=0
  220 std=0.0
      avg=0.0
      ifdp=5
      if(lsqs.ne.0.and.lsqs.lt.6) ifdp=lsqs
      sym1=flsq(ifdp)
      write (iout,230) datum,sym1
  230 format (/,1x,f10.1,' gamma is added to ',a1,'-output',
     1 /, 12x,'x',10x,'y',9x,'z',9x,'del-x',
     2 6x,'del-y',6x,'del-z',6x,'del-h',6x,'del-t',5x,'error')
      do 300 k=1,kk
      if (mtest .eq. mm ) go to 240
      vk1=v(k,1)
      vk2=v(k,2)
      vk3=v(k,3)
      vk5=v(k,4)
      vk6=v(k,5)
      vk4=-(vk1+vk6)
      t(1)=tx*vk1+ty*vk2+tz*vk3
      t(2)=tx*vk2+ty*vk4+tz*vk5
      t(3)=tx*vk3+ty*vk5+tz*vk6
      go to 250
  240 t(1)=q(k,1)
      t(2)=q(k,2)
      t(3)=q(k,3)
  250 t(4)=0.0
      t(5)=0.0
      call compon(3,0,0.0,0.0,t(4),t(5),t(1),t(2),t(3))
      fx=xp(k)
      fy=yp(k)
      zfk=zp(k)*zfac
      conout=t(ifdp)  +datum
      t(ifdp) =conout
c save for contour plot
      a(k)=conout
  260 if ( lsqs .eq. 0) go to 280
      errk=r(k)-conout+datum
c save for residual contour plot
      b(k)=errk
      std=std+errk*errk
      avg=avg+errk
      write (iout,270) k,fy,fx,zfk,t(2),t(1),(t(j),j=3,5),errk
  270 format (i5,f10.3,2f11.3,6f11.2)
      go to 300
  280 write (iout,290) k,fy,fx,zfk,t(2),t(1),(t(j),j=3,5)
  290 format (i5,f10.3,2f11.3,5f11.2)
  300 continue
      if (lsqs .eq. 0) go to 330
  310 if (nk1 .lt. 3) go to 330
      std=sqrt(std/(nk1-2))
      avg=avg/kk
      write (iout,320) std,avg
  320 format (//,' standard deviation=',e13.5,' and average difference='
     1 e13.5,' gammas')
  330 mtest=0
      if(name.ne.'byrow') go to 340
      write(iout1) dum,(a(k),k=1,kk)
      if(lsqs.eq.1) write(iout2) dum,(b(k),k=1,kk)
      go to 31
c
c     contour plots
c
  340 title='calculated values'
c     if (icalc.ne.0.and.iplotr.eq.9.and.name.eq.'gridded')
c    1 call cplot (dc,a,kk,1,title,*365)
      title='residual values'
c     if (ires.ne.0.and.iplotr.eq.9.and.lsqs.ne.0.and.name.eq.'gridded')
c    1 call cplot (dc,b,kk,1,title,*365)
      if(icalc.ne.0.and.name.eq.'gridded')
     1 call stdout(1,a)
      if(ires.ne.0.and.lsqs.ne.0.and.name.eq.'gridded')
     1 call stdout(2,b)
  350 if (mtest .eq. mm .and. lsqs.ge.6) go to 160
      if(lsqs.lt.6) go to 130
  360 continue
      go to 390
  365 write(iprint,375)
  375 format(' contour interval(dc) is 0.')
  370 write(iprint,380)
  380 format(' error - program abort')
  390 close(jin)
      close(iout)
      write(iprint,400) pfile
  400 format(/,' output is in file ',a50)
  420 format(1h+,'Working on row ',i5,' of ',i5)
      stop
      end
c***************
      subroutine minp(*)
c
c     namelist input parameter subroutine
c
c     r.godson - 3/81
c
      character*4 unix,uniz
      character*8 name
      character*50 cafile,rfile,pfile
      character*40 idplot
      character*50 ifile,ifile2,cfile
      character*80 ifmt
      common/st/cafile,rfile,pfile,ifile,ifile2,name,ifmt,
     1 idplot,unix,uniz
      common/maginp/iplotr,ibody,lsqs,dec,fincl,
     1 datum,xscale,dc,height,
     2 sus,rem,d,zi,iobs,icalc,ires
      common/contg/nrow,ncol,dy,dx,xmin,ymin,xmax,ymax
      common /vec/cde,cie,sde,sie,ex,ey,ez,earth,cdci,sdci
      common /trans/ inpt,iout,jin,iprint
      common/widect/naxcol
c     namelist/parms/iplotr,ibody,ifile,ifile2,lsqs,earth,dec,fincl,
c    1 datum,xscale,dc,unix,uniz,name,height,
c    2 ifmt,idplot,iobs,icalc,ires,naxcol
      inpt=5
      iout=8
      iprint=6
      jin=9
      iplotr=9
      naxcol=130
      ibody=0
      icalc=0
      iobs=0
      ires=0
      idplot=' '
      ifile=' '
      ifile2=' '
      lsqs=0
      earth=0.
      dec=0.
      fincl=0.
      datum=0.
      xscale=0.
      dc=0.
      unix='kilm'
      uniz='feet'
      name='gridded'
      height=0.
      write(iprint,10)
   10 format(' enter command file name: '$)
      read(inpt,20) cfile
   20 format(a)
      open(jin,file=cfile,status='old')
c     read(jin,parms,end=30)
      call namemp(jin)
      if(lsqs.eq.0) ires=0
      if(ifile.eq.' ') go to 50
c
c     determine output file names
c
      call names(cfile,cafile,rfile,pfile)
      open(iout,file=pfile,status='unknown')
      return
   50 write(iprint,60)
   60 format(' no fieldpoint data file name specified')
   70 return 1
      end
c***************
      subroutine names(cfile,cafile,rfile,pfile)
      character*(*) cfile,cafile,rfile,pfile
      ic=index(cfile,'.')-1
      pfile=cfile(:ic)//'.prt'
      cafile=cfile(:ic)//'.cal'
      rfile=cfile(:ic)//'.res'
      return
      end
c***************
      subroutine bodm(facz,nbod,ncorn,zdel,tx,ty,tz,ierr)
c
c      reads body corner coordinates
c
c      r.godson - 4/81
c
      character*4 unix,uniz
      character*8 name
      character*40 idplot
      character*50 ifile,ifile2
      character*50 cafile,rfile,pfile
      character*80 ifmt
      common/mplot/ymin,xmin,ymax,xmax
      common/contg/nrow,ncol,dy,dx,yo,xo,xma,yma
      common/st/cafile,rfile,pfile,ifile,ifile2,name,ifmt,
     1 idplot,unix,uniz
      common/maginp/iplotr,ibody,lsqs,dec,fincl,
     1 datum,xscale,dc,height,
     2 sus,rem,d,zi,iobs,icalc,ires
      common /bod/ dstnc(900),ifprnt(50),zee(50),zm(50)
c see ncorn dimension limit in main program
      common /side/xb(900),yb(900), ll(50)
      common /bodk/ mm,nlim,jbody,short,nend
      common /tjcal/tjx(50),tjy(50),tjz(50),mtest
      common /trans/ inpt,iout,jin,iprint
      nend=0
      mm=0
      xmin=1.0e38
      xmax=-1.0e38
      ymin=1.0e38
      ymax=-1.0e38
   10 read(jin,*,end=200)num,ifprnt(mm+1),sus,rem,d,zi,top,bottom
      if(num.eq.0) go to 190
      nst=nend+1
      nend=nend+num
      mm=mm+1
      if(mm.gt.nbod) go to 140
      if(nend.gt.ncorn) go to 70
      ll(mm)=num
      zee(mm)=top*facz
      zm(mm)=bottom*facz
      write(iout,30) mm,top,bottom,num
   30 format (' body',i3,' between depths',f11.3,' and',f11.3,
     1 ' with',i4,' corners.')
   40 if (top .gt. bottom) go to 90
   50 write (iprint,60) top,bottom
   60 format ('0top of body is below bottom',2f7.2)
      go to 180
   70 write(iprint,80) mm,nend,ncorn
   80 format (' body',i3,' includes sufficient corners,',i4,',to exceed'
     1 , ' limitation of',i5,' for all bodies.')
      go to 180
c coordinates of corners of body east,north (x,y) pairs.
c body has at least 3 corners
   90 read(jin,*,end=160) (yb(j),xb(j),j=nst,nend)
      write(iout,100) (yb(j),xb(j),j=nst,nend)
  100 format (5(5x,2f9.3))
      do 105 j=nst,nend
      if(xb(j).lt.xmin) xmin=xb(j)
      if(xb(j).gt.xmax) xmax=xb(j)
      if(yb(j).lt.ymin) ymin=yb(j)
      if(yb(j).gt.ymax) ymax=yb(j)
  105 continue
      if (lsqs.ge.6) go to 120
      ier=0
      call compon(2,ier,d,zi,sus,rem,tx,ty,tz)
      if (ier .gt. 0) go to 110
      if(ier.eq.-1) go to 180
      go to 120
  110 mtest=mtest+1
      tjx(mm)=tx
      tjy(mm)=ty
      tjz(mm)=tz
  120 np=nend+1
      xb(np)=xb(nst)
      yb(np)=yb(nst)
      do 130 j=nst,nend
      dist=(xb(j+1)-xb(j))**2+(yb(j+1)-yb(j))**2
      if (dist .eq. 0.0) go to 130
      dist=sqrt(dist)
      if (dist .lt. short) short=dist
  130 dstnc(j)=dist
      zd=top-bottom
      if ( zd  .lt.  zdel) zdel=zd
      go to 10
  140 write (iprint,150) mm
  150 format (5x,'number of bodies=',i3,' exceeds limitation.')
      go to 180
  160 write(iprint,170) mm
  170 format(' end of file reading x,y pairs for body no. ',i3)
  180 mm=0
  190 return
  200 ierr=1
      return
      end
c***************
      subroutine compon (num,ier,dec,fin,sus,rem,tx,ty,tz)
c various components of earth and magnetization fields
      common /vec/cde,cie,sde,sie,ex,ey,ez,earth,cdci,sdci
      common /trans/ inpt,iout,jin,iprint
      go to (10,90,190),num
c earths field
   10 arth=abs(earth)
      if (earth .eq. 0.0) go to 50
   20 write (iout,30) dec,fin
   30 format (5x,'declination=',f7.1,' and inclination=',f6.1,
     1 ' degrees for earths field')
      if (arth .ne. 0.0) write (iout,40) arth
   40 format (9x,'of',f7.0,' gammas.')
      d=1.745329e-2*dec
      z=1.745329e-2*fin
      cde= cos(d)
      cie= cos(z)
      sie= sin(z)
      sde= sin(d)
      cdci=cde*cie
      sdci=sde*cie
      if (earth .le. 0.0) go to 70
      eh=earth*cie
      ex=earth*cdci
      ey=earth*sdci
      ez=earth*sie
      return
   50 write (iout,60)
   60 format (5x,'magnetic susceptibilities are assumed to be negligible
     1 .',//)
      go to 20
   70 write (iout,80)
   80 format (5x,'magnetic anomaly is assumed to be small compared to ',
     1 'earths field.',//)
      return
c magnetization parameters
   90 tx=0.0
      ty=0.0
      tz=0.0
c     if (sus .gt. 0.0) go to 6
      if (sus .ne. 0.0) go to 130
      if (rem .eq. 0.0) go to 140
      write (iout,100) rem,dec,fin
  100 format (5x,'intensity of total magnetization is',f9.2,'e-5 ',
     1 'with a declination=',f7.1,' and inclination=',f6.1,
     2 ' degrees.',//)
  110 d=1.745329e-2*dec
      z=1.745329e-2*fin
      cdr= cos(d)
      cir= cos(z)
      sir= sin(z)
      sdr= sin(d)
      ehr=rem*cir
      tx=tx+ehr*cdr
      ty=ty+ehr*sdr
      tz=tz+rem*sir
  120 ier= 1
      return
  130 if (earth .ne. 0.0) go to 160
  140 ier=-1
      write (iprint,150)
  150 format (5x,'error. normal field of earth or magnetization not ',
     1 'supplied')
      return
  160 write (iout,170) sus
  170 format (5x,'magnetic susceptibility is',f9.2,'e-5 cgs.',//)
      sus=arth*sus*1.0e-5
      tx=sus*cdci
      ty=sus*sdci
      tz=sus*sie
      if (rem .eq. 0.0) go to 120
      write (iout,180) rem,dec,fin
  180 format (5x,'intensity of remanent magnetization is',f9.2,'e-5 ',
     1 'with a declination=',f7.1,' and inclination=',f6.1,
     2 ' degrees.',//)
      go to 110
  190 if (earth .le. 0.0) go to 200
      delh=(ex+tx  )**2+(ey+ty  )**2
      delt=(ez+tz  )**2+delh
      sus=sqrt(delh)-eh
      rem=sqrt(delt)-earth
      return
  200 sus=tx  *cde+ty  *sde
      rem=tx*cdci+ty*sdci+tz*sie
      return
      end
c***************
      subroutine tlsqs(azi,ad,tm,nk1,lsqs,datum)
      real*8 fn,sx,sy,sxx,sxy,syy,sz,sxz,syz,d,a2,a1
      real*8 st,stt,stx,sty,stz,szz,a3,b1,b2,b3,c1,c2,c3,det
      real*8 ptx,pty,ptz,pyz,pxz,pxy,sgx,sgy,sgz,sig,p
      integer nm(8)/'x','y','z','t','n','d','r','s'/
      common /vec/cde,cie,sde,sie,ex,ey,ez,earth,cdci,sdci
      common/field/v(1000,5),r(1000),q(1000,3)
      common /trans/ inpt,iout,jin,iprint
c      dimension a(1),b(1)
      p(fn,sxy,sx,sy)=(sxy-sx*sy/fn)/fn
      det(a1,a2,a3,b1,b2,b3,c1,c2,c3)=a1*(b2*c3-b3*c2)
     1 +a2*(b3*c1-b1*c3)+a3*(b1*c2-b2*c1)
      sig(fn,sx,sxx)=(sxx-sx*sx/fn)/fn
c observed value r
      kk=nk1-1
      fn=kk
      sx =0.0d0
      sy =0.0d0
      syy=0.0d0
      sxy=0.0d0
      sxx=0.0d0
      l=5
      if (lsqs .eq. 6) go to 40
      if (kk .lt. 3) go to 180
      eyy=sdci*sdci
      e1=cdci*cdci-eyy
      e6=sie*sie-eyy
      e2=2.0*cdci*sdci
      e3=2.0*cdci*sie
      e5=2.0*sdci*sie
c determination of best-fit susceptibility
      do 10 k=1,kk
      rk=r(k)
      ak =     e1*v(k,1)+e2*v(k,2)+e3*v(k,3)+e5*v(k,4)+e6*v(k,5)
      sx =sx +ak
      sy =sy +rk
      syy=syy+rk*rk
      sxy=sxy+ak*rk
   10 sxx=sxx+ak*ak
      xx=sxx-sx*sx/fn
      l=1
      if (xx .le. 0.0) go to 180
      l=4
      xy=sxy-sx*sy/fn
      yy=syy-sy*sy/fn
      tm=xy/xx
      tb=(sy-tm*sx)/fn
      coer=9999.0
      if (yy .le. 0.0) go to 20
      coer=xy/sqrt(xx*yy)
   20 eyy=abs(earth)*1.0e-5
      if (eyy .eq. 0.0) eyy=0.5
      stdm=yy*(1.0-coer*coer)/(xx*(kk-2))
      l=8
      if (stdm .lt. 0.0) go to 180
      stdm=sqrt(stdm)
      tm=tm/eyy
      std=sqrt(yy /(kk -1))/eyy
      tb=tb+datum
      write (iout,30) tm,std,coer,tb,stdm
   30 format (//,' least-squares susceptibility is',f9.1,'e-5 emu. std',
     1 ' dev=',f7.0,' gammas. correl coef=',f6.2,' zero',f6.0,' gammas',
     2 /,32x,f9.1)
c     if (tm .le. 0.1) go to 7
      if (tm .gt. 0.7e5) go to 180
      azi=0.0
      ad=0.0
      return
c least squares determination of total magnetization vector
   40 if (kk .lt. 5) go to 180
      sz =0.0d0
      sxz=0.0d0
      syz=0.0d0
      st =0.0d0
      stt=0.0d0
      stx=0.0d0
      sty=0.0d0
      stz=0.0d0
      szz=0.0d0
      do 50 k=1,kk
      vk1=v(k,1)
      vk2=v(k,2)
      vk3=v(k,3)
      vk5=v(k,4)
      vk6=v(k,5)
      t=r(k)
      x= cdci*vk1+sdci*vk2+sie*vk3
      y= cdci*vk2-sdci*(vk1+vk6)+sie*vk5
      z= cdci*vk3+sdci*vk5+sie*vk6
      st =st +t
      sx =sx +x
      sy =sy +y
      sz =sz +z
      stt=stt+t*t
      stx=stx+t*x
      sty=sty+t*y
      stz=stz+t*z
      szz=szz+z*z
      sxx=sxx+x*x
      sxy=sxy+x*y
      sxz=sxz+x*z
      syz=syz+y*z
      syy=syy+y*y
   50 continue
      ptx=p(fn,stx,st,sx)
      pty=p(fn,sty,st,sy)
      ptz=p(fn,stz,st,sz)
      pyz=p(fn,syz,sy,sz)
      pxz=p(fn,sxz,sx,sz)
      pxy=p(fn,sxy,sx,sy)
      sgx=sig(fn,sx,sxx)
      sgy=sig(fn,sy,syy)
      sgz=sig(fn,sz,szz)
      sgt=sig(fn,st,stt)
      d=det(sgx,pxy,pxz,pxy,sgy,pyz,pxz,pyz,sgz)
      l=6
      if (d  .eq. 0.0d0) go to 180
      x=det(ptx,pty,ptz,pxy,sgy,pyz,pxz,pyz,sgz)/d
      y=det(sgx,pxy,pxz,ptx,pty,ptz,pxz,pyz,sgz)/d
      z=det(sgx,pxy,pxz,pxy,sgy,pyz,ptx,pty,ptz)/d
      a0=-(x*sx+y*sy+z*sz-st)/fn
      tb=a0+datum
      std=x*ptx+y*pty+z*ptz
      coer=std/sgt
      l=7
      if (coer .lt. 0.0) go to 180
      coer=sqrt(coer)
      std=sgt-std
      l=8
      if (std  .lt. 0.0) go to 180
      std =sqrt(std )
      write(iout,60)y,x,z,std,coer,tb
   60 format (//,' least squares determined values. jx=',f8.1,'e-5  jy='
     1,f8.1,'e-5, jz=',f8.1,'e-5',/,5x,' (assuming anomaly is negligibl'
     2, 'e compared to total field of earth)',/,5x,'std',
     3 ' dev=',f7.0,' gammas. correl coef=',f6.2,' zero',f6.0,' gammas')
      th=x*x+y*y
      tm=sqrt(th+z*z)
      th=sqrt(th)
      if (th .gt. 0.0) go to 80
      if (z  .gt. 0.0) go to 70
      azi=-90.
      go to 170
   70 azi=90.
      go to 170
   80 azi=atan (z/th)*57.29578
      if(x) 120,90,130
   90 if(y)100,100,110
  100 ad=-90.
      go to 150
  110 ad=90.
      go to 150
  120 gac=180.
      go to 140
  130 gac=0.0
  140 ad=atan (y/x)*57.29578+gac
  150 write         (iout,160)tm,th,azi,ad
  160 format (/,'   jt=',f8.1,'e-5,  jh=',f8.1,'e-5, zi=',f6.1,', dec=',
     1 f7.1)
      l=4
c     if (tm .le. 0.1) go to 7
      return
  170 ad=0.0
      go to 150
  180 write (iprint,190) nm(l)
  190 format ('0stop. ',a1,'-component of summed field or amplitude ',
     1 'incorrect.')
      nk1=-nk1
      return
      end
c***************
      subroutine calc(kk,mm,short,zdel)
c calculate exact 3d magnetic values
      common /tjcal/tjx(50),tjy(50),tjz(50),mtest
      common/fpts/xp(1001),yp(1001),zp(1000)
      common/field/v(1000,5),r(1000),q(1000,3)
      common /trans/ inpt,iout,jin,iprint
      common /bod/ dstnc(900),ifprnt(50),zee(50),zm(50)
      common /side/xb(900),yb(900), ll(50)
      dimension s(6)
c  part  of shortest length of body edges and thinnest layer.
      stest=short*2.0e-4
      ztest= zdel*2.0e-4
      ztem= 10.0*ztest
      do 250 k=1,kk
      xf =xp(k)
      yf =yp(k)
      zf=zp(k)
      delx=0.0
      dely=0.0
      delz=0.0
      do 10 j=1,5
   10 s(j)=0.0
      nend=0
c number of bodies is mm
      do 230 m =1,mm
      zf1=zee(m)
      zf2= zm(m)
      z1=zf1-zf
      z2=zf2-zf
      if ( abs(z1) .lt. ztest) z1=0.0
      if ( abs(z2) .lt. ztest) z2=0.0
      fact= 1.0
      zz=z1+z2
      z1= abs(z1)
      z2= abs(z2)
      if (zz) 20,230,30
c case where most of mass is above fieldpoint
   20 fact=-1.0
   30 if (z1 .lt. z2) go to 40
      dtest=z2
      z2=z1
      z1=dtest
   40 lk=ll(m)
      ns=nend+1
      nend=ns+lk-1
      z1s=z1*z1
      z2s=z2*z2
      aa1=0.0
      aa2=0.0
      aa3=0.0
      aa5=0.0
      aa6=0.0
      n=ns
      x1=xb(n)-xf
      y1=yb(n)-yf
      if ( abs(x1) .lt. stest) x1=0.0
      if ( abs(y1) .lt. stest) y1=0.0
      r1s=x1*x1+y1*y1
      r1z1= sqrt(r1s+z1s)
      r1z2= sqrt(r1s+z2s)
      r1sf=r1s
      rz1f=r1z1
      rz2f=r1z2
      x1f=x1
      y1f=y1
c number of line segments (vertical sides) or vertices =lk-1
      do 190 l=1,lk
      nm=n
      if (l .ne. lk) go to 50
      x2=x1f
      y2=y1f
      r2s=r1sf
      r2z1=rz1f
      r2z2=rz2f
      go to 60
   50 n=n+1
      x2=xb(n)-xf
      y2=yb(n)-yf
      if ( abs(x2) .lt. stest) x2=0.0
      if ( abs(y2) .lt. stest) y2=0.0
      r2s=x2*x2+y2*y2
      r2z1= sqrt(r2s+z1s)
      r2z2= sqrt(r2s+z2s)
   60 dx= x2- x1
      dy= y2- y1
      dds= dx* dx+ dy* dy
      if (dds .le. 0.0  ) go to 180
      ds=dstnc(nm)
      cost=dy/ds
      sint=dx/ds
      sncos=sint*cost
      coss=cost*cost
c     sins=sint*sint
      p=(x1*y2-x2*y1)/ds
c test if perpendicular distance to line exceeds 1/5000 segment length
      if ( abs(p) .gt. stest) go to 70
      p=0.0
   70 d1=  x1*sint+y1*cost
      d2=  x2*sint+y2*cost
      dd1= abs(d1)
      dd2= abs(d2)
      if (dd1      .lt. stest) d1=0.0
      if (dd2      .lt. stest) d2=0.0
      d12=d1*d2
c  log(r+d) terms
      if (0.0   .eq. p) go to 140
   80 ph=p*p+z1s
      rdz1=r1z1+d1
      if (d1 .ge. 0.0  ) go to 100
      if (d2 .gt. 0.0  ) go to 90
      gl=(r2z1-d2)/(r1z1-d1)
      go to 160
   90 dtest=ph/(d1*d1)
      if (dtest .lt. 0.01  ) rdz1=dd1*0.5  *dtest*(1.0  -0.25  *dtest)
  100 rdz2=r2z1+d2
      if (d2 .ge. 0.0  ) go to 110
      dtest=ph/(d2*d2)
      if (dtest .lt. 0.01  ) rdz2=dd2*0.5  *dtest*(1.0  -0.25  *dtest)
  110 gl=rdz1/rdz2
      go to 160
  120 write (iout,130) yf,xf,ztem
  130 format (' fieldpoint at x=',f10.2,', y=',f10.2,' is raised',f9.5,
     1 ' distance units (on edge of body)' )
      z1=ztem
      z1s=z1*z1
      r1z1=sqrt(r1s+z1s)
      r2z1=sqrt(r2s+z1s)
      if (l .eq. 1) rz1f=r1z1
      go to 80
  140 if (0.0   .ne. z1) go to 80
      if (0.0   .ge. d12) go to 120
      if (0.0   .gt. d1) go to 150
c both positive, exterior, in plane of edge, depth zero
      gl=d1/d2
      go to 160
c both negative, exterior, in plane of edge, depth zero
  150 gl=d2/d1
  160 gl=alog(gl*(r2z2+d2)/(r1z2+d1))
c  log(r+z) terms
      fl=(z2+r2z2)/(z2+r1z2)
      if (0.0   .lt. r1z1) fl=fl*(z1+r1z1)
      if (0.0   .lt. r2z1) fl=fl/(z1+r2z1)
      fl=alog(fl)
c  arctangent(zd/pr) terms
      t12=0.0
      if (0.0   .eq. p) go to 170
      call tanget(p,z2*d2,r2z2,z1*d2,r2z1,z1*d1,r1z1,z2*d1,r1z2,t12)
  170 aa1=aa1- coss*t12+sncos*fl
      aa2=aa2+sncos*t12+ coss*fl
      aa3=aa3+      cost*gl
c     aa4=aa4- sins*t12-sncos*fl
      aa5=aa5-      sint*gl
      aa6=aa6+t12
  180 x1=x2
      y1=y2
      r1s=r2s
      r1z1=r2z1
      r1z2=r2z2
c end, corner loop
  190 continue
      s(1)=s(1)+fact*aa1
      s(2)=s(2)+fact*aa2
      s(3)=s(3)+fact*aa3
      s(4)=s(4)+fact*aa5
      s(5)=s(5)+fact*aa6
      aa4=-(aa1+aa6)
      if ( mtest    .eq. 0) go to 210
      tjxm=tjx(m)
      tjym=tjy(m)
      tjzm=tjz(m)
      delx=delx+fact*(tjxm  *aa1+tjym  *aa2+tjzm  *aa3)
      dely=dely+fact*(tjxm  *aa2+tjym  *aa4+tjzm  *aa5)
      delz=delz+fact*(tjxm  *aa3+tjym  *aa5+tjzm  *aa6)
      if (ifprnt(m) .eq. 0) go to 230
      delh=0.0
      delt=0.0
      tx=delx
      ty=dely
      tz=delz
      call compon(3,0,0.0,0.0,delh,delt,tx,ty,tz)
      write (iout,200) k,yf,xf,zf,dely,delx,delz,delh,delt,ifprnt(m)
  200 format (i5,3f11.4,5f11.2,i3)
  210 if (ifprnt(m) .lt. 0)
     1 write (iout,220) k,aa1,aa2,aa3,aa4,aa5,aa6,m,p,dy,dx
  220 format(i6,6e15.8, i4,f8.2,2f7.3)
c end, body loop
  230 continue
      do 240 n=1,5
  240  v(k,n)=s(n)
      if ( mtest    .eq. 0) go to 250
      q(k,1)=delx
      q(k,2)=dely
      q(k,3)=delz
c end, fieldpoint loop
  250 continue
      return
      end
c***************
      subroutine unit(unix,uniz,facx,facz)
c conversion of horizontal (x) and vertical (z) distances-3d magnetics
      character*4 unix,uniz,feet,kilm,mile,kilf,metr
c     1 lunx,lunz
      common /trans/ inpt,iout,jin,iprint
      data kilm,  feet,  mile,  kilf,  metr
     1   /'kilm','feet','mile','kilf','metr'/
      write (iout,10) unix,uniz
   10 format (5x,'units of distance measurement are in ',a4,
     1 ' and heights are in ',a4)
   20 if (unix .eq. mile) facx=1.609347
      if (uniz .eq. mile) facz=1.609347
      if (unix .eq. feet) facx=3.048006e-4
      if (uniz .eq. feet) facz=3.048006e-4
      if (unix .eq. kilf) facx=3.048006e-1
      if (uniz .eq. kilf) facz=3.048006e-1
      if (unix .eq. kilm) facx=1.0
      if (uniz .eq. kilm) facz=1.0
      if (unix .eq. metr) facx=1.0e-3
      if (uniz .eq. metr) facz=1.0e-3
      if (facx .ne. 0.0) go to 40
      write(iprint,30)
   30 format(' unix specification is incorrect')
      return
   40 if (facz .ne. 0.0) go to 50
      facx=0.0
      return
c multiplying factor for heights to convert to distance units.
c negative sign to compensate for depths + downward later 'calc'.
   50 facz=-facz/facx
      return
      end
c***************
      subroutine tanget(p,zd1,r1,zd2,r2,zd3,r3,zd4,r4,t12)
c arctangent addition with alternating terms negative
      dimension zd(4),r(4)
      r(1)=r1
      r(2)=r2
      r(3)=r3
      r(4)=r4
      zd(1)=zd1
      zd(2)=zd2
      zd(3)=zd3
      zd(4)=zd4
      i=1
      sum=0.0
       tg=0.0
      do 60 j=1,4
      fnum=zd(j)
      if (fnum .eq. 0.0  ) go to 60
c subroutine by-passed for p=0
      f=fnum/(p*r(j))
      if (i .eq. -1) f=-f
      fnum=tg+f
      den=1.0d0-tg*f
      if (den  .eq. 0.0  ) go to 30
       tg=fnum /den
      if (den  .gt. 0.0  ) go to 60
      if (fnum ) 10, 60,20
   10 sum   =    sum-3.141593
      go to 60
   20 sum   =    sum+3.141593
      go to 60
   30 if (fnum .lt. 0.0  ) go to 40
      sum=sum+1.570796
      go to 50
   40 sum=sum-1.570796
   50 tg=0.0
   60 i=-i
      t12= atan(tg)+sum
      return
      end
c***************
      subroutine ptsm(npts,*)
c
c     input data subroutine for user formatted or gridded data
c     r.godson - 3/81
c
      character*4 unix,uniz
      character*8 name
      character*40 idplot
      character*50 ifile,ifile2
      character*64 id,id1
      character*50 cafile,rfile,pfile
      character*80 ifmt
      common/mplot/ymin,xmin,ymax,xmax
      common/fpts/xf(1001),yf(1001),zf(1000)
      common/st/cafile,rfile,pfile,ifile,ifile2,name,ifmt,
     1 idplot,unix,uniz
      common/maginp/iplotr,ibody,lsqs,dec,fincl,
     1 datum,xscale,dc,height,
     2 sus,rem,d,zi,iobs,icalc,ires
      common/contg/nrow,ncol,dy,dx,yo,xo,xma,yma
      common/field/v(1000,5),r(1000),q(1000,3)
      common /trans/ inpt,iout,jin,iprint
      dimension z(1001),z2(1001)
      data in/10/,in2/11/,nrcmax/1000/
      write(iout,10)
   10 format(//,14x,'observed values',/,
     1 7x,'x',10x,'y',10x,'z',10x,'m',/)
      if(name.eq.'gridded') go to 50
c
c     user formatted data
c
      open(in,file=ifile,status='old')
      i=1
   20 read(in,ifmt,end=40) yf(i),xf(i),r(i),zf(i)
      write(iout,30) yf(i),xf(i),zf(i),r(i)
   30 format(1x,4f11.4)
      if(xf(i).lt.xmin) xmin=xf(i)
      if(xf(i).gt.xmax) xmax=xf(i)
      if(yf(i).lt.ymin) ymin=yf(i)
      if(yf(i).gt.ymax) ymax=yf(i)
      i=i+1
      go to 20
   40 k=i-1
      go to 230
c
c     gridded data
c
   50 open(in,file=ifile,status='old',form='unformatted')
      read(in,end=190) id,ncol,nrow,nz,xo,dx,yo,dy
      if((ncol*nrow).gt.nrcmax) name='byrow'
      if(ncol.gt.nrcmax) go to 110
      if(dx.eq.0..or.dy.eq.0.) go to 130
      if(nz.ne.1) go to 170
      if(ifile2.eq.' ') go to 51
      open(in2,file=ifile2,status='old',form='unformatted')
      read(in2,end=190) id1,ncol1,nrow1,nz1,xo1,dx1,yo1,dy1
      if(ncol.ne.ncol1.or.nrow.ne.nrow1) go to 90
      if(dx.ne.dx1.or.dy.ne.dy1.or.xo.ne.xo1.or.yo.ne.yo1) go to 90
      if(nz.ne.nz1) go to 90
   51 npts=0
      if(name.eq.'byrow') go to 280
   60 mcol=ncol+1
      xo1=xo-(2*dx)
      yo1=yo-dy
      k=0
      do 80 i=1,nrow
      call rowin(z,mcol,in)
      if(ifile2.ne.' ') call rowin(z2,mcol,in2)
      do 70 j=2,mcol
      k=k+1
      r(k)=z(j)
      if(ifile2.ne.' ') height=z2(j)
      zf(k)=height
      yf(k)=xo1+(j*dx)
      xf(k)=yo1+(i*dy)
      write(iout,30) yf(k),xf(k),zf(k),r(k)
   70 continue
   80 continue
      if(xf(1).lt.xmin) xmin=xf(1)
      if(xf(k).gt.xmax) xmax=xf(k)
      if(yf(1).lt.ymin) ymin=yf(1)
      if(yf(k).gt.ymax) ymax=yf(k)
      xma=xf(k)
      yma=yf(k)
      go to 230
   90 write(iprint,100)
  100 format(' parameters on two input standard files do not match')
      go to 260
  110 write(iprint,120) nrcmax
c  120 format(' number of grid points .gt. allowed maximum of ',i4)
  120 format(' number of columns .gt. allowed maximum of ',i4)
      go to 260
  130 write(iprint,140)
  140 format(' present implementation does not allow unequal',
     1 ' grid intervals')
      go to 260
  150 write(iprint,160) nrcmax
  160 format(' the number of columns * number of rows requested is',
     1 ' greater than the maximum allowed of ',i5)
      go to 260
  170 write(iprint,180)
  180 format(' only single precision values allowed')
      go to 260
  190 write(iprint,200)
  200 format(' end of file while reading input data')
      go to 260
  210 write(iprint,220) k,nrcmax
  220 format(' number of formatted data values',i5,
     1 ' is greater than allowed number of ',i5)
      go to 260
  230 npts=k
      if(npts.lt.1) go to 240
      ierr=0
      go to 270
  240 write(iprint,250)
  250 format(' number of field points less than 1')
  260 close(in)
      if(ifile2.ne.' ') close(in2)
      return 1
  270 close(in)
      if(ifile2.ne.' ') close(in2)
  280 return
      end
c***************
      subroutine rowin(z,mcol,in)
      dimension z(mcol)
      read(in) z
      return
      end
c***************
      subroutine mag3d(npts,*)
c
c     plotting routine for body corners & fieldpoints
c     r.godson - 4/81
c
      character*4 unix,uniz
      character*1 icross
      character*8 name
      character*40 idplot
      character*50 ifile,ifile2
      character*50 cafile,rfile,pfile
      character*80 ifmt
      integer sym(50),iddplot(10)
      equivalence(idplot,iddplot(1))
      equivalence(icross,icr)
      common/mplot/xmin,ymin,xmax,ymax
      common/fpts/xf(1001),yf(1001),zf(1000)
      common/st/cafile,rfile,pfile,ifile,ifile2,name,ifmt,
     1 idplot,unix,uniz
      common/maginp/iplotr,ibody,lsqs,dec,fincl,
     1 datum,xscale,dc,height,
     2 sus,rem,d,zi,iobs,icalc,ires
      common /bodk/ mm,nlim,jbody,short,nend
      common /side/xb(900),yb(900), ll(50)
      common /trans/ inpt,iout,jin,iprint
      common/contg/nrow,ncol,dy,dx,yo,xo,yma,xma
      dimension xd(2),yd(2),xs(4),ys(4),jdim(2),
c     1 ifmt1(2),nscale(5),xxb(900),yyb(900)
     1 xxb(900),yyb(900)
      data icross/'+'/
c      data jdim(1)/0/,nscale/1,0,2,0,2/,ifmt1/'(f10','.1) '/,
c     1 nshl/6/,nsvl/10/,nvl/12/,
c     2 nx/2/,ny/2/,vexag/1./,
      data jdim(1)/0/,
     2 vexag/1./,
     3 sym/'1','2','3','4','5','6','7','8','9','a','b','c','d','e',
     4 'f','g','h','i','j','k','l','m','n','o','p','q','r','s',
     5 't','u','v','w','x','y','z',15*'*'/
c
c
c         denver plotting system
c
   10 if(iplotr.ne.5) go to 20
      if(xscale.eq.0.) go to 170
      nxl=(xmax-xmin)/xscale*2.+.5
      if(nxl.lt.2) nxl=2
      nyl=nxl
      go to 30
   20 nxl=24
      nyl=17
   30 call scale3(xmin,xmax,nxl,xminp,xmaxp,xdel)
      call scale3(ymin,ymax,nyl,yminp,ymaxp,ydel)
      xd(1)=xminp
      xd(2)=xmaxp
      yd(1)=yminp
      yd(2)=ymaxp
      if(iplotr.eq.5) then
      xs(1)=(xmaxp-xminp)/xscale
      ys(1)=(ymaxp-yminp)/xscale
      if(ys(1).gt.18.) ys(1)=18.
      else
      xs(1)=8.
      ys(1)=7.
      endif
      go to 60
   40 nyl=(ymax-ymin)/xscale*vexag*4.+.5
      if(nyl.lt.2) nyl=2
   50 call scale3(ymin,ymax,nyl,yminp,ymaxp,ydel)
      yd(1)=ymaxp
      yd(2)=yminp
      xs(1)=(xmaxp-xminp)/xscale
      ys(1)=(ymaxp-yminp)/xscale*vexag
   60 call pltset(iplotr,xs(4),ys(4),jdim)
      yyss=ys(4)
      xs(2)=0.
      ys(2)=0.
      xs(3)=1.
      ys(3)=1.
      if((xs(1)+xs(3)+.2).gt.xs(4)) go to 70
      if((ys(1)+ys(3)+.2).gt.ys(4)) go to 70
      go to 100
   70 xs(4)=amin1(xs(4),10.)
   80 ys(4)=amin1(ys(4),8.)
   90 xs(3)=1.0
      ys(3)=1.0
      xs(1)=xs(4)-xs(3)-.2
      ys(1)=ys(4)-ys(3)-.2
      go to 110
  100 xs(4)=xs(1)+xs(3)+.2
      ys(4)=ys(1)+ys(3)+.2
  110 call scale(xd,yd,xs,ys,4,icode)
      if(icode.lt.0) go to 150
      call neatl
      call xaxis(xd,yd,xs,xdel,2,.08,'(f10.2)',10)
      call yaxis(yd,xd,ys,ydel,2,.08,'(f10.2)',10)
      call vchar(2.,.15,iddplot,40,3,.15,0.,0.,0.)
  120 ie=0
c
c     plot body corners
c
      do 140 i=1,mm
      is=ie+1
      ie=ie+ll(i)
      call vchar(yb(is),xb(is),sym(i),1,200,.08,0.,0.,0.)
      m=0
      do 130 k=is,ie
      m=m+1
      yyb(m)=yb(k)
      xxb(m)=xb(k)
  130 continue
      xxb(m+1)=xb(is)
      yyb(m+1)=yb(is)
      call line(yyb,xxb,m+1,0,200)
  140 continue
c
c     plot fieldpoints
c
      do 145 i=1,npts
      call vchar(yf(i),xf(i),icr,1,300,.08,0.,0.,0.)
  145 continue
      jdim(2)=2
      call endpt(jdim)
      return
  150 write(iprint,160)
  160 format(' unable to scale plotter')
      go to 190
  170 write(iprint,180)
  180 format(' xscale is zero')
  190 return 1
      end
c***************
      subroutine stdout(isw,z)
c
c     subroutine to create standard file
c     r.godson - 4/81
c
      character*4 unix,uniz
      character*8 name
      character*32 id1(2),id2(2)
      character*40 idplot
      character*50 file1,file2,file3,ifile,ifile2
      character*80 ifmt
      dimension z(1)
      common/st/file1,file2,file3,ifile,ifile2,name,ifmt,
     1 idplot,unix,uniz
      common/maginp/iplotr,ibody,lsqs,dec,fincl,
     1 datum,xscale,dc,height,
     2 sus,rem,d,zi,iobs,icalc,ires
      common/contg/nrow,ncol,dy,dx,xmin,ymin,xmax,ymax
      data iout/11/,nz/1/,dum/0./,
     1 id1/'magpoly calculated values from m',
     2 'odel                            '/,
     3 id2/'magpoly residual values(observed',
     4 '-calculated)                    '/
      if(isw.eq.1)
     & open(iout,file=file1,status='unknown',form='unformatted')
      if(isw.eq.2)
     & open(iout,file=file2,status='unknown',form='unformatted')
      if(isw.eq.1) write(iout) id1,ncol,nrow,nz,ymin,dx,xmin,dy
      if(isw.eq.2) write(iout) id2,ncol,nrow,nz,ymin,dx,xmin,dy
      is=1
      ie=ncol
      do 10 i=1,nrow
      write(iout) dum,(z(j),j=is,ie)
      is=ie+1
      ie=ie+ncol
   10 continue
      close(iout)
      return
      end
      subroutine scale3(xmin,xmax,n,xminp,xmaxp,dist)
c--see comm. acm, vol.16, no.10, oct 1973, pp.639-640.
      dimension vint(5)
      data vint(1),vint(2),vint(3),vint(4),vint(5)/1.,2.,5.,10.,20./
      if(xmin.lt.xmax.and.n.gt.1) go to 10
      write(6,99999)
99999 format(34h improper input supplied to scale3)
      return
   10 del=.00002
      fn=n
      a=(xmax-xmin)/fn
      al=alog10(a)
      nal=al
      if(a.lt.1.) nal=nal-1
      b=a/10.**nal
      do 20 i=1,3
      if(b.lt.(vint(i)+del)) go to 30
   20 continue
      i=4
   30 dist=vint(i)*10.**nal
      fm1=xmin/dist
      m1=fm1
      if(fm1.lt.0.) m1=m1-1
      if(abs(float(m1)+1.-fm1).lt.del) m1=m1+1
      xminp=dist*float(m1)
      fm2=xmax/dist
      m2=fm2+1.
      if(fm2.lt.(-1.)) m2=m2-1
      if(abs(fm2+1.-float(m2)).lt.del) m2=m2-1
      xmaxp=dist*float(m2)
      np=m2-m1
      if(np.le.n) go to 40
      i=i+1
      go to 30
c
c         next 3 statements comented out by r.godson,usgs,denver,co.
c
c   40 nx=(n-np)/2
c      xminp=xminp-float(nx)*dist
c      xmaxp=xminp+float(n)*dist
40      if(xminp.gt.xmin) xminp=xmin
      if(xmaxp.lt.xmax) xmaxp=xmax
      return
      end
      subroutine namemp(icmd)
c
c     namelist simulator for pc's with no namelist statement
c     this subroutine is program independent but it
c     calls subroutine check which is program dependent
c     The program does not handle the repeat indicator *
c     nor does it handle subscripted array variables or
c     substrings of character variables.
c     It also does not handle complex or logical variables.
c     nn = number of characters in value name(tvar)
c     chv = character variable indicator(logical)
c     inum = an array index number
c     nvar = a number returned by subroutine check to
c     determine whether a variable is an array
c     numa = a number returned by subroutine check
c     which is compared against nvar
c     r.godson,usgs,denver,co., 11/87
c
      parameter(ivar=80)
      character var(ivar),comma,apos,blank,equal,amp,dollar
      character*6 pvar
      character*56 tvar
      logical chv,vset
      data comma/','/,apos/#27/,blank/' '/,equal/'='/
      data amp/'&'/,dollar/'$'/,pvar/' '/,tvar/' '/
c
c     get namelist start name(e.g.&parms)
c
    5 read(icmd,1000,end=910) var
      do 7 i=1,ivar
      if(var(i).ne.amp.and.var(i).ne.dollar) go to 7
      do 6 m=i+1,ivar
      if(var(m).eq.blank) go to 15
    6 continue
    7 continue
      go to 5
c
c     start processing variables
c
   10 read(icmd,1000,end=900) var
 1000 format(80a1)
      m=1
   15 do 20 i=m,ivar
      if(var(i).eq.amp.or.var(i).eq.dollar) go to 900
      if(var(i).ne.blank.and.var(i).ne.comma) go to 30
   20 continue
      go to 10
   30 mm=0
c
c     check for continuation of array values on succeding lines
c
      item=ichar(var(i))
      if(item.eq.43.or.item.eq.45.or.item
     & .eq.46.or.(item.gt.47.and.item.lt.58).or.var(i).eq.apos) then
      k=i
      go to 95
      endif
c
c     get program variable name
c
   35 inum=1
      pvar=blank
      vset=.false.
      do 40 j=i,ivar
      if(var(j).eq.blank.or.var(j).eq.equal) go to 50
      mm=mm+1
      pvar(mm:mm)=var(j)
   40 continue
c
c     get variable value
c
     

   50 do 80 k=j+1,ivar
      if(var(k).eq.blank.or.var(k).eq.equal) go to 80
      if(var(k).eq.apos) then
c
c     character variable
      chv=.true.
      vset=.true.
      nn=0
      do 60 l=k+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
   60 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      else
c
c     not a character variable
      chv=.false.
      nn=0
      do 70 l=k,ivar
      if(var(l).eq.blank.or.var(l).eq.comma
     & .or.var(l).eq.amp.or.var(l).eq.dollar) go to 90
c
c     check for non numeric character
c
      item=ichar(var(l))
      if((item.lt.48.and.(item.ne.46
     & .and.item.ne.45.and.item.ne.43)) .or. (item .gt. 57 .and.
     & (item .ne. 69 .and. item .ne. 101 .and. item.ne.68
     &  .and. item.ne.100))) then
      write(*,*) 'missing first apostropy or a non numeric',
     & ' character in namelist variable ',pvar
      stop
      endif
      nn=nn+1
      tvar(nn:nn)=var(l)
   70 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no delimiter'
      stop
      endif
   80 continue
      write(*,*) 'error in namelist variable ',pvar,' no  value'
      stop
c
c     call program dependent subroutine to assign values
c
   90 call check(pvar,tvar,nn,chv,nvar,numa,inum)
      tvar=blank
      k=l+1
c
c     check to see if variable is an array
c     the variable numa used is program dependent
c
      if(nvar.lt.numa) go to 110
   95 if(.not.chv) then
c
c     array variable
c
      nn=0
      do 100 l=k,ivar
      if(var(l).eq.blank.and.nn.eq.0) go to 100
      if(var(l-1).eq.blank.and.var(l).eq.comma.and.nn.eq.0) go to 100
      if((var(l).eq.amp.or.var(l).eq.dollar).and.nn.eq.0) go to 900
c
c     check for consecutive commas
c
      if(var(l).eq.comma.and.vset) then
      vset=.false.
      inum=inum+1
      go to 90
      else if(var(l).eq.comma) then
      inum=inum + 1
      go to 100
      endif
c
      if(var(l).eq.blank
     & .or.var(l).eq.amp.or.var(l).eq.dollar) then
      inum=inum+1
      go to 90
      endif
      if(ichar(var(l)).gt.57.and.nn.eq.0) go to 120
      nn=nn+1
      tvar(nn:nn)=var(l)
      vset=.true.
  100 continue
      else
c
c     character array variable
c
      do 105 j=k,ivar
      if(var(j).eq.blank) go to 105
      if(var(j).eq.amp.or.var(j).eq.dollar) go to 900
      if(ichar(var(j)).gt.57) then
      l=j
      go to 120
      endif
c
c     check for consecutive commas
c
      if(var(j).eq.comma.and.vset) then
      vset=.false.
      go to 105
      else if(var(j).eq.comma) then
      inum=inum + 1
      go to 105
      endif
      if(var(j).eq.apos) then
      nn=0
      inum=inum+1
      vset=.true.
      do 103 l=j+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
  103 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      endif
  105 continue
c
c     end of array values
c
      endif
  110 m=l+1
      if(var(l).eq.amp.or.var(l).eq.dollar) go to 900
      go to 130
  120 m=l
  130 if(m.lt.81) then
      go to 15
      else
      go to 10
      endif
  900 return
  910 write(*,*) ' error in namelist;no beginning & or $'
      stop
      end
c*********
      subroutine check(pvar,tvar,nn,chv,nvar,numa,inum)
c
c     assigns values to proper variable
c     variables are passed to program minc through common blocks
c     numr=position in the array var where real variables start
c     numa=position in the array var where arrays start
c     nnvar=number of variables in program minc
c
      parameter(nnvar=21,numr=8)
      character*6 pvar,var(nnvar)
      character*56 tvar,kvar,cfmt
      logical chv
      character*4 unix,uniz
      character*8 name
      character*40 idplot
      character*50 ifile,ifile2,rfile,cafile,pfile
      character*80 ifmt
      common/st/cafile,rfile,pfile,ifile,ifile2,name,
     & ifmt,idplot,unix,uniz
      common/maginp/iplotr,ibody,lsqs,
     & dec,fincl,datum,xscale,dc,height,
     & sus,rem,d,zi,iobs,icalc,ires
      common/vec/cde,cie,sde,sie,ex,ey,ez,earth,cdci,sdci
      common/widect/naxcol
      data var/'iplotr','ibody','iobs','icalc','ires','naxcol',
     & 'lsqs','datum','xscale','dc','height','ifile','ifile2',
     & 'unix','uniz','name','ifmt','idplot','earth','dec',
     & 'fincl'/
      inum=1
      numa=22
      do 190 i=1,nnvar
      if(pvar.ne.var(i)) go to 190
      if(.not.chv) then
c
c     noncharacter value
c     right justify then number in variable kvar
c
      m=57-nn
      im=m-1
      kvar(m:56)=tvar(1:nn)
      if(i.lt.numr) then
c
c     integer value
c
      write(cfmt,50) im,nn
   50 format('(',i2,'x,i',i2,')')
      read(kvar,cfmt) jvar
      else
c
c     real value
c
      write(cfmt,60) im,nn
   60 format('(',i2,'x,g',i2,'.0)')
      read(kvar,cfmt) xvar
      endif
      endif
      go to (101,102,103,104,105,106,107,108,109,110,111,112,
     & 113,114,115,116,117,118,119,120,121),i
  101 iplotr=jvar
      go to 200
  102 ibody=jvar
      go to 200
  103 iobs=jvar
      go to 200
  104 icalc=jvar
      go to 200
  105 ires=jvar
      go to 200
  106 naxcol=jvar
      go to 200
  107 lsqs=jvar
      go to 200
  108 datum=xvar
      go to 200
  109 xscale=xvar
      go to 200
  110 dc=xvar
      go to 200
  111 height=xvar
      go to 200
  112 ifile=tvar(1:nn)
      go to 200
  113 ifile2=tvar(1:nn)
      go to 200
  114 unix=tvar(1:nn)
      go to 200
  115 uniz=tvar(1:nn)
      go to 200
  116 name=tvar(1:nn)
      go to 200
  117 ifmt=tvar(1:nn)
      go to 200
  118 idplot=tvar(1:nn)
      go to 200
  119 earth=xvar
      go to 200
  120 dec=xvar
      go to 200
  121 fincl=xvar
      go to 200
  190 continue
      write(*,*) ' error in namelist - ',pvar,' variable not included'
      stop
  200 nvar=i
      return
      end
