c dipole.for
c
c  magnetic field of a dipole
c  profile or grid case, allowing for induced plus remanent magnetization
c  Note that a specified profile azimuth is associated with a profile, 
c  whereas the grid is assumed columns are north-south, rows are east-west.
c  Lin Cordell
c
      character*50 ofile
      dimension f(2000),id(14),pgm(2)
      data pi/3.1415927/
c      dval2=1.e20
      rad=pi/180.
      pi=3.1415927
      print*,'Enter output file name:'
      print 1
1     format(' *'$)
      read(5,100)ofile
100   format(a50)
      pgm(1)='dipo'
      pgm(2)='le  '
      print*,'Enter ncol, nrow (for profile, nrow=1):'
      print 1
      read*,ncol, nrow
      if(ncol.gt.2000)go to 900
      nz=1
      if(nrow.eq.1)then
        print*,'Enter x0, dx:'
        print 1
        read*,x0,dx
        y0=x0
        dy=dx
        print*,'Enter profile azimuth, + clockwise from true north:'
        print 1
        read*,paz
        print*,'Enter dipole coordinates x,z (z=depth > 0):'
        print 1
        read*,xi,zeta
        xi=xi-x0
        print*,xi
      else
        print*,'Enter x0,dx,y0,dy:'
        print 1
        read*,x0,dy,y0,dx
        print*,'Enter dipole coordinates x,y,z :' 
        print*,'x is + east; y is + north, and z=depth; z>0):'
        print 1
        read*, ce,cn,zeta
c  Note internal frame is x north, y east, z down.
        xi=cn-y0
        eta=ce-x0
      endif
      print*,'Enter radius of sphere (km):'
      print 1
      read*,vol
      vol=4.18879*(vol**3)
      print*,'Enter inclination, declination of inducing field:'
      print 1
      read*,rinc,dec
      rinc=rinc*rad
      dec=dec*rad
      si=sin(rinc)
      ci=cos(rinc)
      if(nrow.eq.1)then
        paz=paz*rad
        pc=cos(paz-dec)
        a=sqrt((ci*pc)**2+si*si)
        t1=(ci*pc)/a
        t2=0.
        t3=si/a
      else
        sd=sin(dec)
        cd=cos(dec)
        t1=ci*cd
        t2=ci*sd
        t3=si
      endif
      print*,'Enter inducing-field strength (nT):'
      print 1
      read*,ho
      print*,'Enter susceptibility:'
      print 1
      read*,susp
      rmagi=ho*susp
      print*,'Enter remanent magnetization intensity:'
      print 1
      read*,rem
      rem1=0.
      rem2=0.
      rem3=0.
      if(rem.ne.0.)then
        print*,'Enter remanent inclination, declination:'
        print 1
        read*,rinc,dec
        rinc=rinc*rad
        dec=dec*rad
        si=sin(rinc)
        ci=cos(rinc)
        if(nrow.eq.1)then
          paz=paz*rad
          pc=cos(paz-dec)
          a=sqrt((ci*pc)**2+si*si)
          r1=(ci*pc)/a
          r2=0.
          r3=si/a
        else
          sd=sin(dec)
          cd=cos(dec)
          r1=ci*cd
          r2=ci*sd
          r3=si
        endif
      endif
      c1=vol*susp*ho
      c2=vol*rem
      rj1=c1*t1+c2*r1
      rj2=c1*t2+c2*r2
      rj3=c1*t3+c2*r3
      open(11,file=ofile,form='unformatted')
      write(6,202)
202   format(' Enter id:'/' *'$)
      read(5,203)(id(i),i=1,14)
203   format(14a4)
      write(11)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      zsq=zeta*zeta
      if(nrow.eq.1)go to 50
      do 40 j=1,nrow
      x=xi-(j-1)*dy
      xsq=x*x
      xz=x*zeta
      do 30 i=1,ncol
      y=eta-(i-1)*dy
      ysq=y*y
      xy=x*y
      yz=y*zeta
      d=sqrt(xsq+ysq+zsq)
      d5=1./d**5
      d3=1./d**3
      sxx=3.*xsq*d5-d3
      syy=3.*ysq*d5-d3
      szz=3.*zsq*d5-d3
      sxy=3.*xy*d5
      sxz=3.*xz*d5
      syz=3.*yz*d5
      f(i)=(rj1*sxx+rj2*sxy+rj3*sxz)*t1+(rj1*sxy+rj2*syy+rj3*syz)*t2
     & +(rj1*sxz+rj2*syz+rj3*szz)*t3
30    continue
      call rowio(f,y,1,ncol,11,2)
40    continue
      go to 998
50    continue
      c1=vol*susp*ho
      c2=vol*rem
      do 60 i=1,ncol
      x=xi-(i-1)*dx
      xsq=x*x
      xz=x*zeta
      d=sqrt(xsq+zsq)
      d3=1./d**3
      d5=1./d**5
      sxx=3.*xsq*d5-d3
      szz=3.*zsq*d5-d3
      sxz=3.*xz*d5
      p=t1*sxx+t3*sxz
      q=t1*sxz+t3*szz
      f(i)=c1*(t1*p+t3*q)+c2*r1*p+c2*r3*q
60    continue
      call rowio(f,y,1,ncol,11,2)
      go to 998
900   write(6,901)
901   format(' Can''t handle it.')
998   close(11)
      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
