c***************
c program to input exact talwani polygonal prisms for 3d gravity-plouff
c modified by r.godson for denver multics 5/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,idplot,
     1 unix,uniz
      common/grvin/iplotr,ibody,lsqs,
     1 datum,xscale,dc,height,
     2 iobs,icalc,ires
      common /field/r(1000)
      common/contg/nrow,ncol,dy,dx,yo,xo,xma,yma
      common /fpts/ xf(1001),yf(1001),zf(1000)
      common /bodk/ mm,nlim,jbody,short,nend
      common /trans/ inpt,iout,jin,iprint
      dimension a(1000),b(1000),dstnc(900)
c dimension limits for 3 parameters.  corners, bodies, fieldpoints.
c one more than dimension of total number of fieldpoints
      data in1/10/,in2/11/,iout1/12/,iout2/13/,nz/1/,dum/0./,
     1 id1/'gravpoly calculated values from ',
     2 'model                           '/,
     3 id2/'gravpoly residual values(observe',
     4 'd-calculated                    '/
      data ncorn,nbod/900,50/
      data ddval/1.e30/
c
c     get input parameters
c
      call gravin(*220)
   10 write (iout,20) idplot
   20 format (' gravity effect of polygonal prisms',/,
     1 1x,                            ' for ',5x,a40)
      facx=0.0
      facz=0.0
      call unit(unix,uniz,facx,facz)
      if (facx .eq. 0.0) go to 200
      zfac=1.0/facz
      gfac=6.670*facx
   30 short=1.0e20
c
c     read body coordinates
c
      call bodg(facz,nbod,dstnc,ncorn,gfac)
c mm (common) is total number of layers
      if (mm .eq. 0) go to 200
c
c     read fieldpoint locations and maybe values
c
      call ptsm(nk,*220)
      title='body corners'
      if(ibody.ne.0) call gr3dpt(nk,*220)
      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 gravity 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 220
      write(6,240) i,nrow
      nk=0
      call rowin(xf,mcol,in1)
      if(ifile2.ne.' ') call rowin(yf,mcol,in2)
      do 32 j=2,mcol
      if(lsqs.ne.0) then
              if(xf(j).ge.ddval.or.yf(j).ge.ddval) stop 'dvals found'
      endif
      nk=nk+1
      r(nk)=xf(j)
      if(ifile2.ne.' ') height=yf(j)
      zf(nk)=height
      yf(nk)=xo1+(j*dx)
      xf(nk)=yo1+(i*dy)
c      write(iout,30) yf(k),xf(k),zf(k),r(k)
   32 continue
   33 if (nk .lt. 1) go to 200
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 gr3dpt(nk,*220)
c
c     check if plot of observed values desired
c
      title='observed values'
c     if(iobs.ne.0.and.iplotr.eq.9.and.name.eq.'gridded')
c    1 call cplot(dc,r,nk,1,title,*365)
      do 40 k=1,nk
   40 zf(k)=zf(k)*facz
   50 call gtalw(nk,mm,short,dstnc,a)
      if(datum .ne. 0.) write (iout,60) datum
   60 format (/,1x,f10.2,' mgal is added to gravity values')
      if (lsqs .eq. 0) go to 130
      do 70 k=1,nk
   70 r(k)=r(k)-datum
      if (lsqs .ne. 7) go to 90
      ier=0
      call glsqs(nk,a,ier,datum)
      if (ier .ne. 0) go to 200
   90 std=0.0
      avg=0.0
      do 100 k=1,nk
      errk=r(k)-a(k)
      b(k)=errk
      avg=avg+errk
  100 std=std+errk*errk
  110 if (nk  .lt. 2) go to 130
      std=sqrt(std/(nk -1))
      avg=avg/nk
      write (iout,120) std,avg
  120 format (///,5x,'standard deviation=',
     1 f8.2,' and average difference=',f8.2,' milligals')
  130 write (iout,140)
  140 format (/,15x,'x',9x,'y',9x,'z',9x,'g',4x,'residual')
      do 180 k=1,nk
      x=xf(k)
      y=yf(k)
      zfk=zf(k)*zfac
      grav=a(k)+datum
      a(k)=grav
      if (lsqs .eq. 0) go to 160
      write (iout,150) k,y,x,zfk,  grav,b(k)
  150 format (i6,3x,3f10.3,2f9.2)
      go to 180
  160 write (iout,170) k,y,x,zfk,  grav
  170 format (i6,3x,3f10.3,f9.2)
  180 continue
      if(name.ne.'byrow') go to 340
      write(iout1) dum,(a(k),k=1,nk)
      if(lsqs.eq.1) write(iout2) dum,(b(k),k=1,nk)
      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,nk,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,nk,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)
      go to 220
  365 write(iprint,375)
  375 format(' contour interval(dc) is 0.')
  200 write(iprint,210)
  210 format(' error - program abort')
  220 close(jin)
      close(iout)
      write(iprint,230) pfile
  230 format(/,' output is in file ',a50)
  240 format(1h+,'Working on row ',i5,' of ',i5)
      stop
      end
c***************
      subroutine gravin(*)
c
c     namelist input parameter subroutine
c
c     r.godson - 3/81
c
      character*4 unix,uniz
      character*8 name
      character*40 idplot
      character*50 ifile,ifile2,cfile
      character*50 cafile,rfile,pfile
      character*80 ifmt
      common/st/cafile,rfile,pfile,ifile,ifile2,name,ifmt,idplot,
     1 unix,uniz
      common/grvin/iplotr,ibody,lsqs,
     1 datum,xscale,dc,height,
     2 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
      inpt=5
      iout=8
      iprint=6
      jin=9
      iplotr=9
      naxcol=130
      ibody=0
      iobs=0
      icalc=0
      ires=0
      idplot=' '
      ifile=' '
      ifile2=' '
      lsqs=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 namegp(jin)
      if(lsqs.eq.0) ires=0
      if(ifile.eq.' ') go to 50
c
c     determine output file names
      call names(cfile,pfile,cafile,rfile)
c
      open(iout,file=pfile,status='unknown')
      return
   50 write(iprint,60)
   60 format(' no observed data file name specified')
   70 return 1
      end
c***************
      subroutine names(cfile,pfile,cafile,rfile)
      character*(*) cfile,pfile,cafile,rfile
      ic=index(cfile,'.')-1
      pfile=cfile(:ic)//'.prt'
      cafile=cfile(:ic)//'.cal'
      rfile=cfile(:ic)//'.res'
      return
      end
c***************
      subroutine bodg(facz,nbod,dstnc,ncorn,gfac)
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
      dimension dstnc(1)
      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,idplot,
     1 unix,uniz
      common/grvin/iplotr,ibody,lsqs,
     1 datum,xscale,dc,height,
     2 iobs,icalc,ires
      common /bod/ rho(50),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 /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=180)num,ifprnt(mm+1),dencon,top,bottom
      if(lsqs.eq.7) dencon=1.
      nst=nend+1
      nend=nend+num
      mm=mm+1
      if(mm.gt.nbod) go to 130
      if(nend.gt.ncorn) go to 70
      ll(mm)=num
      zee(mm)=top*facz
      zm(mm)=bottom*facz
      rho(mm)=dencon*gfac
      write(iout,30) mm,top,bottom,num,dencon
   30 format (' body',i3,' between depths',f11.3,' and',f11.3,
     1 ' and has',i4,' vertices. density = ',f8.3)
   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 170
   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 170
c coordinates of corners of body east,north (x,y) pairs.
c body has at least 3 corners
   90 read(jin,*,end=150) (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
  110 np=nend+1
      xb(np)=xb(nst)
      yb(np)=yb(nst)
      do 120 j=nst,nend
      dist=(xb(j+1)-xb(j))**2+(yb(j+1)-yb(j))**2
      if (dist .eq. 0.0) go to 120
      dist=sqrt(dist)
      if (dist .lt. short) short=dist
  120 dstnc(j)=dist
      go to 10
  130 write (iprint,140) mm
  140 format (5x,'number of bodies=',i3,' exceeds limitation.')
      go to 170
  150 write(iprint,160) mm
  160 format(' end of file reading x,y pairs for body no. ',i3)
  170 mm=0
  180 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
      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 gtalw(nk,mm,short,dstnc,a)
c exact solution for talwani 3d gravity in single precision (7 figures)
      common /side/xb(900),yb(900), ll(50)
      common /bod/ rho(50),ifprnt(50),zee(50),zm(50)
      common /fpts/ xp(1001),yp(1001),zp(1000)
      common /trans/ inpt,iout,jdisk,page
      dimension a(1),dstnc(1)
c 1/1000 length of shortest non-zero side
      test=short*1.0e-3
      do 240 k=1,nk
      x=xp(k)
      y=yp(k)
      z=zp(k)
      grav=0.0
      ne  =0
c loop of bodies
      do 230 m=1,mm
c number of vertices
      ns =ne  +1
      lk= ll(m)
      z1=zee(m)-z
      z2=zm (m)-z
      ne  =ns +lk-1
      gf=rho(m)
      jfprnt=ifprnt(m)
      g=0.0
c to give clockwise positive
      h1=abs(z1)
      h2=abs(z2)
      zz=z1+z2
      if (zz)10,230,20
c case where most of mass is above fieldpoint (minus g)
   10 gf=-gf
   20 if (h1 .lt. h2) go to 30
      t=h2
      h2=h1
      h1=t
   30 x1=xb(ns)-x
      y1=yb(ns)-y
c last point is repeat of first
      xf=x1
      yf=y1
      r1s=x1*x1+y1*y1
      r1sf=r1s
      hh1=h1*h1
      hh2=h2*h2
      dh=h2-h1
      r11= sqrt(r1s+hh1)
      r12= sqrt(r1s+hh2)
      r11f=r11
      r12f=r12
      ang=0.0
      t1=0.0
      t2=0.0
      w1=0.0
      w2=0.0
      w3=0.0
      n=ns
c loop of all sides
      do 210 l=1,lk
      nm=n
      if (l .ne. lk) go to 40
      x2=xf
      y2=yf
      r2s=r1sf
      r21=r11f
      r22=r12f
      go to 50
   40 n=n+1
      x2=xb(n)-x
      y2=yb(n)-y
      r2s=x2*x2+y2*y2
      r22= sqrt(r2s+hh2)
      r21= sqrt(r2s+hh1)
   50 ds=dstnc(nm)
      if (test .gt. ds) go to 190
c     same sign as definition used for magnetics. opposite last gravity
      p=(x1*y2-x2*y1)/ds
      dx=x2-x1
      dy=y2-y1
c test if perpendicular distance to line segment is less than 1/1000
c of shortest segment length
      ap=abs(p)
      if (ap .lt. test) go to 190
      d1= (x1*dx+y1*dy)/ds
      d2= (x2*dx+y2*dy)/ds
      dd1= abs(d1)
      dd2= abs(d2)
      if (dd1 .lt.  test) d1=0.0
      if (dd2 .lt.  test) d2=0.0
      fnum=ap*ds
      den=x1*x2+y1*y2
      if (den) 60,100,110
   60 gl=3.141593
      if (p .gt. 0.0) go to 80
      gl=-gl
   70 den=-den
   80 call tandet(fnum/den,w3,ang)
   90 ang=ang+gl
      go to 120
  100 gl=1.570796
      if (p .lt. 0.0) gl=-1.570796
      go to 90
  110 gl=0.0
      if (p .gt. 0.0) go to 80
      go to 70
  120 if (d1 .eq. 0.0) go to 130
      call tandet(-d1*h1/(p*r11),w1,t1)
      call tandet( d1*h2/(p*r12),w2,t2)
  130 if (d2 .eq. 0.0) go to 140
      call tandet( d2*h1/(p*r21),w1,t1)
      call tandet(-d2*h2/(p*r22),w2,t2)
  140 ph=p*p+hh1
      rz1=r11+d1
      if (d1 .ge. 0.0) go to 160
      if (d2 .gt. 0.0) go to 150
      gl=(r21-d2)/(r11-d1)
      go to 180
  150 dtest=ph/(d1*d1)
      if (dtest .lt. 0.01) rz1=dd1*0.5*dtest*(1.0-0.25*dtest)
  160 rz2=r21+d2
      if (d2 .ge. 0.0) go to 170
      dtest=ph/(d2*d2)
      if (dtest .lt. 0.01) rz2=dd2*0.5*dtest*(1.0-0.25*dtest)
  170 gl=rz1/rz2
  180 gl=p*alog(gl*(r22+d2)/(r12+d1))
      g=g-gl
  190 if (jfprnt .lt. 0) write (iout,200) l,p,d1,d2,ang,gl,w1,w2
  200 format (i4,7e16.8)
      x1=x2
      y1=y2
      r1s=r2s
      r12=r22
      r11=r21
  210 continue
      t1=t1+ atan(w1)
      t2=t2+ atan(w2)
      ang=  ang+atan(w3)
      g=gf*(ang*dh+g+h1*t1+h2*t2)
      if (jfprnt .ne. 0) write (iout,220) k,g,ang ,t1,t2
  220 format (i4,f8.2,' mgal',15x,3f9.3)
  230 grav=grav+g
  240 a(k)=grav
      return
      end
c***************
      subroutine tandet(f,zze,sume)
c program for arctangent addition. negative f for subtraction.
      q1=zze+f
      q2=1.0d0-zze*f
      if (q2 .ne. 0.0  ) go to 40
      if (q1) 10,20,20
   10 sume=sume-1.570796
      go to 30
   20 sume=sume+1.570796
   30 zze=0.0
      return
   40 zze=q1/q2
      if (q2 .gt. 0.0  ) return
      if (q1) 50,60,70
   50 sume=sume-3.141593
   60 return
   70 sume=sume+3.141593
      return
      end
c***************
      subroutine glsqs(kk,a,ier,datum)
c least-squares calculation of density regression
      double precision fn,sx,sy,sxx,sxy,syy
      integer nm(4)/'n','x','s','y'/
      common /field/r(1000)
      common /trans/ inpt,iout,jdisk,page
      dimension a(1)
      l=1
      if (kk .lt. 3) go to 40
      fn=kk
      sx =0.0d0
      sy =0.0d0
      syy=0.0d0
      sxy=0.0d0
      sxx=0.0d0
      do 10 k=1,kk
      rk=r(k)
      ak=a(k)
      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=2
      if (xx .le. 0.0) go to 40
      xy=sxy-sx*sy/fn
      yy=syy-sy*sy/fn
      tm=xy/xx
      tb=(sy-tm*sx)/fn
      l=4
      if (yy .le. 0.0) go to 40
      coer=xy/sqrt(xx*yy)
      stdm=yy*(1.0-coer*coer)/(xx*(kk-2))
      l=3
      if (stdm .lt. 0.0) go to 40
      stdm=sqrt(stdm)
      tb=tb+datum
      write (iout,20) tm,stdm,coer,tb
   20 format (//,' least-squares density is',f7.3,' (std',f6.3,') with',
     1 ' correlation coefficient=',           f6.2,3x,'datum=',f7.2,
     2 ' mgal.')
      do 30 k=1,kk
   30 a(k)=tm*a(k)
      return
   40 ier=1
      write (iout,50) nm(l)
   50 format ('0stop. density can be calculated. cutoff ',a1)
      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,idplot,
     1 unix,uniz
      common/grvin/iplotr,ibody,lsqs,
     1 datum,xscale,dc,height,
     2 iobs,icalc,ires
      common/contg/nrow,ncol,dy,dx,yo,xo,xma,yma
      common/field/r(1000)
      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,'g',/)
      if(name.eq.'gridded') go to 50
c
c     user formatted data
c
      open(in,file=ifile,status='old')
      i=1
      nz=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 the 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 gr3dpt(npts,*)
c
c     plotting routine for body corners & fieldpoints
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
      integer sym(50),iddplot(10)
      equivalence(idplot,iddplot(1))
      common/mplot/xmin,ymin,xmax,ymax
      common/fpts/xf(1001),yf(1001),zf(1000)
      common/st/cafile,rfile,pfile,ifile,ifile2,name,ifmt,idplot,
     1 unix,uniz
      common/grvin/iplotr,ibody,lsqs,
     1 datum,xscale,dc,height,
     2 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
      character*1 icross
      dimension xd(2),yd(2),xs(4),ys(4),jdim(2),
     1 xxb(900),yyb(900)
      equivalence(icross,icr)
      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
      xxb(m)=xb(k)
      yyb(m)=yb(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*40 idplot
      character*32 id1(2),id2(2)
      character*50 file1,file2,file3,ifile,ifile2
      character*80 ifmt
      dimension z(1)
      common/st/file1,file2,file3,ifile,ifile2,name,ifmt,idplot,
     1 unix,uniz
      common/grvin/iplotr,ibody,lsqs,
     1 datum,xscale,dc,height,
     2 iobs,icalc,ires
      common/contg/nrow,ncol,dy,dx,xmin,ymin,xmax,ymax
      data iout/11/,nz/1/,dum/0./,
     1 id1/'gravpoly calculated values from ',
     2 'model                           '/,
     3 id2/'gravpoly residual values(observe',
     4 'd-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 namegp(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=18,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/grvin/iplotr,ibody,lsqs,
     & datum,xscale,dc,height,
     & iobs,icalc,ires
      common/widect/naxcol
      data var/'iplotr','ibody','iobs','icalc','ires','naxcol',
     & 'lsqs','datum','xscale','dc','height','ifile','ifile2',
     & 'unix','uniz','name','ifmt','idplot'/
      inum=1
      numa=19
      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),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
  190 continue
      write(*,*) ' error in namelist - ',pvar,' variable not included'
      stop
  200 nvar=i
      return
      end

v

