
c PFS3.FOR
c  Potential-field sounding, 3D.  To solve for point-mass representation
c  of gravity grid using iterative Euler method  -- point masses to be
c  converted to 3D density grid by associated program BUBBLE3.  Input is
c  gravity g and horizontal-derivative gx,gy and vertical-derivative gz
c  files (4 files, in grid format).
c  Output is an ascii x,y,z,mass file.  Z is positive-down.
c
c  Lin Cordell 17 September 1992, 7 January 1993
c
      character*50 ifile,ofile
      dimension g(10000),gx(10000),gy(10000),gz(10000),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,y) 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
      ndim=ncol*nrow
      if(ndim.gt.10000)go to 900
      write(6,2)
2     format(' Enter input gx(x,y) file name:'/' *'$)
      read(5,100)ifile
      open(11,file=ifile,form='unformatted',status='old')
      read(11)id,pgm,ncol2,nrow2,nz,x02,dx2,y02,dy2
      if(ncol2.ne.ncol.or.x02.ne.x0.or.dx2.ne.dx.or.
     & nrow2.ne.nrow.or.y02.ne.y0.or.dy2.ne.dy)then
        print*,'Profile specs must be compatible.'
        stop
      endif
      write(6,3)
3     format(' Enter input gy(x,y) file name:'/' *'$)
      read(5,100)ifile
      open(12,file=ifile,form='unformatted',status='old')
      read(12)id,pgm,ncol2,nrow2,nz,x02,dx2,y02,dy2
      if(ncol2.ne.ncol.or.x02.ne.x0.or.dx2.ne.dx.or.
     & nrow2.ne.nrow.or.y02.ne.y0.or.dy2.ne.dy)then
        print*,'Profile specs ncol, x0, dx must be compatible.'
        stop
      endif
      write(6,4)
4     format(' Enter input gz(x,y) file name:'/' *'$)
      read(5,100)ifile
      open(13,file=ifile,form='unformatted',status='old')
      read(13)id,pgm,ncol2,nrow2,nz,x02,dx2,y02,dy2
      if(ncol2.ne.ncol.or.x02.ne.x0.or.dx2.ne.dx.or.
     & nrow2.ne.nrow.or.y02.ne.y0.or.dy2.ne.dy)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(14,file=ofile,form='formatted',status='new')
      write(6,5)
5     format(1x,'Enter epsilon:'/' *'$)
      read*,eps
c  Read data arrays, locate max|g|
      gmax=-1.e20
      gmin=1.e20
      do 7 j=1,nrow
      jj=(j-1)*ncol
      k1=jj+1
      k2=jj+ncol
      read(10)y,(g(k),k=k1,k2)
      do 6 i=1,ncol
      k=jj+i
      if(g(k).gt.gmax)then
        gmax=g(k)
        igmax=i
        jgmax=j
      endif
      if(g(k).lt.gmin)then
        gmin=g(k)
        igmin=i
        jgmin=j
      endif
6     continue
      read(11)y,(gx(k),k=k1,k2)
      read(12)y,(gy(k),k=k1,k2)
      read(13)y,(gz(k),k=k1,k2)
7     continue
c  Solve for source at (imax,jmax), where abs(g) is max.
8     imax=igmax
      jmax=jgmax
      if(abs(gmin).gt.abs(gmax))then
        imax=igmin
        jmax=jgmin
      endif
c  Move in, if max |g| is on an edge.
      if(imax.eq.1)imax=2
      if(imax.eq.ncol)imax=ncol-1
      if(jmax.eq.1)jmax=2
      if(jmax.eq.nrow)jmax=nrow-1
      k=(jmax-1)*ncol+imax
      gmax=abs(g(k))
      print*,'|Gmax| =',gmax
      if(abs(g(k)).lt.eps)go to 30
      gxp=g(k+1)
      if(abs(gxp).gt.gmax)gxp=g(k)
      gxm=g(k-1)
      if(abs(gxm).gt.gmax)gxm=g(k)
      denom=gxp+gxm-2.*g(k)
      xp=0.
      if(denom.ne.0.)xp=((gxm-gxp)*dx)/denom
      xi=x0+(imax-1)*dx+xp
      gyp=g(k+ncol)
      if(abs(gyp).gt.gmax)gyp=g(k)
      gym=g(k-ncol)
      if(abs(gym).gt.gmax)gym=g(k)
      denom=gyp+gym-2.*g(k)
      yp=0.
      if(denom.ne.0.)yp=((gym-gyp)*dx)/denom
      eta=y0+(jmax-1)*dy+yp
c   Solution for zeta, 5 points only.
      s1=gxm*gz(k-1)+2.*g(k)*gz(k)+gxp*gz(k+1)
     & +gym*gz(k-ncol)+gyp*gz(k+ncol)
      s2=(xp+dx)*gx(k-1)*gz(k-1)+2.*xp*gx(k)*gz(k)+(xp-dx)*gx(k+1)
     & *gz(k+1)+xp*(gx(k-ncol)*gz(k-ncol)+gx(k+ncol)*gz(k+ncol))
      s3=yp*(gy(k-1)*gz(k-1)+gy(k+1)*gz(k+1))+(yp+dy)*gy(k-ncol)
     & *gz(k-ncol)+2.*yp*gy(k)*gz(k)+(yp-dy)*gy(k+ncol)*gz(k+ncol)
      s4=gz(k-1)*gz(k-1)+2.*gz(k)*gz(k)+gz(k+1)*gz(k+1)
     & +gz(k-ncol)*gz(k-ncol)+gz(k+ncol)*gz(k+ncol)
      zeta=(2.*s1-s2-s3)/s4
      if(zeta.le.0.)zeta=.1*dx
      zetasq=zeta*zeta
c   Solution for alpha=gamma*mass.
      alpha=(g(k)*((xp*xp+yp*yp+zetasq)**1.5))/zeta
c   Output this source.
      rmass=alpha/6.67
      write(14,23)xi,eta,zeta,rmass
23    format(4e16.6)
c  Remove effect of this source from g, gx, gy, and gz, locate |gmax|,
c  and, if |gmax| > eps,  recycle.
      gmax=-1.e20
      gmin=1.e20
      do 21 j=1,nrow
      y=y0+(j-1)*dy
      yy=eta-y
      do 20 i=1,ncol
      x=x0+(i-1)*dx
      k=(j-1)*ncol+i
      xx=xi-x
      b=xx*xx+yy*yy+zetasq
      b3=b*sqrt(b)
      b5=b3*b
      g(k)=g(k)-((alpha*zeta)/b3)
      if(g(k).gt.gmax)then
        if(j.gt.1.and.j.lt.nrow-1.and.i.gt.1.and.i.lt.ncol-1)then
          gmax=g(k)
          igmax=i
          jgmax=j
        endif
      endif
      if(g(k).lt.gmin)then
        if(j.gt.1.and.j.lt.nrow-1.and.i.gt.1.and.i.lt.ncol-1)then
          gmin=g(k)
          igmin=i
          jgmin=j
        endif
      endif
      b0=(3.*alpha*zeta)/b5
      gx(k)=gx(k)-(xx*b0)
      gy(k)=gy(k)-(yy*b0)
      gz(k)=gz(k)-(alpha*(((3.*zetasq)/b5)-1./b3))
20    continue
21    continue
      go to 8
c  Iterations completed.
30    continue
      go to 999
900   write(6,901)
901   format(' Can''t handle it.')
999   continue
      stop
      end
