magcalc.f90 Source File

Program MagCalc
!! Need program statement for FORD
!! THIS IS THE MAIN PROGRAM FOR COMPUTING MAGNETIC FIELDS
!! FROM OUTPUT FROM A SIMULATIONS DONE BY GEMINI3D.
!! THIS PROGRAM VERY MUCH MIRRORS THE SETUP OF THE MAIN GEMINI.F90 CODE.

use, intrinsic :: iso_fortran_env, only : stderr=>error_unit
use, intrinsic :: ieee_arithmetic, only : ieee_is_nan

use phys_consts, only : pi,mu0, wp, re, debug
use grid, only : lx1, lx2, lx3, read_grid, clear_grid, lx2all,lx3all,grid_size
use mesh, only : curvmesh
use timeutils, only : dateinc
use io, only : read_configfile,input_plasma_currents,create_outdir_mag,output_magfields
use mpimod, only: mpi_sum, mpi_comm_world, &
mpisetup, mpibreakdown, mpigrid, mpi_manualgrid, halo_end, &
lid, lid2, lid3, myid, myid2, myid3, mpi_realprec, &
tagdv, tagjx, tagjy, tagjz, tagrcubed, tagrx, tagry, tagrz

implicit none

!> VARIABLES READ IN FROM CONFIG.DAT FILE
integer, dimension(3) :: ymd
!! year,month,day of simulation
real(wp) :: UTsec
!! UT (s)
real(wp) :: UTsec0
!! UT start time of simulation (s)
real(wp) :: tdur
!! duration of simulation
real(wp), dimension(3) :: activ
!! f10.7a,f10.7,ap
real(wp) :: tcfl
!! target CFL number
real(wp) :: Teinf
!! exospheric temperature
integer :: potsolve
!! what type of potential solve
integer :: flagperiodic
!! toggles whether or not the grid is treated as periodic in the x3 dimension (affects some of the message passing)
integer :: flagoutput
!! what type of output to do (1 - everything; 2 - avg'd parms.; 3 - ne only)
integer :: flagcap
!! internal capacitance?

!> INPUT AND OUTPUT FILES
character(:), allocatable :: infile    !command line argument input file
character(:), allocatable :: outdir    !" " output directory
character(:), allocatable :: indatsize,indatgrid    !grid size and data filenames
character(:), allocatable :: fieldpointfile

!! GRID STRUCTURE
type(curvmesh) :: x    !structure containg grid locations, finite differences, etc.:  see grid module for details

!STATE VARIABLES
real(wp), dimension(:,:,:), allocatable :: J1,J2,J3      !electrodynamic state variables

!TEMPORAL VARIABLES
real(wp) :: t=0._wp,dt      !time from beginning of simulation (s) and time step (s)
real(wp) :: tout,dtout    !time for next output and time between outputs
real(wp) :: tstart,tfin   !temp. vars. for measuring performance of code blocks
integer :: it,isp        !time and species loop indices

!WORK ARRAYS
integer :: flag2D
real(wp), dimension(:,:,:), allocatable :: xp,yp,zp     !source coordinates computed from simulation sub-grid
real(wp), dimension(:), allocatable  :: xf,yf,zf        !field point coordinates (flat list)
real(wp), dimension(:), allocatable :: r,theta,phi
real(wp), dimension(:,:,:), allocatable :: dV
integer :: ipoints,lpoints    !size of field point arrays
real(wp), dimension(:,:,:), allocatable :: proj_e1er,proj_e2er,proj_e3er
real(wp), dimension(:,:,:), allocatable :: proj_e1etheta,proj_e2etheta,proj_e3etheta
real(wp), dimension(:,:,:), allocatable :: proj_e1ephi,proj_e2ephi,proj_e3ephi
real(wp), dimension(:,:,:), allocatable :: Jx,Jy,Jz
real(wp), dimension(:,:,:), allocatable :: Rx,Ry,Rz,Rcubed
real(wp), dimension(:,:,:), allocatable :: integrand,integrandavg
real(wp), dimension(:,:,:), allocatable :: alt
real(wp), dimension(:), allocatable :: Br,Btheta,Bphi
real(wp), dimension(:), allocatable :: Brall,Bthetaall,Bphiall
real(wp), dimension(:,:), allocatable :: Jxend,Jyend,Jzend,Rxend,Ryend,Rzend,Rcubedend,dVend
real(wp), dimension(:,:), allocatable :: integrandend
real(wp), dimension(:,:), allocatable :: integrandavgend

real(wp), dimension(:,:), allocatable :: Jxtop,Jytop,Jztop,Rxtop,Rytop,Rztop,Rcubedtop,dVtop
real(wp), dimension(:,:), allocatable :: integrandtop
real(wp), dimension(:,:), allocatable :: integrandavgtop

integer :: ix1,ix2,ix3
real(wp) :: rmean,thetamean

!NEUTRAL PERTURBATION VARIABLES (UNUSED)
integer :: flagdneu                  !toggles neutral perturbations (0 - none; 1 - file-based neutral inputs)
integer :: interptype                !toggles whether the neutral input data are interpreted (0 - Cartesian; 1 - axisymmetric)
real(wp) :: dxn,drhon,dzn                 !finite differences for the neutral input data in the horizontal and vertical directions
real(wp) :: sourcemlat,sourcemlon     !mag. lat./long for the neutral source location
character(:), allocatable :: sourcedir          !directory where neutral input data are located
real(wp) :: dtneu                     !time interval [s] in between neutral inputs

!PRECIPITATION FILE INPUT VARIABLES (UNUSED)
integer :: flagprecfile              ! flag toggling precipitation file input (0 - no; 1 - yes)
real(wp) :: dtprec                    ! time interval between precip. inputs
character(:), allocatable :: precdir ! directory containing precip. input files

!! ELECTRIC FIELD FILE INPUT VARIABLES (UNUSED)
integer :: flagE0file                !< flag toggling electric field (potential BCs) file input (0 - no; 1 - yes)
real(wp) :: dtE0                      !< time interval between electric field file inputs
character(:), allocatable :: E0dir   !< directory containing electric field file input data

!! GLOW MODULE INPUT VARIABLES
integer :: flagglow                     !flag toggling GLOW module run (include aurora) (0 - no; 1 - yes)
real(wp) :: dtglow                      !time interval between GLOW runs (s)
real(wp) :: dtglowout                   !time interval between GLOW auroral outputs (s)

!! FOR SPECIFYING THE PROCESS GRID
integer :: lid2in,lid3in

!! FOR HANDLING INPUT
integer :: argc, ierr
character(256) :: argv
logical :: file_exists

!! REGULATOR FOR 1/R^3
!real(wp), parameter :: R3min=1d11     !works well for 192x192 in situ
real(wp), parameter :: R3min=1d9


real(wp), parameter :: Rmin=5d3


!! ## MAIN PROGRAM
debug=.true.     !FIXME:  hardcode this in for now, needs to be set based on user input...

argc = command_argument_count()
if (argc < 2) error stop 'magcalc.f90 --> must specify .ini file to configure simulation and location of'// &
        'output data from simulation (i.e. plasma parameters and fields) and input field point file'


!! INITIALIZE MESSING PASSING VARIABLES, IDS ETC.
call mpisetup()
print *, 'magcalc.f90 --> Process:  ',myid,' of:  ',lid-1,' online...'


!READ CONFIG FILE FROM OUTPUT DIRECTORY
call get_command_argument(1,argv)
outdir=trim(argv)
infile=outdir//'/inputs/config.ini'
!infile = trim(argv)
inquire(file=infile,exist=file_exists)    !needed to deal with ini vs. nml inputs...
if ( .not. file_exists) then
  infile=outdir//'/inputs/config.nml'
end if
if (myid==0) then
  print *, 'Simulation data directory:  ',outdir
  print *, 'Input config file:  ',infile
end if
call read_configfile(infile,ymd,UTsec0,tdur,dtout,activ,tcfl,Teinf,potsolve,flagperiodic,flagoutput,flagcap, &
                     indatsize,indatgrid,flagdneu,interptype,sourcemlat,sourcemlon,dtneu,dxn,drhon,dzn,sourcedir,flagprecfile, &
                     dtprec,precdir,flagE0file,dtE0,E0dir,flagglow,dtglow,dtglowout)


!ESTABLISH A PROCESS GRID
!call grid_size(indatsize)
!call mpigrid(lx2all,lx3all)    !following grid_size these are in scope
!!CHECK THE GRID SIZE AND ESTABLISH A PROCESS GRID
call grid_size(indatsize)
if (argc > 2) then   !user specified process grid
  call get_command_argument(3,argv)
  read(argv,*) lid2in
  call get_command_argument(4,argv)
  read(argv,*) lid3in
  call mpi_manualgrid(lx2all,lx3all,lid2in,lid3in)
else     !try to decide the process grid ourself
  call mpigrid(lx2all,lx3all)    !following grid_size these are in scope
end if


!> LOAD UP THE GRID STRUCTURE/MODULE VARS. FOR THIS SIMULATION - THIS ALSO PERMUTES DIMENSIONS OF 2D GRID, IF NEEDED
if (myid==0) then
  print*, 'Process grid established; reading in full grid file now...'
end if
call read_grid(indatsize,indatgrid,flagperiodic,x)
!! read in a previously generated grid from filename listed in input file, distribute subgrids to individual workers
if (lx2==1) then
  flag2D=1
  !! a 2D simulation was done, which changes how the integrations go...
else
  flag2D=0
end if


! FIXME:  need to copy the input grid file into the output directory


!GRAB THE INFO FOR WHERE THE OUTPUT CALCULATIONS ARE STORED
!call get_command_argument(2,argv)
!outdir = trim(argv)    !this should be the base output directory that the simulation results have been stored in
call get_command_argument(2,argv)
fieldpointfile=trim(argv)
!! this file contains the field points at which we are computing magnetic perturbations, it will be copied into the output directory


!SET UP DIRECTORY TO STORE OUTPUT FILES
if (myid==0) then
  print*, 'Creating output directory...'
  call create_outdir_mag(outdir,fieldpointfile)
end if


!ALLOCATE ARRAYS (AT THIS POINT ALL SIZES ARE SET FOR EACH PROCESS SUBGRID)
allocate(J1(lx1,lx2,lx3),J2(lx1,lx2,lx3),J3(lx1,lx2,lx3))
allocate(Jx(lx1,lx2,lx3),Jy(lx1,lx2,lx3),Jz(lx1,lx2,lx3))


!NOW DEAL WITH THE UNMPRIMED COORDINATES
block
  integer :: u
  open(newunit=u,file=fieldpointfile,status='old',form='unformatted',access='stream',action='read')
  read(u) lpoints    !size of coordinates for field points
  if (myid==0) print *, 'magcalc.f90 --> Number of field points:  ',lpoints
  allocate(r(lpoints),theta(lpoints),phi(lpoints))
  read(u) r,theta,phi
  close(u)
end block

if (myid==0) print *, 'magcalc.f90 --> Range of r,theta,phi',minval(r),maxval(r),minval(theta), &
                           maxval(theta),minval(phi),maxval(phi)
rmean=sum(r)/size(r)
thetamean=sum(theta)/size(theta)
allocate(xf(lpoints),yf(lpoints),zf(lpoints))
xf(:)=r(:)
!yf(:)=r(:)*theta(:)
!zf(:)=r(:)*sin(theta(:))*phi(:)
yf(:)=rmean*theta(:)
zf(:)=rmean*sin(thetamean)*phi(:)


!GET POSITIONS (CARTESIAN) SET UP FOR MAGNETIC COMPUTATIONS.  THESE ARE PRIMED COORDINATES (SOURCE COORDS, I.E. THE SIM GRID)
allocate(xp(lx1,lx2,lx3),yp(lx1,lx2,lx3),zp(lx1,lx2,lx3))
xp(:,:,:)=x%alt(:,:,:)+Re                               !radial distance from Earth's center
!yp(:,:,:)=xp(:,:,:)*x%theta(:,:,:)                      !southward distance (in the direction of the theta spherical coordinate)
!zp(:,:,:)=xp(:,:,:)*sin(x%theta(:,:,:))*x%phi(:,:,:)    !eastward distance
yp(:,:,:)=rmean*x%theta(:,:,:)
!! the integrations are being treated as Cartesian so flatten out the local spherical coordinates into cartesian, as well
zp(:,:,:)=rmean*sin(thetamean)*x%phi(:,:,:)


!print*, myid2,myid3,'--> field point min/max data:  ',minval(xp),maxval(xp),minval(yp),maxval(yp),minval(zp),maxval(zp)


!COMPUTE A SOURCE DIFFERENTIAL VOLUME FOR INTEGRALS
allocate(dV(lx1,lx2,lx3))
allocate(dVend(lx1,lx2),Jxend(lx1,lx2),Jyend(lx1,lx2),Jzend(lx1,lx2))
allocate(Rxend(lx1,lx2),Ryend(lx1,lx2),Rzend(lx1,lx2),Rcubedend(lx1,lx2))
allocate(integrandend(lx1,lx2),integrandavgend(lx1-1,max(lx2-1,1)))

!should these be size lx3???
allocate(dVtop(lx1,lx3),Jxtop(lx1,lx3),Jytop(lx1,lx3),Jztop(lx1,lx3))
allocate(Rxtop(lx1,lx3),Rytop(lx1,lx3),Rztop(lx1,lx3),Rcubedtop(lx1,lx3))
allocate(integrandtop(lx1,lx3),integrandavgtop(lx1-1,max(lx3-1,1)))

if (flag2D/=1) then   !3D differential volume
  do ix3=1,lx3
    do ix2=1,lx2
      do ix1=1,lx1
        dV(ix1,ix2,ix3)=x%h1(ix1,ix2,ix3)*x%h2(ix1,ix2,ix3)*x%h3(ix1,ix2,ix3)*x%dx1(ix1)*x%dx2(ix2)*x%dx3(ix3)
      end do
    end do
  end do
else         !plane geometry assumption
  do ix3=1,lx3
    do ix2=1,lx2
      do ix1=1,lx1
        dV(ix1,ix2,ix3)=x%h1(ix1,ix2,ix3)*x%h3(ix1,ix2,ix3)*x%dx1(ix1)*x%dx3(ix3)
      end do
    end do
  end do
end if


!> GET END VOLUMES SO THE INTEGRALS ARE 'COMPLETE'
call halo_end(dV,dVend,dVtop,tagdV)
!! need to define the differential volume on the edge of this x3-slab in


!print*, myid2,myid3,'--> dV vals.',minval(dV),maxval(dV),minval(dVend),maxval(dVend),minval(dVtop),maxval(dVtop)


!> COMPUTE NEEDED PROJECTIONS
allocate(proj_e1er(lx1,lx2,lx3),proj_e2er(lx1,lx2,lx3),proj_e3er(lx1,lx2,lx3))
allocate(proj_e1etheta(lx1,lx2,lx3),proj_e2etheta(lx1,lx2,lx3),proj_e3etheta(lx1,lx2,lx3))
allocate(proj_e1ephi(lx1,lx2,lx3),proj_e2ephi(lx1,lx2,lx3),proj_e3ephi(lx1,lx2,lx3))
allocate(alt(lx1,lx2,lx3))
alt(:,:,:)=x%alt
proj_e1er(:,:,:)=sum(x%e1*x%er,4)
!! fourth dimension of unit vectors is the 3 Cartesian components of each vector
proj_e2er(:,:,:)=sum(x%e2*x%er,4)
proj_e3er(:,:,:)=sum(x%e3*x%er,4)
proj_e1etheta(:,:,:)=sum(x%e1*x%etheta,4)
proj_e2etheta(:,:,:)=sum(x%e2*x%etheta,4)
proj_e3etheta(:,:,:)=sum(x%e3*x%etheta,4)
proj_e1ephi(:,:,:)=sum(x%e1*x%ephi,4)
proj_e2ephi(:,:,:)=sum(x%e2*x%ephi,4)
proj_e3ephi(:,:,:)=sum(x%e3*x%ephi,4)


!> DEALLOCATE GRID MODULE VARIABLES TO SAVE MEMORY (PROGRAM DOESN'T ACTUALLY NEED THESE ONCE X,Y,Z CREATED
call clear_grid(x)
deallocate(r,theta,phi)


!> STORAGE FOR MAGNETIC FIELD CALCULATIONS
allocate(Rx(lx1,lx2,lx3),Ry(lx1,lx2,lx3),Rz(lx1,lx2,lx3))
allocate(Rcubed(lx1,lx2,lx3))
allocate(integrand(lx1,lx2,lx3),integrandavg(lx1-1,max(lx2-1,1),lx3-1))
!! latter is cell centered hence -1 in size, max is needed to prevent zero sized array
allocate(Br(lpoints),Btheta(lpoints),Bphi(lpoints))
allocate(Brall(lpoints),Bthetaall(lpoints),Bphiall(lpoints))
!! only used by root, but I think workers need to have space allocated for this

!! MAIN LOOP
UTsec=UTsec0; it=1; t=0d0; tout=t;
call dateinc(dtout,ymd,UTsec)     !skip first file
it=3                              !don't trigger any special adaptations to filename
do while (t<tdur)
  !TIME STEP CALCULATION
  dt=dtout    !only compute magnetic field at times when we've done output


  !READ IN THE FULL PLASMA AND FIELD DATA FROM THE OUTPUT FILE (NOTE THAT WE NEED TO KNOW OUTPUT TYPE DONE)
  if (it==1) then
    UTsec=UTsec+0.000001d0
  end if
  if (it==2) then
    UTsec=UTsec-0.000001d0
  end if
  call input_plasma_currents(outdir,flagoutput,ymd,UTsec,J1,J2,J3)    !now everyone has their piece of data


  !FORCE PARALLEL CURRENTS TO ZERO BELOW SOME ALTITUDE LIMIT
  if(myid==0) print *, 'Zeroing out low altitude currents (these are basically always artifacts)...'
  where (alt<75d3)
    J1=0d0
    J2=0d0
    J3=0d0
  end where


  !DEAL WITH THE WEIRD EDGE ARTIFACTS THAT WE GET IN THE PARALLEL CURRENT
  !SOMETIMES, THIS IS FOR THE X2 DIRECTION
  if(myid==0) print *, 'Fixing current edge artifacts...'
  if (myid3==lid3-1) then
    if (lx3>2) then    !do a ZOH
      J1(:,:,lx3-1)=J1(:,:,lx3-2)
      J1(:,:,lx3)=J1(:,:,lx3-2)
    else
      J1(:,:,lx3-1)=0d0
      J1(:,:,lx3)=0d0
    end if
  end if
  if (myid3==0) then
    if (lx3>2) then    !do a ZOH
      J1(:,:,1)=J1(:,:,3)
      J1(:,:,2)=J1(:,:,3)
    else
      J1(:,:,1)=0d0
      J1(:,:,2)=0d0
    end if
  end if

  !X3 EDGES
  if (myid2==lid2-1) then
    if (lx2>2) then    !do a ZOH
      J1(:,lx2-1,:)=J1(:,lx2-2,:)
      J1(:,lx2,:)=J1(:,lx2-2,:)
    else
      J1(:,lx2-1,:)=0d0
      J1(:,lx2-1,:)=0d0
    end if
  end if
  if (myid2==0) then
    if (lx2>2) then    !do a ZOH
      J1(:,1,:)=J1(:,3,:)
      J1(:,2,:)=J1(:,3,:)
    else
      J1(:,1,:)=0d0
      J1(:,2,:)=0d0
    end if
  end if



  !ROTATE MAGNETIC FIELDS INTO VERTICAL,SOUTH,EAST COMPONENTS
  if (myid==0) then
    print *, 'magcalc.f90 --> Rotating currents into geomagnetic coordinates...'
  end if
  Jx=J1*proj_e1er+J2*proj_e2er+J3*proj_e3er                 !vertical
  Jy=J1*proj_e1etheta+J2*proj_e2etheta+J3*proj_e3etheta     !south
  Jz=J1*proj_e1ephi+J2*proj_e2ephi+J3*proj_e3ephi           !east

!    print *, myid2,myid3,'  --> Min/max values of current',minval(Jx),maxval(Jx),minval(Jy),maxval(Jy), &
!                                               minval(Jz),maxval(Jz)


  !GATHER THE END DATA SO WE DON'T LEAVE OUT A POINT IN THE INTEGRATION
  call halo_end(Jx,Jxend,Jxtop,tagJx)
  call halo_end(Jy,Jyend,Jytop,tagJy)
  call halo_end(Jz,Jzend,Jztop,tagJz)

    !print *, myid2,myid3,'  --> Min/max values of end current',minval(Jxend),maxval(Jxend),minval(Jyend),maxval(Jyend), &
    !                                           minval(Jzend),maxval(Jzend)
    !print *, myid2,myid3,'  --> Min/max values of top current',minval(Jxtop),maxval(Jxtop),minval(Jytop),maxval(Jytop), &
    !                                           minval(Jztop),maxval(Jztop)

  !COMPUTE MAGNETIC FIELDS
  do ipoints=1,lpoints
    if (myid == 0) then
      print *, 'magcalc.f90 --> Computing magnetic field for field point:  ',ipoints,' out of:  ',lpoints
      print *, '            --> ...for time:  ',ymd,UTsec
    end if


    Rx(:,:,:)=xf(ipoints)-xp(:,:,:)
    Ry(:,:,:)=yf(ipoints)-yp(:,:,:)
    Rz(:,:,:)=zf(ipoints)-zp(:,:,:)
    call halo_end(Rx,Rxend,Rxtop,tagRx)
    call halo_end(Ry,Ryend,Rytop,tagRy)
    call halo_end(Rz,Rzend,Rztop,tagRz)
!
!        do ix2=1,lx2
!          do ix1=1,lx1
!            if (ieee_is_nan(Rxend(ix1,ix2)) .or. abs(Rxend(ix1,ix2))<Rmin) then
!              if(myid2/=11) then
!                print*,'Rx end:  ',myid2,myid3,Rxend(ix1,ix2),ix1,ix2
!              end if
!            end if
!          end do
!        end do
!
!        do ix3=1,lx3
!          do ix1=1,lx1
!            if (ieee_is_nan(Rxtop(ix1,ix3)) .or. abs(Rxtop(ix1,ix3))<Rmin) then
!              if (myid3/=5) then
!                print*,'Rx top:  ',myid2,myid3,Rxtop(ix1,ix3),ix1,ix3
!              end if
!            end if
!          end do
!        end do

!
!    !enforce a regulator on the distance variables to avoid div by zero
!    where (Rx<Rmin .and. Rx>=0._wp)
!      Rx=Rmin
!    end where
!    where (Ry<Rmin .and. Ry>=0._wp)
!      Ry=Rmin
!    end where
!    where (Rz<Rmin .and. Rz>=0._wp)
!      Rz=Rmin
!    end where
!    where (Rx>-1*Rmin .and. Rx<0._wp)
!      Rx=-1*Rmin
!    end where
!    where (Ry>-1*Rmin .and. Ry<0._wp)
!      Ry=-1*Rmin
!    end where
!    where (Rz>-1*Rmin .and. Rz<0._wp)
!      Rz=-1*Rmin
!    end where
!
!    where (Rxend<Rmin .and. Rxend>=0._wp)
!      Rxend=Rmin
!    end where
!    where (Ryend<Rmin .and. Ryend>=0._wp)
!      Ryend=Rmin
!    end where
!    where (Rzend<Rmin .and. Rzend>=0._wp)
!      Rzend=Rmin
!    end where
!    where (Rxend>-1*Rmin .and. Rxend<0._wp)
!      Rxend=-1*Rmin
!    end where
!    where (Ryend>-1*Rmin .and. Ryend<0._wp)
!      Ryend=-1*Rmin
!    end where
!    where (Rzend>-1*Rmin .and. Rzend<0._wp)
!      Rzend=-1*Rmin
!    end where
!
!    where (Rxtop<Rmin .and. Rxtop>=0._wp)
!      Rxtop=Rmin
!    end where
!    where (Rytop<Rmin .and. Rytop>=0._wp)
!      Rytop=Rmin
!    end where
!    where (Rztop<Rmin .and. Rztop>=0._wp)
!      Rztop=Rmin
!    end where
!    where (Rxtop>-1*Rmin .and. Rxtop<0._wp)
!      Rxtop=-1*Rmin
!    end where
!    where (Rytop>-1*Rmin .and. Rytop<0._wp)
!      Rytop=-1*Rmin
!    end where
!    where (Rztop>-1*Rmin .and. Rztop<0._wp)
!      Rztop=-1*Rmin
!    end where

!    print*,myid,minval(abs(Rx)),minval(abs(Ry)),minval(abs(Rz)),minval(abs(Rxend)),minval(abs(Ryend)), &
!            minval(abs(Rzend)),minval(abs(Rxtop)),minval(abs(Rytop)),minval(abs(Rztop))

    if (flag2D/=1) then
!      Rcubed(:,:,:)=(Rx**2._wp+Ry**2._wp+Rz**2._wp)**(3._wp/2._wp)   !this really is R**3
      Rcubed(:,:,:)=(Rx**2._wp+Ry**2._wp+Rz**2._wp+Rmin**2._wp)**(3._wp/2._wp)   !this really is R**3

!      do ix3=1,lx3
!        do ix2=1,lx2
!          do ix1=1,lx1
!            if (ieee_is_nan(Rcubed(ix1,ix2,ix3)) .or. abs(Rcubed(ix1,ix2,ix3))<R3min) then
!              print*,myid2,myid3,Rcubed(ix1,ix2,ix3),ix1,ix2,ix3
!            end if
!          end do
!        end do
!      end do
!      print*,myid,minval(Rcubed)

!      where(Rcubed<R3min)
!        Rcubed=R3min
!      end where
      call halo_end(Rcubed,Rcubedend,Rcubedtop,tagRcubed)
      if(myid3==lid3-1) Rcubedend=R3min     !< avoids div by zero on the end which is set by the haloing
      if(myid2==lid2-1) Rcubedtop=R3min

        do ix2=1,lx2
          do ix1=1,lx1
            if (ieee_is_nan(Rcubedend(ix1,ix2)) .or. abs(Rcubedend(ix1,ix2))<R3min) then
!              if ( myid2/=11 ) then
                print*,'end:  ',myid2,myid3,Rcubedend(ix1,ix2),ix1,ix2
!              end if
            end if
          end do
        end do

        do ix3=1,lx3
          do ix1=1,lx1
            if (ieee_is_nan(Rcubedtop(ix1,ix3)) .or. abs(Rcubedtop(ix1,ix3))<R3min) then
!              if (myid3/=5) then
                print*,'top:  ',myid2,myid3,Rcubedtop(ix1,ix3),ix1,ix3
!              end if
            end if
          end do
        end do

!      where(Rcubedtop<R3min)
!        Rcubedtop=R3min
!      end where
!      where(Rcubedend<R3min)
!        Rcubedend=R3min
!      end where

      !print*, myid2,myid3,'--> Rcubed:  ',minval(Rcubed),maxval(Rcubed),minval(Rcubedend), &
      !                                   maxval(Rcubedend),minval(Rcubedtop),minval(Rcubedtop)

      !! FIXME: MAY BE MISSING A CORNER POINT HERE???  NO I THINK IT'S OKAY BASED ON SOME SQUARES I DREW, haha...


      !Bx calculation
      integrand(:,:,:)=mu0/4._wp/pi*(Jy*Rz-Jz*Ry)/Rcubed
      integrandend(:,:)=mu0/4._wp/pi*(Jyend*Rzend-Jzend*Ryend)/Rcubedend
      integrandtop(:,:)=mu0/4._wp/pi*(Jytop*Rztop-Jztop*Rytop)/Rcubedtop
      integrandavg(:,:,:)=1d0/8d0*( integrand(1:lx1-1,1:lx2-1,1:lx3-1) + integrand(2:lx1,1:lx2-1,1:lx3-1) + &
                           integrand(1:lx1-1,2:lx2,1:lx3-1) + integrand(2:lx1,2:lx2,1:lx3-1) + &
                           integrand(1:lx1-1,1:lx2-1,2:lx3) + integrand(2:lx1,1:lx2-1,2:lx3) + &
                           integrand(1:lx1-1,2:lx2,2:lx3) + integrand(2:lx1,2:lx2,2:lx3) )
      integrandavgend(:,:)=1d0/8d0*( integrand(1:lx1-1,1:lx2-1,lx3) + integrand(2:lx1,1:lx2-1,lx3) + &
                           integrand(1:lx1-1,2:lx2,lx3) + integrand(2:lx1,2:lx2,lx3) + &
                           integrandend(1:lx1-1,1:lx2-1) + integrandend(2:lx1,1:lx2-1) + &
                           integrandend(1:lx1-1,2:lx2) + integrandend(2:lx1,2:lx2) )
      integrandavgtop(:,:)=1d0/8d0*( integrand(1:lx1-1,lx2,1:lx3-1) + integrand(2:lx1,lx2,1:lx3-1) + &
                           integrand(1:lx1-1,lx2,2:lx3) + integrand(2:lx1,lx2,2:lx3) + &
                           integrandtop(1:lx1-1,1:lx3-1) + integrandtop(2:lx1,1:lx3-1) + &
                           integrandtop(1:lx1-1,2:lx3) + integrandtop(2:lx1,2:lx3) )
      Br(ipoints)=sum(integrandavg*dV(2:lx1,2:lx2,2:lx3))+sum(integrandavgend*dVend(2:lx1,2:lx2))+ &
                    sum(integrandavgtop*dVtop(2:lx1,2:lx3))

      !By
      integrand(:,:,:)=-mu0/4d0/pi*(Jx*Rz-Jz*Rx)/Rcubed
      integrandend(:,:)=-mu0/4d0/pi*(Jxend*Rzend-Jzend*Rxend)/Rcubedend
      integrandtop(:,:)=-mu0/4d0/pi*(Jxtop*Rztop-Jztop*Rxtop)/Rcubedtop
      integrandavg(:,:,:)=1d0/8d0*( integrand(1:lx1-1,1:lx2-1,1:lx3-1) + integrand(2:lx1,1:lx2-1,1:lx3-1) + &
                           integrand(1:lx1-1,2:lx2,1:lx3-1) + integrand(2:lx1,2:lx2,1:lx3-1) + &
                           integrand(1:lx1-1,1:lx2-1,2:lx3) + integrand(2:lx1,1:lx2-1,2:lx3) + &
                           integrand(1:lx1-1,2:lx2,2:lx3) + integrand(2:lx1,2:lx2,2:lx3) )
      integrandavgend(:,:)=1d0/8d0*( integrand(1:lx1-1,1:lx2-1,lx3) + integrand(2:lx1,1:lx2-1,lx3) + &
                           integrand(1:lx1-1,2:lx2,lx3) + integrand(2:lx1,2:lx2,lx3) + &
                           integrandend(1:lx1-1,1:lx2-1) + integrandend(2:lx1,1:lx2-1) + &
                           integrandend(1:lx1-1,2:lx2) + integrandend(2:lx1,2:lx2) )
      integrandavgtop(:,:)=1d0/8d0*( integrand(1:lx1-1,lx2,1:lx3-1) + integrand(2:lx1,lx2,1:lx3-1) + &
                           integrand(1:lx1-1,lx2,2:lx3) + integrand(2:lx1,lx2,2:lx3) + &
                           integrandtop(1:lx1-1,1:lx3-1) + integrandtop(2:lx1,1:lx3-1) + &
                           integrandtop(1:lx1-1,2:lx3) + integrandtop(2:lx1,2:lx3) )
      Btheta(ipoints)=sum(integrandavg*dV(2:lx1,2:lx2,2:lx3))+sum(integrandavgend*dVend(2:lx1,2:lx2))+ &
                        sum(integrandavgtop*dVtop(2:lx1,2:lx3))


      !Bz
      integrand(:,:,:)=mu0/4d0/pi*(Jx*Ry-Jy*Rx)/Rcubed
      integrandend(:,:)=mu0/4d0/pi*(Jxend*Ryend-Jyend*Rxend)/Rcubedend
      integrandtop(:,:)=mu0/4d0/pi*(Jxtop*Rytop-Jytop*Rxtop)/Rcubedtop
      integrandavg(:,:,:)=1d0/8d0*( integrand(1:lx1-1,1:lx2-1,1:lx3-1) + integrand(2:lx1,1:lx2-1,1:lx3-1) + &
                           integrand(1:lx1-1,2:lx2,1:lx3-1) + integrand(2:lx1,2:lx2,1:lx3-1) + &
                           integrand(1:lx1-1,1:lx2-1,2:lx3) + integrand(2:lx1,1:lx2-1,2:lx3) + &
                           integrand(1:lx1-1,2:lx2,2:lx3) + integrand(2:lx1,2:lx2,2:lx3) )
      integrandavgend(:,:)=1d0/8d0*( integrand(1:lx1-1,1:lx2-1,lx3) + integrand(2:lx1,1:lx2-1,lx3) + &
                           integrand(1:lx1-1,2:lx2,lx3) + integrand(2:lx1,2:lx2,lx3) + &
                           integrandend(1:lx1-1,1:lx2-1) + integrandend(2:lx1,1:lx2-1) + &
                           integrandend(1:lx1-1,2:lx2) + integrandend(2:lx1,2:lx2) )
      integrandavgtop(:,:)=1d0/8d0*( integrand(1:lx1-1,lx2,1:lx3-1) + integrand(2:lx1,lx2,1:lx3-1) + &
                           integrand(1:lx1-1,lx2,2:lx3) + integrand(2:lx1,lx2,2:lx3) + &
                           integrandtop(1:lx1-1,1:lx3-1) + integrandtop(2:lx1,1:lx3-1) + &
                           integrandtop(1:lx1-1,2:lx3) + integrandtop(2:lx1,2:lx3) )
      Bphi(ipoints)=sum(integrandavg*dV(2:lx1,2:lx2,2:lx3))+sum(integrandavgend*dVend(2:lx1,2:lx2))+ &
                        sum(integrandavgtop*dVtop(2:lx1,2:lx3))
    else
      Rcubed(:,:,:)=Rx**2+Ry**2    !not really R**3 in 2D, just the denominator of the integrand
!      where(Rcubed<R3min)
!        Rcubed=R3min    !should be R**2???
!      end where
      call halo_end(Rcubed,Rcubedend,Rcubedtop,tagRcubed)
      !! DO WE NEED TO CHECK HERE FOR DIV BY ZERO???
      !! ALSO IN 2D WE KNOW THAT WE ARE ONLY DIVIDED IN THE 3 DIMENSION SO THERE IS NO NEED TO WORRY ABOUT ADDING A 'TOP' ETC.

      !Bx
      integrand(:,:,:)=mu0/4d0/pi*(-2d0*Jz*Ry)/Rcubed
      integrandend(:,:)=mu0/4d0/pi*(-2d0*Jzend*Ryend)/Rcubedend
      integrandavg(:,:,:)=1d0/4d0*( integrand(1:lx1-1,:,1:lx3-1) + integrand(2:lx1,:,1:lx3-1) + &
                           integrand(1:lx1-1,:,2:lx3) + integrand(2:lx1,:,2:lx3) )
      integrandavgend(:,:)=1d0/4d0*( integrand(1:lx1-1,:,lx3) + integrand(2:lx1,:,lx3) + &
                           integrandend(1:lx1-1,:) + integrandend(2:lx1,:) )
      Br(ipoints)=sum(integrandavg*dV(2:lx1,:,2:lx3))+sum(integrandavgend*dVend(2:lx1,:))

      !By
      integrand(:,:,:)=mu0/4/pi*(2*Jz*Rx)/Rcubed
      integrandend(:,:)=mu0/4/pi*(2*Jzend*Rxend)/Rcubedend
      integrandavg(:,:,:)=1d0/4d0*( integrand(1:lx1-1,:,1:lx3-1) + integrand(2:lx1,:,1:lx3-1) + &
                           integrand(1:lx1-1,:,2:lx3) + integrand(2:lx1,:,2:lx3) )
      integrandavgend(:,:)=1d0/4d0*( integrand(1:lx1-1,:,lx3) + integrand(2:lx1,:,lx3) + &
                           integrandend(1:lx1-1,:) + integrandend(2:lx1,:) )
      Btheta(ipoints) = sum(integrandavg*dV(2:lx1,:,2:lx3))+sum(integrandavgend*dVend(2:lx1,:))
      !! without dim= input it just sums everything which is what we want

      !Bz
      integrand(:,:,:)=mu0/4d0/pi*2d0*(Jx*Ry-Jy*Rx)/Rcubed
      integrandend(:,:)=mu0/4d0/pi*2d0*(Jxend*Ryend-Jyend*Rxend)/Rcubedend
      integrandavg(:,:,:)=1d0/4d0*( integrand(1:lx1-1,:,1:lx3-1) + integrand(2:lx1,:,1:lx3-1) + &
                           integrand(1:lx1-1,:,2:lx3) + integrand(2:lx1,:,2:lx3) )
      integrandavgend(:,:)=1d0/4d0*( integrand(1:lx1-1,:,lx3) + integrand(2:lx1,:,lx3) + &
                           integrandend(1:lx1-1,:) + integrandend(2:lx1,:) )
      Bphi(ipoints) = sum(integrandavg*dV(2:lx1,:,2:lx3))+sum(integrandavgend*dVend(2:lx1,:))
      !! without dim= input it just sums everything which is what we want
    end if
  end do
    !print *, myid2,myid3,'  --> Min/max values of field',minval(Br),maxval(Br),minval(Btheta),maxval(Btheta), &
    !                                           minval(Bphi),maxval(Bphi)


  !A REDUCE OPERATION IS NEEDED HERE TO COMBINE MAGNETIC FIELDS (LINEAR SUPERPOSITION) FROM ALL WORKERS
  if (myid ==0) then
    if(debug) print *, 'Attempting reduction of magnetic field...'
  end if
  call mpi_reduce(Br,Brall,lpoints,mpi_realprec,MPI_SUM,0,MPI_COMM_WORLD,ierr)
  call mpi_reduce(Btheta,Bthetaall,lpoints,mpi_realprec,MPI_SUM,0,MPI_COMM_WORLD,ierr)
  call mpi_reduce(Bphi,Bphiall,lpoints,mpi_realprec,MPI_SUM,0,MPI_COMM_WORLD,ierr)
  if (myid == 0) then
    if(debug) print *, 'magcalc.f90 --> Reduced magnetic field...'
    if(debug) print *, '  --> Min/max values of reduced field',minval(Brall),maxval(Brall),minval(Bthetaall),maxval(Bthetaall), &
                                               minval(Bphiall),maxval(Bphiall)
  end if


  !OUTPUT SHOULD BE DONE FOR EVERY INPUT FILE THAT HAS BEEN READ IN
  if (myid==0) then
    call cpu_time(tstart)
    call output_magfields(outdir,ymd,UTsec,Brall,Bthetaall,Bphiall)   !mag field data already reduced so just root needs to output
    call cpu_time(tfin)
   if(debug) print *, 'magcalc.f90 --> Output done for time step:  ',t,' in cpu_time of:  ',tfin-tstart
  end if
  tout=tout+dtout


  !NOW OUR SOLUTION IS FULLY UPDATED SO UPDATE TIME VARIABLES TO MATCH...
  it=it+1; t=t+dt;
  if (myid==0 .and. debug) print *, 'magcalc.f90 --> Moving on to time step (in sec):  ',t,'; end time of simulation:  ',tdur
  call dateinc(dt,ymd,UTsec)
  if (myid==0) then
    print *, 'magcalc.f90 --> Current date',ymd,'Current UT time:  ',UTsec
  end if
end do


!! DEALLOCATE MAIN PROGRAM DATA
deallocate(J1,J2,J3)
deallocate(Jx,Jy,Jz)
deallocate(xp,yp,zp)
deallocate(xf,yf,zf)
deallocate(dV)
deallocate(Rx,Ry,Rz)
deallocate(proj_e1er,proj_e2er,proj_e3er)
deallocate(proj_e1etheta,proj_e2etheta,proj_e3etheta)
deallocate(proj_e1ephi,proj_e2ephi,proj_e3ephi)
deallocate(Rcubed)
deallocate(Rcubedend)
deallocate(Rcubedtop)
deallocate(integrand,integrandavg)
deallocate(Br,Btheta,Bphi)
deallocate(dVend,Jxend,Jyend,Jzend,Rxend,Ryend,Rzend)
deallocate(integrandend,integrandavgend)
deallocate(dVtop,Jxtop,Jytop,Jztop,Rxtop,Rytop,Rztop)
deallocate(integrandtop,integrandavgtop)

!! SHUT DOWN MPI
ierr = mpibreakdown()

if (ierr /= 0) then
  write(stderr, *) 'MAGCALC: abnormal MPI shutdown code', ierr
  error stop
endif

print '(/,A)', 'MAGCALC: simulation complete'

end program