read.f90 Source File

submodule (grid) grid_read

implicit none

interface ! readgrid_{raw,hdf5,nc4}.f90
module subroutine get_grid3(path, flagswap, x, g1all,g2all,g3all, altall,glatall,glonall,Bmagall,Incall,nullptsall,&
  e1all,e2all,e3all,erall,ethetaall,ephiall,rall,thetaall,phiall)
character(*), intent(in) :: path
integer, intent(in) :: flagswap
type(curvmesh), intent(inout) :: x
real(wp), intent(out), dimension(:,:,:) :: g1all,g2all,g3all,altall,glatall,glonall,Bmagall,nullptsall,rall,thetaall,phiall
real(wp), intent(out), dimension(:,:) :: Incall
real(wp), intent(out), dimension(:,:,:,:) :: e1all,e2all,e3all,erall,ethetaall,ephiall
end subroutine get_grid3
end interface

contains


module procedure read_grid
! subroutine read_grid(indatsize,indatgrid,flagperiodic,x)
!! READS GRID INFORMATION FROM A BBINARY FILE AND ALLOCATES
!! GRID STRUCTURE.  NOTE THAT THE INPUT FILE DOES NOT,
!! BY DEFAULT INCLUDE THE GHOST CELLS.  THIS FUNCTION
!! IS A WRAPPER WHICH PASSES WORK OFF TO APPROPRIATE
!! SUBROUTINE DEPENDING ON WHETHER IT IS CALLED BY ROOT
!! OR NOT.   THIS SUBRoutINE ALSO CLASSIFIES THE GRID.



!> READ IN THE GRID DATA
if(myid==0) then
  call read_grid_root(indatsize,indatgrid,x)
else
  call read_grid_workers(x)
end if


!FLAG THE GRID AS PERIODIC, IF REQUESTED; IF SO PERIODICITY WILL BE ASSUMED
!IN THE X3-DIRECTION, NOTE BOTH ROOT AND WORKERS MUST DO THIS!!!
if (flagperiodic==1) then
  x%flagper=.true.
else
  x%flagper=.false.
end if


!DETERMINE THE TYPE OF GRID WE HAVE AND SET AN APPROPRIATE FLAG
!FIXME:  this needs to be done once, globally and also needs to be more
!robust...
if (abs(x%alt(1,1,1)-x%alt(lx1,1,1))<100d3) then    !closed dipole grid
  gridflag=0
else if (x%alt(1,1,1)>x%alt(2,1,1)) then    !open dipole grid with inverted structure wrt altitude
  gridflag=1
else    !something different (viz. non-inverted - lowest altitudes at the logical bottom of the grid)
  gridflag=2
end if

end procedure read_grid


subroutine read_grid_root(indatsize,indatgrid,x)

!------------------------------------------------------------
!--------SOME CODE DUPLICATION WITH WORKER VERSION - CAN WE
!--------CREATE A COMMON SUBROUTINE THAT ALLOCATES SOME VARS?
!--------THE GRID SIZES MUST ALREADY BE DEFINED IN TEH MODULE
!--------, PROBABLY FROM CALLING GRID_SIZE
!------------------------------------------------------------

character(*), intent(in) :: indatsize,indatgrid
type(curvmesh), intent(inout) :: x    !does this need to be inout?  I think if anything is unallocated, it does...

integer u,ierr,iid,ix1,ix2,ix3,icount,icomp!,itell

!NOTE THAT HAVING THESE ARE LOCAL (TEMPORARY) VARS. PREVENTS ROOT FROM WRITING ENTIRE GRID TO FILE AT SOME LATER POINT...
real(wp), dimension(:,:,:), allocatable :: g1all,g2all,g3all   !to temporarily store input data to be distributed
real(wp), dimension(:,:,:), allocatable :: altall,glatall,glonall
real(wp), dimension(:,:,:), allocatable :: rall,thetaall,phiall
real(wp), dimension(:,:,:), allocatable :: Bmagall
real(wp), dimension(:,:), allocatable :: Incall
real(wp), dimension(:,:,:), allocatable :: nullptsall
real(wp), dimension(:,:,:,:), allocatable :: e1all,e2all,e3all,erall,ethetaall,ephiall    !might be best to have a tmp vector array...
real(wp), dimension(:,:,:), allocatable :: mpisendbuf
real(wp), dimension(:,:,:), allocatable :: mpirecvbuf

!print*, 'Entering read_grid_root', lx1,lx2all,lx3all,lid2,lid3,lx2all/lid2


!DETERMINE NUMBER OF SLABS AND CORRESPONDING SIZE FOR EACH WORKER
!NOTE THAT WE WILL ASSUME THAT THE GRID SIZE IS DIVISIBLE BY NUMBER OF WORKERS AS THIS HAS BEEN CHECKED A FEW LINES BELOW
x%lx1=lx1; x%lx2all=lx2all; x%lx3all=lx3all;


!> ADJUST THE SIZES OF THE VARIABLES IF LX3ALL==1, SO THAT THE ALLOCATIONS ARE THE CORRECT SIZE
if (lx3all==1) then
  print *, '2D run: **SWAP** x2 and x3 dims'
  lx3all=lx2all; x%lx3all=lx2all;
  lx2=1
  lx2all=1; x%lx2all=1;
  lx3=lx3all/lid
  !! use total number of processes since we only divide in one direction here...
  flagswap=1
else
  if(lx2all==1) then
    print *, '2D run: do not swap x2 and x3 dims.'
    lx2=1
    lx3=lx3all/lid
  else
    print *, '3D run'
    lx2=lx2all/lid2
    !! should divide evenly if generated from mpigrid
    lx3=lx3all/lid3
  endif
  flagswap=0
end if
x%lx2=lx2; x%lx3=lx3
print '(A,3I6)', 'Grid slab size:  ',lx1,lx2,lx3


!> COMMUNICATE THE GRID SIZE TO THE WORKERS SO THAT THEY CAN ALLOCATE SPACE
ierr=0
!! if this is a root-only simulation we don't want to error out
do iid=1,lid-1
  call mpi_send(lx2,1,MPI_INTEGER,iid,taglx2,MPI_COMM_WORLD,ierr)
  !! need to also pass the lx2all size to all workers to they know
  !if (ierr/=0) error stop 'grid:read_grid_root failed mpi_send lx2'

  call mpi_send(lx3,1,MPI_INTEGER,iid,taglx3,MPI_COMM_WORLD,ierr)
  !if (ierr/=0) error stop 'grid:read_grid_root failed mpi_send lx3'
end do

if (ierr/=0) error stop 'grid:read_grid_root failed mpi_send grid size'

!TELL WORKERS IF WE'VE SWAPPED DIMENSIONS
ierr=0
do iid=1,lid-1
  call mpi_send(flagswap,1,MPI_INTEGER,iid,tagswap,MPI_COMM_WORLD,ierr)
  !if (ierr/=0) error stop 'grid:read_grid_root failed mpi_send flagswap'
end do

if (ierr/=0) error stop 'grid:read_grid_root failed mpi_send flagswap'

!ALLOCATE SPACE FOR ROOTS SLAB OF DATA
allocate(x%x1(-1:lx1+2))
allocate(x%dx1i(lx1),x%x1i(lx1+1),x%dx1(0:lx1+2))


!FULL GRID X2 VARIABLE
allocate(x%x2all(-1:lx2all+2))
allocate(x%dx2all(0:lx2all+2))
allocate(x%dx2iall(lx2all),x%x2iall(lx2all+1))


!FULL-GRID X3-VARIABLE
allocate(x%x3all(-1:lx3all+2))
allocate(x%dx3all(0:lx3all+2))
allocate(x%x3iall(lx3all+1),x%dx3iall(lx3all))

allocate(x%h1all(-1:lx1+2,-1:lx2all+2,-1:lx3all+2),x%h2all(-1:lx1+2,-1:lx2all+2,-1:lx3all+2), &
         x%h3all(-1:lx1+2,-1:lx2all+2,-1:lx3all+2))    !do we need the ghost cell values?  Yes the divergence in compression term is computed using one ghost cell in each direction!!!!
allocate(x%h1x1iall(1:lx1+1,1:lx2all,1:lx3all),x%h2x1iall(1:lx1+1,1:lx2all,1:lx3all), &
           x%h3x1iall(1:lx1+1,1:lx2all,1:lx3all))
allocate(x%h1x2iall(1:lx1,1:lx2all+1,1:lx3all),x%h2x2iall(1:lx1,1:lx2all+1,1:lx3all), &
           x%h3x2iall(1:lx1,1:lx2all+1,1:lx3all))
allocate(x%h1x3iall(1:lx1,1:lx2all,1:lx3all+1),x%h2x3iall(1:lx1,1:lx2all,1:lx3all+1), &
           x%h3x3iall(1:lx1,1:lx2all,1:lx3all+1))

allocate(x%rall(1:lx1,1:lx2all,1:lx3all),x%thetaall(1:lx1,1:lx2all,1:lx3all),x%phiall(1:lx1,1:lx2all,1:lx3all))


!DETERMINE AND ALLOCATE SPACE NEEDED FOR ROOT SUBGRIDS (WORKERS USE SIMILAR ALLOCATE STATEMENTS)
allocate(x%x3(-1:lx3+2))
allocate(x%dx3i(lx3),x%x3i(lx3+1),x%dx3(0:lx3+2))

allocate(x%x2(-1:lx2+2))
allocate(x%dx2i(lx2),x%x2i(lx2+1),x%dx2(0:lx2+2))

allocate(x%h1(-1:lx1+2,-1:lx2+2,-1:lx3+2),x%h2(-1:lx1+2,-1:lx2+2,-1:lx3+2),x%h3(-1:lx1+2,-1:lx2+2,-1:lx3+2))
!    allocate(x%h1(1:lx1,1:lx2,1:lx3),x%h2(1:lx1,1:lx2,1:lx3),x%h3(1:lx1,1:lx2,1:lx3))
allocate(x%h1x1i(1:lx1+1,1:lx2,1:lx3),x%h2x1i(1:lx1+1,1:lx2,1:lx3),x%h3x1i(1:lx1+1,1:lx2,1:lx3))
allocate(x%h1x2i(1:lx1,1:lx2+1,1:lx3),x%h2x2i(1:lx1,1:lx2+1,1:lx3),x%h3x2i(1:lx1,1:lx2+1,1:lx3))
allocate(x%h1x3i(1:lx1,1:lx2,1:lx3+1),x%h2x3i(1:lx1,1:lx2,1:lx3+1),x%h3x3i(1:lx1,1:lx2,1:lx3+1))

allocate(x%glat(1:lx1,1:lx2,1:lx3),x%glon(1:lx1,1:lx2,1:lx3),x%alt(1:lx1,1:lx2,1:lx3))
allocate(x%r(1:lx1,1:lx2,1:lx3),x%theta(1:lx1,1:lx2,1:lx3),x%phi(1:lx1,1:lx2,1:lx3))

allocate(x%Bmag(1:lx1,1:lx2,1:lx3))
allocate(x%I(1:lx2,1:lx3))
allocate(x%nullpts(1:lx1,1:lx2,1:lx3))

allocate(x%e1(1:lx1,1:lx2,1:lx3,1:3),x%e2(1:lx1,1:lx2,1:lx3,1:3),x%e3(1:lx1,1:lx2,1:lx3,1:3))
allocate(x%er(1:lx1,1:lx2,1:lx3,1:3),x%etheta(1:lx1,1:lx2,1:lx3,1:3),x%ephi(1:lx1,1:lx2,1:lx3,1:3))


!NOW WE NEED TO READ IN THE FULL GRID DATA AND PUT IT INTO THE STRUCTURE.
!IF WE HAVE DONE ANY DIMENSION SWAPPING HERE WE NEED TO TAKE THAT INTO ACCOUNT IN THE VARIABLES THAT ARE BEING READ IN
print *, 'Starting grid input from file: ',indatgrid

allocate(g1all(lx1,lx2all,lx3all), g2all(lx1,lx2all,lx3all), g3all(lx1,lx2all,lx3all))
allocate(altall(lx1,lx2all,lx3all), glatall(lx1,lx2all,lx3all), glonall(lx1,lx2all,lx3all))
allocate(Bmagall(lx1,lx2all,lx3all))
allocate(Incall(lx2all,lx3all))
allocate(nullptsall(lx1,lx2all,lx3all))
allocate(e1all(lx1,lx2all,lx3all,3), e2all(lx1,lx2all,lx3all,3), e3all(lx1,lx2all,lx3all,3))
allocate(erall(lx1,lx2all,lx3all,3), ethetaall(lx1,lx2all,lx3all,3), ephiall(lx1,lx2all,lx3all,3))
allocate(rall(lx1,lx2all,lx3all), thetaall(lx1,lx2all,lx3all), phiall(lx1,lx2all,lx3all))

if (flagswap/=1) then
  !! normal (i.e. full 3D) grid ordering, or a 2D grid with 1 element naturally in the second dimension
else
  !! this is apparently a 2D grid, so the x2 and x3 dimensions have been/need to be swapped
  !! MZ - may need to change lx2-->lx2all???
  !! Note there that the fortran arrays are the correct size, but the input data are not!!!
  !! This means tmp variable and permutes...

  print *, '2D grid: **PERMUTE** x2 and x3 dimensions'
end if

call get_grid3(indatgrid, flagswap, x, g1all,g2all,g3all, altall,glatall,glonall,Bmagall,Incall,nullptsall,&
  e1all,e2all,e3all,erall,ethetaall,ephiall,rall,thetaall,phiall)

!! ALLOCATE SPACE FOR ROOTS SUBGRID GRAVITATIONAL FIELD
allocate(g1(1:lx1,1:lx2,1:lx3),g2(1:lx1,1:lx2,1:lx3),g3(1:lx1,1:lx2,1:lx3))


!! SEND FULL X1 AND X2 GRIDS TO EACH WORKER (ONLY X3-DIM. IS INVOLVED IN THE MPI
print*, 'Exchanging grid spacings...'
do iid=1,lid-1
  call mpi_send(x%x1,lx1+4,mpi_realprec,iid,tagx1,MPI_COMM_WORLD,ierr)
  !if (ierr/=0) error stop 'failed mpi_send x1'
  call mpi_send(x%x2all,lx2all+4,mpi_realprec,iid,tagx2all,MPI_COMM_WORLD,ierr)
  !if (ierr/=0) error stop 'failed mpi_send x2all'
  call mpi_send(x%x3all,lx3all+4,mpi_realprec,iid,tagx3all,MPI_COMM_WORLD,ierr)
  !! workers may need a copy of this, e.g. for boudnary conditions
  !if (ierr/=0) error stop 'failed mpi_send x3all'

  call mpi_send(x%dx1,lx1+3,mpi_realprec,iid,tagx1,MPI_COMM_WORLD,ierr)
  !if (ierr/=0) error stop 'failed mpi_send dx1'
  call mpi_send(x%x1i,lx1+1,mpi_realprec,iid,tagx1,MPI_COMM_WORLD,ierr)
  !if (ierr/=0) error stop 'failed mpi_send x1i'
  call mpi_send(x%dx1i,lx1,mpi_realprec,iid,tagx1,MPI_COMM_WORLD,ierr)
  !if (ierr/=0) error stop 'failed mpi_send dx1i'
end do

if (ierr/=0) error stop 'failed mpi_send x grid'

!! NOW SEND THE INFO THAT DEPENDS ON X3 SLAB SIZE
print*, 'Computing subdomain spacing...'
call bcast_send1D_3(x%x3all,tagx3,x%x3)
x%dx3=x%x3(0:lx3+2)-x%x3(-1:lx3+1)     !computing these avoids extra message passing (could be done for other coordinates, as well)
x%x3i(1:lx3+1)=0.5*(x%x3(0:lx3)+x%x3(1:lx3+1))
x%dx3i=x%x3i(2:lx3+1)-x%x3i(1:lx3)


call bcast_send1D_2(x%x2all,tagx2,x%x2)
x%dx2=x%x2(0:lx2+2)-x%x2(-1:lx2+1)     !computing these avoids extra message passing (could be done for other coordinates
x%x2i(1:lx2+1)=0.5*(x%x2(0:lx2)+x%x2(1:lx2+1))
x%dx2i=x%x2i(2:lx2+1)-x%x2i(1:lx2)


print*, 'Dealing with metric factors...'
allocate(mpisendbuf(-1:lx1+2,-1:lx2all+2,-1:lx3all+2),mpirecvbuf(-1:lx1+2,-1:lx2+2,-1:lx3+2))

mpisendbuf=x%h1all    !since metric factors are pointers they are not gauranteed to be contiguous in memory so pack them into a buffer that is...
call bcast_send3D_ghost(mpisendbuf,tagh1,mpirecvbuf)    !special broadcast subroutine to handle 3D arrays with ghost cells
x%h1=mpirecvbuf       !store roots slab of metric factors in its grid structure
mpisendbuf=x%h2all
call bcast_send3D_ghost(mpisendbuf,tagh2,mpirecvbuf)
x%h2=mpirecvbuf
mpisendbuf=x%h3all
call bcast_send3D_ghost(mpisendbuf,tagh3,mpirecvbuf)
x%h3=mpirecvbuf

deallocate(mpisendbuf,mpirecvbuf)    !we need different sized buffers below

call bcast_send(x%h1x1iall,tagh1,x%h1x1i)    !do the weird sizes here (ie. lx1+1) give issues with MPI?  probably...  No because bcast reads the size off of the variable...
call bcast_send(x%h2x1iall,tagh2,x%h2x1i)
call bcast_send(x%h3x1iall,tagh3,x%h3x1i)

call bcast_send3D_x2i(x%h1x2iall,tagh1,x%h1x2i)
call bcast_send3D_x2i(x%h2x2iall,tagh2,x%h2x2i)
call bcast_send3D_x2i(x%h3x2iall,tagh3,x%h3x2i)

call bcast_send3D_x3i(x%h1x3iall,tagh1,x%h1x3i)
call bcast_send3D_x3i(x%h2x3iall,tagh2,x%h2x3i)
call bcast_send3D_x3i(x%h3x3iall,tagh3,x%h3x3i)


print *, 'Sending gravity, etc...'
call bcast_send(g1all,tagh1,g1)
call bcast_send(g2all,tagh2,g2)
call bcast_send(g3all,tagh3,g3)

call bcast_send(glatall,tagglat,x%glat)
call bcast_send(glonall,tagglon,x%glon)
call bcast_send(altall,tagalt,x%alt)

call bcast_send(Bmagall,tagBmag,x%Bmag)
call bcast_send(Incall,taginc,x%I)
call bcast_send(nullptsall,tagnull,x%nullpts)

print *, 'Now sending unit vectors...'
allocate(mpisendbuf(1:lx1,1:lx2all,1:lx3all),mpirecvbuf(1:lx1,1:lx2,1:lx3))    !why is buffering used/needed here???
do icomp=1,3
  mpisendbuf=e1all(:,:,:,icomp)
  call bcast_send(mpisendbuf,tageunit1,mpirecvbuf)
  x%e1(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  mpisendbuf=e2all(:,:,:,icomp)
  call bcast_send(mpisendbuf,tageunit2,mpirecvbuf)
  x%e2(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  mpisendbuf=e3all(:,:,:,icomp)
  call bcast_send(mpisendbuf,tageunit3,mpirecvbuf)
  x%e3(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  mpisendbuf=erall(:,:,:,icomp)
  call bcast_send(mpisendbuf,tager,mpirecvbuf)
  x%er(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  mpisendbuf=ethetaall(:,:,:,icomp)
  call bcast_send(mpisendbuf,tagetheta,mpirecvbuf)
  x%etheta(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  mpisendbuf=ephiall(:,:,:,icomp)
  call bcast_send(mpisendbuf,tagephi,mpirecvbuf)
  x%ephi(:,:,:,icomp)=mpirecvbuf
end do
deallocate(mpisendbuf,mpirecvbuf)

call bcast_send(rall,tagr,x%r)
call bcast_send(thetaall,tagtheta,x%theta)
call bcast_send(phiall,tagphi,x%phi)
print *,  'Done sending slabbed variables to workers...'


!COUNT THE NUMBER OF NULL GRID POINTS AND GENERATE A LIST OF NULL INDICES FOR LATER USE
x%lnull=0;
do ix3=1,lx3
  do ix2=1,lx2
    do ix1=1,lx1
      if(x%nullpts(ix1,ix2,ix3) > 0.5_wp) then
        x%lnull=x%lnull+1
      end if
    end do
  end do
end do
allocate(x%inull(1:x%lnull,1:3))

icount=1
do ix3=1,lx3
  do ix2=1,lx2
    do ix1=1,lx1
      if(x%nullpts(ix1,ix2,ix3) > 0.5_wp) then
        x%inull(icount,:)=[ix1,ix2,ix3]
        icount=icount+1
      end if
    end do
  end do
end do
print *, 'Done computing null grid points...  Process:  ',myid,' has:  ',x%lnull


!COMPUTE DIFFERENTIAL DISTANCES ALONG EACH DIRECTION (TO BE USED IN TIME STEP DETERMINATION...
block
real(wp), dimension(1:lx1,1:lx2,1:lx3) :: tmpdx
allocate(x%dl1i(1:lx1,1:lx2,1:lx3),x%dl2i(1:lx1,1:lx2,1:lx3),x%dl3i(1:lx1,1:lx2,1:lx3))

tmpdx=spread(spread(x%dx1i,2,lx2),3,lx3)
x%dl1i=tmpdx*x%h1(1:lx1,1:lx2,1:lx3)
tmpdx=spread(spread(x%dx2i,1,lx1),3,lx3)
x%dl2i=tmpdx*x%h2(1:lx1,1:lx2,1:lx3)
tmpdx=spread(spread(x%dx3i,1,lx1),2,lx2)
x%dl3i=tmpdx*x%h3(1:lx1,1:lx2,1:lx3)
end block

!    !THIS CODE BLOCK HAS ROOT "PARROT" THE GRID TO A FILE FOR DEBUGGING PURPOSES.
!    !THIS BLOCK MUST BE HERE DUE TO THE FACT MANY OF HTE ALL-GRID VARIABLES ARE ONLY
!    !IN SCOPE INSIDE THIS FUNCITON BEFORE THE DEALLOCATE STATEMENT BELOW
!    open(newunit=outunit,file='testsize.dat',status='replace',form='unformatted', &
!         access='stream', action='write')
!    write(outunit) lx1,lx2,lx3all    !note that these sizes exxclude ghost cells
!    close(outunit)
!
!    open(outunit,file='testgrid.dat',status='replace',form='unformatted', &
!         access='stream')
!    write(outunit) x%x1,x%x1i,x%dx1,x%dx1i    !I guess the number of elements is obtained from the size of the arrays???
!    write(outunit) x%x2,x%x2i,x%dx2,x%dx2i
!    write(outunit) x%x3all,x%x3iall,x%dx3all,x%dx3iall
!    write(outunit) x%h1all,x%h2all,x%h3all
!    write(outunit) x%h1x1iall,x%h2x1iall,x%h3x1iall
!    write(outunit) x%h1x2iall,x%h2x2iall,x%h3x2iall
!    write(outunit) x%h1x3iall,x%h2x3iall,x%h3x3iall
!    write(outunit) g1all,g2all,g3all
!    write(outunit) altall,glatall,glonall
!    write(outunit) Bmagall
!    write(outunit) Incall
!    write(outunit) nullptsall
!    close(outunit)
!    print *, 'Done creating copy of grid...'
!


!    !FLAG THE GRID AS PERIODIC, IF REQUESTED; IF SO PERIODICITY WILL BE ASSUMED IN THE X3-DIRECTION
!    if (flagperiodic==1) then
!      x%flagper=1
!    else
!      x%flagper=0
!    end if
!

!DEALLOCATE ANY FULL GRID VARIABLES THAT ARE NO LONGER NEEDED
deallocate(g1all,g2all,g3all,altall,glatall,glonall,Bmagall,Incall,nullptsall,rall,thetaall,phiall)

end subroutine read_grid_root


subroutine read_grid_workers(x)

type(curvmesh), intent(inout) :: x

integer :: ix1,ix2,ix3,icount,icomp, ierr

!GET ROOTS MESSAGE WITH THE SIZE OF THE GRID WE ARE TO RECEIVE
call mpi_recv(lx2,1,MPI_INTEGER,0,taglx2,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv lx2'
call mpi_recv(lx3,1,MPI_INTEGER,0,taglx3,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv lx3'

x%lx1=lx1; x%lx2=lx2; x%lx2all=lx2all; x%lx3all=lx3all; x%lx3=lx3

!ROOT NEEDS TO TELL US WHETHER WE'VE SWAPPED DIMENSIONS SINCE THIS AFFECTS HOW CURRENTS ARE COMPUTED
call mpi_recv(flagswap,1,MPI_INTEGER,0,tagswap,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv flagswap'

if (flagswap==1) then
  x%lx3all=lx2all
  x%lx2all=lx3all
  lx2all=x%lx2all
  lx3all=x%lx3all
end if


!ALLOCATE SPACE FOR MY SLAB OF DATA
allocate(x%x1(-1:lx1+2))
allocate(x%dx1i(lx1),x%x1i(lx1+1),x%dx1(0:lx1+2))

allocate(x%x2(-1:lx2+2))
allocate(x%dx2i(lx2),x%x2i(lx2+1),x%dx2(0:lx2+2))
allocate(x%x2all(-1:lx2all+2))

!DETERMINE AND ALLOCATE SPACE NEEDED FOR WORKERS SUBGRIDS
allocate(x%x3(-1:lx3+2))
allocate(x%dx3i(lx3),x%x3i(lx3+1),x%dx3(0:lx3+2))
allocate(x%x3all(-1:lx3all+2))

allocate(x%h1(-1:lx1+2,-1:lx2+2,-1:lx3+2),x%h2(-1:lx1+2,-1:lx2+2,-1:lx3+2),x%h3(-1:lx1+2,-1:lx2+2,-1:lx3+2))
allocate(x%h1x1i(1:lx1+1,1:lx2,1:lx3),x%h2x1i(1:lx1+1,1:lx2,1:lx3),x%h3x1i(1:lx1+1,1:lx2,1:lx3))
allocate(x%h1x2i(1:lx1,1:lx2+1,1:lx3),x%h2x2i(1:lx1,1:lx2+1,1:lx3),x%h3x2i(1:lx1,1:lx2+1,1:lx3))
allocate(x%h1x3i(1:lx1,1:lx2,1:lx3+1),x%h2x3i(1:lx1,1:lx2,1:lx3+1),x%h3x3i(1:lx1,1:lx2,1:lx3+1))

allocate(x%glat(1:lx1,1:lx2,1:lx3),x%glon(1:lx1,1:lx2,1:lx3),x%alt(1:lx1,1:lx2,1:lx3))
allocate(x%r(1:lx1,1:lx2,1:lx3),x%theta(1:lx1,1:lx2,1:lx3),x%phi(1:lx1,1:lx2,1:lx3))

allocate(x%Bmag(1:lx1,1:lx2,1:lx3))
allocate(x%I(1:lx2,1:lx3))
allocate(x%nullpts(1:lx1,1:lx2,1:lx3))

allocate(x%e1(1:lx1,1:lx2,1:lx3,1:3),x%e2(1:lx1,1:lx2,1:lx3,1:3),x%e3(1:lx1,1:lx2,1:lx3,1:3))
allocate(x%er(1:lx1,1:lx2,1:lx3,1:3),x%etheta(1:lx1,1:lx2,1:lx3,1:3),x%ephi(1:lx1,1:lx2,1:lx3,1:3))


!ALLOCATE SPACE FOR WORKER'S GRAVITATIONAL FIELD
allocate(g1(1:lx1,1:lx2,1:lx3),g2(1:lx1,1:lx2,1:lx3),g3(1:lx1,1:lx2,1:lx3))


!RECEIVE GRID DATA FROM ROOT
call mpi_recv(x%x1,lx1+4,mpi_realprec,0,tagx1,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv x1'
call mpi_recv(x%x2all,lx2all+4,mpi_realprec,0,tagx2all,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv x2all'
call mpi_recv(x%x3all,lx3all+4,mpi_realprec,0,tagx3all,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv x3all'
call mpi_recv(x%dx1,lx1+3,mpi_realprec,0,tagx1,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv dx1'
call mpi_recv(x%x1i,lx1+1,mpi_realprec,0,tagx1,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv x1i'
call mpi_recv(x%dx1i,lx1,mpi_realprec,0,tagx1,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (ierr/=0) error stop 'failed mpi_recv dx1i'


call bcast_recv1D_3(x%x3,tagx3)
x%dx3=x%x3(0:lx3+2)-x%x3(-1:lx3+1)     !computing these avoids extra message passing (could be done for other coordinates
x%x3i(1:lx3+1)=0.5*(x%x3(0:lx3)+x%x3(1:lx3+1))
x%dx3i=x%x3i(2:lx3+1)-x%x3i(1:lx3)


call bcast_recv1D_2(x%x2,tagx2)
x%dx2=x%x2(0:lx2+2)-x%x2(-1:lx2+1)
x%x2i(1:lx2+1)=0.5*(x%x2(0:lx2)+x%x2(1:lx2+1))
x%dx2i=x%x2i(2:lx2+1)-x%x2i(1:lx2)


block
real(wp), dimension(-1:lx1+2,-1:lx2+2,-1:lx3+2) :: mpirecvbuf

call bcast_recv3D_ghost(mpirecvbuf,tagh1)
x%h1=mpirecvbuf
call bcast_recv3D_ghost(mpirecvbuf,tagh2)
x%h2=mpirecvbuf
call bcast_recv3D_ghost(mpirecvbuf,tagh3)
x%h3=mpirecvbuf
end block

call bcast_recv(x%h1x1i,tagh1)
call bcast_recv(x%h2x1i,tagh2)
call bcast_recv(x%h3x1i,tagh3)

call bcast_recv3D_x2i(x%h1x2i,tagh1)
call bcast_recv3D_x2i(x%h2x2i,tagh2)
call bcast_recv3D_x2i(x%h3x2i,tagh3)

call bcast_recv3D_x3i(x%h1x3i,tagh1)
call bcast_recv3D_x3i(x%h2x3i,tagh2)
call bcast_recv3D_x3i(x%h3x3i,tagh3)

call bcast_recv(g1,tagh1)
call bcast_recv(g2,tagh2)
call bcast_recv(g3,tagh3)

call bcast_recv(x%glat,tagglat)
call bcast_recv(x%glon,tagglon)
call bcast_recv(x%alt,tagalt)

call bcast_recv(x%Bmag,tagBmag)
call bcast_recv(x%I,taginc)           !only time that we need to exchange 2D array data, I think
call bcast_recv(x%nullpts,tagnull)    !note that this is not an integer array

block
real(wp), dimension(1:lx1,1:lx2,1:lx3) :: mpirecvbuf
do icomp=1,3
  call bcast_recv(mpirecvbuf,tageunit1)
  x%e1(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  call bcast_recv(mpirecvbuf,tageunit2)
  x%e2(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  call bcast_recv(mpirecvbuf,tageunit3)
  x%e3(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  call bcast_recv(mpirecvbuf,tager)
  x%er(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  call bcast_recv(mpirecvbuf,tagetheta)
  x%etheta(:,:,:,icomp)=mpirecvbuf
end do
do icomp=1,3
  call bcast_recv(mpirecvbuf,tagephi)
  x%ephi(:,:,:,icomp)=mpirecvbuf
end do
end block

call bcast_recv(x%r,tagr)
call bcast_recv(x%theta,tagtheta)
call bcast_recv(x%phi,tagphi)


!COUNT THE NUMBER OF NULL GRID POINTS AND GENERATE A LIST OF NULL INDICES FOR LATER USE
x%lnull=0;
do ix3=1,lx3
  do ix2=1,lx2
    do ix1=1,lx1
      if(x%nullpts(ix1,ix2,ix3)>0.5d0) then
        x%lnull=x%lnull+1
      end if
    end do
  end do
end do
allocate(x%inull(1:x%lnull,1:3))

icount=1
do ix3=1,lx3
  do ix2=1,lx2
    do ix1=1,lx1
      if(x%nullpts(ix1,ix2,ix3)>0.5d0) then
        x%inull(icount,:)=[ix1,ix2,ix3]
        icount=icount+1
      end if
    end do
  end do
end do
print *, 'Done computing null grid points...  Process:  ',myid,' has:  ',x%lnull


!! COMPUTE DIFFERENTIAL DISTANCES ALONG EACH DIRECTION
block
real(wp), dimension(1:lx1,1:lx2,1:lx3) :: tmpdx
allocate(x%dl1i(1:lx1,1:lx2,1:lx3),x%dl2i(1:lx1,1:lx2,1:lx3),x%dl3i(1:lx1,1:lx2,1:lx3))

tmpdx=spread(spread(x%dx1i,2,lx2),3,lx3)
x%dl1i=tmpdx*x%h1(1:lx1,1:lx2,1:lx3)
tmpdx=spread(spread(x%dx2i,1,lx1),3,lx3)
x%dl2i=tmpdx*x%h2(1:lx1,1:lx2,1:lx3)
tmpdx=spread(spread(x%dx3i,1,lx1),2,lx2)
x%dl3i=tmpdx*x%h3(1:lx1,1:lx2,1:lx3)
end block

end subroutine read_grid_workers

end submodule grid_read