Amie file mods, some bugs!

Config.pl

@@ -109,7 +109,7 @@
 
 &set_planet if $Planet and $Planet ne $PlanetOrig;
 
-&show_settings if $Show; 
+&show_settings_ if $Show; 
 
 print "Config.pl -g=$GridSize\n" if $ShowGridSize and not $Show;
 

README.md

@@ -49,21 +49,21 @@ cd GITM
 ```
 
 3\. Configure the Fortran compiler and download external electrodynamics library
-(gfortran versions 10+)
 
 ```shell
-./Config.pl -install -earth -compiler=gfortran10
+./Config.pl -install -earth -compiler=gfortran
 ```
 
-The biggest issue with the above command is that it assumes that you
-have the gfortran (version 10+) compiler and things like mpif90 work
-ok.  If you don't have gfortran and mpif90, then you need to get these
-things for your computer.  If you have version 9 or before for gfortran,
-you can do:
+The above command is that it assumes that you have a working gfortran compiler and
+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
+> to use different compilation commands than versions 9 and below.
+> If you notice that `Config.pl` does not detect the correct gfoetran version, set
+> `-compiler=gfortran` if you have gfortran <10 or `-compiler=gfortran10` for 
+> gfortran>=10
 
-```shell
-./Config.pl -install -earth -compiler=gfortran
-```
 
 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

share/Scripts/Config.pl

@@ -433,6 +433,18 @@ sub install_code_{
 	&shell_command("echo COMPILER=$Compiler >> $MakefileDef");
 	&shell_command("cat $MakefileDefOrig >> $MakefileDef");
 
+	# Check gfortran version and automatically select gfortran10 if >= 10
+	if ($Compiler eq 'gfortran') {
+	    my $version_output = `gfortran --version 2>/dev/null | head -n1`;
+	    if ($version_output =~ /(\d+)\.\d+\.\d+/) {
+		my $major_version = $1;
+		if ($major_version >= 10) {
+		    $Compiler = 'gfortran10';
+		}
+		print "Detected gfortran version $major_version, using $Compiler makefile\n";
+	    }
+	}
+
 	my $Makefile = "$MakefileConfOrig.$OS.$Compiler";
 	if(-f $Makefile){
 	    &shell_command("cat $Makefile > $MakefileConf");

src/ModGITM.f90

@@ -171,6 +171,8 @@ module ModGITM
 
   real :: SubsolarLatitude, SubsolarLongitude
   real :: MagneticPoleColat, MagneticPoleLon
+  real :: PotentialMax_North, PotentialMin_North
+  real :: PotentialMax_South, PotentialMin_South
   real :: HemisphericPowerNorth, HemisphericPowerSouth
   real :: HemisphericPowerNorth_diffuse, HemisphericPowerSouth_diffuse
   real :: HemisphericPowerNorth_wave, HemisphericPowerSouth_wave

src/ModInputs.f90

@@ -323,6 +323,9 @@ module ModInputs
   ! Lower limit on ion density
   real :: MinIonDensity = 100.0
   real :: MinIonDensityAdvect = 1e5
+  !Lower limits on neutral density
+  real :: MinNeutralDensity = 200.0
+  real :: MinNeutralDensityAdvect = 1e5
 
   !\
   ! Methods for completing chemistry

src/advance_horizontal.f90

@@ -361,8 +361,13 @@ subroutine advance_horizontal(iBlock)
           do iSpecies = 1, nSpecies
             if (NewNum_CV(iLon, iLat, iSpecies) < 0.0) then
               write(*, *) "Species : ", iLon, iLat, iAlt, iSpecies, iBlock
-              call stop_gitm("neg ion density! stopping in advance_horizontal")
-              NewNum_CV(iLon, iLat, iSpecies) = 1.0
+              ! call stop_gitm("neg ion density! stopping in advance_horizontal")
+              if (iSpecies <= nSpecies) then
+                NewNum_CV(iLon, iLat, iSpecies) = MinNeutralDensityAdvect
+              else
+                NewNum_CV(iLon, iLat, iSpecies) = MinNeutralDensity
+              endif
+
               IsFound = .true.
             endif
           enddo

src/aurora.Earth.f90

@@ -379,20 +379,31 @@ subroutine NormalizeDiffuseAuroraToHP
     call calculate_HP(HPn, HPs)
 
     ! Collect all of the powers by summing them together
-    LocalVar = HemisphericPowerNorth/1.0e9
+    LocalVar = HPn/1.0e9
     call MPI_REDUCE(LocalVar, HPn, 1, MPI_REAL, MPI_SUM, &
                     0, iCommGITM, iError)
 
-    LocalVar = HemisphericPowerSouth/1.0e9
+    LocalVar = HPs/1.0e9
     call MPI_REDUCE(LocalVar, HPs, 1, MPI_REAL, MPI_SUM, &
                     0, iCommGITM, iError)
 
     ! Average north and south together
     avepower = (HPn + HPs)/2.0
 
+    if (avepower == 0) then
+      call set_error("AvePower is zero!")
+      avepower = 1.0
+    endif
+
     ! If we are only have one hemisphere or the other, assign to avepower
-    if (HPs < 0.1*HPn) avepower = HPn
-    if (HPn < 0.1*HPs) avepower = HPs
+    if (HPs < 0.1*HPn) then
+      avepower = HPn
+      HPs = HPn
+    endif
+    if (HPn < 0.1*HPs) then
+      avepower = HPs
+      HPn = HPs
+    endif
 
     call MPI_Bcast(avepower, 1, MPI_Real, 0, iCommGITM, ierror)
     call MPI_Bcast(Hps, 1, MPI_Real, 0, iCommGITM, ierror)
@@ -411,41 +422,40 @@ subroutine NormalizeDiffuseAuroraToHP
     else
 
       call get_hpi(CurrentTime, Hpi, iError)
-      if (avepower == 0) then
-        call set_error("AvePower is zero!")
-        avepower = 1.0
-      endif
-      ratio_NH = Hpi/avepower
-      ratio_SH = Hpi/avepower
+      ratio_NH = Hpi/HPn
+      ratio_SH = Hpi/HPs
 
     endif
 
     if (iDebugLevel >= 0) then
       if ((iDebugLevel == 0) .and. IsFirstTime(iBlock) .and. .not. doSeparateHPI) then
         write(*, *) '---------------------------------------------------'
         write(*, *) 'Using auroral normalizing ratios!!! '
-        write(*, *) 'no longer reporting!'
+        write(*, *) '  -> forcing the HP to be the same in both hemispheres!'
+        write(*, *) '  -> no longer reporting!'
         write(*, *) '---------------------------------------------------'
       elseif ((iDebugLevel == 0) .and. IsFirstTime(iBlock) .and. doSeparateHPI) then
         write(*, *) '---------------------------------------------------'
         write(*, *) 'Using HPI from each hemisphere to normalize aurora!!'
-        write(*, *) 'no longer reporting!'
+        write(*, *) '  -> Normalizing to mean, then calculating ratios seperately'
+        write(*, *) '  -> no longer reporting!'
         write(*, *) '---------------------------------------------------'
       endif
-      if (iDebugLevel >= 2) then
+      if (iDebugLevel >= 1) then
         if (doSeparateHPI) then
           write(*, *) 'Auroral normalizing ratios: '
-          write(*, *) 'Hpi(NH)  HPI(SH)  HPI(NH-modeled)  HPI(SH-modeled)  ratio_n  ratio_s'
+          write(*, *) 'HP(N)  HP(S)  HP(N-model)  HP(S-model)  ratio_n  ratio_s'
           write(*, *) Hpi_NH, Hpi_SH, HPn, HPs, ratio_NH, ratio_sh
         else
-          write(*, *) 'Auroral normalizing ratio: ', Hpi_NH, ratio_SH, ratio_NH
+          write(*, *) 'Auroral normalizing ratio (hp, hp model averaged, ratios (S/N)): ', &
+            Hpi, avepower, ratio_SH, ratio_NH
         endif
       endif
 
     endif
     do i = 1, nLats
       do j = 1, nLons
-        if (ElectronEnergyFluxDiffuse(j, i) > 0.1) then
+        if (ElectronEnergyFluxDiffuse(j, i) > 0.001) then
           if (latitude(i, iBlock) < 0.0) then
             ElectronEnergyFluxDiffuse(j, i) = ElectronEnergyFluxDiffuse(j, i)*ratio_sh
           else
@@ -463,12 +473,15 @@ subroutine calculate_HP(HPn, HPs)
     real, intent(out) :: HPn, HPs
     real :: eflx_ergs, eflux
 
+    HPn = 0.0
+    HPs = 0.0
+
     do i = 1, nLats
       do j = 1, nLons
 
         eflx_ergs = ElectronEnergyFluxDiffuse(j, i) !/ (1.0e-7 * 100.0 * 100.0)
 
-        if (eflx_ergs > 0.1) then
+        if (eflx_ergs > 0.001) then
           eflux = eflx_ergs*6.242e11  ! ergs/cm2/s -> eV/cm2/s
 
           !(eV/cm2/s -> J/m2/s)
@@ -477,9 +490,9 @@ subroutine calculate_HP(HPn, HPs)
                   dLonDist_FB(j, i, nAlts, iBlock)
 
           if (latitude(i, iBlock) < 0.0) then
-            HemisphericPowerSouth = HemisphericPowerSouth + power
+            HPs = HPs + power
           else
-            HemisphericPowerNorth = HemisphericPowerNorth + power
+            HPn = HPn + power
           endif
 
         endif

src/finalize.f90

@@ -15,8 +15,8 @@ subroutine finalize_gitm
 
   call start_timing("Finalize")
 
-  if (.not. Is1D) &
-    call UA_calc_electrodynamics(nMLTsTmp, nLatsTmp)
+  ! if (.not. Is1D) &
+  !   call UA_calc_electrodynamics(nMLTsTmp, nLatsTmp)
 
   do iOutputType = 1, nOutputTypes
     do iBlock = 1, nBlocks

src/get_potential.f90

@@ -256,6 +256,21 @@ subroutine get_potential(iBlock)
 
   endif
 
+  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
+  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
+  endif
+
   if (iDebugLevel >= 1) &
     write(*, *) "==> Min, Max, CPC Potential : ", &
     int(minval(Potential(:, :, :, iBlock))/1000.0), &

src/logfile.f90

@@ -112,6 +112,7 @@ 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
 
@@ -162,9 +163,10 @@ subroutine logfile(dir)
     write(iLogFileUnit_, '(a)') " "
     write(iLogFileUnit_, '(a)') "#START"
     write(iLogFileUnit_, '(a)') &
-      "   iStep yyyy mm dd hh mm ss  ms      dt "// &
+      "   iStep year month day hour min sec  ms      dt "// &
       "min(T) max(T) mean(T) min(VV) max(VV) mean(VV) F107 F107A "// &
       "By Bz Vx HP HPn HPs HPn_diff HPs_diff HPn_w HPs_w HPn_m HPs_m "// &
+      "CPCPn CPCPs "// &
       "SubsolarLon SubsolarLat SubsolarVTEC"
   endif
 
@@ -200,6 +202,25 @@ 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)
@@ -250,13 +271,15 @@ subroutine logfile(dir)
 
     if (Is1D) SSVTEC = -1.0
 
-    write(iLogFileUnit_, "(i8,i5,5i3,i4,f8.4,6f13.5,8f9.1,10f10.5,10f10.5,10f8.3)") &
-      iStep, iTimeArray, dt, minTemp, maxTemp, AverageTemp, &
-      minVertVel, maxVertVel, AverageVertVel, &
-      f107, f107A, By, Bz, Vx, Hpi, &
-      HPn/1.0e9, HPs/1.0e9, HPn_d/1.0e9, HPs_d/1.0e9, &
-      HPn_w/1.0e9, HPs_w/1.0e9, HPn_m/1.0e9, HPs_m/1.0e9, &
-      SSLon, SSLat, SSVTEC
+    write(iLogFileUnit_, "(i8,i5,5i3,i4,f8.4,6f9.1,5f7.1,9f8.1,2f7.1,3f8.3)") &
+      iStep, iTimeArray, &  ! i8, i5, 5i3, i4
+      dt, &  ! f8.4
+      minTemp, maxTemp, AverageTemp, minVertVel, maxVertVel, AverageVertVel, & ! 6f13.5
+      f107, f107A, By, Bz, Vx, &  ! 5f9.1
+      Hpi, HPn/1.0e9, HPs/1.0e9, &
+      HPn_d/1.0e9, HPs_d/1.0e9, HPn_w/1.0e9, HPs_w/1.0e9, HPn_m/1.0e9, HPs_m/1.0e9, & ! 9f10.1
+      CPCPn, CPCPs, & ! 2f10.2
+      SSLon, SSLat, SSVTEC   ! 3f8.3
 
     call flush_unit(iLogFileUnit_)
   endif

src/set_inputs.f90

@@ -733,7 +733,7 @@ subroutine set_inputs
         endif
 
       case ("#AMIENORTH")
-        nAMIENorth = 0 
+        nAMIENorth = 0
         call read_in_int(nAMIENorth, iError)
         if (iError == 0 .and. nAMIENorth > 0) then
           allocate(cAMIEListNorth(nAMIENorth))
@@ -755,7 +755,7 @@ subroutine set_inputs
         endif
 
       case ("#AMIESOUTH")
-        nAMIESouth = 0 
+        nAMIESouth = 0
         call read_in_int(nAMIESouth, iError)
         if (iError == 0 .and. nAMIESouth > 0) then
           allocate(cAMIEListSouth(nAMIESouth))
@@ -776,7 +776,6 @@ subroutine set_inputs
           IsDone = .true.
         endif
 
-
       case ("#AURORAMODS")
         HasSetAuroraMods = .true.
         call read_in_logical(NormalizeAuroraToHP, iError)
@@ -967,6 +966,18 @@ subroutine set_inputs
           IsDone = .true.
         endif
 
+      case ("#NEUTRALLIMITS")
+        call read_in_real(MinNeutralDensity, iError)
+        call read_in_real(MinNeutralDensityAdvect, iError)
+        if (iError /= 0) then
+          write(*, *) 'Incorrect format for #NEUTRALLIMITS:'
+          write(*, *) ''
+          write(*, *) '#IONLIMITS'
+          write(*, *) "MinNeutralDensity    (real, default=100 m^-3)"
+          write(*, *) "MinNeutralDensityAdvect    (real, default=1e5 m^-3)"
+          IsDone = .true.
+        endif
+
       case ("#PHOTOELECTRON")
         if (RCMROutType /= 'PHOTOELECTRON') then
           call read_in_real(PhotoElectronHeatingEfficiency, iError)
@@ -1848,7 +1859,7 @@ subroutine set_inputs
           write(*, *) "----------------------------------------------"
           write(*, *) "GITM allows the use of a 10% seasonal "
           write(*, *) "correction factor in the conversion of AE to HP"
-          write(*, *) "Put another T after the onset file if you want this,"
+          write(*, *) "Put another T after convert to HP if you want this,"
           write(*, *) "or put an F to silence this message."
           write(*, *) "----------------------------------------------"
           iError = 0