subroutine dyn_mpas_init_phase4(self, coupling_time_interval)
! Module(s) from MPAS.
use atm_core, only: atm_mpas_init_block
use atm_time_integration, only: mpas_atm_dynamics_init
use mpas_atm_dimensions, only: mpas_atm_set_dims
use mpas_atm_halos, only: atm_build_halo_groups, exchange_halo_group
use mpas_atm_threading, only: mpas_atm_threading_init
use mpas_attlist, only: mpas_modify_att
use mpas_constants, only: mpas_constants_compute_derived
use mpas_derived_types, only: field0dreal, field2dreal, mpas_pool_type, mpas_time_type, &
mpas_start_time
use mpas_field_routines, only: mpas_allocate_scratch_field
use mpas_pool_routines, only: mpas_pool_get_field, mpas_pool_initialize_time_levels
use mpas_string_utils, only: mpas_string_replace
use mpas_timekeeping, only: mpas_get_clock_time, mpas_get_time
class(mpas_dynamical_core_type), intent(inout) :: self
integer, intent(in) :: coupling_time_interval ! Set the time interval, in seconds, over which MPAS dynamical core
! should integrate each time it is called to run.
character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase4'
character(strkind) :: date_time
character(strkind), pointer :: initial_time_1, initial_time_2
character(strkind), pointer :: xtime
integer :: ierr
integer, pointer :: nvertlevels, maxedges, maxedges2, num_scalars
logical, pointer :: config_do_restart
real(rkind), pointer :: config_dt
type(field0dreal), pointer :: field_0d_real
type(field2dreal), pointer :: field_2d_real
type(mpas_pool_type), pointer :: mpas_pool
type(mpas_time_type) :: mpas_time
call self % debug_print(log_level_debug, subname // ' entered')
nullify(initial_time_1, initial_time_2)
nullify(xtime)
nullify(nvertlevels, maxedges, maxedges2, num_scalars)
nullify(config_do_restart)
nullify(config_dt)
nullify(field_0d_real)
nullify(field_2d_real)
nullify(mpas_pool)
if (coupling_time_interval <= 0) then
call self % model_error('Invalid coupling time interval ' // stringify([real(coupling_time_interval, rkind)]), &
subname, __LINE__)
end if
call self % get_variable_pointer(config_dt, 'cfg', 'config_dt')
if (config_dt <= 0.0_rkind) then
call self % model_error('Invalid time step ' // stringify([config_dt]), &
subname, __LINE__)
end if
! `config_dt` in MPAS is a floating-point number. Testing floating-point numbers for divisibility is not trivial and
! should be done carefully.
if (.not. almost_divisible(real(coupling_time_interval, rkind), config_dt)) then
call self % model_error('Coupling time interval ' // stringify([real(coupling_time_interval, rkind)]) // &
' must be divisible by time step ' // stringify([config_dt]), subname, __LINE__)
end if
self % coupling_time_interval = coupling_time_interval
self % number_of_time_steps = 0
call self % debug_print(log_level_info, 'Coupling time interval is ' // &
stringify([real(self % coupling_time_interval, rkind)]) // ' seconds')
call self % debug_print(log_level_info, 'Time step is ' // &
stringify([config_dt]) // ' seconds')
nullify(config_dt)
! Compute derived constants.
call mpas_constants_compute_derived()
! Set up OpenMP threading.
call self % debug_print(log_level_info, 'Setting up OpenMP threading')
call mpas_atm_threading_init(self % domain_ptr % blocklist, ierr=ierr)
if (ierr /= 0) then
call self % model_error('OpenMP threading setup failed for core ' // trim(self % domain_ptr % core % corename), &
subname, __LINE__)
end if
! Set up inner dimensions used by arrays in optimized dynamics subroutines.
call self % debug_print(log_level_info, 'Setting up dimensions')
call self % get_variable_pointer(nvertlevels, 'dim', 'nVertLevels')
call self % get_variable_pointer(maxedges, 'dim', 'maxEdges')
call self % get_variable_pointer(maxedges2, 'dim', 'maxEdges2')
call self % get_variable_pointer(num_scalars, 'dim', 'num_scalars')
call mpas_atm_set_dims(nvertlevels, maxedges, maxedges2, num_scalars)
nullify(nvertlevels, maxedges, maxedges2, num_scalars)
! Build halo exchange groups and set the `exchange_halo_group` procedure pointer, which is used to
! exchange the halo layers of all fields in the named group.
nullify(exchange_halo_group)
call atm_build_halo_groups(self % domain_ptr, ierr=ierr)
if (ierr /= 0) then
call self % model_error('Failed to build halo exchange groups', subname, __LINE__)
end if
if (.not. associated(exchange_halo_group)) then
call self % model_error('Failed to build halo exchange groups', subname, __LINE__)
end if
! Variables in MPAS "state" pool have more than one time level. Copy the values from the first time level of
! such variables into all subsequent time levels to initialize them.
call self % get_variable_pointer(config_do_restart, 'cfg', 'config_do_restart')
if (.not. config_do_restart) then
! Run type is initial run.
call self % debug_print(log_level_info, 'Initializing time levels')
call self % get_pool_pointer(mpas_pool, 'state')
call mpas_pool_initialize_time_levels(mpas_pool)
nullify(mpas_pool)
end if
nullify(config_do_restart)
call exchange_halo_group(self % domain_ptr, 'initialization:u', ierr=ierr)
if (ierr /= 0) then
call self % model_error('Failed to exchange halo layers for group "initialization:u"', subname, __LINE__)
end if
! Initialize atmospheric variables (e.g., momentum, thermodynamic... variables in governing equations)
! as well as various aspects of time in MPAS.
call self % debug_print(log_level_info, 'Initializing atmospheric variables')
! Controlled by `config_start_time` in namelist.
mpas_time = mpas_get_clock_time(self % domain_ptr % clock, mpas_start_time, ierr=ierr)
if (ierr /= 0) then
call self % model_error('Failed to get time for "mpas_start_time"', subname, __LINE__)
end if
call mpas_get_time(mpas_time, datetimestring=date_time, ierr=ierr)
if (ierr /= 0) then
call self % model_error('Failed to get time for "mpas_start_time"', subname, __LINE__)
end if
! Controlled by `config_dt` in namelist.
call self % get_pool_pointer(mpas_pool, 'mesh')
call self % get_variable_pointer(config_dt, 'cfg', 'config_dt')
call atm_mpas_init_block(self % domain_ptr % dminfo, self % domain_ptr % streammanager, self % domain_ptr % blocklist, &
mpas_pool, config_dt)
nullify(mpas_pool)
nullify(config_dt)
call self % get_variable_pointer(xtime, 'state', 'xtime', time_level=1)
xtime = date_time
nullify(xtime)
! Initialize `initial_time` in the second time level. We need to do this manually because initial states
! are read into time level 1, and if we write anything from time level 2, `initial_time` will be invalid.
call self % get_variable_pointer(initial_time_1, 'state', 'initial_time', time_level=1)
call self % get_variable_pointer(initial_time_2, 'state', 'initial_time', time_level=2)
initial_time_2 = initial_time_1
! Set time units to CF-compliant "seconds since <date and time>".
call self % get_pool_pointer(mpas_pool, 'state')
call mpas_pool_get_field(mpas_pool, 'Time', field_0d_real, timelevel=1)
if (.not. associated(field_0d_real)) then
call self % model_error('Failed to find variable "Time"', subname, __LINE__)
end if
call mpas_modify_att(field_0d_real % attlists(1) % attlist, 'units', &
'seconds since ' // mpas_string_replace(initial_time_1, '_', ' '), ierr=ierr)
if (ierr /= 0) then
call self % model_error('Failed to set time units', subname, __LINE__)
end if
nullify(initial_time_1, initial_time_2)
nullify(mpas_pool)
nullify(field_0d_real)
call exchange_halo_group(self % domain_ptr, 'initialization:pv_edge,ru,rw', ierr=ierr)
if (ierr /= 0) then
call self % model_error('Failed to exchange halo layers for group "initialization:pv_edge,ru,rw"', subname, __LINE__)
end if
call self % debug_print(log_level_info, 'Initializing dynamics')
! Prepare dynamics for time integration.
call mpas_atm_dynamics_init(self % domain_ptr)
! Some additional "scratch" fields are needed for interoperability with CAM-SIMA, but they are not initialized by
! `mpas_atm_dynamics_init`. Initialize them below.
call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_uzonal', field_2d_real, timelevel=1)
call mpas_allocate_scratch_field(field_2d_real)
nullify(field_2d_real)
call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_umerid', field_2d_real, timelevel=1)
call mpas_allocate_scratch_field(field_2d_real)
nullify(field_2d_real)
call self % debug_print(log_level_debug, subname // ' completed')
call self % debug_print(log_level_info, 'Successful initialization of MPAS dynamical core')
contains
!> Test if `a` is divisible by `b`, where `a` and `b` are both reals.
!> (KCW, 2024-05-25)
pure function almost_divisible(a, b)
real(rkind), intent(in) :: a, b
logical :: almost_divisible
real(rkind) :: error_tolerance
error_tolerance = epsilon(1.0_rkind) * max(abs(a), abs(b))
if (almost_equal(mod(abs(a), abs(b)), 0.0_rkind, absolute_tolerance=error_tolerance) .or. &
almost_equal(mod(abs(a), abs(b)), abs(b), absolute_tolerance=error_tolerance)) then
almost_divisible = .true.
return
end if
almost_divisible = .false.
end function almost_divisible
!> Test `a` and `b` for approximate equality, where `a` and `b` are both reals.
!> (KCW, 2024-05-25)
pure function almost_equal(a, b, absolute_tolerance, relative_tolerance)
real(rkind), intent(in) :: a, b
real(rkind), optional, intent(in) :: absolute_tolerance, relative_tolerance
logical :: almost_equal
real(rkind) :: error_tolerance
if (present(relative_tolerance)) then
error_tolerance = relative_tolerance * max(abs(a), abs(b))
else
error_tolerance = epsilon(1.0_rkind) * max(abs(a), abs(b))
end if
if (present(absolute_tolerance)) then
error_tolerance = max(absolute_tolerance, error_tolerance)
end if
if (abs(a - b) <= error_tolerance) then
almost_equal = .true.
return
end if
almost_equal = .false.
end function almost_equal
end subroutine dyn_mpas_init_phase4