
c PFS2.FOR
c  Potential-field sounding, 2-D.  To solve for point-mass representation
c  of gravity profile using iterative Euler method  -- point masses to be
c  converted to density grid by associated program BUBBLE.  Input is gravity
c  and horizontal and vertical derivatives (3 files, in 1-row-grid format).
c  Output is and ascii x, z, Mass(x,z) file.  Z is made positive-down for
c  compatibility with the standard-file grid representing the density cross-
c  section (as output from program BUBBLE) having its origin in the lower
c  left corner.
c  Lin Cordell 17 September 1992.
c
      character*50 ifile,ofile
      dimension g(1000),gx(1000),gz(1000),id(14),pgm(2)
      data dval/0.170412e+39/
      dval2=1.e20
c  Read input files and parameters.
      write(6,1)
1     format(' Enter input g(x) file name:'/' *'$)
      read(5,100)ifile
100   format(a50)
      open(10,file=ifile,form='unformatted',status='old')
      read(10)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      if(ncol.gt.1000.or.nz.ne.1.or.dx.eq.0.e0) go to 900
      write(6,2)
2     format(' Enter input gx(x) file name:'/' *'$)
      read(5,100)ifile
      open(11,file=ifile,form='unformatted',status='old')
      read(11)id,pgm,ncolx,nrow,nz,x0x,dxx,y0,dy
      if(ncolx.ne.ncol.or.x0x.ne.x0.or.dx.ne.dxx)then
        print*,'Profile specs ncol, x0, dx must be compatible.'
        stop
      endif
      write(6,3)
3     format(' Enter input gz(x) file name:'/' *'$)
      read(5,100)ifile
      open(12,file=ifile,form='unformatted',status='old')
      read(12)id,pgm,ncolx,nrow,nz,x0x,dxx,y0,dy
      if(ncolx.ne.ncol.or.x0x.ne.x0.or.dx.ne.dxx)then
        print*,'Profile specs ncol, x0, dx must be compatible.'
        stop
      endif
      write(6,201)
201   format(1x,'Enter output file name:'/' *'$)
      read(5,100)ofile
      open(13,form='formatted',status='scratch')
      write(6,4)
4     format(1x,'Enter epsilon:'/' *'$)
      read*,eps
      call rowio(g,y,1,ncol,10,1)
      call rowio(gx,y,1,ncol,11,1)
      call rowio(gz,y,1,ncol,12,1)
c  0th source a dummy, note alpha=0.
      xi=x0+.5*ncol*dx
      zeta=xi
      zetasq=zeta**2
      alpha=0
c  Iterations.
50    continue
      gmax=-1.e20
      gmin=1.e20
c  Remove effect of last source, and find gmax.
      do 40 i=2,ncol-1
      x=x0+(i-1)*dx
      xx=xi-x
      b=xx*xx+zetasq
      bsq=b*b
      calc=((alpha*zeta)/b)
      g(i)=g(i)-calc
      calc=(2.*xx*zeta)/bsq
      gx(i)=gx(i)-calc
      calc=((2.*zetasq)/bsq)-(1./b)
      gz(i)=gz(i)-calc
      if(g(i).gt.gmax) then
        gmax=g(i)
        imax=i
      endif
      if(g(i).lt.gmin) then
        gmin=g(i)
        imin=i
      endif
40    continue
      ix=imax
      if(abs(gmin).gt.abs(gmax))ix=imin
      xgmax=x0+(ix-1)*dx
      print*,'Gmax =',g(ix),', at x =',xgmax
c  Quit if |Gmax| < epsilon.
      if(abs(g(ix)).lt.eps)go to 998
c  Solve for x position of source by three-point interpolating to
c  dg/dx = 0; solve for source depth by Euler, and then mass term
c  by from Gmax.
      xio=(dx*(g(ix-1)-g(ix+1)))/(2.*(g(ix+1)+g(ix-1)-2.*g(ix)))
      s1=0.
      s2=0.
      s3=0.
      do 10 i=-1,1
      s1=s1+g(ix+i)*gz(ix+i)
      s2=s2+(xio-i*dx)*gx(ix+i)*gz(ix+i)
10    s3=s3+gz(ix+i)*gz(ix+i)
      zeta=(s1-s2)/s3
      zetasq=zeta*zeta
      alpha=(g(ix)*(xio*xio+zetasq))/zeta
      xi=x0+(ix-1)*dx+xio
      if(zeta.eq.0.)zeta=.1*dx
      zetaout=-zeta
      write(13,30)xi,zetaout,alpha
30    format(3e16.6)
      go to 50
      go to 998
900   write(6,901)
901   format(' Can''t handle it.')
998   continue
      rewind 13
c  Convert alpha=2*gamma*Mass to mass, change sign of z and mass if z<0,
c  and switch sign of z for grid convention (y0 at base of grid).
      cnst=2.*6.67
      rewind 13
      open(14,file=ofile,form='formatted',status='new')
200   read(13,30,end=202)xi,zeta,alpha
      if(zeta.gt.0)then
        zeta=-zeta
        alpha=-alpha
      endif
      alpha=alpha/cnst
      write(14,30)xi,zeta,alpha
      go to 200
202   continue
999   continue
      stop
      end
      subroutine rowio(a,y,nz,ncol,ld,key)
      dimension a(nz,ncol)
      go to (1,2),key
1     read(ld) y,a
      go to 90
2     write(ld)y,a
90    return
      end
