module grid

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

use mesh, only: curvmesh

use phys_consts, only: Gconst,Me,Re,wp,red,black

use reader, only: get_simsize3

use mpimod, only: mpi_integer, mpi_comm_world, mpi_status_ignore, &

myid, lid, lid2, lid3, &

tagx1, tagx2, tagx3, tagtheta, tagr, tagphi, tagnull, taglx1, taglx2, taglx3, taglx3all, taginc, &

tagh1, tagh2, tagh3, tagglat, tagglon, tageunit1, tageunit2, tageunit3, tagetheta, tager, &

tagalt, tagbmag, tagephi, tagswap, &

mpi_realprec, taglx2all, tagx3all, tagx2all, &

bcast_recv, bcast_send, bcast_recv3D_ghost, bcast_send3D_ghost, bcast_recv3D_x3i, bcast_send3D_x3i, &

bcast_send3D_x2i,bcast_recv3D_x2i, bcast_send1D_2, bcast_recv1D_2, bcast_send1D_3, bcast_recv1D_3

implicit none

private

integer, protected :: lx1,lx2,lx3,lx2all,lx3all

!! this is a useful shorthand for most program units using this module,

!! occassionally a program unit needs to define its own size in which case an only statement

!! is required when using this module.

real(wp), dimension(:,:,:), allocatable, protected :: g1,g2,g3 !gravity, not to be modified by a procedure outside this module

integer, protected :: gridflag !for cataloguing the type of grid that we are using, open, closed, inverted, etc.

integer :: flagswap !have the x2 and x3 dimensions been swapped?

interface ! read.f90

module subroutine read_grid(indatsize,indatgrid,flagperiodic,x)

character(*), intent(in) :: indatsize,indatgrid

integer, intent(in) :: flagperiodic

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

end subroutine read_grid

end interface

public :: lx1,lx2,lx3, lx2all,lx3all, gridflag, flagswap, clear_unitvecs, g1,g2,g3, &

read_grid, clear_grid, grid_size

contains

subroutine grid_size(indatsize)

!! CHECK THE SIZE OF THE GRID TO BE LOADED AND SET SIZES IN THIS MODULE (NOT IN STRUCTURE THOUGH)

character(*), intent(in) :: indatsize

integer :: iid, ierr

logical exists

if (myid==0) then !root must physically read the size info and pass to workers

!! DETERMINE THE SIZE OF THE GRID TO BE LOADED

call get_simsize3(indatsize, lx1, lx2all, lx3all)

if (lx1 < 1 .or. lx2all < 1 .or. lx3all < 1) then

write(stderr,*) 'ERROR: reading ' // indatsize

error stop 'grid.f90: grid size must be strictly positive'

endif

!! check correct number of MPI images. Help avoid confusing errors or bad simulations

if (lx2all > 1) then

if (modulo(lx2all, lid2) /= 0) then

write(stderr,'(/,A,I6,A,I6,/)') 'ERROR: Number of MPI images along x2', lid2, ' is not an integer factor of lx2all: ', lx2all

error stop

endif

endif

if (lx3all > 1) then

if (modulo(lx3all, lid3) /= 0) then

write(stderr,'(/,A,I6,A,I6,/)') 'ERROR: Number of MPI images along x3', lid3, ' is not an integer factor of lx3all: ', lx3all

error stop

endif

endif

do iid=1,lid-1

call mpi_send(lx1,1,MPI_INTEGER,iid,taglx1,MPI_COMM_WORLD,ierr)

if (ierr/=0) error stop 'grid:grid_size lx1 failed mpi_send'

call mpi_send(lx2all,1,MPI_INTEGER,iid,taglx2all,MPI_COMM_WORLD,ierr)

if (ierr/=0) error stop 'grid:grid_size lx2all failed mpi_send'

call mpi_send(lx3all,1,MPI_INTEGER,iid,taglx3all,MPI_COMM_WORLD,ierr)

if (ierr/=0) error stop 'grid:grid_size lx3all failed mpi_send'

end do

print *, 'grid:grid_size reporting full grid size: ',lx1,lx2all,lx3all

else

call mpi_recv(lx1,1,MPI_INTEGER,0,taglx1,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

if (ierr/=0) error stop 'grid:grid_size lx1 failed mpi_send'

call mpi_recv(lx2all,1,MPI_INTEGER,0,taglx2all,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

if (ierr/=0) error stop 'grid:grid_size lx2all failed mpi_send'

call mpi_recv(lx3all,1,MPI_INTEGER,0,taglx3all,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

if (ierr/=0) error stop 'grid:grid_size lx3all failed mpi_send'

end if

end subroutine grid_size

subroutine clear_grid(x)

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

!------------------------------------------------------------

!-------DEALLOCATES GRID VARIABLES.

!------------------------------------------------------------

deallocate(x%x3all,x%x2all)

deallocate(x%x1,x%x2,x%x3)

deallocate(x%dx1i,x%x1i,x%dx1)

deallocate(x%dx2i,x%x2i,x%dx2)

deallocate(x%dx3i,x%x3i,x%dx3)

deallocate(x%glat,x%glon,x%alt)

deallocate(x%r,x%theta,x%phi)

deallocate(x%h1,x%h2,x%h3)

deallocate(x%h1x1i,x%h2x1i,x%h3x1i)

deallocate(x%h1x2i,x%h2x2i,x%h3x2i)

deallocate(x%h1x3i,x%h2x3i,x%h3x3i)

deallocate(x%I,x%Bmag,x%nullpts)

! deallocate(x%e1,x%e2,x%e3) !handled by clear_unitvecs (assuming netural perturbations are used)

! deallocate(x%er,x%etheta,x%ephi)

deallocate(x%dl1i,x%dl2i,x%dl3i)

if (myid == 0) then

deallocate(x%x2iall,x%dx2all,x%dx2iall)

deallocate(x%x3iall,x%dx3all,x%dx3iall)

deallocate(x%h1all,x%h2all,x%h3all)

deallocate(x%h1x1iall,x%h2x1iall,x%h3x1iall)

deallocate(x%h1x2iall,x%h2x2iall,x%h3x2iall)

deallocate(x%h1x3iall,x%h2x3iall,x%h3x3iall)

deallocate(x%rall,x%thetaall,x%phiall)

end if

!THIS NEEDS TO DEALLOCATE BOTH GRAVITY AND THE MAGNETIC FIELD MAGNITUDE

call clear_grav()

end subroutine clear_grid

subroutine clear_unitvecs(x)

!! DEALLOCATE GRID UNIT VECTORS, WHICH TAKE UP A LOT OF MEMORY

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

deallocate(x%e1,x%e2,x%e3,x%er,x%etheta,x%ephi)

end subroutine clear_unitvecs

subroutine load_grav(alt)

!! LOAD UP GRAV. FIELD ARRAY. IT IS EXPECTED THAT

!! GHOST CELLS WILL HAVE BEEN TRIMMED FROM ARRAYS BEFORE THEY ARE PASSED INTO THIS ROUTINE.

real(wp), dimension(:,:,:), intent(in) :: alt

allocate(g1(lx1,lx2,lx3),g2(lx1,lx2,lx3),g3(lx1,lx2,lx3))

g1 = -1 * Gconst * Me / (Re + alt)**2

end subroutine load_grav

subroutine clear_grav()

!! DEALLOCATE GRAV. FIELD ARRAY.

deallocate(g1,g2,g3)

end subroutine clear_grav

end module grid