!> mpi-related gemini functionality

module gemini3d_mpi

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

use, intrinsic :: iso_c_binding, only : c_ptr

use phys_consts, only: wp,debug

use mpimod, only: mpi_manualgrid, process_grid_auto, mpi_cfg, mpibreakdown, mpisetup, tag=>gemini_mpi, halo

use meshobj, only: curvmesh

use gemini3d_config, only: gemini_cfg

use io, only: output_plasma,output_aur,find_milestone,input_plasma,create_outdir

use potential_comm, only: get_BGEfields,velocities

use grid, only: lx1,lx2,lx3, grid_drift, read_grid, calc_subgrid_size

use collisions, only: conductivities

use potentialBCs_mumps, only: init_Efieldinput

use potential_comm,only : pot2perpfield, electrodynamics, BGfields_boundaries_root, BGfields_boundaries_worker

use neutral_perturbations, only: init_neutral_perturb,neutral_perturb,clear_neutral_perturb

use neutral_background, only: init_neutral_background, neutral_background_fileinput, neutral_background_empirical, &

clear_neutral_background_fileinput

use neutral, only: neutral_aggregate, neutral_wind_aggregate

use temporal_mpi, only : dt_comm

use advec_mpi, only: halo_interface_vels_allspec

use multifluid_mpi, only: halo_allparams, halo_fluidvars

use sources_mpi, only: RK2_prep_mpi_allspec, RK2_global_boundary_allspec

use ionization_mpi, only: get_gavg_Tinf

use gemini3d, only: fluidvar_pointers,fluidauxvar_pointers, electrovar_pointers, gemini_work, &

v2grid, v3grid, setv2v3, set_start_timefromcfg, init_precipinput_in, precip_perturb_in, &

solflux_perturb_in, init_solfluxinput_in, &

get_it, tneuBG, user_populate

use sanity_check, only : check_finite_perturb

!=======

! v2grid, v3grid, setv2v3, set_start_timefromcfg

!use gemini_work_def, only: gemini_work

!>>>>>>> 969ab4efaa53ca11dd996ac250d15542314487b1

implicit none (type, external)

private

public :: init_procgrid, outdir_fullgridvaralloc, read_grid_in, get_initial_state, &

BGfield_Lagrangian, check_dryrun, check_fileoutput, &

get_initial_drifts, init_inputdata_in, init_Efieldinput_in, pot2perpfield_in, init_neutralperturb_in, dt_select, &

neutral_perturb_in, electrodynamics_in, init_neutralbackground_in, &

halo_interface_vels_allspec_in, halo_allparams_in, &

RK2_prep_mpi_allspec_in,get_gavg_Tinf_in, clear_neutral_perturb_in, clear_neutral_background_in, &

mpisetup_in, mpiparms, calc_subgrid_size_in, &

RK2_global_boundary_allspec_in, halo_fluidvars_in, efield_perturb_in, inputdata_perturb_in

real(wp), parameter :: dtscale=2 ! controls how rapidly the time step is allowed to change

contains

!> call the mpi setup from our module

subroutine mpisetup_in()

call mpisetup()

end subroutine mpisetup_in

!> my id and the total number of processes for this run

subroutine mpiparms(myid,lid)

integer, intent(inout) :: myid,lid

myid=mpi_cfg%myid

lid=mpi_cfg%lid

end subroutine mpiparms

!> create output directory and allocate full grid potential storage

subroutine outdir_fullgridvaralloc(cfg,intvars,lx1,lx2all,lx3all)

type(gemini_cfg), intent(in) :: cfg

type(gemini_work), intent(inout) :: intvars

integer, intent(in) :: lx1,lx2all,lx3all

!> create a place, if necessary, for output datafiles

if (mpi_cfg%myid==0) then

call create_outdir(cfg)

allocate(intvars%Phiall(-1:lx1+2,-1:lx2all+2,-1:lx3all+2))

end if

end subroutine outdir_fullgridvaralloc

!> read in the grid and distribute to workers

subroutine read_grid_in(cfg,x)

type(gemini_cfg), intent(in) :: cfg

class(curvmesh), pointer, intent(inout) :: x

call read_grid(cfg%indatsize,cfg%indatgrid,cfg%flagperiodic, x)

!! read in a previously generated grid from filenames listed in input file

end subroutine read_grid_in

!> interface for setting simulation subgrid sizes for a particular worker in grid module

subroutine calc_subgrid_size_in(lx2all,lx3all)

integer, intent(in) :: lx2all, lx3all

call calc_subgrid_size(lx2all,lx3all)

end subroutine calc_subgrid_size_in

! !> load initial conditions and check if this is a restart run; set time variables accordingly

! subroutine get_initial_state(cfg,fluidvars,electrovars,intvars,x,UTsec,ymd,tdur)

! type(gemini_cfg), intent(inout) :: cfg

! real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

! real(wp), dimension(:,:,:,:), pointer, intent(inout) :: electrovars

! type(gemini_work), intent(inout) :: intvars

! class(curvmesh), intent(in) :: x

! real(wp), intent(inout) :: UTsec

! integer, dimension(3), intent(inout) :: ymd

! real(wp), intent(inout) :: tdur

!

! real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

! real(wp), dimension(:,:,:), pointer :: E1,E2,E3,J1,J2,J3,Phi

! integer, dimension(3) :: ymdtmp

! real(wp) :: UTsectmp,ttmp

! character(:), allocatable :: filetmp

!

! call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

! call electrovar_pointers(electrovars,E1,E2,E3,J1,J2,J3,Phi)

!

! call find_milestone(cfg, ttmp, ymdtmp, UTsectmp, filetmp)

! if ( ttmp > 0 ) then

! !! restart scenario

! if (mpi_cfg%myid==0) then

! print*, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'

! print*, '! Restarting simulation from time: ',ymdtmp,UTsectmp

! print*, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'

! end if

!

! !! Set start variables accordingly and read in the milestone

! UTsec=UTsectmp

! ymd=ymdtmp

!

! ! FIXME: instead keep tdur and just adjust the start time of the simulation to be closer to endtime

! tdur=cfg%tdur-ttmp ! subtract off time that has elapsed to milestone

! ! FIXME: need to feed in the t variable and overwrite it if we are restarting.

!

! if (mpi_cfg%myid==0) then

! print*, 'Treating the following file as initial conditions: ',filetmp

! print*, ' full duration: ',cfg%tdur,'; remaining simulation time: ',tdur ! FIXME: should be tdur-t once t adjusted

! end if

!

! if (tdur <= 1e-6_wp .and. mpi_cfg%myid==0) error stop 'Cannot restart simulation from the final time step!'

!

! cfg%tdur=tdur ! just to insure consistency

! call input_plasma(cfg%outdir, x%x1,x%x2all,x%x3all,cfg%indatsize,filetmp,ns,vs1,Ts,Phi,intvars%Phiall)

! else !! start at the beginning

! ! UTsec = cfg%UTsec0

! ! ymd = cfg%ymd0

! ! tdur = cfg%tdur

! call set_start_timefromcfg(cfg,ymd,UTsec,tdur)

!

! if (tdur <= 1e-6_wp .and. mpi_cfg%myid==0) error stop 'Simulation is of zero time duration'

! print*, 'Starting from beginning of simulation...'

! call input_plasma(cfg%outdir, x%x1,x%x2all,x%x3all,cfg%indatsize,cfg%indatfile,ns,vs1,Ts,Phi,intvars%Phiall)

! end if

! end subroutine get_initial_state

!> load initial conditions and check if this is a restart run; set time variables accordingly

subroutine get_initial_state(cfg,fluidvars,electrovars,intvars,x,UTsec,ymd,tdur,t,tmilestone)

type(gemini_cfg), intent(inout) :: cfg

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: electrovars

type(gemini_work), intent(inout) :: intvars

class(curvmesh), intent(in) :: x

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

integer, dimension(3), intent(inout) :: ymd

real(wp), intent(inout) :: tdur,t,tmilestone

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

real(wp), dimension(:,:,:), pointer :: E1,E2,E3,J1,J2,J3,Phi

integer, dimension(3) :: ymdtmp

real(wp) :: UTsectmp

character(:), allocatable :: filetmp

real(wp) :: tremaining

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call electrovar_pointers(electrovars,E1,E2,E3,J1,J2,J3,Phi)

call find_milestone(cfg, t, ymdtmp, UTsectmp, filetmp)

if ( t > 0 ) then

!! Set start variables accordingly and read in the milestone

UTsec=UTsectmp

ymd=ymdtmp

tremaining=cfg%tdur-t ! subtract off time that has elapsed to milestone

tdur=cfg%tdur

tmilestone=t

if (mpi_cfg%myid==0) then

print*, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'

print*, '! Restarting simulation from time: ',ymdtmp,UTsectmp

print*, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'

print*, 'Treating the following file as initial conditions: ',filetmp

print*, ' full duration: ',cfg%tdur,'; remaining simulation time: ',tremaining

print*, ' simulation start time: ',t

if (tremaining<=1e-6_wp) error stop 'Cannot restart simulation from the final time step!'

end if

!cfg%tdur=tdur ! just to insure consistency

call input_plasma(cfg%outdir, x%x1,x%x2all,x%x3all,cfg%indatsize,filetmp,ns,vs1,Ts,Phi,intvars%Phiall)

else !! start at the beginning

! UTsec = cfg%UTsec0

! ymd = cfg%ymd0

! tdur = cfg%tdur

call set_start_timefromcfg(cfg,ymd,UTsec,tdur)

t=0._wp

tdur=cfg%tdur

if (tdur <= 1e-6_wp .and. mpi_cfg%myid==0) error stop 'Simulation is of zero time duration'

print*, 'Starting from beginning of simulation...'

call input_plasma(cfg%outdir, x%x1,x%x2all,x%x3all,cfg%indatsize,cfg%indatfile,ns,vs1,Ts,Phi,intvars%Phiall)

end if

end subroutine get_initial_state

!> check whether file output should be done and complete it

subroutine check_fileoutput(cfg,fluidvars,electrovars,intvars,t,tout,tglowout,tmilestone,flagoutput,ymd,UTsec)

type(gemini_cfg), intent(in) :: cfg

real(wp), dimension(:,:,:,:), pointer, intent(in) :: fluidvars

real(wp), dimension(:,:,:,:), pointer, intent(in) :: electrovars

type(gemini_work), intent(inout) :: intvars

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

real(wp), intent(inout) :: tout,tglowout,tmilestone

integer, intent(inout) :: flagoutput

integer, dimension(3), intent(in) :: ymd

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

real(wp) :: tstart,tfin

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

real(wp), dimension(:,:,:), pointer :: E1,E2,E3,J1,J2,J3,Phi

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call electrovar_pointers(electrovars,E1,E2,E3,J1,J2,J3,Phi)

if (abs(t-tout) < 1d-5) then

tout = tout + cfg%dtout

if (cfg%nooutput ) then

if (mpi_cfg%myid==0) write(stderr,*) 'WARNING: skipping file output at sim time (sec)',t

return

endif

!! close enough to warrant an output now...

if (mpi_cfg%myid==0 .and. debug) call cpu_time(tstart)

!! We may need to adjust flagoutput if we are hitting a milestone

flagoutput=cfg%flagoutput

call user_populate(fluidvars,electrovars,intvars) ! custom user data

if (cfg%mcadence>0 .and. abs(t-tmilestone) < 1d-5) then

flagoutput=1 !force a full output at the milestone

call output_plasma(cfg%outdir,flagoutput,ymd, &

UTsec,vs2,vs3,ns,vs1,Ts,intvars%Phiall,J1,J2,J3, &

cfg%out_format,intvars%user_output)

tmilestone = t + cfg%dtout * cfg%mcadence

if(mpi_cfg%myid==0) print*, 'Milestone output triggered.'

else

call output_plasma(cfg%outdir,flagoutput,ymd, &

UTsec,vs2,vs3,ns,vs1,Ts,intvars%Phiall,J1,J2,J3, &

cfg%out_format,intvars%user_output)

end if

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

call cpu_time(tfin)

print *, 'Plasma output done for time step: ',t,' in cpu_time of: ',tfin-tstart

endif

end if

!> GLOW file output

if ((cfg%flagglow /= 0) .and. (abs(t-tglowout) < 1d-5)) then !same as plasma output

call cpu_time(tstart)

call output_aur(cfg%outdir, cfg%flagglow, ymd, UTsec, intvars%iver, cfg%out_format)

if (mpi_cfg%myid==0) then

call cpu_time(tfin)

print *, 'Auroral output done for time step: ',t,' in cpu_time of: ',tfin-tstart

end if

tglowout = tglowout + cfg%dtglowout

end if

end subroutine

!> check whether a dryrun simulation was done

subroutine check_dryrun(cfg)

type(gemini_cfg), intent(in) :: cfg

character(8) :: date

character(10) :: time

integer :: ierr

if (cfg%dryrun) then

ierr = mpibreakdown()

if (ierr /= 0) error stop 'Gemini dry run MPI shutdown failure'

call date_and_time(date,time)

print '(/,A)', 'DONE: ' // date(1:4) // '-' // date(5:6) // '-' // date(7:8) // 'T' &

// time(1:2) // ':' // time(3:4) // ':' // time(5:)

stop "OK: Gemini dry run"

endif

end subroutine check_dryrun

!> prep simulation for use of Lagrangian grid, if needed

subroutine BGfield_Lagrangian(cfg,x,electrovars,intvars)

type(gemini_cfg), intent(in) :: cfg

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: electrovars

type(gemini_work), intent(inout) :: intvars

real(wp) :: v2gridtmp,v3gridtmp

real(wp), dimension(:,:,:), allocatable :: E01,E02,E03

real(wp), dimension(:,:,:), pointer :: E1,E2,E3,J1,J2,J3,Phi

call electrovar_pointers(electrovars,E1,E2,E3,J1,J2,J3,Phi) ! bind pointers for EM variables

allocate(E01(lx1,lx2,lx3),E02(lx1,lx2,lx3),E03(lx1,lx2,lx3)) ! allocate without ghost cells

E01=0; E02=0; E03=0;

if (cfg%flagE0file==1) then

call get_BGEfields(x,E01,E02,E03,intvars%efield)

end if

if (cfg%flaglagrangian) then ! Lagrangian (moving) grid; compute from input background electric fields

call grid_drift(x,E02,E03,v2gridtmp,v3gridtmp)

call setv2v3(v2gridtmp,v3gridtmp)

! must recompute wind data since we've adjust grid drift

call neutral_wind_aggregate(v2grid,v3grid,intvars%atmos,intvars%atmosperturb,.true.)

if (mpi_cfg%myid==0) print*, mpi_cfg%myid,' using Lagrangian grid moving at: ',v2grid,v3grid

else ! stationary grid

v2gridtmp=0._wp

v3gridtmp=0._wp

call setv2v3(v2gridtmp,v3gridtmp)

! no neutral adjustment needed if grid isn't moving

E1(1:lx1,1:lx2,1:lx3) = E1(1:lx1,1:lx2,1:lx3) + E01

E2(1:lx1,1:lx2,1:lx3) = E2(1:lx1,1:lx2,1:lx3) + E02

E3(1:lx1,1:lx2,1:lx3) = E3(1:lx1,1:lx2,1:lx3) + E03

end if

deallocate(E01,E02,E03)

if (mpi_cfg%myid==0) then

print*, 'Recomputed initial fields including background (BGfieldLagrangian):'

print*, ' ',minval(E1(1:lx1,1:lx2,1:lx3)),maxval(E1(1:lx1,1:lx2,1:lx3))

print*, ' ',minval(E2(1:lx1,1:lx2,1:lx3)),maxval(E2(1:lx1,1:lx2,1:lx3))

print*, ' ',minval(E3(1:lx1,1:lx2,1:lx3)),maxval(E3(1:lx1,1:lx2,1:lx3))

end if

end subroutine BGfield_Lagrangian

!> initial drifts at the start of the simulation

subroutine get_initial_drifts(cfg,x,fluidvars,fluidauxvars,electrovars,intvars)

type(gemini_cfg), intent(in) :: cfg

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

real(wp), dimension(:,:,:,:), pointer, intent(in) :: fluidauxvars

real(wp), dimension(:,:,:,:), pointer, intent(in) :: electrovars

type(gemini_work), intent(in) :: intvars

real(wp), dimension(:,:,:), allocatable :: sig0,sigP,sigH,sigPgrav,sigHgrav

real(wp), dimension(:,:,:,:), allocatable :: muP,muH,nusn

integer :: lx1,lx2,lx3,lsp

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

real(wp), dimension(:,:,:,:), pointer :: rhovs1,rhoes

real(wp), dimension(:,:,:), pointer :: rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom

real(wp), dimension(:,:,:),pointer :: E1,E2,E3,J1,J2,J3,Phi

! bind pointers

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call fluidauxvar_pointers(fluidauxvars,rhovs1,rhoes,rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom)

call electrovar_pointers(electrovars,E1,E2,E3,J1,J2,J3,Phi)

! calculate drifts

lx1=x%lx1; lx2=x%lx2; lx3=x%lx3; lsp=size(ns,4);

allocate(sig0(lx1,lx2,lx3),sigP(lx1,lx2,lx3),sigH(lx1,lx2,lx3),sigPgrav(lx1,lx2,lx3),sigHgrav(lx1,lx2,lx3))

allocate(muP(lx1,lx2,lx3,lsp),muH(lx1,lx2,lx3,lsp),nusn(lx1,lx2,lx3,lsp))

call conductivities(intvars%atmos%nn,intvars%atmos%Tn,ns,Ts,vs1,B1,sig0,sigP,sigH,muP,muH,nusn,sigPgrav,sigHgrav)

call velocities(muP,muH,nusn,E2,E3,intvars%atmos%vn2,intvars%atmos%vn3,ns,Ts,x, &

cfg%flaggravdrift,cfg%flagdiamagnetic,vs2,vs3)

deallocate(sig0,sigP,sigH,muP,muH,nusn,sigPgrav,sigHgrav)

if(mpi_cfg%myid==0) then

print*, 'Recomputed initial drifts: '

print*, ' ',minval(vs2(1:lx1,1:lx2,1:lx3,1:lsp)),maxval(vs2(1:lx1,1:lx2,1:lx3,1:lsp))

print*, ' ',minval(vs3(1:lx1,1:lx2,1:lx3,1:lsp)),maxval(vs3(1:lx1,1:lx2,1:lx3,1:lsp))

end if

end subroutine get_initial_drifts

!> initialize the process gridf for this simulation

subroutine init_procgrid(lx2all,lx3all,lid2in,lid3in)

integer, intent(in) :: lx2all,lx3all,lid2in,lid3in

if (lid2in==-1) then

call process_grid_auto(lx2all, lx3all)

!! grid_size defines lx2all and lx3all

else

call mpi_manualgrid(lx2all, lx3all, lid2in, lid3in)

endif

print '(A, I0, A1, I0)', 'process grid (Number MPI processes) x2, x3: ',mpi_cfg%lid2, ' ', mpi_cfg%lid3

print '(A, I0, A, I0, A1, I0)', 'Process:',mpi_cfg%myid,' at process grid location: ',mpi_cfg%myid2,' ',mpi_cfg%myid3

end subroutine init_procgrid

!> initialize all inputdata classes (internals will check whether present in config file)

subroutine init_inputdata_in(cfg,x,dt,t,ymd,UTsec,intvars)

type(gemini_cfg), intent(in) :: cfg

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

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

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

integer, dimension(3), intent(in) :: ymd

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

type(gemini_work), intent(inout) :: intvars

call init_neutralbackground_in(dt,cfg,x,intvars,ymd,UTsec)

call init_precipinput_in(cfg,x,dt,t,ymd,UTsec,intvars)

call init_Efieldinput_in(cfg,x,dt,intvars,ymd,UTsec)

call init_neutralperturb_in(dt,cfg,x,intvars,ymd,UTsec)

call init_solfluxinput_in(cfg,x,dt,t,ymd,UTsec,intvars)

!> add other inputdata streams here...

end subroutine init_inputdata_in

!> initialize electric field input data

subroutine init_Efieldinput_in(cfg,x,dt,intvars,ymd,UTsec)

type(gemini_cfg), intent(in) :: cfg

class(curvmesh), intent(in) :: x

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

type(gemini_work), intent(inout) :: intvars

integer, dimension(3), intent(in) :: ymd

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

call init_Efieldinput(dt,cfg,ymd,UTsec,x,intvars%efield)

end subroutine init_Efieldinput_in

!> convert potential to electric field by differentiating

subroutine pot2perpfield_in(x,electrovars)

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: electrovars

real(wp), dimension(:,:,:), pointer :: E1,E2,E3,J1,J2,J3,Phi

call electrovar_pointers(electrovars,E1,E2,E3,J1,J2,J3,Phi)

E1 = 0

call pot2perpfield(Phi,x,E2,E3)

if(mpi_cfg%myid==0) then

print '(A)', 'Recomputed initial dist. fields via pot2perpfield:'

print*, ' gemini ',minval(E1(1:lx1,1:lx2,1:lx3)),maxval(E1(1:lx1,1:lx2,1:lx3))

print*, ' gemini ',minval(E2(1:lx1,1:lx2,1:lx3)),maxval(E2(1:lx1,1:lx2,1:lx3))

print*, ' gemini ',minval(E3(1:lx1,1:lx2,1:lx3)),maxval(E3(1:lx1,1:lx2,1:lx3))

print*, ' gemini ',minval(Phi(1:lx1,1:lx2,1:lx3)),maxval(Phi(1:lx1,1:lx2,1:lx3))

end if

end subroutine pot2perpfield_in

!> initialize neutral perturbations

subroutine init_neutralperturb_in(dt,cfg,x,intvars,ymd,UTsec)

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

type(gemini_cfg), intent(in) :: cfg

class(curvmesh), intent(in) :: x

type(gemini_work), intent(inout) :: intvars

integer, dimension(3), intent(in) :: ymd

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

call init_neutral_perturb(cfg,x,dt,ymd,UTsec,intvars%atmosperturb)

call neutral_aggregate(v2grid,v3grid,intvars%atmos,intvars%atmosperturb)

end subroutine init_neutralperturb_in

!> initialize neutral background information

subroutine init_neutralbackground_in(dt,cfg,x,intvars,ymd,UTsec)

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

type(gemini_cfg), intent(in) :: cfg

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

type(gemini_work), intent(inout) :: intvars

integer, dimension(3), intent(in) :: ymd

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

call init_neutral_background(dt,cfg,ymd,UTsec,x,v2grid,v3grid,intvars%atmos,intvars%atmosbackground)

call neutral_aggregate(v2grid,v3grid,intvars%atmos,intvars%atmosperturb) ! some problems can occur if no data are populated into neutral arrays used by rest of code

end subroutine init_neutralbackground_in

!> select time step and throttle if changing too rapidly

subroutine dt_select(cfg,x,fluidvars,fluidauxvars,t,tout,tglowout,dt)

type(gemini_cfg), intent(in) :: cfg

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(in) :: fluidvars

real(wp), dimension(:,:,:,:), pointer, intent(in) :: fluidauxvars

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

real(wp), intent(inout) :: dt

real(wp) :: dtprev

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

real(wp), dimension(:,:,:,:), pointer :: rhovs1,rhoes

real(wp), dimension(:,:,:), pointer :: rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom

! bind pointers for fluid and auxiliary variables

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call fluidauxvar_pointers(fluidauxvars,rhovs1,rhoes,rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom)

!> save prior time step

dtprev = dt

!> time step calculation, requires workers to report their most stringent local stability constraint

call dt_comm(t,tout,tglowout,cfg,ns,Ts,vs1,vs2,vs3,B1,B2,B3,x,dt)

!> do not allow the time step to change too rapidly

if (get_it()>1) then

if(dt/dtprev > dtscale) then

!! throttle how quickly we allow dt to increase

dt=dtscale*dtprev

if (mpi_cfg%myid == 0) then

print '(A,EN14.3)', 'Throttling dt to: ',dt

end if

end if

end if

end subroutine dt_select

!> update all perturbed inputdata quantities

subroutine inputdata_perturb_in(cfg,intvars,x,dt,t,ymd,UTsec)

type(gemini_cfg), intent(in) :: cfg

type(gemini_work), intent(inout) :: intvars

class(curvmesh), intent(inout) :: x ! unit vectors could be deallocated in this procedure

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

integer, dimension(3), intent(in) :: ymd

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

call neutral_background_in(cfg,intvars,x,dt,t,ymd,UTsec)

call neutral_perturb_in(cfg,intvars,x,dt,t,ymd,UTsec)

call efield_perturb_in(cfg,intvars,x,dt,t,ymd,UTsec)

call precip_perturb_in(cfg,intvars,x,dt,t,ymd,UTsec)

call solflux_perturb_in(cfg,intvars,x,dt,t,ymd,UTsec)

end subroutine inputdata_perturb_in

!> update the solar flux inputdata if present

subroutine neutral_background_in(cfg,intvars,x,dt,t,ymd,UTsec)

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

type(gemini_cfg), intent(in) :: cfg

integer, dimension(3), intent(in) :: ymd

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

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

type(gemini_work), intent(inout) :: intvars

if (cfg%flagneutralBGfile==1) then

!print*, 'File-based neutral background...'

call neutral_background_fileinput(dt,t,cfg,ymd,UTsec,x,intvars%atmos,intvars%atmosbackground) ! load into base array variables

call neutral_aggregate(v2grid,v3grid,intvars%atmos,intvars%atmosperturb)

else

!> get neutral background

! if ( get_it()/=1 .and. cfg%flagneuBG .and. t>tneuBG) then

if ( get_it()==1 .or. (cfg%flagneuBG .and. t>tneuBG) ) then

!^we dont' throttle for tneuBG so we have to do things this way to not skip over...

!call cpu_time(tstart)

call neutral_background_empirical(cfg,ymd,UTsec,x,v2grid,v3grid,intvars%atmos) ! load background states from empirical models into base array variables

call neutral_aggregate(v2grid,v3grid,intvars%atmos,intvars%atmosperturb) ! apply to variables in this program unit

tneuBG=tneuBG+cfg%dtneuBG

!if (myid==0) then

! call cpu_time(tfin)

! print *, 'Neutral background at time: ',t,' calculated in time: ',tfin-tstart

!end if

end if

end if

end subroutine neutral_background_in

!> compute neutral perturbations and apply to main code variables

subroutine neutral_perturb_in(cfg,intvars,x,dt,t,ymd,UTsec)

type(gemini_cfg), intent(in) :: cfg

type(gemini_work), intent(inout) :: intvars

class(curvmesh), intent(inout) :: x ! unit vectors could be deallocated in this procedure

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

integer, dimension(3), intent(in) :: ymd

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

if (cfg%flagdneu==1) then

call neutral_perturb(cfg,dt,t,ymd,UTsec,x,v2grid,v3grid,intvars%atmos,intvars%atmosperturb)

call check_finite_perturb(cfg%outdir, t, mpi_cfg%myid, intvars%atmos%nn, intvars%atmos%Tn, &

intvars%atmos%vn1, intvars%atmos%vn2, intvars%atmos%vn3)

call neutral_aggregate(v2grid,v3grid,intvars%atmos,intvars%atmosperturb) ! extra step to tell GEMINI to agg the perturbations and BG

end if

!> no default for neutral perturbations -- aggregate functions will ignore if not flagged.

end subroutine neutral_perturb_in

!> compute boundary conditions for electric field solutions

subroutine efield_perturb_in(cfg,intvars,x,dt,t,ymd,UTsec)

type(gemini_cfg), intent(in) :: cfg

type(gemini_work), intent(inout) :: intvars

class(curvmesh), intent(inout) :: x ! unit vectors could be deallocated in this procedure

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

integer, dimension(3), intent(in) :: ymd

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

!> assign values to the background fields, etc., irrespective of whether or not we do a potential solve

if (mpi_cfg%myid/=0) then

call BGfields_boundaries_worker(intvars%flagdirich,intvars%E01,intvars%E02,intvars%E03, &

intvars%Vminx1slab,intvars%Vmaxx1slab)

else

call BGfields_boundaries_root(dt,t,ymd,UTsec,cfg,x,intvars%efield, &

intvars%flagdirich,intvars%Vminx1,intvars%Vmaxx1,intvars%Vminx2,intvars%Vmaxx2,&

intvars%Vminx3,intvars%Vmaxx3,intvars%E01,intvars%E02,intvars%E03, &

intvars%Vminx1slab,intvars%Vmaxx1slab)

end if

end subroutine efield_perturb_in

!> call electrodynamics solution

subroutine electrodynamics_in(cfg,fluidvars,fluidauxvars,electrovars,intvars,x,t,dt,ymd,UTsec)

type(gemini_cfg), intent(in) :: cfg

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

real(wp), dimension(:,:,:,:), pointer, intent(in) :: fluidauxvars

real(wp), dimension(:,:,:,:), pointer, intent(in) :: electrovars

type(gemini_work), intent(inout) :: intvars

class(curvmesh), intent(in) :: x

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

integer, dimension(3), intent(in) :: ymd

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

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

real(wp), dimension(:,:,:,:), pointer :: rhovs1,rhoes

real(wp), dimension(:,:,:), pointer :: rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom

real(wp), dimension(:,:,:),pointer :: E1,E2,E3,J1,J2,J3,Phi

! bind pointers

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call fluidauxvar_pointers(fluidauxvars,rhovs1,rhoes,rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom)

call electrovar_pointers(electrovars,E1,E2,E3,J1,J2,J3,Phi)

! E&M solves

call electrodynamics(get_it(),t,dt,intvars%atmos%nn,intvars%atmos%vn2, &

intvars%atmos%vn3,intvars%atmos%Tn, &

cfg,ns,Ts,vs1,B1,vs2,vs3,x,intvars%efield,E1,E2,E3,J1,J2,J3, &

intvars%Phiall,intvars%flagdirich, &

intvars%Vminx1,intvars%Vmaxx1,intvars%Vminx2,intvars%Vmaxx2, &

intvars%Vminx3,intvars%Vmaxx3,intvars%Vminx1slab,intvars%Vmaxx1slab, &

intvars%E01,intvars%E02,intvars%E03, &

ymd,UTsec,intvars%sig0,intvars%sigP,intvars%sigH,intvars%sigNCP,intvars%sigNCH)

end subroutine electrodynamics_in

! FIXME: deprecated; is easier/better just to halo all params prior to advection and interface vels

! calculation...

!> haloing for computing cell interface velocities

subroutine halo_interface_vels_allspec_in(x,fluidvars,lsp)

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

integer, intent(in) :: lsp

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call halo_interface_vels_allspec(x%flagper,vs2,vs3,lsp)

end subroutine halo_interface_vels_allspec_in

!> halo all ***advected*** parameters

subroutine halo_allparams_in(x,fluidvars,fluidauxvars)

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidauxvars

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

real(wp), dimension(:,:,:,:), pointer :: rhovs1,rhoes

real(wp), dimension(:,:,:), pointer :: rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom

! bind pointers

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call fluidauxvar_pointers(fluidauxvars,rhovs1,rhoes,rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom)

! halo all fluid parameters

call halo_allparams(ns,rhovs1,rhoes,x%flagper)

end subroutine halo_allparams_in

!> halo all parameters, including velocities

subroutine halo_fluidvars_in(x,fluidvars,fluidauxvars)

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidauxvars

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

real(wp), dimension(:,:,:,:), pointer :: rhovs1,rhoes

real(wp), dimension(:,:,:), pointer :: rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom

! bind pointers

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call fluidauxvar_pointers(fluidauxvars,rhovs1,rhoes,rhov2,rhov3,B1,B2,B3,v1,v2,v3,rhom)

! halo all fluid parameters

call halo_fluidvars(ns,rhovs1,rhoes,vs2,vs3,x%flagper)

end subroutine halo_fluidvars_in

!> prepare and then halo data for compression substep

subroutine RK2_prep_mpi_allspec_in(x,fluidvars)

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

real(wp), dimension(-1:x%lx1+2,-1:x%lx2+2,-1:x%lx3+2) :: param

integer :: isp,lsp

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

lsp=size(ns,4)

call RK2_prep_mpi_allspec(vs1,vs2,vs3,x%flagper)

end subroutine RK2_prep_mpi_allspec_in

!> prepare data for compression substep deal with global boundaries

subroutine RK2_global_boundary_allspec_in(x,fluidvars)

class(curvmesh), intent(in) :: x

real(wp), dimension(:,:,:,:), pointer, intent(inout) :: fluidvars

real(wp), dimension(:,:,:,:), pointer :: ns,vs1,vs2,vs3,Ts

call fluidvar_pointers(fluidvars,ns,vs1,vs2,vs3,Ts)

call RK2_global_boundary_allspec(vs1,vs2,vs3,x%flagper)

end subroutine RK2_global_boundary_allspec_in

!> agree on average value of gravity and exospheric temp

subroutine get_gavg_Tinf_in(intvars)

type(gemini_work), intent(inout) :: intvars

call get_gavg_Tinf(intvars%atmos,intvars%gavg,intvars%Tninf)

end subroutine get_gavg_Tinf_in

!> deallocate module storage for neutral perturbations

subroutine clear_neutral_perturb_in(intvars)

type(gemini_work), intent(inout) :: intvars

call clear_neutral_perturb(intvars%atmosperturb)

end subroutine clear_neutral_perturb_in

!> deallocate module storage for neutral background

subroutine clear_neutral_background_in(intvars)

type(gemini_work), intent(inout) :: intvars

call clear_neutral_background_fileinput(intvars%atmosbackground)

end subroutine clear_neutral_background_in

end module gemini3d_mpi