diff --git a/doc/Configuration.rst b/doc/Configuration.rst index 8898bdd..519d515 100644 --- a/doc/Configuration.rst +++ b/doc/Configuration.rst @@ -215,17 +215,45 @@ output fields retain their existing single-precision NetCDF representation. Number of grid points around lfn=0 where ``fire_viscosity_bg`` is used ``fmoist_run``: *logical* (Default: ``.false.``) - Runs moisture model on the atmospheric grid, outputting the result as a variable named ``fmc_gc`` + Enables evolving class moisture on the fire grid. The first moisture call + initializes the four dead classes from equilibrium and the live class from + ``fuelmc_g_live``, using the first supplied atmospheric forcing. Initialization + occurs independently of the update cadence and copies that forcing into the + previous-forcing state. Missing, nonfinite, or invalid forcing stops the run. + Effective moisture ``fmc_g`` is used by the spread model; class state + ``fmc_gc`` is stored in checkpoints, not ordinary output. + + The effective initialization time is the first call's ``itimestep*dt`` + timestamp, normally ``dt`` seconds from the original run origin. Initialization + applies zero elapsed moisture time. Later updates use the interval since + initialization or the last update. This retains the existing component's + timestamp convention without introducing an artificial initial advance. ``fmoist_freq``: *integer* (Default: ``0``) - Frequency to run moisture model. + Frequency of moisture updates after initialization. 0: use ``fmoist_dt`` k>0: every "k" timesteps ``fmoist_dt``: *real* (Default: ``600.0``) [Units: s] - Time step of moisture model (only used if ``fmoist_freq=0``) + Numerical update interval of the moisture model (only used if ``fmoist_freq=0``), + measured from its effective initialization time. It is distinct from a fuel + class's physical response time. Step-based updates retain their phase relative + to the original run origin. For example, with ``dt=4`` and ``fmoist_freq=3``, + initialization occurs at 4 s, the first update at 12 s uses 8 s elapsed time, + and subsequent updates use 12 s intervals. + + Checkpoints save initialization status/time, class moisture, previous forcing, + last/next update times, and the moisture activation and cadence settings. + A resumed run restores this state without reinitialization and requires the + same activation and cadence. New checkpoints cannot be written before + moisture initialization. PR #39 time-based checkpoints are accepted when their + positive next-update time establishes initialization and their saved interval + matches ``fmoist_dt``. Legacy step-based files are rejected because their + metadata cannot distinguish initialized state from the old delayed-initialization + defect. Prescribed-moisture legacy files remain readable. This compatibility + rule does not resolve the separate humidity-convention audit. ``fire_wind_height``: *real* (Default: ``6.096``) [Units: m] diff --git a/nuopc/fire_behavior_nuopc.F90 b/nuopc/fire_behavior_nuopc.F90 index 91fa869..6e26c23 100644 --- a/nuopc/fire_behavior_nuopc.F90 +++ b/nuopc/fire_behavior_nuopc.F90 @@ -545,6 +545,8 @@ subroutine Realize(model, rc) if (NUOPC_IsConnected(importState, fieldName="mean_prec_rate")) then imp_rainrte = .TRUE. + ! Start rate accumulation once; retain rainfall restored by Advertise. + if (.not. config_flags%restart) grid%fire_rain = 0.0 ! importable field on Grid: mean_prec_rate field = ESMF_FieldCreate(name="mean_prec_rate", grid=fire_grid, & typekind=ESMF_TYPEKIND_R8, rc=rc) diff --git a/physics/fmc_wrffire_mod.F90 b/physics/fmc_wrffire_mod.F90 index a7dab51..b5bd5c9 100644 --- a/physics/fmc_wrffire_mod.F90 +++ b/physics/fmc_wrffire_mod.F90 @@ -13,7 +13,9 @@ module fmc_wrffire_mod ! 50% Entering dormancy, coloration starting, some leaves may have dropped from stem ! 30% Completely cured, treat as dead fuel - use, intrinsic :: iso_fortran_env, only : REAL32 + use, intrinsic :: iso_fortran_env, only : REAL32, INT32 + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite + use netcdf, only : NF90_FILL_FLOAT use stderrout_mod, only : Stop_simulation use netcdf_mod, only : Add_netcdf_att, Add_netcdf_dim, Add_netcdf_var_mpi, NAME_DIM_X, NAME_DIM_Y use ros_mod, only : ros_t @@ -25,6 +27,7 @@ module fmc_wrffire_mod private public :: fmc_wrffire_t, NAME_VAR_FMC_GC, NAME_ATT_FMOIST_LASTTIME, NAME_ATT_FMOIST_NEXTTIME + public :: NAME_ATT_FMC_INITIALIZED, NAME_ATT_FMC_INITTIME integer, parameter :: MOISTURE_CLASSES = 5, NUM_FMEP = 2 real, parameter :: FMEP_DECAY_TLAG = 999999 ! time constant of assimilated adjustments of equilibria decay @@ -32,9 +35,15 @@ module fmc_wrffire_mod character (len = *), parameter :: NAME_VAR_FMC_GC = 'fmc_gc' character (len = *), parameter :: NAME_ATT_FMOIST_LASTTIME = 'fmc_fmoist_lasttime' character (len = *), parameter :: NAME_ATT_FMOIST_NEXTTIME = 'fmc_fmoist_nexttime' + character (len = *), parameter :: NAME_ATT_FMC_INITIALIZED = 'fmc_initialized' + character (len = *), parameter :: NAME_ATT_FMC_INITTIME = 'fmc_initialization_time' type, extends(fmc_t) :: fmc_wrffire_t real :: fmoist_lasttime, fmoist_nexttime, dt_moisture + real :: initialization_time = 0.0 + ! Component-local state: all five moisture classes and previous forcing + ! are initialized by the first advance or restored from a valid restart. + logical :: initialized = .false. real, dimension(MOISTURE_CLASSES) :: rec_drying_lag_sec, rec_wetting_lag_sec, fmc_gc_initial_value logical :: run_advance_moisture real, dimension(:, :, :), allocatable :: fmc_gc ! fuel moisture contents by class @@ -92,6 +101,7 @@ subroutine Add_restart_vars (this, file_restart, cfbm_comm, nx, ny, ifps, ifpe, if (.not. allocated (this%fmc_gc)) call Stop_simulation ('FMC restart write requires allocated fmc_gc') + if (.not. this%initialized) call Stop_simulation ('FMC restart write requires initialized moisture and forcing history') n_moisture_classes = size (this%fmc_gc, 2) if (rank == 0) call Add_netcdf_dim (file_restart, NAME_DIM_MOISTURE_CLASS, n_moisture_classes) @@ -107,6 +117,8 @@ subroutine Add_restart_vars (this, file_restart, cfbm_comm, nx, ny, ifps, ifpe, ifps, ifpe, jfps, jfpe, NAME_VAR_FMC_GC, this%fmc_gc(ifps:ifpe, 1:n_moisture_classes, jfps:jfpe)) if (rank == 0) then + call Add_netcdf_att (file_restart, 'global', NAME_ATT_FMC_INITIALIZED, int (1, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', NAME_ATT_FMC_INITTIME, real (this%initialization_time, kind = REAL32)) call Add_netcdf_att (file_restart, 'global', NAME_ATT_FMOIST_LASTTIME, real (this%fmoist_lasttime, kind = REAL32)) call Add_netcdf_att (file_restart, 'global', NAME_ATT_FMOIST_NEXTTIME, real (this%fmoist_nexttime, kind = REAL32)) end if @@ -141,6 +153,9 @@ subroutine Init_fmc_wrffire (this, fuels, fuelmc_g, fuelmc_g_live, ifms, ifme, j this%fmoist_lasttime = itimestep * dt this%fmoist_nexttime = this%fmoist_lasttime + this%initialized = .false. + this%initialization_time = 0.0 + this%dt_moisture = 0.0 fgi_c = 0.0 fgi_c(1:mfuelcats, 1) = fuels%fgi_1h @@ -202,10 +217,12 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms, integer :: ij real :: time_start, moisture_time + logical :: initialize time_start = itimestep * dt moisture_time = time_start + initialize = .not. this%initialized this%run_advance_moisture = .false. if (fmoist_freq > 0) then @@ -214,8 +231,17 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms, if (.not. moisture_time < this%fmoist_nexttime) this%run_advance_moisture = .true. end if + ! Establish class state at the first call's model timestamp, even when + ! the first scheduled update is later. Initialization has no elapsed time. + if (initialize) this%run_advance_moisture = .true. if (this%run_advance_moisture) then - this%dt_moisture = moisture_time - this%fmoist_lasttime + if (initialize) then + this%dt_moisture = 0.0 + this%initialization_time = moisture_time + else + this%dt_moisture = moisture_time - this%fmoist_lasttime + if (this%dt_moisture <= 0.0) call Stop_simulation ('FMC advance requires increasing moisture time') + end if this%fmoist_lasttime = moisture_time if (fmoist_freq > 0)then continue @@ -228,10 +254,11 @@ subroutine Advance_fmc_model (this, fmoist_freq, fmoist_dt, itimestep, dt, ifms, !$OMP PARALLEL DO & !$OMP PRIVATE (ij) do ij = 1, num_tiles - call this%Advance_moisture_classes (itimestep == 1, ifms, ifme, jfms, jfme, i_start(ij), i_end(ij), j_start(ij), j_end(ij), & + call this%Advance_moisture_classes (initialize, ifms, ifme, jfms, jfme, i_start(ij), i_end(ij), j_start(ij), j_end(ij), & fire_rain, fire_t2, fire_q2, fire_psfc, fire_rain_old, fire_t2_old, fire_q2_old, fire_psfc_old, fire_rh_fire, fuelmc_g) end do !$OMP END PARALLEL DO + this%initialized = .true. !$OMP PARALLEL DO & !$OMP PRIVATE (ij) @@ -301,7 +328,14 @@ subroutine Advance_moisture_classes (this, initialize, ifms, ifme, jfms, jfme, i real :: epsilon, Pws, Pw - if (initialize) call Copy2old () + ! Reject unavailable forcing before copying it into the persistent history + ! or evaluating logarithms and powers in the equilibrium calculation. + call Validate_forcing (rain, t2, q2, psfc) + if (initialize) then + call Copy2old () + else + call Validate_forcing (rain_old, t2_old, q2_old, psfc_old) + end if rhmax = -huge (rhmax) rhmin = huge (rhmin) @@ -442,6 +476,25 @@ subroutine Advance_moisture_classes (this, initialize, ifms, ifme, jfms, jfme, i contains + subroutine Validate_forcing (rain_in, t_in, q_in, p_in) + + real, intent (in) :: rain_in(ifms:ifme, jfms:jfme), t_in(ifms:ifme, jfms:jfme), & + q_in(ifms:ifme, jfms:jfme), p_in(ifms:ifme, jfms:jfme) + integer :: i, j + + do j = jfts, jfte + do i = ifts, ifte + if (.not. all (ieee_is_finite ([rain_in(i,j), t_in(i,j), q_in(i,j), p_in(i,j)]))) & + call Stop_simulation ('FMC requires finite atmospheric forcing') + if (any ([rain_in(i,j), t_in(i,j), q_in(i,j), p_in(i,j)] >= NF90_FILL_FLOAT)) & + call Stop_simulation ('FMC atmospheric forcing is unavailable') + if (t_in(i,j) <= 0.0 .or. p_in(i,j) <= 0.0 .or. q_in(i,j) < 0.0 .or. rain_in(i,j) < 0.0) & + call Stop_simulation ('FMC requires positive temperature and pressure and nonnegative humidity and rain') + end do + end do + + end subroutine Validate_forcing + subroutine Copy2old () implicit none diff --git a/state/state_restart_submod.F90 b/state/state_restart_submod.F90 index 48d41fc..c1bf51d 100644 --- a/state/state_restart_submod.F90 +++ b/state/state_restart_submod.F90 @@ -3,10 +3,13 @@ #ifdef DM_PARALLEL use mpi #endif - use fmc_wrffire_mod, only : fmc_wrffire_t, NAME_VAR_FMC_GC, NAME_ATT_FMOIST_LASTTIME, NAME_ATT_FMOIST_NEXTTIME - use mpi_mod, only : Do_halo_exchange_with_corners - use netcdf_mod, only : Add_netcdf_att, Get_netcdf_att, Get_netcdf_var_mpi, Is_netcdf_file_present + use netcdf, only : nf90_open, nf90_close, nf90_inquire_attribute, nf90_nowrite, nf90_global, nf90_noerr, nf90_enotatt use, intrinsic :: iso_fortran_env, only : INT32, REAL32 + use, intrinsic :: ieee_arithmetic, only : ieee_is_finite + use fmc_wrffire_mod, only : fmc_wrffire_t, NAME_VAR_FMC_GC, NAME_ATT_FMOIST_LASTTIME, NAME_ATT_FMOIST_NEXTTIME, & + NAME_ATT_FMC_INITIALIZED, NAME_ATT_FMC_INITTIME + use mpi_mod, only : Do_halo_exchange_with_corners + use netcdf_mod, only : Add_netcdf_att, Get_netcdf_att, Get_netcdf_var_mpi, Is_netcdf_file_present, Is_netcdf_var_present implicit none @@ -380,6 +383,18 @@ module subroutine Read_restart (this, config_flags) call Validate_restart_integer (file_restart, 'fuel_opt', config_flags%fuel_opt) call Validate_restart_integer (file_restart, 'ros_opt', config_flags%ros_opt) call Validate_restart_integer (file_restart, 'fmc_opt', config_flags%fmc_opt) + if (Has_restart_attribute (file_restart, 'fmoist_run')) then + call Validate_restart_integer (file_restart, 'fmoist_run', merge (1, 0, config_flags%fmoist_run)) + if (config_flags%fmoist_run) then + call Validate_restart_integer (file_restart, 'fmoist_freq', config_flags%fmoist_freq) + call Validate_restart_real (file_restart, 'fmoist_dt', config_flags%fmoist_dt) + end if + else + if (Is_netcdf_var_present (file_restart, NAME_VAR_FMC_GC) .neqv. config_flags%fmoist_run) & + call Stop_simulation ('Legacy restart moisture activation differs from fmoist_run') + if (Has_restart_attribute (file_restart, NAME_ATT_FMC_INITIALIZED)) & + call Stop_simulation ('Restart moisture cadence metadata is incomplete') + end if call Validate_restart_integer (file_restart, 'emis_opt', config_flags%emis_opt) call Validate_restart_integer (file_restart, 'fire_upwinding', config_flags%fire_upwinding) call Validate_restart_integer (file_restart, 'fire_upwinding_reinit', config_flags%fire_upwinding_reinit) @@ -637,8 +652,8 @@ subroutine Read_restart_fmc (file_name) character (len = *), intent (in) :: file_name real (kind = REAL32) :: att_real32 - real :: fmoist_lasttime, fmoist_nexttime - integer :: n_moisture_classes + real :: fmoist_lasttime, fmoist_nexttime, initialization_time + integer :: n_moisture_classes, initialized, i, j if (.not. allocated (this%fmc_param)) call Stop_simulation ('FMC restart read requires allocated fmc_param') @@ -654,16 +669,65 @@ subroutine Read_restart_fmc (file_name) fmoist_lasttime = 0.0 fmoist_nexttime = 0.0 + initialization_time = 0.0 + initialized = 0 if (rank == 0) then call Get_netcdf_att (file_name, 'global', NAME_ATT_FMOIST_LASTTIME, att_real32) fmoist_lasttime = att_real32 call Get_netcdf_att (file_name, 'global', NAME_ATT_FMOIST_NEXTTIME, att_real32) fmoist_nexttime = att_real32 + if (Has_restart_attribute (file_name, NAME_ATT_FMC_INITIALIZED)) then + call Get_netcdf_att (file_name, 'global', NAME_ATT_FMC_INITIALIZED, initialized) + call Get_netcdf_att (file_name, 'global', NAME_ATT_FMC_INITTIME, att_real32) + initialization_time = att_real32 + else + if (Has_restart_attribute (file_name, NAME_ATT_FMC_INITTIME) .or. & + Has_restart_attribute (file_name, 'fmoist_run')) & + call Stop_simulation ('Restart moisture initialization metadata is incomplete') + ! PR #39 advances nexttime only for time-based moisture. A zero + ! nexttime cannot distinguish initialized from defective step runs. + if (config_flags%fmoist_freq /= 0 .or. fmoist_nexttime <= fmoist_lasttime) & + call Stop_simulation ('Legacy restart cannot establish moisture initialization for step-based cadence') + if (abs ((fmoist_nexttime - fmoist_lasttime) - config_flags%fmoist_dt) > & + RESTART_INTERVAL_TOL * max (1.0, abs (config_flags%fmoist_dt))) & + call Stop_simulation ('Legacy restart moisture interval differs from fmoist_dt') + initialization_time = this%dt + initialized = 1 + end if + if (initialized /= 1) call Stop_simulation ('Restart moisture state is not initialized') + if (.not. all (ieee_is_finite ([initialization_time, fmoist_lasttime, fmoist_nexttime]))) & + call Stop_simulation ('Restart moisture times must be finite') + if (initialization_time < this%dt .or. fmoist_lasttime < initialization_time .or. & + fmoist_lasttime > this%itimestep * this%dt .or. fmoist_nexttime < 0.0) & + call Stop_simulation ('Restart moisture times are inconsistent with the fire clock') + if (config_flags%fmoist_freq == 0 .and. fmoist_nexttime <= fmoist_lasttime) & + call Stop_simulation ('Restart next moisture update must follow the last update') end if call Broadcast_restart_real (this, fmoist_lasttime) call Broadcast_restart_real (this, fmoist_nexttime) + call Broadcast_restart_real (this, initialization_time) + call Broadcast_restart_integer (this, initialized) fmc_param%fmoist_lasttime = fmoist_lasttime fmc_param%fmoist_nexttime = fmoist_nexttime + fmc_param%initialization_time = initialization_time + fmc_param%initialized = initialized == 1 + fmc_param%dt_moisture = 0.0 + + do j = this%jfps, this%jfpe + do i = this%ifps, this%ifpe + if (.not. all (ieee_is_finite (fmc_param%fmc_gc(i,:,j))) .or. & + any (fmc_param%fmc_gc(i,:,j) < 0.0) .or. any (fmc_param%fmc_gc(i,:,j) >= NF90_FILL_FLOAT)) & + call Stop_simulation ('Restart class moisture must be finite, available and nonnegative') + if (.not. all (ieee_is_finite ([this%fire_t2_old(i,j), this%fire_psfc_old(i,j), & + this%fire_q2_old(i,j), this%fire_rain_old(i,j)]))) & + call Stop_simulation ('Restart moisture forcing history must be finite') + if (this%fire_t2_old(i,j) <= 0.0 .or. this%fire_psfc_old(i,j) <= 0.0 .or. & + this%fire_q2_old(i,j) < 0.0 .or. this%fire_rain_old(i,j) < 0.0 .or. & + any ([this%fire_t2_old(i,j), this%fire_psfc_old(i,j), this%fire_q2_old(i,j), & + this%fire_rain_old(i,j)] >= NF90_FILL_FLOAT)) & + call Stop_simulation ('Restart moisture forcing history is unavailable or invalid') + end do + end do class default call Stop_simulation ('Restart read is not implemented for selected FMC component') @@ -704,6 +768,19 @@ subroutine Read_restart_field_2d (file_name, var_name, nx, ny, ifms, ifme, jfms, end subroutine Read_restart_field_2d + logical function Has_restart_attribute (file_name, att_name) result (present_attribute) + character (len = *), intent (in) :: file_name, att_name + integer :: ncid, status, close_status + + status = nf90_open (file_name, nf90_nowrite, ncid) + if (status /= nf90_noerr) call Stop_simulation ('Cannot open restart metadata: '//file_name) + status = nf90_inquire_attribute (ncid, nf90_global, att_name) + close_status = nf90_close (ncid) + if (close_status /= nf90_noerr) call Stop_simulation ('Cannot close restart metadata: '//file_name) + if (status /= nf90_noerr .and. status /= nf90_enotatt) call Stop_simulation ('Cannot inspect restart attribute: '//att_name) + present_attribute = status == nf90_noerr + end function Has_restart_attribute + subroutine Validate_restart_integer (file_name, att_name, expected_value) implicit none @@ -815,6 +892,11 @@ module subroutine Write_restart (this, config_flags) call Add_netcdf_att (file_restart, 'global', 'fuel_opt', int (config_flags%fuel_opt, kind = INT32)) call Add_netcdf_att (file_restart, 'global', 'ros_opt', int (config_flags%ros_opt, kind = INT32)) call Add_netcdf_att (file_restart, 'global', 'fmc_opt', int (config_flags%fmc_opt, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'fmoist_run', int (merge (1, 0, config_flags%fmoist_run), kind = INT32)) + if (config_flags%fmoist_run) then + call Add_netcdf_att (file_restart, 'global', 'fmoist_freq', int (config_flags%fmoist_freq, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'fmoist_dt', real (config_flags%fmoist_dt, kind = REAL32)) + end if call Add_netcdf_att (file_restart, 'global', 'emis_opt', int (config_flags%emis_opt, kind = INT32)) call Add_netcdf_att (file_restart, 'global', 'fire_upwinding', int (config_flags%fire_upwinding, kind = INT32)) call Add_netcdf_att (file_restart, 'global', 'fire_upwinding_reinit', int (config_flags%fire_upwinding_reinit, kind = INT32)) diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 6aa31ee..d228ab7 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -17,7 +17,7 @@ endfunction() foreach(test IN ITEMS tiles_unit tile_namelist_unit observed_perimeter_unit namelist_validation_unit time_interval_unit coordinate_precision_unit projection_precision_unit bilinear_precision_unit - fuel_precision_unit wind_order_unit wind_state_unit) + fuel_precision_unit wind_order_unit wind_state_unit fmc_initialization_unit) add_executable(${test} ${test}.F90) target_link_libraries(${test} PRIVATE firelib) endforeach() @@ -128,3 +128,24 @@ add_cfbm_unit_test(NAME wind_state_10m_unit COMMAND wind_state_unit u10m LABELS wind interpolation EXPECT "state wind checks passed") add_cfbm_unit_test(NAME wind_state_invalid_unit COMMAND wind_state_unit invalid LABELS wind interpolation EXPECT "STOP:Error: wrong wind_vinterp_opt") + +foreach(scenario IN ITEMS each delayed time late_first_call) + add_cfbm_unit_test(NAME fmc_${scenario}_unit COMMAND fmc_initialization_unit ${scenario} + LABELS moisture restart EXPECT "moisture initialization checks passed") + if(OPENMP) + set_tests_properties(fmc_${scenario}_unit PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=4;OMP_DYNAMIC=FALSE") + endif() +endforeach() +foreach(scenario IN ITEMS temperature pressure humidity rain) + add_cfbm_unit_test(NAME fmc_${scenario}_unit COMMAND fmc_initialization_unit ${scenario} + LABELS moisture restart EXPECT "STOP:FMC requires positive temperature and pressure and nonnegative humidity and rain") +endforeach() +foreach(scenario IN ITEMS nan old_forcing) + add_cfbm_unit_test(NAME fmc_${scenario}_unit COMMAND fmc_initialization_unit ${scenario} + LABELS moisture restart EXPECT "STOP:FMC requires finite atmospheric forcing") +endforeach() +add_cfbm_unit_test(NAME fmc_missing_unit COMMAND fmc_initialization_unit missing + LABELS moisture restart EXPECT "STOP:FMC atmospheric forcing is unavailable") + +add_cfbm_unit_test(NAME fmc_duplicate_unit COMMAND fmc_initialization_unit duplicate + LABELS moisture restart EXPECT "STOP:FMC advance requires increasing moisture time") diff --git a/tests/unit/fmc_initialization_unit.F90 b/tests/unit/fmc_initialization_unit.F90 new file mode 100644 index 0000000..43403bd --- /dev/null +++ b/tests/unit/fmc_initialization_unit.F90 @@ -0,0 +1,152 @@ +program fmc_initialization_unit + + ! Exercise the moisture lifecycle through the production advance routine. + ! moisture%initialized covers all five moisture classes and their previous + ! atmospheric forcing; allocating storage with Init must leave it false. + ! Check first-call initialization for each-step, delayed-step and time-based + ! cadence, including a first call after step 1. Two owned tiles surrounded by + ! unavailable halos test that forcing validation respects the tile bounds. + ! Verify zero initial elapsed time, unchanged state between updates, and the + ! first drying increment against separately calculated equilibrium values. + ! CTest expects explicit failures for invalid current/previous forcing and a + ! repeated update time. Restart-file restoration is covered by regression cases. + use, intrinsic :: ieee_arithmetic, only : ieee_value, ieee_quiet_nan, ieee_is_finite + use fmc_wrffire_mod, only : fmc_wrffire_t + use fuel_anderson_mod, only : fuel_anderson_t + use ros_wrffire_mod, only : ros_wrffire_t + use netcdf, only : NF90_FILL_FLOAT +#ifdef DM_PARALLEL + use mpi +#endif + implicit none + + type (fmc_wrffire_t) :: moisture + type (fuel_anderson_t) :: fuels + type (ros_wrffire_t) :: ros + real :: rain(0:5,0:3), t(0:5,0:3), q(0:5,0:3), p(0:5,0:3) + real :: rain_old(0:5,0:3), t_old(0:5,0:3), q_old(0:5,0:3), p_old(0:5,0:3) + real :: rh(0:5,0:3), fmc(0:5,0:3), category(0:5,0:3) + real :: initial(4,5,2), initial_fmc(4,2), interval, elapsed, expected, d, w, change + integer :: first_step, next_step, frequency, step, ierr + integer :: i_start(2) = [1,1], i_end(2) = [4,4], j_start(2) = [1,2], j_end(2) = [1,2] + character (len = 32) :: scenario + +#ifdef DM_PARALLEL + call MPI_Init (ierr) +#endif + call get_command_argument (1, scenario) + call fuels%Initialization (0.8) + call moisture%Init (fuels, 0.08, 0.30, 0, 5, 0, 3, 0, 4.0) + call ros%Init (0, 5, 0, 3) + if (moisture%initialized) error stop 'allocation initialized moisture before forcing' + moisture%fmc_gc = NF90_FILL_FLOAT + t = NF90_FILL_FLOAT + q = NF90_FILL_FLOAT + p = NF90_FILL_FLOAT + rain = NF90_FILL_FLOAT + t(1:4,1:2) = 300.0 + q(1:4,1:2) = 0.008 + p(1:4,1:2) = 90000.0 + rain(1:4,1:2) = 7.0 + t_old = NF90_FILL_FLOAT + q_old = NF90_FILL_FLOAT + p_old = NF90_FILL_FLOAT + rain_old = NF90_FILL_FLOAT + category = 1.0 + frequency = 3 + interval = 16.0 + first_step = 1 + next_step = 3 + select case (scenario) + case ('each') + frequency = 1 + next_step = 2 + case ('time') + frequency = 0 + next_step = 5 + case ('late_first_call') + first_step = 7 + next_step = 9 + case ('missing') + t(4,2) = NF90_FILL_FLOAT + case ('nan') + q(4,2) = ieee_value (0.0, ieee_quiet_nan) + case ('temperature') + t(4,2) = 0.0 + case ('pressure') + p(4,2) = 0.0 + case ('humidity') + q(4,2) = -0.001 + case ('rain') + rain(4,2) = -1.0 + end select + + call Advance (first_step) + if (.not. moisture%initialized) error stop 'delayed cadence skipped initialization' + if (moisture%initialization_time /= first_step * 4.0) error stop 'wrong initialization time' + if (moisture%fmoist_lasttime /= first_step * 4.0) error stop 'wrong first moisture time' + if (moisture%dt_moisture /= 0.0) error stop 'artificial elapsed interval at initialization' + if (.not. all (ieee_is_finite (moisture%fmc_gc(1:4,:,1:2)))) error stop 'invalid initial class state' + if (any (moisture%fmc_gc(1:4,5,1:2) /= 0.30)) error stop 'initial live moisture advanced' + call Equilibria (300.0, 0.008, d, w) + expected = 0.5 * (d + w) + if (any (abs (moisture%fmc_gc(1:4,1:4,1:2) - expected) > 2.0e-7)) error stop 'wrong initial equilibrium' + if (any (t_old(1:4,1:2) /= t(1:4,1:2)) .or. any (q_old(1:4,1:2) /= q(1:4,1:2)) .or. & + any (p_old(1:4,1:2) /= p(1:4,1:2)) .or. any (rain_old(1:4,1:2) /= rain(1:4,1:2))) & + error stop 'initial forcing history was not copied' + initial = moisture%fmc_gc(1:4,:,1:2) + initial_fmc = fmc(1:4,1:2) + if (scenario == 'old_forcing') q_old(4,2) = ieee_value (0.0, ieee_quiet_nan) + + ! Dry the forcing after initialization. Classes and previous forcing must + ! remain at their initialized values until the next scheduled moisture update. + t(1:4,1:2) = 310.0 + q(1:4,1:2) = 0.002 + do step = first_step + 1, next_step - 1 + call Advance (step) + if (any (moisture%fmc_gc(1:4,:,1:2) /= initial)) error stop 'off-cadence class advance' + if (any (fmc(1:4,1:2) /= initial_fmc)) error stop 'off-cadence effective moisture advance' + if (any (t_old(1:4,1:2) /= 300.0)) error stop 'off-cadence forcing history changed' + end do + call Advance (next_step) + elapsed = (next_step - first_step) * 4.0 + if (moisture%dt_moisture /= elapsed) error stop 'wrong elapsed interval after initialization' + call Equilibria (305.0, 0.005, d, w) + change = elapsed / 3600.0 + expected = initial(1,1,1) + (max (min (initial(1,1,1), d), w) - initial(1,1,1)) * change * (1.0 - 0.5 * change) + if (any (abs (moisture%fmc_gc(1:4,1,1:2) - expected) > 2.0e-7)) error stop 'first drying increment differs' + if (any (moisture%fmc_gc(1:4,1,1:2) >= initial(:,1,:))) error stop 'drying did not advance moisture' + if (moisture%initialization_time /= first_step * 4.0) error stop 'initialization was repeated' + if (frequency == 0 .and. moisture%fmoist_nexttime /= next_step * 4.0 + interval) error stop 'time cadence changed' + if (scenario == 'duplicate') call Advance (next_step) + print *, 'moisture initialization checks passed: ', trim (scenario) +#ifdef DM_PARALLEL + call MPI_Finalize (ierr) +#endif + +contains + + subroutine Advance (n) + integer, intent (in) :: n + call moisture%Advance_fmc_model (frequency, interval, n, 4.0, 0, 5, 0, 3, & + i_start, i_end, j_start, j_end, 2, rain, t, q, p, rain_old, t_old, q_old, p_old, rh, & + 0.08, fmc, category, fuels, ros) + end subroutine Advance + + subroutine Equilibria (temperature, humidity, drying, wetting) + real, intent (in) :: temperature, humidity + real, intent (out) :: drying, wetting + real :: vapor_pressure, saturation_pressure, percent_rh + ! Retain the existing thermodynamic expression here. P01 owns its correction. + vapor_pressure = humidity * 90000.0 / (0.622 + 0.378 * humidity) + saturation_pressure = exp (54.842763 - 6763.22 / temperature - 4.210 * log (temperature) + & + 0.000367 * temperature + tanh (0.0415 * (temperature - 218.8)) * & + (53.878 - 1331.22 / temperature - 9.44523 * log (temperature) + 0.014025 * temperature)) + percent_rh = min (vapor_pressure / saturation_pressure, 1.0) * 100.0 + drying = (0.942 * percent_rh ** 0.679 + 0.000499 * exp (0.1 * percent_rh) + & + 0.18 * (294.25 - temperature) * (1.0 - exp (-0.115 * percent_rh))) * 0.01 + wetting = (0.618 * percent_rh ** 0.753 + 0.000454 * exp (0.1 * percent_rh) + & + 0.18 * (294.25 - temperature) * (1.0 - exp (-0.115 * percent_rh))) * 0.01 + end subroutine Equilibria + +end program fmc_initialization_unit