c es_fwd.for
c
c  To calculate field of a point_sources file on a grid.
c  (The forward version of es.for)
c  Input from post (.pos) data file.
c  Lin Cordell, 1 Mar 91.
c  revised Jeff Phillips, 27 Oct 92.
      character*50 ifile,ofile
      character*1 ans
      character*4 id1,id2
      dimension id(14),pgm(2)
      dimension g(3000),e(3000)
      dimension x(8000),y(8000),z(8000),c(8000)
      data dval/0.1e+39/,ddval/1.e+30/,e/3000*0.0/
      print*,'Enter input sources file name:'
      print 1
1     format(1x,'*'$)
      read(5,800)ifile
800   format(a)
      open(10,file=ifile,status='old',form='unformatted')
2     print*,'Do you want to compute the field on a constant level?'
      print 1
      read(5,800) ans
      if(ans.eq.'y'.or.ans.eq.'Y') then
        print*,'Enter z (up is negative) in km.'
        print 1
        read*,zlevel
      else if(ans.eq.'n'.or.ans.eq.'N') then
        print*,'Enter input surface grid:'
        print 1
        read(5,800) ifile
        open(12,file=ifile,status='old',form='unformatted')
        print*,'Enter factor by which surface units are to be'
        print*,'multiplied in order to make them negative up in km:'
        read(5,*) efac
      else
        print*,'Please answer "y" or "n"'
        go to 2
      endif
      do 400 k=1,8000
      read(10,end=404)id1,id2,x(k),y(k),z(k),c(k),v,v,v,v
      if(ans.eq.'y'.or.ans.eq.'Y') then
        z(k)=(z(k)-zlevel)**2
      endif
400   continue
      print*,'WARNING - using only the first 8000 sources'
404   kmax=k-1
      close(10)
      print*,kmax,' sources.'
      print*,'Enter output grid file name:'
      print 1
      read(5,800)ofile
      open(11,file=ofile,status='new',form='unformatted')
      if(ans.eq.'y'.or.ans.eq.'Y') then
        print*,'Enter x0,y0,dx,dy,ncol,nrow'
        print 1
        read*,x0,y0,dx,dy,ncol,nrow
      else
        read(12) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      endif
      if(ncol.gt.3000) stop 'too many columns'
      print*,'Enter id:'
      print 1
      read(5,4)(id(i),i=1,14)
4     format(14a4)
      pgm(1)='es_f'
      pgm(2)='wd  '
c  Calculate grid.
      print*,'Enter datum to be added to calculated field:'
      print 1
      read*,datum
      nz=1
      write(11)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      do 75 j=1,nrow
      yg=y0+(j-1)*dy
      if(ans.eq.'n'.or.ans.eq.'N') read(12)xg,(e(i),i=1,ncol)
      do 74 i=1,ncol
      xg=x0+(i-1)*dx
      if(e(i).gt.ddval) then
        g(i)=dval
        go to 74
      else
        g(i)=datum
      endif
      do 73 k=1,kmax
      if(ans.eq.'y'.or.ans.eq.'Y') then
        rsg=(xg-x(k))**2+(yg-y(k))**2+z(k)
      else
        rsg=(xg-x(k))**2+(yg-y(k))**2+(e(i)*efac-z(k))**2
      endif
      rsg=c(k)/sqrt(rsg)
73    g(i)=g(i)+rsg
74    continue
      write(11)yg,(g(i),i=1,ncol)
75    continue
999   continue
      close(11)
      stop
      end
   
