module neutral

use, intrinsic :: iso_fortran_env, only: sp => real32

use, intrinsic :: ieee_arithmetic, only: ieee_is_finite

use mpi, only: mpi_integer, mpi_comm_world, mpi_status_ignore

use phys_consts, only: wp, lnchem, pi, re, debug

use grid, only: curvmesh, lx1, lx2, lx3, clear_unitvecs, gridflag

use interpolation, only : interp2, interp3

use timeutils, only : doy_calc,dateinc, date_filename

use mpimod, only: myid, lid, taglrho, taglz, mpi_realprec, tagdno, tagdnn2, tagdno2, tagdtn, tagdvnrho, tagdvnz, tagly, &

taglx,tagxn,tagxnrange,tagynrange,tagyn,tagzn,tagdvnx

! also links gtd7 from vendor/msis00/

implicit none

private

!! ALL ARRAYS THAT FOLLOW ARE USED WHEN INCLUDING NEUTRAL PERTURBATIONS FROM ANOTHER MODEL

!! ARRAYS TO STORE THE NEUTRAL GRID INFORMATION

!! as long as the neutral module is in scope these persist and do not require a "save"; this variable only used by the axisymmetric interpolation

real(wp), dimension(:), allocatable, private :: rhon !used for axisymmetric 2D simulations

real(wp), dimension(:), allocatable, private :: yn !used in cartesian 2D and 3D interpolation

real(wp), dimension(:), allocatable, private :: zn

real(wp), dimension(:), allocatable, private :: xn !for 3D cartesian interpolation

integer, private :: lrhon,lzn,lyn,lxn

!! STORAGE FOR NEUTRAL SIMULATION DATA.

! These will be singleton in the second dimension (longitude) in the case of 2D interpolation...

!! THESE ARE INCLUDED AS MODULE VARIATIONS TO AVOID HAVING TO REALLOCATE AND DEALLOCIATE EACH TIME WE NEED TO INTERP

real(wp), dimension(:,:,:), allocatable, private :: dnO,dnN2,dnO2,dvnrho,dvnz,dvnx,dTn

!!full grid parameters for root to store input from files.

real(wp), dimension(:), allocatable, private :: xnall

real(wp), dimension(:), allocatable, private :: ynall

integer, private :: lxnall,lynall

real(wp), dimension(:,:,:), allocatable, private :: dnOall,dnN2all,dnO2all,dvnrhoall,dvnzall,dvnxall,dTnall

!ARRAYS TO STORE NEUTRAL DATA THAT HAS BEEN INTERPOLATED

real(wp), dimension(:,:,:), allocatable, private :: dnOiprev,dnN2iprev,dnO2iprev,dvnrhoiprev,dvnziprev,dTniprev, &

dvn1iprev,dvn2iprev,dvn3iprev,dvnxiprev

real(wp), private :: tprev

integer, dimension(3), private :: ymdprev

!! denoted time corresponding to "prev" interpolated data

real(wp), private :: UTsecprev

real(wp), dimension(:,:,:), allocatable, private :: dnOinext,dnN2inext,dnO2inext,dvnrhoinext,dvnzinext, &

dTninext,dvn1inext,dvn2inext,dvn3inext,dvnxinext

real(wp), private :: tnext

integer, dimension(3), private :: ymdnext

real(wp), private :: UTsecnext

!SPACE TO STORE PROJECTION FACTORS

real(wp), dimension(:,:,:), allocatable, private :: proj_erhop_e1,proj_ezp_e1,proj_erhop_e2,proj_ezp_e2,proj_erhop_e3,proj_ezp_e3 !these projections are used in the axisymmetric interpolation

real(wp), dimension(:,:,:), allocatable, private :: proj_eyp_e1,proj_eyp_e2,proj_eyp_e3 !these are for Cartesian projections

real(wp), dimension(:,:,:), allocatable, private :: proj_exp_e1,proj_exp_e2,proj_exp_e3

!PLASMA GRID ZI AND RHOI LOCATIONS FOR INTERPOLATIONS

real(wp), dimension(:), allocatable, private :: zi,yi,xi,rhoi !this is to be a flat listing of sites on the, rhoi only used in axisymmetric and yi only in cartesian

!USED FOR 3D INTERPOLATION WHERE WORKER DIVISIONS ARE COMPLICATED (note that the first dim starts at zero so it matches mpi ID)

real(wp), dimension(:,:), private, allocatable :: extents !roots array that is used to store min/max x,y,z of each works

integer, dimension(:,:), private, allocatable :: indx !roots array that contain indices for each workers needed piece of the neutral data

integer, dimension(:,:), private, allocatable :: slabsizes

!! BASE MSIS ATMOSPHERIC STATE ON WHICH TO APPLY PERTURBATIONS

real(wp), dimension(:,:,:,:), allocatable, protected :: nnmsis

real(wp), dimension(:,:,:), allocatable, protected :: Tnmsis

real(wp), dimension(:,:,:), allocatable, protected :: vn1base,vn2base,vn3base

public :: Tnmsis, neutral_atmos, make_dneu, clear_dneu, neutral_perturb

contains

subroutine neutral_atmos(ymd,UTsecd,glat,glon,alt,activ,nn,Tn)

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

!-------CALL NRL-MSISE-00 AND ORGANIZE THE RESULTS. APPEND

!-------OTHER AUXILIARY NEUTRAL DENSITY DATA USED BY MAIN

!-------CODE

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

integer, dimension(3) :: ymd

real(wp) :: UTsecd

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

real(wp), dimension(3) :: activ

real(wp), dimension(1:size(alt,1),1:size(alt,2),1:size(alt,3),lnchem), intent(out) :: nn

real(wp), dimension(1:size(alt,1),1:size(alt,2),1:size(alt,3)), intent(out) :: Tn

integer :: ix1,ix2,ix3,lx1,lx2,lx3

integer :: iyd,mass=48

integer :: dom,month,year,doy,yearshort

real :: sec,f107a,f107,ap(7),stl,ap3

real :: altnow,latnow,lonnow

real :: d(9),t(2)

! real(wp), dimension(1:size(alt,1),1:size(alt,2),1:size(alt,3)) :: nnow

! real(wp), dimension(1:size(alt,1),1:size(alt,2),1:size(alt,3)) :: altalt !an alternate altitude variable which fixes below ground values to 1km

lx1=size(alt,1)

lx2=size(alt,2)

lx3=size(alt,3)

!! CONVERT DATE INFO INTO EXPECTED FORM AND KIND

f107a=real(activ(1),sp)

f107=real(activ(2),sp)

ap=real(activ(3),sp)

ap3=real(activ(3),sp)

dom=ymd(3)

month=ymd(2)

year=ymd(1)

doy=doy_calc(year, month, dom)

yearshort=mod(year,100)

iyd=yearshort*1000+doy

sec=floor(UTsecd)

ap(2)=ap3 !superfluous for now

!! ITERATED LAT, LON, ALT DATA

call meters(.true.) !switch to mksa units

do ix3=1,lx3

do ix2=1,lx2

do ix1=1,lx1

altnow=real(alt(ix1,ix2,ix3)/1d3,sp)

if (altnow<0.0) then

altnow = 1._sp !so that MSIS does not get called with below ground values and so that we set them to something sensible that won't mess up the conductance calculations

end if

latnow=real(glat(ix1,ix2,ix3),sp)

lonnow=real(glon(ix1,ix2,ix3),sp)

stl=sec/3600.0+lonnow/15.0

call gtd7(iyd,sec,altnow,latnow,lonnow,stl,f107a,f107,ap,mass,d,t)

nn(ix1,ix2,ix3,1)=real(d(2),wp)

nn(ix1,ix2,ix3,2)=real(d(3),wp)

nn(ix1,ix2,ix3,3)=real(d(4),wp)

nn(ix1,ix2,ix3,4)=real(d(7),wp)

nn(ix1,ix2,ix3,5)=real(d(8),wp)

Tn(ix1,ix2,ix3)=real(t(2),wp)

nn(ix1,ix2,ix3,6)=4d-1*exp(-3700.0/Tn(ix1,ix2,ix3))*nn(ix1,ix2,ix3,3)+ &

5d-7*nn(ix1,ix2,ix3,1) !Mitra, 1968

end do

end do

end do

!UPDATE THE REFERENCE ATMOSPHERE VALUES

nnmsis=nn; Tnmsis=Tn; vn1base=0d0; vn2base=0d0; vn3base=0d0;

end subroutine neutral_atmos

!THIS IS WRAPPER FOR THE NEUTRAL PERTURBATION CODES THAT DO EITHER

!AXISYMMETRIC OR CARTESIAN OR 3D INTERPOLATION

subroutine neutral_perturb(interptype,dt,dtneu,t,ymd,UTsec,neudir,dxn,drhon,dzn,meanlat,meanlong,x,nn,Tn,vn1,vn2,vn3)

integer, intent(in) :: interptype

real(wp), intent(in) :: dt,dtneu

real(wp), intent(in) :: t

integer, dimension(3), intent(in) :: ymd !date for which we wish to calculate perturbations

real(wp), intent(in) :: UTsec

real(wp), intent(in) :: dxn,drhon,dzn !neutral grid spacing

real(wp), intent(in) :: meanlat, meanlong !neutral source center location

character(*), intent(in) :: neudir !directory where neutral simulation data is kept

type(curvmesh), intent(inout) :: x !grid structure (inout becuase we want to be able to deallocate unit vectors once we are done with them)

real(wp), dimension(:,:,:,:), intent(out) :: nn !neutral params interpolated to plasma grid at requested time

real(wp), dimension(:,:,:), intent(out) :: Tn,vn1,vn2,vn3

if (interptype==0) then !cartesian interpolation drho inputs (radial distance) will be interpreted as dy (horizontal distance)

call neutral_perturb_cart(dt,dtneu,t,ymd,UTsec,neudir,drhon,dzn,meanlat,meanlong,x,nn,Tn,vn1,vn2,vn3)

else if (interptype==1) then !axisymmetric interpolation

call neutral_perturb_axisymm(dt,dtneu,t,ymd,UTsec,neudir,drhon,dzn,meanlat,meanlong,x,nn,Tn,vn1,vn2,vn3)

else if (interptype==3) then !3D interpolation drhon is takent to be dyn (northward distance)

call neutral_perturb_3D(dt,dtneu,t,ymd,UTsec,neudir,dxn,drhon,dzn,meanlat,meanlong,x,nn,Tn,vn1,vn2,vn3)

else

error stop '...Invalid interpolation type specified from input file...'

end if

end subroutine neutral_perturb

!FOR CONSISTENCY I'D LIKE TO STRUCTURE NEUTRAL PERTURB OPERATIONS LIKE GRAVITY IS HANDLED IN THE GRID MODULE, I.E. HAVE AN EXPLICIT CONSTRUCTORS/DESTRUCTOR TYPE ROUTINE THAT HANDLES ALLOCATION AND DEALLOCATION, WHICH WILL CLEAN UP THE NEUTRAL_PERTURB SUBROUTINE, I.E. REMOVE ALLOCATES OF PERSISTENT MODULE VARIABLES.

subroutine neutral_perturb_axisymm(dt,dtneu,t,ymd,UTsec,neudir,drhon,dzn,meanlat,meanlong,x,nn,Tn,vn1,vn2,vn3)

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

!-------COMPUTE NEUTRAL PERTURBATIONS FOR THIS TIME STEP. ADD

!-------THEM TO MSIS PERTURBATIONS TO GET ABSOLUTE VALUES FOR

!-------EACH PARAMETER

!-------

!-------THIS VERSION ASSUMES THE INPUT NEUTRAL DATA ARE IN

!-------CYLINDRICAL COORDINATES.

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

real(wp), intent(in) :: dt,dtneu

real(wp), intent(in) :: t

integer, dimension(3), intent(in) :: ymd !date for which we wish to calculate perturbations

real(wp), intent(in) :: UTsec

real(wp), intent(in) :: drhon,dzn !neutral grid spacing

real(wp), intent(in) :: meanlat, meanlong !neutral source center location

character(*), intent(in) :: neudir !directory where neutral simulation data is kept

type(curvmesh), intent(inout) :: x !grid structure (inout becuase we want to be able to deallocate unit vectors once we are done with them)

real(wp), dimension(:,:,:,:), intent(out) :: nn !neutral params interpolated to plasma grid at requested time

real(wp), dimension(:,:,:), intent(out) :: Tn,vn1,vn2,vn3

integer :: ix1,ix2,ix3,iid!,irhon,izn

integer, dimension(3) :: ymdtmp

real(wp) :: UTsectmp

real(wp), dimension(size(nn,1),size(nn,2),size(nn,3)) :: dnOinow,dnN2inow,dnO2inow,dTninow,dvn1inow,dvn2inow,dvn3inow !current time step perturbations (centered in time)

!CHECK WHETHER WE NEED TO LOAD A NEW FILE

if (t+dt/2d0>=tnext .or. t<=0d0) then !negative time means that we need to load the first frame

!IF FIRST LOAD ATTEMPT CREATE A NEUTRAL GRID AND COMPUTE GRID SITES FOR IONOSPHERIC GRID. Since this needs an input file, I'm leaving it under this condition here

if (.not. allocated(zn)) then !means this is the first tiem we've tried to load neutral simulation data, should we check for a previous neutral file to load??? or just assume everything starts at zero? This needs to somehow check for an existing file under certain conditiosn, maybe if it==1??? Actually we don't even need that we can just check that the neutral grid is allocated (or not)

!initialize dates

ymdprev=ymd

UTsecprev=UTsec

ymdnext=ymdprev

UTsecnext=UTsecprev

!Create a neutral grid, do some allocations and projections

call gridproj_dneu2D(drhon,dzn,meanlat,meanlong,neudir,.false.,x) !set false to denote not Cartesian...

end if

!Read in neutral data from a file

call read_dneu2D(tprev,tnext,t,dtneu,dt,neudir,ymdtmp,UTsectmp,.false.)

!Spatial interpolatin for the frame we just read in

if (myid==0 .and. debug) then

print *, 'Spatial interpolation and rotation of vectors for date: ',ymdtmp,' ',UTsectmp

end if

call spaceinterp_dneu2D(.false.)

!UPDATE OUR CONCEPT OF PREVIOUS AND NEXT TIMES

tprev=tnext

UTsecprev=UTsecnext

ymdprev=ymdnext

tnext=tprev+dtneu

UTsecnext=UTsectmp

ymdnext=ymdtmp

end if !done loading frame data...

!Interpolation in time

call timeinterp_dneu(t,dt,dNOinow,dnN2inow,dnO2inow,dvn1inow,dvn2inow,dvn3inow,dTninow)

!Add interpolated perturbations to reference atmosphere arrays

nn(:,:,:,1)=nnmsis(:,:,:,1)+dnOinow

nn(:,:,:,2)=nnmsis(:,:,:,2)+dnN2inow

nn(:,:,:,3)=nnmsis(:,:,:,3)+dnO2inow

nn(:,:,:,1)=max(nn(:,:,:,1),1._wp)

nn(:,:,:,2)=max(nn(:,:,:,2),1._wp)

nn(:,:,:,3)=max(nn(:,:,:,3),1._wp)

Tn=Tnmsis+dTninow

Tn=max(Tn,50._wp)

vn1=vn1base+dvn1inow

vn2=vn2base+dvn2inow

vn3=vn3base+dvn3inow

end subroutine neutral_perturb_axisymm

!! THIS SHARES SO MUCH CODE WITH THE AXISYMMETRIC VERSION THAT THEY SHOULD PROBABLY BE COMBINED

subroutine neutral_perturb_cart(dt,dtneu,t,ymd,UTsec,neudir,dyn,dzn,meanlat,meanlong,x,nn,Tn,vn1,vn2,vn3)

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

!-------COMPUTE NEUTRAL PERTURBATIONS FOR THIS TIME STEP. ADD

!-------THEM TO MSIS PERTURBATIONS TO GET ABSOLUTE VALUES FOR

!-------EACH PARAMETER.

!-------

!-------THIS VERSION ASSUMES THE INPUT NEUTRAL DATA ARE IN

!-------CARTESIAN COORDINATES.

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

real(wp), intent(in) :: dt,dtneu

real(wp), intent(in) :: t

integer, dimension(3), intent(in) :: ymd !date for which we wish to calculate perturbations

real(wp), intent(in) :: UTsec

real(wp), intent(in) :: dyn,dzn !neutral grid spacing

real(wp), intent(in) :: meanlat, meanlong !neutral source center location

character(*), intent(in) :: neudir !directory where neutral simulation data is kept

type(curvmesh), intent(inout) :: x !grid structure (inout becuase we want to be able to deallocate unit vectors once we are done with them)

real(wp), dimension(:,:,:,:), intent(out) :: nn !neutral params interpolated to plasma grid at requested time

real(wp), dimension(:,:,:), intent(out) :: Tn,vn1,vn2,vn3

integer :: ix1,ix2,ix3,iid

integer, dimension(3) :: ymdtmp

real(wp) :: UTsectmp

real(wp), dimension(size(nn,1),size(nn,2),size(nn,3)) :: dnOinow,dnN2inow,dnO2inow,dTninow,dvn1inow,dvn2inow,dvn3inow !current time step perturbations (centered in time)

!CHECK WHETHER WE NEED TO LOAD A NEW FILE

if (t+dt/2d0>=tnext .or. t<=0d0) then

!IF FIRST LOAD ATTEMPT CREATE A NEUTRAL GRID AND COMPUTE GRID SITES FOR IONOSPHERIC GRID. Since this needs an input file, I'm leaving it under this condition here

if (.not. allocated(zn)) then !means this is the first tiem we've tried to load neutral simulation data, should we check for a previous neutral file to load???

!initialize dates

ymdprev=ymd

UTsecprev=UTsec

ymdnext=ymdprev

UTsecnext=UTsecprev

!Create a neutral grid, do some allocations and projections

call gridproj_dneu2D(dyn,dzn,meanlat,meanlong,neudir,.true.,x) !set true to denote Cartesian...

end if

!Read in neutral data from a file

call read_dneu2D(tprev,tnext,t,dtneu,dt,neudir,ymdtmp,UTsectmp,.true.)

!Spatial interpolatin for the frame we just read in

if (myid==0 .and. debug) then

print *, 'Spatial interpolation and rotation of vectors for date: ',ymdtmp,' ',UTsectmp

end if

call spaceinterp_dneu2D(.true.)

!UPDATE OUR CONCEPT OF PREVIOUS AND NEXT TIMES

tprev=tnext

UTsecprev=UTsecnext

ymdprev=ymdnext

tnext=tprev+dtneu

UTsecnext=UTsectmp

ymdnext=ymdtmp

end if

!Interpolation in time

call timeinterp_dneu(t,dt,dNOinow,dnN2inow,dnO2inow,dvn1inow,dvn2inow,dvn3inow,dTninow)

!NOW UPDATE THE PROVIDED NEUTRAL ARRAYS

nn(:,:,:,1)=nnmsis(:,:,:,1)+dnOinow

nn(:,:,:,2)=nnmsis(:,:,:,2)+dnN2inow

nn(:,:,:,3)=nnmsis(:,:,:,3)+dnO2inow

nn(:,:,:,1)=max(nn(:,:,:,1),1._wp)

nn(:,:,:,2)=max(nn(:,:,:,2),1._wp)

nn(:,:,:,3)=max(nn(:,:,:,3),1._wp)

Tn=Tnmsis+dTninow

Tn=max(Tn,51._wp)

vn1=vn1base+dvn1inow

vn2=vn2base+dvn2inow

vn3=vn3base+dvn3inow

end subroutine neutral_perturb_cart

subroutine neutral_perturb_3D(dt,dtneu,t,ymd,UTsec,neudir,dxn,dyn,dzn,meanlat,meanlong,x,nn,Tn,vn1,vn2,vn3)

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

!-------COMPUTE NEUTRAL PERTURBATIONS FOR THIS TIME STEP. ADD

!-------THEM TO MSIS PERTURBATIONS TO GET ABSOLUTE VALUES FOR

!-------EACH PARAMETER.

!-------

!-------THIS VERSION ASSUMES THE INPUT NEUTRAL DATA ARE IN

!-------CARTESIAN COORDINATES AND IN THREE DIMENSIONS.

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

real(wp), intent(in) :: dt,dtneu

real(wp), intent(in) :: t

integer, dimension(3), intent(in) :: ymd !date for which we wish to calculate perturbations

real(wp), intent(in) :: UTsec

real(wp), intent(in) :: dxn,dyn,dzn !neutral grid spacing

real(wp), intent(in) :: meanlat, meanlong !neutral source center location

character(*), intent(in) :: neudir !directory where neutral simulation data is kept

type(curvmesh), intent(inout) :: x !grid structure (inout becuase we want to be able to deallocate unit vectors once we are done with them)

real(wp), dimension(:,:,:,:), intent(out) :: nn !neutral params interpolated to plasma grid at requested time

real(wp), dimension(:,:,:), intent(out) :: Tn,vn1,vn2,vn3

integer :: ix1,ix2,ix3,iid

integer, dimension(3) :: ymdtmp

real(wp) :: UTsectmp

real(wp), dimension(size(nn,1),size(nn,2),size(nn,3)) :: dnOinow,dnN2inow,dnO2inow,dTninow,dvn1inow,dvn2inow,dvn3inow !current time step perturbations (centered in time)

real(wp) :: starttime,endtime

!CHECK WHETHER WE NEED TO LOAD A NEW FILE

if (t+dt/2d0>=tnext .or. t<=0d0) then

!IF FIRST LOAD ATTEMPT CREATE A NEUTRAL GRID AND COMPUTE GRID SITES FOR IONOSPHERIC GRID. Since this needs an input file, I'm leaving it under this condition here

if (.not. allocated(zn)) then !means this is the first tiem we've tried to load neutral simulation data, should we check for a previous neutral file to load???

!initialize dates

ymdprev=ymd

UTsecprev=UTsec

ymdnext=ymdprev

UTsecnext=UTsecprev

!Create a neutral grid, do some allocations and projections

if (myid==0 .and. debug) then

print*, 'Creating a neutral grid...'

end if

call gridproj_dneu3D(dxn,dyn,dzn,meanlat,meanlong,neudir,x) !set true to denote Cartesian...

end if

!Read in neutral data from a file

if (myid==0 .and. debug) then

print*, 'Reading in data from neutral file'

call cpu_time(starttime)

end if

call read_dneu3D(tprev,tnext,t,dtneu,dt,neudir,ymdtmp,UTsectmp)

if (myid==0 .and. debug) then

call cpu_time(endtime)

print*, 'Neutral data input required time: ',endtime-starttime

end if

!Spatial interpolatin for the frame we just read in

if (myid==0 .and. debug) then

print *, 'Spatial interpolation and rotation of vectors for date: ',ymdtmp,' ',UTsectmp

call cpu_time(starttime)

end if

call spaceinterp_dneu3D()

if (myid==0 .and. debug) then

call cpu_time(endtime)

print*, 'Spatial interpolation in 3D took time: ',endtime-starttime

end if

!UPDATE OUR CONCEPT OF PREVIOUS AND NEXT TIMES

tprev=tnext

UTsecprev=UTsecnext

ymdprev=ymdnext

tnext=tprev+dtneu

UTsecnext=UTsectmp

ymdnext=ymdtmp

end if

!Interpolation in time

if (myid==0 .and. debug) then

print*, 'Interpolating in time'

end if

call timeinterp_dneu(t,dt,dNOinow,dnN2inow,dnO2inow,dvn1inow,dvn2inow,dvn3inow,dTninow)

!NOW UPDATE THE PROVIDED NEUTRAL ARRAYS

nn(:,:,:,1)=nnmsis(:,:,:,1)+dnOinow

nn(:,:,:,2)=nnmsis(:,:,:,2)+dnN2inow

nn(:,:,:,3)=nnmsis(:,:,:,3)+dnO2inow

nn(:,:,:,1)=max(nn(:,:,:,1),1._wp)

nn(:,:,:,2)=max(nn(:,:,:,2),1._wp)

nn(:,:,:,3)=max(nn(:,:,:,3),1._wp)

Tn=Tnmsis+dTninow

Tn=max(Tn,51._wp)

vn1=vn1base+dvn1inow

vn2=vn2base+dvn2inow

vn3=vn3base+dvn3inow

end subroutine neutral_perturb_3D

subroutine gridproj_dneu2D(dhorzn,dzn,meanlat,meanlong,neudir,flagcart,x)

!Read in the grid for the neutral data and project unit vectors into the appropriiate directions.

!Also allocate module-scope variables for storing neutral perturbations read in from input files.

real(wp), intent(in) :: dhorzn,dzn !neutral grid spacing in horizontal "rho or y" and vertical directions

real(wp), intent(in) :: meanlat, meanlong !neutral source center location

character(*), intent(in) :: neudir !< directory where neutral simulation data is kept

logical, intent(in) :: flagcart !whether or not the input data are to be interpreted as Cartesian

type(curvmesh), intent(inout) :: x !inout to allow deallocation of unit vectors once we are done with them, should consider exporting this to another functino to be called from main program to avoid having x writeable...

integer :: lhorzn

real(wp) :: meanyn

character(:), allocatable :: filename

real(wp) :: theta1,phi1,theta2,phi2,gammarads,theta3,phi3,gamma1,gamma2,phip

real(wp) :: xp,yp

real(wp), dimension(3) :: erhop,ezp,eyp,tmpvec

real(wp) :: tmpsca

integer :: ix1,ix2,ix3,ihorzn,izn,iid,ierr

real(wp), dimension(x%lx1,x%lx2,x%lx3) :: zimat,rhoimat,yimat

!Establish the size of the grid based on input file and distribute to workers

if (myid==0) then !root

filename = neudir // '/simsize.dat'

print '(A,/,A)', 'Inputting neutral size from file: ',filename

block

integer :: u

open(newunit=u, file=filename,status='old',form='unformatted',access='stream')

read(u) lhorzn,lzn

close(u)

end block

print *, 'Neutral data has lhorzn,lz size: ',lhorzn,lzn,' with spacing dhorzn,dz',dhorzn,dzn

if (lhorzn < 1 .or. lzn < 1) error stop 'grid size must be strictly positive'

do iid=1,lid-1

call mpi_send(lhorzn,1,MPI_INTEGER,iid,taglrho,MPI_COMM_WORLD,ierr)

call mpi_send(lzn,1,MPI_INTEGER,iid,taglz,MPI_COMM_WORLD,ierr)

end do

else !workers

call mpi_recv(lhorzn,1,MPI_INTEGER,0,taglrho,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

call mpi_recv(lzn,1,MPI_INTEGER,0,taglz,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

end if

!Everyone must allocate space for the grid of input data

allocate(zn(lzn)) !these are module-scope variables

if (flagcart) then

allocate(rhon(1)) !not used in Cartesian code so just set to something

allocate(yn(lhorzn))

lyn=lhorzn

else

allocate(rhon(lhorzn))

allocate(yn(1)) !not used in the axisymmetric code so just initialize to something

lrhon=lhorzn

end if

!Note that the second dimension ("longitude") is singleton so that we are able to also use these vars for 3D input

allocate(dnO(lzn,1,lhorzn),dnN2(lzn,1,lhorzn),dnO2(lzn,1,lhorzn),dvnrho(lzn,1,lhorzn),dvnz(lzn,1,lhorzn),dTn(lzn,1,lhorzn))

!Define a grid (input data) by assuming that the spacing is constant

if (flagcart) then !Cartesian neutral simulation

yn=[ ((real(ihorzn,8)-1._wp)*dhorzn, ihorzn=1,lhorzn) ]

meanyn=sum(yn,1)/size(yn,1)

yn=yn-meanyn !the neutral grid should be centered on zero for a cartesian interpolation

else

rhon=[ ((real(ihorzn,8)-1._wp)*dhorzn, ihorzn=1,lhorzn) ]

end if

zn=[ ((real(izn,8)-1._wp)*dzn, izn=1,lzn) ]

if (myid==0) then

if (flagcart) then

print *, 'Creating neutral grid with y,z extent:',minval(yn),maxval(yn),minval(zn),maxval(zn)

else

print *, 'Creating neutral grid with rho,z extent: ',minval(rhon),maxval(rhon),minval(zn),maxval(zn)

end if

end if

!Neutral source locations specified in input file, here referenced by spherical magnetic coordinates.

phi1=meanlong*pi/180d0

theta1=pi/2d0-meanlat*pi/180d0

!Convert plasma simulation grid locations to z,rho values to be used in interoplation. altitude ~ zi; lat/lon --> rhoi. Also compute unit vectors and projections

if (myid==0) then

print *, 'Computing alt,radial distance values for plasma grid and completing rotations'

end if

zimat=x%alt !vertical coordinate

do ix3=1,lx3

do ix2=1,lx2

do ix1=1,lx1

!INTERPOLATION BASED ON GEOMAGNETIC COORDINATES

theta2=x%theta(ix1,ix2,ix3) !field point zenith angle

if (lx2/=1) then

phi2=x%phi(ix1,ix2,ix3) !field point azimuth, full 3D calculation

else

phi2=phi1 !assume the longitude is the samem as the source in 2D, i.e. assume the source epicenter is in the meridian of the grid

end if

!COMPUTE DISTANCES

gammarads=cos(theta1)*cos(theta2)+sin(theta1)*sin(theta2)*cos(phi1-phi2) !this is actually cos(gamma)

if (gammarads>1._wp) then !handles weird precision issues in 2D

gammarads=1._wp

else if (gammarads<-1._wp) then

gammarads=-1._wp

end if

gammarads=acos(gammarads) !angle between source location annd field point (in radians)

rhoimat(ix1,ix2,ix3)=Re*gammarads !rho here interpreted as the arc-length defined by angle between epicenter and ``field point''

!we need a phi locationi (not spherical phi, but azimuth angle from epicenter), as well, but not for interpolation - just for doing vector rotations

theta3=theta2

phi3=phi1

gamma1=cos(theta2)*cos(theta3)+sin(theta2)*sin(theta3)*cos(phi2-phi3)

if (gamma1>1._wp) then !handles weird precision issues in 2D

gamma1=1._wp

else if (gamma1<-1._wp) then

gamma1=-1._wp

end if

gamma1=acos(gamma1)

gamma2=cos(theta1)*cos(theta3)+sin(theta1)*sin(theta3)*cos(phi1-phi3)

if (gamma2>1._wp) then !handles weird precision issues in 2D

gamma2=1._wp

else if (gamma2<-1._wp) then

gamma2=-1._wp

end if

gamma2=acos(gamma2)

xp=Re*gamma1

yp=Re*gamma2 !this will likely always be positive, since we are using center of earth as our origin, so this should be interpreted as distance as opposed to displacement

!COMPUTE COORDIANTES FROM DISTANCES

if (theta3>theta1) then !place distances in correct quadrant, here field point (theta3=theta2) is is SOUTHward of source point (theta1), whreas yp is distance northward so throw in a negative sign

yp=-1._wp*yp !do we want an abs here to be safe

end if

if (phi2<phi3) then !assume we aren't doing a global grid otherwise need to check for wrapping, here field point (phi2) less than soure point (phi3=phi1)

xp=-1._wp*xp

end if

phip=atan2(yp,xp)

if(flagcart) then

yimat(ix1,ix2,ix3)=yp

end if

!PROJECTIONS FROM NEUTURAL GRID VECTORS TO PLASMA GRID VECTORS

!projection factors for mapping from axisymmetric to dipole (go ahead and compute projections so we don't have to do it repeatedly as sim runs

ezp=x%er(ix1,ix2,ix3,:)

tmpvec=ezp*x%e2(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_ezp_e2(ix1,ix2,ix3)=tmpsca

tmpvec=ezp*x%e1(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_ezp_e1(ix1,ix2,ix3)=tmpsca

tmpvec=ezp*x%e3(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec) !should be zero, but leave it general for now

proj_ezp_e3(ix1,ix2,ix3)=tmpsca

if (flagcart) then

eyp=-1._wp*x%etheta(ix1,ix2,ix3,:)

tmpvec=eyp*x%e1(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_eyp_e1(ix1,ix2,ix3)=tmpsca

tmpvec=eyp*x%e2(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_eyp_e2(ix1,ix2,ix3)=tmpsca

tmpvec=eyp*x%e3(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_eyp_e3(ix1,ix2,ix3)=tmpsca

else

erhop=cos(phip)*x%e3(ix1,ix2,ix3,:)+(-1._wp)*sin(phip)*x%etheta(ix1,ix2,ix3,:) !unit vector for azimuth (referenced from epicenter - not geocenter!!!) in cartesian geocentric-geomagnetic coords.

tmpvec=erhop*x%e1(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_erhop_e1(ix1,ix2,ix3)=tmpsca

tmpvec=erhop*x%e2(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_erhop_e2(ix1,ix2,ix3)=tmpsca

tmpvec=erhop*x%e3(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_erhop_e3(ix1,ix2,ix3)=tmpsca

end if

end do

end do

end do

!Assign values for flat lists of grid points

zi=pack(zimat,.true.) !create a flat list of grid points to be used by interpolation ffunctions

if (flagcart) then

yi=pack(yimat,.true.)

else

rhoi=pack(rhoimat,.true.)

end if

!GRID UNIT VECTORS NO LONGER NEEDED ONCE PROJECTIONS ARE CALCULATED...

call clear_unitvecs(x)

!PRINT OUT SOME BASIC INFO ABOUT THE GRID THAT WE'VE LOADED

if (myid==0 .and. debug) then

if (flagcart) then

print *, 'Min/max yn,zn values',minval(yn),maxval(yn),minval(zn),maxval(zn)

print *, 'Min/max yi,zi values',minval(yi),maxval(yi),minval(zi),maxval(zi)

else

print *, 'Min/max rhon,zn values',minval(rhon),maxval(rhon),minval(zn),maxval(zn)

print *, 'Min/max rhoi,zi values',minval(rhoi),maxval(rhoi),minval(zi),maxval(zi)

end if

print *, 'Source lat/long: ',meanlat,meanlong

print *, 'Plasma grid lat range: ',minval(x%glat(:,:,:)),maxval(x%glat(:,:,:))

print *, 'Plasma grid lon range: ',minval(x%glon(:,:,:)),maxval(x%glon(:,:,:))

end if

end subroutine gridproj_dneu2D

subroutine gridproj_dneu3D(dxn,dyn,dzn,meanlat,meanlong,neudir,x)

!Read in the grid for the neutral data and project unit vectors into the appropriiate directions.

!Also allocate module-scope variables for storing neutral perturbations read in from input files.

real(wp), intent(in) :: dxn,dyn,dzn !neutral grid spacing in horizontal "rho or y" and vertical directions

real(wp), intent(in) :: meanlat, meanlong !neutral source center location

character(*), intent(in) :: neudir !directory where neutral simulation data is kept

type(curvmesh), intent(inout) :: x !inout to allow deallocation of unit vectors once we are done with them, should consider exporting this to another functino to be called from main program to avoid having x writeable...

real(wp) :: meanyn

real(wp) :: meanxn

character(:), allocatable :: filename

real(wp) :: theta1,phi1,theta2,phi2,gammarads,theta3,phi3,gamma1,gamma2,phip

real(wp) :: xp,yp

real(wp), dimension(3) :: erhop,ezp,eyp,tmpvec,exprm

real(wp) :: tmpsca

integer :: ix1,ix2,ix3,iyn,izn,ixn,iid,ierr

real(wp), dimension(x%lx1,x%lx2,x%lx3) :: zimat,rhoimat,yimat,ximat

real(wp) :: maxzn

real(wp), dimension(2) :: xnrange,ynrange

integer, dimension(6) :: indices

!Neutral source locations specified in input file, here referenced by spherical magnetic coordinates.

phi1=meanlong*pi/180._wp

theta1=pi/2._wp-meanlat*pi/180._wp

!Convert plasma simulation grid locations to z,rho values to be used in interoplation. altitude ~ zi; lat/lon --> rhoi. Also compute unit vectors and projections

if (myid==0) then

print *, 'Computing alt,radial distance values for plasma grid and completing rotations'

end if

zimat=x%alt !vertical coordinate

do ix3=1,lx3

do ix2=1,lx2

do ix1=1,lx1

!INTERPOLATION BASED ON GEOMAGNETIC COORDINATES

theta2=x%theta(ix1,ix2,ix3) !field point zenith angle

if (lx2/=1) then

phi2=x%phi(ix1,ix2,ix3) !field point azimuth, full 3D calculation

else

phi2=phi1 !assume the longitude is the samem as the source in 2D, i.e. assume the source epicenter is in the meridian of the grid

end if

!!COMPUTE DISTANCES - ZZZ possibly superfluous for 3D case???

!gammarads=cos(theta1)*cos(theta2)+sin(theta1)*sin(theta2)*cos(phi1-phi2) !this is actually cos(gamma)

!if (gammarads>1._wp) then !handles weird precision issues in 2D

! gammarads=1._wp

!else if (gammarads<-1._wp) then

! gammarads=-1._wp

!end if

!gammarads=acos(gammarads) !angle between source location annd field point (in radians)

!rhoimat(ix1,ix2,ix3)=Re*gammarads !rho here interpreted as the arc-length defined by angle between epicenter and ``field point''

!! ZZZ end possibly superfluous block of code...

!we need a phi locationi (not spherical phi, but azimuth angle from epicenter), as well, but not for interpolation - just for doing vector rotations

theta3=theta2

phi3=phi1

gamma1=cos(theta2)*cos(theta3)+sin(theta2)*sin(theta3)*cos(phi2-phi3)

if (gamma1>1._wp) then !handles weird precision issues in 2D

gamma1=1._wp

else if (gamma1<-1._wp) then

gamma1=-1._wp

end if

gamma1=acos(gamma1)

gamma2=cos(theta1)*cos(theta3)+sin(theta1)*sin(theta3)*cos(phi1-phi3)

if (gamma2>1._wp) then !handles weird precision issues in 2D

gamma2=1._wp

else if (gamma2<-1._wp) then

gamma2=-1._wp

end if

gamma2=acos(gamma2)

xp=Re*gamma1

yp=Re*gamma2 !this will likely always be positive, since we are using center of earth as our origin, so this should be interpreted as distance as opposed to displacement

!COMPUTE COORDIANTES FROM DISTANCES

if (theta3>theta1) then !place distances in correct quadrant, here field point (theta3=theta2) is is SOUTHward of source point (theta1), whreas yp is distance northward so throw in a negative sign

yp=-1._wp*yp !do we want an abs here to be safe

end if

if (phi2<phi3) then !assume we aren't doing a global grid otherwise need to check for wrapping, here field point (phi2) less than soure point (phi3=phi1)

xp=-1._wp*xp

end if

!phip=atan2(yp,xp)

ximat(ix1,ix2,ix3)=xp !eastward distance

yimat(ix1,ix2,ix3)=yp !northward distance

!PROJECTIONS FROM NEUTURAL GRID VECTORS TO PLASMA GRID VECTORS

!projection factors for mapping from axisymmetric to dipole (go ahead and compute projections so we don't have to do it repeatedly as sim runs

ezp=x%er(ix1,ix2,ix3,:)

tmpvec=ezp*x%e2(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_ezp_e2(ix1,ix2,ix3)=tmpsca

tmpvec=ezp*x%e1(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_ezp_e1(ix1,ix2,ix3)=tmpsca

tmpvec=ezp*x%e3(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec) !should be zero, but leave it general for now

proj_ezp_e3(ix1,ix2,ix3)=tmpsca

eyp=-1._wp*x%etheta(ix1,ix2,ix3,:)

tmpvec=eyp*x%e1(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_eyp_e1(ix1,ix2,ix3)=tmpsca

tmpvec=eyp*x%e2(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_eyp_e2(ix1,ix2,ix3)=tmpsca

tmpvec=eyp*x%e3(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_eyp_e3(ix1,ix2,ix3)=tmpsca

exprm=x%ephi(ix1,ix2,ix3,:) !for 3D interpolation need to have a unit vector/projection onto x-direction (longitude)

tmpvec=exprm*x%e1(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_exp_e1(ix1,ix2,ix3)=tmpsca

tmpvec=exprm*x%e2(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_exp_e2(ix1,ix2,ix3)=tmpsca

tmpvec=exprm*x%e3(ix1,ix2,ix3,:)

tmpsca=sum(tmpvec)

proj_exp_e3(ix1,ix2,ix3)=tmpsca

end do

end do

end do

!Assign values for flat lists of grid points

zi=pack(zimat,.true.) !create a flat list of grid points to be used by interpolation functions

yi=pack(yimat,.true.)

xi=pack(ximat,.true.)

!GRID UNIT VECTORS NO LONGER NEEDED ONCE PROJECTIONS ARE CALCULATED, so go ahead and free some space

if (myid==0) then

print*, '...Clearing out unit vectors (after projections)...'

end if

call clear_unitvecs(x)

if(myid==0) then

print*, 'Projection checking: ',minval(proj_exp_e1),maxval(proj_exp_e1),minval(proj_exp_e2),maxval(proj_exp_e2), &

minval(proj_exp_e3),maxval(proj_exp_e3)

end if

!Establish the size of the grid based on input file and distribute to workers

if (myid==0) then !root

filename = neudir // '/simsize.dat'

print '(A,/,A)', 'Inputting neutral size from file:', filename

block

integer :: u

open(newunit=u,file=filename,status='old',form='unformatted',access='stream')

read(u) lxnall,lynall,lzn

close(u)

end block

print *, 'Neutral data has lx,ly,lz size: ',lxnall,lynall,lzn,' with spacing dx,dy,dz',dxn,dyn,dzn

if (lxnall < 1 .or. lynall < 1 .or. lzn < 1) error stop 'grid size must be strictly positive'

!root must allocate space for the entire grid of input data - this might be doable one parameter at a time???

allocate(zn(lzn)) !the z coordinate is never split up in message passing - want to use full altitude range...

allocate(rhon(1)) !not used in Cartesian or 3D code so just set to something; could likely be left unallocated

allocate(xnall(lxnall))

allocate(ynall(lynall))

allocate(dnOall(lzn,lxnall,lynall),dnN2all(lzn,lxnall,lynall),dnO2all(lzn,lxnall,lynall),dvnrhoall(lzn,lxnall,lynall), &

dvnzall(lzn,lxnall,lynall),dvnxall(lzn,lxnall,lynall),dTnall(lzn,lxnall,lynall)) !ZZZ - note that these might be deallocated after each read to clean up memory management a bit...

!calculate the z grid (same for all) and distribute to workers so we can figure out their x-y slabs

print*, '...creating vertical grid and sending to workers...'

zn=[ ((real(izn,8)-1._wp)*dzn, izn=1,lzn) ] !root calculates and distributes but this is the same for all workers - assmes that the max neutral grid extent in altitude is always less than the plasma grid (should almost always be true)

maxzn=maxval(zn)

do iid=1,lid-1

call mpi_send(lzn,1,MPI_INTEGER,iid,taglz,MPI_COMM_WORLD,ierr)

call mpi_send(zn,lzn,mpi_realprec,iid,tagzn,MPI_COMM_WORLD,ierr)

end do

!Define a global neutral grid (input data) by assuming that the spacing is constant

ynall=[ ((real(iyn,8)-1._wp)*dyn, iyn=1,lynall) ]

meanyn=sum(ynall,1)/size(ynall,1)

ynall=ynall-meanyn !the neutral grid should be centered on zero for a cartesian interpolation

xnall=[ ((real(ixn,8)-1._wp)*dxn, ixn=1,lxnall) ]

meanxn=sum(xnall,1)/size(xnall,1)

xnall=xnall-meanxn !the neutral grid should be centered on zero for a cartesian interpolation

print *, 'Created full neutral grid with y,z extent:',minval(xnall),maxval(xnall),minval(ynall), &

maxval(ynall),minval(zn),maxval(zn)

!calculate the extent of my piece of the grid using max altitude specified for the neutral grid

call slabrange(maxzn,ximat,yimat,zimat,meanlat,xnrange,ynrange)

allocate(extents(0:lid-1,6),indx(0:lid-1,6),slabsizes(0:lid-1,2))

extents(0,1:6)=[0._wp,maxzn,xnrange(1),xnrange(2),ynrange(1),ynrange(2)]

!receive extents of each of the other workers: extents(lid,6)

print*, 'Receiving xn and yn ranges from workers...'

do iid=1,lid-1

call mpi_recv(xnrange,2,mpi_realprec,iid,tagxnrange,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

call mpi_recv(ynrange,2,mpi_realprec,iid,tagynrange,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

extents(iid,1:6)=[0._wp,maxzn,xnrange(1),xnrange(2),ynrange(1),ynrange(2)] !need to store values as xnrange overwritten for each worker

print*, 'Subgrid extents: ',iid,extents(iid,:)

end do

!find index into into neutral arrays for each worker: indx(lid,6)

print*, 'Root grid check: ',ynall(1),ynall(lynall)

print*, 'Converting ranges to indices...'

do iid=0,lid-1

call range2inds(extents(iid,1:6),zn,xnall,ynall,indices)

indx(iid,1:6)=indices

print*, 'Subgrid indices',iid,indx(iid,:)

end do

!send each worker the sizes for their particular chunk (all different) and send worker that grid chunk

print*,'Sending sizes and xn,yn subgrids to workers...'

do iid=1,lid-1

lxn=indx(iid,4)-indx(iid,3)+1

lyn=indx(iid,6)-indx(iid,5)+1

slabsizes(iid,1:2)=[lxn,lyn]

call mpi_send(lyn,1,MPI_INTEGER,iid,taglrho,MPI_COMM_WORLD,ierr)

call mpi_send(lxn,1,MPI_INTEGER,iid,taglx,MPI_COMM_WORLD,ierr)

allocate(xn(lxn),yn(lyn))

xn=xnall(indx(iid,3):indx(iid,4))

yn=ynall(indx(iid,5):indx(iid,6))

call mpi_send(xn,lxn,mpi_realprec,iid,tagxn,MPI_COMM_WORLD,ierr)

call mpi_send(yn,lyn,mpi_realprec,iid,tagyn,MPI_COMM_WORLD,ierr)

deallocate(xn,yn)

end do

!have root store its part to the full neutral grid

print*, 'Root is picking out its own subgrid...'

lxn=indx(0,4)-indx(0,3)+1

lyn=indx(0,6)-indx(0,5)+1

slabsizes(0,1:2)=[lxn,lyn]

allocate(xn(lxn),yn(lyn))

xn=xnall(indx(0,3):indx(0,4))

yn=ynall(indx(0,5):indx(0,6))

else !workers

!get teh z-grid from root so we know what the max altitude we have to deal with will be

call mpi_recv(lzn,1,MPI_INTEGER,0,taglz,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

allocate(zn(lzn))

call mpi_recv(zn,lzn,mpi_realprec,0,tagzn,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

maxzn=maxval(zn)

!calculate the extent of my grid

call slabrange(maxzn,ximat,yimat,zimat,meanlat,xnrange,ynrange)

!send ranges to root

call mpi_send(xnrange,2,mpi_realprec,0,tagxnrange,MPI_COMM_WORLD,ierr)

call mpi_send(ynrange,2,mpi_realprec,0,tagynrange,MPI_COMM_WORLD,ierr)

!receive my sizes from root, allocate then receive my pieces of the grid

call mpi_recv(lxn,1,MPI_INTEGER,0,taglx,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

call mpi_recv(lyn,1,MPI_INTEGER,0,taglrho,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

allocate(xn(lxn),yn(lyn))

call mpi_recv(xn,lxn,mpi_realprec,0,tagxn,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

call mpi_recv(yn,lyn,mpi_realprec,0,tagyn,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)

end if

!AT THIS POINT WE CAN ALLOCATE THE SUBGRID SIZES

allocate(dnO(lzn,lxn,lyn),dnN2(lzn,lxn,lyn),dnO2(lzn,lxn,lyn),dvnrho(lzn,lxn,lyn), &

dvnz(lzn,lxn,lyn),dvnx(lzn,lxn,lyn),dTn(lzn,lxn,lyn))

!PRINT OUT SOME BASIC INFO ABOUT THE GRID THAT WE'VE LOADED

if (debug) then

print *, 'Min/max zn,xn,yn values',myid,minval(zn),maxval(zn),minval(xn),maxval(xn),minval(yn),maxval(yn)

print *, 'Min/max zi,xi,yi values',myid,minval(zi),maxval(zi),minval(xi),maxval(xi),minval(yi),maxval(yi)

!print *, 'Source lat/long: ',myid,meanlat,meanlong

!print *, 'Plasma grid lat range: ',myid,minval(x%glat(:,:,:)),maxval(x%glat(:,:,:))

!print *, 'Plasma grid lon range: ',myid,minval(x%glon(:,:,:)),maxval(x%glon(:,:,:))

end if