WIP: Tracking CC heat flux

src/base/types/raijuTypes.F90

@@ -403,8 +403,10 @@ module raijutypes
 
         ! Loss-related things
         type(raijuLPHolder_T), dimension(:), allocatable :: lps
-        !class(baseRaijuLoss_T), dimension(:), allocatable :: lps
-            !! Collection of loss processes
+            !! Array of loss processes
+        integer :: lp_cc_idx = -1
+            !! If CC active, this is the index of CC object within lps array
+        
         ! (Ngi, Ngj, Nk) Varibles coming from RAIJU
         real(rp), dimension(:,:,:), allocatable :: lossRates
             !! (Ngi, Ngj, Nk) [1/s] Loss rates for each grid and lambda point. Generally stays the same over coupling time so we store them all here
@@ -417,7 +419,9 @@ module raijutypes
         real(rp), dimension(:,:,:), allocatable :: precipEFlux
             !! (Ngi, Ngj, Nk) [erg/cm^2/s] Precipitation energy fluxes
         real(rp), dimension(:,:,:), allocatable :: dEta_dt
-            !! (Ngi, Ngj, Ngk) [eta units/s] Average 
+            !! (Ngi, Ngj, Nk) [eta units/s] Average 
+        real(rp), dimension(:,:,:), allocatable :: CCHeatFlux
+            !! (Ngi, Ngj, Nk) [erg/cm^2/s] Heat flux from RC ions to plasmasphere electrons due to Coulumb collisions
         
         ! (Ngi, Ngj, Nspc+1) (First Nspc index is bulk) Moments
         ! Last dimension will be D/P of different populations (not necessarily same as species)

src/raiju/raijuAdvancer.F90

@@ -150,6 +150,7 @@ module raijuAdvancer
                 State%precipNFlux(:,:,k) = State%precipNFlux(:,:,k)/State%dt
                 State%precipEFlux(:,:,k) = State%precipEFlux(:,:,k)/State%dt
                 State%dEta_dt(:,:,k) = State%dEta_dt(:,:,k)/State%dt
+                State%CCHeatFlux(:,:,k) = State%CCHeatFlux(:,:,k)/State%dt
             endif
 
             State%nStepk(k) = State%nStepk(k) + n

src/raiju/raijuIO.F90

@@ -171,7 +171,7 @@ module raijuIO
         logical , dimension(:,:,:), allocatable :: tmpGood3D
         !real(rp), dimension(:,:), allocatable :: outActiveShell, outEnt
         !real(rp), dimension(:,:,:), allocatable :: outDen, outIntensity
-        real(rp), dimension(:,:,:), allocatable :: outPrecipN, outPrecipE, outPrecipAvgE
+        real(rp), dimension(:,:,:), allocatable :: outPrecipN, outPrecipE, outPrecipAvgE, outCCHeatFlux
         !real(rp), dimension(:,:,:), allocatable :: outPEff  ! effective Potential
 
         if (present(doGhostsO)) then
@@ -328,11 +328,14 @@ module raijuIO
             allocate(outPrecipN   (is:ie,js:je,Grid%nSpc))
             allocate(outPrecipE   (is:ie,js:je,Grid%nSpc))
             allocate(outPrecipAvgE(is:ie,js:je,Grid%nSpc))
+            allocate(outCCHeatFlux(is:ie,js:je,Grid%nSpc))
+            
             do s=1,Grid%nSpc
                 ks = Grid%spc(s)%kStart
                 ke = Grid%spc(s)%kEnd
-                outPrecipN(:,:,s) = sum(State%precipNFlux(is:ie,js:je,kS:kE), dim=3)
-                outPrecipE(:,:,s) = sum(State%precipEFlux(is:ie,js:je,kS:kE), dim=3)
+                outPrecipN(:,:,s)    = sum(State%precipNFlux(is:ie,js:je,kS:kE), dim=3)
+                outPrecipE(:,:,s)    = sum(State%precipEFlux(is:ie,js:je,kS:kE), dim=3)
+                outCCHeatFlux(:,:,s) = sum(State%CCHeatFlux (is:ie,js:je,kS:kE), dim=3)
 
                 where (outPrecipN(:,:,s) > TINY)
                     outPrecipAvgE(:,:,s) = outPrecipE(:,:,s)/outPrecipN(:,:,s) * erg2kev  ! Avg E [keV]
@@ -343,9 +346,11 @@ module raijuIO
             call AddOutVar(IOVars, "precipNFlux", outPrecipN   , uStr="#/cm^2/s")
             call AddOutVar(IOVars, "precipEFlux", outPrecipE   , uStr="erg/cm^2/s")
             call AddOutVar(IOVars, "precipAvgE" , outPrecipAvgE, uStr="keV")
+            call AddOutVar(IOVars, "CCHeatFlux" , outCCHeatFlux, uStr="erg/cm^2/s")
             deallocate(outPrecipN)
             deallocate(outPrecipE)
             deallocate(outPrecipAvgE)
+            deallocate(outCCHeatFlux)
         endif
 
         

src/raiju/raijuPreAdvancer.F90

@@ -43,6 +43,7 @@ module raijuPreAdvancer
         State%precipType_ele = 0.0
         State%precipNFlux = 0.0
         State%precipEFlux = 0.0
+        State%CCHeatFlux = 0.0
 
         ! Moments to etas, initial active shell calculation
         call Tic("BCs")

src/raiju/raijuStarter.F90

@@ -402,6 +402,7 @@ module raijustarter
             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) )
+            allocate( State%CCHeatFlux     (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%Nk) )
             ! Coupling output data
             allocate( State%Den  (sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%nSpc+1) )
             allocate( State%Press(sh%isg:sh%ieg, sh%jsg:sh%jeg, Grid%nSpc+1) )

src/raiju/raijulosses.F90

@@ -32,6 +32,7 @@ module raijulosses
         State%lossRates = 0.0
         State%lossRatesPrecip = 0.0
         State%dEta_dt = 0.0
+        State%CCHeatFlux = 0.0
 
         if (.not. Model%doLosses) then
             return
@@ -56,6 +57,7 @@ module raijulosses
                 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(initSS==1) then
                 if (allocated(State%lps(iLP)%p)) deallocate(State%lps(iLP)%p)
                 allocate( raiLoss_SS_T :: State%lps(iLP)%p )
@@ -184,7 +186,8 @@ module raijulosses
             !! Time delta [s]
 
         integer :: i,j
-        real(rp) :: deleta, eta0, pNFlux
+        real(rp) :: deleta, eta0, pNFlux, tau
+
         
         ! ! !$OMP PARALLEL DO default(shared) collapse(1) &
         ! ! !$OMP schedule(dynamic) &
@@ -205,9 +208,22 @@ module raijulosses
                 pNFlux = deleta2NFlux(deleta, Model%planet%rp_m, Grid%Brcc(i,j), dt)
                 State%precipNFlux(i,j,k) = State%precipNFlux(i,j,k) + pNFlux
                 State%precipEFlux(i,j,k) = State%precipEFlux(i,j,k) + nFlux2EFlux(pNFlux, Grid%alamc(k), State%bVol_cc(i,j))
+
+                ! Do special stuff for Coulomb collision effects
+                if (Model%doCC) then
+                    ! We can estimate heat transfer to plasmasphere electrons by energy lost from RC species to CC
+                    ! So we can follow same prodecure as above, by using just CC tau and dividing later by coupling dt to get average heat flux
+                    ! Treating this separately from precipication since its not actually precipitating ions
+                    tau = max(TINY, State%lps(State%lp_cc_idx)%p%calcTau(Model, Grid, State, i,j,k))
+                    deleta = eta0*(1.0 - exp(-dt/tau))
+                    pNFlux = deleta2NFlux(deleta, Model%planet%rp_m, Grid%Brcc(i,j), dt)
+                    State%CCHeatFlux(i,j,k) = State%CCHeatFlux(i,j,k) + nFlux2EFlux(pNFlux, Grid%alamc(k), State%bvol_cc(i,j))
+                endif
             enddo
         enddo
 
+        
+
     end subroutine applyLosses