dyn_mpas_init_phase3 Subroutine

    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 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
        integer :: mesh_format
        integer, pointer :: num_scalars
        type(mpas_pool_type), pointer :: mpas_pool

        call self % debug_print(log_level_debug, subname // ' entered')

        nullify(mpas_pool)
        nullify(num_scalars)

        ! 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.
        self % number_of_constituents = max(1, number_of_constituents)

        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`.
        call mpas_pool_add_config(self % domain_ptr % configs, 'cam_pcnst', self % number_of_constituents)

        ! 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 (num_scalars /= self % number_of_constituents) then
            call self % model_error('Failed to allocate constituents', subname, __LINE__)
        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