diff --git a/CMakeLists.txt b/CMakeLists.txt index e67f82d..8d9e5fb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -67,6 +67,7 @@ endif() # State list(APPEND _state_files state/state_mod.F90 + state/state_restart_submod.F90 state/ignition_line_mod.F90 state/tiles_mod.F90) diff --git a/driver/fire_behavior.F90 b/driver/fire_behavior.F90 index 9b53fb2..f519af3 100644 --- a/driver/fire_behavior.F90 +++ b/driver/fire_behavior.F90 @@ -8,6 +8,7 @@ program fire_behavior use initialize_mod, only : Init_fire_state, Init_atm_state use advance_mod, only : Advance_state use wrfdata_mod, only : wrfdata_t + use stderrout_mod, only : Print_message use, intrinsic :: iso_fortran_env, only : ERROR_UNIT, OUTPUT_UNIT implicit none @@ -71,16 +72,26 @@ program fire_behavior end select - if (DEBUG_LOCAL) write (OUTPUT_UNIT, *) ' Saving fire state...' - call grid%Save_state () + if (.not. config_flags%restart) then + if (DEBUG_LOCAL) call Print_message (' Saving fire state') + call grid%Save_state () + end if + + call grid%Handle_restart (config_flags, initialize = .true.) if (DEBUG_LOCAL) write (OUTPUT_UNIT, *) ' Starting temporal loop...' do while (grid%datetime_now < grid%datetime_end) call Advance_state (grid, config_flags) + + call grid%Check_model_clock () + ! Save the forcing used for the completed interval, as in the NUOPC cap. ! Refresh atmospheric fields afterward for the next fire advance. + call grid%Handle_output (config_flags) if (config_flags%ideal_opt == 0) call grid%Handle_wrfdata_update (atm_state, config_flags) + + call grid%Handle_restart (config_flags) end do if (DEBUG_LOCAL) write (OUTPUT_UNIT, *) ' Completed temporal loop' diff --git a/driver/initialize_mod.F90 b/driver/initialize_mod.F90 index 1f8be08..6ec51d0 100644 --- a/driver/initialize_mod.F90 +++ b/driver/initialize_mod.F90 @@ -63,46 +63,52 @@ subroutine Init_fire_state (grid, config_flags, wrf) if (DEBUG_LOCAL) call Print_message (' Initialization...') if (config_flags%ideal_opt == 0) then ! Real world - if (DEBUG_LOCAL) call Print_message (' Reading geogrid file') + if (.not. config_flags%restart) then + if (DEBUG_LOCAL) call Print_message (' Reading geogrid file') #ifdef DM_PARALLEL - if (grid%is_cfbm_comm_set) then - call Convert_mpi_comm_to_f08 (grid%cfbm_comm, cfbm_comm_f08) - call Mpi_comm_rank (cfbm_comm_f08, rank, ierr) - else - call Stop_simulation ('The MPI communicator cfbm_comm has not been set') - end if + if (grid%is_cfbm_comm_set) then + call Convert_mpi_comm_to_f08 (grid%cfbm_comm, cfbm_comm_f08) + call Mpi_comm_rank (cfbm_comm_f08, rank, ierr) + else + call Stop_simulation ('The MPI communicator cfbm_comm has not been set') + end if #else - rank = 0 + rank = 0 #endif - if (rank == 0) geogrid = geogrid_t (file_name = 'geo_em.d01.nc') + if (rank == 0) geogrid = geogrid_t (file_name = 'geo_em.d01.nc') #ifdef DM_PARALLEL - call MPI_Bcast (geogrid%cen_lon, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) - call MPI_Bcast (geogrid%cen_lat, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) - call MPI_Bcast (geogrid%dx, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) - call MPI_Bcast (geogrid%dy, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) - call MPI_Bcast (geogrid%true_lat_1, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) - call MPI_Bcast (geogrid%true_lat_2, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) - call MPI_Bcast (geogrid%stand_lon, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) - - call MPI_Bcast(geogrid%map_proj, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - call MPI_Bcast(geogrid%sr_x, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - call MPI_Bcast(geogrid%sr_y, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - - call MPI_Bcast(geogrid%ifds, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - call MPI_Bcast(geogrid%ifde, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - call MPI_Bcast(geogrid%jfds, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - call MPI_Bcast(geogrid%jfde, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - - call MPI_Bcast(geogrid%ids, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - call MPI_Bcast(geogrid%ide, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - call MPI_Bcast(geogrid%jds, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) - call MPI_Bcast(geogrid%jde, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast (geogrid%cen_lon, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) + call MPI_Bcast (geogrid%cen_lat, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) + call MPI_Bcast (geogrid%dx, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) + call MPI_Bcast (geogrid%dy, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) + call MPI_Bcast (geogrid%true_lat_1, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) + call MPI_Bcast (geogrid%true_lat_2, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) + call MPI_Bcast (geogrid%stand_lon, 1, MPI_REAL, 0, cfbm_comm_f08, ierr) + + call MPI_Bcast(geogrid%map_proj, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast(geogrid%sr_x, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast(geogrid%sr_y, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + + call MPI_Bcast(geogrid%ifds, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast(geogrid%ifde, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast(geogrid%jfds, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast(geogrid%jfde, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + + call MPI_Bcast(geogrid%ids, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast(geogrid%ide, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast(geogrid%jds, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Bcast(geogrid%jde, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) #endif + end if if (DEBUG_LOCAL) call Print_message (' Initializing fire state') - call grid%Initialization (config_flags, geogrid) + if (config_flags%restart) then + call grid%Initialization (config_flags) + else + call grid%Initialization (config_flags, geogrid) + end if else ! Ideal @@ -112,11 +118,16 @@ subroutine Init_fire_state (grid, config_flags, wrf) call Init_fire_components (grid, config_flags) - if (present (wrf)) then + if (present (wrf) .and. .not. config_flags%restart) then if (DEBUG_LOCAL) call Print_message (' Initializing atmospheric state') call grid%Handle_wrfdata_update (wrf, config_flags) end if + if (config_flags%restart) then + if (DEBUG_LOCAL) call Print_message (' Reading restart state') + call grid%Read_restart (config_flags) + end if + if (DEBUG_LOCAL) then ! print lat/lons open (newunit = unit_out, file = 'latlons_c.dat') @@ -187,6 +198,8 @@ subroutine Init_fire_state_within_wrf (state, config_flags, & call Init_fire_components (state, config_flags) + if (config_flags%restart) call state%Read_restart (config_flags) + if (DEBUG_LOCAL) call Print_message (' Leaving subroutine Init_fire_state_within_wrf...') end subroutine Init_fire_state_within_wrf diff --git a/io/namelist_mod.F90 b/io/namelist_mod.F90 index 033512e..1befbc0 100644 --- a/io/namelist_mod.F90 +++ b/io/namelist_mod.F90 @@ -38,6 +38,10 @@ module namelist_mod integer :: interval_output = -1 ! Frequency to save the output [s] real :: dt = 2.0 ! Time step of the fire model [s] + ! Restart controls + logical :: restart = .false. ! Read initial dynamic fire state from fire_restart_.nc + integer :: restart_interval = -1 ! Restart output interval [s]; -1 disables restart writes + integer :: num_tiles = 1 ! Number of tiles in each patch integer :: tile_strategy = 0 ! Strategy for the tile decomposition: 0) ... @@ -216,6 +220,8 @@ subroutine Broadcast_nml (this, mpi_comm_cfbm) call Broadcast_integer (this%ideal_opt) call Broadcast_integer (this%devel_opt) + call Broadcast_logical (this%restart) + call Broadcast_integer (this%restart_interval) call Broadcast_integer (this%fuel_opt) call Broadcast_integer (this%ros_opt) call Broadcast_integer (this%fmc_opt) @@ -410,6 +416,9 @@ subroutine Check_nml (this, require_atm_interval) call Stop_simulation ('perimeter ignition time must align with a fire timestep') end if + if (this%restart_interval /= -1 .and. this%restart_interval <= 0) & + call Stop_simulation ('restart_interval must be positive or -1') + end subroutine Check_nml subroutine Check_time_intervals (this, require_atm_interval) @@ -829,15 +838,16 @@ subroutine Init_time_block (this, file_name) integer :: start_year, start_month, start_day, start_hour, start_minute, start_second, & end_year, end_month, end_day, end_hour, end_minute, end_second, interval_output, & - num_tiles, tile_strategy + restart_interval, num_tiles, tile_strategy real :: dt + logical :: restart character (len = :), allocatable :: msg integer :: unit_nml, io_stat namelist /time/ start_year, start_month, start_day, start_hour, start_minute, start_second, & end_year, end_month, end_day, end_hour, end_minute, end_second, dt, interval_output, & - num_tiles, tile_strategy + restart, restart_interval, num_tiles, tile_strategy ! Set default values @@ -855,6 +865,8 @@ subroutine Init_time_block (this, file_name) end_second = this%end_second dt = this%dt interval_output = this%interval_output + restart = this%restart + restart_interval = this%restart_interval num_tiles = this%num_tiles tile_strategy = this%tile_strategy @@ -886,6 +898,8 @@ subroutine Init_time_block (this, file_name) this%end_second = end_second this%dt = dt this%interval_output = interval_output + this%restart = restart + this%restart_interval = restart_interval this%num_tiles = num_tiles this%tile_strategy = tile_strategy diff --git a/io/netcdf_mod.F90 b/io/netcdf_mod.F90 index b7a22b7..eeea7d1 100644 --- a/io/netcdf_mod.F90 +++ b/io/netcdf_mod.F90 @@ -1,8 +1,9 @@ module netcdf_mod + use, intrinsic :: iso_fortran_env, only : REAL64 use netcdf, only : NF90_FILL_FLOAT use stderrout_mod, only : Stop_simulation, Print_message - use mpi_mod, only : Gather_var2d + use mpi_mod, only : Gather_var2d, Distribute_global_var2d implicit none @@ -10,26 +11,42 @@ module netcdf_mod character (len = 2), parameter :: NAME_DIM_X = 'nx', NAME_DIM_Y = 'ny' public :: Get_netcdf_var, Get_netcdf_att, Get_netcdf_dim, Create_netcdf_file, Add_netcdf_dim, Add_netcdf_var, & - Is_netcdf_file_present, Is_netcdf_var_present, Add_netcdf_var_mpi, NAME_DIM_X, NAME_DIM_Y, NF90_FILL_FLOAT + Add_netcdf_att, Is_netcdf_file_present, Is_netcdf_var_present, Add_netcdf_var_mpi, Get_netcdf_var_mpi, & + NAME_DIM_X, NAME_DIM_Y, NF90_FILL_FLOAT + + interface Add_netcdf_att + module procedure Add_netcdf_att_int32 + module procedure Add_netcdf_att_real32 + end interface Add_netcdf_att interface Add_netcdf_var module procedure Add_netcdf_var_real32_2d + module procedure Add_netcdf_var_real64_2d module procedure Add_netcdf_var_real32_3d end interface Add_netcdf_var interface Add_netcdf_var_mpi module procedure Add_netcdf_var_real32_2d_mpi + module procedure Add_netcdf_var_real64_2d_mpi + module procedure Add_netcdf_var_real32_3d_mpi end interface Add_netcdf_var_mpi interface Get_netcdf_var module procedure Get_netcdf_var_char_1d module procedure Get_netcdf_var_int32_3d module procedure Get_netcdf_var_real32_2d + module procedure Get_netcdf_var_real64_2d module procedure Get_netcdf_var_real32_3d module procedure Get_netcdf_var_real64_3d module procedure Get_netcdf_var_real32_4d end interface Get_netcdf_var + interface Get_netcdf_var_mpi + module procedure Get_netcdf_var_real32_2d_mpi + module procedure Get_netcdf_var_real64_2d_mpi + module procedure Get_netcdf_var_real32_3d_mpi + end interface Get_netcdf_var_mpi + interface Get_netcdf_att module procedure Get_netcdf_att_int32 module procedure Get_netcdf_att_real32 @@ -37,6 +54,84 @@ module netcdf_mod contains + subroutine Add_netcdf_att_int32 (file_name, var_name, att_name, att_value) + + use netcdf + + use, intrinsic :: iso_fortran_env, only : INT32 + + implicit none + + character (len = *), intent(in) :: file_name, var_name, att_name + integer (kind = INT32), intent(in) :: att_value + + integer :: status, ncid, varid + + + status = nf90_open (trim(file_name), NF90_WRITE, ncid) + call Check_status (status) + + status = nf90_redef (ncid) + call Check_status (status) + + if (var_name == 'global') then + status = nf90_put_att (ncid, NF90_GLOBAL, att_name, att_value) + call Check_status (status) + else + status = nf90_inq_varid (ncid, var_name, varid) + call Check_status (status) + + status = nf90_put_att (ncid, varid, att_name, att_value) + call Check_status (status) + end if + + status = nf90_enddef (ncid) + call Check_status (status) + + status = nf90_close (ncid) + call Check_status (status) + + end subroutine Add_netcdf_att_int32 + + subroutine Add_netcdf_att_real32 (file_name, var_name, att_name, att_value) + + use netcdf + + use, intrinsic :: iso_fortran_env, only : REAL32 + + implicit none + + character (len = *), intent(in) :: file_name, var_name, att_name + real (kind = REAL32), intent(in) :: att_value + + integer :: status, ncid, varid + + + status = nf90_open (trim(file_name), NF90_WRITE, ncid) + call Check_status (status) + + status = nf90_redef (ncid) + call Check_status (status) + + if (var_name == 'global') then + status = nf90_put_att (ncid, NF90_GLOBAL, att_name, att_value) + call Check_status (status) + else + status = nf90_inq_varid (ncid, var_name, varid) + call Check_status (status) + + status = nf90_put_att (ncid, varid, att_name, att_value) + call Check_status (status) + end if + + status = nf90_enddef (ncid) + call Check_status (status) + + status = nf90_close (ncid) + call Check_status (status) + + end subroutine Add_netcdf_att_real32 + subroutine Add_netcdf_dim (file_name, name_dim, val_dim) use netcdf @@ -123,6 +218,63 @@ subroutine Add_netcdf_var_real32_2d (file_name, name_dims, varname, var, fill_va end subroutine Add_netcdf_var_real32_2d + subroutine Add_netcdf_var_real64_2d (file_name, name_dims, varname, var, fill_value, units, long_name) + + use netcdf + + implicit none + + character (len = *), intent (in) :: file_name, varname + character (len = *), dimension(:), intent (in) :: name_dims + real (kind = REAL64), dimension(:, :), intent (in) :: var + real (kind = REAL64), intent (in), optional :: fill_value + character (len = *), intent (in), optional :: units, long_name + + integer, parameter :: N_DIMS = 2 + integer :: ncidout, status, varid, n + integer :: dimids(N_DIMS) + + + status = nf90_open (file_name, NF90_WRITE, ncidout) + call Check_status (status) + + status = nf90_redef (ncidout) + call Check_status (status) + + do n = 1, N_DIMS + status = nf90_inq_dimid(ncidout, name_dims(n), dimids(n)) + call Check_status (status) + end do + + status = nf90_inq_varid (ncidout, varname, varid) + if (status == NF90_ENOTVAR) then + status = nf90_def_var (ncidout, varname, NF90_DOUBLE, dimids, varid) + call Check_status (status) + if (present (fill_value)) then + status = nf90_put_att (ncidout, varid, '_FillValue', fill_value) + call Check_status (status) + end if + if (present (units)) then + status = nf90_put_att (ncidout, varid, 'units', units) + call Check_status (status) + end if + if (present (long_name)) then + status = nf90_put_att (ncidout, varid, 'long_name', long_name) + call Check_status (status) + end if + end if + + status = nf90_enddef (ncidout) + call Check_status (status) + + status = nf90_put_var (ncidout, varid, var) + call Check_status (status) + + status = nf90_close (ncidout) + call Check_status (status) + + end subroutine Add_netcdf_var_real64_2d + subroutine Add_netcdf_var_real32_2d_mpi (file_name, cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var_name, var2d_local, & fill_value, units, long_name) @@ -168,6 +320,145 @@ subroutine Add_netcdf_var_real32_2d_mpi (file_name, cfbm_comm, nx, ny, ifps, ifp end subroutine Add_netcdf_var_real32_2d_mpi + subroutine Add_netcdf_var_real64_2d_mpi (file_name, cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var_name, var2d_local, & + fill_value, units, long_name) + +#ifdef DM_PARALLEL + use mpi +#endif + implicit none + + integer, intent (in) :: cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe + character (len = *), intent (in) :: file_name, var_name + real (kind = REAL64), dimension(ifps:ifpe, jfps:jfpe), intent (in) :: var2d_local + real (kind = REAL64), intent (in), optional :: fill_value + character (len = *), intent (in), optional :: units, long_name + + real (kind = REAL64), dimension(nx, ny) :: var2d + integer :: rank, ierr + logical, parameter :: DEBUG_LOCAL = .false. + + + if (DEBUG_LOCAL) call Print_message ('Entering Add_netcdf_var_real64_2d_mpi...') + + if (DEBUG_LOCAL) call Print_message ('file name = ' // trim (file_name)) + if (DEBUG_LOCAL) call Print_message ('var name = ' // trim (var_name)) + if (DEBUG_LOCAL) call Print_message ('dim X name = ' // trim (NAME_DIM_X)) + if (DEBUG_LOCAL) call Print_message ('dim Y name = ' // trim (NAME_DIM_Y)) + +#ifdef DM_PARALLEL + call Mpi_comm_rank (cfbm_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') + + call Gather_var2d (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_local(ifps:ifpe, jfps:jfpe), var2d) + + if (rank == 0) then + call Add_netcdf_var (file_name, [NAME_DIM_X, NAME_DIM_Y], var_name, var2d(1:nx, 1:ny), & + fill_value, units, long_name) + end if +#else + call Add_netcdf_var (file_name, [NAME_DIM_X, NAME_DIM_Y], var_name, var2d_local(1:nx, 1:ny), & + fill_value, units, long_name) +#endif + + if (DEBUG_LOCAL) call Print_message ('Leaving Add_netcdf_var_real64_2d_mpi...') + + end subroutine Add_netcdf_var_real64_2d_mpi + + subroutine Get_netcdf_var_real32_2d_mpi (file_name, cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var_name, var2d_local) + +#ifdef DM_PARALLEL + use mpi +#endif + implicit none + + character (len = *), intent (in) :: file_name, var_name + integer, intent (in) :: cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe + real, dimension(ifps:ifpe, jfps:jfpe), intent (out) :: var2d_local + + real, dimension(:, :), allocatable :: var2d_global + character (len = :), allocatable :: msg + integer :: rank, ierr + logical, parameter :: DEBUG_LOCAL = .false. + + + if (DEBUG_LOCAL) call Print_message ('Entering Get_netcdf_var_real32_2d_mpi...') + +#ifdef DM_PARALLEL + call Mpi_comm_rank (cfbm_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') + + if (rank == 0) then + call Get_netcdf_var (file_name, var_name, var2d_global) + if (size (var2d_global, 1) /= nx .or. size (var2d_global, 2) /= ny) then + msg = 'NetCDF variable has unexpected dimensions: '//trim (var_name) + call Stop_simulation (msg) + end if + end if + + call Distribute_global_var2d (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_global, var2d_local) +#else + call Get_netcdf_var (file_name, var_name, var2d_global) + if (size (var2d_global, 1) /= nx .or. size (var2d_global, 2) /= ny) then + msg = 'NetCDF variable has unexpected dimensions: '//trim (var_name) + call Stop_simulation (msg) + end if + var2d_local(1:nx, 1:ny) = var2d_global(1:nx, 1:ny) +#endif + + if (allocated (var2d_global)) deallocate (var2d_global) + + if (DEBUG_LOCAL) call Print_message ('Leaving Get_netcdf_var_real32_2d_mpi...') + + end subroutine Get_netcdf_var_real32_2d_mpi + + subroutine Get_netcdf_var_real64_2d_mpi (file_name, cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var_name, var2d_local) + +#ifdef DM_PARALLEL + use mpi +#endif + implicit none + + character (len = *), intent (in) :: file_name, var_name + integer, intent (in) :: cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe + real (kind = REAL64), dimension(ifps:ifpe, jfps:jfpe), intent (out) :: var2d_local + + real (kind = REAL64), dimension(:, :), allocatable :: var2d_global + character (len = :), allocatable :: msg + integer :: rank, ierr + logical, parameter :: DEBUG_LOCAL = .false. + + + if (DEBUG_LOCAL) call Print_message ('Entering Get_netcdf_var_real64_2d_mpi...') + +#ifdef DM_PARALLEL + call Mpi_comm_rank (cfbm_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') + + if (rank == 0) then + call Get_netcdf_var (file_name, var_name, var2d_global) + if (size (var2d_global, 1) /= nx .or. size (var2d_global, 2) /= ny) then + msg = 'NetCDF variable has unexpected dimensions: '//trim (var_name) + call Stop_simulation (msg) + end if + end if + + call Distribute_global_var2d (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_global, var2d_local) +#else + call Get_netcdf_var (file_name, var_name, var2d_global) + if (size (var2d_global, 1) /= nx .or. size (var2d_global, 2) /= ny) then + msg = 'NetCDF variable has unexpected dimensions: '//trim (var_name) + call Stop_simulation (msg) + end if + var2d_local(1:nx, 1:ny) = var2d_global(1:nx, 1:ny) +#endif + + if (allocated (var2d_global)) deallocate (var2d_global) + + if (DEBUG_LOCAL) call Print_message ('Leaving Get_netcdf_var_real64_2d_mpi...') + + end subroutine Get_netcdf_var_real64_2d_mpi + subroutine Add_netcdf_var_real32_3d (file_name, name_dims, varname, var) use netcdf @@ -211,6 +502,104 @@ subroutine Add_netcdf_var_real32_3d (file_name, name_dims, varname, var) end subroutine Add_netcdf_var_real32_3d + subroutine Add_netcdf_var_real32_3d_mpi (file_name, name_dims, cfbm_comm, nx, ny, n2, ifps, ifpe, jfps, jfpe, varname, var3d_local) + +#ifdef DM_PARALLEL + use mpi +#endif + implicit none + + character (len = *), intent (in) :: file_name, varname + character (len = *), dimension(:), intent (in) :: name_dims + integer, intent (in) :: cfbm_comm, nx, ny, n2, ifps, ifpe, jfps, jfpe + real, dimension(ifps:ifpe, n2, jfps:jfpe), intent (in) :: var3d_local + + real, dimension(ifps:ifpe, jfps:jfpe) :: var2d_local + real, dimension(nx, ny) :: var2d_global + real, dimension(:, :, :), allocatable :: var3d_global + integer :: rank, ierr, k + logical, parameter :: DEBUG_LOCAL = .false. + + + if (DEBUG_LOCAL) call Print_message ('Entering Add_netcdf_var_real32_3d_mpi...') + +#ifdef DM_PARALLEL + call Mpi_comm_rank (cfbm_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') + + if (rank == 0) allocate (var3d_global(nx, n2, ny)) + do k = 1, n2 + var2d_local(ifps:ifpe, jfps:jfpe) = var3d_local(ifps:ifpe, k, jfps:jfpe) + call Gather_var2d (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_local, var2d_global) + if (rank == 0) var3d_global(1:nx, k, 1:ny) = var2d_global(1:nx, 1:ny) + end do + + if (rank == 0) call Add_netcdf_var (file_name, name_dims, varname, var3d_global) +#else + call Add_netcdf_var (file_name, name_dims, varname, var3d_local(1:nx, 1:n2, 1:ny)) +#endif + + if (allocated (var3d_global)) deallocate (var3d_global) + + if (DEBUG_LOCAL) call Print_message ('Leaving Add_netcdf_var_real32_3d_mpi...') + + end subroutine Add_netcdf_var_real32_3d_mpi + + subroutine Get_netcdf_var_real32_3d_mpi (file_name, cfbm_comm, nx, ny, n2, ifps, ifpe, jfps, jfpe, var_name, var3d_local) + +#ifdef DM_PARALLEL + use mpi +#endif + implicit none + + character (len = *), intent (in) :: file_name, var_name + integer, intent (in) :: cfbm_comm, nx, ny, n2, ifps, ifpe, jfps, jfpe + real, dimension(ifps:ifpe, n2, jfps:jfpe), intent (out) :: var3d_local + + real, dimension(:, :), allocatable :: var2d_global + real, dimension(:, :, :), allocatable :: var3d_global + character (len = :), allocatable :: msg + integer :: rank, ierr, k + logical, parameter :: DEBUG_LOCAL = .false. + + + if (DEBUG_LOCAL) call Print_message ('Entering Get_netcdf_var_real32_3d_mpi...') + +#ifdef DM_PARALLEL + call Mpi_comm_rank (cfbm_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') + + if (rank == 0) then + call Get_netcdf_var (file_name, var_name, var3d_global) + if (size (var3d_global, 1) /= nx .or. size (var3d_global, 2) /= n2 .or. size (var3d_global, 3) /= ny) then + msg = 'NetCDF variable has unexpected dimensions: '//trim (var_name) + call Stop_simulation (msg) + end if + end if + + do k = 1, n2 + if (rank == 0) then + allocate (var2d_global(nx, ny)) + var2d_global(1:nx, 1:ny) = var3d_global(1:nx, k, 1:ny) + end if + call Distribute_global_var2d (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_global, var3d_local(ifps:ifpe, k, jfps:jfpe)) + if (allocated (var2d_global)) deallocate (var2d_global) + end do +#else + call Get_netcdf_var (file_name, var_name, var3d_global) + if (size (var3d_global, 1) /= nx .or. size (var3d_global, 2) /= n2 .or. size (var3d_global, 3) /= ny) then + msg = 'NetCDF variable has unexpected dimensions: '//trim (var_name) + call Stop_simulation (msg) + end if + var3d_local(1:nx, 1:n2, 1:ny) = var3d_global(1:nx, 1:n2, 1:ny) +#endif + + if (allocated (var3d_global)) deallocate (var3d_global) + + if (DEBUG_LOCAL) call Print_message ('Leaving Get_netcdf_var_real32_3d_mpi...') + + end subroutine Get_netcdf_var_real32_3d_mpi + subroutine Check_status (status) @@ -510,6 +899,61 @@ subroutine Get_netcdf_var_real32_2d (file_name, var_name, output) end subroutine Get_netcdf_var_real32_2d + subroutine Get_netcdf_var_real64_2d (file_name, var_name, output) + + use netcdf + use, intrinsic :: iso_fortran_env, only : REAL64 + + implicit none + + character (len = *), intent(in) :: file_name, var_name + real (kind = REAL64), dimension (:, :), allocatable :: output + + integer :: status, ncid, ivar, nf_type, nvdims + integer :: i, n1, n2, io_stat + integer, dimension(NF90_MAX_VAR_DIMS) :: dimids, idims + character (len = :), allocatable :: msg + + + ! Open file + status = nf90_open (trim(file_name), NF90_NOWRITE, ncid) + call Check_status (status) + + ! Get var + status = nf90_inq_varid (ncid, trim (var_name), ivar) + call Check_status (status) + status = nf90_inquire_variable (ncid, ivar, xtype = nf_type, ndims = nvdims, dimids = dimids) + call Check_status (status) + + if (nf_type == NF90_FLOAT .or. nf_type == NF90_DOUBLE) then + if (nvdims == 2) then + do i = 1, nvdims + status = nf90_inquire_dimension (ncid, dimids(i), len = idims(i)) + call Check_status (status) + end do + + n1 = idims(1) + n2 = idims(2) + allocate (output(n1, n2), stat = io_stat) + if (io_stat /= 0) Call Stop_simulation ('Problems allocating output variable') + + status = nf90_get_var (ncid, ivar, output) + call Check_status (status) + else + msg = 'FATAL ERROR: ' // trim (var_name) // ' is not a 2D array' + call Stop_simulation (msg) + end if + else + msg = 'FATAL ERROR: ' // trim (var_name) // ' is not a floating-point variable' + call Stop_simulation (msg) + end if + + ! Closing file + status = nf90_close (ncid) + call Check_status (status) + + end subroutine Get_netcdf_var_real64_2d + subroutine Get_netcdf_var_real32_3d (file_name, var_name, output) use netcdf diff --git a/physics/fmc_wrffire_mod.F90 b/physics/fmc_wrffire_mod.F90 index 7110ac9..a7dab51 100644 --- a/physics/fmc_wrffire_mod.F90 +++ b/physics/fmc_wrffire_mod.F90 @@ -13,10 +13,10 @@ 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 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 - use state_mod, only: state_fire_t - use namelist_mod, only : namelist_t use fmc_mod, only : fmc_t use fuel_mod, only : fuel_t @@ -24,10 +24,14 @@ module fmc_wrffire_mod private - public :: fmc_wrffire_t + public :: fmc_wrffire_t, NAME_VAR_FMC_GC, NAME_ATT_FMOIST_LASTTIME, NAME_ATT_FMOIST_NEXTTIME integer, parameter :: MOISTURE_CLASSES = 5, NUM_FMEP = 2 real, parameter :: FMEP_DECAY_TLAG = 999999 ! time constant of assimilated adjustments of equilibria decay + character (len = *), parameter :: NAME_DIM_MOISTURE_CLASS = 'moisture_class' + 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' type, extends(fmc_t) :: fmc_wrffire_t real :: fmoist_lasttime, fmoist_nexttime, dt_moisture @@ -42,6 +46,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 :: Add_restart_vars => Add_restart_vars procedure, public :: Init => Init_fmc_wrffire end type fmc_wrffire_t @@ -64,6 +69,50 @@ module fmc_wrffire_mod contains + subroutine Add_restart_vars (this, file_restart, cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, rank) + +#ifdef DM_PARALLEL + use mpi +#endif + implicit none + + class (fmc_wrffire_t), intent (in) :: this + character (len = *), intent (in) :: file_restart + integer, intent (in) :: cfbm_comm + integer, intent (in) :: nx + integer, intent (in) :: ny + integer, intent (in) :: ifps + integer, intent (in) :: ifpe + integer, intent (in) :: jfps + integer, intent (in) :: jfpe + integer, intent (in) :: rank + + character (len = 32), dimension(3) :: dim_names_fmc_gc + integer :: ierr, n_moisture_classes + + + if (.not. allocated (this%fmc_gc)) call Stop_simulation ('FMC restart write requires allocated fmc_gc') + + n_moisture_classes = size (this%fmc_gc, 2) + if (rank == 0) call Add_netcdf_dim (file_restart, NAME_DIM_MOISTURE_CLASS, n_moisture_classes) +#ifdef DM_PARALLEL + call MPI_Barrier (cfbm_comm, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with MPI_Barrier for FMC restart write') +#endif + + dim_names_fmc_gc(1) = NAME_DIM_X + dim_names_fmc_gc(2) = NAME_DIM_MOISTURE_CLASS + dim_names_fmc_gc(3) = NAME_DIM_Y + call Add_netcdf_var_mpi (file_restart, dim_names_fmc_gc, cfbm_comm, nx, ny, n_moisture_classes, & + 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_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 + + end subroutine Add_restart_vars + subroutine Init_fmc_wrffire (this, fuels, fuelmc_g, fuelmc_g_live, ifms, ifme, jfms, jfme, itimestep, dt) implicit none diff --git a/share/fmc_mod.F90 b/share/fmc_mod.F90 index 5500968..165ee5f 100644 --- a/share/fmc_mod.F90 +++ b/share/fmc_mod.F90 @@ -14,10 +14,25 @@ module fmc_mod type, abstract :: fmc_t contains procedure (Advance_fmc_model), deferred :: Advance_fmc_model + procedure (Add_restart_vars), deferred :: Add_restart_vars procedure (Init), deferred :: Init end type fmc_t abstract interface + subroutine Add_restart_vars (this, file_restart, cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, rank) + import :: fmc_t + class (fmc_t), intent (in) :: this + character (len = *), intent (in) :: file_restart + integer, intent (in) :: cfbm_comm + integer, intent (in) :: nx + integer, intent (in) :: ny + integer, intent (in) :: ifps + integer, intent (in) :: ifpe + integer, intent (in) :: jfps + integer, intent (in) :: jfpe + integer, intent (in) :: rank + end subroutine Add_restart_vars + subroutine Init (this, fuels, fuelmc_g, fuelmc_g_live, ifms, ifme, jfms, jfme, itimestep, dt) import :: fmc_t, fuel_t class (fmc_t), intent (in out) :: this diff --git a/share/mpi_mod.F90 b/share/mpi_mod.F90 index 34291f8..283dea6 100644 --- a/share/mpi_mod.F90 +++ b/share/mpi_mod.F90 @@ -1,6 +1,6 @@ module mpi_mod - use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT + use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT, REAL64 implicit none @@ -8,7 +8,7 @@ module mpi_mod public :: Calc_tasks_in_x_and_y, Calc_patch_dims, Gather_var2d, Do_halo_exchange, Do_halo_exchange_with_corners, & Max_across_mpi_tasks, Sum_across_mpi_tasks, Distribute_var2d, Min_across_mpi_tasks, Convert_mpi_comm_to_f08, & - Print_cart_info, topology_dim_order + Distribute_global_var2d, Print_cart_info, topology_dim_order ! --- Topology Ordering Flag --- ! 0 = Default (Fortran order: X is Dim 0, Y is Dim 1) @@ -16,6 +16,18 @@ module mpi_mod integer :: topology_dim_order = 0 ! ------------------------------ + interface Distribute_global_var2d + module procedure Distribute_global_var2d_real32, Distribute_global_var2d_real64 + end interface Distribute_global_var2d + + interface Gather_var2d + module procedure Gather_var2d_real32, Gather_var2d_real64 + end interface Gather_var2d + + interface Do_halo_exchange_with_corners + module procedure Do_halo_exchange_with_corners_real32, Do_halo_exchange_with_corners_real64 + end interface Do_halo_exchange_with_corners + contains subroutine Calc_patch_dims (nx, ny, px, py, coords, istart, iend, jstart, jend) @@ -204,6 +216,146 @@ subroutine Distribute_var2d (local_field, ips, ipe, jps, jpe, comm) end subroutine Distribute_var2d + subroutine Distribute_global_var2d_real32 (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_global, var2d_local) + +#ifdef DM_PARALLEL + use mpi +#endif + implicit none + + integer, intent (in) :: cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe + real, dimension(:, :), allocatable, intent (in out) :: var2d_global + real, dimension(ifps:ifpe, jfps:jfpe), intent (out) :: var2d_local + + integer :: i, ierr, j, kbuf, nprocs, nx_local, ny_local, r, rank + integer, dimension(:), allocatable :: all_ifps, all_ifpe, all_jfps, all_jfpe, displs, sendcounts + real, dimension(:), allocatable :: sendbuf + + +#ifdef DM_PARALLEL + call Mpi_comm_rank (cfbm_comm, rank, ierr) + call Mpi_comm_size (cfbm_comm, nprocs, ierr) + + nx_local = ifpe - ifps + 1 + ny_local = jfpe - jfps + 1 + + if (rank == 0) then + allocate (all_ifps(nprocs), all_ifpe(nprocs), all_jfps(nprocs), all_jfpe(nprocs)) + else + allocate (all_ifps(1), all_ifpe(1), all_jfps(1), all_jfpe(1)) + end if + + call MPI_Gather (ifps, 1, MPI_INTEGER, all_ifps, 1, MPI_INTEGER, 0, cfbm_comm, ierr) + call MPI_Gather (ifpe, 1, MPI_INTEGER, all_ifpe, 1, MPI_INTEGER, 0, cfbm_comm, ierr) + call MPI_Gather (jfps, 1, MPI_INTEGER, all_jfps, 1, MPI_INTEGER, 0, cfbm_comm, ierr) + call MPI_Gather (jfpe, 1, MPI_INTEGER, all_jfpe, 1, MPI_INTEGER, 0, cfbm_comm, ierr) + + if (rank == 0) then + allocate (sendcounts(nprocs), displs(nprocs)) + do r = 1, nprocs + sendcounts(r) = (all_ifpe(r) - all_ifps(r) + 1) * (all_jfpe(r) - all_jfps(r) + 1) + end do + + displs(1) = 0 + do r = 2, nprocs + displs(r) = displs(r - 1) + sendcounts(r - 1) + end do + + allocate (sendbuf(sum (sendcounts))) + kbuf = 0 + do r = 1, nprocs + do j = all_jfps(r), all_jfpe(r) + do i = all_ifps(r), all_ifpe(r) + kbuf = kbuf + 1 + sendbuf(kbuf) = var2d_global(i, j) + end do + end do + end do + else + allocate (sendcounts(1), displs(1), sendbuf(1)) + end if + + call MPI_Scatterv (sendbuf, sendcounts, displs, MPI_REAL, var2d_local, nx_local * ny_local, MPI_REAL, 0, cfbm_comm, ierr) + + deallocate (all_ifps, all_ifpe, all_jfps, all_jfpe, sendcounts, displs, sendbuf) +#else + var2d_local(1:nx, 1:ny) = var2d_global(1:nx, 1:ny) +#endif + + end subroutine Distribute_global_var2d_real32 + + subroutine Distribute_global_var2d_real64 (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_global, var2d_local) + +#ifdef DM_PARALLEL + use mpi_f08 +#endif + implicit none + +#ifdef DM_PARALLEL + type (MPI_Comm) :: cfbm_comm_f08 +#endif + integer, intent (in) :: cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe + real (kind = REAL64), dimension(:, :), allocatable, intent (in out) :: var2d_global + real (kind = REAL64), dimension(ifps:ifpe, jfps:jfpe), intent (out) :: var2d_local + + integer :: i, ierr, j, kbuf, nprocs, nx_local, ny_local, r, rank + integer, dimension(:), allocatable :: all_ifps, all_ifpe, all_jfps, all_jfpe, displs, sendcounts + real (kind = REAL64), dimension(:), allocatable :: sendbuf + + +#ifdef DM_PARALLEL + call Convert_mpi_comm_to_f08 (cfbm_comm, cfbm_comm_f08) + call Mpi_comm_rank (cfbm_comm_f08, rank, ierr) + call Mpi_comm_size (cfbm_comm_f08, nprocs, ierr) + + nx_local = ifpe - ifps + 1 + ny_local = jfpe - jfps + 1 + + if (rank == 0) then + allocate (all_ifps(nprocs), all_ifpe(nprocs), all_jfps(nprocs), all_jfpe(nprocs)) + else + allocate (all_ifps(1), all_ifpe(1), all_jfps(1), all_jfpe(1)) + end if + + call MPI_Gather (ifps, 1, MPI_INTEGER, all_ifps, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Gather (ifpe, 1, MPI_INTEGER, all_ifpe, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Gather (jfps, 1, MPI_INTEGER, all_jfps, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + call MPI_Gather (jfpe, 1, MPI_INTEGER, all_jfpe, 1, MPI_INTEGER, 0, cfbm_comm_f08, ierr) + + if (rank == 0) then + allocate (sendcounts(nprocs), displs(nprocs)) + do r = 1, nprocs + sendcounts(r) = (all_ifpe(r) - all_ifps(r) + 1) * (all_jfpe(r) - all_jfps(r) + 1) + end do + + displs(1) = 0 + do r = 2, nprocs + displs(r) = displs(r - 1) + sendcounts(r - 1) + end do + + allocate (sendbuf(sum (sendcounts))) + kbuf = 0 + do r = 1, nprocs + do j = all_jfps(r), all_jfpe(r) + do i = all_ifps(r), all_ifpe(r) + kbuf = kbuf + 1 + sendbuf(kbuf) = var2d_global(i, j) + end do + end do + end do + else + allocate (sendcounts(1), displs(1), sendbuf(1)) + end if + + call MPI_Scatterv (sendbuf, sendcounts, displs, MPI_DOUBLE_PRECISION, var2d_local, nx_local * ny_local, MPI_DOUBLE_PRECISION, 0, cfbm_comm_f08, ierr) + + deallocate (all_ifps, all_ifpe, all_jfps, all_jfpe, sendcounts, displs, sendbuf) +#else + var2d_local(1:nx, 1:ny) = var2d_global(1:nx, 1:ny) +#endif + + end subroutine Distribute_global_var2d_real64 + subroutine Do_halo_exchange (patch, ims, ime, jms, jme, ips, ipe, jps, jpe, nghost, cart_comm) ! It does not update the corners in the halo @@ -364,7 +516,7 @@ subroutine Do_halo_exchange (patch, ims, ime, jms, jme, ips, ipe, jps, jpe, ngho end subroutine Do_halo_exchange - subroutine Do_halo_exchange_with_corners (patch, ims, ime, jms, jme, ips, ipe, jps, jpe, nghost, cart_comm) + subroutine Do_halo_exchange_with_corners_real32 (patch, ims, ime, jms, jme, ips, ipe, jps, jpe, nghost, cart_comm) #ifdef DM_PARALLEL use mpi @@ -520,9 +672,173 @@ subroutine Do_halo_exchange_with_corners (patch, ims, ime, jms, jme, ips, ipe, j deallocate (sendbuf_up, recvbuf_down, sendbuf_down, recvbuf_up) #endif - end subroutine Do_halo_exchange_with_corners + end subroutine Do_halo_exchange_with_corners_real32 + + subroutine Do_halo_exchange_with_corners_real64 (patch, ims, ime, jms, jme, ips, ipe, jps, jpe, nghost, cart_comm) + +#ifdef DM_PARALLEL + use mpi_f08 +#endif + + implicit none + +#ifdef DM_PARALLEL + type (MPI_Comm) :: cart_comm_f08 +#endif + integer, intent (in) :: cart_comm, ims, ime, jms, jme, ips, ipe, jps, jpe, nghost + real (kind = REAL64), dimension(ims:ime, jms:jme), intent (inout) :: patch + + integer :: ierr, nbr_left, nbr_right, nbr_up, nbr_down, tag_base + integer :: nx, ny, rank, i, j, k +#ifdef DM_PARALLEL + type (MPI_Request), dimension(8) :: reqs +#endif + integer, dimension(2) :: coords + real (kind = REAL64), dimension(:), allocatable :: sendbuf_right, recvbuf_left, sendbuf_left, recvbuf_right + real (kind = REAL64), dimension(:), allocatable :: sendbuf_up, recvbuf_down, sendbuf_down, recvbuf_up + + ! New variables to handle the topology mapping + integer :: dim_x, dim_y + +#ifdef DM_PARALLEL + call Convert_mpi_comm_to_f08 (cart_comm, cart_comm_f08) + call MPI_Comm_rank(cart_comm_f08, rank, ierr) + call MPI_Cart_coords(cart_comm_f08, rank, 2, coords, ierr) + + nx = ipe - ips + 1 + ny = jpe - jps + 1 + tag_base = 1000 + + ! --------------------------------------------------------- + ! ADAPT DIMENSIONS BASED ON MODULE FLAG + ! --------------------------------------------------------- + if (topology_dim_order == 0) then + ! Default internal behavior [X, Y] + dim_x = 0 + dim_y = 1 + else + ! External framework behavior [Y, X] + dim_x = 1 + dim_y = 0 + end if + + ! Neighbor ranks using mapped dimensions + call MPI_Cart_shift (cart_comm_f08, dim_x, 1, nbr_left, nbr_right, ierr) + call MPI_Cart_shift (cart_comm_f08, dim_y, 1, nbr_down, nbr_up, ierr) + ! --------------------------------------------------------- + + ! Halo exchange in X direction + allocate (sendbuf_right(ny * nghost), recvbuf_left(ny * nghost)) + allocate (sendbuf_left(ny * nghost), recvbuf_right(ny * nghost)) + + ! Pack right send buffer + k = 0 + do j = jps, jpe + do i = ipe - nghost + 1, ipe + k = k + 1 + sendbuf_right(k) = patch(i, j) + end do + end do + + ! Pack left send buffer + k = 0 + do j = jps, jpe + do i = ips, ips + nghost - 1 + k = k + 1 + sendbuf_left(k) = patch(i, j) + end do + end do + + ! Exchange left/right + call MPI_Irecv (recvbuf_left, ny * nghost, MPI_DOUBLE_PRECISION, nbr_left, tag_base + 0, cart_comm_f08, reqs(1), ierr) + call MPI_Isend (sendbuf_right, ny * nghost, MPI_DOUBLE_PRECISION, nbr_right, tag_base + 0, cart_comm_f08, reqs(2), ierr) + call MPI_Irecv (recvbuf_right, ny * nghost, MPI_DOUBLE_PRECISION, nbr_right, tag_base + 1, cart_comm_f08, reqs(3), ierr) + call MPI_Isend (sendbuf_left, ny * nghost, MPI_DOUBLE_PRECISION, nbr_left, tag_base + 1, cart_comm_f08, reqs(4), ierr) + + call MPI_Waitall (4, reqs, MPI_STATUSES_IGNORE, ierr) + + ! Unpack left halo + if (nbr_left /= MPI_PROC_NULL) then + k = 0 + do j = jps, jpe + do i = ips - nghost, ips - 1 + k = k + 1 + patch(i, j) = recvbuf_left(k) + end do + end do + end if + + ! Unpack right halo + if (nbr_right /= MPI_PROC_NULL) then + k = 0 + do j = jps, jpe + do i = ipe + 1, ipe + nghost + k = k + 1 + patch(i, j) = recvbuf_right(k) + end do + end do + end if + + deallocate (sendbuf_right, recvbuf_left, sendbuf_left, recvbuf_right) + + ! Exchange in Y direction (including halos in X direction) + allocate (sendbuf_up((nx + 2 * nghost) * nghost), recvbuf_down((nx + 2 * nghost) * nghost)) + allocate (sendbuf_down((nx + 2 * nghost) * nghost), recvbuf_up((nx + 2 * nghost) * nghost)) + + ! Pack up send buffer + k = 0 + do j = jpe - nghost + 1, jpe + do i = ips - nghost, ipe + nghost + k = k + 1 + sendbuf_up(k) = patch(i, j) + end do + end do + + ! Pack down send buffer + k = 0 + do j = jps, jps + nghost - 1 + do i = ips - nghost, ipe + nghost + k = k + 1 + sendbuf_down(k) = patch(i, j) + end do + end do + + ! Exchange up/down + call MPI_Irecv (recvbuf_down, (nx + 2 * nghost) * nghost, MPI_DOUBLE_PRECISION, nbr_down, tag_base + 2, cart_comm_f08, reqs(1), ierr) + call MPI_Isend (sendbuf_up, (nx + 2 * nghost) * nghost, MPI_DOUBLE_PRECISION, nbr_up, tag_base + 2, cart_comm_f08, reqs(2), ierr) + call MPI_Irecv (recvbuf_up, (nx + 2 * nghost) * nghost, MPI_DOUBLE_PRECISION, nbr_up, tag_base + 3, cart_comm_f08, reqs(3), ierr) + call MPI_Isend (sendbuf_down, (nx + 2 * nghost) * nghost, MPI_DOUBLE_PRECISION, nbr_down, tag_base + 3, cart_comm_f08, reqs(4), ierr) + + call MPI_Waitall(4, reqs, MPI_STATUSES_IGNORE, ierr) + + ! Unpack down halo + if (nbr_down /= MPI_PROC_NULL) then + k = 0 + do j = jps - nghost, jps - 1 + do i = ips - nghost, ipe + nghost + k = k + 1 + patch(i, j) = recvbuf_down(k) + end do + end do + end if + + ! Unpack up halo + if (nbr_up /= MPI_PROC_NULL) then + k = 0 + do j = jpe + 1, jpe + nghost + do i = ips - nghost, ipe + nghost + k = k + 1 + patch(i, j) = recvbuf_up(k) + end do + end do + end if + + deallocate (sendbuf_up, recvbuf_down, sendbuf_down, recvbuf_up) +#endif + + end subroutine Do_halo_exchange_with_corners_real64 - subroutine Gather_var2d (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_local, var2d_global) + subroutine Gather_var2d_real32 (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_local, var2d_global) #ifdef DM_PARALLEL use mpi_f08 @@ -578,7 +894,65 @@ subroutine Gather_var2d (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_local, end if #endif - end subroutine Gather_var2d + end subroutine Gather_var2d_real32 + + subroutine Gather_var2d_real64 (cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, var2d_local, var2d_global) + +#ifdef DM_PARALLEL + use mpi_f08 +#endif + + implicit none + + integer, intent (in) :: cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe + real (kind = REAL64), dimension(ifps:ifpe, jfps:jfpe), intent (in) :: var2d_local + real (kind = REAL64), dimension(nx, ny), intent (out) :: var2d_global + + real (kind = REAL64), dimension(:, :), allocatable :: buf + integer :: rank, ierr, ntasks, src, i_start, j_start, nx_local, ny_local, nxl, nyl + + integer, dimension(4) :: metadata +#ifdef DM_PARALLEL + type(MPI_Comm) :: cfbm_comm_f08 +#endif + + +#ifdef DM_PARALLEL + call Convert_mpi_comm_to_f08 (cfbm_comm, cfbm_comm_f08) + + call Mpi_comm_size (cfbm_comm_f08, ntasks, ierr) + call Mpi_comm_rank (cfbm_comm_f08, rank, ierr) + + nx_local = ifpe - ifps + 1 + ny_local = jfpe - jfps + 1 + if (rank /= 0) then + metadata = [ifps, jfps, nx_local, ny_local] + call Mpi_Send (metadata, 4, MPI_INTEGER, 0, 200, cfbm_comm_f08, ierr) + call Mpi_Send (var2d_local, nx_local * ny_local, MPI_DOUBLE_PRECISION, 0, 100, cfbm_comm_f08, ierr) + else + ! Rank 0 places its own data + var2d_global(ifps:ifpe, jfps:jfpe) = var2d_local + + ! Receive from others + do src = 1, ntasks - 1 + ! Receive metadata + call Mpi_Recv (metadata, 4, MPI_INTEGER, src, 200, cfbm_comm_f08, MPI_STATUS_IGNORE, ierr) + i_start = metadata(1) + j_start = metadata(2) + nxl = metadata(3) + nyl = metadata(4) + allocate (buf(nxl, nyl)) + call Mpi_Recv (buf, nxl * nyl, MPI_DOUBLE_PRECISION, src, 100, cfbm_comm_f08, MPI_STATUS_IGNORE, ierr) + + ! Place in global array + var2d_global(i_start:i_start + nxl - 1, j_start:j_start + nyl - 1) = buf + + deallocate (buf) + end do + end if +#endif + + end subroutine Gather_var2d_real64 subroutine Max_across_mpi_tasks (local_max, cart_comm, global_max) diff --git a/state/state_mod.F90 b/state/state_mod.F90 index 05804c1..190bb00 100644 --- a/state/state_mod.F90 +++ b/state/state_mod.F90 @@ -13,8 +13,8 @@ module state_mod use ignition_line_mod, only : ignition_line_t use interp_mod, only : VINTERP_WINDS_FROM_3D_WINDS, VINTERP_WINDS_FROM_10M_WINDS use namelist_mod, only : namelist_t - use netcdf_mod, only : Create_netcdf_file, Add_netcdf_dim, Add_netcdf_var, Add_netcdf_var_mpi, NAME_DIM_X, NAME_DIM_Y, & - NF90_FILL_FLOAT + use netcdf_mod, only : Create_netcdf_file, Add_netcdf_dim, Add_netcdf_var, Add_netcdf_var_mpi, & + NAME_DIM_X, NAME_DIM_Y, NF90_FILL_FLOAT use proj_lc_mod, only : proj_lc_t use ros_mod, only : ros_t use stderrout_mod, only : Stop_simulation, Print_message @@ -28,7 +28,7 @@ module state_mod public :: state_fire_t, N_POINTS_IN_HALO - integer, parameter :: INIT_MODE_NONE = 0, INIT_MODE_GEOGRID = 1, INIT_MODE_WRF = 2, INIT_MODE_IDEAL = 3 + integer, parameter :: INIT_MODE_NONE = 0, INIT_MODE_GEOGRID = 1, INIT_MODE_WRF = 2, INIT_MODE_IDEAL = 3, INIT_MODE_RESTART = 4 integer, parameter :: N_POINTS_IN_HALO = 5, N_DIMS = 2 logical, dimension(2), parameter :: PERIODS = [ .false., .false. ] @@ -114,6 +114,8 @@ module state_mod ! Output integer :: output_level + integer :: restart_step_interval = -1 + logical :: is_restart_output_initialized = .false. ! Resolved initialization path; supplied grid data can override ideal_opt. integer :: init_mode = INIT_MODE_NONE logical :: fire_perimeter_ignited = .false. @@ -127,8 +129,10 @@ module state_mod contains procedure, public :: Allocate_vars => Allocate_vars procedure, public :: Apply_wafs => Apply_wafs + procedure, public :: Check_model_clock => Check_model_clock procedure, public :: Convert_sb_to_ander => Convert_scottburgan_to_anderson procedure, public :: Handle_output => Handle_output + procedure, public :: Handle_restart => Handle_restart procedure, public :: Handle_wrfdata_update => Handle_wrfdata_update procedure, public :: Init_fuel_vars => Init_fuel_vars procedure, public :: Initialization => Init_domain @@ -141,12 +145,54 @@ module state_mod procedure :: Interpolate_vars_atm_to_fire => Interpolate_vars_atm_to_fire procedure, public :: Print => Print_domain ! private procedure, public :: Print_tiles => Print_tiles + procedure, public :: Read_restart => Read_restart procedure, public :: Save_state => Save_state procedure, public :: Set_vars_to_default => Set_vars_to_default procedure, public :: Set_mpi_comm_cfbm => Set_mpi_comm_cfbm procedure, public :: Set_time_stamps => Set_time_stamps + procedure, public :: Write_restart => Write_restart end type state_fire_t + interface + module function Build_restart_file_name (restart_datetime) result (file_restart) + character (len = *), intent (in) :: restart_datetime + character (len = :), allocatable :: file_restart + end function Build_restart_file_name + + module subroutine Init_restart_dimensions (this, restart_file, ids0, ide0, jds0, jde0, ips, ipe, jps, jpe) + class (state_fire_t), intent (in out) :: this + character (len = *), intent (in) :: restart_file + integer, intent (out) :: ids0, ide0, jds0, jde0, ips, ipe, jps, jpe + end subroutine Init_restart_dimensions + + module subroutine Init_restart_projection (this, restart_file, proj) + class (state_fire_t), intent (in out) :: this + character (len = *), intent (in) :: restart_file + type (proj_lc_t), intent (out) :: proj + end subroutine Init_restart_projection + + module subroutine Read_restart_static_fields (this, restart_file) + class (state_fire_t), intent (in out) :: this + character (len = *), intent (in) :: restart_file + end subroutine Read_restart_static_fields + + module subroutine Read_restart (this, config_flags) + class (state_fire_t), intent (in out) :: this + type (namelist_t), intent (in) :: config_flags + end subroutine Read_restart + + module subroutine Handle_restart (this, config_flags, initialize) + class (state_fire_t), intent (in out) :: this + type (namelist_t), intent (in) :: config_flags + logical, intent (in), optional :: initialize + end subroutine Handle_restart + + module subroutine Write_restart (this, config_flags) + class (state_fire_t), intent (in) :: this + type (namelist_t), intent (in) :: config_flags + end subroutine Write_restart + end interface + contains subroutine Allocate_vars (this, ifms, ifme, jfms, jfme) @@ -260,6 +306,21 @@ subroutine Convert_scottburgan_to_anderson (this) end subroutine Convert_scottburgan_to_anderson + subroutine Check_model_clock (this) + + implicit none + + class (state_fire_t), intent (in) :: this + + type (datetime_t) :: datetime_check + + + datetime_check = this%datetime_start + call datetime_check%Add_seconds (this%itimestep * this%dt) + if (datetime_check /= this%datetime_now) call Stop_simulation ('Model clock is inconsistent with itimestep and dt') + + end subroutine Check_model_clock + subroutine Handle_output (this, config_flags) implicit none @@ -331,26 +392,40 @@ subroutine Init_domain (this, config_flags, geogrid, & real, intent (in), optional :: cen_lat, cen_lon, truelat1, truelat2, stand_lon, dx, dy real, dimension(:, :), intent (in), optional :: nfuel_cat, zsf, dzdxf, dzdyf + type (datetime_t) :: datetime_restart type (proj_lc_t) :: proj logical, parameter :: DEBUG_LOCAL = .false. - integer :: ids0, ide0, jds0, jde0, i, j, px, py, ntasks, ierr, cart_comm, rank, ips, ipe, jps, jpe, is_lfn_init_allocated + logical :: has_wrf_metadata + integer :: ids0, ide0, jds0, jde0, i, j, px, py, ntasks, ierr, cart_comm, rank, ips, ipe, jps, jpe, & + is_lfn_init_allocated integer, dimension(2) :: coords character (len = 300) :: msg + character (len = :), allocatable :: file_restart if (DEBUG_LOCAL) call Print_message ('Entering Init_domain...') + has_wrf_metadata = present (ifds) .and. present (ifde) .and. present (ifms) .and. present (ifme) .and. & + present (ifps) .and. present (ifpe) .and. present (jfds) .and. present (jfde) .and. present (jfms) .and. & + present (jfme) .and. present (jfps) .and. present (jfpe) .and. present (kfds) .and. present (kfde) .and. & + present (kfms) .and. present (kfme) .and. present (kfps) .and. present (kfpe) .and. present (kfts) .and. & + present (kfte) .and. present (ide) .and. present (jde) .and. present (cen_lat) .and. present (cen_lon) .and. & + present (truelat1) .and. present (truelat2) .and. present (stand_lon) .and. present (dx) .and. present (dy) .and. & + present (sr_x) .and. present (sr_y) .and. present (nfuel_cat) .and. present (zsf) .and. present (dzdxf) .and. & + present (dzdyf) + + if (config_flags%ideal_opt == 0 .and. config_flags%restart .and. .not. present (geogrid) .and. & + .not. has_wrf_metadata) then + datetime_restart = datetime_t (config_flags%start_year, config_flags%start_month, config_flags%start_day, & + config_flags%start_hour, config_flags%start_minute, config_flags%start_second) + file_restart = Build_restart_file_name (datetime_restart%datetime) + end if + this%init_mode = INIT_MODE_NONE if (config_flags%ideal_opt == 1) this%init_mode = INIT_MODE_IDEAL if (present (geogrid)) this%init_mode = INIT_MODE_GEOGRID - if (present (ifds) .and. present (ifde) .and. present (ifms) .and. present (ifme) .and. present (ifps) .and. present (ifpe) .and. & - present (jfds) .and. present (jfde) .and. present (jfms) .and. present (jfme) .and. present (jfps) .and. present (jfpe) .and. & - present (kfds) .and. present (kfde) .and. present (kfms) .and. present (kfme) .and. present (kfps) .and. present (kfpe) .and. & - present (kfts) .and. present (kfte) .and. present (ide) .and. present (jde) .and. & - present (cen_lat) .and. present (cen_lon) .and. present (truelat1) .and. present (truelat2) .and. present (stand_lon) .and. & - present (dx) .and. present (dy) .and. present (sr_x) .and. present (sr_y) .and. present (nfuel_cat) .and. present (zsf) .and. & - present (dzdxf) .and. present (dzdyf)) & - this%init_mode = INIT_MODE_WRF + if (allocated (file_restart)) this%init_mode = INIT_MODE_RESTART + if (has_wrf_metadata) this%init_mode = INIT_MODE_WRF if (this%init_mode == INIT_MODE_NONE) & call Stop_simulation ('Not enough information to initialize domain') @@ -358,7 +433,7 @@ subroutine Init_domain (this, config_flags, geogrid, & ! Set dimensions if (DEBUG_LOCAL) call Print_message (' Setting dimensions...') Set_dims: select case (this%init_mode) - case (INIT_MODE_GEOGRID, INIT_MODE_IDEAL) + case (INIT_MODE_GEOGRID, INIT_MODE_IDEAL, INIT_MODE_RESTART) if (this%init_mode == INIT_MODE_GEOGRID) then @@ -453,6 +528,9 @@ subroutine Init_domain (this, config_flags, geogrid, & jps = jds0 jpe = jde0 #endif + else if (this%init_mode == INIT_MODE_RESTART) then + + call Init_restart_dimensions (this, trim (file_restart), ids0, ide0, jds0, jde0, ips, ipe, jps, jpe) end if this%ifds = ids0 @@ -576,6 +654,9 @@ subroutine Init_domain (this, config_flags, geogrid, & call this%Init_latlons (proj) + case (INIT_MODE_RESTART) + call Init_restart_projection (this, trim (file_restart), proj) + case default call Stop_simulation ('Not ready to complete fire state initialization 2') @@ -626,12 +707,15 @@ subroutine Init_domain (this, config_flags, geogrid, & if (config_flags%fire_is_real_perim) & call Stop_simulation ('Not ready to initialize from fire perimeter in idealized mode') + case (INIT_MODE_RESTART) + call Read_restart_static_fields (this, trim (file_restart)) + case default call Stop_simulation ('Not ready to complete fire state initialization 3') end select Set_topo_fuels - if (config_flags%fuel_opt == FUEL_ANDERSON) call this%Convert_sb_to_ander () + if (config_flags%fuel_opt == FUEL_ANDERSON .and. this%init_mode /= INIT_MODE_RESTART) call this%Convert_sb_to_ander () ! Set clock if (DEBUG_LOCAL) call Print_message (' Setting clock...') diff --git a/state/state_restart_submod.F90 b/state/state_restart_submod.F90 new file mode 100644 index 0000000..48d41fc --- /dev/null +++ b/state/state_restart_submod.F90 @@ -0,0 +1,907 @@ + submodule (state_mod) state_restart_submod + +#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, intrinsic :: iso_fortran_env, only : INT32, REAL32 + + implicit none + + integer, parameter :: RESTART_IO_ROOT = 0 + real, parameter :: RESTART_INTERVAL_TOL = 1.0e-6 + + contains + + module function Build_restart_file_name (restart_datetime) result (file_restart) + + implicit none + + character (len = *), intent (in) :: restart_datetime + character (len = :), allocatable :: file_restart + + file_restart = 'fire_restart_'//trim (restart_datetime)//'.nc' + + end function Build_restart_file_name + + module subroutine Handle_restart (this, config_flags, initialize) + + implicit none + + class (state_fire_t), intent (in out) :: this + type (namelist_t), intent (in) :: config_flags + logical, intent (in), optional :: initialize + + real :: restart_interval_seconds, restart_interval_error + character (len = 256) :: msg + logical :: do_initialize + + + do_initialize = .false. + if (present (initialize)) do_initialize = initialize + + if (do_initialize .or. .not. this%is_restart_output_initialized) then + this%restart_step_interval = -1 + if (config_flags%restart_interval > 0) then + this%restart_step_interval = nint (real (config_flags%restart_interval) / this%dt) + restart_interval_seconds = this%restart_step_interval * this%dt + restart_interval_error = abs (restart_interval_seconds - real (config_flags%restart_interval)) + + if (this%restart_step_interval <= 0 .or. & + restart_interval_error > max (RESTART_INTERVAL_TOL, abs (real (config_flags%restart_interval)) * RESTART_INTERVAL_TOL)) then + write (msg, '(a, i0, a, f12.6)') 'restart_interval must map to an integer number of time steps: restart_interval = ', & + config_flags%restart_interval, ', dt = ', this%dt + call Stop_simulation (msg) + end if + end if + this%is_restart_output_initialized = .true. + end if + + if (do_initialize) return + + if (this%restart_step_interval > 0) then + if (mod (this%itimestep, this%restart_step_interval) == 0) call this%Write_restart (config_flags) + end if + + end subroutine Handle_restart + + module subroutine Init_restart_dimensions (this, restart_file, ids0, ide0, jds0, jde0, ips, ipe, jps, jpe) + + implicit none + + class (state_fire_t), intent (in out) :: this + character (len = *), intent (in) :: restart_file + integer, intent (out) :: ids0, ide0, jds0, jde0, ips, ipe, jps, jpe + + integer :: restart_nx, restart_ny + integer (kind = INT32) :: att_int32 + integer :: ntasks, px, py, ierr, cart_comm, rank + integer, dimension(2) :: coords + character (len = 300) :: msg + + + restart_nx = 0 + restart_ny = 0 + if (Is_restart_io_root (this)) then + call Get_netcdf_att (trim (restart_file), 'global', 'nx', att_int32) + restart_nx = att_int32 + call Get_netcdf_att (trim (restart_file), 'global', 'ny', att_int32) + restart_ny = att_int32 + end if + call Broadcast_restart_integer (this, restart_nx) + call Broadcast_restart_integer (this, restart_ny) + + ids0 = 1 + ide0 = restart_nx + jds0 = 1 + jde0 = restart_ny + +#ifdef DM_PARALLEL + if (.not. this%is_cfbm_comm_set) call Stop_simulation ('The MPI CFBM communicator has not been set') + + call Mpi_comm_size (this%cfbm_comm, ntasks, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems getting the number of MPI tasks') + this%ntasks = ntasks + + call Calc_tasks_in_x_and_y (this%ntasks, restart_nx, restart_ny, px, py) + this%px = px + this%py = py + write (msg, '(a25, 2(1x, i5))') 'MPI TASKS in x and y =', this%px, this%py + call Print_message (msg) + + call Mpi_cart_create (this%cfbm_comm, N_DIMS, [this%px, this%py], PERIODS, REORDER, cart_comm, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_cart_create') + this%cart_comm = cart_comm + + call Mpi_comm_rank (this%cart_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') + + call Mpi_cart_coords (this%cart_comm, rank, N_DIMS, coords, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_cart_coords') + + call Calc_patch_dims (restart_nx, restart_ny, this%px, this%py, coords, ips, ipe, jps, jpe) +#else + ips = ids0 + ipe = ide0 + jps = jds0 + jpe = jde0 +#endif + + end subroutine Init_restart_dimensions + + module subroutine Init_restart_projection (this, restart_file, proj) + + implicit none + + class (state_fire_t), intent (in out) :: this + character (len = *), intent (in) :: restart_file + type (proj_lc_t), intent (out) :: proj + + integer :: restart_map_proj, restart_sr_x, restart_sr_y + integer (kind = INT32) :: att_int32 + real (kind = REAL32) :: att_real32 + real :: restart_dx, restart_dy, restart_cen_lat, restart_cen_lon, restart_stand_lon, restart_true_lat_1, restart_true_lat_2 + + + restart_dx = 0.0 + restart_dy = 0.0 + restart_map_proj = 0 + restart_sr_x = 0 + restart_sr_y = 0 + restart_cen_lat = 0.0 + restart_cen_lon = 0.0 + restart_stand_lon = 0.0 + restart_true_lat_1 = 0.0 + restart_true_lat_2 = 0.0 + + if (Is_restart_io_root (this)) then + call Get_netcdf_att (trim (restart_file), 'global', 'dx', att_real32) + restart_dx = att_real32 + call Get_netcdf_att (trim (restart_file), 'global', 'dy', att_real32) + restart_dy = att_real32 + call Get_netcdf_att (trim (restart_file), 'global', 'map_proj', att_int32) + restart_map_proj = att_int32 + call Get_netcdf_att (trim (restart_file), 'global', 'sr_x', att_int32) + restart_sr_x = att_int32 + call Get_netcdf_att (trim (restart_file), 'global', 'sr_y', att_int32) + restart_sr_y = att_int32 + call Get_netcdf_att (trim (restart_file), 'global', 'cen_lat', att_real32) + restart_cen_lat = att_real32 + call Get_netcdf_att (trim (restart_file), 'global', 'cen_lon', att_real32) + restart_cen_lon = att_real32 + call Get_netcdf_att (trim (restart_file), 'global', 'stand_lon', att_real32) + restart_stand_lon = att_real32 + call Get_netcdf_att (trim (restart_file), 'global', 'true_lat_1', att_real32) + restart_true_lat_1 = att_real32 + call Get_netcdf_att (trim (restart_file), 'global', 'true_lat_2', att_real32) + restart_true_lat_2 = att_real32 + end if + + call Broadcast_restart_real (this, restart_dx) + call Broadcast_restart_real (this, restart_dy) + call Broadcast_restart_integer (this, restart_map_proj) + call Broadcast_restart_integer (this, restart_sr_x) + call Broadcast_restart_integer (this, restart_sr_y) + call Broadcast_restart_real (this, restart_cen_lat) + call Broadcast_restart_real (this, restart_cen_lon) + call Broadcast_restart_real (this, restart_stand_lon) + call Broadcast_restart_real (this, restart_true_lat_1) + call Broadcast_restart_real (this, restart_true_lat_2) + + this%dx = restart_dx + this%dy = restart_dy + + if (restart_map_proj /= 1) call Stop_simulation ('Restart map projection is not supported') + + if (restart_sr_x <= 0 .or. restart_sr_y <= 0) call Stop_simulation ('Restart subgrid ratios must be positive') + if (mod (this%nx, restart_sr_x) /= 0 .or. mod (this%ny, restart_sr_y) /= 0) & + call Stop_simulation ('Restart subgrid ratios do not divide fire-grid dimensions') + if (this%nx / restart_sr_x <= 1 .or. this%ny / restart_sr_y <= 1) & + call Stop_simulation ('Restart subgrid ratios imply invalid atmospheric dimensions') + + this%cen_lat = restart_cen_lat + this%cen_lon = restart_cen_lon + + proj = proj_lc_t (cen_lat = this%cen_lat , cen_lon = this%cen_lon, dx = this%dx * restart_sr_x, & + dy = this%dy * restart_sr_y, & + standard_lon = restart_stand_lon, true_lat_1 = restart_true_lat_1, & + true_lat_2 = restart_true_lat_2, nx = this%nx / restart_sr_x - 1, ny = this%ny / restart_sr_y - 1) + + call this%Init_latlons (proj, srx = restart_sr_x, sry = restart_sr_y) + + end subroutine Init_restart_projection + + module subroutine Read_restart_static_fields (this, restart_file) + + implicit none + + class (state_fire_t), intent (in out) :: this + character (len = *), intent (in) :: restart_file + + call Read_restart_field_2d (trim (restart_file), 'zsf', this%nx, this%ny, & + this%ifms, this%ifme, this%jfms, this%jfme, this%ifps, this%ifpe, this%jfps, this%jfpe, this%cfbm_comm, this%zsf) + call Read_restart_field_2d (trim (restart_file), 'dzdxf', this%nx, this%ny, & + this%ifms, this%ifme, this%jfms, this%jfme, this%ifps, this%ifpe, this%jfps, this%jfpe, this%cfbm_comm, this%dzdxf) + call Read_restart_field_2d (trim (restart_file), 'dzdyf', this%nx, this%ny, & + this%ifms, this%ifme, this%jfms, this%jfme, this%ifps, this%ifpe, this%jfps, this%jfpe, this%cfbm_comm, this%dzdyf) + call Read_restart_field_2d (trim (restart_file), 'nfuel_cat', this%nx, this%ny, & + this%ifms, this%ifme, this%jfms, this%jfme, this%ifps, this%ifpe, this%jfps, this%jfpe, this%cfbm_comm, this%nfuel_cat) + call Read_restart_field_2d (trim (restart_file), 'fz0', this%nx, this%ny, & + this%ifms, this%ifme, this%jfms, this%jfme, this%ifps, this%ifpe, this%jfps, this%jfpe, this%cfbm_comm, this%fz0) + + end subroutine Read_restart_static_fields + + subroutine Require_restart_mpi_comm (this) + + implicit none + + class (state_fire_t), intent (in) :: this + + +#ifdef DM_PARALLEL + if (.not. this%is_cfbm_comm_set) call Stop_simulation ('The MPI CFBM communicator has not been set') +#endif + + end subroutine Require_restart_mpi_comm + + function Is_restart_io_root (this) result (am_root) + + implicit none + + class (state_fire_t), intent (in) :: this + logical :: am_root + + integer :: rank, ierr + + +#ifdef DM_PARALLEL + call Require_restart_mpi_comm (this) + call Mpi_comm_rank (this%cfbm_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') + am_root = rank == RESTART_IO_ROOT +#else + am_root = .true. +#endif + + end function Is_restart_io_root + + subroutine Restart_io_barrier (this) + + implicit none + + class (state_fire_t), intent (in) :: this + + integer :: ierr + + +#ifdef DM_PARALLEL + call Require_restart_mpi_comm (this) + call MPI_Barrier (this%cfbm_comm, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with MPI_Barrier ') +#endif + + end subroutine Restart_io_barrier + + subroutine Broadcast_restart_integer (this, value) + + implicit none + + class (state_fire_t), intent (in) :: this + integer, intent (in out) :: value + + integer :: ierr + + +#ifdef DM_PARALLEL + call Require_restart_mpi_comm (this) + call MPI_Bcast (value, 1, MPI_INTEGER, RESTART_IO_ROOT, this%cfbm_comm, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with MPI_Bcast for restart integer metadata') +#endif + + end subroutine Broadcast_restart_integer + + subroutine Broadcast_restart_real (this, value) + + implicit none + + class (state_fire_t), intent (in) :: this + real, intent (in out) :: value + + integer :: packed_value + + +#ifdef DM_PARALLEL + packed_value = 0 + if (Is_restart_io_root (this)) packed_value = transfer (value, packed_value) + call Broadcast_restart_integer (this, packed_value) + value = transfer (packed_value, value) +#endif + + end subroutine Broadcast_restart_real + + module subroutine Read_restart (this, config_flags) + + implicit none + + class (state_fire_t), intent (in out) :: this + type (namelist_t), intent (in) :: config_flags + + type (datetime_t) :: datetime_restart, datetime_check + character (len = :), allocatable :: file_restart + integer (kind = INT32) :: att_int, restart_year, restart_month, restart_day, restart_hour, restart_minute, restart_second, & + start_year, start_month, start_day, start_hour, start_minute, start_second + integer (kind = INT32) :: perimeter_ignited + integer :: ij, expected_sr_x, expected_sr_y, ierr, rank + logical, parameter :: DEBUG_LOCAL = .false. + + + if (DEBUG_LOCAL) call Print_message ('Entering Read_restart...') + +#ifdef DM_PARALLEL + call Require_restart_mpi_comm (this) + call Mpi_comm_rank (this%cfbm_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') +#else + rank = 0 +#endif + + if (config_flags%ideal_opt /= 0 .and. config_flags%ideal_opt /= 1) & + call Stop_simulation ('Read_restart is implemented for standalone idealized and real runs only') + + datetime_restart = datetime_t (config_flags%start_year, config_flags%start_month, config_flags%start_day, & + config_flags%start_hour, config_flags%start_minute, config_flags%start_second) + file_restart = Build_restart_file_name (datetime_restart%datetime) + + if (rank == 0) then + call Is_netcdf_file_present (file_restart) + + call Validate_restart_integer (file_restart, 'restart_year', config_flags%start_year) + call Validate_restart_integer (file_restart, 'restart_month', config_flags%start_month) + call Validate_restart_integer (file_restart, 'restart_day', config_flags%start_day) + call Validate_restart_integer (file_restart, 'restart_hour', config_flags%start_hour) + call Validate_restart_integer (file_restart, 'restart_minute', config_flags%start_minute) + call Validate_restart_integer (file_restart, 'restart_second', config_flags%start_second) + + call Validate_restart_integer (file_restart, 'nx', this%nx) + call Validate_restart_integer (file_restart, 'ny', this%ny) + call Validate_restart_real (file_restart, 'dt', this%dt) + call Validate_restart_real (file_restart, 'dx', this%dx) + call Validate_restart_real (file_restart, 'dy', this%dy) + call Validate_restart_real (file_restart, 'cen_lat', this%cen_lat) + call Validate_restart_real (file_restart, 'cen_lon', this%cen_lon) + expected_sr_x = nint (this%proj%dx / this%dx) + expected_sr_y = nint (this%proj%dy / this%dy) + call Validate_restart_integer (file_restart, 'map_proj', 1) + call Validate_restart_integer (file_restart, 'sr_x', expected_sr_x) + call Validate_restart_integer (file_restart, 'sr_y', expected_sr_y) + call Validate_restart_integer (file_restart, 'ideal_opt', config_flags%ideal_opt) + 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) + 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) + call Validate_restart_integer (file_restart, 'fire_lsm_reinit_iter', config_flags%fire_lsm_reinit_iter) + call Validate_restart_real (file_restart, 'fire_viscosity', config_flags%fire_viscosity) + call Validate_restart_real (file_restart, 'fire_viscosity_bg', config_flags%fire_viscosity_bg) + call Validate_restart_integer (file_restart, 'fire_viscosity_ngp', config_flags%fire_viscosity_ngp) + call Validate_restart_real (file_restart, 'fire_viscosity_band', config_flags%fire_viscosity_band) + call Validate_restart_real (file_restart, 'reinit_pseudot_coef', config_flags%reinit_pseudot_coef) + + select case (config_flags%ideal_opt) + case (0) + call Validate_restart_real (file_restart, 'stand_lon', this%proj%standard_lon) + call Validate_restart_real (file_restart, 'true_lat_1', this%proj%true_lat_1) + call Validate_restart_real (file_restart, 'true_lat_2', this%proj%true_lat_2) + + case (1) + call Validate_restart_real (file_restart, 'stand_lon', config_flags%stand_lon) + call Validate_restart_real (file_restart, 'true_lat_1', config_flags%true_lat_1) + call Validate_restart_real (file_restart, 'true_lat_2', config_flags%true_lat_2) + end select + + call Get_netcdf_att (file_restart, 'global', 'start_year', start_year) + call Get_netcdf_att (file_restart, 'global', 'start_month', start_month) + call Get_netcdf_att (file_restart, 'global', 'start_day', start_day) + call Get_netcdf_att (file_restart, 'global', 'start_hour', start_hour) + call Get_netcdf_att (file_restart, 'global', 'start_minute', start_minute) + call Get_netcdf_att (file_restart, 'global', 'start_second', start_second) + + call Get_netcdf_att (file_restart, 'global', 'restart_year', restart_year) + call Get_netcdf_att (file_restart, 'global', 'restart_month', restart_month) + call Get_netcdf_att (file_restart, 'global', 'restart_day', restart_day) + call Get_netcdf_att (file_restart, 'global', 'restart_hour', restart_hour) + call Get_netcdf_att (file_restart, 'global', 'restart_minute', restart_minute) + call Get_netcdf_att (file_restart, 'global', 'restart_second', restart_second) + + call Get_netcdf_att (file_restart, 'global', 'fire_perimeter_ignited', perimeter_ignited) + if (perimeter_ignited /= 0 .and. perimeter_ignited /= 1) & + call Stop_simulation ('Restart fire_perimeter_ignited must be 0 or 1') + + call Get_netcdf_att (file_restart, 'global', 'itimestep', att_int) + end if + + call Broadcast_restart_integer (this, start_year) + call Broadcast_restart_integer (this, start_month) + call Broadcast_restart_integer (this, start_day) + call Broadcast_restart_integer (this, start_hour) + call Broadcast_restart_integer (this, start_minute) + call Broadcast_restart_integer (this, start_second) + call Broadcast_restart_integer (this, restart_year) + call Broadcast_restart_integer (this, restart_month) + call Broadcast_restart_integer (this, restart_day) + call Broadcast_restart_integer (this, restart_hour) + call Broadcast_restart_integer (this, restart_minute) + call Broadcast_restart_integer (this, restart_second) + call Broadcast_restart_integer (this, att_int) + this%itimestep = att_int + call Broadcast_restart_integer (this, perimeter_ignited) + this%fire_perimeter_ignited = perimeter_ignited == 1 + + this%datetime_start = datetime_t (start_year, start_month, start_day, start_hour, start_minute, start_second) + this%datetime_now = datetime_t (restart_year, restart_month, restart_day, restart_hour, restart_minute, restart_second) + + datetime_check = this%datetime_start + call datetime_check%Add_seconds (this%itimestep * this%dt) + if (datetime_check /= this%datetime_now) call Stop_simulation ('Restart clock is inconsistent with itimestep and dt') + + call Set_next_datetime_after_restart (this%datetime_next_output, config_flags%interval_output) + call Set_next_datetime_after_restart (this%datetime_next_atm_update, config_flags%interval_atm) + + call Read_restart_field (file_restart, 'lfn', this%lfn) + call Read_restart_field (file_restart, 'lfn_hist', this%lfn_hist) + call Read_restart_field (file_restart, 'lfn_0', this%lfn_0) + call Read_restart_field (file_restart, 'lfn_1', this%lfn_1) + call Read_restart_field (file_restart, 'lfn_2', this%lfn_2) + call Read_restart_field (file_restart, 'lfn_s0', this%lfn_s0) + call Read_restart_field (file_restart, 'lfn_s1', this%lfn_s1) + call Read_restart_field (file_restart, 'lfn_s2', this%lfn_s2) + call Read_restart_field (file_restart, 'lfn_s3', this%lfn_s3) + call Read_restart_field (file_restart, 'lfn_out', this%lfn_out) + call Read_restart_field (file_restart, 'tign_g', this%tign_g) + call Get_netcdf_var_mpi (file_restart, this%cfbm_comm, this%nx, this%ny, & + this%ifps, this%ifpe, this%jfps, this%jfpe, 'fuel_frac', & + this%fuel_frac(this%ifps:this%ifpe, this%jfps:this%jfpe)) + call Read_restart_field (file_restart, 'fire_area', this%fire_area) + call Get_netcdf_var_mpi (file_restart, this%cfbm_comm, this%nx, this%ny, & + this%ifps, this%ifpe, this%jfps, this%jfpe, 'fuel_frac_burnt_dt', & + this%fuel_frac_burnt_dt(this%ifps:this%ifpe, this%jfps:this%jfpe)) + call Read_restart_field (file_restart, 'fgrnhfx', this%fgrnhfx) + call Read_restart_field (file_restart, 'fgrnqfx', this%fgrnqfx) + call Read_restart_field (file_restart, 'fcanhfx', this%fcanhfx) + call Read_restart_field (file_restart, 'fcanqfx', this%fcanqfx) + call Read_restart_field (file_restart, 'flame_length', this%flame_length) + call Read_restart_field (file_restart, 'ros', this%ros) + call Read_restart_field (file_restart, 'ros_front', this%ros_front) + call Read_restart_field (file_restart, 'emis_smoke', this%emis_smoke) + call Read_restart_field (file_restart, 'fmc_g', this%fmc_g) + call Read_restart_field (file_restart, 'fuel_load_g', this%fuel_load_g) + call Read_restart_field (file_restart, 'fuel_time', this%fuel_time) + call Read_restart_field (file_restart, 'zsf', this%zsf) + call Read_restart_field (file_restart, 'dzdxf', this%dzdxf) + call Read_restart_field (file_restart, 'dzdyf', this%dzdyf) + call Read_restart_field (file_restart, 'nfuel_cat', this%nfuel_cat) + call Read_restart_field (file_restart, 'uf', this%uf) + call Read_restart_field (file_restart, 'vf', this%vf) + call Read_restart_field (file_restart, 'fz0', this%fz0) + + if (config_flags%ideal_opt == 0) then + call Read_restart_field (file_restart, 'fire_t2', this%fire_t2) + call Read_restart_field (file_restart, 'fire_q2', this%fire_q2) + call Read_restart_field (file_restart, 'fire_psfc', this%fire_psfc) + call Read_restart_field (file_restart, 'fire_rain', this%fire_rain) + if (config_flags%fmoist_run) then + call Read_restart_field (file_restart, 'fire_t2_old', this%fire_t2_old) + call Read_restart_field (file_restart, 'fire_q2_old', this%fire_q2_old) + call Read_restart_field (file_restart, 'fire_psfc_old', this%fire_psfc_old) + call Read_restart_field (file_restart, 'fire_rain_old', this%fire_rain_old) + end if + end if + if (config_flags%fmoist_run) call Read_restart_fmc (file_restart) + + call Exchange_restart_halos () + + if (allocated (this%ros_param) .and. allocated (this%fuels)) then + do ij = 1, this%num_tiles + 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 + end if + + if (DEBUG_LOCAL) call Print_message ('Leaving Read_restart...') + + contains + + subroutine Exchange_restart_halos () + + implicit none + + +#ifdef DM_PARALLEL + call Exchange_restart_field (this%lfn) + call Exchange_restart_field (this%lfn_hist) + call Exchange_restart_field (this%lfn_0) + call Exchange_restart_field (this%lfn_1) + call Exchange_restart_field (this%lfn_2) + call Exchange_restart_field (this%lfn_s0) + call Exchange_restart_field (this%lfn_s1) + call Exchange_restart_field (this%lfn_s2) + call Exchange_restart_field (this%lfn_s3) + call Exchange_restart_field (this%lfn_out) + call Exchange_restart_field (this%tign_g) + call Do_halo_exchange_with_corners (this%fuel_frac, & + this%ifms, this%ifme, this%jfms, this%jfme, & + this%ifps, this%ifpe, this%jfps, this%jfpe, N_POINTS_IN_HALO, this%cart_comm) + call Exchange_restart_field (this%fire_area) + call Do_halo_exchange_with_corners (this%fuel_frac_burnt_dt, & + this%ifms, this%ifme, this%jfms, this%jfme, & + this%ifps, this%ifpe, this%jfps, this%jfpe, N_POINTS_IN_HALO, this%cart_comm) + call Exchange_restart_field (this%fgrnhfx) + call Exchange_restart_field (this%fgrnqfx) + call Exchange_restart_field (this%fcanhfx) + call Exchange_restart_field (this%fcanqfx) + call Exchange_restart_field (this%flame_length) + call Exchange_restart_field (this%ros) + call Exchange_restart_field (this%ros_front) + call Exchange_restart_field (this%emis_smoke) + call Exchange_restart_field (this%fmc_g) + call Exchange_restart_field (this%fuel_load_g) + call Exchange_restart_field (this%fuel_time) + call Exchange_restart_field (this%zsf) + call Exchange_restart_field (this%dzdxf) + call Exchange_restart_field (this%dzdyf) + call Exchange_restart_field (this%nfuel_cat) + call Exchange_restart_field (this%uf) + call Exchange_restart_field (this%vf) + call Exchange_restart_field (this%fz0) + + if (config_flags%ideal_opt == 0) then + call Exchange_restart_field (this%fire_t2) + call Exchange_restart_field (this%fire_q2) + call Exchange_restart_field (this%fire_psfc) + call Exchange_restart_field (this%fire_rain) + if (config_flags%fmoist_run) then + call Exchange_restart_field (this%fire_t2_old) + call Exchange_restart_field (this%fire_q2_old) + call Exchange_restart_field (this%fire_psfc_old) + call Exchange_restart_field (this%fire_rain_old) + end if + end if + + if (config_flags%fmoist_run) call Exchange_restart_fmc_halos () +#endif + + end subroutine Exchange_restart_halos + + subroutine Exchange_restart_field (var) + + implicit none + + real, dimension(this%ifms:this%ifme, this%jfms:this%jfme), intent (in out) :: var + + +#ifdef DM_PARALLEL + call Do_halo_exchange_with_corners (var, this%ifms, this%ifme, this%jfms, this%jfme, & + this%ifps, this%ifpe, this%jfps, this%jfpe, N_POINTS_IN_HALO, this%cart_comm) +#endif + + end subroutine Exchange_restart_field + + subroutine Exchange_restart_fmc_halos () + + implicit none + + integer :: k, n_moisture_classes + + +#ifdef DM_PARALLEL + if (.not. allocated (this%fmc_param)) call Stop_simulation ('FMC restart halo exchange requires allocated fmc_param') + + select type (fmc_param => this%fmc_param) + type is (fmc_wrffire_t) + if (.not. allocated (fmc_param%fmc_gc)) call Stop_simulation ('FMC restart halo exchange requires allocated fmc_gc') + + n_moisture_classes = size (fmc_param%fmc_gc, 2) + do k = 1, n_moisture_classes + call Do_halo_exchange_with_corners (fmc_param%fmc_gc(this%ifms:this%ifme, k, this%jfms:this%jfme), & + this%ifms, this%ifme, this%jfms, this%jfme, & + this%ifps, this%ifpe, this%jfps, this%jfpe, N_POINTS_IN_HALO, this%cart_comm) + end do + + class default + call Stop_simulation ('Restart halo exchange is not implemented for selected FMC component') + end select +#endif + + end subroutine Exchange_restart_fmc_halos + + subroutine Read_restart_field (file_name, var_name, var) + + implicit none + + character (len = *), intent (in) :: file_name, var_name + real, dimension(this%ifms:this%ifme, this%jfms:this%jfme), intent (in out) :: var + + + call Get_netcdf_var_mpi (file_name, this%cfbm_comm, this%nx, this%ny, & + this%ifps, this%ifpe, this%jfps, this%jfpe, var_name, var(this%ifps:this%ifpe, this%jfps:this%jfpe)) + + end subroutine Read_restart_field + + subroutine Read_restart_fmc (file_name) + + implicit none + + character (len = *), intent (in) :: file_name + + real (kind = REAL32) :: att_real32 + real :: fmoist_lasttime, fmoist_nexttime + integer :: n_moisture_classes + + + if (.not. allocated (this%fmc_param)) call Stop_simulation ('FMC restart read requires allocated fmc_param') + + select type (fmc_param => this%fmc_param) + type is (fmc_wrffire_t) + if (.not. allocated (fmc_param%fmc_gc)) call Stop_simulation ('FMC restart read requires allocated fmc_gc') + + n_moisture_classes = size (fmc_param%fmc_gc, 2) + call Get_netcdf_var_mpi (file_name, this%cfbm_comm, this%nx, this%ny, n_moisture_classes, & + 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)) + + fmoist_lasttime = 0.0 + fmoist_nexttime = 0.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 + end if + call Broadcast_restart_real (this, fmoist_lasttime) + call Broadcast_restart_real (this, fmoist_nexttime) + fmc_param%fmoist_lasttime = fmoist_lasttime + fmc_param%fmoist_nexttime = fmoist_nexttime + + class default + call Stop_simulation ('Restart read is not implemented for selected FMC component') + end select + + end subroutine Read_restart_fmc + + subroutine Set_next_datetime_after_restart (datetime_next, interval_seconds) + + implicit none + + type (datetime_t), intent (in out) :: datetime_next + integer, intent (in) :: interval_seconds + + datetime_next = this%datetime_now + if (interval_seconds > 0) then + datetime_next = this%datetime_start + do while (datetime_next <= this%datetime_now) + call datetime_next%Add_seconds (interval_seconds) + end do + end if + + end subroutine Set_next_datetime_after_restart + + end subroutine Read_restart + + subroutine Read_restart_field_2d (file_name, var_name, nx, ny, ifms, ifme, jfms, jfme, ifps, ifpe, jfps, jfpe, cfbm_comm, var) + + implicit none + + character (len = *), intent (in) :: file_name, var_name + integer, intent (in) :: nx, ny, ifms, ifme, jfms, jfme, ifps, ifpe, jfps, jfpe, cfbm_comm + real, dimension(ifms:ifme, jfms:jfme), intent (in out) :: var + + + call Get_netcdf_var_mpi (file_name, cfbm_comm, nx, ny, ifps, ifpe, jfps, jfpe, & + var_name, var(ifps:ifpe, jfps:jfpe)) + + end subroutine Read_restart_field_2d + + subroutine Validate_restart_integer (file_name, att_name, expected_value) + + implicit none + + character (len = *), intent (in) :: file_name, att_name + integer, intent (in) :: expected_value + + integer (kind = INT32) :: att_value + character (len = :), allocatable :: msg + + + call Get_netcdf_att (file_name, 'global', att_name, att_value) + if (att_value /= expected_value) then + msg = 'Restart metadata mismatch: '//trim (att_name) + call Stop_simulation (msg) + end if + + end subroutine Validate_restart_integer + + subroutine Validate_restart_real (file_name, att_name, expected_value) + + implicit none + + character (len = *), intent (in) :: file_name, att_name + real, intent (in) :: expected_value + + real (kind = REAL32) :: att_value, tolerance + character (len = :), allocatable :: msg + + + call Get_netcdf_att (file_name, 'global', att_name, att_value) + tolerance = max (1.0e-6_REAL32, abs (real (expected_value, kind = REAL32)) * 1.0e-6_REAL32) + if (abs (att_value - real (expected_value, kind = REAL32)) > tolerance) then + msg = 'Restart metadata mismatch: '//trim (att_name) + call Stop_simulation (msg) + end if + + end subroutine Validate_restart_real + + module subroutine Write_restart (this, config_flags) + + implicit none + + class (state_fire_t), intent (in) :: this + type (namelist_t), intent (in) :: config_flags + + character (len = :), allocatable :: file_restart + integer :: start_year, start_month, start_day, start_hour, start_minute, start_second, & + restart_year, restart_month, restart_day, restart_hour, restart_minute, restart_second + integer :: restart_sr_x, restart_sr_y, rank, ierr + logical, parameter :: DEBUG_LOCAL = .false. + + + if (DEBUG_LOCAL) call Print_message ('Entering Write_restart...') + +#ifdef DM_PARALLEL + if (.not. this%is_cfbm_comm_set) call Stop_simulation ('The MPI CFBM communicator has not been set') + call Mpi_comm_rank (this%cfbm_comm, rank, ierr) + if (ierr /= MPI_SUCCESS) call Stop_simulation ('Problems with Mpi_comm_rank ') +#else + rank = 0 +#endif + + if (config_flags%ideal_opt /= 0 .and. config_flags%ideal_opt /= 1) & + call Stop_simulation ('Write_restart is implemented for standalone idealized and real runs only') + + file_restart = Build_restart_file_name (this%datetime_now%datetime) + + call this%datetime_start%Get_datetime_as_ints (start_year, start_month, start_day, start_hour, start_minute, start_second) + call this%datetime_now%Get_datetime_as_ints (restart_year, restart_month, restart_day, restart_hour, restart_minute, restart_second) + restart_sr_x = nint (this%proj%dx / this%dx) + restart_sr_y = nint (this%proj%dy / this%dy) + if (restart_sr_x <= 0 .or. restart_sr_y <= 0) call Stop_simulation ('Restart subgrid ratios must be positive') + + if (rank == 0) then + call Create_netcdf_file (file_name = file_restart) + call Add_netcdf_dim (file_restart, NAME_DIM_X, this%nx) + call Add_netcdf_dim (file_restart, NAME_DIM_Y, this%ny) + + call Add_netcdf_att (file_restart, 'global', 'start_year', int (start_year, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'start_month', int (start_month, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'start_day', int (start_day, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'start_hour', int (start_hour, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'start_minute', int (start_minute, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'start_second', int (start_second, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'restart_year', int (restart_year, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'restart_month', int (restart_month, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'restart_day', int (restart_day, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'restart_hour', int (restart_hour, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'restart_minute', int (restart_minute, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'restart_second', int (restart_second, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'itimestep', int (this%itimestep, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'fire_perimeter_ignited', & + int (merge (1, 0, this%fire_perimeter_ignited), kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'nx', int (this%nx, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'ny', int (this%ny, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'dt', real (this%dt, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'dx', real (this%dx, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'dy', real (this%dy, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'cen_lat', real (this%cen_lat, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'cen_lon', real (this%cen_lon, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'map_proj', int (1, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'sr_x', int (restart_sr_x, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'sr_y', int (restart_sr_y, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'stand_lon', real (this%proj%standard_lon, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'true_lat_1', real (this%proj%true_lat_1, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'true_lat_2', real (this%proj%true_lat_2, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'ideal_opt', int (config_flags%ideal_opt, kind = INT32)) + 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', '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)) + call Add_netcdf_att (file_restart, 'global', 'fire_lsm_reinit_iter', int (config_flags%fire_lsm_reinit_iter, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'fire_viscosity', real (config_flags%fire_viscosity, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'fire_viscosity_bg', real (config_flags%fire_viscosity_bg, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'fire_viscosity_ngp', int (config_flags%fire_viscosity_ngp, kind = INT32)) + call Add_netcdf_att (file_restart, 'global', 'fire_viscosity_band', real (config_flags%fire_viscosity_band, kind = REAL32)) + call Add_netcdf_att (file_restart, 'global', 'reinit_pseudot_coef', real (config_flags%reinit_pseudot_coef, kind = REAL32)) + end if + + call Restart_io_barrier (this) + + call Add_restart_field ('lfn', this%lfn) + call Add_restart_field ('lfn_hist', this%lfn_hist) + call Add_restart_field ('lfn_0', this%lfn_0) + call Add_restart_field ('lfn_1', this%lfn_1) + call Add_restart_field ('lfn_2', this%lfn_2) + call Add_restart_field ('lfn_s0', this%lfn_s0) + call Add_restart_field ('lfn_s1', this%lfn_s1) + call Add_restart_field ('lfn_s2', this%lfn_s2) + call Add_restart_field ('lfn_s3', this%lfn_s3) + call Add_restart_field ('lfn_out', this%lfn_out) + call Add_restart_field ('tign_g', this%tign_g) + call Add_netcdf_var_mpi (file_restart, this%cfbm_comm, this%nx, this%ny, & + this%ifps, this%ifpe, this%jfps, this%jfpe, 'fuel_frac', & + this%fuel_frac(this%ifps:this%ifpe, this%jfps:this%jfpe)) + call Add_restart_field ('fire_area', this%fire_area) + call Add_netcdf_var_mpi (file_restart, this%cfbm_comm, this%nx, this%ny, & + this%ifps, this%ifpe, this%jfps, this%jfpe, 'fuel_frac_burnt_dt', & + this%fuel_frac_burnt_dt(this%ifps:this%ifpe, this%jfps:this%jfpe)) + call Add_restart_field ('fgrnhfx', this%fgrnhfx) + call Add_restart_field ('fgrnqfx', this%fgrnqfx) + call Add_restart_field ('fcanhfx', this%fcanhfx) + call Add_restart_field ('fcanqfx', this%fcanqfx) + call Add_restart_field ('flame_length', this%flame_length) + call Add_restart_field ('ros', this%ros) + call Add_restart_field ('ros_front', this%ros_front) + call Add_restart_field ('emis_smoke', this%emis_smoke) + call Add_restart_field ('fmc_g', this%fmc_g) + call Add_restart_field ('fuel_load_g', this%fuel_load_g) + call Add_restart_field ('fuel_time', this%fuel_time) + call Add_restart_field ('zsf', this%zsf) + call Add_restart_field ('dzdxf', this%dzdxf) + call Add_restart_field ('dzdyf', this%dzdyf) + call Add_restart_field ('nfuel_cat', this%nfuel_cat) + call Add_restart_field ('uf', this%uf) + call Add_restart_field ('vf', this%vf) + call Add_restart_field ('fz0', this%fz0) + + if (config_flags%ideal_opt == 0) then + call Add_restart_field ('fire_t2', this%fire_t2) + call Add_restart_field ('fire_q2', this%fire_q2) + call Add_restart_field ('fire_psfc', this%fire_psfc) + call Add_restart_field ('fire_rain', this%fire_rain) + if (config_flags%fmoist_run) then + call Add_restart_field ('fire_t2_old', this%fire_t2_old) + call Add_restart_field ('fire_q2_old', this%fire_q2_old) + call Add_restart_field ('fire_psfc_old', this%fire_psfc_old) + call Add_restart_field ('fire_rain_old', this%fire_rain_old) + end if + end if + if (config_flags%fmoist_run) then + if (.not. allocated (this%fmc_param)) call Stop_simulation ('FMC restart write requires allocated fmc_param') + call this%fmc_param%Add_restart_vars (file_restart, this%cfbm_comm, this%nx, this%ny, & + this%ifps, this%ifpe, this%jfps, this%jfpe, rank) + end if + + call Restart_io_barrier (this) + + if (DEBUG_LOCAL) call Print_message ('Leaving Write_restart...') + + contains + + subroutine Add_restart_field (var_name, var) + + implicit none + + character (len = *), intent (in) :: var_name + real, dimension(this%ifms:this%ifme, this%jfms:this%jfme), intent (in) :: var + + + call Add_netcdf_var_mpi (file_restart, this%cfbm_comm, this%nx, this%ny, this%ifps, this%ifpe, this%jfps, this%jfpe, & + var_name, var(this%ifps:this%ifpe, this%jfps:this%jfpe)) + + end subroutine Add_restart_field + + end subroutine Write_restart + + end submodule state_restart_submod