c#######################################################################
      module read_coord_interface
      interface
        subroutine read_coord (fname,n,x)
        character(*) :: fname
        integer :: n
        real(8), dimension(:), pointer :: x
        end
      end interface
      end module
c#######################################################################
      program DIPOLE_PLUS_BIPOLE
c
c-----------------------------------------------------------------------
c
c ****** Get the normal magnetic field Br at r=R0 for a model with
c ****** a large-scale dipole and a sub-surface bipole that models
c ****** an active region.
c
c-----------------------------------------------------------------------
c
      use read_coord_interface
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
      real(8) :: mr,mt,mp
      real(8) :: rd,td,pd
      real(8) :: r0=1._8
c
      integer :: i,j,ierr
      real(8) :: xd,yd,zd,mx,my,mz,b0
      real(8) :: st,ct,sp,cp
      character(512) :: tfile,pfile,brfile
c
      integer :: nt,np
      real(8), dimension(:), pointer :: t,p
      real(8), dimension(:,:), pointer :: br
c
c-----------------------------------------------------------------------
c
c ****** Read the theta and phi meshes.
c
      write (*,'(a,$)') 'Enter the t mesh file name: '
      read (*,'(a)') tfile
      write (*,'(a,$)') 'Enter the p mesh file name: '
      read (*,'(a)') pfile
c
      call read_coord (tfile,nt,t)
      call read_coord (pfile,np,p)
c
c ****** Allocate the Br array.
c
      allocate (br(nt,np))
c
c ****** Get the bipole parameters.
c
      write (*,'(a,$)') 'Enter the bipole position [r,t,p]: '
      read (*,*) rd,td,pd
      write (*,'(a,$)') 'Enter the bipole orientation [mr,mt,mp]: '
      read (*,*) mr,mt,mp
c
c ****** Convert the bipole position to Cartesian coordinates.
c
      call s2c (rd,td,pd,xd,yd,zd)
c
c ****** Convert the bipole direction to Cartesian coordinates.
c
      st=sin(td)
      ct=cos(td)
      sp=sin(pd)
      cp=cos(pd)
c
      mx=mr*st*cp+mt*ct*cp-mp*sp
      my=mr*st*sp+mt*ct*sp+mp*cp
      mz=mr*ct   -mt*st
c
c ****** Evaluate the magnetic field at r=R0 due to the bipole.
c
      do j=1,np
        do i=1,nt
          br(i,j)=0.
          call add_bipole (xd,yd,zd,mx,my,mz,
     &                     r0,t(i),p(j),
     &                     br(i,j))
        enddo
      enddo
c
c ****** Read in the strength of the dipole.
c
      write (*,'(a,$)') 'Enter the dipole polar Br field: '
      read (*,*) b0
c
c ****** Add in the global magnetic field.
c
      do j=1,np
        do i=1,nt
          br(i,j)=br(i,j)+b0*cos(t(i))
        enddo
      enddo
c
c ****** Write the magnetic field to an HDF file.
c
      write (*,'(a,$)') 'Enter the Br HDF file name: '
      read (*,'(a)') brfile
c
      call wrhdf_2d (brfile,.true.,nt,np,br,t,p,.true.,ierr)
c
      end
c#######################################################################
      subroutine add_bipole (xd,yd,zd,mx,my,mz,r,t,p,br)
c
c-----------------------------------------------------------------------
c
c ****** Add the magnetic field at position (R,T,P) in spherical
c ****** coordinates for a bipole centered at (XD,YD,ZD) with magnetic
c ****** moment (MX,MY,MZ).
c
c ****** The radial component of the magnetic field is added to BR.
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
      real(8), parameter :: one=1._8
      real(8), parameter :: three=3._8
c
c-----------------------------------------------------------------------
c
      real(8) :: xd,yd,zd
      real(8) :: mx,my,mz
      real(8) :: r,t,p
      real(8) :: br
c
c-----------------------------------------------------------------------
c
      real(8) :: xp,yp,zp
      real(8) :: ax,ay,az
      real(8) :: bx,by,bz
      real(8) :: bt,bp
      real(8) :: r2i,r3i,mdotx
      real(8) :: x,y,z
      real(8) :: st,ct,sp,cp
c
c-----------------------------------------------------------------------
c
c ****** Get the evaluation point in Cartesian coordinates.
c
      call s2c (r,t,p,x,y,z)
c
      xp=x-xd
      yp=y-yd
      zp=z-zd
c
      r2i=one/(xp**2+yp**2+zp**2)
      r3i=one/sqrt(xp**2+yp**2+zp**2)**3
c
c ****** The components of the vector potential are given by:
c
      ax=(my*zp-mz*yp)*r3i
      ay=(mz*xp-mx*zp)*r3i
      az=(mx*yp-my*xp)*r3i
c
c ****** The components of the magnetic field are given by:
c
      mdotx=(mx*xp+my*yp+mz*zp)
c
      bx=(-mx+three*mdotx*xp*r2i)*r3i
      by=(-my+three*mdotx*yp*r2i)*r3i
      bz=(-mz+three*mdotx*zp*r2i)*r3i
c
      st=sin(t)
      ct=cos(t)
      sp=sin(p)
      cp=cos(p)
c
      br= bx*st*cp+by*st*sp+bz*ct
      bt= bx*ct*cp+by*ct*sp-bz*st
      bp=-bx*sp   +by*cp
c
      return
      end
c#######################################################################
      subroutine s2c (r,t,p,x,y,z)
c
c-----------------------------------------------------------------------
c
c ****** Convert from spherical coordinates (R,T,P) 
c ****** to Cartesian coordinates (X,Y,Z).
c
c ****** This routine assumes that T and P are in radians.
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
      real(8) :: r,t,p,x,y,z
c
c-----------------------------------------------------------------------
c
      real(8) :: st
c
c-----------------------------------------------------------------------
c
      st=sin(t)
      x=r*st*cos(p)
      y=r*st*sin(p)
      z=r*cos(t)
c
      return
      end
c#######################################################################
      subroutine read_coord (fname,n,x)
c
c-----------------------------------------------------------------------
c
c ****** Read a coordinate array.
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
      character(*) :: fname
      integer :: n
      real(8), dimension(:), pointer :: x
c
c-----------------------------------------------------------------------
c
      integer :: ierr,i
      real(8) :: dum
c
c-----------------------------------------------------------------------
c
      call ffopen (1,fname,'r',ierr)
      if (ierr.ne.0) go to 910
c
c ****** Skip the header.
c
      read (1,*,err=900,end=900)
c
      n=0
      do
        read (1,*,err=900,end=100)
        n=n+1
      enddo
  100 continue
c
      close (1)
c
c ****** Allocate memory for the coordinate.
c
      allocate (x(n))
c
c ****** Read the coordinates.
c
      call ffopen (1,fname,'r',ierr)
      if (ierr.ne.0) go to 910
c
c ****** Skip the header.
c
      read (1,*,err=900,end=900)
c
      do i=1,n
        read (1,*,err=900,end=900) dum,x(i),dum,dum
      enddo
c
      close (1)
c
      write (*,*)
      write (*,*) '### File name: ',trim(fname)
      write (*,*) '### Number of points read = ',n
c
      return
c
c ****** Error returns.
c
  900 continue
c
      write (*,*)
      write (*,*) '### ERROR in READ_COORD:'
      write (*,*) '### The format of the coordinate file is not valid.'
      write (*,*) 'File name: ',trim(fname)
      call exit (1)
c
  910 continue
      write (*,*)
      write (*,*) '### ERROR in READ_COORD:'
      write (*,*) '### Could not open the coordinate file.'
      write (*,*) 'File name: ',trim(fname)
      call exit (1)
c
      end