From 138f9d5f327fdeb860e3ff6fde9fbaca9b93ed7c Mon Sep 17 00:00:00 2001 From: Maria Frediani Date: Wed, 21 Oct 2026 10:54:19 -0600 Subject: [PATCH 1/2] Test reaction and flame quantities across moisture extinction --- tests/unit/CMakeLists.txt | 30 ++++ tests/unit/moisture_extinction_unit.F90 | 192 ++++++++++++++++++++++++ 2 files changed, 222 insertions(+) create mode 100644 tests/unit/moisture_extinction_unit.F90 diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 6aa31ee..1b3aabd 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -22,6 +22,36 @@ 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:Rothermel 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:Rothermel 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:Rothermel 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:Rothermel fuelmce must be finite") +endforeach() + # 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) diff --git a/tests/unit/moisture_extinction_unit.F90 b/tests/unit/moisture_extinction_unit.F90 new file mode 100644 index 0000000..cc9cdee --- /dev/null +++ b/tests/unit/moisture_extinction_unit.F90 @@ -0,0 +1,192 @@ +program moisture_extinction_unit + + use, intrinsic :: iso_fortran_env, only : REAL64 + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite, ieee_value, ieee_quiet_nan, & + ieee_positive_inf, ieee_negative_inf + use constants_mod, only : CMBCNST, CONVERT_J_PER_KG_TO_BTU_PER_POUND + use fuel_anderson_mod, only : fuel_anderson_t + use ros_wrffire_mod, only : ros_wrffire_t + use fire_physics_mod, only : Calc_flame_length +#ifdef DM_PARALLEL + use mpi +#endif + + implicit none + + type(fuel_anderson_t) :: fuels + type(ros_wrffire_t) :: params + real :: category(0:15,0:3), moisture(0:15,0:3), original(1:14,1:2) + real :: spread(0:15,0:3), flame(0:15,0:3), front(0:15,0:3), area(0:15,0:3) + real :: ir, expected_ir, expected_r0, expected_iboros, factor, nan, wind, slope + integer :: i, j, tile, sample, direction, ierr + character(len=40) :: scenario + +#ifdef DM_PARALLEL + call MPI_Init(ierr) +#endif + call get_command_argument(1, scenario) + call fuels%Initialization(1.0) + call params%Init(0,15,0,3) + nan = ieee_value(0.0, ieee_quiet_nan) + ! Unavailable halos and the nonburnable cell must not enter fuel physics. + category = nan + moisture = nan + do i = 1, 14 + category(i,1:2) = real(i) + end do + area = 0.5 + + do sample = 1, 4 + select case (trim(scenario)) + case ('below') + factor = 0.25 * (sample - 1) + case ('edge') + factor = 1.0 + case ('at') + factor = 1.0 + case ('above') + factor = real(sample + 1) + case default + factor = 0.5 + end select + do i = 1, 13 + moisture(i,1:2) = factor * fuels%fuelmce(i) + if (scenario == 'edge') moisture(i,1:2) = nearest(fuels%fuelmce(i), -1.0) + if (scenario == 'above' .and. sample == 1) moisture(i,1:2) = nearest(fuels%fuelmce(i), 1.0) + end do + select case (trim(scenario)) + case ('negative_moisture') + moisture(1,1) = -0.01 + case ('nan_moisture') + moisture(1,1) = nan + case ('inf_moisture') + moisture(1,1) = ieee_value(0.0, ieee_positive_inf) + case ('neginf_moisture') + moisture(1,1) = ieee_value(0.0, ieee_negative_inf) + case ('zero_extinction') + fuels%fuelmce(1) = 0.0 + case ('negative_extinction') + fuels%fuelmce(1) = -0.1 + case ('nan_extinction') + fuels%fuelmce(1) = nan + case ('inf_extinction') + fuels%fuelmce(1) = ieee_value(0.0, ieee_positive_inf) + case ('neginf_extinction') + fuels%fuelmce(1) = ieee_value(0.0, ieee_negative_inf) + case ('below', 'edge', 'at', 'above') + case default + error stop 'unknown extinction scenario' + end select + original = moisture(1:14,1:2) + ! Exercise disjoint owned tiles through the same routine used on restart. + !$omp parallel do private(tile) + do tile = 1, 2 + call params%Set_params(0,15,0,3,1,14,tile,tile,fuels,category,moisture) + end do + !$omp end parallel do + if (index(scenario, 'moisture') > 0 .or. index(scenario, 'extinction') > 0) & + error stop 'invalid input was accepted' + if (any(moisture(1:13,1:2) /= original(1:13,1:2))) error stop 'physical moisture was modified' + + do i = 1, 13 + call Reference_reaction(i, moisture(i,1), expected_ir, expected_r0, expected_iboros) + do j = 1, 2 + ! Recover reaction intensity from its stored residence-time product. + ir = params%iboros(i,j) / (1055.0 / (0.3048**2 * 60.0) * 1.e-3 * (60.0 * 384.0 / fuels%savr(i))) + call Check_value(ir, expected_ir, 'reaction intensity') + call Check_value(params%iboros(i,j), expected_iboros, 'iboros') + call Check_value(params%r_0(i,j), expected_r0, 'base spread') + if (scenario == 'at' .or. scenario == 'above') then + if (ir /= 0.0 .or. params%iboros(i,j) /= 0.0 .or. params%r_0(i,j) /= 0.0) & + error stop 'nonzero reaction at or above extinction' + end if + end do + ! Hexadecimal values permit an exact before/after comparison below extinction. + if (scenario == 'below') write(*,'(a,2i3,2(1x,z8.8))') 'below signature ', sample, i, & + transfer(params%iboros(i,1),0), transfer(params%r_0(i,1),0) + end do + if (any(params%iboros(14,1:2) /= 0.0) .or. any(params%r_0(14,1:2) /= 0.0)) & + error stop 'nonburnable cell has a reaction' + + do direction = 1, 3 + wind = 0.0 + slope = 0.0 + if (direction > 1) wind = 3.0 + if (direction == 3) slope = 0.2 + do j = 1, 2 + do i = 1, 14 + spread(i,j) = params%Calc_ros(0,15,0,3,i,j,1.0,0.0,wind,0.0,slope,0.0) + end do + end do + call Calc_flame_length(1,14,1,2,0,15,0,3,spread,params%iboros,flame,front,area) + if (any(.not. ieee_is_finite(spread(1:14,1:2))) .or. any(spread(1:14,1:2) < 0.0)) & + error stop 'invalid spread rate' + if (any(.not. ieee_is_finite(flame(1:14,1:2))) .or. any(flame(1:14,1:2) < 0.0)) & + error stop 'invalid flame length' + if (any(front(1:14,1:2) /= spread(1:14,1:2))) error stop 'front spread differs' + do i = 1, 13 + if (moisture(i,1) >= fuels%fuelmce(i)) then + if (any(flame(i,1:2) /= 0.0)) error stop 'nonzero extinction flame length' + ! Anderson 4 uses the independent empirical chaparral spread law. + if (fuels%ichap(i) == 0 .and. any(spread(i,1:2) /= 0.0)) error stop 'Rothermel spread at extinction' + else + if (any(flame(i,1:2) <= 0.0 .and. spread(i,1:2) > 0.0)) error stop 'missing below-extinction flame' + end if + if (fuels%ichap(i) /= 0) call Check_value(spread(i,1), & + min(max(1.2974 * wind**1.41, 0.03333) * wind / (wind + tiny(wind)), 6.0), 'chaparral spread') + end do + if (any(flame(14,1:2) /= 0.0) .or. any(spread(14,1:2) /= 0.0)) error stop 'nonburnable propagation' + end do + end do + print *, 'moisture extinction checks passed: ', trim(scenario) +#ifdef DM_PARALLEL + call MPI_Finalize(ierr) +#endif + +contains + + subroutine Reference_reaction(k, m, intensity, base_ros, integrated) + integer, intent(in) :: k + real, intent(in) :: m + real, intent(out) :: intensity, base_ros, integrated + real(REAL64) :: load, depth, packing, optimum, sigma, a, rate, damping, flux, heat + real :: ratio + ! Evaluate the dimensional Rothermel expressions in double precision. + ! Retain single-precision damping to test the existing polynomial even + ! one representable moisture value below extinction, where it cancels. + ratio = m / fuels%fuelmce(k) + damping = 0.0_REAL64 + if (m < fuels%fuelmce(k)) damping = 1.0 - 2.59 * ratio + 5.11 * ratio**2 - 3.52 * ratio**3 + sigma = fuels%savr(k) + load = real(fuels%fgi(k), REAL64) / (1.0_REAL64 + m) * 0.3048_REAL64**2 * 2.205_REAL64 + depth = real(fuels%fueldepthm(k), REAL64) / 0.3048_REAL64 + packing = load / (depth * fuels%fueldens(k)) + optimum = 3.348_REAL64 * sigma**(-0.8189_REAL64) + a = 1.0_REAL64 / (4.774_REAL64 * sigma**0.1_REAL64 - 7.27_REAL64) + rate = sigma**1.5_REAL64 / (495.0_REAL64 + 0.0594_REAL64 * sigma**1.5_REAL64) * & + (packing / optimum)**a * exp(a * (1.0_REAL64 - packing / optimum)) + heat = CMBCNST * CONVERT_J_PER_KG_TO_BTU_PER_POUND + intensity = rate * load / (1.0_REAL64 + fuels%st(k)) * heat * damping * & + 0.174_REAL64 * real(fuels%se(k), REAL64)**(-0.19_REAL64) + integrated = intensity * 1055.0_REAL64 / (0.3048_REAL64**2 * 60.0_REAL64) * & + 1.e-3_REAL64 * (60.0_REAL64 * 384.0_REAL64 / sigma) + flux = exp((0.792_REAL64 + 0.681_REAL64 * sqrt(sigma)) * (packing + 0.1_REAL64)) / & + (192.0_REAL64 + 0.2595_REAL64 * sigma) + base_ros = intensity * flux / (load / depth * exp(-138.0_REAL64 / sigma) * & + (250.0_REAL64 + 1116.0_REAL64 * m)) + end subroutine Reference_reaction + + subroutine Check_value(value, expected, quantity) + real, intent(in) :: value, expected + character(len=*), intent(in) :: quantity + if (.not. ieee_is_finite(value)) error stop 'nonfinite reaction quantity' + if (value < 0.0) error stop 'negative reaction quantity' + ! Bound accumulated single-precision rounding against the independent + ! double-precision calculation, without a nonzero absolute-error floor. + if (abs(value - expected) > 64.0 * epsilon(value) * abs(expected)) then + print *, quantity, value, expected + error stop 'reaction reference mismatch' + end if + end subroutine Check_value + +end program moisture_extinction_unit From e210e158994db47508fd17994caf79aa4c84d9e7 Mon Sep 17 00:00:00 2001 From: Maria Frediani Date: Wed, 21 Oct 2026 10:55:01 -0600 Subject: [PATCH 2/2] Set Rothermel reaction damping to zero at moisture extinction --- doc/Configuration.rst | 22 ++++++++++++++++++++++ physics/ros_wrffire_mod.F90 | 30 +++++++++++++++++++++++++++--- 2 files changed, 49 insertions(+), 3 deletions(-) diff --git a/doc/Configuration.rst b/doc/Configuration.rst index 8898bdd..d40623b 100644 --- a/doc/Configuration.rst +++ b/doc/Configuration.rst @@ -264,6 +264,28 @@ 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. Nonburnable cells do + not evaluate these fuel quantities. + + 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) diff --git a/physics/ros_wrffire_mod.F90 b/physics/ros_wrffire_mod.F90 index b8cc25f..d110181 100644 --- a/physics/ros_wrffire_mod.F90 +++ b/physics/ros_wrffire_mod.F90 @@ -1,10 +1,13 @@ module ros_wrffire_mod + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite + use constants_mod, only : CMBCNST, CONVERT_J_PER_KG_TO_BTU_PER_POUND use fuel_mod, only : fuel_t use namelist_mod, only : namelist_t use ros_mod, only : ros_t use state_mod, only : state_fire_t + use stderrout_mod, only : Stop_simulation implicit none @@ -140,6 +143,22 @@ subroutine Set_ros_parameters_wrffire (this, ifms, ifme, jfms, jfme, ifts, ifte, this%r_0(i, j) = 0.0 this%iboros(i, j) = 0.0 else + ! Moisture is water mass per dry fuel mass; wetter-than-extinction + ! fuel is valid, but invalid inputs must not be hidden by damping. + if (.not. ieee_is_finite (fmc_g(i, j))) then + write (msg, '(a,2(1x,i0))') 'Rothermel 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,2(1x,i0))') 'Rothermel fmc_g must be nonnegative at i,j:', i, j + call Stop_simulation (msg) + end if + if (.not. ieee_is_finite (fuels%fuelmce(k))) then + write (msg, '(a,1x,i0)') 'Rothermel fuelmce must be finite for fuel category:', k + call Stop_simulation (msg) + else if (fuels%fuelmce(k) <= 0.0) then + write (msg, '(a,1x,i0)') 'Rothermel fuelmce must be positive for fuel category:', k + call Stop_simulation (msg) + end if this%ischap(i, j) = fuels%ichap(k) ! Settings of fire spread parameters from Rothermel ! No need to recalculate if FMC does not change @@ -176,9 +195,14 @@ 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. + ! 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