diff --git a/Makefile b/Makefile index 4676738..6714465 100755 --- a/Makefile +++ b/Makefile @@ -236,14 +236,10 @@ test1_run: cd ${TESTDIR1}; ${MPIRUN} ./ram_scb.exe | tee runlog test1_check: - ${SCRIPTDIR}/DiffNum.pl -b -a=1e-9 \ - ${TESTDIR1}/output_ram/log_d20130317_t000000.log \ - ${IMDIR}/output/test1/log.ref \ - > test1.diff ${SCRIPTDIR}/DiffNum.pl -b -a=1e-9 \ ${TESTDIR1}/output_ram/pressure_d20130317_t001500.dat \ ${IMDIR}/output/test1/pressure.ref \ - >> test1.diff + > test1.diff ncdump -v "Flux_H","B_xyz" \ ${TESTDIR1}/output_ram/sat1_d20130317_t000000.nc \ | sed -e '1,/data:/d' > \ @@ -450,7 +446,38 @@ testEMIC_check: >> testEMIC.diff @echo "Test Successful!" +testSpecies: + @echo "starting..." > testSpecies.diff + @echo "testSpecies_compile..." >> testSpecies.diff + make testSpecies_compile + @echo "testSpecies_rundir..." >> testSpecies.diff + make testSpecies_rundir PARAMFILE=PARAM.in.testSpecies + @echo "testSpecies_run..." >> testSpecies.diff + make testSpecies_run MPIRUN= + @echo "testSpecies_check..." >> testSpecies.diff + make testSpecies_check + +testSpecies_compile: + make +testSpecies_rundir: + rm -rf ${TESTDIRC} + make rundir RUNDIR=${TESTDIRC} STANDALONE="YES" + cp Param/${PARAMFILE} ${TESTDIRC}/PARAM.in + cp input/sat*.dat ${TESTDIRC}/ + +testSpecies_run: + cd ${TESTDIRC}; ${MPIRUN} ./ram_scb.exe | tee runlog; + +testSpecies_check: + ${SCRIPTDIR}/DiffNum.pl -b -a=1e-9 \ + ${TESTDIRC}/output_ram/log_d20130317_t000000.log \ + ${IMDIR}/output/testSpecies/log.ref \ + > testSpecies.diff + ${SCRIPTDIR}/DiffNum.pl -b -a=1e-9 \ + ${TESTDIRC}/output_ram/pressure_d20130317_t001500.dat \ + ${IMDIR}/output/testSpecies/pressure.ref \ + >> testSpecies.diff #TEST SCE---------------------------------- testSCE: @@ -485,4 +512,5 @@ testSCE_check: ${TESTDIRC}/output_ram/pressure_d20130317_t001000.dat \ ${IMDIR}/output/testSCE/pressure.ref \ >> testSCE.diff + @echo "Test Successful!" diff --git a/input/initialization.nc b/input/initialization.nc index 4fbc762..ab3ddf8 100644 Binary files a/input/initialization.nc and b/input/initialization.nc differ diff --git a/src/ModRamBoundary.f90 b/src/ModRamBoundary.f90 index c2dd7f1..6cc871a 100644 --- a/src/ModRamBoundary.f90 +++ b/src/ModRamBoundary.f90 @@ -116,8 +116,8 @@ subroutine get_geomlt_flux(NameParticleIn, fluxOut_II) ! outside of the files energy range. rE = 0 lE = 0 - if (eGrid_SI(iSpec,1).gt.EkeV(1)) lE = 1 - if (eGrid_SI(iSpec,NEL_).lt.EkeV(nE)) rE = 1 + if (eGrid_SI(iSpec,1).gt.EkeV(iSpec, 1)) lE = 1 + if (eGrid_SI(iSpec,NEL_).lt.EkeV(iSpec, nE)) rE = 1 pE = rE + lE allocate(flux_II(0:NTL,NEL_+pE), logFlux_II(nT,NEL_+pE), logELan(NEL_+pE), logERam(nE)) flux_II = 0.0; logFlux_II = 0.0; logELan = 0.0; logERam = 0.0 @@ -201,7 +201,7 @@ subroutine get_geomlt_flux(NameParticleIn, fluxOut_II) logELan(1+lE:NEL_+le) = log10(eGrid_SI(iSpec,1:NEL_)) if (lE.eq.1) logELan(1) = log10(0.1000) if (rE.eq.1) logELan(NEL_+pE) = log10(1000.00) - logERam = log10(Ekev) + logERam = log10(Ekev(iSpec,:)) ! Interpolate/Extrapolate in energy space; return to normal units. ! Interpolation in Log space isn't really needed with the GSL interpolation diff --git a/src/ModRamCoul.f90 b/src/ModRamCoul.f90 index e13ebdc..55ac9ee 100644 --- a/src/ModRamCoul.f90 +++ b/src/ModRamCoul.f90 @@ -87,8 +87,8 @@ SUBROUTINE COULPARA(S) !...Collisions with plasmaspheric ions COULI(S,K,1)=-CCI*VBND(S,K)*CCO*GRBND(S,K)**2 !...Energy deposition rate [eV/cm2/s] in 1 hemisphere, check: - CEDR(S,K,1)=EDRCO*EKEV(K)*EDRE*(GREL(S,K)+1)/GREL(S,K)**2/V(S,K)**2 - CIDR(S,K,1)=EDRCO*EKEV(K)*EDRI*(GREL(S,K)+1)/GREL(S,K)**2/V(S,K)**2 + CEDR(S,K,1)=EDRCO*EKEV(S,K)*EDRE*(GREL(S,K)+1)/GREL(S,K)**2/V(S,K)**2 + CIDR(S,K,1)=EDRCO*EKEV(S,K)*EDRI*(GREL(S,K)+1)/GREL(S,K)**2/V(S,K)**2 !...Coulomb scattering CCDE=CCD*CDE*GREL(S,K)/(GREL(S,K)**2-1)**(1.5) CCDI=CCD*CDI*GREL(S,K)/(GREL(S,K)**2-1)**(1.5) @@ -152,7 +152,7 @@ SUBROUTINE COULEN(S) T=TimeRamElapsed - EZERO=EKEV(1)-WE(1) + EZERO=EKEV(S,1)-WE(S,1) GRZERO(S)=1.+EZERO*1000.*Q/RMAS(S)/CS/CS F(NE+1)=0. F(NE+2)=0. @@ -201,14 +201,14 @@ SUBROUTINE COULEN(S) IF (R.LE.0) FBND(K)=FUP IF (R.GT.0) THEN LIMITER=MAX(MIN(BetaLim*R,1.),MIN(R,BetaLim)) - CORR=-0.5*(CccolE(K)/DE(K)-ISIGN)*X + CORR=-0.5*(CccolE(K)/DE(S,K)-ISIGN)*X FBND(K)=FUP+LIMITER*CORR END IF END IF END DO DO K=2,NE - F2(S,I,J,K,L)=F2(S,I,J,K,L)-CccolE(K)/WE(K)*FBND(K)+CccolE(K-1)/WE(K)*FBND(K-1) + F2(S,I,J,K,L)=F2(S,I,J,K,L)-CccolE(K)/WE(S,K)*FBND(K)+CccolE(K-1)/WE(S,K)*FBND(K-1) if (f2(s,i,j,k,l) < 0._dp) then ! write(*,*) 'in COULEN f2<0 ', S,i,j,k,l, f2(S,i,j,k,l) f2(S,i,j,k,l)=1E-15 diff --git a/src/ModRamCouple.f90 b/src/ModRamCouple.f90 index 3b68777..56687b6 100755 --- a/src/ModRamCouple.f90 +++ b/src/ModRamCouple.f90 @@ -484,8 +484,6 @@ subroutine generate_flux ! Convert from moments to fluxes: do iE=1, nE - eCent = 1000.0 * Ekev(iE) ! Center of energy bin in eV. - factor=4.0E6*Ekev(iE) ! Unit conversion factor * eCent in KeV. do iT=1, nT select case(TypeMhd) case('single') @@ -493,6 +491,8 @@ subroutine generate_flux ! Assuming a maxwellian distribution about MHD temperature, divide ! densities into energy bins and calculate fluxes. do iS = 1, nS + eCent = 1000.0 * Ekev(iS,iE) ! Center of energy bin in eV. + factor=4.0E6*Ekev(iS,iE) ! Unit conversion factor * eCent in KeV. FluxBats_IIS(iE,iT,iS) = 1./sqrt(species(iS)%s_mass) & *exp(-1.0*eCent/MhdDensPres_VII(2,iT,iS)) & *factor*MhdDensPres_VII(1,iT,iS) & diff --git a/src/ModRamDrift.f90 b/src/ModRamDrift.f90 index e84f228..cb37570 100644 --- a/src/ModRamDrift.f90 +++ b/src/ModRamDrift.f90 @@ -68,7 +68,7 @@ SUBROUTINE DRIFTPARA(S) P1(I)=DTs/DPHI/2/MDR/RLZ(I) END DO DO K=1,NE - P2(I,K)=DTs*EKEV(K)*1000*(GREL(S,K)+1)/GREL(S,K)/RLZ(I)**2/DPHI/QS + P2(I,K)=DTs*EKEV(S,K)*1000*(GREL(S,K)+1)/GREL(S,K)/RLZ(I)**2/DPHI/QS END DO END DO @@ -80,7 +80,7 @@ SUBROUTINE DRIFTPARA(S) END DO MUDOT(I,NPA)=0. DO K=1,NE - EDOT(I,K)=EBND(K)*DTs/RLZ(I)*(GRBND(S,K)+1)/GRBND(S,K)/2. + EDOT(I,K)=EBND(S,K)*DTs/RLZ(I)*(GRBND(S,K)+1)/GRBND(S,K)/2. END DO END DO @@ -128,7 +128,7 @@ SUBROUTINE DRIFTR(S) ! Gradient Curvature Radial Drift DO K=1,NE - P4=DTs*EKEV(K)*1000.0*(GREL(S,K)+1)/GREL(S,K)/DPHI/MDR/QS + P4=DTs*EKEV(S,K)*1000.0*(GREL(S,K)+1)/GREL(S,K)/DPHI/MDR/QS DO L=1,NPA DO J=1,NT F(1:NR) = F2(S,:,J,K,L) @@ -307,7 +307,7 @@ SUBROUTINE DRIFTE(S) DtDriftE(S)=10000.0 OME=7.3E-5 - EZERO=EKEV(1)-WE(1) + EZERO=EKEV(S,1)-WE(S,1) GRZERO(S)=1.+EZERO*1000.*Q/RMAS(S)/CS/CS F(NE+1)=0. F(NE+2)=0. @@ -334,7 +334,7 @@ SUBROUTINE DRIFTE(S) F(1) = F(2)*GREL(S,1)/GREL(S,2)*SQRT((GREL(S,2)**2-1)/(GREL(S,1)**2-1)) F(0) = F(1)*GRZERO(S)/GREL(S,1)*SQRT((GREL(S,1)**2-1)/(GRZERO(S)**2-1)) DO K=1,NE - EDT1 = EBND(K)*1e3*(GRBND(S,K)+1)/2/GRBND(S,K)/FNHS(I,J,L)/RLZ(I)/BNES(I,J)/QS + EDT1 = EBND(S,K)*1e3*(GRBND(S,K)+1)/2/GRBND(S,K)/FNHS(I,J,L)/RLZ(I)/BNES(I,J)/QS DRDT = DRD1+EDT1*DRD2*RLZ(I) DPDT = DPD1-EDT1*DPD2 dBdt1 = dBdt(I,J)*(1.-FNIS(I,J,L)/2./FNHS(I,J,L))*RLZ(I)/BNES(I,J) @@ -342,7 +342,7 @@ SUBROUTINE DRIFTE(S) CDriftE(I,J,K,L) = EDOT(I,K)*((GPR1+GPR2+GPR3)*DRDT+(GPP1+GPP2)*DPDT+dBdt1+dIdt1) if (outsideMGNP(i,j) == 0) then ctemp = max(abs(CDriftE(I,J,K,L)),1E-10) - DtDriftE(S) = min(DtDriftE(S), FracCFL*DTs*DE(K)/ctemp) + DtDriftE(S) = min(DtDriftE(S), FracCFL*DTs*DE(S,K)/ctemp) endif sgn = 1 if (CDriftE(I,J,K,L).lt.0) sgn = -1 @@ -355,13 +355,13 @@ SUBROUTINE DRIFTE(S) IF (R.LE.0) FBND(K)=FUP IF (R.GT.0) THEN LIMITER=MAX(MIN(BetaLim*R,1.),MIN(R,BetaLim)) - CORR=-0.5*(CDriftE(I,J,K,L)/DE(K)-sgn)*X + CORR=-0.5*(CDriftE(I,J,K,L)/DE(S,K)-sgn)*X FBND(K)=FUP+LIMITER*CORR END IF END IF END DO DO K=2,NE - F2(S,I,J,K,L)=F2(S,I,J,K,L)-CDriftE(I,J,K,L)/WE(K)*FBND(K)+CDriftE(I,J,K-1,L)/WE(K)*FBND(K-1) + F2(S,I,J,K,L)=F2(S,I,J,K,L)-CDriftE(I,J,K,L)/WE(S,K)*FBND(K)+CDriftE(I,J,K-1,L)/WE(S,K)*FBND(K-1) if (f2(s,i,j,k,l).lt.0) then ! write(*,*) 'in DRIFTE f2<0 ', S,i,j,k,l, f2(S,i,j,k,l) f2(S,i,j,k,l)=1E-15 @@ -421,7 +421,7 @@ SUBROUTINE DRIFTMU(S) GMR3 = (BOUNIS(I+1,J,L)-BOUNIS(I-1,J,L))/2/MDR/BOUNIS(I,J,L) GMP1 = (BNES(I,J1)-BNES(I,J0))/4/DPHI/BNES(I,J) GMP2 = (BOUNIS(I,J1,L)-BOUNIS(I,J0,L))/2/DPHI/BOUNIS(I,J,L) - EDT = EKEV(K)*1e3*(GREL(S,K)+1)/2/GREL(S,K)/BOUNHS(I,J,L)/RLZ(I)/BNES(I,J)/QS + EDT = EKEV(S,K)*1e3*(GREL(S,K)+1)/2/GREL(S,K)/BOUNHS(I,J,L)/RLZ(I)/BNES(I,J)/QS DRM2 = (BOUNIS(I,J1,L)-BOUNIS(I,J0,L))/2/DPHI+(BOUNIS(I,J,L)-2*BOUNHS(I,J,L)) & *(BNES(I,J1)-BNES(I,J0))/4/BNES(I,J)/DPHI DPM2 = BOUNIS(I,J,L)+(BOUNIS(I+1,J,L)-BOUNIS(I-1,J,L))*RLZ(I)/2/MDR & diff --git a/src/ModRamFunctions.f90 b/src/ModRamFunctions.f90 index c2976d1..91ebbe5 100755 --- a/src/ModRamFunctions.f90 +++ b/src/ModRamFunctions.f90 @@ -59,7 +59,7 @@ subroutine get_ramdst(dstOut) if(l.ge.uPa(i))cycle do j=1, nT-1 if (outsideMGNP(i,j) == 0) then - sumEnergy=sumEnergy+f2(s,i,j,k,l)*wE(k)*wMu(L)*eKeV(k) + sumEnergy=sumEnergy+f2(s,i,j,k,l)*wE(S,k)*wMu(L)*eKeV(s,k) endif end do end do; end do; end do; end do diff --git a/src/ModRamIO.f90 b/src/ModRamIO.f90 index 5475689..8e71687 100644 --- a/src/ModRamIO.f90 +++ b/src/ModRamIO.f90 @@ -561,14 +561,11 @@ subroutine read_initial integer :: i, j, k, l, iS, GSLerr, iRDim, iTDim, iEDim, iPaDim, & iR, iT, iE, iPa, S - integer :: iFluxEVar, iFluxHVar, iFluxHeVar, iFluxOVar, & - iFileID, iStatus, iPParTVar, iPPerTVar - real(DP), allocatable :: iF2(:,:,:,:,:), iFNHS(:,:,:), iFNIS(:,:,:), iBOUNHS(:,:,:), & - iBOUNIS(:,:,:), iEIR(:,:), iEIP(:,:), iBNES(:,:), & - iHDNS(:,:,:), idBdt(:,:), idIdt(:,:,:), idIbndt(:,:,:), & - iLz(:), iMLT(:), iEkeV(:), iPaVar(:), radGrid(:,:), & - angleGrid(:,:), iPParT(:,:,:), iPPerT(:,:,:), & - tPParT(:,:,:), tPPerT(:,:,:) + integer :: iFluxVar, iFluxEVar, iFluxHVar, iFluxHeVar, iFluxOVar, & + iFileID, iStatus, iEkeVVar + real(DP), allocatable :: iF2(:,:,:,:,:), iLz(:), iMLT(:), iEkeV(:,:), iPaVar(:), & + radGrid(:,:), angleGrid(:,:), energyGrid(:,:), pitchGrid(:,:), & + jF2(:,:,:,:) real(DP) :: DL1, DPHI character(len=100) :: NameFile @@ -597,74 +594,71 @@ subroutine read_initial iStatus = nf90_inquire_dimension(iFileID, iPaDim, len = iPa) ! ALLOCATE ARRAYS FOR DATA - ALLOCATE(iF2(nS,iR,iT,iE,iPa), iLz(iR+1), iMLT(iT), iEkeV(iE), iPaVar(iPa)) - ALLOCATE(iFNHS(iR+1,iT,iPa), iBOUNHS(iR+1,iT,iPa), idBdt(iR+1,iT), iHDNS(iR+1,iT,iPa), & - iFNIS(iR+1,iT,iPa), iBOUNIS(iR+1,iT,iPa), idIdt(iR+1,iT,iPa), iBNES(iR+1,iT), & - idIbndt(iR+1,iT,iPa), iEIR(iR+1,iT), iEIP(iR+1,iT), tPParT(4,iR,iT), & - tPPerT(4,iR,iT),iPParT(nS,iR,iT),iPPerT(nS,iR,iT)) + ALLOCATE(iF2(nS,iR,iT,iE,iPa), iLz(iR+1), iMLT(iT), iEkeV(nS, iE), iPaVar(iPa)) + ! ALLOCATE ARRAYS FOR INTERPOLATION + ALLOCATE(radGrid(nR,nT),angleGrid(nR,nT),energyGrid(nE,nPa),pitchGrid(nE,nPa),jF2(iR,iT,nE,nPa)) + !ALLOCATE(jF2(iR,iT,nE,nPa)) ! GET VARIABLE IDS - !! FLUXES - iStatus = nf90_inq_varid(iFileID, 'FluxE', iFluxEVar) - iStatus = nf90_inq_varid(iFileID, 'FluxH', iFluxHVar) - iStatus = nf90_inq_varid(iFileID, 'FluxHe', iFluxHeVar) - iStatus = nf90_inq_varid(iFileID, 'FluxO', iFluxOVar) - - !! PRESSURES - iStatus = nf90_inq_varid(iFileID, 'PParT', iPParTVar) - iStatus = nf90_inq_varid(iFileID, 'PPerT', iPPerTVar) - - ! READ DATA - !! FLUXES - do i = 1, nS - select case(species(i)%s_name) - case('Electron') - iStatus = nf90_get_var(iFileID, iFluxEVar, iF2(i,:,:,:,:)) - case('Hydrogen') - iStatus = nf90_get_var(iFileID, iFluxHVar, iF2(i,:,:,:,:)) - case('HeliumP1') - iStatus = nf90_get_var(iFileID, iFluxHeVar, iF2(i,:,:,:,:)) - case('OxygenP1') - iStatus = nf90_get_var(iFileID, iFluxOVar, iF2(i,:,:,:,:)) - iF2(i,:,:,:,:) = (1. - OfracN)*iF2(i,:,:,:,:) - case('Nitrogen') - ! If we want to initialize some nitrogen, we assume that a percentage - ! of the oxygen in the initialization file is actually nitrogen - iStatus = nf90_get_var(iFileID, iFluxOVar, iF2(i,:,:,:,:)) - iF2(i,:,:,:,:) = OfracN*iF2(i,:,:,:,:) - case default - iF2(i,:,:,:,:) = 0._dp - end select - enddo - - !! PRESSURES - iStatus = nf90_get_var(iFileID, iPParTVar, tPParT(:,:,:)) - iStatus = nf90_get_var(iFileID, iPPerTVar, tPPerT(:,:,:)) - do i = 1, nS - select case(species(i)%s_name) - case('Electron') - iPPerT(i,:,:) = tPPerT(1,:,:) - iPParT(i,:,:) = tPParT(1,:,:) - case('Hydrogen') - iPPerT(i,:,:) = tPPerT(2,:,:) - iPParT(i,:,:) = tPParT(2,:,:) - case('HeliumP1') - iPPerT(i,:,:) = tPPerT(3,:,:) - iPParT(i,:,:) = tPParT(3,:,:) - case('OxygenP1') - iPPerT(i,:,:) = tPPerT(4,:,:) - iPParT(i,:,:) = tPParT(4,:,:) - case default - iPPerT(i,:,:) = 0._dp - iPParT(i,:,:) = 0._dp - end select - enddo - - + iStatus = nf90_inq_varid(iFileID, 'FluxE', iFluxVar) + if (iStatus == nf90_noerr) then + iStatus = nf90_inq_varid(iFileID, 'FluxE', iFluxEVar) + iStatus = nf90_inq_varid(iFileID, 'FluxH', iFluxHVar) + iStatus = nf90_inq_varid(iFileID, 'FluxHe', iFluxHeVar) + iStatus = nf90_inq_varid(iFileID, 'FluxO', iFluxOVar) + iStatus = nf90_inq_varid(iFileID, 'EnergyGrid', iEkeVVar) + do i = 1, nS + iStatus = nf90_get_var(iFileID, iEkeVVar, iEkeV(i, :)) + select case(species(i)%s_name) + case('Electron') + iStatus = nf90_get_var(iFileID, iFluxEVar, iF2(i,:,:,:,:)) + case('Hydrogen') + iStatus = nf90_get_var(iFileID, iFluxHVar, iF2(i,:,:,:,:)) + case('HeliumP1') + iStatus = nf90_get_var(iFileID, iFluxHeVar, iF2(i,:,:,:,:)) + case('OxygenP1') + iStatus = nf90_get_var(iFileID, iFluxOVar, iF2(i,:,:,:,:)) + iF2(i,:,:,:,:) = (1. - OfracN)*iF2(i,:,:,:,:) + case('Nitrogen') + ! If we want to initialize some nitrogen, we assume that a percentage + ! of the oxygen in the initialization file is actually nitrogen + iStatus = nf90_get_var(iFileID, iFluxOVar, iF2(i,:,:,:,:)) + iF2(i,:,:,:,:) = OfracN*iF2(i,:,:,:,:) + case default + iF2(i,:,:,:,:) = 0._dp + end select + enddo + else + do iS = 1, nS + iStatus = nf90_inq_varid(iFileID, 'EnergyGrid_'//species(iS)%s_name, iEkeVVar) + if (iStatus /= nf90_noerr) then + iEkeV(iS, :) = 0.0 + else + iStatus = nf90_get_var(iFileID, iEkeVVar, iEkeV(iS, :)) + endif + + iStatus = nf90_inq_varid(iFileID, 'Flux_'//species(iS)%s_name, iFluxVar) + if (iStatus /= nf90_noerr) then + iF2(iS,:,:,:,:) = 0.0 + else + iStatus = nf90_get_var(iFileID, iFluxVar, iF2(iS,:,:,:,:)) + endif + + select case(species(iS)%s_name) + case('OxygenP1') + iF2(iS,:,:,:,:) = (1. - OfracN)*iF2(iS,:,:,:,:) + case('Nitrogen') + iF2(iS,:,:,:,:) = OfracN*iF2(iS,:,:,:,:) + end select + enddo + endif ! CLOSE INITIALIZATION FILE iStatus = nf90_close(iFileID) call ncdf_check(iStatus, NameSub) + ! In the future we may want to be able to interpolate onto a different RadiusMax/RadiusMin + ! but for now we are assuming the same sized spatial grid, just potentially different number + ! of grid points. DL1 = (RadiusMax - RadiusMin)/(iR - 1) DO I=1,iR+1 iLz(I)=2.+(I-2)*DL1 @@ -675,50 +669,41 @@ subroutine read_initial iMLT(J)=(J-1)*DPHI END DO - iEkeV = EkeV - iPaVar = Pa + ! This is a hack for now, it will cause errors if the dimension of pitch angle bins isn't the same + iPaVar = cos(Pa*PI_d/180.) - ! Now we need to check the dimensions of the initialization file and - ! interpolate if they are different - if ((nR.eq.iR).and.(nT.eq.iT)) then - F2 = iF2 - PParT = iPParT - PPerT = iPPerT - else - ! Interpolate spatially for each energy and pitch angle (if required) - ALLOCATE(radGrid(nR+1,nT),angleGrid(nR+1,nT)) - DO i=1,NR+1 + do iS=1,nS + DO i=1,NR radGrid(i,:) = Lz(i) ENDDO DO j=1,NT angleGrid(:,j) = MLT(j)*2*PI_d/24 ENDDO - do iS=1,nS - CALL GSL_Interpolation_2D(iLz(1:iR), iMLT, iPParT(iS,:,:), radGrid(1:nR,:), & - angleGrid(1:nR,:), PParT(iS,:,:), GSLerr) - CALL GSL_Interpolation_2D(iLz(1:iR), iMLT, iPPerT(iS,:,:), radGrid(1:nR,:), & - angleGrid(1:nR,:), PPerT(iS,:,:), GSLerr) + DO k=1,nE + energyGrid(k,:) = EkeV(iS,k) + ENDDO + DO l=1,nPa + pitchGrid(:,l) = cos(Pa(l)*PI_d/180.0) + ENDDO + jF2(:,:,:,:) = 0.0 + do i=1,iR + do j=1,iT + CALL GSL_Interpolation_2D(iEkeV(iS,:), iPaVar, iF2(iS,i,j,:,:), energyGrid, & + pitchGrid, jF2(i,j,:,:), GSLerr) + enddo enddo - do l=1,nPa - do k=1,nE - do iS=1,nS - CALL GSL_Interpolation_2D(iLz(1:iR), iMLT, iF2(iS,:,:,k,l), radGrid(1:nR,:), & - angleGrid(1:nR,:), F2(iS,1:nR,:,k,l), GSLerr) - enddo + do k=1,nE + do l=1,nPa + CALL GSL_Interpolation_2D(iLz(1:iR), iMLT, jF2(:,:,k,l), radGrid(1:nR,:), angleGrid(1:nR,:), & + F2(iS,1:nR,:,k,l), GSLerr) enddo enddo - DEALLOCATE(radGrid,angleGrid) - - ! Interpolate across pitch angle - if ((nPa.ne.iPa).or.(nE.ne.iE)) then - call CON_stop('Changing pitch angle and energy resolution not currently supported') - endif - endif + enddo F2(:,:,:,1,:) = F2(:,:,:,2,:) DEALLOCATE(iF2, iEkeV, iMLT, iLz, iPaVar) - DEALLOCATE(iFNHS, iBOUNHS, iFNIS, iBOUNIS, iBNES, iHDNS, iEIR, iEIP, & - idBdt, idIbndt, idIdt, iPParT, iPPerT) + DEALLOCATE(radGrid, angleGrid, energyGrid, pitchGrid, jF2) + !DEALLOCATE(jF2) end subroutine read_initial @@ -759,10 +744,6 @@ subroutine write2DFlux iStatus = nf90_def_dim(iFileID, 'nE', nE, nEDim) iStatus = nf90_def_dim(iFileID, 'nPa', nPa, nPaDim) - iStatus = nf90_def_var(iFileID, 'EnergyGrid', nf90_double, & - (/nEDim/), iGridVar) - iStatus = nf90_put_var(iFileID, iGridVar, EKEV(:)) - iStatus = nf90_def_var(iFileID, 'PitchAngleGrid', nf90_double, & (/nPaDim/), iGridVar) iStatus = nf90_put_var(iFileID, iGridVar, PA(:)) @@ -778,6 +759,10 @@ subroutine write2DFlux !! FLUXES allocate(F(nR,nT,nE,nPa)) do S = 1, nS + iStatus = nf90_def_var(iFileID, 'EnergyGrid'//species(S)%s_name, nf90_double, & + (/nEDim/), iGridVar) + iStatus = nf90_put_var(iFileID, iGridVar, EKEV(S,:)) + iStatus = nf90_def_var(iFileID, 'Flux'//species(S)%s_name, nf90_double, & (/nRDim,nTDim,nEDim,nPaDim/), iFluxVar) iStatus = nf90_def_var_deflate(iFileID, iFluxVar, 0, yDeflate, iDeflate) @@ -848,13 +833,13 @@ subroutine writeLosses DO L=1,NPA DO J=1,NT-1 IF (L.LT.UPA(I)) THEN - WEIGHT=F2(S,I,J,K,L)*WE(K)*WMU(L)/FFACTOR(S,I,K,L)/FNHS(I,J,L) + WEIGHT=F2(S,I,J,K,L)*WE(S,K)*WMU(L)/FFACTOR(S,I,K,L)/FNHS(I,J,L) IF (MLT(J).LE.6.OR.MLT(J).GE.18.) THEN XNN(S,I)=XNN(S,I)+WEIGHT - ENERN(S,I)=ENERN(S,I)+EKEV(K)*WEIGHT + ENERN(S,I)=ENERN(S,I)+EKEV(S,K)*WEIGHT ELSE XND(S,I)=XND(S,I)+WEIGHT - ENERD(S,I)=ENERD(S,I)+EKEV(K)*WEIGHT + ENERD(S,I)=ENERD(S,I)+EKEV(S,K)*WEIGHT ENDIF ENDIF END DO @@ -1066,13 +1051,13 @@ subroutine ram_hour_write(S) !Previously WRESULT in ram_all if (f2(S,i,j,k,l).lt.1E-5) f2(S,i,j,k,l)=1E-5 if (f(i,j,k,l).lt.1E-5) f(i,j,k,l)=1E-5 IF (L.LT.UPA(I)) THEN - WEIGHT=F2(S,I,J,K,L)*WE(K)*WMU(L) + WEIGHT=F2(S,I,J,K,L)*WE(S,K)*WMU(L) IF (MLT(J).LE.6.OR.MLT(J).GE.18.) THEN XNN(S,I)=XNN(S,I)+WEIGHT - ENERN(S,I)=ENERN(S,I)+EKEV(K)*WEIGHT + ENERN(S,I)=ENERN(S,I)+EKEV(S,K)*WEIGHT ELSE XND(S,I)=XND(S,I)+WEIGHT - ENERD(S,I)=ENERD(S,I)+EKEV(K)*WEIGHT + ENERD(S,I)=ENERD(S,I)+EKEV(S,K)*WEIGHT ENDIF ENDIF END DO @@ -1117,7 +1102,7 @@ subroutine ram_hour_write(S) !Previously WRESULT in ram_all WRITE(UNITTMP_,32) StringDate,LZ(I),KP,MLT(J), species(S)%s_code if (outsideMGNP(I,J) == 1) F(I,J,:,:) = 1e-31 DO 27 K=4,NE-1 -27 WRITE(UNITTMP_,30) EKEV(K),(F(I,J,K,L),L=2,NPA-2) +27 WRITE(UNITTMP_,30) EKEV(S,K),(F(I,J,K,L),L=2,NPA-2) 25 CONTINUE END DO close(UNITTMP_) @@ -1140,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 @@ -1163,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 @@ -1313,7 +1298,7 @@ subroutine write_dsbnd(S) OPEN(UNIT=UNITTMP_,FILE=NameFluxFile, STATUS='UNKNOWN') WRITE(UNITTMP_,*)'EKEV FGEOSB [1/cm2/s/ster/keV] T=',TimeRamElapsed/3600,Kp,F107 DO K=2,NE - WRITE(UNITTMP_,*) EKEV(K),(FGEOS(S,J,K,2)/FFACTOR(S,NR,K,2),J=1,NT) + WRITE(UNITTMP_,*) EKEV(S,K),(FGEOS(S,J,K,2)/FFACTOR(S,NR,K,2),J=1,NT) END DO CLOSE(UNITTMP_) diff --git a/src/ModRamInit.f90 b/src/ModRamInit.f90 index 5923f5f..7ba8eb4 100644 --- a/src/ModRamInit.f90 +++ b/src/ModRamInit.f90 @@ -95,7 +95,7 @@ subroutine ram_allocate ! ModRamInit Variables ALLOCATE(RMAS(NS), V(NS,NE), VBND(NS,NE), GREL(NS,NE), GRBND(NS,NE), FACGR(NS,NE), & - EPP(NS,NE), ERNH(NS,NE), UPA(NR), WE(NE), DE(NE), EKEV(NE), EBND(NE), & + EPP(NS,NE), ERNH(NS,NE), UPA(NR), WE(NS,NE), DE(NS,NE), EKEV(NS,NE), EBND(NS,NE), & PHI(NT), LT(NT), MLT(NT), MU(NPA), DMU(NPA), WMU(NPA), PAbn(NPA), LZ(NR+1), & RLZ(NR+1), AMLA(Slen), BE(NR+1,Slen), GridExtend(NRExtend), ZRPabn(NR,NPA,Slen), & FFACTOR(NS,NR,NE,NPA), PA(NPA)) @@ -419,44 +419,49 @@ SUBROUTINE ARRAYS(S) PHI(J)=(J-1)*DPHI ! Magnetic local time in radian MLT(J)=PHI(J)*12./PI ! Magnetic local time in hour END DO - IP1=(MLT(2)-MLT(1))/0.5 + IP=(MLT(2)-MLT(1))/0.5 - DO iS=1,nS - RMAS(iS)=MP*species(iS)%s_mass ! rest mass of each species (kg) - END DO + RMAS(S)=MP*species(S)%s_mass ! rest mass of each species (kg) ! Calculate Kinetic Energy EKEV [keV] at cent, RW depends on NE - ELB=EnergyMin ! Lower limit of energy in keV - IF (abs(ELB-0.01).le.1e-9) THEN - WE(1)=2.8E-3 ! |_._|___.___|____.____|______.______| - RW=1.36 ! . < DE > < WE > - END IF ! EKEV EBND - IF (abs(ELB-0.1).le.1e-9) THEN ! relativistic - WE(1)=3E-2 - RW=1.27 - END IF - IF (abs(ELB-1.0).le.1e-9) THEN - WE(1)=0.31 - RW=1.16 - END IF - - EKEV(1)=ELB+0.5*WE(1) - GREL(S,1)=1.+EKEV(1)*1000.*Q/RMAS(S)/CS/CS + ELB=species(S)%Emin + if (species(S)%Emax < 0.0) then + RW = 1.27 + WE(S,1) = 3E-2 + else + RW = EXP(LOG(ELB) + (LOG(species(S)%Emax) - LOG(ELB))/nE)/ELB + WE(S,1) = ELB*RW - ELB + endif + !IF (abs(ELB-0.01).le.1e-9) THEN + ! WE(1)=2.8E-3 ! |_._|___.___|____.____|______.______| + ! RW=1.36 ! . < DE > < WE > + !END IF ! EKEV EBND + !IF (abs(ELB-0.1).le.1e-9) THEN ! relativistic + ! WE(1)=3E-2 + ! RW=1.27 + !END IF + !IF (abs(ELB-1.0).le.1e-9) THEN + ! WE(1)=0.31 + ! RW=1.16 + !END IF + + EKEV(S,1)=ELB+0.5*WE(S,1) + GREL(S,1)=1.+EKEV(S,1)*1000.*Q/RMAS(S)/CS/CS V(S,1)=CS*SQRT(GREL(S,1)**2-1.)/GREL(S,1) - EBND(1)=ELB+WE(1) - GRBND(S,1)=1.+EBND(1)*1000.*Q/RMAS(S)/CS/CS + EBND(S,1)=ELB+WE(S,1) + GRBND(S,1)=1.+EBND(S,1)*1000.*Q/RMAS(S)/CS/CS VBND(S,1)=CS*SQRT(GRBND(S,1)**2-1.)/GRBND(S,1) DO K=1,NE-1 - WE(K+1)=WE(K)*RW ! WE(K) [keV] is a power series - EBND(K+1)=EBND(K)+WE(K+1) ! E[keV] at bound of grid - DE(K)=0.5*(WE(K)+WE(K+1)) - EKEV(K+1)=EKEV(K)+DE(K) ! E[keV] at cent of grid - GREL(S,K+1)=1.+EKEV(K+1)*1000.*Q/RMAS(S)/CS/CS + WE(S,K+1)=WE(S,K)*RW ! WE(K) [keV] is a power series + EBND(S,K+1)=EBND(S,K)+WE(S,K+1) ! E[keV] at bound of grid + DE(S,K)=0.5*(WE(S,K)+WE(S,K+1)) + EKEV(S,K+1)=EKEV(S,K)+DE(S,K) ! E[keV] at cent of grid + GREL(S,K+1)=1.+EKEV(S,K+1)*1000.*Q/RMAS(S)/CS/CS V(S,K+1)=CS*SQRT(GREL(S,K+1)**2-1.)/GREL(S,K+1) ! Veloc [m/s] at cent - GRBND(S,K+1)=1.+EBND(K+1)*1000.*Q/RMAS(S)/CS/CS + GRBND(S,K+1)=1.+EBND(S,K+1)*1000.*Q/RMAS(S)/CS/CS VBND(S,K+1)=CS*SQRT(GRBND(S,K+1)**2-1.)/GRBND(S,K+1) ! Veloc [m/s] at bound END DO - DE(NE)=0.5*WE(NE)*(1.+RW) + DE(S,NE)=0.5*WE(S,NE)*(1.+RW) ! CONE - pitch angle loss cone in degree ! Dipole loss cone. This needs to be changed in future version to support a @@ -571,8 +576,8 @@ SUBROUTINE ARRAYS(S) ! Energy factors used in RAM equations DO K=1,NE - ERNH(S,K)=WE(K)*GREL(S,K)/SQRT((GREL(S,K)-1.)*(GREL(S,K)+1.)) ! [1/cm3] - EPP(S,K)=ERNH(S,K)*EKEV(K) + ERNH(S,K)=WE(S,K)*GREL(S,K)/SQRT((GREL(S,K)-1.)*(GREL(S,K)+1.)) ! [1/cm3] + EPP(S,K)=ERNH(S,K)*EKEV(S,K) FACGR(S,K)=GREL(S,K)*SQRT((GREL(S,K)-1.)*(GREL(S,K)+1.)) END DO @@ -694,6 +699,7 @@ SUBROUTINE init_input .and.(trim(species(i)%Initialization).ne.'na')) then call load_injection_file(i, trim(species(i)%Initialization), F2(i,:,:,:,:)) endif + call ANISCH(i) enddo ! Reset dipole initialization for SWMF runs diff --git a/src/ModRamInjection.f90 b/src/ModRamInjection.f90 index ce5c2fc..84fe608 100644 --- a/src/ModRamInjection.f90 +++ b/src/ModRamInjection.f90 @@ -78,7 +78,7 @@ subroutine load_injection_file(S, fname, fluxOut) iE = 1 do i = 2,nE - if ((energy > EkeV(i-1)).and.(energy < EkeV(i))) then + if ((energy > EkeV(S,i-1)).and.(energy < EkeV(S,i))) then iE = i exit endif @@ -93,8 +93,8 @@ subroutine load_injection_file(S, fname, fluxOut) enddo ! Add the number density to the grid assuming a maxwellian distribution in energy - factor = 4.0E6*EkeV(iE) - den = 1._dp/sqrt(species(S)%s_mass)*exp(-1._dp*EkeV(iE)/energy)*factor*den*energy**(-1.5) + factor = 4.0E6*EkeV(S,iE) + den = 1._dp/sqrt(species(S)%s_mass)*exp(-1._dp*EkeV(S,iE)/energy)*factor*den*energy**(-1.5) !den = FFactor(S,iR,iE,iPa)*FNHS(iR,iT,iPa)*den/dE(iE)/dMu(iPa) num_density(iR,iT,iE,iPa) = num_density(iR,iT,iE,iPa) + den*FFactor(S,iR,iE,iPa)*FNHS(iR,iT,iPa) enddo read_file diff --git a/src/ModRamLoss.f90 b/src/ModRamLoss.f90 index 3524f32..2692b58 100644 --- a/src/ModRamLoss.f90 +++ b/src/ModRamLoss.f90 @@ -41,7 +41,7 @@ SUBROUTINE CEPARA(S) DO K=2,NE DO I=2,NR DO J=1,NT - X=LOG10(EKEV(K)) + X=LOG10(EKEV(S,K)) IF (X.LT.-2.) X=-2. Y=-18.767-0.11017*X-3.8173e-2*X**2-0.1232*X**3-5.0488e-2*X**4 ALPHA=10.**Y*V(S,K)*HDNS(I,J,L)*DTs @@ -56,7 +56,7 @@ SUBROUTINE CEPARA(S) DO K=2,NE DO I=2,NR DO J=1,NT - X=LOG10(EKEV(K)) + X=LOG10(EKEV(S,K)) IF(X.LT.-2.) X=-2. Y=-20.789+0.92316*X-0.68017*X**2+0.66153*X**3-0.20998*X**4 ALPHA=10.**Y*V(S,K)*HDNS(I,J,L)*DTs @@ -71,7 +71,7 @@ SUBROUTINE CEPARA(S) DO K=2,NE DO I=2,NR DO J=1,NT - X=LOG10(EKEV(K)) + X=LOG10(EKEV(S,K)) IF(X.LT.-2.) X=-2. Y=-18.987-0.10613*X-5.4841E-3*X**2-1.6262E-2*X**3-7.0554E-3*X**4 ALPHA=10.**Y*V(S,K)*HDNS(I,J,L)*DTs @@ -438,7 +438,7 @@ subroutine PARA_FLC(S) if(DoWriteFLCDiffCoeff)then do l=1, nPA-1 - write(unittmp_, '(2(2x, f11.3),5(2x,e11.3))') EkeV(k), acos(MU(l)+0.5*WMU(l)), r_gyro(i,j,k), & + write(unittmp_, '(2(2x, f11.3),5(2x,e11.3))') EkeV(S,k), acos(MU(l)+0.5*WMU(l)), r_gyro(i,j,k), & r_curvEq(i,j)/REarth, epsil(i,j,k),tau_bounce(l),FLC_coef(S,i,j,k,l) end do end if diff --git a/src/ModRamRestart.f90 b/src/ModRamRestart.f90 old mode 100755 new mode 100644 index b10f690..4e984ac --- a/src/ModRamRestart.f90 +++ b/src/ModRamRestart.f90 @@ -6,9 +6,8 @@ module ModRamRestart implicit none - - contains - !========================================================================== + contains +!========================================================================== subroutine write_restart use ModRamParams, ONLY: TimedRestarts !!!! Module Variables @@ -16,10 +15,10 @@ subroutine write_restart use ModRamFunctions, ONLY: RamFileName use ModRamTiming, ONLY: TimeRamElapsed, TimeRamStart, TimeRamNow, DtsNext, & TOld - use ModRamGrids, ONLY: NR, NT, NE, NPA + use ModRamGrids, ONLY: NR, NT, NE, NPA, nS use ModRamVariables, ONLY: F2, PParT, PPerT, FNHS, FNIS, BOUNHS, BOUNIS, & BNES, HDNS, EIR, EIP, dBdt, dIdt, dIbndt, VTN, & - VTOL, VT, EIR, EIP, EKEV, PA, NECR, TAU + VTOL, VT, EIR, EIP, EKEV, PA, NECR, TAU, species use ModScbGrids, ONLY: nthe, npsi, nzeta use ModScbParams, ONLY: constZ, constTheta use ModScbVariables, ONLY: x, y, z, alphaVal, psiVal, chiVal, xpsiout, & @@ -34,8 +33,7 @@ subroutine write_restart implicit none - integer :: iFluxEVar, iFluxHVar, iFluxHeVar, iFluxOVar, iPParTVar, & - iPPerTVar, iHVar, iBHVar, iEGridVar, iPaGridVar, & + integer :: iFluxVar, iPParTVar, iPPerTVar, iHVar, iBHVar, iEGridVar, iPaGridVar, & iIVar, iBIVar, iBNESVar, iHDNSVar, iEIRVar, iEIPVar, & iDtVar, iXVar, iYVar, iZVar, iBxVar, iByVar, iBzVar, & iVTNVar, iAValVar, iBValVar, iVTOLVar, & @@ -45,7 +43,7 @@ subroutine write_restart iNECRVar, itauVar integer :: nRDim, nTDim, nEDim, nPaDim, nSDim, nThetaDim, nPsiDim, & - nZetaDim + nZetaDim, iS integer, parameter :: iDeflate = 2, yDeflate = 1 character(len=999) :: NameFile @@ -97,44 +95,32 @@ subroutine write_restart ! START DEFINE MODE !! TEMP STUFF - iStatus = nf90_def_var(iFileID, 'EnergyGrid', nf90_double, & - (/nEDim/), iEGridVar) + do iS = 1, nS + iStatus = nf90_def_var(iFileID, 'EnergyGrid_'//species(iS)%s_name, nf90_double, & + (/nEDim/), iEgridVar) + enddo iStatus = nf90_def_var(iFileID, 'PitchAngleGrid', nf90_double, & (/nPaDim/), iPaGridVar) !! FLUXES - !!! Electron Flux - iStatus = nf90_def_var(iFileID, 'FluxE', nf90_double, & - (/nRdim,nTDim,nEDim,nPaDim/), iFluxEVar) - iStatus = nf90_def_var_deflate(iFileID, iFluxEVar, 0, yDeflate, iDeflate) - - iStatus = nf90_put_att(iFileID, iFluxEVar, 'title', & - 'This is an example title') - - !!! Proton Flux - iStatus = nf90_def_var(iFileID, 'FluxH', nf90_double, & - (/nRdim,nTDim,nEDim,nPaDim/), iFluxHVar) - iStatus = nf90_def_var_deflate(iFileID, iFluxHVar, 0, yDeflate, iDeflate) - - !!! Oxygen Flux - iStatus = nf90_def_var(iFileID, 'FluxO', nf90_double, & - (/nRdim,nTDim,nEDim,nPaDim/), iFluxOVar) - iStatus = nf90_def_var_deflate(iFileID, iFluxOVar, 0, yDeflate, iDeflate) - - !!! Helium Flux - iStatus = nf90_def_var(iFileID, 'FluxHe', nf90_double, & - (/nRdim,nTDim,nEDim,nPaDim/), iFluxHeVar) - iStatus = nf90_def_var_deflate(iFileID, iFluxHeVar, 0, yDeflate, iDeflate) + do iS = 1, nS + iStatus = nf90_def_var(iFileID, 'Flux_'//species(iS)%s_name, nf90_double, & + (/nRDim, nTDim, nEDim, nPaDim/), iFluxVar) + iStatus = nf90_def_var_deflate(iFileID, iFluxVar, 0, yDeflate, iDeflate) + enddo !! PRESSURES !!! Total Parallel Pressures - iStatus = nf90_def_var(iFileID, 'PParT', nf90_double, & - (/nSdim,nRDim,nTDim/), iPParTVar) - iStatus = nf90_def_var_deflate(iFileID, iPParTVar, 0, yDeflate, iDeflate) - - !!! Total Perpendicular Pressures - iStatus = nf90_def_var(iFileID, 'PPerT', nf90_double, & - (/nSdim,nRDim,nTDim/), iPPerTVar) - iStatus = nf90_def_var_deflate(iFileID, iPPerTVar, 0, yDeflate, iDeflate) + do iS = 1, nS + !!! Total Parallel Pressures + iStatus = nf90_def_var(iFileID, 'PParT_'//species(iS)%s_name, nf90_double, & + (/nRDim,nTDim/), iPParTVar) + iStatus = nf90_def_var_deflate(iFileID, iPParTVar, 0, yDeflate, iDeflate) + + !!! Total Perpendicular Pressures + iStatus = nf90_def_var(iFileID, 'PPerT_'//species(iS)%s_name, nf90_double, & + (/nRDim,nTDim/), iPPerTVar) + iStatus = nf90_def_var_deflate(iFileID, iPPerTVar, 0, yDeflate, iDeflate) + enddo !! MAGNETIC FIELD !! Bx @@ -309,18 +295,25 @@ subroutine write_restart ! START WRITE MODE !! TEMP STUFF - iStatus = nf90_put_var(iFileID, iEGridVar, EKEV(:)) + do iS = 1, nS + iStatus = nf90_inq_varid(iFileID, 'EnergyGrid_'//species(iS)%s_name, iEGridVar) + iStatus = nf90_put_var(iFileID, iEGridVar, EKEV(iS,:)) + enddo iStatus = nf90_put_var(iFileID, iPaGridVar, PA(:)) !! FLUXES - iStatus = nf90_put_var(iFileID, iFluxEVar, F2(1,:,:,:,:)) - iStatus = nf90_put_var(iFileID, iFluxHVar, F2(2,:,:,:,:)) - iStatus = nf90_put_var(iFileID, iFluxHeVar, F2(3,:,:,:,:)) - iStatus = nf90_put_var(iFileID, iFluxOVar, F2(4,:,:,:,:)) + do iS = 1, nS + iStatus = nf90_inq_varid(iFileID, 'Flux_'//species(iS)%s_name, iFluxVar) + iStatus = nf90_put_var(iFileID, iFluxVar, F2(iS,:,:,:,:)) + enddo !! PRESSURES - iStatus = nf90_put_var(iFileID, iPParTVar, PParT(:,:,:)) - iStatus = nf90_put_var(iFileID, iPPerTVar, PPerT(:,:,:)) + do iS = 1, nS + iStatus = nf90_inq_varid(iFileID, 'PParT_'//species(iS)%s_name, iPParTVar) + iStatus = nf90_put_var(iFileID, iPParTVar, PParT(iS,:,:)) + iStatus = nf90_inq_varid(iFileID, 'PPerT_'//species(iS)%s_name, iPPerTVar) + iStatus = nf90_put_var(iFileID, iPPerTVar, PPerT(iS,:,:)) + enddo !! MAGNETIC FIELD iStatus = nf90_put_var(iFileID, iBxVar, bX(:,:,:)) @@ -386,14 +379,14 @@ end subroutine write_restart subroutine read_restart !!!! Module Variables use ModRamMain, ONLY: PathRestartIn - use ModRamGrids, ONLY: nPa, nT, nR + use ModRamGrids, ONLY: nPa, nT, nR, nS, nE use ModRamFunctions, ONLY: RamFileName use ModRamTiming, ONLY: DtsNext, TOld use ModRamVariables, ONLY: F2, PParT, PPerT, FNHS, FNIS, BOUNHS, BOUNIS, & BNES, HDNS, EIR, EIP, dBdt, dIdt, dIbndt, VTN, & VTOL, VT, EIR, EIP, PParH, PPerH, PParO, PAbn, & PPerO, PParHe, PPerHe, PParE, PPerE, LZ, MU, & - NECR, tau + NECR, tau, EKEV, species use ModScbGrids, ONLY: npsi, nzeta use ModScbParams, ONLY: constTheta, constZ use ModScbVariables, ONLY: x, y, z, alphaVal, psiVal, chiVal, xpsiout, & @@ -408,8 +401,8 @@ subroutine read_restart use nrtype, ONLY: DP, pi_d, twopi_d implicit none - integer :: j, L - integer :: iFluxEVar, iFluxHVar, iFluxHeVar, iFluxOVar, iPParTVar, & + integer :: j, L, iS + integer :: iFluxVar, iPParTVar, iEGridVar, & iPPerTVar, iHVar, iBHVar, iBxVar, iByVar, iBzVar, & iIVar, iBIVar, iBNESVar, iHDNSVar, iEIRVar, iEIPVar, & iDtVar, iXVar, iYVar, iZVar, & @@ -433,16 +426,6 @@ subroutine read_restart call ncdf_check(iStatus, NameSub) ! GET VARIABLE IDS - !! FLUXES - iStatus = nf90_inq_varid(iFileID, 'FluxE', iFluxEVar) - iStatus = nf90_inq_varid(iFileID, 'FluxH', iFluxHVar) - iStatus = nf90_inq_varid(iFileID, 'FluxHe', iFluxHeVar) - iStatus = nf90_inq_varid(iFileID, 'FluxO', iFluxOVar) - - !! PRESSURES - iStatus = nf90_inq_varid(iFileID, 'PParT', iPParTVar) - iStatus = nf90_inq_varid(iFileID, 'PPerT', iPPerTVar) - !! MAGNETIC FIELD iStatus = nf90_inq_varid(iFileID, 'Bx', iBxVar) iStatus = nf90_inq_varid(iFileID, 'By', iByVar) @@ -494,24 +477,41 @@ subroutine read_restart iStatus = nf90_inq_varid(iFileID, 'necr', iNECRVar) iStatus = nf90_inq_varid(iFileID, 'tau', itauVar) - ! READ DATA - !! FLUXES - iStatus = nf90_get_var(iFileID, iFluxEVar, F2(1,:,:,:,:)) - iStatus = nf90_get_var(iFileID, iFluxHVar, F2(2,:,:,:,:)) - iStatus = nf90_get_var(iFileID, iFluxHeVar, F2(3,:,:,:,:)) - iStatus = nf90_get_var(iFileID, iFluxOVar, F2(4,:,:,:,:)) + ! Check which version we are reading + iStatus = nf90_inq_varid(iFileID, 'FluxE', iFluxVar) + if (iStatus == nf90_noerr) then + iStatus = nf90_inq_varid(iFileID, 'FluxE', iFluxVar) + iStatus = nf90_get_var(iFileID, iFluxVar, F2(1,:,:,:,:)) - !! PRESSURES - iStatus = nf90_get_var(iFileID, iPParTVar, PParT(:,:,:)) - iStatus = nf90_get_var(iFileID, iPPerTVar, PPerT(:,:,:)) - PPerO = PPerT(4,:,:) - PParO = PParT(4,:,:) - PPerHe = PPerT(3,:,:) - PParHe = PParT(3,:,:) - PPerE = PPerT(1,:,:) - PParE = PParT(1,:,:) - PPerH = PPerT(2,:,:) - PParH = PParT(2,:,:) + iStatus = nf90_inq_varid(iFileID, 'FluxH', iFluxVar) + iStatus = nf90_get_var(iFileID, iFluxVar, F2(2,:,:,:,:)) + + iStatus = nf90_inq_varid(iFileID, 'FluxHe', iFluxVar) + iStatus = nf90_get_var(iFileID, iFluxVar, F2(3,:,:,:,:)) + + iStatus = nf90_inq_varid(iFileID, 'FluxO', iFluxVar) + iStatus = nf90_get_var(iFileID, iFluxVar, F2(4,:,:,:,:)) + + iStatus = nf90_inq_varid(iFileID, 'PParT', iPParTVar) + iStatus = nf90_get_var(iFileID, iPParTVar, PParT(:,:,:)) + + iStatus = nf90_inq_varid(iFileID, 'PPerT', iPPerTVar) + iStatus = nf90_get_var(iFileID, iPPerTVar, PPerT(:,:,:)) + else + do iS = 1, nS + iStatus = nf90_inq_varid(iFileID, 'EnergyGrid_'//species(iS)%s_name, iEGridVar) + iStatus = nf90_get_var(iFileID, iEGridVar, EKEV(iS,:)) + + iStatus = nf90_inq_varid(iFileID, 'Flux_'//species(iS)%s_name, iFluxVar) + iStatus = nf90_get_var(iFileID, iFluxVar, F2(iS,:,:,:,:)) + + iStatus = nf90_inq_varid(iFileID, 'PParT_'//species(iS)%s_name, iPParTVar) + iStatus = nf90_get_var(iFileID, iPParTVar, PParT(iS,:,:)) + + iStatus = nf90_inq_varid(iFileID, 'PPerT_'//species(iS)%s_name, iPPerTVar) + iStatus = nf90_get_var(iFileID, iPPerTVar, PPerT(iS,:,:)) + enddo + endif !! MAGNETIC FIELD iStatus = nf90_get_var(iFileID, iBxVar, bX(:,:,:)) diff --git a/src/ModRamRun.f90 b/src/ModRamRun.f90 index f9811d6..6b7d8b5 100644 --- a/src/ModRamRun.f90 +++ b/src/ModRamRun.f90 @@ -247,8 +247,8 @@ SUBROUTINE SUMRC(S) DO K=2,NE DO L=2,NPA DO J=1,NT-1 - WEIGHT=F2(S,I,J,K,L)*WE(K)*WMU(L) - SETRC(S)=SETRC(S)+EKEV(K)*WEIGHT + WEIGHT=F2(S,I,J,K,L)*WE(S,K)*WMU(L) + SETRC(S)=SETRC(S)+EKEV(S,K)*WEIGHT END DO END DO END DO @@ -379,7 +379,7 @@ SUBROUTINE ANISCH(S) PPER=PPER+EPP(S,K)*SUME PPAR=PPAR+EPP(S,K)*SUMA RNHT=RNHT+ERNH(S,K)*SUMN - EDEN=EDEN+ERNH(S,K)*EKEV(K)*SUMN + EDEN=EDEN+ERNH(S,K)*EKEV(S,K)*SUMN ENDDO ANIS=PPER/2./PPAR-1. EPAR=2*PPAR/RNHT @@ -497,7 +497,7 @@ SUBROUTINE ANISCH(S) ENDDO ENDDO DO K=2,NE - ER1=LOG10(EKEV(K)) + ER1=LOG10(EKEV(S, K)) CALL GSL_Interpolation_2D(ALENOR,fpofc,DUMP,ER1,xfrl,Y,GSLerr) DWAVE(L)=10.**Y*fnorm/GREL(S,K)**2*(1.-MUBOUN*MUBOUN)/MUBOUN ! denorm pa [1/s] ATAW(I,J,K,L)=DWAVE(L) ! call WPADIF twice, implicit @@ -567,7 +567,7 @@ SUBROUTINE ANISCH(S) END DO END DO DO K=2,NE - ER1 = log10(EKEV(K)) + ER1 = log10(EKEV(S,K)) ! EkeV_emic: 0.1 keV to 1000 keV call GSL_Interpolation_2D(logEkeV_emic, fp2c_emic, & DUMP2_h, ER1, xfrl, Y, GSLerr) @@ -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 diff --git a/src/ModRamSats.f90 b/src/ModRamSats.f90 index aaa182f..9ce88e5 100644 --- a/src/ModRamSats.f90 +++ b/src/ModRamSats.f90 @@ -406,20 +406,22 @@ subroutine create_sat_file(FileNameIn) iStatus = nf90_put_att(iFileID, iOFluxVar, 'units', '1/cm2/s/keV') enddo - ! ENERGY GRID - iStatus = nf90_def_var(iFileID, 'energy_grid', nf90_float, & - iEnDim, iEgridVar) - iStatus = nf90_put_att(iFileID, iEgridVar, 'title', & - 'Energy grid as values at window centers.') - iStatus = nf90_put_att(iFileID, iEgridVar, 'units', 'KeV') - - iStatus = nf90_def_var(iFileID, 'energy_width', nf90_float, & - iEnDim, iEwidVar) - iStatus = nf90_put_att(iFileID, iEwidVar, 'title', & - 'Widths of each energy bin whose centers are listed in energy_grid.') - iStatus = nf90_put_att(iFileID, iEwidVar, 'units', 'KeV') - - ! PITCH-ANGLE GRID + ! ENERGY GRID + do iS = 1, nS + iStatus = nf90_def_var(iFileID, 'energy_grid_'//species(iS)%s_code, nf90_float, & + iEnDim, iEgridVar) + iStatus = nf90_put_att(iFileID, iEgridVar, 'title', & + 'Energy grid as values at window centers.') + iStatus = nf90_put_att(iFileID, iEgridVar, 'units', 'KeV') + + iStatus = nf90_def_var(iFileID, 'energy_width_'//species(iS)%s_code, nf90_float, & + iEnDim, iEwidVar) + iStatus = nf90_put_att(iFileID, iEwidVar, 'title', & + 'Widths of each energy bin whose centers are listed in energy_grid.') + iStatus = nf90_put_att(iFileID, iEwidVar, 'units', 'KeV') + enddo + + ! PITCH-ANGLE GRID iStatus = nf90_def_var(iFileID, 'pa_grid', nf90_float, & iPaDim, iPgridVar) iStatus = nf90_put_att(iFileID, iPgridVar, 'title', & @@ -439,8 +441,13 @@ subroutine create_sat_file(FileNameIn) iStatus = nf90_enddef(iFileID) ! Write static (no time dimension) variables. - iStatus = nf90_put_var(iFileID, iEgridVar, Ekev) - iStatus = nf90_put_var(iFileID, iEwidVar, wE) + do iS = 1, nS + iStatus = nf90_inq_varid(iFileID, 'energy_grid_'//species(iS)%s_code, iEgridVar) + iStatus = nf90_put_var(iFileID, iEgridVar, Ekev(iS,:)) + + iStatus = nf90_inq_varid(iFileID, 'energy_width_'//species(iS)%s_code, iEwidVar) + iStatus = nf90_put_var(iFileID, iEwidVar, wE(iS,:)) + enddo iStatus = nf90_put_var(iFileID, iPgridVar, Mu) iStatus = nf90_put_var(iFileID, iPwidVar, wMu) iStatus = nf90_put_var(iFileID, iDtVar, DtWriteSat) diff --git a/src/ModRamSce.f90 b/src/ModRamSce.f90 index 36345eb..47e2a83 100644 --- a/src/ModRamSce.f90 +++ b/src/ModRamSce.f90 @@ -119,7 +119,7 @@ subroutine calculate_precip_flux_jr(IS, nTheta, nPhi, rIono, energy_flux, ave_e, ave_flux(i,j,k) = ave_flux(i,j,k)+f(i,j,k,l)*WMU(l) end do ave_flux(i,j,k) = ave_flux(i,j,k)/(mu(NPA)-mu(UPA(i))) - num_fluxeq(i,j) = num_fluxeq(i,j) + pi_d*ave_flux(i,j,k)*WE(k) + num_fluxeq(i,j) = num_fluxeq(i,j) + pi_d*ave_flux(i,j,k)*WE(IS,k) end do end do end do @@ -228,18 +228,22 @@ subroutine calculate_precip_flux_jr(IS, nTheta, nPhi, rIono, energy_flux, ave_e, energy_flux_iono(i,j) = 0.0 num_flux_iono(i,j) = 0.0 ave_e_iono(i,j) = 0.0 + do k=1, nE + energy_flux_iono(i,j) = energy_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*EKEV(IS,k)*WE(IS,k) + num_flux_iono(i,j) = num_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*WE(IS,k) + end do if(conductance_model .eq. 9)then ! GLOW conductance do k=1,nE - if (EKEV(k) .ge. 0.5e-3 .and. EKEV(k) .le. 46)then ! GLOW energy range: 0.5 eV to 46 keV - energy_flux_iono(i,j) = energy_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*EKEV(k)*WE(k) - num_flux_iono(i,j) = num_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*WE(k) + if (EKEV(iS,k) .ge. 0.5e-3 .and. EKEV(iS,k) .le. 46)then ! GLOW energy range: 0.5 eV to 46 keV + energy_flux_iono(i,j) = energy_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*EKEV(iS,k)*WE(iS,k) + num_flux_iono(i,j) = num_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*WE(iS,k) end if end do else do k=1, nE - if (EKEV(k) .ge. 0.5 .and. EKEV(k) .le. 50)then ! Robinson energy range: 500eV to 50 keV - energy_flux_iono(i,j) = energy_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*EKEV(k)*WE(k) - num_flux_iono(i,j) = num_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*WE(k) + if (EKEV(iS,k) .ge. 0.5 .and. EKEV(iS,k) .le. 50)then ! Robinson energy range: 500eV to 50 keV + energy_flux_iono(i,j) = energy_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*EKEV(iS,k)*WE(iS,k) + num_flux_iono(i,j) = num_flux_iono(i,j) + pi_d*ave_fluxEQ(i,j,k)*WE(iS,k) end if end do end if @@ -264,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 diff --git a/src/ModRamSpecies.f90 b/src/ModRamSpecies.f90 index 6974cea..c64d6e2 100644 --- a/src/ModRamSpecies.f90 +++ b/src/ModRamSpecies.f90 @@ -27,9 +27,10 @@ MODULE ModRamSpecies real(DP), allocatable :: CEX_nO(:) real(DP), allocatable :: CEX_nN(:) real(DP) :: plasmasphereRatio + real(DP) :: Emin, Emax ! Min and Max energy range of species end type SpeciesType - integer, parameter :: nSpecies = 6 + integer, parameter :: nSpecies = 7 type(SpeciesType), dimension(nSpecies) :: RAMSpecies contains @@ -54,6 +55,8 @@ subroutine DefineSpecies RAMSpecies(1)%CEX_species = 'na' RAMSpecies(1)%Initialization = 'InitializationFile' RAMSpecies(1)%plasmasphereRatio = 1.0 + RAMSpecies(1)%Emin = 0.1 + RAMSpecies(1)%Emax = -1.0 ! Protons RAMSpecies(2)%s_name = "Hydrogen" @@ -70,6 +73,8 @@ subroutine DefineSpecies RAMSpecies(2)%CEX_species = 'na' RAMSpecies(2)%Initialization = 'InitializationFile' RAMSpecies(2)%plasmasphereRatio = 0.77 + RAMSpecies(2)%Emin = 0.1 + RAMSpecies(2)%Emax = -1.0 ! Helium +1 RAMSpecies(3)%s_name = "HeliumP1" @@ -86,6 +91,8 @@ subroutine DefineSpecies RAMSpecies(3)%CEX_species = 'na' RAMSpecies(3)%Initialization = 'InitializationFile' RAMSpecies(3)%plasmasphereRatio = 0.2 + RAMSpecies(3)%Emin = 0.1 + RAMSpecies(3)%Emax = -1.0 ! Oxygen +1 RAMSpecies(4)%s_name = "OxygenP1" @@ -102,6 +109,8 @@ subroutine DefineSpecies RAMSpecies(4)%CEX_species = 'na' RAMSpecies(4)%Initialization = 'InitializationFile' RAMSpecies(4)%plasmasphereRatio = 0.03 + RAMSpecies(4)%Emin = 0.1 + RAMSpecies(4)%Emax = -1.0 ! Nitrogen +1 RAMSpecies(5)%s_name = "Nitrogen" @@ -118,6 +127,8 @@ subroutine DefineSpecies RAMSpecies(5)%CEX_species = 'nH' RAMSpecies(5)%Initialization = 'na' RAMSpecies(5)%plasmasphereRatio = 0.0 + RAMSpecies(5)%Emin = 0.1 + RAMSpecies(5)%Emax = -1.0 ! Strontium +1 RAMSpecies(6)%s_name = "Strontium" @@ -134,6 +145,26 @@ subroutine DefineSpecies RAMSpecies(6)%CEX_species = 'nH nO nN' RAMSpecies(6)%Initialization = 'StrontiumPlusOneIons.dat' RAMSpecies(6)%plasmasphereRatio = 0.0 + RAMSpecies(6)%Emin = 0.1 + RAMSpecies(6)%Emax = 500.0 + + ! High Energy Electrons + RAMSpecies(7)%s_name = "RadBeltE" + RAMSpecies(7)%s_code = "Re" + RAMSpecies(7)%s_mass = 5.4462E-4 + RAMSpecies(7)%s_comp = 1.0 + RAMSpecies(7)%s_charge = -1 + RAMSpecies(7)%SCB = .false. + RAMSpecies(7)%WPI = .true. + RAMSpecies(7)%CEX = .false. + RAMSpecies(7)%FLC = .false. + RAMSpecies(7)%EMIC = .false. + RAMSpecies(7)%CEX_file = 'na' + RAMSpecies(7)%CEX_species = 'na' + RAMSpecies(7)%Initialization = 'InitializationFile' + RAMSpecies(7)%plasmasphereRatio = 1.0 + RAMSpecies(7)%Emin = 100.0 + RAMSpecies(7)%Emax = 10000.0 end subroutine DefineSpecies diff --git a/src/ModRamVariables.f90 b/src/ModRamVariables.f90 index 27b67a9..b327b51 100644 --- a/src/ModRamVariables.f90 +++ b/src/ModRamVariables.f90 @@ -41,8 +41,8 @@ Module ModRamVariables ! ModRamInit variables real(DP), ALLOCATABLE :: RMAS(:), V(:,:), VBND(:,:), GREL(:,:), GRBND(:,:), & - FACGR(:,:), EPP(:,:), ERNH(:,:), UPA(:), WE(:), DE(:), & - EKEV(:), EBND(:), PHI(:), LT(:), MLT(:), MU(:), DMU(:), & + FACGR(:,:), EPP(:,:), ERNH(:,:), UPA(:), WE(:,:), DE(:,:), & + EKEV(:,:), EBND(:,:), PHI(:), LT(:), MLT(:), MU(:), DMU(:), & WMU(:), PAbn(:), LZ(:), RLZ(:), AMLA(:), BE(:,:), & GridExtend(:), ZRPabn(:,:,:), FFACTOR(:,:,:,:), PA(:) real(DP) :: PHIOFS, IR1, DL1, MDR, IP1, CONF1, CONF2, RFACTOR diff --git a/src/ModRamWPI.f90 b/src/ModRamWPI.f90 index ffb32ef..aba88c9 100644 --- a/src/ModRamWPI.f90 +++ b/src/ModRamWPI.f90 @@ -74,7 +74,7 @@ SUBROUTINE WAVEPARA1(S) ! Calculates the losses due to the w-p interaction DO K=2,NE DO II=2,NR - xE=EKEV(K)/1000. + xE=EKEV(S,K)/1000. xL=LZ(II) if (xL.ge.1.65.and.xL.le.5.0) then do i=8,2,-1 @@ -150,7 +150,7 @@ SUBROUTINE WAVEPARA2(S) DO K=2,NE DO I=2,NR - EMEV=EKEV(K)*0.001 + EMEV=EKEV(S,K)*0.001 R1=0.08*EMEV**(-1.32) R2=0.4*10.**(2.*LZ(I)-6.+0.4*log10(29.*EMEV)) tau_wave=min(R1,R2) @@ -612,16 +612,16 @@ SUBROUTINE WAVELO(S) IF (LZ(I).LE.RLpp(J)) THEN TAU_LIF=WALOS1(I,K)*((10./Bw)**2) ELSEIF (LZ(I).GT.RLpp(J)) THEN - IF (EKEV(K).LE.1000.) THEN + IF (EKEV(S,K).LE.1000.) THEN TAU_LIF=WALOS2(I,K)*(1+WALOS3(I,K)/WALOS2(I,K)) - if (ekev(k).le.1.1) then - tau_lif=tau_lif*37.5813*exp(-1.81255*ekev(k)) - elseif (ekev(k).gt.1.1.and.ekev(k).le.5.) then - tau_lif=tau_lif*(7.5-1.15*ekev(k)) + if (ekev(S,k).le.1.1) then + tau_lif=tau_lif*37.5813*exp(-1.81255*ekev(S,k)) + elseif (ekev(S,k).gt.1.1.and.ekev(S,k).le.5.) then + tau_lif=tau_lif*(7.5-1.15*ekev(S,k)) else tau_lif=tau_lif endif - ELSEIF(EKEV(K).GT.1000.) THEN + ELSEIF(EKEV(S,K).GT.1000.) THEN TAU_LIF=5.*3600*24/KP ENDIF ENDIF diff --git a/src/ModSceIono.f90 b/src/ModSceIono.f90 index ece117a..01939c9 100644 --- a/src/ModSceIono.f90 +++ b/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(:,:), & diff --git a/src/ModSceRun.f90 b/src/ModSceRun.f90 index 8438997..3d2b02b 100644 --- a/src/ModSceRun.f90 +++ b/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, &