Updates to make energy by species work with SCE

There have been several updates since I started this branch, these
changes bring us to a working point with the current upstream changes

Author: Miles A. Engel

src/ModRamIO.f90

@@ -1125,7 +1125,7 @@ subroutine ram_hour_write(S) !Previously WRESULT in ram_all
               AVEFL(I,J,K)=AVEFL(I,J,K)+F(I,J,K,L)*WMU(L)
             ENDDO
             AVEFL(I,J,K)=AVEFL(I,J,K)/(MU(NPA)-MU(UPA(I)))
-            PRECFL=PRECFL+AVEFL(I,J,K)*PI*WE(K)
+            PRECFL=PRECFL+AVEFL(I,J,K)*PI*WE(S,K)
           ENDDO
           WRITE(UNITTMP_,70) LZ(I),PHI(J),PRECFL
         END DO
@@ -1148,7 +1148,7 @@ subroutine ram_hour_write(S) !Previously WRESULT in ram_all
         if (outsideMGNP(I,J) == 1) F(I,J,:,:) = 1e-31
           DO 822 K=2,NE
 !            WRITE(20,31) T/3600.,LZ(I),KP,MLT(J),EKEV(K),F(I,J,K,2)
-           WRITE(30,31) T/3600.,LZ(I),KP,MLT(J),EKEV(K),F(I,J,K,27)
+           WRITE(30,31) T/3600.,LZ(I),KP,MLT(J),EKEV(S,K),F(I,J,K,27)
 !	    WRITE(40,31) T/3600.,LZ(I),KP,MLT(J),EKEV(K),F(I,J,K,UPA(I))
 822         CONTINUE
       ENDDO

src/ModRamRun.f90

@@ -588,7 +588,7 @@ SUBROUTINE ANISCH(S)
                    taudaa_he = ATAW_emic_he(I,J,K,L)/MUBOUN/BOUNHS(I,J,L)/(1.-MUBOUN*MUBOUN)
 
                    IF(DoSaveLwgr .and. (MOD(INT(T),3600).EQ.0.)) THEN
-                   write(24,556)ekev(k), PABn(L), taudaa_h, taudaa_he, &
+                   write(24,556)ekev(S,k), PABn(L), taudaa_h, taudaa_he, &
                         1./taudaa_h, 1./taudaa_he
                    ENDIF
 

src/ModRamSce.f90

@@ -268,7 +268,7 @@ subroutine calculate_precip_flux_jr(IS, nTheta, nPhi, rIono, energy_flux, ave_e,
           do i=1, npsi
              do j=2, nzeta
                 do k=1, nE
-                   write(UnitTmp_,'(f8.4,1x,f8.4,1x,f8.4,1x,E14.6)')rGrid(i,j), aGrid(i,j), EKEV(k), ave_fluxEQ(i,j,k)                   
+                   write(UnitTmp_,'(f8.4,1x,f8.4,1x,f8.4,1x,E14.6)')rGrid(i,j), aGrid(i,j), EKEV(iS,k), ave_fluxEQ(i,j,k) 
                 end do
              end do
           end do

src/ModSceIono.f90

@@ -9,7 +9,7 @@ MODULE ModSceIono
 
   contains
 !==================================================================================================
-  subroutine ionosphere_fluxes(iModel)
+  subroutine ionosphere_fluxes(iS, iModel)
 
     !\
     ! Combine the diffusive and discrete precipitating fluxes
@@ -26,7 +26,7 @@ subroutine ionosphere_fluxes(iModel)
 
     implicit none
 
-    integer, intent(in) :: iModel
+    integer, intent(in) :: iS, iModel
     real(DP), allocatable :: discrete_ef(:,:), discrete_ae(:,:), &
                              diffuse_ef(:,:), diffuse_ae(:,:)
     !---------------------------------------------------------------------------
@@ -76,7 +76,7 @@ subroutine ionosphere_fluxes(iModel)
   end subroutine ionosphere_fluxes
 
 !==================================================================================================
-  subroutine ionosphere_conductance(Sigma0, SigmaH, SigmaP, SigmaThTh, SigmaThPs, &
+  subroutine ionosphere_conductance(iS, Sigma0, SigmaH, SigmaP, SigmaThTh, SigmaThPs, &
                                     SigmaPsPs, dSigmaThTh_dTheta, dSigmaThPs_dTheta, &
                                     dSigmaPsPs_dTheta, dSigmaThTh_dPsi, dSigmaThPs_dPsi, &
                                     dSigmaPsPs_dPsi, Eflux, Ave_E, Eflux_diff, Theta, Psi, nTheta, &
@@ -111,7 +111,7 @@ subroutine ionosphere_conductance(Sigma0, SigmaH, SigmaP, SigmaThTh, SigmaThPs,
 
     implicit none
 
-    integer, intent(in)     :: iModel
+    integer, intent(in)     :: iS, iModel
     integer, intent(in)     :: nTheta, nPsi
     real(DP), intent(inout) :: Sigma0(:,:), SigmaH(:,:), SigmaP(:,:), SigmaThTh(:,:), &
                                SigmaThPs(:,:), SigmaPsPs(:,:), dSigmaThTh_dTheta(:,:), &
@@ -279,7 +279,7 @@ subroutine ionosphere_conductance(Sigma0, SigmaH, SigmaP, SigmaThTh, SigmaThPs,
              if (EFlux(i,j)*1000. > 0.001 .and. Ave_E(i,j)/2.*1000. >1)then
                 call glow_aurora_conductance(idate, real(ut,DP), glat(i,j), glong(i,j), &
                      SigmaP_Glow(i,j), SigmaH_Glow(i,j), &
-                     EFlux(i,j)*1000., Ave_E(i,j)/2., EFlux_Diff(i,j,:), EkeV(:), nE, &
+                     EFlux(i,j)*1000., Ave_E(i,j)/2., EFlux_Diff(i,j,:), EkeV(iS,:), nE, &
                      real(ap,DP), real(f107r,DP), real(f107p,DP),&
                      real(f107a,DP), DoUseFullSpec, &
                      zz(i,j,1:nzGlow), ionrate(i,j,1:nzGlow), eDen(i,j,1:nzGlow), &
@@ -448,7 +448,7 @@ subroutine ionosphere_conductance(Sigma0, SigmaH, SigmaP, SigmaThTh, SigmaThPs,
   end subroutine ionosphere_conductance
 
 !==================================================================================================
-  subroutine ionosphere_solver(Jr, SigmaThTh, SigmaThPs, SigmaPsPs, &
+  subroutine ionosphere_solver(iS, Jr, SigmaThTh, SigmaThPs, SigmaPsPs, &
                                dSigmaThTh_dTheta, dSigmaThPs_dTheta, &
                                dSigmaPsPs_dTheta, dSigmaThTh_dPsi, &
                                dSigmaThPs_dPsi, dSigmaPsPs_dPsi, &
@@ -516,6 +516,7 @@ subroutine ionosphere_solver(Jr, SigmaThTh, SigmaThPs, SigmaPsPs, &
 
     implicit none
 
+    integer, intent(in) :: iS
     real(DP), intent(inout) :: Jr(:,:), SigmaThTh(:,:), SigmaThPs(:,:), SigmaPsPs(:,:), &
                                dSigmaThTh_dTheta(:,:), dSigmaThPs_dTheta(:,:), &
                                dSigmaPsPs_dTheta(:,:), dSigmaThTh_dPsi(:,:), &

src/ModSceRun.f90

@@ -82,8 +82,8 @@ subroutine sce_solve
     iono_north_im_eflux_diff(:,Iono_nPsi,:) = iono_north_im_eflux_diff(:,1,:)
 
     IONO_NORTH_JR  = iono_north_im_jr
-    call ionosphere_fluxes(Conductance_Model)
-    call ionosphere_conductance(IONO_NORTH_Sigma0, IONO_NORTH_SigmaH, IONO_NORTH_SigmaP, &
+    call ionosphere_fluxes(S, Conductance_Model)
+    call ionosphere_conductance(S, IONO_NORTH_Sigma0, IONO_NORTH_SigmaH, IONO_NORTH_SigmaP, &
                                 IONO_NORTH_SigmaThTh, IONO_NORTH_SigmaThPs, IONO_NORTH_SigmaPsPs, &
                                 IONO_NORTH_dSigmaThTh_dTheta, IONO_NORTH_dSigmaThPs_dTheta, &
                                 IONO_NORTH_dSigmaPsPs_dTheta, IONO_NORTH_dSigmaThTh_dPsi, &
@@ -92,7 +92,7 @@ subroutine sce_solve
                                 IONO_NORTH_Theta, IONO_NORTH_Psi, IONO_nTheta, IONO_nPsi, &
                                 dTheta_North, dPsi_North, IONO_NORTH_X, IONO_NORTH_Y, IONO_NORTH_Z,&
                                 Conductance_Model)
-    call ionosphere_solver(IONO_NORTH_JR, IONO_NORTH_SigmaThTh, IONO_NORTH_SigmaThPs, &
+    call ionosphere_solver(S, IONO_NORTH_JR, IONO_NORTH_SigmaThTh, IONO_NORTH_SigmaThPs, &
                            IONO_NORTH_SigmaPsPs, IONO_NORTH_dSigmaThTh_dTheta, &
                            IONO_NORTH_dSigmaThPs_dTheta, IONO_NORTH_dSigmaPsPs_dTheta, &
                            IONO_NORTH_dSigmaThTh_dPsi, IONO_NORTH_dSigmaThPs_dPsi, &