More amie types

src/ModInputs.f90

@@ -89,6 +89,10 @@ module ModInputs
   real :: CuspMltHalfWidth = 1.5
   real :: CuspLatHalfWidth = 1.0
 
+  logical :: UsePolarRain = .false.
+  real :: polarRainEFlux = 0.1 ! 0.1 ergs/cm2/s is a small value
+  real :: polarRainAveE = 0.3 ! 300 eV is a nominal polar rain value
+
   logical :: DoOverwriteIonosphere = .false.
   logical :: DoOverwriteWithIRI = .true.
   logical :: DoOverwriteWithSami = .false.

src/aurora.Earth.f90

@@ -75,12 +75,15 @@ subroutine aurora(iBlock)
   endif
 
   if (iProc == 0 .and. AveEFactor /= 1.0 .and. IsFirstTime(iBlock)) then
-    write(*, *) "Auroral Experiments!!!!"
-    write(*, *) "AveEFactor : ", AveEFactor
+    write(*, *) "Reminder / Warning!!!"
+    write(*, *) " --> Auroral Experiment, AveEFactor is set"
+    write(*, *) " --> AveEFactor : ", AveEFactor
+    write(*, *) " --> This factor is currently applied in get_potential!"
   endif
   if (iProc == 0 .and. IsKappaAurora .and. IsFirstTime(iBlock)) then
-    write(*, *) "Auroral Experiments!!!!"
-    write(*, *) "kappa : ", AuroraKappa
+    write(*, *) "Reminder / Warning!!!"
+    write(*, *) " --> Auroral Experiments, AuroraKappa is set"
+    write(*, *) " --> kappa : ", AuroraKappa
   endif
 
   do i = 1, nLats
@@ -96,19 +99,18 @@ subroutine aurora(iBlock)
       HasSomeAurora = .false.
 
       ! For diffuse auroral models (default)
-      if (ElectronEnergyFluxDiffuse(j, i) > 0.1 &
-          .and. ElectronAverageEnergyDiffuse(j, i) > 0.1 &
+      if (ElectronEnergyFluxDiffuse(j, i) > 0.001 &
+          .and. ElectronAverageEnergyDiffuse(j, i) > 0.01 &
           .and. ElectronAverageEnergyDiffuse(j, i) < MaxAveEAurora &
           .and. UseDiffuseAurora &
           ) then
         call do_diffuse_aurora(ElectronEnergyFluxDiffuse(j, i), &
-                               ! aveEFactor is 1 unless set in #auroramods
-                               ElectronAverageEnergyDiffuse(j, i)*AveEFactor, &
+                               ElectronAverageEnergyDiffuse(j, i), &
                                e_diffuse_ED_flux)
         HasSomeAurora = .true.
       endif
 
-      if (IonEnergyFlux(j, i) > 0.1 &
+      if (IonEnergyFlux(j, i) > 0.001 &
           .and. IonAverageEnergy(j, i) > 0.1 &
           .and. IonAverageEnergy(j, i) < MaxAveEAurora &
           .and. UseIonAurora &
@@ -121,7 +123,7 @@ subroutine aurora(iBlock)
 
       ! Monoenergetic aurora
       if (UseMonoAurora &
-          .and. ElectronAverageEnergyMono(j, i) > 1. &
+          .and. ElectronAverageEnergyMono(j, i) > 0.1 &
           .and. ElectronEnergyFluxMono(j, i) > 0.1 &
           .and. ElectronAverageEnergyMono(j, i) < MaxAveEAurora &
           ) then
@@ -149,7 +151,7 @@ subroutine aurora(iBlock)
         ED_Ion_EnergyFlux(n) = i_diffuse_ED_flux(n) ! for consistency
       enddo
 
-      if (maxval(ED_EnergyFlux) > 0.1 .or. maxval(ED_Ion_EnergyFlux) > 0.1) &
+      if (HasSomeAurora) &
         call calc_fang_rates(j, i, iBlock, AuroralBulkIonRate)
 
     enddo
@@ -304,8 +306,8 @@ subroutine do_mono_aurora(eflux_ergs, av_kev, Mono_ED_EnergyFlux)
       HemisphericPowerNorth_mono = HemisphericPowerNorth_mono + power
     endif
 
-    ! Mono is treated as a (skinny) gaussian:
-    if (eflux > 0.1 .and. avee > 0.1) &
+    ! Mono is treated as a (skinny) gaussian, width of 0.1 keV
+    if (avee > 0.0) &
       call calc_gaussian(avee, Mono_ED_EnergyFlux, 0.1)
 
     ! Normalize to E-Flux
@@ -339,8 +341,8 @@ subroutine do_wave_aurora(eflux_ergs, av_kev, wave_EDEnergyFlux)
       HemisphericPowerNorth_wave = HemisphericPowerNorth_wave + power
     endif
 
-    ! Wave is a gaussian, centered at aveE
-    if (eflux > 0.1 .and. avee > 0.1) &
+    ! Wave is a gaussian, centered at aveE, width of 0.4 keV
+    if (avee > 0.1) &
       call calc_gaussian(avee, tmp_wavflux, 0.4)
 
     ! But only in a few bins!

src/finalize.f90

@@ -13,6 +13,8 @@ subroutine finalize_gitm
   integer :: iError, iBlock, iOutputType
   integer :: nMLTsTmp, nLatsTmp
 
+  call start_timing("Finalize")
+
   if (.not. Is1D) &
     call UA_calc_electrodynamics(nMLTsTmp, nLatsTmp)
 
@@ -31,6 +33,12 @@ subroutine finalize_gitm
     close(iOutputUnit_)
   endif
 
+  if (iProc == 0) then
+    call report_errors
+    call report_warnings
+  endif
+
+  call end_timing("Finalize")
   call end_timing("GITM")
 
   if (iDebugLevel >= 0) call report_timing("all")
@@ -49,9 +57,4 @@ subroutine finalize_gitm
   ! cleanup mpi
   if (.not. IsFrameWork) call MPI_FINALIZE(iError)
 
-  if (iProc == 0) then
-    call report_errors
-    call report_warnings
-  endif
-
 end subroutine finalize_gitm

src/get_potential.f90

@@ -61,19 +61,13 @@ subroutine init_get_potential
 
   ! Now run some checks on user's settings:
   if (iProc == 0) then
-  if (IEModel_%iAurora_ == iOvationPrime_) then
-    if (NormalizeAuroraToHP .and. (iProc == 0)) &
-      call raise_warning("You probably should not use NormalizeAuroraToHP and Ovation")
+    if (IEModel_%iAurora_ == iOvationPrime_) then
+      if (NormalizeAuroraToHP) &
+        call raise_warning("You probably should not use NormalizeAuroraToHP and Ovation")
 
-    if (UseIonAurora .and. IsKappaAurora) &
-      call raise_warning("Kappa aurora & ion precipitation cannot be used simultaneously, yet.")
-  endif
-
-  if (IEModel_%iAurora_ /= iOvationPrime_ .and. UseMonoAurora) &
-    call set_error("Mono-energetic aurora can only be used with ovation currently")
-
-  if (IEModel_%iAurora_ /= iOvationPrime_ .and. UseWaveAurora) &
-    call set_error("Wave aurora can only be used with ovation currently")
+      if (UseIonAurora .and. IsKappaAurora) &
+        call raise_warning("Kappa aurora & ion precipitation cannot be used simultaneously, yet.")
+    endif
   endif
 
   if (IEModel_%iAurora_ /= iFRE_ .and. NormalizeAuroraToHP) &
@@ -127,6 +121,7 @@ subroutine get_potential(iBlock)
   real, dimension(-1:nLons + 2, -1:nLats + 2, 2) :: TempPotential, AMIEPotential
   real, dimension(-1:nLons + 2, -1:nLats + 2) :: Grid, dynamo, SubMLats, SubMLons
   real, dimension(-1:nLons + 2, -1:nLats + 2) :: lats, mlts, EFlux
+  real, dimension(-1:nLons + 2, -1:nLats + 2) :: polarCap
   real :: by, bz, CuspLat, CuspMlt
 
   logical :: UAl_UseGridBasedEIE
@@ -245,6 +240,12 @@ subroutine get_potential(iBlock)
 
   endif
 
+  if (iDebugLevel >= 1) &
+    write(*, *) "==> Min, Max, CPC Potential : ", &
+    minval(Potential(:, :, :, iBlock))/1000.0, &
+    maxval(Potential(:, :, :, iBlock))/1000.0, &
+    (maxval(Potential(:, :, :, iBlock)) - minval(Potential(:, :, :, iBlock)))/1000.0
+
   ! -----------------------------------------------------
   ! Now get the aurora.
   ! This assumes that the field lines are basically
@@ -268,26 +269,29 @@ subroutine get_potential(iBlock)
       call stop_gitm("Stopping in get_potential")
     endif
 
-    ! Sometimes, in AMIE, things get messed up in the
-    ! Average energy, so go through and fix some of these.
-    ! (also checked in aurora.f90)
+    ! Sometimes, in AMIE, things get messed up in the Average energy,
+    ! so go through and fix some of these (Also checked in aurora.f90.)
     if (iDebugLevel > 1) then
-    do iLat = -1, nLats + 2
-      do iLon = -1, nLons + 2
-        if (ElectronAverageEnergyDiffuse(iLon, iLat) < 0.0) then
-          ElectronAverageEnergyDiffuse(iLon, iLat) = 0.1
-          write(*, *) "ave e i,j Negative : ", iLon, iLat, &
-            ElectronAverageEnergyDiffuse(iLon, iLat)
-        endif
-        if (ElectronAverageEnergyDiffuse(iLon, iLat) > 100.0) then
-          write(*, *) "ave e i,j Positive : ", iLon, iLat, &
-            ElectronAverageEnergyDiffuse(iLon, iLat)
-          ElectronAverageEnergyDiffuse(iLon, iLat) = 0.1
-        endif
+      do iLat = -1, nLats + 2
+        do iLon = -1, nLons + 2
+          if (ElectronAverageEnergyDiffuse(iLon, iLat) < 0.0) then
+            ElectronAverageEnergyDiffuse(iLon, iLat) = 0.1
+            write(*, *) "ave e i,j Negative : ", iLon, iLat, &
+              ElectronAverageEnergyDiffuse(iLon, iLat)
+          endif
+          if (ElectronAverageEnergyDiffuse(iLon, iLat) > 100.0) then
+            write(*, *) "ave e i,j Positive : ", iLon, iLat, &
+              ElectronAverageEnergyDiffuse(iLon, iLat)
+            ElectronAverageEnergyDiffuse(iLon, iLat) = 0.1
+          endif
+        enddo
       enddo
-    enddo
     endif
 
+    ! Adjust the Average Energy of the Diffuse Aurora, if desired:
+    if (iDebugLevel > 1) write(*, *) '=> Adjusting average energy of the aurora : ', AveEFactor
+    ElectronAverageEnergyDiffuse = ElectronAverageEnergyDiffuse*AveEFactor
+
     ! -----------------------------
     ! Ion, Wave- & Mono- aurora
     ! -----------------------------
@@ -335,6 +339,12 @@ subroutine get_potential(iBlock)
           endif
         endif
 
+        ! Typical values are:
+        !   - electron average energy = 30-100 eV (< 220 eV)
+        !   - electron energy flux = 1-2 ergs/cm2/s
+        !   - ion average energy = 300 - 3000 keV
+        !   - ion energy flux = ~ 6x less than electrons
+
         EFlux = CuspEFlux* &
                 exp(-abs(lats - CuspLat)/CuspLatHalfWidth)* &
                 exp(-abs(mlts - CuspMlt)/CuspMltHalfWidth)
@@ -361,11 +371,31 @@ subroutine get_potential(iBlock)
 
   endif
 
-  if (iDebugLevel >= 1) &
-    write(*, *) "==> Min, Max, CPC Potential : ", &
-    minval(Potential(:, :, :, iBlock))/1000.0, &
-    maxval(Potential(:, :, :, iBlock))/1000.0, &
-    (maxval(Potential(:, :, :, iBlock)) - minval(Potential(:, :, :, iBlock)))/1000.0
+  ! -----------------------------
+  ! Polar Rain, if desired
+  ! -----------------------------
+
+  if (UsePolarRain) then
+
+    ! Get the polar cap - this variable is 1 if it is the polar cap and 0 is not
+    call IEModel_%get_polarcap(polarCap)
+
+    ! Then, if we are in the polar cap, fill in the energy flux with values
+    ! entered by the user in the UAM.in file
+    ! Typical values of polar rain are [Newell et al., 1990]:
+    ! aveE ~ 80 - 110 eV
+    ! eflux ~ 0.01 - 0.1 ergs/cm2/s for nominal events; 1.0 for extreme events
+    do iLat = -1, nLats + 2
+      do iLon = -1, nLons + 2
+        if (polarCap(iLon, iLat) > 0 .and. &
+            ElectronEnergyFluxDiffuse(iLon, iLat) < 0.1) then
+          ElectronEnergyFluxDiffuse(iLon, iLat) = polarRainEFlux
+          ElectronAverageEnergyDiffuse(iLon, iLat) = polarRainAveE
+        endif
+      enddo
+    enddo
+
+  endif
 
   call end_timing("get_potential")
 

src/logfile.f90

@@ -274,6 +274,7 @@ subroutine write_code_information(dir)
   use ModIoUnit, ONLY: io_unit_new
   use ModUtilities, ONLY: flush_unit
   use ModRCMR
+  use ModGITMVersion
 
   implicit none
 
@@ -291,6 +292,10 @@ subroutine write_code_information(dir)
 
     write(iCodeInfoFileUnit_, *) "GITM2 Run Information"
     write(iCodeInfoFileUnit_, *) "---------------------"
+    write(iCodeInfoFileUnit_, *) "GITM Version ", trim(GitmVersion)
+    write(iCodeInfoFileUnit_, *) "GITM Full Version ", trim(GitmVersionFull)
+    write(iCodeInfoFileUnit_, *) "Electrodynamics Version ", trim(ElectrodynamicsVersionFull)
+    write(iCodeInfoFileUnit_, *) "---------------------"
     write(iCodeInfoFileUnit_, *) ""
 
     write(iCodeInfoFileUnit_, *) "nSpecies", nSpecies

src/output_common.f90

@@ -644,15 +644,15 @@ subroutine output_header
 
     if (cType(3:5) == "THM") then
 
-      write(iOutputUnit_, "(I7,A1,a)") 4, " ", "EUV Heating"
-      write(iOutputUnit_, "(I7,A1,a)") 5, " ", "Conduction"
-      write(iOutputUnit_, "(I7,A1,a)") 6, " ", "Molecular Conduction"
-      write(iOutputUnit_, "(I7,A1,a)") 7, " ", "Eddy Conduction"
-      write(iOutputUnit_, "(I7,A1,a)") 8, " ", "Eddy Adiabatic Conduction"
-      write(iOutputUnit_, "(I7,A1,a)") 9, " ", "Chemical Heating"
-      write(iOutputUnit_, "(I7,A1,a)") 11, " ", "Joule Heating"
-      write(iOutputUnit_, "(I7,A1,a)") 12, " ", "NO Cooling"
-      write(iOutputUnit_, "(I7,A1,a)") 13, " ", "O Cooling"
+      write(iOutputUnit_, "(I7,A1,a)") 4, " ", "EUV Heating (K/s)"
+      write(iOutputUnit_, "(I7,A1,a)") 5, " ", "Conduction (K/s)"
+      write(iOutputUnit_, "(I7,A1,a)") 6, " ", "Molecular Conduction (K/s)"
+      write(iOutputUnit_, "(I7,A1,a)") 7, " ", "Eddy Conduction (K/s)"
+      write(iOutputUnit_, "(I7,A1,a)") 8, " ", "Eddy Adiabatic Conduction (K/s)"
+      write(iOutputUnit_, "(I7,A1,a)") 9, " ", "Chemical Heating (K/s)"
+      write(iOutputUnit_, "(I7,A1,a)") 11, " ", "Joule Heating (K/s)"
+      write(iOutputUnit_, "(I7,A1,a)") 12, " ", "NO Cooling (K/s)"
+      write(iOutputUnit_, "(I7,A1,a)") 13, " ", "O Cooling (K/s)"
       write(iOutputUnit_, "(I7,A1,a)") 14, " ", "Total Abs EUV"
       if (cType(1:2) == "1D") then
         do iSpecies = 1, nSpeciesTotal
@@ -1346,16 +1346,16 @@ subroutine output_3dthm(iBlock)
           Longitude(iLon, iBlock), &
           Latitude(iLat, iBlock), &
           Altitude_GB(iLon, iLat, iAlt, iBlock), &
-          EuvHeating(iiLon, iiLat, iiAlt, iBlock)*dt*TempUnit(iiLon, iiLat, iiAlt), &
-          Conduction(iiLon, iiLat, iiAlt)*TempUnit(iiLon, iiLat, iiAlt), &
+          EuvHeating(iiLon, iiLat, iiAlt, iBlock)*TempUnit(iiLon, iiLat, iiAlt), &
+          Conduction(iiLon, iiLat, iiAlt)*TempUnit(iiLon, iiLat, iiAlt)/dt, &
           MoleConduction(iiLon, iiLat, iiAlt), &
           EddyCond(iiLon, iiLat, iiAlt), &
           EddyCondAdia(iiLon, iiLat, iiAlt), &
-          ChemicalHeatingRate(iiLon, iiLat, iiAlt)*TempUnit(iiLon, iiLat, iiAlt), &
-          JouleHeating(iiLon, iiLat, iiAlt)*dt*TempUnit(iiLon, iiLat, iiAlt), &
-          -NOCooling(iiLon, iiLat, iiAlt)*dt*TempUnit(iiLon, iiLat, iiAlt), &
-          -OCooling(iiLon, iiLat, iiAlt)*dt*TempUnit(iiLon, iiLat, iiAlt), &
-          EuvTotal(iiLon, iiLat, iiAlt, iBlock)*dt, &
+          ChemicalHeatingRate(iiLon, iiLat, iiAlt)*TempUnit(iiLon, iiLat, iiAlt)/dt, &
+          JouleHeating(iiLon, iiLat, iiAlt)*TempUnit(iiLon, iiLat, iiAlt), &
+          -NOCooling(iiLon, iiLat, iiAlt)*TempUnit(iiLon, iiLat, iiAlt), &
+          -OCooling(iiLon, iiLat, iiAlt)*TempUnit(iiLon, iiLat, iiAlt), &
+          EuvTotal(iiLon, iiLat, iiAlt, iBlock), &
           cp(iiLon, iiLat, iiAlt, iBlock), &
           rho(iiLon, iiLat, iiAlt, iBlock), &
           sqrt(sum(EField(iLon, iLat, iAlt, :)**2)), & ! magnitude of E.F.

src/set_inputs.f90

@@ -796,6 +796,21 @@ subroutine set_inputs
           IsDone = .true.
         endif
 
+      case ("#USEPOLARRAIN")
+        call read_in_logical(UsePolarRain, iError)
+        call read_in_real(polarRainAveE, iError)
+        call read_in_real(polarRainEFlux, iError)
+        if (iError /= 0) then
+          write(*, *) 'Incorrect format for #USECUSP'
+          write(*, *) 'This is for specifying a cusp.'
+          write(*, *) ''
+          write(*, *) '#USEPOLARRAIN'
+          write(*, *) 'UsePolarRain        (logical)'
+          write(*, *) 'PolarRainAveE       (real)'
+          write(*, *) 'PolarRainEFlux      (real)'
+          IsDone = .true.
+        endif
+
       case ("#USEREGIONALAMIE")
         call read_in_logical(UseRegionalAMIE, iError)
         call read_in_logical(UseTwoAMIEPotentials, iError)