subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file)
! Module(s) from external libraries.
use pio, only: file_desc_t, pio_file_is_open
! Module(s) from MPAS.
use dyn_mpas_procedures, only: stringify
use mpas_bootstrapping, only: mpas_bootstrap_framework_phase1, mpas_bootstrap_framework_phase2
use mpas_derived_types, only: mpas_io_pnetcdf, mpas_pool_type
use mpas_pool_routines, only: mpas_pool_add_config, mpas_pool_add_dimension, mpas_pool_get_dimension
class(mpas_dynamical_core_type), intent(inout) :: self
integer, intent(in) :: number_of_constituents
type(file_desc_t), pointer, intent(in) :: pio_file
character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase3'
character(strkind) :: mesh_filename
character(strkind), pointer :: config_les_model
integer :: mesh_format
integer, pointer :: num_scalars
type(mpas_pool_type), pointer :: mpas_pool
call self % debug_print(log_level_debug, subname // ' entered')
nullify(config_les_model)
nullify(mpas_pool)
nullify(num_scalars)
call self % get_variable_pointer(config_les_model, 'cfg', 'config_les_model')
self % les_model = (trim(adjustl(config_les_model)) /= 'none')
nullify(config_les_model)
! In MPAS, there must be at least one constituent, `qv`, which denotes water vapor mixing ratio.
! Because MPAS has some hard-coded array accesses through the `index_qv` index, it will crash
! (i.e., segmentation fault due to invalid memory access) if `qv` is not allocated.
!
! Since version 8.4.0, if the large eddy simulation (LES) capability is enabled, two additional constituents,
! `qc` and `tke`, must also exist.
!
! Such capability is implemented provided that:
! 1. `qc` must be one of the constituents just like `qv`;
! 2. `tke` is treated as an extra "phantom" constituent that is internal to MPAS only.
if (.not. self % les_model) then
self % number_of_constituents = max(1, number_of_constituents) ! "1" is for `qv`.
else
self % number_of_constituents = max(2, number_of_constituents) ! "2" is for `qv` and `qc`.
end if
call self % debug_print(log_level_info, 'Number of constituents is ' // stringify([self % number_of_constituents]))
! Adding a config named "cam_pcnst" with the number of constituents will indicate to MPAS that
! it is operating as a dynamical core, and therefore it needs to allocate scalars separately
! from other Registry-defined fields. The special logic is located in `atm_setup_block`.
! This must be done before calling `mpas_bootstrap_framework_phase1`.
if (.not. self % les_model) then
! No need to add an extra "phantom" constituent, `tke`.
call mpas_pool_add_config(self % domain_ptr % configs, 'cam_pcnst', self % number_of_constituents)
else
! Need to add an extra "phantom" constituent, `tke`. There is additional logic
! for it in `dyn_mpas_define_scalar`.
call mpas_pool_add_config(self % domain_ptr % configs, 'cam_pcnst', self % number_of_constituents + 1)
end if
! Not actually used because a PIO file descriptor is directly supplied.
mesh_filename = 'external mesh'
mesh_format = mpas_io_pnetcdf
call self % debug_print(log_level_info, 'Checking PIO file descriptor')
if (.not. associated(pio_file)) then
call self % model_error('Invalid PIO file descriptor', subname, __LINE__)
end if
if (.not. pio_file_is_open(pio_file)) then
call self % model_error('Invalid PIO file descriptor', subname, __LINE__)
end if
call self % debug_print(log_level_info, 'Bootstrapping MPAS framework (Phase 1/2)')
! Finish setting up blocks.
call mpas_bootstrap_framework_phase1(self % domain_ptr, mesh_filename, mesh_format, pio_file_desc=pio_file)
call self % debug_print(log_level_info, 'Bootstrapping MPAS framework (Phase 2/2)')
! Finish setting up fields.
call mpas_bootstrap_framework_phase2(self % domain_ptr, pio_file_desc=pio_file)
! "num_scalars" is a dimension variable, but it only exists in MPAS "state" pool.
! Fix this inconsistency by also adding it to MPAS "dimension" pool.
call self % get_pool_pointer(mpas_pool, 'state')
call mpas_pool_get_dimension(mpas_pool, 'num_scalars', num_scalars)
if (.not. associated(num_scalars)) then
call self % model_error('Failed to find variable "num_scalars"', subname, __LINE__)
end if
! While we are at it, check if its value is consistent.
if (.not. self % les_model) then
if (num_scalars /= self % number_of_constituents) then
call self % model_error('Failed to allocate constituents', subname, __LINE__)
end if
else
if (num_scalars /= self % number_of_constituents + 1) then
call self % model_error('Failed to allocate constituents', subname, __LINE__)
end if
end if
call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'num_scalars', num_scalars)
nullify(mpas_pool)
nullify(num_scalars)
! At this point, what follows next depends on the specific use case. In no particular order:
! * Use `dyn_mpas_define_scalar` to define the names of constituents at run-time.
! * Use `dyn_mpas_read_write_stream` to read mesh variables.
! * Follow up with a call to `dyn_mpas_compute_unit_vector` immediately. This is by design.
! * For setting analytic initial condition, use `get_variable_pointer` to inject data directly into MPAS memory.
! * Use `dyn_mpas_compute_edge_wind` where appropriate.
! * Use `dyn_mpas_exchange_halo` where appropriate.
! * Use `dyn_mpas_read_write_stream` to read initial condition or restart.
! * Finally, use `dyn_mpas_init_phase4` to conclude the initialization of MPAS dynamical core.
call self % debug_print(log_level_debug, subname // ' completed')
end subroutine dyn_mpas_init_phase3