Enforce nonnegative reaction quantities at moisture extinction

doc/Configuration.rst

@@ -292,6 +292,36 @@ output fields retain their existing single-precision NetCDF representation.
    Rate of Spread (ROS) parameterization option.
      0: Rothermel model (only option currently implemented)
 
+   For burnable cells, effective ground moisture ``fmc_g`` must be finite and
+   nonnegative, and fuel-table extinction moisture ``fuelmce`` must be finite
+   and positive. Both use kg water per kg dry fuel. Invalid inputs stop with
+   the cell indices or fuel category in the diagnostic. Fuel-table values are
+   checked during initialization, and prescribed ``fuelmc_g`` is checked in
+   the namelist (zero moisture is allowed). The moisture component validates
+   fixed averaging weights during initialization and checks computed classes
+   before averaging and updating spread parameters. The restart reader checks
+   restored class and effective moisture separately: ``fmc_g`` is read directly,
+   not recomputed from ``fmc_gc``, and prescribed-moisture restarts have no class
+   state. Spatial moisture checks run outside OpenMP tile regions and exclude
+   nonburnable cells.
+   Existing restart checks also require available class state and forcing history.
+
+   Aggregate moisture damping retains the existing polynomial below extinction
+   and is zero at or above extinction. Reaction intensity, ``iboros`` (reaction
+   intensity times residence time), and the Rothermel base spread rate are
+   therefore zero at extinction. The stored physical moisture is not capped.
+   The independent empirical chaparral spread law for Anderson category 4 is
+   unchanged; its reaction-based flame length is zero at or above extinction.
+
+   Restart files retain physical moisture. The reader reconstructs spread
+   parameters from that state using the current damping calculation, so no
+   new restart field or namelist option is required. Existing PR #39 files
+   remain readable when their inputs satisfy the same validation. A legacy
+   checkpoint can contain flame diagnostics calculated with the old damping;
+   these are replaced at the next fire advance. Readability does not imply
+   numerical equivalence across the extinction correction. Continuity must be
+   assessed with continuous and restarted runs of the same model revision.
+
 ``fmc_opt``: *integer* (Default: ``-1``)
    :term:`FMC` model
      -1 = Constant fuel moisture (only option currently implemented)

io/namelist_mod.F90

@@ -396,6 +396,9 @@ subroutine Check_nml (this, require_atm_interval)
 
       if (this%ideal_opt /= 0 .and. this%fmoist_run) &
           call Stop_simulation ('ideal runs do not support a FMC model')
+      ! Prescribed moisture may describe dry fuel, but must be a physical mass ratio.
+      if (.not. ieee_is_finite (this%fuelmc_g)) call Stop_simulation ('fuelmc_g must be finite')
+      if (this%fuelmc_g < 0.0) call Stop_simulation ('fuelmc_g must be nonnegative')
       if (this%fire_is_real_perim .and. this%fire_num_ignitions /= 1) &
           call Stop_simulation ('perimeter mode requires one ignition record for its ignition time')
       if (.not. this%fire_is_real_perim .and. this%fire_num_ignitions <= 0) &

physics/fmc_wrffire_mod.F90

@@ -55,6 +55,7 @@ module fmc_wrffire_mod
       procedure, public :: Advance_fmc_model => Advance_fmc_model
       procedure, public :: Advance_moisture_classes => Advance_moisture_classes
       procedure, public :: Average_moisture_classes => Average_moisture_classes
+      procedure, public :: Check_moisture_classes => Check_moisture_classes
       procedure, public :: Add_restart_vars => Add_restart_vars
       procedure, public :: Init => Init_fmc_wrffire
     end type fmc_wrffire_t
@@ -170,11 +171,20 @@ subroutine Init_fmc_wrffire (this, fuels, fuelmc_g, fuelmc_g_live, ifms, ifme, j
       this%fmc_gc_initial_value(4) = fuelmc_g
       this%fmc_gc_initial_value(5) = fuelmc_g_live
 
-        ! Calc averaging weights
+      ! Fixed weights must form a nonnegative, normalized mean for each loaded
+      ! fuel category. Validate once here; subsequent updates check class state.
       do i = 1, mfuelcats
         fgi_t(i) = 0.0
         do k = 1, MOISTURE_CLASSES
+          if (.not. ieee_is_finite (fgi_c(i, k))) then
+            write (msg, '(a,2(1x,i0))') 'Moisture averaging fuel load must be finite for category,class:', i, k
+            call Stop_simulation (msg)
+          end if
           if (fgi_c(i, k) >= 0.0) then
+            if (fgi_c(i, k) > huge (fgi_t(i)) - fgi_t(i)) then
+              write (msg, '(a,1x,i0)') 'Moisture averaging total fuel load exceeds representable range for category:', i
+              call Stop_simulation (msg)
+            end if
             fgi_t(i) = fgi_t(i) + fgi_c(i, k)
           else
             write (msg, *) 'fuel load in category', i, ' fuel class ', k, ' is ', fgi_c(i, k), ',must be nonegative.'
@@ -189,6 +199,17 @@ subroutine Init_fmc_wrffire (this, fuels, fuelmc_g, fuelmc_g_live, ifms, ifme, j
             this%fmc_gw(i, k) = 0.0
           end if
         end do
+        if (.not. all (ieee_is_finite (this%fmc_gw(i,:))) .or. any (this%fmc_gw(i,:) < 0.0)) then
+          write (msg, '(a,1x,i0)') 'Moisture averaging weights must be finite and nonnegative for category:', i
+          call Stop_simulation (msg)
+        end if
+        if (fgi_t(i) > 0.0) then
+          ! Allow rounding from summing MOISTURE_CLASSES default-real fractions.
+          if (abs (sum (this%fmc_gw(i,:)) - 1.0) > MOISTURE_CLASSES * epsilon (1.0)) then
+            write (msg, '(a,1x,i0)') 'Moisture averaging weights must sum to one for category:', i
+            call Stop_simulation (msg)
+          end if
+        end if
       end do
 
         ! moisture model derived scalars
@@ -260,6 +281,13 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms,
         !$OMP END PARALLEL DO
         this%initialized = .true.
 
+        ! Validate the completed update on the initial thread, before averaging.
+        ! Fatal input diagnostics must not race across OpenMP tiles in MPI_Abort.
+        do ij = 1, num_tiles
+          call this%Check_moisture_classes (ifms, ifme, jfms, jfme, i_start(ij), i_end(ij), j_start(ij), j_end(ij), &
+              nfuel_cat, fuels%no_fuel_cat, 'Computed')
+        end do
+
         !$OMP PARALLEL DO   &
         !$OMP PRIVATE (ij)
         do ij = 1, num_tiles
@@ -278,6 +306,35 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms,
 
     end subroutine Advance_fmc_model
 
+    subroutine Check_moisture_classes (this, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, nfuel_cat, no_fuel_cat, source)
+
+      implicit none
+
+      class (fmc_wrffire_t), intent (in) :: this
+      integer, intent (in) :: ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, no_fuel_cat
+      real, dimension(ifms:ifme, jfms:jfme), intent (in) :: nfuel_cat
+      character (len = *), intent (in) :: source
+      integer :: i, j, k
+      character (len = 160) :: msg
+
+
+      do j = jfts, jfte
+        do k = 1, MOISTURE_CLASSES
+          do i = ifts, ifte
+            if (int (nfuel_cat(i, j)) == no_fuel_cat) cycle
+            if (.not. ieee_is_finite (this%fmc_gc(i, k, j))) then
+              write (msg, '(a,a,3(1x,i0))') source, ' fmc_gc must be finite at i,class,j:', i, k, j
+              call Stop_simulation (msg)
+            else if (this%fmc_gc(i, k, j) < 0.0) then
+              write (msg, '(a,a,3(1x,i0))') source, ' fmc_gc must be nonnegative at i,class,j:', i, k, j
+              call Stop_simulation (msg)
+            end if
+          end do
+        end do
+      end do
+
+    end subroutine Check_moisture_classes
+
     subroutine Average_moisture_classes (this, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, nfuel_cat, fmc_g)
 
       implicit none

physics/fuel_anderson_mod.F90

@@ -117,6 +117,8 @@ subroutine Init_anderson_fuel_model (this, fuelmc_c)
         end if
       end do
 
+      call this%Check_fuels ()
+
       ! Initialize wind adjustment factor
       call this%Calc_wind_adjustment_factor()
 

physics/ros_wrffire_mod.F90

@@ -176,9 +176,16 @@ subroutine Set_ros_parameters_wrffire (this, ifms, ifme, jfms, jfme, ifts, ifte,
 
              ! net fuel loading, lb/ft^2
             wn = fuelload/(1 + fuels%st(k))
-            rtemp1 = fmc_g(i, j) / fuels%fuelmce(k)
               ! moist damp coef
-            etam = 1.0 - 2.59 * rtemp1 + 5.11 * rtemp1 ** 2 - 3.52 * rtemp1 ** 3
+              ! Zero reaction at extinction without capping physical moisture.
+              ! Polynomial cancellation can leave a positive residual at equality.
+              ! Retain the existing polynomial and arithmetic below extinction.
+            if (fmc_g(i, j) >= fuels%fuelmce(k)) then
+              etam = 0.0
+            else
+              rtemp1 = fmc_g(i, j) / fuels%fuelmce(k)
+              etam = 1.0 - 2.59 * rtemp1 + 5.11 * rtemp1 ** 2 - 3.52 * rtemp1 ** 3
+            end if
               ! mineral damping coef
             etas = 0.174 * fuels%se(k) ** (-0.19)
               !rxn intensity,btu/ft^2 min

share/fmc_mod.F90

@@ -1,7 +1,9 @@
   module fmc_mod
 
+    use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
     use fuel_mod, only : fuel_t
     use ros_mod, only : ros_t
+    use stderrout_mod, only : Stop_simulation
 
     implicit none
 
@@ -16,6 +18,7 @@ module fmc_mod
       procedure (Advance_fmc_model), deferred :: Advance_fmc_model
       procedure (Add_restart_vars), deferred :: Add_restart_vars
       procedure (Init), deferred :: Init
+      procedure, public :: Check_fmc => Check_fmc
     end type fmc_t
 
     abstract interface
@@ -57,4 +60,37 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms,
       end subroutine Advance_fmc_model
     end interface
 
+  contains
+
+    subroutine Check_fmc (this, ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, nfuel_cat, no_fuel_cat, fmc_g, source)
+
+      implicit none
+
+      class (fmc_t), intent (in) :: this
+      integer, intent (in) :: ifms, ifme, jfms, jfme, ifts, ifte, jfts, jfte, no_fuel_cat
+      real, dimension(ifms:ifme, jfms:jfme), intent (in) :: nfuel_cat, fmc_g
+      character (len = *), intent (in) :: source
+      integer :: i, j
+      character (len = 160) :: msg
+
+
+      ! Check computed or restored moisture before it reaches spread physics.
+      ! Prescribed moisture has an FMC object but no initialized class arrays;
+      ! this common check therefore uses only the supplied effective moisture.
+      ! Only owned burnable cells require moisture; halo values may be unavailable.
+      do j = jfts, jfte
+        do i = ifts, ifte
+          if (int (nfuel_cat(i, j)) == no_fuel_cat) cycle
+          if (.not. ieee_is_finite (fmc_g(i, j))) then
+            write (msg, '(a,a,2(1x,i0))') source, ' fmc_g must be finite at i,j:', i, j
+            call Stop_simulation (msg)
+          else if (fmc_g(i, j) < 0.0) then
+            write (msg, '(a,a,2(1x,i0))') source, ' fmc_g must be nonnegative at i,j:', i, j
+            call Stop_simulation (msg)
+          end if
+        end do
+      end do
+
+    end subroutine Check_fmc
+
   end module fmc_mod

share/fuel_mod.F90

@@ -1,5 +1,8 @@
   module fuel_mod
 
+    use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
+    use stderrout_mod, only : Stop_simulation
+
     implicit none
 
     private
@@ -40,6 +43,7 @@ module fuel_mod
       real, dimension(:), allocatable :: waf
     contains
       procedure (Initialization), deferred :: Initialization
+      procedure, public :: Check_fuels => Check_fuels
       procedure, public :: Calc_wind_adjustment_factor => Calc_wind_adjustment_factor
     end type fuel_t
 
@@ -53,6 +57,29 @@ end subroutine Initialization
 
   contains
 
+    subroutine Check_fuels (this)
+
+      implicit none
+
+      class (fuel_t), intent (in) :: this
+      integer :: k
+      character (len = 128) :: msg
+
+
+      ! Extinction moisture is a fixed table property, checked once at initialization.
+      do k = 1, this%n_fuel_cat
+        if (k == this%no_fuel_cat) cycle
+        if (.not. ieee_is_finite (this%fuelmce(k))) then
+          write (msg, '(a,1x,i0)') 'fuelmce must be finite for fuel category:', k
+          call Stop_simulation (msg)
+        else if (this%fuelmce(k) <= 0.0) then
+          write (msg, '(a,1x,i0)') 'fuelmce must be positive for fuel category:', k
+          call Stop_simulation (msg)
+        end if
+      end do
+
+    end subroutine Check_fuels
+
     pure function Crosswalk_from_scottburgan_to_anderson (fuel_cat) result (return_value)
 
       implicit none

state/state_restart_submod.F90

@@ -520,6 +520,11 @@ module subroutine Read_restart (this, config_flags)
 
       if (allocated (this%ros_param) .and. allocated (this%fuels)) then
         do ij = 1, this%num_tiles
+          ! fmc_g is read independently from the checkpoint, not recomputed from
+          ! fmc_gc. Valid classes cannot establish its validity, and prescribed
+          ! moisture restarts have no class state; check it before rebuilding ROS.
+          call this%fmc_param%Check_fmc (this%ifms, this%ifme, this%jfms, this%jfme, this%i_start(ij), this%i_end(ij), &
+              this%j_start(ij), this%j_end(ij), this%nfuel_cat, this%fuels%no_fuel_cat, this%fmc_g, 'Restart')
           call this%ros_param%Set_params (this%ifms, this%ifme, this%jfms, this%jfme, this%i_start(ij), this%i_end(ij), &
               this%j_start(ij), this%j_end(ij), this%fuels, this%nfuel_cat, this%fmc_g)
         end do
@@ -667,6 +672,9 @@ subroutine Read_restart_fmc (file_name)
                 this%ifps, this%ifpe, this%jfps, this%jfpe, NAME_VAR_FMC_GC, &
                 fmc_param%fmc_gc(this%ifps:this%ifpe, 1:n_moisture_classes, this%jfps:this%jfpe))
 
+            call fmc_param%Check_moisture_classes (this%ifms, this%ifme, this%jfms, this%jfme, &
+                this%ifps, this%ifpe, this%jfps, this%jfpe, this%nfuel_cat, this%fuels%no_fuel_cat, 'Restart')
+
             fmoist_lasttime = 0.0
             fmoist_nexttime = 0.0
             initialization_time = 0.0

tests/unit/CMakeLists.txt

@@ -22,6 +22,91 @@ foreach(test IN ITEMS tiles_unit tile_namelist_unit observed_perimeter_unit
   target_link_libraries(${test} PRIVATE firelib)
 endforeach()
 
+add_executable(moisture_extinction_unit moisture_extinction_unit.F90)
+target_link_libraries(moisture_extinction_unit PRIVATE firelib)
+foreach(scenario IN ITEMS below edge at above)
+  add_cfbm_unit_test(NAME moisture_extinction_${scenario}_unit
+    COMMAND moisture_extinction_unit ${scenario} LABELS fuel moisture ros
+    EXPECT "moisture extinction checks passed: ${scenario}")
+  if(OPENMP)
+    set_tests_properties(moisture_extinction_${scenario}_unit PROPERTIES
+      ENVIRONMENT "OMP_NUM_THREADS=2;OMP_DYNAMIC=FALSE")
+  endif()
+endforeach()
+add_cfbm_unit_test(NAME moisture_extinction_negative_moisture_unit
+  COMMAND moisture_extinction_unit negative_moisture LABELS fuel moisture ros
+  EXPECT "STOP:Restart fmc_g must be nonnegative")
+foreach(scenario IN ITEMS nan_moisture inf_moisture neginf_moisture)
+  add_cfbm_unit_test(NAME moisture_extinction_${scenario}_unit
+    COMMAND moisture_extinction_unit ${scenario} LABELS fuel moisture ros
+    EXPECT "STOP:Restart fmc_g must be finite")
+endforeach()
+foreach(scenario IN ITEMS zero_extinction negative_extinction)
+  add_cfbm_unit_test(NAME moisture_extinction_${scenario}_unit
+    COMMAND moisture_extinction_unit ${scenario} LABELS fuel moisture ros
+    EXPECT "STOP:fuelmce must be positive")
+endforeach()
+foreach(scenario IN ITEMS nan_extinction inf_extinction neginf_extinction)
+  add_cfbm_unit_test(NAME moisture_extinction_${scenario}_unit
+    COMMAND moisture_extinction_unit ${scenario} LABELS fuel moisture ros
+    EXPECT "STOP:fuelmce must be finite")
+endforeach()
+
+add_cfbm_unit_test(NAME namelist_fuelmc_zero_unit
+  COMMAND namelist_validation_unit fuelmc_zero LABELS namelist moisture
+  EXPECT "namelist validation accepted: fuelmc_zero")
+add_cfbm_unit_test(NAME namelist_fuelmc_negative_unit
+  COMMAND namelist_validation_unit fuelmc_negative LABELS namelist moisture
+  EXPECT "STOP:fuelmc_g must be nonnegative")
+foreach(scenario IN ITEMS nan inf neginf)
+  add_cfbm_unit_test(NAME namelist_fuelmc_${scenario}_unit
+    COMMAND namelist_validation_unit fuelmc_${scenario} LABELS namelist moisture
+    EXPECT "STOP:fuelmc_g must be finite")
+endforeach()
+
+add_executable(fmc_validation_unit fmc_validation_unit.F90)
+target_link_libraries(fmc_validation_unit PRIVATE firelib)
+foreach(scenario IN ITEMS valid zero_load restored_zero restored_wet)
+  add_cfbm_unit_test(NAME fmc_validation_${scenario}_unit
+    COMMAND fmc_validation_unit ${scenario} LABELS fuel moisture
+    EXPECT "moisture validation checks passed: ${scenario}")
+endforeach()
+foreach(quantity IN ITEMS class restored)
+  foreach(scenario IN ITEMS negative nan inf neginf)
+    if(quantity STREQUAL "class")
+      set(field fmc_gc)
+      set(source Computed)
+    else()
+      set(field fmc_g)
+      set(source Restart)
+    endif()
+    if(scenario STREQUAL "negative")
+      set(requirement nonnegative)
+    else()
+      set(requirement finite)
+    endif()
+    add_cfbm_unit_test(NAME fmc_validation_${quantity}_${scenario}_unit
+      COMMAND fmc_validation_unit ${quantity}_${scenario} LABELS fuel moisture
+      EXPECT "STOP:${source} ${field} must be ${requirement}")
+  endforeach()
+endforeach()
+add_cfbm_unit_test(NAME fmc_validation_load_negative_unit
+  COMMAND fmc_validation_unit load_negative LABELS fuel moisture
+  EXPECT "STOP:.*fuel load in category.*must be nonegative")
+foreach(scenario IN ITEMS nan inf neginf)
+  add_cfbm_unit_test(NAME fmc_validation_load_${scenario}_unit
+    COMMAND fmc_validation_unit load_${scenario} LABELS fuel moisture
+    EXPECT "STOP:Moisture averaging fuel load must be finite for category,class:")
+endforeach()
+add_cfbm_unit_test(NAME fmc_validation_load_overflow_unit
+  COMMAND fmc_validation_unit load_overflow LABELS fuel moisture
+  EXPECT "STOP:Moisture averaging total fuel load exceeds representable range for category:")
+if(OPENMP)
+  get_property(fmc_tests DIRECTORY PROPERTY TESTS)
+  list(FILTER fmc_tests INCLUDE REGEX "^fmc_validation_")
+  set_tests_properties(${fmc_tests} PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=2;OMP_DYNAMIC=FALSE")
+endif()
+
 # Check the production wind routine's array bounds even in Release builds.
 # Keep its test module separate from firelib's normal compilation products.
 target_sources(wind_order_unit PRIVATE ${PROJECT_SOURCE_DIR}/io/coupling_mod.F90)