First cut at plasmasphere refilling

src/raiju/raijuEtautils.F90

@@ -427,17 +427,21 @@ module raijuetautils
         type(raijuState_T), intent(inout) :: State
 
         integer :: i,j,k0
-        real(rp) :: maxX
+        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
+        logical :: isGood
         
-        maxX = 2.0 !Max over-filling relative to target, i.e. don't go above maxX x den-target
         
         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,isGood,xeq,yeq,rad,dppT,cc2eta) &
-        !$OMP private()
-        
+        !$OMP private(dpsph,etaT,deta,dndt)
+
         do j=Grid%shGrid%jsg,Grid%shGrid%jeg
             do i=Grid%shGrid%isg,Grid%shGrid%ieg
                 isGood = (State%active(i,j) == RAIJUACTIVE)
@@ -448,7 +452,7 @@ module raijuetautils
                 rad = sqrt(xeq**2.0 + yeq**2.0)
 
                 !Closed field line, calculate Gallagher w/ current Kp to get target density
-                dppT = GallagherXY(xmin(i,j),ymin(i,j),NowKp)
+                dppT = GallagherXY(xeq,yeq,NowKp)
                 cc2eta = State%bvol_cc(i,j)*sclEta
                 eta2cc = 1.0/cc2eta !Convert eta to #/cc
 
@@ -458,11 +462,11 @@ module raijuetautils
                 if (dpsph >= maxX*dppT) cycle !Too much already there
                 etaT = dppT/eta2cc
 
-                if ((rad <= RInMHD) .and. (dppT > dpsph)) then
-                    !If this is inside MHD inner boundary, be at least at target value
-                    State%eta(i,j,k0) = etaT
-                    cycle
-                endif
+                ! if ((rad <= RInMHD) .and. (dppT > dpsph)) then
+                !     !If this is inside MHD inner boundary, be at least at target value
+                !     State%eta(i,j,k0) = etaT
+                !     cycle
+                ! endif
                 
                 !If still here then calculate refilling
                 dndt = 10.0**(3.48-0.331*rad) !cm^-3/day, Denton+ 2012 eqn 1

src/raiju/raijuStarter.F90

@@ -485,7 +485,7 @@ module raijustarter
             endif
 
             State%KpTS%wID = Model%tsF
-            call State%tsKp%initTS("Kp")
+            call State%KpTS%initTS("Kp")
         
         end associate