c esm_fwd.for
c
c  Calculates the magnetic field on a grid of equivalent dipole sources.
c  (The forward version of esmag.for)
c  Jeff Phillips, 27 Oct 92.
c
c  Restricted to 10000 sources, 3000 output columns
      parameter (MAXP = 10000)
      parameter (MAXF = 3000)
      character*50 ifileo,ifilec,ofile
      character*1 ans
      character*4 id1,id2
      character*24 fmt
      character*56 id
      character*8 pgm
      dimension g(MAXF)
      dimension c(MAXP),es(MAXF)
      dimension x(MAXP),y(MAXP),z(MAXP),e(MAXF)
      dimension v(6)
      real mx,my,mz
      common /normal/mx,my,mz
      data dval/0.1e+39/,ddval/1.e+30/,e/MAXF*0.0/
      data pi/3.14159265/
c get the equivalent source file
      print*,'Enter input sources file name:'
      print 1
1     format(1x,'*'$)
      read(5,800)ifileo
800   format(a)
      open(10,file=ifileo,status='old',form='unformatted')
      zmin=1.e38
      zmax=-zmin
      do 400 k=1,MAXP
      read(10,end=404)id1,id2,xs,ys,v
      x(k)=xs
      y(k)=ys
      z(k)=v(1)
      if(z(k).lt.zmin) zmin=z(k)
      if(z(k).gt.zmax) zmax=z(k)
      c(k)=v(2)
400   continue
      print*,'WARNING - Using only the first ',MAXP,' sources.'
404   kmax=k-1
      close(10)
      print*,'Found',kmax,' sources within the following vertical bounds
     1:'
        print*,'        upper = ',zmin
        print*,'        lower = ',zmax
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
        ans='y'
41      print*,'Enter z (up is negative) in units of dx.'
        print 1
        read*,zlevel
        if(zlevel.ge.zmin) then
          print*,'ERROR - Choose a level above the sources.'
          go to 41
        endif
      else if(ans.eq.'n'.or.ans.eq.'N') then
        ans='n'
        print*,'Enter continuation surface grid:'
        print 1
        read(5,800) ifilec
        open(12,file=ifilec,status='old',form='unformatted')
        read(12) id,pgm,nc1,nr1,nz,xo,dx,yo,dy
        nc=nc1
        nr=nr1
        if(nc.gt.MAXF) stop 'ERROR - Too many columns.'
        gmin=1.e38
        gmax=-gmin
        do 20 j=1,nr
        read(12) dlt,(g(i),i=1,nc)
        do 20 i=1,nc
        if(g(i).gt.ddval) go to 20
        if(g(i).lt.gmin) gmin=g(i)
        if(g(i).gt.gmax) gmax=g(i)
   20   continue
        print*,'CONTINUATION SURFACE'
        print*,' title = ',id
        print*,' grid dx/dy = ',dx,dy
        print*,' surface min/max = ',gmin,gmax

        print*,'Enter factor by which surface units are to be'
        print*,'multiplied so that they increase vertically down'
        print*,'and have the same units as dx: '
        print 1
        read(5,*) efac
        if(efac.gt.0.0) then
          gmin=gmin*efac
          gmax=gmax*efac
        else
          temp=gmin
          gmin=gmax*efac
          gmax=temp*efac
        endif
        print*,'The bounds of the continuation surface are:'
        print*,'        upper = ',gmin
        print*,'        lower = ',gmax
      else
        print*,'Please answer "y" or "n".'
        go to 2
      endif
      print*,'Enter mode (0=mag,1=pseudograv):'
      print 1
      read*, mode
      if(mode.eq.0) then
        print*,' Input inclination = ',v(5)
        print*,' Input declination = ',v(6)
        print*,' Enter new inclination, declination:'
        print 1
        read*, ainc,decl
        ainc=ainc*pi/180.
        decl=decl*pi/180.
        mx=cos(ainc)*cos(decl)
        my=cos(ainc)*sin(decl)
        mz=sin(ainc)
      endif
c get filter parameters
      print*,'How should the field be filtered during downward continuat
     1ion?'
      print*,'     0 - no filtering'
      print*,'     1 - slight lowpass'
      print*,'     2 - moderate lowpass'
      print*,'     3 - strong lowpass'
      print*,'     4 - severe lowpass'
   22 print*,'Enter filtering factor: '
      print 1
      read(5,*) ifac
      if(ifac.lt.0.or.ifac.gt.4) then
          print*,'Try again'
          go to 22
      endif
      if(ifac.eq.0.and.ans.eq.'n') go to 24
c get the observation surface grid
        print*,'Enter observation surface grid:'
        print 1
        read(5,800) ifileo
        open(11,file=ifileo,status='old',form='unformatted')
        read(11) id,pgm,nc,nr,nz,xo,dx,yo,dy
        if(nc.gt.MAXF) stop 'ERROR - Too many columns.'
        if(ifac.eq.0) go to 24
        if(ans.eq.'y') go to 40
        if(nc.ne.nc1.or.nr.ne.nr1) stop 'ERROR - Grids not registered.'
40      gmin=1.e38
        gmax=-gmin
        do 10 j=1,nr
        read(11) dlt,(g(i),i=1,nc)
        do 10 i=1,nc
        if(g(i).gt.ddval) go to 10
        if(g(i).lt.gmin) gmin=g(i)
        if(g(i).gt.gmax) gmax=g(i)
   10   continue
        print*,'OBSERVATION SURFACE'
        print*,' title = ',id
        print*,' grid dx/dy = ',dx,dy
        print*,' surface min/max = ',gmin,gmax
        print*,'Enter factor by which surface units are to be'
        print*,'multiplied so that they increase vertically down'
        print*,'and have the same units as dx: '
        print 1
        read(5,*) ofac
        if(ofac.gt.0.0) then
          gmin=gmin*ofac
          gmax=gmax*ofac
        else
          temp=gmin
          gmin=gmax*ofac
          gmax=temp*ofac
        endif
        print*,'The bounds of the observation surface are:'
        print*,'        upper = ',gmin
        print*,'        lower = ',gmax

c adjust the continuation surface to effect the filtering
      icount=0
      if(ans.eq.'y') then
c continuation surface is flat, but observation surface is not
          rewind(11)
          read(11)id
        open(12,file='emag_fwd.tmp',form='unformatted',status='unknown')
          write(12) id,pgm,nc,nr,nz,xo,dx,yo,dy
          do 31 j=1,nr
          read(11) dlt,(g(i),i=1,nc)
          do 30 i=1,nc
          if(g(i)*ofac.lt.zlevel) then
            g(i)=zlevel/ofac-ifac*(zlevel/ofac-g(i))/4.
            ans='n'
            icount=icount+1
          else
            g(i)=zlevel/ofac
          endif
   30     continue
          write(12) dlt,(g(i),i=1,nc)
   31     continue
          close(11)
          efac=ofac
c          ans='n'
          print*,'Continuation surface adjusted at ',icount,' points'
          print*,'     and saved to file emag_fwd.tmp.'
      else
c neither surface is flat
          close(12)
          open(10,file=ifilec,status='old',form='unformatted')
          read(10)id
          rewind(11)
          read(11)id
        open(12,file='emag_fwd.tmp',form='unformatted',status='unknown')
          write(12) id,pgm,nc,nr,nz,xo,dx,yo,dy
          do 35 j=1,nr
          read(11) dlt,(g(i),i=1,nc)
          read(10) dlt,(e(i),i=1,nc)
          do 34 i=1,nc
          if(e(i)*efac.gt.g(i)*ofac) then
            e(i)=e(i)-ifac*(e(i)-g(i)*ofac/efac)/4.
            icount=icount+1
          endif
   34     continue
          write(12) dlt,(e(i),i=1,nc)
   35     continue
          close(10)
          close(11)
          print*,'Continuation surface adjusted at ',icount,' points'
          print*,'     and saved to file emag_fwd.tmp.'
      endif
24    if(ans.eq.'y') go to 25
      rewind(12)
      read(12) id
25    print*,'Enter output grid file name:'
      print 1
      read(5,800)ofile
      open(13,file=ofile,status='unknown',form='unformatted')
      print*,'Enter id:'
      print 1
      read(5,800) id
      pgm='emag_fwd'
c  Calculate grid.
      print*,'Input datum = ',v(4)
106   print*,'Enter datum to be added to calculated field [input datum]:
     1'
      datum=v(4)
      print 1
      read (*,'(a)') ifileo
      leng=len_trim(ifileo)
      if(leng.ne.0) then
        lenm=leng-index(ifileo(1:leng),'.')
        if(lenm.eq.leng) lenm=0
        write(fmt,'(a,i3,a,i3,a)') '(f',leng,'.',lenm,')'
        read(ifileo(1:leng),fmt,err=106) datum
      endif
      print*,datum
c
      nz=1
      write(13)id,pgm,nc,nr,nz,xo,dx,yo,dy
      do 75 j=1,nr
      write(6,100) j,nr
      yg=yo+(j-1)*dy
      if(ans.eq.'n'.or.ans.eq.'N') then
        read(12) dlt,(es(i),i=1,nc)
        do 72 i=1,nc
        if(es(i).lt.ddval) then
          e(i)=es(i)*efac
        endif
72      continue
      endif
      do 74 i=1,nc
      xg=xo+(i-1)*dx
      if(es(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
        if(mode.eq.0) then
          rsg=dipole(y(k),x(k),z(k),c(k),yg,xg,zlevel)
        else
          rsg=pole(y(k),x(k),z(k),c(k),yg,xg,zlevel)
        endif
      else
        if(mode.eq.0) then
          rsg=dipole(y(k),x(k),z(k),c(k),yg,xg,e(i))
        else
          rsg=pole(y(k),x(k),z(k),c(k),yg,xg,e(i))
        endif
      endif
73    g(i)=g(i)+rsg
74    continue
      write(13)  dlt,(g(i),i=1,nc)
75    continue
100   format(1h+,'working on row ',i5,' of ',i5)
999   continue
      close(11)
      stop
      end
      function dipole(xd,yd,zd,cnst,x,y,z)
      real mx,my,mz
      common /normal/mx,my,mz
      rx=(x-xd)
      ry=(y-yd)
      rz=(z-zd)
      r2=rx**2+ry**2+rz**2
      r=sqrt(r2)
      dipole=cnst*(3.*(mx*rx+my*ry+mz*rz)**2/r2-1.)/r**3
      end
      function pole(xd,yd,zd,cnst,x,y,z)
      rx=(x-xd)
      ry=(y-yd)
      rz=(z-zd)
      r2=rx**2+ry**2+rz**2
      pole=cnst/r2
      end
