c  phony.for
c
c  Part of the equivalent-source gridding operation.  This
c  program creates a set of phony stations from the part of
c  the grid plugged using megaplug.  The phony and actual data
c  are to be recycled thru the es program to obtain a set of
c  sources sufficient to produce a grid without recourse to
c  plugging.  Input is the plugged and the unplugged es grids, and
c  maybe a topo grid. Output is a post file. Elevation associated
c  with the phony data is either obtained from an input elevation
c  grid (preferable), or input average elevation.
c
c  Lin Cordell 1 Mar 91.
      common x0,dx,y0,dy,hcnst
      character*50 pfile,ufile,tfile,ofile
      dimension p(10000),u(10000),t(10000)
      dimension id(14),pgm(2)
      print*,'Enter plugged grid file name:'
      print 1
1     format(' *',$)
      read(5,100)pfile
100   format(a)
      print*,'Enter unplugged grid file name:'
      print 1
      read(5,100)ufile
      open(10,file=ufile,form='unformatted',status='old')
      read(10)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      if(ncol.gt.10000.or.nz.ne.1.or.dx.eq.0.e0)then
        print*,'Cant handle it.'
        stop
      endif
      open(11,file=pfile,form='unformatted',status='old')
      read(11)id,pgm,ncol2,nrow2,nz,x02,dx2,y02,dy2
      ier=0
      if(ncol2.ne.ncol.or.nrow2.ne.nrow)ier=1
      if(x02.ne.x0.or.y02.ne.y0)ier=1
      if(dx2.ne.dx.or.dy2.ne.dy)ier=1
      if(ier.eq.1)then
        print*,'Cant handle it.'
        stop
      endif
2     print*,'Topo option: enter 1 for input grid:'
      print*,'             enter 2 for input constant:'
      print 1
      read*, iopt
      hcnst=0
      if(iopt.eq.1)then
        print*,'Enter topo grid file name:'
        read(5,100)tfile
        open(13,file=tfile,form='unformatted',status='old')
        read(13)id,pgm,ncol2,nrow2,nz,x02,dx2,y02,dy2
        ier=0
        if(ncol2.ne.ncol.or.nrow2.ne.nrow)ier=1
        if(x02.ne.x0.or.y02.ne.y0)ier=1
        if(dx2.ne.dx.or.dy2.ne.dy)ier=1
        if(ier.eq.1)then
          print*,'Cant handle it.'
          stop
        endif
      else
        if(iopt.ne.2)go to 2
        print*,'Enter constant elevation, positive up, '
        print*,'in the same units as the real data points:'
        print 1
        read*,hcnst
      endif
      print*,'Enter output file name:'
      print 101
101   format(' *'$)
      read(5,100)ofile
      open(12,file=ofile,form='unformatted',status='new')
      print*,'POST FIELD   GRAVITY             MAGNETICS'
      print*,'    1        free-air anomaly    total field residual'
      print*,'    2        Bouguer  anomaly    total field'
      print*,'    3        elevation           terrain clearance'
      print*,'    4        inner zone TC       barometric altitude'
      print*,'    5        outer zone TC       fiducial number'
      print*,'    6        observed gravity    julien date'
      print*,'Enter post fields for data, elevation:'
      print 101
      read(5,*) ifld, ielv
      call work(p,u,t,ncol,nrow,iopt,ifld,ielv)
      stop
      end
      subroutine work(p,u,t,ncol,nrow,iopt,ifld,ielv)
      character*8 ident
      dimension p(ncol),u(ncol),t(ncol),v(6)
      common x0,dx,y0,dy,hcnst
      icount=0
      ident='        '
      v=0
      do 200 j=1,nrow
      read(10)y,u
      read(11)y,p
      if(iopt.eq.1)read(13)y1,t
      y=y0+(j-1)*dy
      do 210 i=1,ncol
      if(u(i).lt.1.e20)go to 210
      icount=icount+1
      x=x0+(i-1)*dx
      v(ifld)=p(i)
      go to(1,2),iopt
c1     write(12)ident,x,y,v,p(i),t(i),v,v,v
1     v(ielv)=t(i)
      write(12)ident,x,y,v
      go to 210
c2     write(12)ident,x,y,v,p(i),hcnst,v,v,v
2     v(ielv)=hcnst
      write(12)ident,x,y,v
210   continue
200   continue
      print*,'Number of data points is',icount
      return
      end
