mirror of
https://github.com/JHUAPL/kaiju.git
synced 2026-01-08 22:58:05 -05:00
Merge remote-tracking branch 'origin/development' into ewinter-derecho_testing
This commit is contained in:
Eric Winter
@@ -846,7 +846,7 @@ def run_preprocessing_steps(options: dict, args: dict):
|
||||
|
||||
# If needed, create the solar wind file by fetching data from CDAWeb.
|
||||
# NOTE: Assumes kaipy is installed.
|
||||
if options["simulation"]["bcwind_available"] == "N":
|
||||
if options["simulation"]["bcwind_available"].upper() == "N":
|
||||
cmd = "cda2wind"
|
||||
args = [cmd, "-t0", options["simulation"]["start_date"], "-t1",
|
||||
options["simulation"]["stop_date"], "-interp", "-bx",
|
||||
|
||||
@@ -72,10 +72,14 @@ module raijudefs
|
||||
logical, parameter :: def_doUseVelLRs = .true.
|
||||
|
||||
! Domain limits
|
||||
! Buffer not allowed beyond min of maxTail and maxSun
|
||||
real(rp), parameter :: def_maxTail_buffer = 15.0 ! [Rp]
|
||||
real(rp), parameter :: def_maxSun_buffer = 10.0 ! [Rp]
|
||||
! Active not allowed beyond min of maxTail and maxSun
|
||||
real(rp), parameter :: def_maxTail_active = 10.0 ! [Rp]
|
||||
real(rp), parameter :: def_maxSun_active = 10.0 ! [Rp]
|
||||
! Active is forced below activeDomRad
|
||||
real(rp), parameter :: def_activeDomRad = 3.0 ! [Rp]
|
||||
|
||||
! Settings
|
||||
integer, parameter :: raiRecLen = 8
|
||||
@@ -97,14 +101,19 @@ module raijudefs
|
||||
!! [deg] Max allowable angle between any two normals surrounding a cell corner
|
||||
real(rp), parameter :: def_bminThresh = 1.0
|
||||
!! [nT] default allowable bmin strencgh
|
||||
real(rp), parameter :: def_nBounce = 1.0
|
||||
real(rp), parameter :: def_nBounce = 3.0
|
||||
!! Number of Alfven bounces (Tb) required to be considered a "good" flux tube for coupling
|
||||
|
||||
|
||||
real(rp), parameter :: def_lim_vaFrac_soft = 0.6_rp
|
||||
real(rp), parameter :: def_lim_vaFrac_hard = 0.4_rp
|
||||
real(rp), parameter :: def_lim_bmin_soft = 5.0_rp
|
||||
! [nT]
|
||||
!! [nT]
|
||||
real(rp), parameter :: def_lim_bmin_hard = 2.0_rp
|
||||
!! [nT]
|
||||
|
||||
! Plasmasphere stuff
|
||||
real(rp), parameter :: def_psphEvolRad = 2.25
|
||||
!! [Rp] radius under which plasmasphere is not evolved
|
||||
|
||||
end module raijudefs
|
||||
@@ -245,6 +245,8 @@ module raijutypes
|
||||
real(rp) :: psphInitKp
|
||||
logical :: doPsphEvol
|
||||
!! Whether or not to actually evolve the plasmasphere
|
||||
real(rp) :: psphEvolRad
|
||||
!! [Rp] Radius below which plasmasphere is not evolved
|
||||
! TODO: Extra params for refilling rate, determining initial profile, etc.
|
||||
|
||||
! Some constants
|
||||
@@ -261,6 +263,8 @@ module raijutypes
|
||||
!! Maximum tailward extent of the active region
|
||||
real(rp) :: maxSun_active
|
||||
!! Maximum sunward extent of the active region
|
||||
real(rp) :: activeDomRad
|
||||
!! [Rp] Cells are forced to be active below this radius
|
||||
|
||||
! Active shell settings
|
||||
logical :: doActiveShell
|
||||
@@ -390,6 +394,10 @@ module raijutypes
|
||||
type(IOClock_T) :: IO
|
||||
!! Timers for IO operations
|
||||
|
||||
! I feel like philosophically this should be in Grid but that feels weird so its here
|
||||
type(TimeSeries_T) :: KpTS
|
||||
!! Kp timeseries from solar wind file
|
||||
|
||||
! -- Solver values -- !
|
||||
real(rp), dimension(:,:,:), allocatable :: eta
|
||||
!! (Ngi, Ngj, Nk) [#/cc * Rp/T] etas
|
||||
|
||||
@@ -175,6 +175,10 @@ module volttypes
|
||||
!! [s] Time scale over which we ramp up to full IM_TSCL for MHD ingestion
|
||||
logical :: doColdstartCX = .true.
|
||||
!! Whether or not we apply charge exchange to initial coldstart ion profile
|
||||
real(rp) :: tsclSm_dL = 1.0_rp
|
||||
!! [Rp] theta/L-direction smoothing scale for mhd ingestion timescale
|
||||
real(rp) :: tsclSm_dMLT = 1.0_rp
|
||||
!! [hr] phi/MLT-direction smoothing scale for mhd ingestion timescale
|
||||
|
||||
! --- Grid --- !
|
||||
type(ShellGrid_T) :: shGr
|
||||
@@ -219,13 +223,18 @@ module volttypes
|
||||
!! Used to limit active region to tubes that can reasonably be treated as averaged and slowly-evolving
|
||||
type(ShellGridVar_T) :: Tb
|
||||
!! (Ngi, Ngj) Bounce timesale
|
||||
|
||||
type(ShellGridVar_T) :: avgBeta
|
||||
!! (Ngi, Ngj) Cross-sectional flux-tube area - averaged beta (woof)
|
||||
|
||||
type(ShellGridVar_T) :: pot_total
|
||||
!! Total electrostatic potential from (ionosphere + corot) [kV]
|
||||
type(ShellGridVar_T) :: pot_corot
|
||||
!! Just corotation potential, just for output purposes [kV]
|
||||
|
||||
! Post-advance stuff leaving raiju
|
||||
type(ShellGridVar_T) :: tscl_mhdIngest
|
||||
!! [s] Timescale at which we think MHD should ingest our info
|
||||
|
||||
contains
|
||||
|
||||
procedure :: toIMAG => volt2RAIJU
|
||||
|
||||
@@ -90,11 +90,6 @@ module streamline
|
||||
doShue = .false.
|
||||
endif
|
||||
|
||||
if ( (.not. doSH) .and. (.not. doNH) ) then
|
||||
!What are you even asking for?
|
||||
return
|
||||
endif
|
||||
|
||||
!Allocate temp arrays to hold information along each direction
|
||||
!Trailing dimension (1-2) represents -/+ directions
|
||||
allocate(xyzn (0:MaxFL,NDIM,2))
|
||||
|
||||
@@ -111,11 +111,16 @@ module raijuAdvancer
|
||||
odt = State%dtk(k)
|
||||
|
||||
! NOTE: When we actually implement activeShells these will need to be updated within time loop
|
||||
where (State%active .eq. RAIJUACTIVE)
|
||||
isGoodEvol = .true.
|
||||
elsewhere
|
||||
isGoodEvol = .false.
|
||||
end where
|
||||
if (Grid%spc(Grid%k2spc(k))%flav==F_PSPH) then
|
||||
call calcEvolDom_psph(Model, Grid, State, isGoodEvol)
|
||||
else
|
||||
where (State%active .eq. RAIJUACTIVE)
|
||||
isGoodEvol = .true.
|
||||
elsewhere
|
||||
isGoodEvol = .false.
|
||||
end where
|
||||
endif
|
||||
|
||||
where (State%active .ne. RAIJUINACTIVE)
|
||||
isGoodRecon = .true.
|
||||
elsewhere
|
||||
@@ -284,4 +289,33 @@ module raijuAdvancer
|
||||
|
||||
end subroutine calcEtaHalf
|
||||
|
||||
|
||||
subroutine calcEvolDom_psph(Model, Grid, State, isGoodEvol)
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T ), intent(in) :: Grid
|
||||
type(raijuState_T), intent(in) :: State
|
||||
|
||||
logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg, &
|
||||
Grid%shGrid%jsg:Grid%shGrid%jeg), intent(inout) :: isGoodEvol
|
||||
|
||||
integer :: i,j
|
||||
real(rp) :: Req
|
||||
|
||||
isGoodEvol = .false. ! Can't use fillArray yet, no version for logicals
|
||||
|
||||
do j=Grid%shGrid%jsg,Grid%shGrid%jeg
|
||||
do i=Grid%shGrid%isg,Grid%shGrid%ieg
|
||||
Req = sqrt(State%xyzMincc(i,j,XDIR)**2 + State%xyzMincc(i,j,YDIR)**2)
|
||||
if (Req < Model%psphEvolRad) then
|
||||
continue ! Keep points within this radius set to no evol
|
||||
else
|
||||
if (State%active(i,j) == RAIJUACTIVE) then
|
||||
isGoodEvol(i,j) = .true.
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine
|
||||
|
||||
end module raijuAdvancer
|
||||
|
||||
@@ -34,6 +34,10 @@ module raijuDomain
|
||||
State%active = RAIJUACTIVE
|
||||
end where
|
||||
|
||||
! HOWEVER...
|
||||
! Force all cells below X Re to be active, with enough buffer cells, and suffer the consequences if that includes bad flux tubes
|
||||
call forceActiveWithinRadius(Model, Grid, State, State%active)
|
||||
|
||||
end subroutine setActiveDomain
|
||||
|
||||
|
||||
@@ -572,4 +576,33 @@ module raijuDomain
|
||||
|
||||
end subroutine calcCornerNormAngles
|
||||
|
||||
|
||||
subroutine forceActiveWithinRadius(Model, Grid, State, active)
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T ), intent(in) :: Grid
|
||||
type(raijuState_T), intent(in) :: State
|
||||
integer, dimension(Grid%shGrid%isg:Grid%shGrid%ieg,Grid%shGrid%jsg:Grid%shGrid%jeg), intent(inout) :: active
|
||||
|
||||
integer :: i,j
|
||||
real(rp) :: rad_next
|
||||
|
||||
do j=Grid%shGrid%jsg,Grid%shGrid%jeg
|
||||
i = Grid%shGrid%ie
|
||||
rad_next = sqrt(State%xyzMincc(i-1,j,XDIR)**2 + State%xyzMincc(i-1,j,YDIR)**2)
|
||||
do while(rad_next < Model%activeDomRad)
|
||||
i = i - 1
|
||||
rad_next = sqrt(State%xyzMincc(i-1,j,XDIR)**2 + State%xyzMincc(i-1,j,YDIR)**2)
|
||||
enddo
|
||||
! i is now the last cell in this j within activeDomRad
|
||||
if (State%active(i,j) .ne. RAIJUACTIVE) then
|
||||
! Something bad happened, active domain is really small
|
||||
! Throw caution to the wind and make it active anyways, good luck everybody
|
||||
active(i :Grid%shGrid%ieg ,j) = RAIJUACTIVE
|
||||
active(i-Grid%nB-1 :i-1 ,j) = RAIJUBUFFER
|
||||
active(Grid%shGrid%isg:i-Grid%nB-2 ,j) = RAIJUINACTIVE
|
||||
endif
|
||||
enddo
|
||||
|
||||
end subroutine forceActiveWithinRadius
|
||||
|
||||
end module raijuDomain
|
||||
@@ -420,4 +420,61 @@ module raijuetautils
|
||||
|
||||
end function getInitPsphere
|
||||
|
||||
!Do plasmasphere refilling for the interval we're about to advance
|
||||
subroutine plasmasphereRefill(Model,Grid,State)
|
||||
type(raijuModel_T), intent(in) :: Model
|
||||
type(raijuGrid_T) , intent(in) :: Grid
|
||||
type(raijuState_T), intent(inout) :: State
|
||||
|
||||
integer :: i,j,k0
|
||||
real(rp), parameter :: maxX = 2.0 !Max over-filling relative to target, i.e. don't go above maxX x den-target
|
||||
real(rp), parameter :: day2s = 24.0*60.0*60,s2day=1.0/day2s
|
||||
real(rp) :: NowKp
|
||||
real(rp) :: xeq,yeq,rad,cc2eta,eta2cc
|
||||
real(rp) :: dppT,dpsph,etaT,dndt,deta
|
||||
|
||||
|
||||
k0 = Grid%spc(spcIdx(Grid, F_PSPH))%kStart !plasmasphere index
|
||||
NowKp = State%KpTS%evalAt(State%t)
|
||||
|
||||
!$OMP PARALLEL DO default(shared) &
|
||||
!$OMP private(i,j,xeq,yeq,rad,dppT,cc2eta,eta2cc) &
|
||||
!$OMP private(dpsph,etaT,deta,dndt)
|
||||
do j=Grid%shGrid%jsg,Grid%shGrid%jeg
|
||||
do i=Grid%shGrid%isg,Grid%shGrid%ieg
|
||||
if (State%active(i,j) == RAIJUACTIVE) then
|
||||
|
||||
xeq = State%xyzMincc(i,j,XDIR)
|
||||
yeq = State%xyzMincc(i,j,YDIR)
|
||||
rad = sqrt(xeq**2.0 + yeq**2.0)
|
||||
|
||||
!Closed field line, calculate Gallagher w/ current Kp to get target density
|
||||
dppT = GallagherXY(xeq,yeq,NowKp)
|
||||
cc2eta = State%bvol_cc(i,j)*sclEta
|
||||
eta2cc = 1.0/cc2eta !Convert eta to #/cc
|
||||
|
||||
dpsph = eta2cc*State%eta(i,j,k0) !Current plasmasphere density [#/cc]
|
||||
|
||||
!Check for other outs before doing anything
|
||||
if (dpsph >= maxX*dppT) then
|
||||
continue ! Too much already there, don't touch it
|
||||
else
|
||||
etaT = dppT/eta2cc
|
||||
|
||||
if ((rad <= Model%psphEvolRad) .and. (dppT > dpsph)) then
|
||||
!If this is inside MHD inner boundary, be at least at target value
|
||||
State%eta(i,j,k0) = etaT
|
||||
else
|
||||
!If still here then calculate refilling
|
||||
dndt = 10.0**(3.48-0.331*rad) !cm^-3/day, Denton+ 2012 eqn 1
|
||||
deta = (State%dt*s2day)*dndt/eta2cc !Change in eta over dt
|
||||
State%eta(i,j,k0) = State%eta(i,j,k0) + deta
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
enddo !i
|
||||
enddo !j
|
||||
|
||||
end subroutine plasmasphereRefill
|
||||
|
||||
end module raijuetautils
|
||||
|
||||
@@ -78,6 +78,12 @@ module raijuOut
|
||||
State%IO%tRes = State%IO%tRes + State%IO%dtRes
|
||||
State%IO%nRes = State%IO%nRes + 1
|
||||
|
||||
! Handle case where we are coupled to voltron and 5s ahead of everyone
|
||||
! 5s will accumulate in restart file even though we only ever stay 5s ahead
|
||||
if (.not. Model%isSA .and. State%t > State%dt) then
|
||||
State%t = State%t - State%dt ! Current time - coupling time
|
||||
endif
|
||||
|
||||
|
||||
end subroutine raijuResInput
|
||||
|
||||
|
||||
@@ -41,6 +41,9 @@ module raijuPreAdvancer
|
||||
call applyRaijuBCs(Model, Grid, State, doWholeDomainO=State%isFirstCpl) ! If fullEtaMap=True, mom2eta map is applied to the whole domain
|
||||
call Toc("BCs")
|
||||
|
||||
! Handle plasmaasphere refilling for the full step about to happen
|
||||
call plasmasphereRefill(Model,Grid,State)
|
||||
|
||||
! Handle edge cases that may effect the validity of information carried over from last coupling period
|
||||
call prepEtaLast(Grid%shGrid, State, State%isFirstCpl)
|
||||
|
||||
|
||||
@@ -143,6 +143,7 @@ module raijustarter
|
||||
!--- 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
|
||||
@@ -164,11 +165,13 @@ module raijustarter
|
||||
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%doUseVelLRs,'sim/useVelLRs',def_doUseVelLRs)
|
||||
@@ -483,6 +486,9 @@ module raijustarter
|
||||
allocate( State%etaFacePDMR (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
|
||||
allocate( State%etaFlux (sh%isg:sh%ieg+1, sh%jsg:sh%jeg+1, Grid%Nk, 2) )
|
||||
endif
|
||||
|
||||
State%KpTS%wID = Model%tsF
|
||||
call State%KpTS%initTS("Kp")
|
||||
|
||||
end associate
|
||||
|
||||
|
||||
@@ -93,7 +93,7 @@ module raijulosses
|
||||
integer :: k
|
||||
logical, dimension(Grid%shGrid%isg:Grid%shGrid%ieg,Grid%shGrid%jsg:Grid%shGrid%jeg) :: isGood
|
||||
|
||||
where (State%active .eq. RAIJUACTIVE)
|
||||
where (State%active .ne. RAIJUINACTIVE)
|
||||
isGood = .true.
|
||||
elsewhere
|
||||
isGood = .false.
|
||||
|
||||
@@ -167,7 +167,6 @@ module mixmain
|
||||
call dragonking_total(I(h)%aurora,I(h)%G,I(h)%St,I(h)%conductance)
|
||||
|
||||
if (present(gcm)) then
|
||||
write(*,*) 'doGCM!'
|
||||
call conductance_total(I(h)%conductance,I(h)%G,I(h)%St,gcm,h)
|
||||
else
|
||||
!write(*,*) "conductance: total"
|
||||
|
||||
@@ -144,7 +144,7 @@ module mixparams
|
||||
call xmlInp%Set_Val(Params%apply_cap,"conductance/apply_cap",.true.)
|
||||
call xmlInp%Set_Val(Params%gamma,"precipitation/gamma",5.0/3.0)
|
||||
call xmlInp%Set_Val(Params%doEMA,"conductance/doEMA",.true.)
|
||||
call xmlInp%Set_Val(Params%doET,"conductance/doET",.true.)
|
||||
call xmlInp%Set_Val(Params%doET,"conductance/doET",.false.)
|
||||
|
||||
! if this is on, Params%F107 read from xml file above will be overwritten
|
||||
! and F107 will be set from the SW file (i.e., from OMNI data)
|
||||
|
||||
@@ -1,257 +0,0 @@
|
||||
|
||||
!Various routines to handle flux-tube calculations to feed inner magnetosphere models
|
||||
module imagtubes
|
||||
use volttypes
|
||||
use voltcpltypes
|
||||
use streamline
|
||||
use chmpdbz, ONLY : DipoleB0
|
||||
use imaghelper
|
||||
use planethelper
|
||||
use shellGrid
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
!subroutine init_TubeShell(sh, tubeShell, maskO, TioTeO)
|
||||
subroutine init_TubeShell(sh, tubeShell, TioTeO)
|
||||
type(ShellGrid_T), intent(in) :: sh
|
||||
!! Voltron shellgrid
|
||||
type(TubeShell_T), intent(inout) :: tubeShell
|
||||
!! IMAGTubeShell object we are initializing
|
||||
!logical, dimension(:,:), intent(in), optional :: maskO
|
||||
real(rp), intent(in), optional :: TioteO
|
||||
!! Default Ti/Te ratio
|
||||
|
||||
integer :: i
|
||||
real(rp) :: tiote
|
||||
|
||||
if (present(TioteO)) then
|
||||
tiote = TioteO
|
||||
else
|
||||
tiote = 4.0
|
||||
endif
|
||||
|
||||
|
||||
do i=1,NDIM
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%xyz0(i))
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%X_bmin(i))
|
||||
enddo
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%lat0 )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%lon0 )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%latc )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%lonc )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%invlat )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%topo )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%bmin )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%bVol )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%Lb )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%Tb )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%wMAG )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%rCurv )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%avgBeta)
|
||||
do i=0,MAXTUBEFLUIDS
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%avgP(i))
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%avgN(i))
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%stdP(i))
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%stdN(i))
|
||||
enddo
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%losscone )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%lossconec)
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%TioTe0 )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%nTrc )
|
||||
|
||||
tubeShell%TioTe0%data = tiote
|
||||
|
||||
|
||||
end subroutine init_TubeShell
|
||||
|
||||
! Dipole flux tube info
|
||||
subroutine DipoleTube(vApp,lat,lon,ijTube)
|
||||
type(voltApp_T), intent(in) :: vApp
|
||||
real(rp), intent(in) :: lat,lon
|
||||
type(IMAGTube_T), intent(out) :: ijTube
|
||||
|
||||
real(rp) :: L,colat
|
||||
real(rp) :: mdipole
|
||||
|
||||
mdipole = ABS(vapp%planet%magMoment)*G2T ! dipole moment in T
|
||||
colat = PI/2 - lat
|
||||
L = 1.0/(sin(colat)**2.0)
|
||||
!ijTube%Vol = 32./35.*L**4.0/mdipole
|
||||
!Use full dipole FTV formula, convert Rx/nT => Rx/T
|
||||
ijTube%Vol = DipFTV_L(L,vapp%planet%magMoment)*1.0e+9
|
||||
|
||||
ijTube%X_bmin(XDIR) = L*cos(lon)*vApp%planet%rp_m !Rp=>meters
|
||||
ijTube%X_bmin(YDIR) = L*sin(lon)*vApp%planet%rp_m !Rp=>meters
|
||||
ijTube%X_bmin(ZDIR) = 0.0
|
||||
ijTube%bmin = mdipole/L**3.0
|
||||
|
||||
ijTube%topo = 2
|
||||
|
||||
ijTube%beta_average = 0.0
|
||||
ijTube%Pave = 0.0
|
||||
ijTube%Nave = 0.0
|
||||
|
||||
ijTube%latc = -lat
|
||||
ijTube%lonc = lon
|
||||
ijTube%Lb = L !Just lazily using L shell
|
||||
ijTube%Tb = 0.0
|
||||
ijTube%losscone = 0.0
|
||||
ijTube%rCurv = L/3.0
|
||||
|
||||
ijTube%wIMAG = 1.0 !Much imag
|
||||
ijTube%TioTe0 = def_tiote
|
||||
end subroutine DipoleTube
|
||||
|
||||
|
||||
subroutine MHDTube(ebTrcApp,planet,colat,lon,r,ijTube,bTrc,nTrcO,doShiftO)
|
||||
type(ebTrcApp_T), intent(in) :: ebTrcApp
|
||||
type(planet_T), intent(in) :: planet
|
||||
real(rp), intent(in) :: colat, lon, r
|
||||
type(IMAGTube_T), intent(inout) :: ijTube
|
||||
type(magLine_T), intent(inout) :: bTrc
|
||||
integer, intent(in), optional :: nTrcO
|
||||
logical, intent(in), optional :: doShiftO
|
||||
|
||||
integer :: s
|
||||
real(rp) :: t, bMin,bIon
|
||||
real(rp), dimension(NDIM) :: x0, bEq, xyzIon
|
||||
real(rp), dimension(NDIM) :: xyzC,xyzIonC
|
||||
integer :: OCb
|
||||
real(rp) :: bD,bP,dvB,bBeta,rCurv, bDRC, bPRC, bBetaRC, stdP, stdD
|
||||
real(rp) :: VaMKS,CsMKS,VebMKS !Speeds in km/s
|
||||
real(rp) :: TiEV !Temperature in ev
|
||||
|
||||
logical :: doShift,doShue
|
||||
|
||||
if (present(doShiftO)) then
|
||||
doShift = doShiftO
|
||||
else
|
||||
doShift = .false.
|
||||
endif
|
||||
|
||||
! Take seed point in spherical coordinates
|
||||
xyzIon(XDIR) = r*sin(colat)*cos(lon)
|
||||
xyzIon(YDIR) = r*sin(colat)*sin(lon)
|
||||
xyzIon(ZDIR) = r*cos(colat)
|
||||
if (doShift) then
|
||||
!! TODO: Decide if this is how we want to do it
|
||||
x0 = DipoleShift(xyzIon, planet%ri_m/planet%rp_m + TINY)
|
||||
else
|
||||
x0 = xyzIon
|
||||
endif
|
||||
|
||||
bIon = norm2(DipoleB0(xyzIon))*oBScl*1.0e-9 !EB=>T, ionospheric field strength
|
||||
|
||||
|
||||
associate(ebModel=>ebTrcApp%ebModel,ebGr=>ebTrcApp%ebState%ebGr,ebState=>ebTrcApp%ebState)
|
||||
|
||||
!Now do field line trace
|
||||
t = ebState%eb1%time !Time in CHIMP units
|
||||
!!TODO: What do we do when we want times in between steps? Somethign similar to what slice/chop do to output
|
||||
|
||||
if (present(nTrcO)) then
|
||||
call genLine(ebModel,ebState,x0,t,bTrc,nTrcO,doShueO=.false.,doNHO=.true.)
|
||||
else
|
||||
call genLine(ebModel,ebState,x0,t,bTrc, doShueO=.false.,doNHO=.true.)
|
||||
endif
|
||||
|
||||
|
||||
!Topology
|
||||
!OCB = 0 (solar wind), 1 (half-closed), 2 (both ends closed)
|
||||
OCb = FLTop(ebModel,ebGr,bTrc)
|
||||
ijTube%topo = OCb
|
||||
|
||||
if (OCb /= 2) then
|
||||
! If the field line hasn't closed after 15 minutes we are legally allowed to leave
|
||||
ijTube%X_bmin = 0.0
|
||||
ijTube%bmin = 0.0
|
||||
ijTube%Vol = -1.0
|
||||
ijTube%Pave = 0.0
|
||||
ijTube%Nave = 0.0
|
||||
ijTube%Pstd = 0.0
|
||||
ijTube%Nstd = 0.0
|
||||
ijTube%beta_average = 0.0
|
||||
ijTube%latc = 0.0
|
||||
ijTube%lonc = 0.0
|
||||
ijTube%Lb = 0.0
|
||||
ijTube%Tb = 0.0
|
||||
ijTube%losscone = 0.0
|
||||
ijTube%rCurv = 0.0
|
||||
ijTube%wIMAG = 0.0
|
||||
ijTube%Veb = 0.0
|
||||
return
|
||||
endif
|
||||
|
||||
!Get diagnostics from closed field line
|
||||
!Minimal surface (bEq in Rp, bMin in EB)
|
||||
call FLEq(ebModel,bTrc,bEq,bMin)
|
||||
bMin = bMin*oBScl*1.0e-9 !EB=>Tesla
|
||||
bEq = bEq*planet%rp_m ! Rp -> meters
|
||||
|
||||
! Plasma quantities
|
||||
do s=0,ebModel%nSpc
|
||||
!dvB = Flux-tube volume (Re/EB)
|
||||
!write(*,*)"FLThermo, s=",s
|
||||
call FLThermo(ebModel,ebGr,bTrc,bD,bP,dvB,bBeta,s)
|
||||
call FLStdev (ebModel,ebGr,bTrc,bD,bP,stdD,stdP,s)
|
||||
!call FLStdev (ebModel,ebGr,bTrc,stdD,stdP,s)
|
||||
bP = bP*1.0e-9 !nPa=>Pa
|
||||
bD = bD*1.0e+6 !#/cc => #/m3
|
||||
ijTube%Pave(s) = bP
|
||||
ijTube%Nave(s) = bD
|
||||
ijTube%Pstd(s) = stdP * 1.0e-9 ! nPa=>Pa
|
||||
ijTube%Nstd(s) = stdD * 1.0e+6 ! #/cc=>#/m3
|
||||
if (s .eq. RCFLUID) then
|
||||
bPRC = bP
|
||||
bDRC = bD
|
||||
bBetaRC = bBeta
|
||||
endif
|
||||
enddo
|
||||
|
||||
!Converts Rp/EB => Rp/T
|
||||
dvB = dvB/(oBScl*1.0e-9)
|
||||
|
||||
ijTube%X_bmin = bEq
|
||||
ijTube%bmin = bMin
|
||||
ijTube%Vol = dvB ! Taken from last FLThermo call
|
||||
ijTube%beta_average = bBetaRC
|
||||
|
||||
call FLConj(ebModel,ebGr,bTrc,xyzC)
|
||||
if (doShift) then
|
||||
xyzIonC = DipoleShift(xyzC, planet%ri_m)
|
||||
else
|
||||
xyzIonC = xyzC
|
||||
endif
|
||||
ijTube%latc = asin(xyzIonC(ZDIR)/norm2(xyzIonC))
|
||||
ijTube%lonc = modulo( atan2(xyzIonC(YDIR),xyzIonC(XDIR)),2*PI )
|
||||
ijTube%Lb = FLArc(ebModel,ebGr,bTrc)
|
||||
!NOTE: Bounce timescale may be altered to use IMAG hot density
|
||||
ijTube%Tb = FLAlfvenX(ebModel,ebGr,bTrc)
|
||||
!ijTube%Tb = FLFastX(ebModel,ebGr,bTrc)
|
||||
|
||||
!write(*,*) 'Bounce compare: = ', FLFastX(ebModel,ebGr,bTrc)/FLAlfvenX(ebModel,ebGr,bTrc)
|
||||
|
||||
ijTube%losscone = asin(sqrt(bMin/bIon))
|
||||
|
||||
!Get curvature radius and ExB velocity [km/s]
|
||||
call FLCurvRadius(ebModel,ebGr,ebState,bTrc,rCurv,VebMKS)
|
||||
ijTube%rCurv = rCurv
|
||||
ijTube%Veb = VebMKS
|
||||
|
||||
!Get confidence interval
|
||||
!VaMKS = flux tube arc length [km] / Alfven crossing time [s]
|
||||
VaMKS = (ijTube%Lb*planet%rp_m*1.0e-3)/ijTube%Tb
|
||||
!CsMKS = 9.79 x sqrt(5/3 * Ti) km/s, Ti eV
|
||||
!TiEV = (1.0e+3)*DP2kT(bDRC*1.0e-6,bPRC*1.0e+9) !Temp in eV
|
||||
TiEV = (1.0e+3)*DP2kT(ijTube%Nave(BLK)*1.0D-6,ijTube%Pave(BLK)*1.0D9) !Temp in eV
|
||||
CsMKS = 9.79*sqrt((5.0/3)*TiEV)
|
||||
|
||||
ijTube%wIMAG = VaMKS/( sqrt(VaMKS**2.0 + CsMKS**2.0) + VebMKS)
|
||||
|
||||
end associate
|
||||
|
||||
end subroutine MHDTube
|
||||
|
||||
end module imagtubes
|
||||
@@ -134,6 +134,10 @@ submodule (volttypes) raijuCplTypesSub
|
||||
real(rp) :: active_interp
|
||||
real(rp) :: d_cold, t_cold, d_hot, p_hot
|
||||
real(rp) :: tScl, rampC
|
||||
real(rp) :: wMHD,wRAI,alpha,beta,p_mhd
|
||||
|
||||
logical, parameter :: doWolf = .true.
|
||||
real(rp), parameter :: wRAI0 = 0.5
|
||||
|
||||
associate(Model=>App%raiApp%Model, State=>App%raiApp%State, sh=>App%raiApp%Grid%shGrid, spcList=>App%raiApp%Grid%spc)
|
||||
|
||||
@@ -179,22 +183,32 @@ submodule (volttypes) raijuCplTypesSub
|
||||
endif
|
||||
enddo
|
||||
|
||||
|
||||
call InterpShellVar_TSC_pnt(sh, State%Tb, th, ph, tScl)
|
||||
tScl = Model%nBounce*tScl ! [s]
|
||||
!tScl = 10.0_rp ! [s]
|
||||
!Wolf-limit the pressure going back
|
||||
if (doWolf) then
|
||||
!You're sending the wolf? Shit, that's all you had to say
|
||||
call InterpShellVar_TSC_pnt(App%shGr, App%avgBeta , th, ph, beta )
|
||||
call InterpShellVar_TSC_pnt(App%shGr, App%Pavg(BLK), th, ph, p_mhd) !MHD flux-tube averaged pressure (sum over all fluids)
|
||||
alpha = 1.0 + beta*5.0/6.0
|
||||
wRAI = 1.0/alpha !Full wolf-limited weight
|
||||
wRAI = max(wRAI,wRAI0) !Don't allow Raiju contribution to drop below wRCM0
|
||||
call ClampValue(wRAI,0.0_rp,1.0_rp)
|
||||
wMHD = 1 - wRAI ! = (alpha-1)/alpha w/o any clamping
|
||||
imW(IM_P_RING) = imW(IM_P_RING)*wRAI + wMHD*p_mhd
|
||||
endif
|
||||
|
||||
!call InterpShellVar_TSC_pnt(sh, State%Tb, th, ph, tScl)
|
||||
!tScl = Model%nBounce*tScl ! [s]
|
||||
|
||||
! 1/(x)^4 for x from 1 to 0.5 goes from 1 to 16. Higher exponents means stronger ramp-up
|
||||
!tScl = 15.0_rp/(App%vaFrac%data(i0,j0))**4 ! [s]
|
||||
!tScl = 15.0_rp/(App%vaFrac%data(i0,j0))**2 ! [s]
|
||||
|
||||
call InterpShellVar_TSC_pnt(sh, App%tscl_mhdIngest, th, ph, tScl)
|
||||
|
||||
! Adjust IM_TSCL if we wanna ramp up over time
|
||||
if (t < App%startup_blendTscl) then
|
||||
rampC = RampDown(t, 0.0_rp, App%startup_blendTscl)
|
||||
!tScl = sqrt(tScl*App%startup_blendTscl)*rampC + (1-rampC)*tScl ! idk
|
||||
tScl = rampC*30.0_rp*tScl + (1-rampC)*tScl ! No good reason for 30 except for wanting starting tScl to be ~8-10 minutes
|
||||
!if (th > 50*deg2rad .and. th < 55*deg2rad .and. ph > 130*deg2rad .and. ph < 150*deg2rad) then
|
||||
! write(*,*)"--",t,App%startup_blendTscl,rampC,tScl
|
||||
!endif
|
||||
endif
|
||||
|
||||
imW(IM_TSCL) = tScl
|
||||
@@ -373,6 +387,8 @@ submodule (volttypes) raijuCplTypesSub
|
||||
real(rp), intent(in) :: dt
|
||||
|
||||
call App%raiApp%AdvanceModel(dt)
|
||||
|
||||
call raiCpl_PostAdvance(App)
|
||||
App%raiApp%State%t = App%raiApp%State%t + dt
|
||||
App%raiApp%State%ts = App%raiApp%State%ts + 1
|
||||
App%raiApp%State%mjd = T2MJD(dt,App%raiApp%State%mjd)
|
||||
|
||||
@@ -10,8 +10,8 @@ module raijuCplHelper
|
||||
use raijugrids
|
||||
use ioh5
|
||||
use files
|
||||
use arrayutil
|
||||
|
||||
use imagtubes
|
||||
use mixdefs
|
||||
use raijuColdStartHelper, only : initRaijuColdStarter
|
||||
|
||||
@@ -38,6 +38,8 @@ module raijuCplHelper
|
||||
|
||||
! Options
|
||||
call iXML%Set_Val(raiCpl%startup_blendTscl, "cpl/startupTscl", raiCpl%startup_blendTscl)
|
||||
call iXML%Set_Val(raiCpl%tsclSm_dL , "cpl/tsclSm_dL" , raiCpl%tsclSm_dL )
|
||||
call iXML%Set_Val(raiCpl%tsclSm_dMLT, "cpl/tsclSm_dMLT", raiCpl%tsclSm_dMLT)
|
||||
|
||||
! State sub-modules that need coupler settings
|
||||
call initRaijuColdStarter(raiCpl%raiApp%Model, iXML, raiCpl%raiApp%State%coldStarter,tEndO=raiCpl%startup_blendTscl)
|
||||
@@ -79,6 +81,9 @@ module raijuCplHelper
|
||||
enddo
|
||||
call initShellVar(raiCpl%shGr, SHGR_CC, raiCpl%tiote)
|
||||
call initShellVar(raiCpl%shGr, SHGR_CC, raiCpl%Tb)
|
||||
call initShellVar(raiCpl%shGr, SHGR_CC, raiCpl%avgBeta)
|
||||
|
||||
call initShellVar(raiCpl%shGr, SHGR_CC, raiCpl%tscl_mhdIngest)
|
||||
end associate
|
||||
|
||||
! Initial values
|
||||
@@ -146,8 +151,8 @@ module raijuCplHelper
|
||||
do i=1,NDIM
|
||||
call InterpShellVar_TSC_SG(raiCpl%shGr, raiCpl%xyzMin(i), raiCpl%shGr, raiCpl%xyzMincc(i), srcMaskO=topoSrcMask)
|
||||
enddo
|
||||
call InterpShellVar_TSC_SG(voltGrid, tubeShell%Tb, raiCpl%shGr, raiCpl%Tb, srcMaskO=topoSrcMask)
|
||||
|
||||
call InterpShellVar_TSC_SG(voltGrid, tubeShell%Tb , raiCpl%shGr, raiCpl%Tb , srcMaskO=topoSrcMask)
|
||||
call InterpShellVar_TSC_SG(voltGrid, tubeShell%avgBeta, raiCpl%shGr, raiCpl%avgBeta, srcMaskO=topoSrcMask)
|
||||
|
||||
end subroutine tubeShell2RaiCpl
|
||||
|
||||
@@ -219,50 +224,6 @@ module raijuCplHelper
|
||||
end associate
|
||||
end subroutine
|
||||
|
||||
!------
|
||||
! One-way driving from file helpers
|
||||
!------
|
||||
|
||||
subroutine genImagTubes(raiCpl, vApp)
|
||||
class(raijuCoupler_T), intent(inout) :: raiCpl
|
||||
type(voltApp_T), intent(in ) :: vApp
|
||||
|
||||
integer :: i,j
|
||||
real(rp) :: seedR, eqR
|
||||
type(magLine_T) :: magLine
|
||||
! Get field line info and potential from voltron
|
||||
! And put the data into RAIJU's fromV coupling object
|
||||
|
||||
associate(sh=>raiCpl%raiApp%Grid%shGrid , &
|
||||
planet=>raiCpl%raiApp%Model%planet, &
|
||||
ebApp =>vApp%ebTrcApp)
|
||||
|
||||
seedR = planet%ri_m/planet%rp_m
|
||||
! Do field line tracing, populate fromV%ijTubes
|
||||
!$OMP PARALLEL DO default(shared) &
|
||||
!$OMP schedule(dynamic) &
|
||||
!$OMP private(i,j,eqR)
|
||||
do i=sh%isg,sh%ieg+1
|
||||
do j=sh%jsg,sh%jeg+1
|
||||
call CleanLine(raiCpl%magLines(i,j))
|
||||
|
||||
eqR = DipColat2L(raiCpl%raiApp%Grid%thRp(i)) ! Function assumes colat coming from 1 Rp, make sure we use the right theta value
|
||||
if (eqR < raiCpl%opt%mhdRin) then
|
||||
call DipoleTube(vApp, sh%th(i), sh%ph(j), raiCpl%ijTubes(i,j))
|
||||
else
|
||||
call MHDTube(ebApp, planet, & !ebTrcApp, planet
|
||||
sh%th(i), sh%ph(j), seedR, & ! colat, lon, r
|
||||
raiCpl%ijTubes(i,j), raiCpl%magLines(i,j), & ! IMAGTube_T, magLine_T
|
||||
doShiftO=.true.)
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
end associate
|
||||
|
||||
end subroutine genImagTubes
|
||||
|
||||
|
||||
!------
|
||||
! Real-time coupling stuff
|
||||
!------
|
||||
@@ -283,6 +244,177 @@ module raijuCplHelper
|
||||
end subroutine
|
||||
|
||||
|
||||
!------
|
||||
! Post-advance calculations
|
||||
!------
|
||||
|
||||
subroutine raiCpl_PostAdvance(raiCpl)
|
||||
class(raijuCoupler_T), intent(inout) :: raiCpl
|
||||
|
||||
call raiCpl_calcTsclMHD(raiCpl)
|
||||
end subroutine raiCpl_PostAdvance
|
||||
|
||||
|
||||
subroutine raiCpl_calcTsclMHD(raiCpl)
|
||||
class(raijuCoupler_T), intent(inout) :: raiCpl
|
||||
|
||||
integer :: i,j
|
||||
real(rp) :: vaFrac_cc
|
||||
|
||||
associate(tscl=>raiCpl%tscl_mhdIngest, State=>raiCpl%raiApp%State, sh=>raiCpl%shGr)
|
||||
|
||||
! Defaults
|
||||
call fillArray(tscl%data, State%dt)
|
||||
where (State%active /= RAIJUINACTIVE)
|
||||
tscl%mask = .true.
|
||||
elsewhere
|
||||
tscl%mask = .false.
|
||||
end where
|
||||
|
||||
! First, calculate our tscl point-by-point
|
||||
!$OMP PARALLEL DO default(shared) &
|
||||
!$OMP private(i,j,vaFrac_cc)
|
||||
do j=sh%jsg,sh%jeg
|
||||
do i=sh%isg,sh%ieg
|
||||
if (tscl%mask(i,j)) then
|
||||
vaFrac_cc = 0.25*sum(State%vaFrac(i:i+1,j:j+1))
|
||||
tscl%data(i,j) = raiCpl%raiApp%Model%nBounce*State%dt/(vaFrac_cc)**2
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Now do smoothing
|
||||
call smoothRaijuVar_eq(State, sh, raiCpl%tsclSm_dMLT, raiCpl%tsclSm_dL, tscl)
|
||||
|
||||
! We had the mask include buffer cells for smoothing reasons, but now we want anyone interpolating to only use active cells
|
||||
where (State%active == RAIJUACTIVE)
|
||||
tscl%mask = .true.
|
||||
elsewhere
|
||||
tscl%mask = .false.
|
||||
end where
|
||||
|
||||
end associate
|
||||
|
||||
|
||||
contains
|
||||
|
||||
subroutine smoothRaijuVar_eq(State, sh, dMLT, dL, var)
|
||||
!! Smooth a raiju variable with stencil in equatorial projection
|
||||
type(raijuState_T), intent(in) :: State
|
||||
type(ShellGrid_T), intent(in) :: sh
|
||||
real(rp), intent(in) :: dMLT
|
||||
real(rp), intent(in) :: dL
|
||||
type(ShellGridVar_T), intent(inout) :: var
|
||||
!! Variable to smooth. Assumes var%mask can be used to include/exclude points
|
||||
|
||||
integer :: i, j, j_eval, ipnt, jpnt
|
||||
integer :: dj, dj_half, nGood
|
||||
logical :: doIScan
|
||||
real(rp) :: var_sum
|
||||
real(rp) :: L_center, L_pnt
|
||||
real(rp) :: dPhi_rad, dMLT_pnt, dL_pnt
|
||||
real(rp), dimension(:,:), allocatable :: tmp_var_sm
|
||||
|
||||
|
||||
if (.not. sh%isPhiUniform) then
|
||||
write(*,*) "ERROR: raiCpl_calcTsclMHD expects raiCpl ShellGrid to have periodic phi, but it does not"
|
||||
write(*,*) " Goodbye."
|
||||
stop
|
||||
endif
|
||||
|
||||
if (var%loc /= SHGR_CC ) then
|
||||
write(*,*) "ERROR: raiCpl_calcTsclMHD expects raiCpl ShellGridVar to be located at cell center"
|
||||
write(*,*) " Given var with location enum:",var%loc
|
||||
write(*,*) " Goodbye."
|
||||
stop
|
||||
endif
|
||||
|
||||
allocate(tmp_var_sm(var%isv:var%iev, var%jsv:var%jev))
|
||||
call fillArray(tmp_var_sm, 0.0_rp)
|
||||
|
||||
dphi_rad = dMLT * PI/12.0_rp ! Convert delta-MLT to a delta-phi in radians
|
||||
dj_half = 0
|
||||
do while(sh%ph(sh%js + dj_half + 1) - sh%phc(sh%js) < dphi_rad/2.0_rp) ! Increase dj_half until its dphi is half of target dphi_rad (you're welcome)
|
||||
dj_half = dj_half + 1
|
||||
enddo
|
||||
|
||||
!$OMP PARALLEL DO default(shared) &
|
||||
!$OMP schedule(dynamic) &
|
||||
!$OMP private(i, j, dj, doIScan) &
|
||||
!$OMP private(nGood, var_sum, ipnt, jpnt, L_center, L_pnt, dL_pnt)
|
||||
do j=sh%jsg,sh%jeg
|
||||
do i=sh%isg,sh%ieg
|
||||
if (.not. var%mask(i,j)) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
L_center = norm2(State%xyzMincc(i,j,XDIR:YDIR)) ! XY plane to L shell / radius from planet
|
||||
|
||||
nGood = 0
|
||||
var_sum = 0.0_rp
|
||||
! Loop over j stencil range
|
||||
do dj=-dj_half,dj_half
|
||||
jpnt = j + dj
|
||||
if (jpnt < sh%js) then
|
||||
jpnt = jpnt + sh%Np
|
||||
elseif (jpnt > sh%je) then
|
||||
jpnt = jpnt - sh%Np
|
||||
endif
|
||||
|
||||
! Sweep in i+ direction (earthward)
|
||||
doIScan = .true.
|
||||
ipnt = i
|
||||
do while (doIScan)
|
||||
L_pnt = norm2(State%xyzMincc(ipnt,jpnt,XDIR:YDIR))
|
||||
dL_pnt = abs(L_center - L_pnt)
|
||||
|
||||
! First decide if we are gonna keep going after this point
|
||||
if (ipnt >= sh%ieg .or. dL_pnt > dL/2.0_rp) then
|
||||
doIScan = .false.
|
||||
endif
|
||||
|
||||
! Decide if we should include this point
|
||||
if (var%mask(ipnt,jpnt) .and. (dL_pnt < dL/2.0_rp)) then
|
||||
nGood = nGood + 1
|
||||
var_sum = var_sum + var%data(ipnt,jpnt)
|
||||
endif
|
||||
ipnt = ipnt + 1
|
||||
enddo
|
||||
|
||||
! Same thing but in the -i direction
|
||||
doIScan = .true.
|
||||
ipnt = i-1
|
||||
do while(doIScan)
|
||||
L_pnt = norm2(State%xyzMincc(ipnt,jpnt,XDIR:YDIR))
|
||||
dL_pnt = abs(L_center - L_pnt)
|
||||
|
||||
if (ipnt <= sh%isg .or. dL_pnt > dL/2.0_rp) then
|
||||
doIScan = .false.
|
||||
endif
|
||||
|
||||
if (var%mask(ipnt,jpnt) .and. (dL_pnt < dL/2.0_rp)) then
|
||||
nGood = nGood + 1
|
||||
var_sum = var_sum + var%data(ipnt,jpnt)
|
||||
endif
|
||||
ipnt = ipnt - 1
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Now calculate and save our average
|
||||
if (nGood > 0) then
|
||||
tmp_var_sm(i,j) = var_sum/(1.0_rp*nGood)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Put it back into the original variable and we're done
|
||||
var%data(:,:) = tmp_var_sm(:,:)
|
||||
|
||||
end subroutine smoothRaijuVar_eq
|
||||
|
||||
end subroutine raiCpl_calcTsclMHD
|
||||
|
||||
|
||||
!------
|
||||
! Coupler I/O
|
||||
!------
|
||||
@@ -302,23 +434,24 @@ module raijuCplHelper
|
||||
call writeShellGrid(raiCpl%shGr, ResF,"/ShellGrid")
|
||||
|
||||
call ClearIO(IOVars)
|
||||
call AddOutVar(IOVars, "tLastUpdate", raiCpl%tLastUpdate, uStr="s")
|
||||
call AddOutSGV(IOVars, "Pavg" , raiCpl%Pavg , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Davg" , raiCpl%Davg , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Pstd" , raiCpl%Pstd , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Dstd" , raiCpl%Dstd , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Bmin" , raiCpl%Bmin , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "xyzMin" , raiCpl%xyzMin , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "xyzMincc" , raiCpl%xyzMincc , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "topo" , raiCpl%topo , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "thcon" , raiCpl%thcon , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "phcon" , raiCpl%phcon , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "bvol" , raiCpl%bvol , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "bvol_cc" , raiCpl%bvol_cc , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "vaFrac" , raiCpl%vaFrac , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Tb" , raiCpl%Tb , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "pot_total" , raiCpl%pot_total , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "pot_corot" , raiCpl%pot_corot , doWriteMaskO=.true.)
|
||||
call AddOutVar(IOVars, "tLastUpdate" , raiCpl%tLastUpdate , uStr="s")
|
||||
call AddOutSGV(IOVars, "Pavg" , raiCpl%Pavg , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Davg" , raiCpl%Davg , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Pstd" , raiCpl%Pstd , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Dstd" , raiCpl%Dstd , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Bmin" , raiCpl%Bmin , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "xyzMin" , raiCpl%xyzMin , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "xyzMincc" , raiCpl%xyzMincc , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "topo" , raiCpl%topo , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "thcon" , raiCpl%thcon , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "phcon" , raiCpl%phcon , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "bvol" , raiCpl%bvol , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "bvol_cc" , raiCpl%bvol_cc , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "vaFrac" , raiCpl%vaFrac , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "Tb" , raiCpl%Tb , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "pot_total" , raiCpl%pot_total , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "pot_corot" , raiCpl%pot_corot , doWriteMaskO=.true.)
|
||||
call AddOutSGV(IOVars, "tscl_mhdIngest", raiCpl%tscl_mhdIngest, doWriteMaskO=.true.)
|
||||
call WriteVars(IOVars, .false., ResF)
|
||||
call MapSymLink(ResF,lnResF)
|
||||
end subroutine
|
||||
@@ -347,23 +480,24 @@ module raijuCplHelper
|
||||
raiCpl%tLastUpdate = GetIOReal(IOVars, "tLastUpdate")
|
||||
|
||||
! ShellGridVars
|
||||
call ReadInSGV(raiCpl%Pavg ,ResF, "Pavg" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Davg ,ResF, "Davg" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Pstd ,ResF, "Pstd" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Dstd ,ResF, "Dstd" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Bmin ,ResF, "Bmin" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%xyzMin ,ResF, "xyzMin" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%xyzMincc ,ResF, "xyzMincc" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%topo ,ResF, "topo" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%thcon ,ResF, "thcon" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%phcon ,ResF, "phcon" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%bvol ,ResF, "bvol" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%bvol_cc ,ResF, "bvol_cc" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%vaFrac ,ResF, "vaFrac" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Tb ,ResF, "Tb" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%pot_total,ResF, "pot_total", doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%pot_corot,ResF, "pot_corot", doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Pavg ,ResF, "Pavg" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Davg ,ResF, "Davg" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Pstd ,ResF, "Pstd" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Dstd ,ResF, "Dstd" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Bmin ,ResF, "Bmin" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%xyzMin ,ResF, "xyzMin" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%xyzMincc ,ResF, "xyzMincc" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%topo ,ResF, "topo" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%thcon ,ResF, "thcon" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%phcon ,ResF, "phcon" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%bvol ,ResF, "bvol" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%bvol_cc ,ResF, "bvol_cc" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%vaFrac ,ResF, "vaFrac" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%Tb ,ResF, "Tb" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%pot_total ,ResF, "pot_total" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%pot_corot ,ResF, "pot_corot" , doIOpO=.false.)
|
||||
call ReadInSGV(raiCpl%tscl_mhdIngest,ResF, "tscl_mhdIngest", doIOpO=.false.)
|
||||
|
||||
end subroutine
|
||||
|
||||
end module raijuCplHelper
|
||||
end module raijuCplHelper
|
||||
|
||||
@@ -382,6 +382,8 @@ module voltapp_mpi
|
||||
call convertGameraToRemix(vApp%mhd2mix, vApp%gApp, vApp%remixApp)
|
||||
call Toc("G2R")
|
||||
|
||||
call MJDRecalc(vApp%MJD)
|
||||
|
||||
if (vApp%doGCM .and. vApp%time >=0 .and. vApp%gcmCplRank /= -1) then
|
||||
call Tic("GCM2MIX")
|
||||
call coupleGCM2MIX(vApp%gcm,vApp%remixApp%ion,vApp%MJD,vApp%time,vApp%mageCplComm,vApp%gcmCplRank)
|
||||
|
||||
@@ -12,6 +12,59 @@ module tubehelper
|
||||
|
||||
contains
|
||||
|
||||
|
||||
subroutine init_TubeShell(sh, tubeShell, TioTeO)
|
||||
type(ShellGrid_T), intent(in) :: sh
|
||||
!! Voltron shellgrid
|
||||
type(TubeShell_T), intent(inout) :: tubeShell
|
||||
!! IMAGTubeShell object we are initializing
|
||||
!logical, dimension(:,:), intent(in), optional :: maskO
|
||||
real(rp), intent(in), optional :: TioteO
|
||||
!! Default Ti/Te ratio
|
||||
|
||||
integer :: i
|
||||
real(rp) :: tiote
|
||||
|
||||
if (present(TioteO)) then
|
||||
tiote = TioteO
|
||||
else
|
||||
tiote = 4.0
|
||||
endif
|
||||
|
||||
|
||||
do i=1,NDIM
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%xyz0(i))
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%X_bmin(i))
|
||||
enddo
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%lat0 )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%lon0 )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%latc )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%lonc )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%invlat )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%topo )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%bmin )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%bVol )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%Lb )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%Tb )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%wMAG )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%rCurv )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%avgBeta)
|
||||
do i=0,MAXTUBEFLUIDS
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%avgP(i))
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%avgN(i))
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%stdP(i))
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%stdN(i))
|
||||
enddo
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%losscone )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%lossconec)
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%TioTe0 )
|
||||
call initShellVar(sh, SHGR_CORNER, tubeShell%nTrc )
|
||||
|
||||
tubeShell%TioTe0%data = tiote
|
||||
|
||||
|
||||
end subroutine init_TubeShell
|
||||
|
||||
subroutine FakeTube(P,xyz0,bTube)
|
||||
!! Given a seed xyz0 generate a fake tube that seems plausible
|
||||
type(planet_T), intent(in) :: P
|
||||
|
||||
@@ -2,10 +2,10 @@ module voltappHelper
|
||||
|
||||
use kdefs
|
||||
use voltTypes
|
||||
use imagtubes
|
||||
use kai2geo
|
||||
use mixinterfaceutils
|
||||
|
||||
use tubehelper
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
Reference in New Issue
Block a user