From fa72b53562f86b89f4d32cfbfabe9e6f0a9119fc Mon Sep 17 00:00:00 2001 From: scivision Date: Thu, 30 Jul 2026 12:06:48 -0400 Subject: [PATCH 01/11] MSIS2: consistency lint --- cmake/libraries.json | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/libraries.json b/cmake/libraries.json index 35d1c25f..9330caaf 100644 --- a/cmake/libraries.json +++ b/cmake/libraries.json @@ -17,10 +17,10 @@ }, "msis": { "stem": "https://github.com/gemini3d/msis/archive/", - "archive": "1577d765f86f4bb180e589236be3336d1dbbc26f.tar.gz" + "archive": "586c9b96931d80011db417bcb34140580258e786.tar.gz" }, "mumps": { - "stem": "https://github.com/scivision/mumps/archive/", + "stem": "https://github.com/scivision/mumps-superbuild/archive/", "archive": "595ad9b758411da45a2cfd272e2164a55e7eae11.tar.gz" } } From a4c1d7055da8c2480f21f3f4363bee7b8a3f1e68 Mon Sep 17 00:00:00 2001 From: scivision Date: Fri, 31 Jul 2026 20:19:59 -0400 Subject: [PATCH 02/11] h5fortran: FindHDF5 programmatic support library search. HDF5 2.2.0 release --- cmake/libraries.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/libraries.json b/cmake/libraries.json index 9330caaf..a4fe4419 100644 --- a/cmake/libraries.json +++ b/cmake/libraries.json @@ -9,7 +9,7 @@ }, "h5fortran": { "stem": "https://github.com/ECLAIRWaveS/h5fortran/archive/", - "archive": "641a527b1aaf5a9503e7324481c79503e091a6c4.tar.gz" + "archive": "5db3c2efc6425a2219cfb3c1365fcc206b807cfd.tar.gz" }, "hwm14": { "stem": "https://github.com/gemini3d/hwm14/archive/", From 65e13e789f8cd5bb1352d71be6ceaa86eb8f60e9 Mon Sep 17 00:00:00 2001 From: scivision Date: Mon, 3 Aug 2026 11:49:06 -0400 Subject: [PATCH 03/11] add gemini3d_msis2 option to enable MSIS2.x Must also specifiy MSIS 2.x if desired in config.nml --- CMakeLists.txt | 3 +++ options.cmake | 2 ++ 2 files changed, 5 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1ee940a0..8b0f327c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -138,6 +138,9 @@ if(gemini3d_hwm14) endif() # --- MSISE00 / MSIS 2.x +if(gemini3d_msis2) + set(msis2 on) +endif() set(msis_BUILD_UTILS on) string(JSON _stem GET "${json}" "msis" "stem") string(JSON _archive GET "${json}" "msis" "archive") diff --git a/options.cmake b/options.cmake index 844a59fd..66d741ad 100644 --- a/options.cmake +++ b/options.cmake @@ -30,6 +30,8 @@ option(gemini3d_python "Python-based self-checks") # Matlab checks take much longer than Python, and Python covers much more option(gemini3d_matlab "Matlab-based self-checks") +option(gemini3d_msis2 "enable MSIS2 neutral atmosphere model (must specify in config.nml as well)" on) + # append .debug to debug libraries, because the computation speed penalty is so great set(CMAKE_DEBUG_POSTFIX .debug) From 903f8c149a17bfba74c548b77246e756bd591456 Mon Sep 17 00:00:00 2001 From: scivision Date: Mon, 3 Aug 2026 12:11:42 -0400 Subject: [PATCH 04/11] msis_setup: -features outputs MSIS2 if present, and puts msis_setup in top binary dir --- CMakeLists.txt | 1 + cmake/libraries.json | 2 +- cmake/summary.cmake | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8b0f327c..29f0f8fb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -147,6 +147,7 @@ string(JSON _archive GET "${json}" "msis" "archive") set(msis_url "${_stem}${_archive}") FetchContent_Declare(msis URL ${msis_url}) FetchContent_MakeAvailable(msis) +set_target_properties(msis_setup PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}) # --- Gemini3D library add_subdirectory(src) diff --git a/cmake/libraries.json b/cmake/libraries.json index a4fe4419..4754fbfa 100644 --- a/cmake/libraries.json +++ b/cmake/libraries.json @@ -17,7 +17,7 @@ }, "msis": { "stem": "https://github.com/gemini3d/msis/archive/", - "archive": "586c9b96931d80011db417bcb34140580258e786.tar.gz" + "archive": "e8437eb61ff0bdb3ef336d96fd279dd75d9aec98.tar.gz" }, "mumps": { "stem": "https://github.com/scivision/mumps-superbuild/archive/", diff --git a/cmake/summary.cmake b/cmake/summary.cmake index f1e4ac68..01df7ab8 100644 --- a/cmake/summary.cmake +++ b/cmake/summary.cmake @@ -4,6 +4,7 @@ if(gemini3d_IS_TOP_LEVEL) add_feature_info(GLOW gemini3d_glow "airglow / aurora model") add_feature_info(HWM14 gemini3d_hwm14 "HWM14 neutral winds model") +add_feature_info(MSIS2 gemini3d_msis2 "enable MSIS 2.x neutral atmosphere model") add_feature_info(PyGemini gemini3d_python "simulation generation, HPC script generator and plotting") add_feature_info(MatGemini gemini3d_matlab "checks not as extensive as Python, and slow") From e01ce60bc3ab39934fa6faafa78e4eda1de064ea Mon Sep 17 00:00:00 2001 From: scivision Date: Mon, 3 Aug 2026 16:58:30 -0400 Subject: [PATCH 05/11] MSIS2 needs msis21.parm copied to gemini.bin directory --- CMakeLists.txt | 4 ++++ test/sim.cmake | 3 +++ 2 files changed, 7 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 29f0f8fb..ab7d4901 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -148,6 +148,10 @@ set(msis_url "${_stem}${_archive}") FetchContent_Declare(msis URL ${msis_url}) FetchContent_MakeAvailable(msis) set_target_properties(msis_setup PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}) +if(gemini3d_msis2) + set(msis2_parm_file ${msis_SOURCE_DIR}/src/msis2/msis21.parm) + configure_file(${msis2_parm_file} msis21.parm COPYONLY) +endif() # --- Gemini3D library add_subdirectory(src) diff --git a/test/sim.cmake b/test/sim.cmake index 3aea4046..7436d95a 100644 --- a/test/sim.cmake +++ b/test/sim.cmake @@ -24,6 +24,9 @@ else() endif() endif() +if(gemini3d_msis2) + list(APPEND gemini_test_sim_names mini2dns_msis2_fang) +endif() foreach(_s IN LISTS gemini_test_sim_names) setup_gemini_test(${_s}) From d9c65438fb54dfeb66400fbae62a9569972eef31 Mon Sep 17 00:00:00 2001 From: scivision Date: Mon, 3 Aug 2026 17:24:32 -0400 Subject: [PATCH 06/11] ci: add msis 2.x test --- .github/workflows/ci.yml | 20 ++++++++++++++++++- CMakePresets.json | 43 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 61 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e6449485..054afb7e 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -30,7 +30,6 @@ jobs: strategy: fail-fast: false matrix: - hwm14: [false] os: [ubuntu-latest] gcc: [11, 12, 13, 14] @@ -58,6 +57,25 @@ jobs: - &linux-build uses: ./.github/workflows/composite-unix + hwm-msis: + strategy: + matrix: + model: [msis2] + + env: + CMAKE: cmake + CTEST: ctest + + runs-on: ubuntu-latest + + steps: + - *checkout + + - *pkg + + - name: workflow ${{ matrix.model }} + run: cmake --workflow ${{ matrix.model }} + cmake-oldest: timeout-minutes: 30 diff --git a/CMakePresets.json b/CMakePresets.json index 35adcde8..5ec528e4 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -16,7 +16,20 @@ { "name": "gmake", "inherits": "default", "generator": "Unix Makefiles", "binaryDir": "build-${presetName}", "displayName": "GNU generator", - "description": "this helps catch GNU Make-specific bugs that might not appear with Ninja" }, + "description": "this helps catch GNU Make-specific bugs that might not appear with Ninja" +}, +{ "name": "msis2", "inherits": "default", + "displayName": "Use MSIS2 model", + "cacheVariables": { + "gemini3d_msis2": true + } +}, +{ "name": "hwm14", "inherits": "default", + "displayName": "Use HWM14 model", + "cacheVariables": { + "gemini3d_hwm14": true + } +}, { "name": "linux", "inherits": "default", "displayName": "Linux build (for WSL on GHA)", @@ -54,6 +67,8 @@ "buildPresets": [ { "name": "default", "configurePreset": "default", "configuration": "Release" }, { "name": "no-build-hdf5", "configurePreset": "no-build-hdf5", "inherits": "default" }, + { "name": "msis2", "configurePreset": "msis2", "inherits": "default" }, + { "name": "hwm14", "configurePreset": "hwm14", "inherits": "default" }, { "name": "gmake", "inherits": "default", "configurePreset": "gmake", "jobs": 4 }, { "name": "linux", "configurePreset": "linux", "inherits": "default" }, { "name": "release", "inherits": "default", "configurePreset": "release" }, @@ -75,6 +90,18 @@ "stopOnFailure": false } }, +{ "name": "msis2", "configurePreset": "msis2", "inherits": "default", + "displayName": "Tests for MSIS 2.x", + "filter": { + "include": {"name": "msis2"} + } +}, +{ "name": "hwm14", "configurePreset": "hwm14", "inherits": "default", + "displayName": "Tests for HWM14", + "filter": { + "include": {"name": "hwm14"} + } +}, { "name": "download", "inherits": "default", "displayName": "download reference data for subsequent offline tests", @@ -148,6 +175,20 @@ { "type": "build", "name": "default" } ] }, + { "name": "msis2", "displayName": "MSIS 2.x workflow", + "steps": [ + { "type": "configure", "name": "msis2" }, + { "type": "build", "name": "msis2" }, + { "type": "test", "name": "msis2" } + ] + }, + { "name": "hwm14", "displayName": "HWM14 workflow", + "steps": [ + { "type": "configure", "name": "hwm14" }, + { "type": "build", "name": "hwm14" }, + { "type": "test", "name": "hwm14" } + ] + }, { "name": "no-build-hdf5", "displayName": "Only find, don't build HDF5 library from source if not found", "steps": [ From a1356b21676885a96c4f3a2e3d34dc1194fde007 Mon Sep 17 00:00:00 2001 From: zettergm Date: Mon, 3 Aug 2026 18:38:00 -0400 Subject: [PATCH 07/11] add C interface for alternate energy solution method --- include/gemini3d.h | 1 + src/libgemini_c.f90 | 28 +++++++++++++++++++++++++++- 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/include/gemini3d.h b/include/gemini3d.h index 2012ca91..5b7066a7 100644 --- a/include/gemini3d.h +++ b/include/gemini3d.h @@ -85,6 +85,7 @@ extern void VNRicht_artvisc_C(double**, void**); extern void compression_C(double**, double**, void**, int*, void**, double*); extern void rhoe2T_C(double**, double**); extern void energy_diffusion_C(void**, int*, void**, double**, double**, void**, double*); +extern void diffusion_source_loss_energy_C(void**, int*, void**, double**, double**, void**, double*); extern void source_loss_allparams_C(void**, double**, double**, double**, void**, int*, void**, double*); extern void source_loss_mass_C(void**, double**,double**,double**,void**,int*,void**,double*); diff --git a/src/libgemini_c.f90 b/src/libgemini_c.f90 index c7512916..317eb58f 100644 --- a/src/libgemini_c.f90 +++ b/src/libgemini_c.f90 @@ -42,7 +42,7 @@ module gemini3d_C sweep1_allspec_mass_in,sweep1_allspec_momentum_in,sweep1_allspec_energy_in, & sweep2_allspec_mass_in,sweep2_allspec_momentum_in,sweep2_allspec_energy_in, & rhov12v1_in, VNRicht_artvisc_in, compression_in, rhoe2T_in, clean_param_in, & - energy_diffusion_in, source_loss_allparams_in, & + energy_diffusion_in, diffusion_source_loss_energy_in, source_loss_allparams_in, & source_loss_mass_in, source_loss_momentum_in, source_loss_energy_in, & clear_ionization_arrays, impact_ionization_in, solar_ionization_in, & dateinc_in, get_subgrid_size,get_fullgrid_size,get_config_vars, get_species_size, fluidvar_pointers, & @@ -935,6 +935,32 @@ subroutine energy_diffusion_C(cfgC,xtype,xC,fluidvarsC,electrovarsC,intvarsC,dt) end subroutine energy_diffusion_C + !> diffusion of energy + subroutine diffusion_source_loss_energy_C(cfgC,xtype,xC,fluidvarsC,electrovarsC,intvarsC,dt) & + bind(C,name="diffusion_source_loss_energy_C") + type(c_ptr), intent(in) :: cfgC + integer(C_INT), intent(in) :: xtype + type(c_ptr), intent(in) :: xC + type(c_ptr), intent(inout) :: fluidvarsC + type(c_ptr), intent(in) :: electrovarsC + type(c_ptr), intent(in) :: intvarsC + real(wp), intent(in) :: dt + + type(gemini_cfg), pointer :: cfg + class(curvmesh), pointer :: x + real(wp), dimension(:,:,:,:), pointer :: fluidvars + real(wp), dimension(:,:,:,:), pointer :: electrovars + type(gemini_work), pointer :: intvars + + call c_f_pointer(cfgC, cfg) + x=>set_gridpointer_dyntype(xtype, xC) + call c_f_pointer(fluidvarsC,fluidvars,[(lx1+4),(lx2+4),(lx3+4),(5*lsp)]) + call c_f_pointer(electrovarsC,electrovars,[(lx1+4),(lx2+4),(lx3+4),7]) + call c_f_pointer(intvarsC,intvars) + call diffusion_source_loss_energy_in(cfg,x,fluidvars,electrovars,intvars,dt) + end subroutine diffusion_source_loss_energy_C + + !> source/loss numerical solutions subroutine source_loss_allparams_C(cfgC,fluidvarsC,fluidauxvarsC,electrovarsC,intvarsC,xtype,xC,dt) & bind(C, name="source_loss_allparams_C") From c32bdcbefc54a15f2838d1cb6544c3f91ce4e2c0 Mon Sep 17 00:00:00 2001 From: zettergm Date: Thu, 6 Aug 2026 19:11:37 -0400 Subject: [PATCH 08/11] revert altnull to default to 80km --- src/numerical/mesh/meshobj.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/numerical/mesh/meshobj.f90 b/src/numerical/mesh/meshobj.f90 index a7d34861..d6576b6f 100644 --- a/src/numerical/mesh/meshobj.f90 +++ b/src/numerical/mesh/meshobj.f90 @@ -159,7 +159,7 @@ module meshobj !> contains information about where to assign null cells, most uses will assume default value - real(wp) :: altnull=0 !80e3 + real(wp) :: altnull=80e3 !80e3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! type-bound procedures !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! contains From 477889779c4d00eac4fb8d992638c9e64aa3052c Mon Sep 17 00:00:00 2001 From: zettergm Date: Mon, 10 Aug 2026 17:56:42 -0400 Subject: [PATCH 09/11] assume flagnodivJ0 wants background currents to cancel out --- src/numerical/potential/potential_comm_mumps.f90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/numerical/potential/potential_comm_mumps.f90 b/src/numerical/potential/potential_comm_mumps.f90 index fc48ae91..ee1b6886 100644 --- a/src/numerical/potential/potential_comm_mumps.f90 +++ b/src/numerical/potential/potential_comm_mumps.f90 @@ -425,12 +425,20 @@ subroutine potential_sourceterms(incap,sigP,sigH,sigPgrav,sigHgrav,E02,E03,vn2,v J2 = 0 J3 = 0 !! zero everything out to initialize since *accumulating* sources + if (.not. cfg%flagnodivJ0) then call acc_perpBGconductioncurrents(sigP,sigH,E02,E03,J2,J3) !background conduction currents only if (debug .and. mpi_cfg%myid==0) print *, 'Workers have computed background field currents...' + call acc_perpwindcurrents(sigP,sigH,vn2,vn3,B1,J2,J3) ! always include wind effects + if (debug .and. mpi_cfg%myid==0) print *, 'Workers have computed wind currents...' + !^ flagnodivJ0 basically assumes that the background electric field being provided to the code is exactly + ! that needed to balance out the current denstiy from the winds. I.e. there is no need to accumulate + ! current involved in the establishment of the background state from these sources since they naturally + ! yield zero divergence. Or at least that is what we explicitly enforce with this flag. This is relevant + ! especially for EPB simulations where the full region over which the background state is established is not + ! resolved by the model and artificial potentials will form, thus, if one just computes currents explicitly, + ! subject to some boundary conditions that do not include explicity accounting for currents entering and leaving end if - call acc_perpwindcurrents(sigP,sigH,vn2,vn3,B1,J2,J3) ! always include wind effects - if (debug .and. mpi_cfg%myid==0) print *, 'Workers have computed wind currents...' if (cfg%flagdiamagnetic) then call acc_pressurecurrents(muP,muH,ns,Ts,x,J2,J3) if (debug .and. mpi_cfg%myid==0) print *, 'Workers have computed pressure currents...' From 47d4b312a6b062980f6a8bf3c851560e19198a3c Mon Sep 17 00:00:00 2001 From: zettergm Date: Tue, 11 Aug 2026 19:20:23 -0400 Subject: [PATCH 10/11] user control of density fill values; update echo printing of cfg --- docs/Readme_input.md | 8 ++ src/gemini_init.f90 | 119 ++++++++++++++++----- src/io/config_nml.f90 | 68 +++++++----- src/numerical/constants/phys_consts.in.f90 | 9 +- 4 files changed, 149 insertions(+), 55 deletions(-) diff --git a/docs/Readme_input.md b/docs/Readme_input.md index ae89e98f..ca0c31d8 100644 --- a/docs/Readme_input.md +++ b/docs/Readme_input.md @@ -233,6 +233,14 @@ flagJ1ve=.true. ! use parallel current density in calculation of parallel &nightQ flagnightQ=.true. ! use newer photoionization calculation including nighttime Q / + +! Set fill and min values for densities +! (optional - defaults to reasonable values for most applications) +&mindens_user + mindens_userval=1.0e-100 + mindensnull_userval=1.0e-20 + mindensdiv_userval=1.0e-5 +/ ``` diff --git a/src/gemini_init.f90 b/src/gemini_init.f90 index 7706a35a..51355d73 100644 --- a/src/gemini_init.f90 +++ b/src/gemini_init.f90 @@ -3,6 +3,7 @@ module gemini_init use gemini3d_config, only : gemini_cfg use mpimod, only : mpi_cfg use filesystem, only : assert_is_file, assert_is_dir +use phys_consts, only: mindens, mindensnull, mindensdiv implicit none (type, external) @@ -45,13 +46,25 @@ subroutine check_input_files(cfg) call assert_is_file(cfg%indatfile) print *, '******************** input config ****************' + + ! base print '(A)', 'simulation directory: ' // cfg%outdir print '(A51,I6,A1,I0.2,A1,I0.2)', ' start year-month-day: ', cfg%ymd0(1), '-', cfg%ymd0(2),'-', cfg%ymd0(3) print '(A51,F10.3)', 'start time: ',cfg%UTsec0 print '(A51,F10.3)', 'duration: ',cfg%tdur print '(A51,F10.3)', 'output every: ',cfg%dtout + print*, 'F10.7 and geomagnetic indices: ',cfg%activ + print*, 'Top boundary electron temperature: ',cfg%Teinf + + ! flags + print*, 'Potential solve: ', cfg%potsolve + print*, 'Periodic: ', cfg%flagperiodic + print*, 'Output type: ', cfg%flagoutput + + ! files print '(A,/,A,/,A,/,A)', 'gemini.f90: using input data files:', cfg%indatsize, cfg%indatgrid, cfg%indatfile + ! neutral perturb if(cfg%flagdneu==1) then if (.not. (cfg%interptype==5 .or. cfg%interptype==6)) call assert_is_dir(cfg%sourcedir) print *, 'Neutral disturbance mlat,mlon: ',cfg%sourcemlat,cfg%sourcemlon @@ -62,6 +75,7 @@ subroutine check_input_files(cfg) print *, "no neutral disturbance specified." end if + ! precip if (cfg%flagprecfile==1) then call assert_is_dir(cfg%precdir) print '(A,F10.3)', 'Precipitation file input cadence (s): ',cfg%dtprec @@ -70,6 +84,7 @@ subroutine check_input_files(cfg) print *, "no precipitation specified" end if + ! efield if(cfg%flagE0file==1) then call assert_is_dir(cfg%E0dir) print *, 'Electric field file input cadence (s): ',cfg%dtE0 @@ -78,6 +93,7 @@ subroutine check_input_files(cfg) print *, "no Efield specified" end if + ! solflux if(cfg%flagsolfluxfile==1) then call assert_is_dir(cfg%solfluxdir) print *, 'Solar flux file input cadence (s): ',cfg%dtsolflux @@ -86,7 +102,19 @@ subroutine check_input_files(cfg) print *, "no solar flux specified" end if + ! neutral_BG + if (cfg%flagneuBG) then + print*, 'Neutral background updated at: ',cfg%dtneuBG + print*, 'Using MSIS version: ', cfg%msis_version + end if + + ! neutral_BG file + if (cfg%flagneutralBGfile==1) then + print*, 'Neutral background file input cadence (s): ', cfg%dtneuBGfile + print*, 'Neutral background file input directory: ',cfg%neutralBGdir + end if + ! glow if (cfg%flagglow==1) then print *, 'GLOW enabled for auroral emission calculations.' print *, 'GLOW electron transport calculation cadence (s): ', cfg%dtglow @@ -95,88 +123,125 @@ subroutine check_input_files(cfg) print *, "GLOW disabled" end if - if (cfg%msis_version > 0) then - print '(A,f3.1,A)', 'MSIS ', real(cfg%msis_version)/10, ' enabled for neutral atmosphere calculations.' - else - print '(A)', "MSISE00 enabled for neutral atmosphere calculations." - end if - - if (cfg%flagneuBG) then - print*, 'Neutral background updated at: ',cfg%dtneuBG - end if - + ! EIA if (cfg%flagEIA) then print*, 'EIA enables with peok equatorial drift: ',cfg%v0equator else print*, 'EIA disabled' end if - if (cfg%flagneuBG) then - print*, 'Variable background neutral atmosphere enabled at cadence: ',cfg%dtneuBG - else - print*, 'Variable background neutral atmosphere disabled.' - end if - + ! precip_BG print*, 'Background precipitation has total energy flux and energy: ',cfg%PhiWBG,cfg%W0BG + ! Jpar if (cfg%flagJpar) then print*, 'Parallel current calculation enabled.' else print*, 'Parallel current calculation disabled.' end if + ! capacitance print*, 'Inertial capacitance calculation type: ',cfg%flagcap + ! diffusion print*, 'Diffusion solve type: ',cfg%diffsolvetype + ! milestone if (cfg%mcadence > 0) then print*, 'Milestone output selected; cadence (every nth output) of: ',cfg%mcadence else print*, 'Milestone output disabled.' end if + ! gravdrift if (cfg%flaggravdrift) then print*, 'Gravitational drift terms enabled.' else print*, 'Gravitaional drift terms disabled.' end if - if (cfg%flagtwoway) then - print*, 'Two-way coupling enabled.' - else - print*, 'Two-way coupling disabled.' - end if - + ! lagrangian if (cfg%flaglagrangian) then print*, 'Lagrangian grid enabled.' else print*, 'Lagrangian grid disabled' end if - if (cfg%flagmagpole) then - print*, 'Using year-based magnetic pole.' + ! diamagnetic + if (cfg%flagdiamagnetic) then + print*, 'Diamagnetic drift terms enabled.' else - print*, 'Using default magnetic pole.' + print*, 'Diamagnetic drift terms disabled.' + end if + + ! twoway_coupled + if (cfg%flagtwoway) then + print*, 'Two-way coupling enabled.' + else + print*, 'Two-way coupling disabled.' end if + ! nodivJ0 + if (cfg%flagnodivJ0) then + print*, 'Excluding background current divergence.' + else + print*, 'Including background current divergence.' + end if + + ! FBI if (cfg%flagFBI>=1) then print*, 'Anomalous electrojet heating enabled.' else print*, 'Anomalous electrojet heating disabled.' end if - if (cfg%flagFBI==2) then print*, 'Nonlinear current enabled.' else print*, 'Nonlinear current disabled.' end if + ! evibcool if (cfg%flagevibcool==1) then print*, 'Using updated model for electron inelastic collisions.' else print*, 'Using legacy model for electron inelastic collisions.' end if - + + ! magpole + if (cfg%flagmagpole) then + print*, 'Using year-based magnetic pole.' + else + print*, 'Using default magnetic pole.' + end if + + ! J1ve + if (cfg%flagJ1ve) then + print*, 'Computing parallel electron drift from parallel current density.' + else + print*, 'Computing parallel electron drift from ambipolar assumption.' + end if + + ! nightQ + if (cfg%flagnightQ) then + print*, 'Using nighttime ionization calculation.' + else + print*, 'Neglecting nighttime ionization.' + end if + + ! fang + print*, 'Fang impact ionization version: ',cfg%flag_fang + + ! fang_pars + if (cfg%flag_fang==0) then + print*, 'diff_num_flux: ',cfg%diff_num_flux + print*, 'kappa: ', cfg%kappa + print*, 'bimax_frac: ', cfg%bimax_frac + print*, 'W0_char: ', cfg%W0_char + end if + + ! mindens_user + print*, 'Density fill values: ',mindens,mindensnull,mindensdiv + print *, '**************** end input config ***************' end if diff --git a/src/io/config_nml.f90 b/src/io/config_nml.f90 index 893d67a6..06036c2c 100644 --- a/src/io/config_nml.f90 +++ b/src/io/config_nml.f90 @@ -3,6 +3,7 @@ use, intrinsic :: iso_fortran_env, only : stderr => error_unit use gemini3d_sysinfo, only : expand_envvar, get_compiler_vendor use filesystem, only : absolute +use phys_consts, only: mindens, mindensnull, mindensdiv implicit none (type, external) @@ -72,6 +73,11 @@ ! add nightime ionization logical :: flagnightQ = .false. + ! in case the user wants to specify minimum allowed density + real(wp) :: mindens_userval=1.0e-100 + real(wp) :: mindensnull_userval=1.0e-20 + real(wp) :: mindensdiv_userval=1.0e-5 + namelist /base/ ymd, UTsec0, tdur, dtout, activ, tcfl, Teinf namelist /files/ file_format, indat_size, indat_grid, indat_file namelist /flags/ potsolve, flagperiodic, flagoutput @@ -100,6 +106,7 @@ namelist /magpole/ flagmagpole namelist /J1ve/ flagJ1ve namelist /nightQ/ flagnightQ + namelist /mindens_user/ mindens_userval, mindensnull_userval, mindensdiv_userval if(.not. allocated(cfg%outdir)) error stop 'gemini3d:config:config_nml please specify simulation output directory' if(.not. allocated(cfg%infile)) error stop 'gemini3d:config:config_nml please specify simulation configuration file config.nml' @@ -197,6 +204,19 @@ cfg%solfluxdir = "" endif + !> neural background (optional) + if (namelist_exists(u,'neutral_BG')) then + rewind(u) + read(u, nml=neutral_BG, iostat=i) + call check_nml_io(i, cfg%infile, "neutral_BG") + cfg%flagneuBG=flagneuBG + cfg%dtneuBG=dtneuBG + cfg%msis_version = msis_version + else + cfg%flagneuBG=.false. + cfg%msis_version = 0 + end if + if (namelist_exists(u, "neutralBG_file", verbose)) then cfg%flagneutralBGfile = 1 rewind(u) @@ -209,6 +229,17 @@ cfg%neutralBGdir = "" endif + if (namelist_exists(u, "glow", verbose)) then + cfg%flagglow = 1 + rewind(u) + read(u, nml=glow, iostat=i) + call check_nml_io(i, cfg%infile, "glow") + cfg%dtglow = dtglow + cfg%dtglowout = dtglowout + else + cfg%flagglow = 0 + endif + if (namelist_exists(u, "fang", verbose)) then rewind(u) read(u, nml=fang, iostat=i) @@ -234,17 +265,6 @@ cfg%W0_char = 3000._wp ! same as W0BG default endif - if (namelist_exists(u, "glow", verbose)) then - cfg%flagglow = 1 - rewind(u) - read(u, nml=glow, iostat=i) - call check_nml_io(i, cfg%infile, "glow") - cfg%dtglow = dtglow - cfg%dtglowout = dtglowout - else - cfg%flagglow = 0 - endif - !> EIA (optional) if (namelist_exists(u,'EIA')) then rewind(u) @@ -256,19 +276,6 @@ cfg%flagEIA=.false. end if - !> neural background (optional) - if (namelist_exists(u,'neutral_BG')) then - rewind(u) - read(u, nml=neutral_BG, iostat=i) - call check_nml_io(i, cfg%infile, "neutral_BG") - cfg%flagneuBG=flagneuBG - cfg%dtneuBG=dtneuBG - cfg%msis_version = msis_version - else - cfg%flagneuBG=.false. - cfg%msis_version = 0 - end if - !> precip background (optional) if (namelist_exists(u,'precip_BG')) then rewind(u) @@ -416,6 +423,19 @@ cfg%flagnightQ = .false. ! not adding nighttime ionization (default uses the older version) end if + if (namelist_exists(u, 'mindens_user')) then + rewind(u) + read(u, nml=mindens_user, iostat=i) + call check_nml_io(i, cfg%infile, "mindens_user") + mindens = mindens_userval ! this is different from the others since we just directly set the module variable, rather than cfg + mindensnull = mindensnull_userval + mindensdiv = mindensdiv_userval + else + mindens = 1.0e-100 + mindensnull = 1.0e-20_wp + mindensdiv = 1.0e-5_wp + end if + close(u) end procedure read_nml diff --git a/src/numerical/constants/phys_consts.in.f90 b/src/numerical/constants/phys_consts.in.f90 index c52b4ae7..f2f791ab 100644 --- a/src/numerical/constants/phys_consts.in.f90 +++ b/src/numerical/constants/phys_consts.in.f90 @@ -84,10 +84,11 @@ module phys_consts 1304] !wavelength of each auroral line, housekeeping -!> HOUSEKEEPING PARAMETERS for conditioning densities -real(wp), parameter :: mindens = 1.0e-100_wp -real(wp), parameter :: mindensnull = 1.0e-20_wp -real(wp), parameter :: mindensdiv = 1.0e-5_wp +!> HOUSEKEEPING PARAMETERS for conditioning densities, these are writeable and may be adjusted in io +! module if the user specifies them in the nml input file. +real(wp) :: mindens = 1.0e-100_wp +real(wp) :: mindensnull = 1.0e-20_wp +real(wp) :: mindensdiv = 1.0e-5_wp !To control the amount of console output; can be changed by user command line flag "-debug" From dbe9ea9a14133677a083c941de475a73ee6a2403 Mon Sep 17 00:00:00 2001 From: zettergm Date: Tue, 11 Aug 2026 19:44:13 -0400 Subject: [PATCH 11/11] left out some changes from previous commit --- src/io/config_nml.f90 | 50 +++++++++++++++++++++---------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/src/io/config_nml.f90 b/src/io/config_nml.f90 index 06036c2c..b68d5ed1 100644 --- a/src/io/config_nml.f90 +++ b/src/io/config_nml.f90 @@ -240,31 +240,6 @@ cfg%flagglow = 0 endif - if (namelist_exists(u, "fang", verbose)) then - rewind(u) - read(u, nml=fang, iostat=i) - call check_nml_io(i, cfg%infile, "fang") - cfg%flag_fang = flag_fang - else - cfg%flag_fang = 2008 !< legacy default - endif - - if (namelist_exists(u, "fang_pars", verbose)) then - rewind(u) - read(u, nml=fang_pars, iostat=i) - call check_nml_io(i, cfg%infile, "fang_pars") - cfg%flag_fang = 0 ! force fang flag for integrated spectrum - cfg%diff_num_flux = diff_num_flux - cfg%kappa = kappa - cfg%bimax_frac = bimax_frac - cfg%W0_char = W0_char - else - cfg%diff_num_flux = 0 ! Maxwellian, same as Fang et al. 2008 within 5% in most cases - cfg%kappa = 1e4_wp ! close to Maxwellian - cfg%bimax_frac = 1._wp ! Maxwellian - cfg%W0_char = 3000._wp ! same as W0BG default - endif - !> EIA (optional) if (namelist_exists(u,'EIA')) then rewind(u) @@ -423,6 +398,31 @@ cfg%flagnightQ = .false. ! not adding nighttime ionization (default uses the older version) end if + if (namelist_exists(u, "fang", verbose)) then + rewind(u) + read(u, nml=fang, iostat=i) + call check_nml_io(i, cfg%infile, "fang") + cfg%flag_fang = flag_fang + else + cfg%flag_fang = 2008 !< legacy default + endif + + if (namelist_exists(u, "fang_pars", verbose)) then + rewind(u) + read(u, nml=fang_pars, iostat=i) + call check_nml_io(i, cfg%infile, "fang_pars") + cfg%flag_fang = 0 ! force fang flag for integrated spectrum + cfg%diff_num_flux = diff_num_flux + cfg%kappa = kappa + cfg%bimax_frac = bimax_frac + cfg%W0_char = W0_char + else + cfg%diff_num_flux = 0 ! Maxwellian, same as Fang et al. 2008 within 5% in most cases + cfg%kappa = 1e4_wp ! close to Maxwellian + cfg%bimax_frac = 1._wp ! Maxwellian + cfg%W0_char = 3000._wp ! same as W0BG default + endif + if (namelist_exists(u, 'mindens_user')) then rewind(u) read(u, nml=mindens_user, iostat=i)