JHUAPL/kaiju
1
0
Fork 0
mirror of https://github.com/JHUAPL/kaiju.git synced 2026-01-09 18:57:53 -05:00
Code Issues Packages Projects Releases 1 Wiki Activity

Moving wave model code into 'new' raijuLoss format

Browse Source
This commit is contained in:
Anthony
2025-05-01 13:55:10 -06:00
parent c7239e0de4
commit 000c96caef
8 changed files with 531 additions and 390 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
src/base/defs/raijudefs.F90
View File

@@ -86,13 +86,6 @@ module raijudefs
!! Fraction that a lambda channel must contribute to total pressure or density in order to be worthy of being evolved
real(rp), parameter :: pressFracThreshDef = 0.01
!! If fraction of total pressure below lowest lambda bin is greater than this, we complain
! Wave model stuff
real(rp), parameter :: def_NpsphHigh = 100.0 ! [#/cc]
real(rp), parameter :: def_NpsphLow = 10.0 ! [#/cc]
real(rp), parameter :: def_ChorusLmax = 7.0 ! [Re]
real(rp), parameter :: def_PsheetLmin = 8.0 ! [Re]
real(rp), parameter :: def_ChorusEMin = 1.1 ! [keV]
! Coupling defaults
real(rp), parameter :: def_vaFracThresh = 0.10

106
src/base/types/raijuTypes.F90
View File

@@ -78,6 +78,8 @@ module raijutypes
!! Tell someone if this loss process applies to given species
procedure :: calcTau => baseLossCalcTau
!! Report instantaneous loss rate for a given energy and lat/lon
procedure :: doOutput => baseLossDoOutput
!! Output any relevant state info to file
end type baseRaijuLoss_T
type :: raijuLPHolder_T
@@ -88,44 +90,44 @@ module raijutypes
!------
! Precipitation models
!------
type eLossWM_T
!! Parameters used in electron wave model from Dedong Wang and Shanshan Bao
! -- Model -- !
real(rp) :: NpsphHigh = def_NpsphHigh
real(rp) :: NpsphLow = def_NpsphLow
real(rp) :: ChorusLMax = def_ChorusLmax
real(rp) :: PsheetLMin = def_PsheetLmin
real(rp) :: ChorusEMin = def_ChorusEMin
logical :: doOutput = .false.
!! Whether or not we will be asked to provide detailed info in the output file
type(TimeSeries_T) :: KpTS
!! Kp data from wind file
! -- Grid -- !
! Chorus info
integer :: Nkp, Nmlt, Nl, Ne
!! Number of bins for Kp, MLT, L shell, and Energy
real(rp), allocatable, dimension(:) :: Kp1D
!! 1D array of Kp dimension for Tau4D
real(rp), allocatable, dimension(:) :: MLT1D
!! 1D array of MLT dimension for Tau4D
real(rp), allocatable, dimension(:) :: L1D
!! 1D array of L shell dimension for Tau4D [Re]
real(rp), allocatable, dimension(:) :: Energy1D
!! 1D array of energy dimension for Tau4D [MeV]
real(rp), allocatable, dimension(:,:,:,:) :: Tau4D
!! Tau(Kp, MLT, L, E) table electron scattering table [seconds]
! -- State -- !
real(rp), allocatable, dimension(:,:) :: wPS
real(rp), allocatable, dimension(:,:) :: wHISS
real(rp), allocatable, dimension(:,:) :: wCHORUS
end type eLossWM_T
! type eLossWM_T
! !! Parameters used in electron wave model from Dedong Wang and Shanshan Bao
!
! ! -- Model -- !
! real(rp) :: NpsphHigh = def_NpsphHigh
! real(rp) :: NpsphLow = def_NpsphLow
! real(rp) :: ChorusLMax = def_ChorusLmax
! real(rp) :: PsheetLMin = def_PsheetLmin
! real(rp) :: ChorusEMin = def_ChorusEMin
!
! logical :: doOutput = .false.
! !! Whether or not we will be asked to provide detailed info in the output file
!
!
! type(TimeSeries_T) :: KpTS
! !! Kp data from wind file
!
! ! -- Grid -- !
! ! Chorus info
! integer :: Nkp, Nmlt, Nl, Ne
! !! Number of bins for Kp, MLT, L shell, and Energy
! real(rp), allocatable, dimension(:) :: Kp1D
! !! 1D array of Kp dimension for Tau4D
! real(rp), allocatable, dimension(:) :: MLT1D
! !! 1D array of MLT dimension for Tau4D
! real(rp), allocatable, dimension(:) :: L1D
! !! 1D array of L shell dimension for Tau4D [Re]
! real(rp), allocatable, dimension(:) :: Energy1D
! !! 1D array of energy dimension for Tau4D [MeV]
! real(rp), allocatable, dimension(:,:,:,:) :: Tau4D
! !! Tau(Kp, MLT, L, E) table electron scattering table [seconds]
!
! ! -- State -- !
! real(rp), allocatable, dimension(:,:) :: wPS
! real(rp), allocatable, dimension(:,:) :: wHISS
! real(rp), allocatable, dimension(:,:) :: wCHORUS
!
! end type eLossWM_T
!------
! General helpers
@@ -275,14 +277,16 @@ module raijutypes
!! Planet info like radius, mag. moment, etc.
! Losses
logical :: doSS, doCC, doCX, doFLC
!! (Ions) Do strong scattering / coulomb collisions / charge exchange / field-line curvature
integer :: eLossModel
!! Enumerator indicating active electon loss model
procedure(raijuELossRate_T), pointer, nopass :: eLossRateFn => NULL()
!! Pointer to electron loss function
type(eLossWM_T) :: eLossWM
!! Container for electron Wave Model data
logical :: doSS
!! Do strong scattering
logical :: doCC
!! Do coulomb collisions
logical :: doCX
!! Do charge exchange
logical :: doFLC
!! Do field-line curvature
logical :: doEWM
!! Do electron loss model
! Coupling info
integer :: nFluidIn = 0
@@ -656,6 +660,16 @@ module raijutypes
tau = HUGE
end function baseLossCalcTau
subroutine baseLossDoOutput(this, Model, Grid, State, gStr, doGhostsO)
!Import :: baseRaijuLoss_T, raijuModel_T, raijuGrid_T, raijuState_T, rp
class(baseRaijuLoss_T), intent(in) :: this
type(raijuModel_T), intent(in) :: Model
type(raijuGrid_T) , intent(in) :: Grid
type(raijuState_T), intent(in) :: State
character(len=strLen), intent(in) :: gStr
logical, intent(in), optional :: doGhostsO
end subroutine baseLossDoOutput

2
src/raiju/losses/raijuLoss_SS.F90
View File

@@ -27,7 +27,7 @@ module raijuLoss_SS
type(raijuGrid_T) , intent(in) :: Grid
type(XML_Input_T) , intent(in) :: xmlInp
! Idk what planet you are, electrons gonna scatter
! Idc what planet you are, electrons gonna scatter
this%reqsGood = .true.
! And they gonna precipitate
this%isPrecip = .true.

458
src/raiju/losses/raijuLoss_eWM_BW.F90 Normal file
View File

@@ -0,0 +1,458 @@
module raijuLoss_eWM_BW
! Bao-Wang electron Wave Model - based loss and precipitation model
use kdefs, only : mec2, vc_cgs
use raijudefs
use raijutypes
use raijuSpeciesHelper
use raijuLoss_SS
use arrayutil
implicit none
integer, parameter, private :: MAXIOVAR = 50
type, extends(baseRaijuLoss_T) :: raiLoss_eWM_BW_T
logical :: reqsGood = .false.
type(raiLoss_SS_T) :: loss_SS
!! Our personal instance of the strong scattering implementation to apply where needed
! -- Model -- !
real(rp) :: NpsphHigh = 100.0 ! [#/cc]
real(rp) :: NpsphLow = 10.0 ! [#/cc]
real(rp) :: ChorusLMax = 7.0 ! [Re]
real(rp) :: PsheetLMin = 8.0 ! [Re]
real(rp) :: ChorusEMin = 1.1 ! [keV]
type(TimeSeries_T) :: KpTS
!! Kp data from wind file
! -- Grid -- !
! Chorus info
integer :: Nkp, Nmlt, Nl, Ne
!! Number of bins for Kp, MLT, L shell, and Energy
real(rp), allocatable, dimension(:) :: Kp1D
!! 1D array of Kp dimension for Tau4D
real(rp), allocatable, dimension(:) :: MLT1D
!! 1D array of MLT dimension for Tau4D
real(rp), allocatable, dimension(:) :: L1D
!! 1D array of L shell dimension for Tau4D [Re]
real(rp), allocatable, dimension(:) :: Energy1D
!! 1D array of energy dimension for Tau4D [MeV]
real(rp), allocatable, dimension(:,:,:,:) :: Tau4D
!! Tau(Kp, MLT, L, E) table electron scattering table [seconds]
! -- State -- !
real(rp), allocatable, dimension(:,:) :: wPS
real(rp), allocatable, dimension(:,:) :: wHISS
real(rp), allocatable, dimension(:,:) :: wCHORUS
real(rp), allocatable, dimension(:,:,:) :: tauTotal
!! (i,j,k) Final taus considering all processes and weights
contains
procedure :: doInit => eWM_BW_LossInit
procedure :: doUpdate => eWM_BW_DoUpdate
procedure :: isValidSpc => eWM_BW_LossValidSpc
procedure :: calcTau => eWM_BW_CalcTau
procedure :: doOutput => eWM_BW_DoOutput
end type raiLoss_eWM_BW_T
contains
subroutine eWM_BW_LossInit(this, Model, Grid, xmlInp)
class(raiLoss_eWM_BW_T), intent(inout) :: this
type(raijuModel_T), intent(in) :: Model
type(raijuGrid_T) , intent(in) :: Grid
type(XML_Input_T) , intent(in) :: xmlInp
type(raijuSpecies_T) :: spc_ele
type(IOVAR_T), dimension(MAXIOVAR) :: IOVars
type(IOVar_T) :: tauVAR
! Only valid for Earth
if (trim(toUpper(Model%planet%name)) .ne. "EARTH") then
write(*,*)"WARNING in raijuLoss_eWM_BW: Not simulating Earth, idk what waves are like elswhere"
write(*,*)"No electron wave model happening"
return
else
this%reqsGood = .true.
endif
this%isPrecip = .true.
! If we are still here then we are gonna be used, set everything up
spc_ele = Grid%spc(spcIdx(Grid, F_HOTE))
! Init our instance of RAIJU strong scattering model
call this%loss_SS%doInit(Model, Grid, xmlInp)
! We need Kp from wind file, try to get that first
call xmlInp%Set_Val(this%KpTS%wID,"/Kaiju/Gamera/wind/tsfile","NONE")
call this%KpTS%initTS("Kp",doLoudO=.false.)
call ClearIO(IOVars)
if(.not. ioExist(Model%configFName, "Tau", "waveModel")) then
write(*,*) "This config file not structured for RAIJU, use genRAIJU.py. Good day."
stop
endif
!Chorus wave
call AddInVar(IOVars,"Kp")
call AddInVar(IOVars,"MLT")
call AddInVar(IOVars,"L")
call AddInVar(IOVars,"Ek")
call AddInVar(IOVars,"Tau")
call ReadVars(IOVars,.false.,Model%configFname, "waveModel")
this%Nkp = IOVars(FindIO(IOVars, "Kp" ))%N
this%Nmlt = IOVars(FindIO(IOVars, "MLT"))%N
this%Nl = IOVars(FindIO(IOVars, "L" ))%N
this%Ne = IOVars(FindIO(IOVars, "Ek" ))%N
tauVAR = IOVars(FindIO(IOVars, "Tau"))
! Do some dimension checks
if (tauVAR%Nr /= 4) then
write(*,*) "tauDim:",tauVAR%Nr
write(*,*) 'Currently only support tau model files in the form tau(Kp,MLT,L,Ek)'
stop
endif
if (this%Nkp /= tauVAR%dims(1) .or. this%Nmlt /= tauVAR%dims(2) .or. this%Nl /= tauVAR%dims(3) .or. this%Ne /= tauVAR%dims(4)) then
write(*,*) "tauDims:",tauVAR%dims,"Nk:",this%Nkp,"Nm:",this%Nmlt,"Nl:",this%Nl,"Ne:",this%Ne
write(*,*) 'Dimensions of tau arrays are not compatible'
stop
endif
! If still here we're good to allocate and store data
allocate(this%Kp1D (this%Nkp ))
allocate(this%MLT1D (this%Nmlt))
allocate(this%L1D (this%Nl ))
allocate(this%Energy1D(this%Ne ))
allocate(this%Tau4D(this%Nkp,this%Nmlt,this%Nl,this%Ne))
call IOArray1DFill(IOVars,"Kp" , this%Kp1D )
call IOArray1DFill(IOVars,"MLT", this%MLT1D )
call IOArray1DFill(IOVars,"L" , this%L1D )
call IOArray1DFill(IOVars,"Ek" , this%Energy1D)
call IOArray4DFill(IOVars,"Tau", this%Tau4D )
call ClearIO(IOVars)
!Array order check: array is in acsending order
!Chorus
if(this%Kp1D(1) > this%Kp1D(this%Nkp)) then
write(*,*) "Kp: ",this%Kp1D
write(*,*) "reorder wave model so Kp is in ascending order"
stop
end if
if(this%MLT1D(1) > this%MLT1D(this%Nmlt)) then
write(*,*) "MLT: ",this%MLT1D
write(*,*) "reorder wave model so MLT is in ascending order"
stop
end if
if(this%L1D(1) > this%L1D(this%Nl)) then
write(*,*) "L: ",this%L1D
write(*,*) "reorder wave model so L shell is in ascending order"
stop
end if
if(this%Energy1D(1) > this%Energy1D(this%Ne)) then
write(*,*) "Ek: ",this%Energy1D
write(*,*) "reorder wave model so Ek is in ascending order"
stop
end if
allocate(this%wPS (Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg))
allocate(this%wHISS (Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg))
allocate(this%wCHORUS(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg))
allocate(this%tauTotal(Grid%shGrid%isg:Grid%shGrid%ieg, Grid%shGrid%jsg:Grid%shGrid%jeg, spc_ele%kStart:spc_ele%kEnd))
end subroutine eWM_BW_LossInit
subroutine eWM_BW_DoUpdate(this, Model, Grid, State)
!! We are called at the start of every step, once coupling info, boundary conditions, etc. have beenr resolved
!! We calculate all of the weighting, loss taus, for all electrons, and store it here
!! When our calcTau gets called for an i,j,k, we'll just report the right value from our State
class(raiLoss_eWM_BW_T), intent(inout) :: this
type(raijuModel_T), intent(in) :: Model
type(raijuGrid_T) , intent(in) :: Grid
type(raijuState_T), intent(in) :: State
integer :: psphIdx, eleIdx
logical :: isGood
integer :: i,j,k
real(rp) :: NpsphPnt
!! Density [#/cc] of plasmasphere at point i,j
real(rp) :: L, MLT, E, Kp
!! L shell and MLT of given point, channel energy in MeV, current Kp
real(rp) :: wLBlend, wNBlend
!! L-weighting of blending between IMAG and PS. 0=PS
!! Density-weighting between Chorus and Hiss
real(rp) :: tauPS, tauHiss, tauCHORUS
! Zero everyone out in prep for new values
call fillArray(this%wPS , 0.0_rp)
call fillArray(this%wHISS , 0.0_rp)
call fillArray(this%wCHORUS , 0.0_rp)
call fillArray(this%tauTotal, HUGE )
call this%loss_SS%doUpdate(Model, Grid, State) ! Just in case its doUpdate actually does something in the future
! Things that aren't i,j,k-dependent
psphIdx = spcIdx(Grid, F_PSPH)
eleIdx = spcIdx(Grid, F_HOTE)
Kp = this%KpTS%evalAt(State%t)
associate(sh=>Grid%shGrid, spc=>Grid%spc(eleIdx))
!$OMP PARALLEL DO default(shared) &
!$OMP private(i,j,k,isGood,L,MLT,E,NpsphPnt,wNBlend,wLBlend,tauPS,tauHiss,tauCHORUS)
do j=sh%jsg,sh%jeg
do i=sh%isg,sh%ieg
isGood = State%active(i,j)
if (.not. isGood) then
cycle
endif
L = sqrt(State%xyzMincc(i,j,1)**2 + State%xyzMincc(i,j,2)**2) ! [Re]
MLT = atan2(State%xyzMincc(i,j,2),State%xyzMincc(i,j,1))/pi*12.D0+12.D0
E = abs(Grid%alamc(k) * State%bvol_cc(i,j)**(-2./3.)) * 1.0E-6 ! [MeV]
NpsphPnt = State%Den(psphIdx)%data(i,j)
! Calculate blending
wNBlend = dlog(NpsphPnt/this%NpsphLow) / dlog(this%NpsphHigh/this%NpsphLow)
call ClampValue(wNBlend, 0.0_rp, 1.0_rp)
!! 1 => Psphere Hiss, 0 => Other
wLBlend = RampDown(L, this%ChorusLMax, this%PsheetLMin - this%ChorusLMax)
!! 1 => Chorus, 0 => PS
! If psphere density is high enough, we always apply Hiss regardless of L
! This means Hiss model needs to appropriately handle any L value its given
! If psphere density is low, we choose between Chorus and plasma sheet strong-scattering based on L
! Also note: the weighting is not k-dependent, but if its not too costly to do for each k then ¯\_(ツ)_/¯
! Calculate weights
this%wHISS(i,j) = wNBlend
this%wCHORUS(i,j) = (1 - wNBlend)*wLBlend
this%wPS(i,j) = (1 - wNBlend)*(1-wLBlend)
! Now calculate loss taus and accumulate
! Energy-independent taus first
if (this%wHiss(i,j) > TINY) then
tauHiss = calcHissTau(MLT, L, E, Kp)
else
tauHiss = 0.0_rp
endif
if (this%wPS(i,j) > TINY) then
tauPS = this%loss_SS%calcTau(Model, Grid, State, i, j, k)
else
tauPS = 0.0_rp
endif
this%tauTotal(i,j,:) = this%wHISS(i,j)*tauHiss + this%wPS(i,j)*tauPS
if (this%wCHORUS(i,j) > TINY) then
do k=spc%kStart,spc%kEnd
!! Implement chorus tau calculation here
tauCHORUS = 0.0_rp
this%tauTotal(i,j,k) = this%tauTotal(i,j,k) + this%wCHORUS(i,j)*tauCHORUS
enddo
endif
enddo
enddo
end associate
end subroutine eWM_BW_DoUpdate
function eWM_BW_CalcTau(this, Model, Grid, State, i, j, k) result(tau)
!! We are not doing any heavy calculations since we did it in DoUpdate
!! Just give back thr right value after some checks
class(raiLoss_eWM_BW_T), intent(in) :: this
type(raijuModel_T ), intent(in) :: Model
type(raijuGrid_T ), intent(in) :: Grid
type(raijuState_T ), intent(in) :: State
integer, intent(in) :: i, j, k
real(rp) :: tau
associate(spc => Grid%spc(Grid%k2spc(k)))
tau = HUGE
if (.not. this%reqsGood) then
return
endif
if (k < spc%kStart .or. k > spc%kEnd) then
return
endif
! If still here then the requester is right to be calling us, give back our tau
tau = this%tauTotal(i,j,k)
end associate
end function eWM_BW_CalcTau
function eWM_BW_LossValidSpc(this, spc) result(isValid)
class(raiLoss_eWM_BW_T), intent(in) :: this
type(raijuSpecies_T), intent(in) :: spc
logical :: isValid
isValid = .false.
if (.not. this%reqsGood) return
if ( (spc%spcType .eq. RAIJUELE) ) then
isValid = .true.
endif
end function eWM_BW_LossValidSpc
subroutine eWM_BW_DoOutput(this, Model, Grid, State, gStr, doGhostsO)
class(raiLoss_eWM_BW_T), intent(in) :: this
type(raijuModel_T), intent(in) :: Model
type(raijuGrid_T) , intent(in) :: Grid
type(raijuState_T), intent(in) :: State
character(len=strLen), intent(in) :: gStr
!! Group ("Step#X") to append our data to in a sub-group ("eWM")
logical, intent(in), optional :: doGhostsO
character(len=strLen) :: gStrWM
integer :: is, ie, js, je, ks, ke
logical :: doGhosts
type(IOVAR_T), dimension(MAXIOVAR) :: IOVars
if (present(doGhostsO)) then
doGhosts = doGhostsO
else
doGhosts = .false.
endif
if (doGhosts) then
is = Grid%shGrid%isg
ie = Grid%shGrid%ieg
js = Grid%shGrid%jsg
je = Grid%shGrid%jeg
else
is = Grid%shGrid%is
ie = Grid%shGrid%ie
js = Grid%shGrid%js
je = Grid%shGrid%je
endif
write(gStrWM, "(A, A)") trim(gStr), "/eWM"
!Reset IO chain
call ClearIO(IOVars)
! Fat output for precip info
associate(spc=>Grid%spc(spcIdx(Grid, F_HOTE)))
call AddOutVar(IOVars, "wPS" , this%wPS (is:ie,js:je), dStr="0-1 weighting for plasma sheet strong scattering losses")
call AddOutVar(IOVars, "wHISS" , this%wHISS (is:ie,js:je), dStr="0-1 weighting for plasmasphere hiss losses")
call AddOutVar(IOVars, "wCHORUS", this%wCHORUS(is:ie,js:je), dStr="0-1 weighting for inner mag chorus wave losses")
call AddOutVar(IOVars, "tauTotal", this%tauTotal(is:ie,js:je,spc%kStart:spc%kEnd), dStr="[s] combined lifetimes of all contributing loss processes")
call WriteVars(IOVars,.true.,Model%raijuH5, gStrWM)
end associate
end subroutine eWM_BW_DoOutput
function calcHissTau(MLTin, Lin, Ein, Kpin) result(tau)
! Empirical lifetime against plasmaspheric hiss pitch angle diffusion, based on Orlova et al. 2015JA021878.
! Improvements relative to 2014GL060100: 1. Hiss wave intensity distribution model is based on new data
! (O14 was based on single-component E field in CRRES data. O16 used Spasojevic+2015 model based on EMFISIS B data on VAP);
! 2. Wave spectrum is assumed differently (O14 assume Gaussian spectrum based on CRRES data).
! Electron lifetime tau(L,E,MLT,Kp) = tau_av(L,E)/g(MLT)/h(Kp),
! where 1.5<L<5.5, E=log10(Ek[MeV]) for 1 KeV < Ek < 10 MeV.
! log10(tau_av(L,E)) = a1+a2*L+a3*E+...+a20*E^3, when E >= f(L).
! f(L) = 0.1328*L^2-2.1463*L+3.7857.
! g(MLT) = 10^g0(MLT)/G0
! h(Kp) = 10^h0(Kp)/H0
! G0 = int_0^24(10^g0(MLT))dMLT / 24 = 782.3.
! g0(MLT) = b2*MLT^2 + b1*MLT + b0
! H0 = 1315.
! h0(Kp) = c2*Kp^2 + c1*Kp + c0
real(rp), intent(in) :: MLTin, Lin
!! L in Re
real(rp), intent(in) :: Ein, Kpin
!! E in MeV
real(rp) :: tau_av, tau, rateK
!! tau in s, rate in 1/s
real(rp) :: MLT, L, E, Kp
!! Parames we clamp if needed
real(rp) :: L2, L3, L4, fL, E2, E3, E4, E5, LE
!! Polynomials
real(rp), dimension(20) :: le_pol
!! Container for polynomials
real(rp) :: b0, b1, b2, G0, g0_MLT, g_MLT, c0, c1, c2, H0, h0_Kp, h_Kp
!! Coefficients for g_MLT and h_Kp
real(rp), dimension(20), parameter :: a1_20 = [77.323, -92.641, -55.754, 44.497, 48.981, 8.9067, -10.704, &
-15.711, -3.3326, 1.5189, 1.294, 2.2546, 0.31889, -0.85916, &
-0.22182, 0.034318, 0.097248, -0.12192, -0.062765, 0.0063218]
!! Coefficients for polynomials
rateK = 0.0
MLT = MLTin
Kp = Kpin
L = Lin
call ClampValue(L,1.5_rp,5.5_rp)
L2 = L*L
fL = 0.1328*L2 - 2.1463*L + 3.7857
E = log10(Ein)
call ClampValue(E,max(-3.0_rp,fL),1.0_rp) ! 1 keV to 1 Mev
E2 = E*E
E3 = E2*E
E4 = E3*E
E5 = E4*E
L3 = L2*L
L4 = L3*L
LE = L*E
le_pol = (/1.D0,L,E,L2,LE,E2,L3,L2*E,L*E2,E3,L4,L3*E,L2*E2,L*E3,E4,L*E4,L2*E3,L4*E,L2*L3,E5/)
b0 = 2.080
b1 = 0.1773
b2 = -0.007338
G0 = 782.3
g0_MLT = b2*MLT*MLT + b1*MLT + b0
g_MLT = 10**g0_MLT/G0
c0 = 2.598
c1 = 0.2321
c2 = -0.01414
H0 = 1315.0
Kp = min(Kp,5.0_rp) ! 0<Kp<5, 1.2% data in Kp bin of 4.3-7.6
h0_Kp = c2*Kp*Kp + c1*Kp + c0
h_Kp = 10**h0_Kp/H0
tau_av = 10.0**(dot_product(a1_20,le_pol))*86400.D0 ! seconds
tau = tau_av/g_MLT/h_Kp
!rateK = 1.0/tau
!write(*,*) Ein, L, fL, E, tau, rateK, '--'
!stop
end function calcHissTau
end module raijuLoss_eWM_BW

309
src/raiju/raijuELossWM.F90
View File

@@ -1,309 +0,0 @@
module raijuELossWM
use ioh5
use files
use raijudefs
use raijutypes
use raijuspecieshelper, only : spcIdx
implicit none
integer, parameter, private :: MAXIOVAR = 50
contains
subroutine initEWM(eWM, configFname, xmlInp, shGrid)
type(eLossWM_T), intent(inout) :: eWM
character(len=strLen) :: configFname
type(XML_Input_T), intent(in) :: xmlInp
type(ShellGrid_T), intent(in) :: shGrid
type(IOVAR_T), dimension(MAXIOVAR) :: IOVars
type(IOVar_T) :: tauVAR
! Determine if we are expected to provide info to output file
call xmlInp%Set_Val(eWM%doOutput,"/Kaiju/RAIJU/losses/doOutput",eWM%doOutput)
! We need Kp from wind file, try to get that first
call xmlInp%Set_Val(eWM%KpTS%wID,"/Kaiju/Gamera/wind/tsfile","NONE")
call eWM%KpTS%initTS("Kp",doLoudO=.false.)
call ClearIO(IOVars)
if(.not. ioExist(configfname, "Tau", "waveModel")) then
write(*,*) "This config file not structured for RAIJU, use genRAIJU.py. Good day."
stop
endif
!Chorus wave
call AddInVar(IOVars,"Kp")
call AddInVar(IOVars,"MLT")
call AddInVar(IOVars,"L")
call AddInVar(IOVars,"Ek")
call AddInVar(IOVars,"Tau")
call ReadVars(IOVars,.false.,configFname, "waveModel")
eWM%Nkp = IOVars(FindIO(IOVars, "Kp"))%N
eWM%Nmlt = IOVars(FindIO(IOVars, "MLT"))%N
eWM%Nl = IOVars(FindIO(IOVars, "L"))%N
eWM%Ne = IOVars(FindIO(IOVars, "Ek"))%N
tauVAR = IOVars(FindIO(IOVars, "Tau"))
! Do some dimension checks
if (tauVAR%Nr /= 4) then
write(*,*) "tauDim:",tauVAR%Nr
write(*,*) 'Currently only support tau model files in the form tau(Kp,MLT,L,Ek)'
stop
endif
if (eWM%Nkp /= tauVAR%dims(1) .or. eWM%Nmlt /= tauVAR%dims(2) .or. eWM%Nl /= tauVAR%dims(3) .or. eWM%Ne /= tauVAR%dims(4)) then
write(*,*) "tauDims:",tauVAR%dims,"Nk:",eWM%Nkp,"Nm:",eWM%Nmlt,"Nl:",eWM%Nl,"Ne:",eWM%Ne
write(*,*) 'Dimensions of tau arrays are not compatible'
stop
endif
! If still here we're good to allocate and store data
allocate(eWM%Kp1D (eWM%Nkp ))
allocate(eWM%MLT1D (eWM%Nmlt))
allocate(eWM%L1D (eWM%Nl ))
allocate(eWM%Energy1D(eWM%Ne ))
allocate(eWM%Tau4D(eWM%Nkp,eWM%Nmlt,eWM%Nl,eWM%Ne))
call IOArray1DFill(IOVars,"Kp" , eWM%Kp1D )
call IOArray1DFill(IOVars,"MLT", eWM%MLT1D )
call IOArray1DFill(IOVars,"L" , eWM%L1D )
call IOArray1DFill(IOVars,"Ek" , eWM%Energy1D)
call IOArray4DFill(IOVars,"Tau", eWM%Tau4D )
call ClearIO(IOVars)
!Array order check: array is in acsending order
!Chorus
if(eWM%Kp1D(1) > eWM%Kp1D(eWM%Nkp)) then
write(*,*) "Kp: ",eWM%Kp1D
write(*,*) "reorder wave model so Kp is in ascending order"
stop
end if
if(eWM%MLT1D(1) > eWM%MLT1D(eWM%Nmlt)) then
write(*,*) "MLT: ",eWM%MLT1D
write(*,*) "reorder wave model so MLT is in ascending order"
stop
end if
if(eWM%L1D(1) > eWM%L1D(eWM%Nl)) then
write(*,*) "L: ",eWM%L1D
write(*,*) "reorder wave model so L shell is in ascending order"
stop
end if
if(eWM%Energy1D(1) > eWM%Energy1D(eWM%Ne)) then
write(*,*) "Ek: ",eWM%Energy1D
write(*,*) "reorder wave model so Ek is in ascending order"
stop
end if
if (eWM%doOutput) then
allocate(eWM%wPS (shGrid%isg:shGrid%ieg, shGrid%jsg:shGrid%jeg))
allocate(eWM%wHISS (shGrid%isg:shGrid%ieg, shGrid%jsg:shGrid%jeg))
allocate(eWM%wCHORUS(shGrid%isg:shGrid%ieg, shGrid%jsg:shGrid%jeg))
endif
end subroutine initEWM
function calcELossRate_WM(Model, Grid, State, i, j, k) result(lossRate2)
type(raijuModel_T) , intent(inout) :: Model
type(raijuGrid_T) , intent(in) :: Grid
type(raijuState_T) , intent(in) :: State
integer, intent(in) :: i, j, k
real(rp), dimension(2) :: lossRate2
integer :: psphIdx, eleIdx
real(rp) :: NpsphPnt
!! Density [#/cc] of plasmasphere at point i,j
real(rp) :: L, MLT, E, Kp
!! L shell and MLT of given point, channel energy in MeV, current Kp
real(rp) :: wLBlend, wNBlend
!! L-weighting of blending between IMAG and PS. 0=PS
!! Density-weighting between Chorus and Hiss
real(rp) :: wPS, wHISS, wCHORUS
!! Weights for plasma sheet, plasmasphere hiss, and chorus models
associate(eWM=>Model%eLossWM)
lossRate2 = 0.0
!! (1) = value, (2) = type
L = sqrt(State%xyzMin(i,j,1)**2 + State%xyzMin(i,j,2)**2) ! [Re]
MLT = atan2(State%xyzMin(i,j,2),State%xyzMin(i,j,1))/pi*12.D0+12.D0
E = abs(Grid%alamc(k) * State%bvol(i,j)**(-2./3.)) * 1.0E-6 ! [MeV]
Kp = eWM%KpTS%evalAt(State%t)
psphIdx = spcIdx(Grid, F_PSPH)
NpsphPnt = State%Den(psphIdx)%data(i,j)
! Calculate blending
wNBlend = dlog(NpsphPnt/eWM%NpsphLow) / dlog(eWM%NpsphHigh/eWM%NpsphLow)
call ClampValue(wNBlend, 0.0_rp, 1.0_rp)
!! 1 => Psphere Hiss, 0 => Other
wLBlend = RampDown(L, eWM%ChorusLMax, eWM%PsheetLMin - eWM%ChorusLMax)
!! 1 => Chorus, 0 => PS
! If psphere density is high enough, we always apply Hiss regardless of L
! This means Hiss model needs to appropriately handle any L value its given
! If psphere density is low, we choose between Chorus and plasma sheet strong-scattering based on L
! Also note: the weighting is not k-dependent, but if its not too costly to do for each k then ¯\_(ツ)_/¯
! Calculate weights
wHISS = wNBlend
wCHORUS = (1 - wNBlend)*wLBlend
wPS = (1 - wNBlend)*(1-wLBlend)
! Make weighting variable, idk if this will work well
lossRate2(2) = lossRate2(2) + int(wHISS*10)
lossRate2(2) = lossRate2(2) + int(wCHORUS*100)
lossRate2(2) = lossRate2(2) + int(wPS*1000)
! Calculate loss rates and accumulate
if (wHiss > TINY) then
lossRate2(1) = lossRate2(1) + wHISS*calcHissRate(MLT, L, E, Kp)
endif
! Save output info if expected
! All electron k's will use this function, only first k should write something
eleIdx = Grid%k2spc(k)
if (eWM%doOutput .and. k .eq. Grid%spc(eleIdx)%kStart) then
eWM%wPS(i,j) = wPS
eWM%wHISS(i,j) = wHISS
eWM%wCHORUS(i,j) = wCHORUS
endif
end associate
end function calcELossRate_WM
function calcHissRate(MLTin, Lin, Ein, Kpin) result(rateK)
! Empirical lifetime against plasmaspheric hiss pitch angle diffusion, based on Orlova et al. 2015JA021878.
! Improvements relative to 2014GL060100: 1. Hiss wave intensity distribution model is based on new data
! (O14 was based on single-component E field in CRRES data. O16 used Spasojevic+2015 model based on EMFISIS B data on VAP);
! 2. Wave spectrum is assumed differently (O14 assume Gaussian spectrum based on CRRES data).
! Electron lifetime tau(L,E,MLT,Kp) = tau_av(L,E)/g(MLT)/h(Kp),
! where 1.5<L<5.5, E=log10(Ek[MeV]) for 1 KeV < Ek < 10 MeV.
! log10(tau_av(L,E)) = a1+a2*L+a3*E+...+a20*E^3, when E >= f(L).
! f(L) = 0.1328*L^2-2.1463*L+3.7857.
! g(MLT) = 10^g0(MLT)/G0
! h(Kp) = 10^h0(Kp)/H0
! G0 = int_0^24(10^g0(MLT))dMLT / 24 = 782.3.
! g0(MLT) = b2*MLT^2 + b1*MLT + b0
! H0 = 1315.
! h0(Kp) = c2*Kp^2 + c1*Kp + c0
real(rp), intent(in) :: MLTin, Lin
!! L in Re
real(rp), intent(in) :: Ein, Kpin
!! E in MeV
real(rp) :: tau_av, tau, rateK
!! tau in s, rate in 1/s
real(rp) :: MLT, L, E, Kp
!! Parames we clamp if needed
real(rp) :: L2, L3, L4, fL, E2, E3, E4, E5, LE
!! Polynomials
real(rp), dimension(20) :: le_pol
!! Container for polynomials
real(rp) :: b0, b1, b2, G0, g0_MLT, g_MLT, c0, c1, c2, H0, h0_Kp, h_Kp
!! Coefficients for g_MLT and h_Kp
real(rp), dimension(20), parameter :: a1_20 = [77.323, -92.641, -55.754, 44.497, 48.981, 8.9067, -10.704, &
-15.711, -3.3326, 1.5189, 1.294, 2.2546, 0.31889, -0.85916, &
-0.22182, 0.034318, 0.097248, -0.12192, -0.062765, 0.0063218]
!! Coefficients for polynomials
rateK = 0.0
MLT = MLTin
Kp = Kpin
L = Lin
call ClampValue(L,1.5_rp,5.5_rp)
L2 = L*L
fL = 0.1328*L2 - 2.1463*L + 3.7857
E = log10(Ein)
call ClampValue(E,max(-3.0_rp,fL),1.0_rp) ! 1 keV to 1 Mev
E2 = E*E
E3 = E2*E
E4 = E3*E
E5 = E4*E
L3 = L2*L
L4 = L3*L
LE = L*E
le_pol = (/1.D0,L,E,L2,LE,E2,L3,L2*E,L*E2,E3,L4,L3*E,L2*E2,L*E3,E4,L*E4,L2*E3,L4*E,L2*L3,E5/)
b0 = 2.080
b1 = 0.1773
b2 = -0.007338
G0 = 782.3
g0_MLT = b2*MLT*MLT + b1*MLT + b0
g_MLT = 10**g0_MLT/G0
c0 = 2.598
c1 = 0.2321
c2 = -0.01414
H0 = 1315.0
Kp = min(Kp,5.0_rp) ! 0<Kp<5, 1.2% data in Kp bin of 4.3-7.6
h0_Kp = c2*Kp*Kp + c1*Kp + c0
h_Kp = 10**h0_Kp/H0
tau_av = 10.0**(dot_product(a1_20,le_pol))*86400.D0 ! seconds
tau = tau_av/g_MLT/h_Kp
rateK = 1.0/tau
!write(*,*) Ein, L, fL, E, tau, rateK, '--'
!stop
end function calcHissRate
subroutine eWMOutput(Model, Grid, State, gStr, doGhostsO)
type(raijuModel_T), intent(inout) :: Model
type(raijuGrid_T ), intent(in) :: Grid
type(raijuState_T), intent(in) :: State
character(len=strLen), intent(in) :: gStr
logical, optional, intent(in) :: doGhostsO
character(len=strLen) :: gStrWM
integer :: is, ie, js, je, ks, ke
logical :: doGhosts
type(IOVAR_T), dimension(MAXIOVAR) :: IOVars
if (present(doGhostsO)) then
doGhosts = doGhostsO
else
doGhosts = .false.
endif
if (doGhosts) then
is = Grid%shGrid%isg
ie = Grid%shGrid%ieg
js = Grid%shGrid%jsg
je = Grid%shGrid%jeg
else
is = Grid%shGrid%is
ie = Grid%shGrid%ie
js = Grid%shGrid%js
je = Grid%shGrid%je
endif
write(gStrWM, "(A, A)") trim(gStr), "/eWM"
!Reset IO chain
call ClearIO(IOVars)
! Fat output for precip info
call AddOutVar(IOVars, "wPS" , Model%eLossWM%wPS (is:ie,js:je), dStr="0-1 weighting for plasma sheet strong scattering losses")
call AddOutVar(IOVars, "wHISS" , Model%eLossWM%wHISS (is:ie,js:je), dStr="0-1 weighting for plasmasphere hiss losses")
call AddOutVar(IOVars, "wCHORUS", Model%eLossWM%wCHORUS(is:ie,js:je), dStr="0-1 weighting for inner mag chorus wave losses")
call WriteVars(IOVars,.true.,Model%raijuH5, gStrWM)
end subroutine eWMOutput
end module raijuELossWM

6
src/raiju/raijuIO.F90
View File

@@ -7,7 +7,6 @@ module raijuIO
use raijutypes
use raijuetautils
use raijuELossWM, only : eWMOutput
use shellGrid
use shellInterp
use shellGridIO
@@ -423,11 +422,6 @@ module raijuIO
call WriteVars(IOVars,.true.,Model%raijuH5, gStr)
! Let sub-models add stuff if they want
if (Model%doLosses .and. Model%eLossWM%doOutput) then
call eWMOutput(Model, Grid, State, gStr, doGhostsO)
endif
end subroutine WriteRaiju

2
src/raiju/raijuStarter.F90
View File

@@ -14,7 +14,6 @@ module raijustarter
use raijuRecon, only : maxOrderSupported
use raijuout
use raijuICHelpers
use raijuELossWM
use raijulosses, only : initRaiLosses
use raijuColdStartHelper, only : initRaijuColdStarter
@@ -235,6 +234,7 @@ module raijustarter
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",.false.)
!--- Output ---!
call iXML%Set_Val(Model%isLoud , "output/loudConsole",.false.)

31
src/raiju/raijulosses.F90
View File

@@ -10,6 +10,7 @@ module raijulosses
use raijuLoss_CX
use raijuLoss_CC
use raijuLoss_SS
use raijuLoss_eWM_BW
implicit none
@@ -26,7 +27,7 @@ module raijulosses
!! iterator
integer :: numLPs
!! number of loss proccces we're gonna have
integer :: initCX, initCC, initSS, initFLC
integer :: initCX, initCC, initSS, initFLC, initEWM
!! Flag for if we need to allocate and init this LP
State%lossRates = 0.0
@@ -40,26 +41,29 @@ module raijulosses
initCX = merge(1, 0, Model%doCX)
initCC = merge(1, 0, Model%doCC)
initSS = merge(1, 0, Model%doSS)
initFLC = merge(1, 0, Model%doFLC)
numLPs = initCX + initCC + initSS + initFLC
initEWM = merge(1, 0, Model%doEWM)
initSS = merge(1, 0, Model%doSS .and. .not. Model%doEWM) ! TODO: Kinda jank, fix later
numLPs = initCX + initCC + initEWM + initSS + initFLC
allocate(State%lps(numLPs))
do iLP=1,numLPs
if (allocated(State%lps(iLP)%p)) deallocate(State%lps(iLP)%p)
! Determine which loss process is next in line for initting
if (initCX==1) then
if (allocated(State%lps(iLP)%p)) deallocate(State%lps(iLP)%p)
allocate( raiLoss_CX_T :: State%lps(iLP)%p )
initCX = 0
elseif(initCC==1) then
if (allocated(State%lps(iLP)%p)) deallocate(State%lps(iLP)%p)
allocate( raiLoss_CC_T :: State%lps(iLP)%p )
initCC = 0
State%lp_cc_idx = iLP
elseif(initEWM==1) then
allocate( raiLoss_eWM_BW_T :: State%lps(iLP)%p )
initEWM = 0
elseif(initSS==1) then
if (allocated(State%lps(iLP)%p)) deallocate(State%lps(iLP)%p)
allocate( raiLoss_SS_T :: State%lps(iLP)%p )
initSS = 0
endif
@@ -178,19 +182,6 @@ module raijulosses
end subroutine calcChannelLossRates
!------
! High-level calc
!------
subroutine calcElectronLossRate(Model, Grid, State, k)
type(raijuModel_T), intent(inout) :: Model
type(raijuGrid_T), intent(in) :: Grid
type(raijuState_T), intent(inout) :: State
integer, intent(in) :: k
end subroutine calcElectronLossRate
!------
! Apply losses and calc useful info
!------