JHUAPL/kaiju
1
0
Fork 0
mirror of https://github.com/JHUAPL/kaiju.git synced 2026-01-08 02:23:51 -05:00
Code Issues Packages Projects Releases 1 Wiki Activity
Files
22607a408429c23d825cde1f9f6aa1eccfecc7ea
kaiju/src/raiju/raijuStarter.F90

602 lines
26 KiB
Fortran
Raw Blame History

! Initialize kaimag
module raijustarter
! Base
use shellgrid
use xml_input
use planethelper
use arrayutil
! Raiju
use raijudefs
use raijutypes
use raijugrids
use raijuetautils
use raijuRecon, only : maxOrderSupported
use raijuout
use raijuICHelpers
use raijulosses, only : initRaiLosses
use raijuColdStartHelper, only : initRaijuColdStarter
! Cmake points to this
use raijuuseric
implicit none
contains
!------
! Main Initialization Routines
!------
subroutine raijuInit(app, inpXML, shGridO)
class(raijuApp_T), intent(inout) :: app
type(XML_Input_T), intent(in) :: inpXML
type(ShellGrid_T), intent(in), optional :: shGridO
!! If provided, and xml arg is set, we make a child grid for ourselves off of this one
type(XML_Input_T) :: iXML, iXML_volt
character(len=strLen) :: tmpStr
! Make sure root is Kaiju/raiju
call inpXML%GetFileStr(tmpStr)
! Create new XML reader w/ RAIJU as root
iXML = New_XML_Input(trim(tmpStr),'Kaiju/RAIJU',.true.)
! Init model, grid, state
call raijuInitModel(app%Model, iXML)
if(present(shGridO)) then
call raijuInitGrid(app%Model, app%Grid, app%opt, iXML, shGridO)
else
call raijuInitGrid(app%Model, app%Grid, app%opt, iXML)
endif
call raijuInitState(app%Model,app%Grid,app%State,iXML)
! Do losses after everything else has been set, just in case they need something from it
call initRaiLosses(app%Model, app%Grid, app%State, iXML)
! Initialize IOCLOCK
if (app%Model%isSA) then
! If we are driving, we look to raiju section for IO timing
call app%State%IO%init(iXML,app%State%t,app%State%ts)
else
! If we are not driving, we look to voltron section for IO timing
iXML_volt = New_XML_Input(trim(tmpStr),'Kaiju/VOLTRON',.true.)
call app%State%IO%init(iXML_volt,app%State%t,app%State%ts)
endif
! TODO: Add some final checks before returning
! e.g., if we are mapping fluids, are we mapping something to HOTP?
! (totally just thought of it randomly and not because I spent 4 hours realizing this was my problem)
end subroutine raijuInit
! Sets up Model; Grid and State must be set up separately
subroutine raijuInitModel(Model, iXML)
type(raijuModel_T), intent(inout) :: Model
type(XML_Input_T) , intent(in) :: iXML
character(len=strLen) :: tmpStr
character(len=strLen) :: tmpResId
real(rp) :: cfl0
write(*,*) "raijuInitModel is starting"
! Assume we are running with others, but allow for solo run
call iXML%Set_Val(Model%isSA, "driver/isSA", .false.)
! Assuming that if being controlled by e.g. Voltron, someone else will set RunID accordingly
! If we are in SA mode, need to set it ourselves
if (Model%isSA) then
call iXML%Set_Val(Model%RunID, "prob/RunID","raijuSA") ! raiju stand-alone
else
call iXML%Set_Val(Model%RunID, "/Kaiju/Gamera/sim/runid","msphere")
endif
! Timing info, if provided
call iXML%Set_Val(Model%t0 ,'time/T0',0.0)
call iXML%Set_Val(Model%tFin,'time/tFin',60.0)
call iXML%Set_Val(Model%dt ,'time/dt',1.0)
! Config file
call iXML%Set_Val(Model%configFName, "config/fname","raijuconfig.h5")
call CheckFileOrDie(Model%configFName,"RAIJU unable to open config file")
! Solar wind filename
call iXML%Set_Val(Model%tsF, '/Kaiju/Gamera/wind/tsFile', 'bcwind.h5')
call iXML%Set_Val(Model%isMPI, "mpi/isMPI",.false.)
if (Model%isMPI) then
write(*,*) "(RAIJU) MPI not implemented yet, dying."
stop
endif
!--- Restarts ---!
if (Model%isSA) then
tmpStr = "restart/doRes"
else
tmpStr = "/Kaiju/Gamera/restart/doRes"
endif
call iXML%Set_Val(Model%isRestart, trim(tmpStr),.false.)
if (Model%isRestart) then
! Restart number
if (Model%isSA) then
tmpStr = "restart/nRes"
else
tmpStr = "/Kaiju/Gamera/restart/nRes"
endif
call iXML%Set_Val(Model%nResIn, trim(tmpStr), Model%nResIn)
! Restart id (allows it to be different from our current runId)
if (Model%isSA) then
tmpStr = "restart/resId"
else
tmpStr = "/Kaiju/Gamera/restart/resId"
endif
call iXML%Set_Val(Model%resId, trim(tmpStr), Model%RunID)
call genResInFname(Model, Model%ResF, runIdO=Model%resId) ! Determine filename to read from
endif
!--- Plasmasphere ---!
call iXML%Set_Val(Model%doPlasmasphere, "plasmasphere/doPsphere",.true.)
call iXML%Set_Val(Model%doPsphEvol, 'plasmasphere/doEvol',.true.)
call iXML%Set_Val(Model%psphEvolRad, 'plasmasphere/evolRad', def_psphEvolRad)
! Determine number of species. First set default, then read from xml to overwrite if present
if (Model%doPlasmasphere) then
Model%nSpc = 3
call iXML%Set_Val(Model%psphInitKp, "plasmasphere/initKp",3.0)
else
Model%nSpc = 2
endif
call iXML%Set_Val(Model%doExcessToPsph, "prob/doExcessMap",.true.)
!! Allow mapping of excess H+ to plasmasphere channel
if (Model%doExcesstoPsph .and. .not. Model%doPlasmasphere) then
write(*,*)"Error: Can't map excess to psph because psph doesn't exist!"
endif
call iXML%Set_Val(Model%nSpc, "prob/nSpc",Model%nSpc)
!--- Domain ---!
call iXML%Set_Val(Model%n_bndLim , "domain/bndLim" , 3)
call iXML%Set_Val(Model%maxTail_buffer, "domain/tail_buffer", def_maxTail_buffer)
call iXML%Set_Val(Model%maxSun_buffer , "domain/sun_buffer" , def_maxSun_buffer)
call iXML%Set_Val(Model%maxTail_active, "domain/tail_active", def_maxTail_active)
call iXML%Set_Val(Model%maxSun_active , "domain/sun_active" , def_maxSun_active)
call iXML%Set_Val(Model%activeDomRad , "domain/activeRad" , def_activeDomRad)
! Store all distances as positive values, we'll add signs as needed later
Model%maxTail_buffer = abs(Model%maxTail_buffer)
Model%maxSun_buffer = abs(Model%maxSun_buffer)
Model%maxTail_active = abs(Model%maxTail_active)
Model%maxSun_active = abs(Model%maxSun_active)
Model%activeDomRad = abs(Model%activeDomRad)
!---Solver ---!
call iXML%Set_Val(Model%doSmoothGrads,'sim/doSmoothGrads',def_doSmoothGrads)
call iXML%Set_Val(Model%doUseVelLRs,'sim/useVelLRs',def_doUseVelLRs)
call iXML%Set_Val(Model%maxItersPerSec,'sim/maxIter',def_maxItersPerSec)
call iXML%Set_Val(Model%maxOrder,'sim/maxOrder',7)
if (Model%maxOrder > maxOrderSupported) then
write(*,*) "(RAIJU) Too much order, allowable orders <= ",maxOrderSupported
stop
endif
call iXML%Set_Val(Model%PDMB,'sim/pdmb',def_pdmb)
cfl0 = min(0.5/(Model%PDMB+0.5), cflMax) !Set CFL based on PDM
!Set CFL from XML
call iXML%Set_Val(Model%CFL ,'sim/cfl' ,cfl0)
if (Model%CFL > cfl0) then
write(*,*) '-------------------------------------'
write(*,*) '(RAIJU) WARNING, CFL is above critical value!'
write(*,*) 'CFL/Critical/PDMB = ', Model%CFL,cfl0,Model%PDMB
write(*,*) '-------------------------------------'
else
write(*,*) 'CFL # = ', Model%CFL
endif
! Certain physical constants that shouldn't be constants
call iXML%Set_Val(Model%tiote, "prob/tiote",4.0)
! Active shell settings
call iXML%Set_Val(Model%doActiveShell, "activeShell/doAS",.false.)
call iXML%Set_Val(Model%worthyFrac, "activeShell/worthyFrac",fracWorthyDef)
! Corotation
if (Model%isSA) then
call iXML%Set_Val(Model%doOwnCorot, "prob/doCorot",.true.)
else
call iXML%Set_Val(Model%doOwnCorot, "prob/doCorot",.false.)
endif
if (Model%doOwnCorot) then
call iXML%Set_Val(Model%doGeoCorot, "prob/doGeoCorot",.false.)
endif
! Lambda channel settings
call iXML%Set_Val(Model%doDynamicLambdaRanges, "lambdas/dynamicRanges",.false.)
! Determine which DP2eta mapping should be used (e.g. Maxwellian, Kappa, user-defined)
call iXML%Set_Val(tmpStr, "moments/DP2EtaMap","MAXWELL")
select case(tmpStr)
case("MAXWELL")
Model%dp2etaMap => Maxwell2Eta
case("KAPPA")
Model%dp2etaMap => Kappa2Eta
call iXML%Set_Val(Model%kappa, "moments/kappa",6.0)
! TODO: Add user option which will point to whatever is in the chosen IC file
case DEFAULT
write(*,*) "(RAIJU) Received unavailable DP2Eta mapping: ",tmpStr
write(*,*) " Dying."
stop
end select
!--- Losses ---!
call iXML%Set_Val(Model%doLosses, "losses/doLosses",.true.)
!call iXML%Set_Val(Model%doSS , "losses/doSS" ,.true. )
call iXML%Set_Val(Model%doCC , "losses/doCC" ,.true. )
call iXML%Set_Val(Model%doCX , "losses/doCX" ,.true.)
!call iXML%Set_Val(Model%doFLC , "losses/doFLC",.false.)
call iXML%Set_Val(Model%doEWM , "losses/doEWM",.true.)
call iXML%Set_Val(Model%ewmType , "losses/ewmType","SS")
!--- Output ---!
call iXML%Set_Val(Model%isLoud , "output/loudConsole",.false.)
call iXML%Set_Val(Model%writeGhosts , "output/writeGhosts",.false.)
call iXML%Set_Val(Model%doOutput_3DLoss , "output/doLossExtras" ,.false.) ! Several (Ni,Nj,Nk) arrays
call iXML%Set_Val(Model%doOutput_debug , "output/doDebug" ,.false.)
! Model Hax?
call iXML%Set_Val(Model%doHack_rcmEtaPush, "hax/rcmEtaPush",.false.)
! Determine coupling information
call setMHD2RaiInfo(Model, iXML)
! Set planet params
call getPlanetParams(Model%planet, iXML, doLoudO=.true.)
end subroutine raijuInitModel
subroutine setMHD2RaiInfo(Model, iXML)
type(raijuModel_T), intent(inout) :: Model
type(XML_Input_T), intent(in) :: iXML
integer :: nFluids
integer :: f
character(len=strLen) :: xmlFluidTag
! Domain determination
call iXML%Set_Val(Model%lim_vaFrac_soft, "domain/vaFrac_soft", def_lim_vaFrac_soft)
call iXML%Set_Val(Model%lim_vaFrac_hard, "domain/vaFrac_hard", def_lim_vaFrac_hard)
call iXML%Set_Val(Model%lim_bmin_soft , "domain/bmin_soft" , def_lim_bmin_soft )
call iXML%Set_Val(Model%lim_bmin_hard , "domain/bmin_hard" , def_lim_bmin_hard )
!call iXML%Set_Val(Model%vaFracThresh , "cpl/vaFracThresh" , def_vaFracThresh)
!call iXML%Set_Val(Model%bminThresh , "cpl/bminThresh" , def_bminThresh)
!call iXML%Set_Val(Model%PstdThresh , "cpl/Pstd" , def_PstdThresh)
!call iXML%Set_Val(Model%normAngThresh, "cpl/normAngThresh", def_normAngle)
! Note: Not used now within raiju, just passing through. But should probably include as a constraint
call iXML%Set_Val(Model%nBounce, "cpl/nBounce", def_nBounce)
! Convert degrees to dot product value
!Model%normAngThresh = cos(Model%normAngThresh*PI/180.0)
! Fluid mapping
call iXML%Set_Val(nFluids, "cpl/nFluidsIn",0)
Model%nFluidIn = nFluids
if (nFluids > 0) then
allocate(Model%fluidInMaps(nFluids))
do f=1,nFluids
write(xmlFluidTag,'(A,I0,A)')"fluidIn",f,"/imhd"
call iXML%Set_Val(Model%fluidInMaps(f)%idx_mhd, xmlFluidTag, 0)
write(xmlFluidTag,'(A,I0,A)')"fluidIn",f,"/flav"
call iXML%Set_Val(Model%fluidInMaps(f)%flav, xmlFluidTag, -1)
write(xmlFluidTag,'(A,I0,A)')"fluidIn",f,"/excessToPsph"
call iXML%Set_Val(Model%fluidInMaps(f)%doExcessToPsph, xmlFluidTag, .false.)
enddo
endif
end subroutine setMHD2RaiInfo
subroutine raijuInitGrid(Model, Grid, opts, iXML, shGridO)
type(raijuModel_T) , intent(inout) :: Model
type(raijuGrid_T) , intent(inout) :: Grid
type(raijuOptions_T), intent(in) :: opts
type(XML_Input_T), intent(in) :: iXML
type(ShellGrid_T), optional, intent(in) :: shGridO
character(len=strLen) :: tmpStr
type(ShellGrid_T) :: shGrid
logical :: doHack_constB
! Set grid params
call iXML%Set_Val(Grid%nB, "grid/Nbnd", 4 ) ! Number of cells between open boundary and active domain
! If we are restarting, use shellGrid from file
if (Model%isRestart) then
tmpStr = "RESTART"
else
! Otherwise, ask user how we should be making our grid
call iXML%Set_Val(tmpStr, "grid/gType","UNISPH")
endif
! Fill out Grid object depending on chosen method
select case(tmpStr)
case("RESTART")
Grid%gType = RAI_G_SHGRID ! Idk, not important right now
call GenShellGridFromFile(Grid%shGrid, RAI_SG_NAME, Model%ResF)
case("UNISPH")
Grid%gType = RAI_G_UNISPH
! Generate our own grid from scratch
call raijuGenUniSphGrid(Model, Grid, iXML)
case("WARPSPH")
Grid%gType = RAI_G_WARPSPH
! Generate our own grid from scratch
!call raijuGenWarpSphGrid_Fok2021(Model, Grid, iXML)
call raijuGenWarpSphGrid_Shafee2008(Model, Grid, iXML)
case("SHGRID")
Grid%gType = RAI_G_SHGRID
! Then we should be receiving a predefined ShellGrid that Voltron has set up
if(present(shGridO)) then
shGrid = shGridO
call raijuGenGridFromShGrid(Grid%shGrid, shGrid, iXML, opts)
else
write(*,*) "RAIJU expecting a ShellGrid_T but didn't receive one. Dying."
endif
case DEFAULT
write(*,*) "RAIJU Received invalid grid definition: ",tmpStr
write(*,*) " Dying."
stop
end select
! Make sure we aren't going beyond the equator with our ghost cells, that wouldn't make sense
if (any(Grid%shGrid%th(Grid%shGrid%ie+1:) > PI/2.0_rp)) then
write(*,*)"ERROR: RAIJU ghosts going beyond 0 deg lat. Probably you are pushing it in too close to planet's surface"
stop
endif
! Finalize the spatial part of the grid
call finalizeLLGrid(Grid, Model%planet)
! Hax?
call iXML%Set_Val(doHack_constB, "hax/doConstB",.false.)
if (doHack_constB) then
Grid%Bmag = Model%planet%magMoment*G2nT
Grid%Brcc = Model%planet%magMoment*G2nT
Grid%BrFace = Model%planet%magMoment*G2nT
endif
! Now handle lambda grid
Grid%nSpc = Model%nSpc ! Make a copy of nSpc
Grid%nFluidIn = Model%nFluidIn ! Make a copy of nFluidIn
call iXML%Set_Val(Grid%ignoreConfigMismatch, "config/ignoreMismatch",.false.)
call initLambdaGrid(Model, Grid, Model%configFName)
end subroutine raijuInitGrid
subroutine raijuInitState(Model, Grid, State, iXML)
type(raijuModel_T), intent(inout) :: Model
type(raijuGrid_T) , intent(in) :: Grid
type(raijuState_T), intent(inout) :: State
type(XML_Input_T), intent(in) :: iXML
integer :: s
! Allocate arrays
associate(sh=>Grid%shGrid)
! dt for every lambda channel
allocate( State%dtk (Grid%Nk) )
call fillArray(State%dtk, 0.0_rp)
! nSteps for each channel
allocate( State%nStepk(Grid%Nk) )
call fillArray(State%nStepk, 0)
! Where we keep all our stuff
allocate( State%eta (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%eta, 0.0_rp)
! Where we keep all our stuff but a half-step ahead of now
allocate( State%eta_half (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%eta_half, 0.0_rp)
! Where we kept all our stuff one step ago
allocate( State%eta_last (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%eta_last, 0.0_rp)
! Where all the stuff sorta was over the last State%dt
allocate( State%eta_avg (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%eta_avg, 0.0_rp)
! I shells shat should be evolved for each k
allocate( State%activeShells (sh%isg:sh%ieg, Grid%Nk) )
State%activeShells = .false.
! Effective potential (used for output only)
allocate( State%pEff(sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk) )
call fillArray(State%pEff, 0.0_rp)
! Gradient of ionspheric potential
allocate( State%gradPotE (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, 2) )
call fillArray(State%gradPotE, 0.0_rp)
! Gradient of corotation potential
allocate( State%gradPotCorot (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, 2) )
call fillArray(State%gradPotCorot, 0.0_rp)
! Gradient of (flux tube volume ^ -2/3)
allocate( State%gradVM (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, 2) )
call fillArray(State%gradVM, 0.0_rp)
! Interface and cell velocities
allocate( State%gradPotE_cc (sh%isg:sh%ieg, sh%jsg:sh%jeg, 2) )
call fillArray(State%gradPotE_cc, 0.0_rp)
allocate( State%gradPotCorot_cc(sh%isg:sh%ieg, sh%jsg:sh%jeg, 2) )
call fillArray(State%gradPotCorot_cc, 0.0_rp)
allocate( State%gradVM_cc (sh%isg:sh%ieg, sh%jsg:sh%jeg, 2) )
call fillArray(State%gradVM_cc, 0.0_rp)
allocate( State%iVel (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
call fillArray(State%iVel, 0.0_rp)
allocate( State%iVelL(sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
call fillArray(State%iVelL, 0.0_rp)
allocate( State%iVelR(sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
call fillArray(State%iVelR, 0.0_rp)
allocate( State%cVel (sh%isg:sh%ieg , sh%jsg:sh%jeg , Grid%Nk, 2) )
call fillArray(State%cVel, 0.0_rp)
! Coupling input moments
allocate( State%Pavg(sh%isg:sh%ieg , sh%jsg:sh%jeg, 0:Grid%nFluidIn) )
call fillArray(State%Pavg, 0.0_rp)
allocate( State%Davg(sh%isg:sh%ieg , sh%jsg:sh%jeg, 0:Grid%nFluidIn) )
call fillArray(State%Davg, 0.0_rp)
allocate( State%Pstd(sh%isg:sh%ieg , sh%jsg:sh%jeg, 0:Grid%nFluidIn) )
call fillArray(State%Pstd, 0.0_rp)
allocate( State%Dstd(sh%isg:sh%ieg , sh%jsg:sh%jeg, 0:Grid%nFluidIn) )
call fillArray(State%Dstd, 0.0_rp)
allocate( State%domWeights(sh%isg:sh%ieg , sh%jsg:sh%jeg) )
call fillArray(State%domWeights, 0.0_rp)
allocate( State%tiote(sh%isg:sh%ieg , sh%jsg:sh%jeg) )
call fillArray(State%tiote, 0.0_rp)
call initShellVar(Grid%shGrid, SHGR_CC, State%Tb)
State%Tb%data = 0.0
State%Tb%mask = .false.
! Bmin surface
allocate( State%Bmin (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, 3 ) )
call fillArray(State%Bmin, 0.0_rp)
allocate( State%xyzMin (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, 3 ) )
call fillArray(State%xyzMin, 0.0_rp)
allocate( State%xyzMincc(sh%isg:sh%ieg , sh%jsg:sh%jeg , 3 ) )
call fillArray(State%xyzMincc, 0.0_rp)
allocate( State%thcon (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1 ) )
call fillArray(State%thcon, 0.0_rp)
allocate( State%phcon (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1 ) )
call fillArray(State%phcon, 0.0_rp)
! 2D corner quantities
allocate( State%topo (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
call fillArray(State%topo, 0)
allocate( State%espot (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
call fillArray(State%espot, 0.0_rp)
allocate( State%pot_corot(sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
call fillArray(State%pot_corot, 0.0_rp)
allocate( State%bvol (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
call fillArray(State%bvol, 0.0_rp)
allocate( State%bvol_cc (sh%isg:sh%ieg , sh%jsg:sh%jeg ) )
call fillArray(State%bvol_cc, 0.0_rp)
allocate( State%vaFrac (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1) )
call fillArray(State%vaFrac, 0.0_rp)
! 1D cell-centered quantities
allocate( State%bndLoc(sh%jsg:sh%jeg) )
call fillArray(State%bndLoc, 0)
! 2D cell-centered quantities
allocate( State%active (sh%isg:sh%ieg, sh%jsg:sh%jeg) )
call fillArray(State%active, 0)
allocate( State%active_last (sh%isg:sh%ieg, sh%jsg:sh%jeg) )
call fillArray(State%active_last, 0)
allocate( State%OCBDist(sh%isg:sh%ieg, sh%jsg:sh%jeg) )
call fillArray(State%OCBDist, 0)
allocate( State%lossRates (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%lossRates, 0.0_rp)
allocate( State%precipType_ele (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%precipType_ele, 0.0_rp)
allocate( State%lossRatesPrecip(sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%lossRatesPrecip, 0.0_rp)
!allocate( State%precipNFlux (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
!allocate( State%precipEFlux (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
allocate( State%dEta_dt (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%dEta_dt, 0.0_rp)
allocate( State%CCHeatFlux (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
call fillArray(State%CCHeatFlux, 0.0_rp)
! Coupling output data
allocate(State%Den (0:Model%nSpc))
allocate(State%Press(0:Model%nSpc))
allocate(State%Den_avg (0:Model%nSpc))
allocate(State%Press_avg(0:Model%nSpc))
do s=0,Model%nSpc
call initShellVar(Grid%shGrid, SHGR_CC, State%Den (s))
call initShellVar(Grid%shGrid, SHGR_CC, State%Press(s))
State%Den (s)%data = 0.0
State%Press(s)%data = 0.0
State%Den (s)%mask = .false.
State%Press(s)%mask = .false.
call initShellVar(Grid%shGrid, SHGR_CC, State%Den_avg (s))
call initShellVar(Grid%shGrid, SHGR_CC, State%Press_avg(s))
State%Den_avg (s)%data = 0.0
State%Press_avg(s)%data = 0.0
State%Den_avg (s)%mask = .false.
State%Press_avg(s)%mask = .false.
enddo
allocate(State%precipNFlux(Grid%Nk))
allocate(State%precipEFlux(Grid%Nk))
do s=1,Grid%Nk
call initShellVar(Grid%shGrid, SHGR_CC, State%precipNFlux(s))
call initShellVar(Grid%shGrid, SHGR_CC, State%precipEFlux(s))
State%precipNFlux(s)%data = 0.0
State%precipEFlux(s)%data = 0.0
State%precipNFlux(s)%mask = .false.
State%precipEFlux(s)%mask = .false.
enddo
! Only bother allocating persistent versions of debug stuff if we need them
if (Model%doOutput_debug) then
allocate( State%etaFaceReconL(sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
call fillArray(State%etaFaceReconL, 0.0_rp)
allocate( State%etaFaceReconR(sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
call fillArray(State%etaFaceReconR, 0.0_rp)
allocate( State%etaFacePDML (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
call fillArray(State%etaFacePDML, 0.0_rp)
allocate( State%etaFacePDMR (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
call fillArray(State%etaFacePDMR, 0.0_rp)
allocate( State%etaFlux (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
call fillArray(State%etaFlux, 0.0_rp)
endif
State%KpTS%wID = Model%tsF
call State%KpTS%initTS("Kp")
end associate
! For now, just set t to tStart and ts to 0
State%t = Model%t0
State%ts = 0
State%nStepk = 0
! Set all activeShells to true. There should be doActiveShell checks anywhere it is used, but just in case, make sure default has all shells on
State%activeShells = .true.
! Similarly, set vaFrac to safe value in case stand-alone never writes to it
State%vaFrac = 1.0
State%isFirstCpl = .true.
! Init State sub-modules
if (Model%isSA) then
! If we are standalone, this is the place to get coldStarter settings
! Otherwise, raiCpl should call it so we can use other models' info
call initRaijuColdStarter(Model, iXML, State%coldStarter)
endif
! Init state
call iXML%Set_Val(Model%icStr, "prob/IC","DIP")
select case(trim(Model%icStr))
case("DIP")
!! Simple Maxwellian in a dipole field
!Model%initState => initRaijuIC_DIP
call initRaijuIC_DIP(Model, Grid, State, iXML)
case("TTA")
!! Test Theta Advection
!Model%initState => initRaijuIC_testThetaAdvection
call initRaijuIC_testThetaAdvection(Model, Grid, State, iXML)
case("USER")
! Call the IC in the module raijuuseric
! This module is set in cmake via the RAIJUIC variable
!Model%initState => userInitStateFunc
call raijuInitState_useric(Model, Grid, State, iXML)
case DEFAULT
write(*,*)"Invalid IC name to RAIJU, see raijuStarter.F90:raijuInitState. Bye."
stop
end select
!call Model%initState(Model, Grid, State, iXML)
end subroutine raijuInitState
end module raijustarter
Reference in New Issue View Git Blame Copy Permalink