Pre-release Bugfixes & Documentation Update

README.md

@@ -7,27 +7,7 @@ on Linux and MacOS as well as ifort on NASA's Pleiades computer.
 
 For the complete documentation, see [GITM's Read the Docs Page](https://gitm.readthedocs.io).
 
----
-
-<details>
-<summary><b>GITM's default branch recently changed names. To access the latest features, please update your local refs.</b></summary>
-
-From your local clone of this repository, run the following commands to update the name of the default branch.
-
-```sh
-git branch -m master main
-git fetch origin
-git branch -u origin/main main
-git remote set-head origin -a
-```
-
-Optionally, run the following command to remove tracking references to the old branch name.
-
-```sh
-git remote prune origin
-```
-
-</details>
+GITM's stable version is the main branch, which is downloaded by default.  If you want the latest changes, but are ok with possibly unstable code, you can use the 'develop' version, which is described below.  If you don't know what you are doing, please just use the main version (i.e., don't checkout a different branch).
 
 ## Quick Start
 
@@ -48,7 +28,13 @@ Substitute the URL with the `https` link from the "Code" button above if you do
 cd GITM
 ```
 
-3\. Configure the Fortran compiler and download external electrodynamics library
+(if you want/need to change to a different branch, do that here with the command:
+```shell
+git checkout develop
+```
+but, again, we don't recommend this unless you know what you are doing!)
+
+3\. Configure the Fortran compiler and download the external electrodynamics library (the install should do this automagically):
 
 ```shell
 ./Config.pl -install -earth -compiler=gfortran
@@ -58,16 +44,16 @@ The above command is that it assumes that you have a working gfortran compiler a
 things like mpif90 work ok.  If you don't have gfortran and mpif90, then you need
 to get these things for your computer. 
 
-> `Config.pl` should automatically your gfortran version since gfortran>=10 needs
+> `Config.pl` should automatically determine your gfortran version since gfortran>=10 needs
 > to use different compilation commands than versions 9 and below.
-> If you notice that `Config.pl` does not detect the correct gfoetran version, set
+> If you notice that `Config.pl` does not detect the correct gfortran version, set
 > `-compiler=gfortran` if you have gfortran <10 or `-compiler=gfortran10` for 
 > gfortran>=10
 
-
 In theory, Mars, Venus, Titan, and LV-426 should work.  These are in
 various states of completion, so I wouldn't count on them being
-perfect.
+perfect. To configure with one of these, simple use the planet as
+an option (like '-venus', instead of '-earth)'.
 
 If running on Pleiades, you need to have these
 in your start-up script (.cshrc, .bashrc, etc):
@@ -130,13 +116,13 @@ lon. See below for how to set the resolution.
 9\. Go into the output directory:
 
 ```shell
-cd data
+cd UA/data
 ```
 
 10\. Make some plots with an old plotter:
 
 ```shell
-../../../srcPython/plot_model_results.py -var=3 -alt=120 3DALL_t021221_000500.bin
+../../../srcPython/plot_model_results.py -var=3 -alt=300 3DALL_t021221_000500.bin
 ```
 
 Then look at the png file that is created.  You can use a `-h` to see

src/ModGITM.f90

@@ -135,6 +135,8 @@ module ModGITM
     VerticalIonDrag
 
   real, allocatable :: Potential(:, :, :, :), PotentialY(:, :, :, :)
+  ! Separate potentials, leaving nBlocks out:
+  real, allocatable :: MagnetosphericPotential(:, :, :), DynamoPotential(:, :, :)
 
   real, dimension(-1:nLons + 2, -1:nLats + 2, -1:nAlts + 2, 3) :: &
     ExB, EField
@@ -171,8 +173,7 @@ module ModGITM
 
   real :: SubsolarLatitude, SubsolarLongitude
   real :: MagneticPoleColat, MagneticPoleLon
-  real :: PotentialMax_North, PotentialMin_North
-  real :: PotentialMax_South, PotentialMin_South
+  real :: CPCPn, CPCPs
   real :: HemisphericPowerNorth, HemisphericPowerSouth
   real :: HemisphericPowerNorth_diffuse, HemisphericPowerSouth_diffuse
   real :: HemisphericPowerNorth_wave, HemisphericPowerSouth_wave
@@ -254,6 +255,8 @@ subroutine init_mod_gitm
     allocate(cMax_GDB(0:nLons + 1, 0:nLats + 1, 0:nAlts + 1, 3, nBlocks))
     allocate(IonPressureGradient(-1:nLons + 2, -1:nLats + 2, -1:nAlts + 2, 3, nBlocks))
     allocate(Potential(-1:nLons + 2, -1:nLats + 2, -1:nAlts + 2, nBlocks))
+    allocate(MagnetosphericPotential(-1:nLons + 2, -1:nLats + 2, -1:nAlts + 2))
+    allocate(DynamoPotential(-1:nLons + 2, -1:nLats + 2, -1:nAlts + 2))
     allocate(PotentialY(-1:nLons + 2, -1:nLats + 2, -1:nAlts + 2, nBlocks))
     PotentialY = 0.0
     allocate(Velocity(-1:nLons + 2, -1:nLats + 2, -1:nAlts + 2, 3, nBlocks))

src/ModInputs.f90

@@ -75,6 +75,9 @@ module ModInputs
   logical :: UseWaveAurora = .false.
   logical :: UseIonAurora = .false.
 
+  logical :: doStopIfNoAurora = .false.
+  logical :: doStopIfNoPotential = .false.
+
   ! logical :: UseOvationSMEMono = .false.
   ! logical :: UseOvationSMEWave = .false.
   ! logical :: UseOvationSMEIon = .false.

src/get_potential.f90

@@ -130,8 +130,11 @@ subroutine get_potential(iBlock)
   logical :: IsFirstTime = .true.
   logical :: IsFirstPotential(nBlocksMax) = .true.
   logical :: IsFirstAurora(nBlocksMax) = .true.
-  real    :: mP, dis, TempWeight
-  real    ::  LocalSumDiffPot, MeanDiffPot, LatBoundNow
+  real :: mP, dis, TempWeight
+  real :: LocalSumDiffPot, MeanDiffPot, LatBoundNow
+  real :: PotentialMin_North, PotentialMin_South
+  real :: PotentialMax_North, PotentialMax_South
+  real :: maxEnergyFlux
 
   real, dimension(-1:nLons + 2, -1:nLats + 2) :: TempPotential2d
   real, dimension(-1:nLons + 2, -1:nLats + 2, 2) :: TempPotential, AMIEPotential
@@ -150,7 +153,6 @@ subroutine get_potential(iBlock)
   if (index(cPlanet, "Earth") == 0) then
 
     potential = 0.0
-    ElectronEnergyFluxDiffuse = 0.1
     ElectronEnergyFluxDiffuse = 0.0001
     return
 
@@ -169,7 +171,7 @@ subroutine get_potential(iBlock)
     call set_ie_indices(IEModel_, CurrentTime)
 
     Potential(:, :, :, iBlock) = 0.0
-
+    DynamoPotential = 0.0
     do iAlt = -1, nAlts + 2
 
       call iemodel_%grid( &
@@ -190,6 +192,7 @@ subroutine get_potential(iBlock)
 
       call iemodel_%get_potential(TempPotential2d)
       TempPotential(:, :, 1) = TempPotential2d
+      MagnetosphericPotential(:, :, iAlt) = TempPotential2d
 
       if (.not. isOk) then
         call set_error("Error in routine get_potential (getting potential):")
@@ -204,6 +207,7 @@ subroutine get_potential(iBlock)
         call get_dynamo_potential( &
           MLongitude(-1:nLons + 2, -1:nLats + 2, iAlt, iBlock), &
           MLatitude(-1:nLons + 2, -1:nLats + 2, iAlt, iBlock), dynamo)
+        DynamoPotential(:, :, iAlt) = dynamo
 
         ! Set latitude boundary between region of high lat convection
         ! and region of neutral wind dyanmo based on if SWMF potential
@@ -256,26 +260,39 @@ subroutine get_potential(iBlock)
 
   endif
 
+  PotentialMax_North = 0.0
+  PotentialMin_North = 0.0
+  PotentialMax_South = 0.0
+  PotentialMin_South = 0.0
+
   if (maxval(Mlatitude(:, :, :, iBlock)) > 45.0) then
-    PotentialMax_North = maxval(Potential(:, :, :, iBlock))/1000.0
-    PotentialMin_North = minval(Potential(:, :, :, iBlock))/1000.0
-  else
-    PotentialMax_North = 0.0
-    PotentialMin_North = 0.0
+    PotentialMax_North = maxval(MagnetosphericPotential(:, :, :))/1000.0
+    PotentialMin_North = minval(MagnetosphericPotential(:, :, :))/1000.0
   endif
-  if (maxval(Mlatitude(:, :, :, iBlock)) < -45.0) then
-    PotentialMax_South = maxval(Potential(:, :, :, iBlock))/1000.0
-    PotentialMin_South = minval(Potential(:, :, :, iBlock))/1000.0
-  else
-    PotentialMax_South = 0.0
-    PotentialMin_South = 0.0
+  if (minval(Mlatitude(:, :, :, iBlock)) < -45.0) then
+    PotentialMax_South = maxval(MagnetosphericPotential(:, :, :))/1000.0
+    PotentialMin_South = minval(MagnetosphericPotential(:, :, :))/1000.0
   endif
 
+  call get_min_value_across_pes(PotentialMin_North)
+  call get_max_value_across_pes(PotentialMax_North)
+  call get_min_value_across_pes(PotentialMin_South)
+  call get_max_value_across_pes(PotentialMax_South)
+
+  ! Volt -> kV
+  CPCPn = (PotentialMax_North - PotentialMin_North)
+  CPCPs = (PotentialMax_South - PotentialMin_South)
+
   if (iDebugLevel >= 1) &
-    write(*, *) "==> Min, Max, CPC Potential : ", &
-    int(minval(Potential(:, :, :, iBlock))/1000.0), &
-    int(maxval(Potential(:, :, :, iBlock))/1000.0), &
-    int((maxval(Potential(:, :, :, iBlock)) - minval(Potential(:, :, :, iBlock)))/1000.0)
+    write(*, *) "==> CPC Potential (North, South): ", CPCPn, CPCPs
+
+  if ((CPCPn == 0) .and. (CPCPs == 0)) then
+    if (doStopIfNoPotential) then
+      call set_error("CPCP is zero and model is not zero! Must Stop!")
+      call report_errors
+      call stop_gitm("Stopping in get_potential")
+    endif
+  endif
 
   ! -----------------------------------------------------
   ! Now get the aurora.
@@ -332,6 +349,15 @@ subroutine get_potential(iBlock)
       enddo
     endif
 
+    ! Check to see if there is any Diffuse aurora and stop if there is none:
+    maxEnergyFlux = maxval(ElectronEnergyFluxDiffuse)
+    call get_max_value_across_pes(maxEnergyFlux)
+    if ((maxEnergyFlux == 0) .and. (doStopIfNoAurora)) then
+      call set_error("Wave Aurora is zero and it was requested! Must Stop!")
+      call report_errors
+      call stop_gitm("Stopping in get_potential")
+    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
@@ -340,14 +366,41 @@ subroutine get_potential(iBlock)
     ! Ion, Wave- & Mono- aurora
     ! -----------------------------
 
-    if (UseWaveAurora) &
+    if (UseWaveAurora) then
       call IEModel_%get_electron_wave_aurora(ElectronEnergyFluxWave, ElectronAverageEnergyWave)
+      ! Check to see if there is any aurora:
+      maxEnergyFlux = maxval(ElectronEnergyFluxWave)
+      call get_max_value_across_pes(maxEnergyFlux)
+      if (maxEnergyFlux == 0) then
+        call set_error("Wave Aurora is zero and it was requested! Must Stop!")
+        call report_errors
+        call stop_gitm("Stopping in get_potential")
+      endif
+    endif
 
-    if (UseMonoAurora) &
+    if (UseMonoAurora) then
       call IEModel_%get_electron_mono_aurora(ElectronEnergyFluxMono, ElectronAverageEnergyMono)
+      ! Check to see if there is any aurora:
+      maxEnergyFlux = maxval(ElectronEnergyFluxMono)
+      call get_max_value_across_pes(maxEnergyFlux)
+      if (maxEnergyFlux == 0) then
+        call set_error("Mono Aurora is zero and it was requested! Must Stop!")
+        call report_errors
+        call stop_gitm("Stopping in get_potential")
+      endif
+    endif
 
-    if (UseIonAurora) &
+    if (UseIonAurora) then
       call IEModel_%get_ion_diffuse_aurora(IonEnergyFlux, IonAverageEnergy)
+      ! Check to see if there is any aurora:
+      maxEnergyFlux = maxval(IonEnergyFlux)
+      call get_max_value_across_pes(maxEnergyFlux)
+      if (maxEnergyFlux == 0) then
+        call set_error("Ion Aurora is zero and it was requested! Must Stop!")
+        call report_errors
+        call stop_gitm("Stopping in get_potential")
+      endif
+    endif
 
     ! -----------------------------
     ! Cusp, if desired

src/library.f90

@@ -5,6 +5,48 @@
 !
 !---------------------------------------------------------------------------
 
+subroutine get_min_value_across_pes(value)
+
+  use ModGITM, only: iCommGITM
+  use ModMpi
+
+  implicit none
+  real, intent(inout) :: value
+  real :: localVar
+  integer :: iError
+
+  LocalVar = value
+  call MPI_REDUCE(LocalVar, value, 1, MPI_REAL, MPI_MIN, &
+                  0, iCommGITM, iError)
+  call MPI_BCAST(value, 1, MPI_Real, 0, iCommGITM, iError)
+
+end subroutine get_min_value_across_pes
+
+!---------------------------------------------------------------------------
+!
+!---------------------------------------------------------------------------
+
+subroutine get_max_value_across_pes(value)
+
+  use ModGITM, only: iCommGITM
+  use ModMpi
+
+  implicit none
+  real, intent(inout) :: value
+  real :: localVar
+  integer :: iError
+
+  LocalVar = value
+  call MPI_REDUCE(LocalVar, value, 1, MPI_REAL, MPI_MAX, &
+                  0, iCommGITM, iError)
+  call MPI_BCAST(value, 1, MPI_Real, 0, iCommGITM, iError)
+
+end subroutine get_max_value_across_pes
+
+!---------------------------------------------------------------------------
+!
+!---------------------------------------------------------------------------
+
 subroutine report(str, iLevel)
 
   use ModInputs, only: iDebugLevel

src/logfile.f90

@@ -112,7 +112,6 @@ subroutine logfile(dir)
   real    :: minTemp, maxTemp, localVar, minVertVel, maxVertVel
   real    :: AverageTemp, AverageVertVel, TotalVolume, Bx, By, Bz, Vx, Hpi
   real    :: HPn, HPs, HPn_d, HPs_d, HPn_w, HPs_w, HPn_m, HPs_m
-  real    :: PotMaxN, PotMaxS, PotMinN, PotMinS, CPCPn, CPCPs
   real    :: SSLon, SSLat, SSVTEC
   integer :: iError
 
@@ -202,25 +201,6 @@ subroutine logfile(dir)
   call MPI_REDUCE(LocalVar, AverageVertVel, 1, MPI_REAL, MPI_SUM, &
                   0, iCommGITM, iError)
 
-  LocalVar = PotentialMin_North
-  call MPI_REDUCE(LocalVar, PotMinN, 1, MPI_REAL, MPI_MIN, &
-                  0, iCommGITM, iError)
-
-  LocalVar = PotentialMin_South
-  call MPI_REDUCE(LocalVar, PotMinS, 1, MPI_REAL, MPI_MIN, &
-                  0, iCommGITM, iError)
-
-  LocalVar = PotentialMax_North
-  call MPI_REDUCE(LocalVar, PotMaxN, 1, MPI_REAL, MPI_MAX, &
-                  0, iCommGITM, iError)
-
-  LocalVar = PotentialMax_South
-  call MPI_REDUCE(LocalVar, PotMaxS, 1, MPI_REAL, MPI_MAX, &
-                  0, iCommGITM, iError)
-
-  CPCPn = PotMaxN - PotMinN
-  CPCPs = PotMaxS - PotMinS
-
   LocalVar = HemisphericPowerNorth
   call MPI_REDUCE(LocalVar, HPn, 1, MPI_REAL, MPI_SUM, &
                   0, iCommGITM, iError)

src/set_inputs.f90

@@ -709,13 +709,23 @@ subroutine set_inputs
           write(*, *) '(and low-latitude) electrodynamic drivers, such as'
           write(*, *) 'the potential and the aurora.'
           write(*, *) '#ELECTRODYNAMICS'
-          write(*, *) 'AuroralModel     [fta], fre, pem, ovation, hpi/ihp, amie, etc.'
+          write(*, *) 'AuroralModel     [fta], fre, pem, ovation, hpi/ihp, amie, zero, etc.'
           write(*, *) 'DtAurora    (real, seconds [60.0])'
-          write(*, *) 'PotentialModel   [weimer05], hepmay, amie, etc.'
+          write(*, *) 'PotentialModel   [weimer05], hepmay, amie, zero, etc.'
           write(*, *) 'DtPotential (real, seconds [60.0])'
           IsDone = .true.
         endif
 
+        call lower_case(cAuroralModel)
+        call lower_case(cPotentialModel)
+
+        ! This is saying if auroral model is NOT zero, then stop if there is no aurora
+        if (index(cAuroralModel, 'zero') == 0) &
+          doStopIfNoAurora = .true.
+        ! This is saying if potential model is NOT zero, then stop if there is no potential
+        if (index(cPotentialModel, 'zero') == 0) &
+          doStopIfNoPotential = .true.
+
       case ("#AMIEFILES")
         call read_in_string(cAMIEFileNorth, iError)
         call read_in_string(cAMIEFileSouth, iError)