dyn_mpas_define_scalar Subroutine

    subroutine dyn_mpas_define_scalar(self, constituent_name, input_alias, is_water_species)
        ! Module(s) from MPAS.
        use dyn_mpas_procedures, only: index_unique, stringify
        use mpas_derived_types, only: mpas_pool_type
        use mpas_pool_routines, only: mpas_pool_add_dimension

        class(mpas_dynamical_core_type), intent(inout) :: self
        character(*), intent(in) :: constituent_name(:)
        character(*), intent(in) :: input_alias(:)
        logical, intent(in) :: is_water_species(:)

        !> Possible CCPP standard names of `qv`, which denotes water vapor mixing ratio.
        !> They are hard-coded here because MPAS needs to know where `qv` is.
        !> Index 1 is exactly what MPAS wants. Others also work, but need to be converted.
        character(*), parameter :: mpas_scalar_qv_standard_name(*) = [character(strkind) :: &
            'water_vapor_mixing_ratio_wrt_dry_air', &
            'water_vapor_mixing_ratio_wrt_moist_air', &
            'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water' &
        ]

        !> Possible CCPP standard names of `qc`, which denotes cloud liquid water mixing ratio.
        !> They are hard-coded here because MPAS needs to know where `qc` is.
        !> Index 1 is exactly what MPAS wants. Others also work, but need to be converted.
        character(*), parameter :: mpas_scalar_qc_standard_name(*) = [character(strkind) :: &
            'cloud_liquid_water_mixing_ratio_wrt_dry_air', &
            'cloud_liquid_water_mixing_ratio_wrt_moist_air', &
            'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' &
        ]

        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_define_scalar'
        character(strkind) :: cerr
        integer :: i
        integer :: ierr
        integer :: index_qv, index_qc, index_tke, index_water_start, index_water_end
        type(mpas_pool_type), pointer :: mpas_pool

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

        nullify(mpas_pool)

        ! `dyn_mpas_init_phase3` must be called before this subroutine.
        if (self % number_of_constituents == 0) then
            call self % model_error('Constituents must be allocated before being defined', subname, __LINE__)
        end if

        allocate(self % constituent_name(self % number_of_constituents), errmsg=cerr, stat=ierr)

        if (ierr /= 0) then
            call self % model_error('Failed to allocate constituent_name' // new_line('') // &
                'Allocation returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
                subname, __LINE__)
        end if

        allocate(self % input_alias(self % number_of_constituents), errmsg=cerr, stat=ierr)

        if (ierr /= 0) then
            call self % model_error('Failed to allocate input_alias' // new_line('') // &
                'Allocation returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
                subname, __LINE__)
        end if

        allocate(self % is_water_species(self % number_of_constituents), errmsg=cerr, stat=ierr)

        if (ierr /= 0) then
            call self % model_error('Failed to allocate is_water_species' // new_line('') // &
                'Allocation returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
                subname, __LINE__)
        end if

        ! Input sanitization.

        if (size(constituent_name) /= size(input_alias)) then
            call self % model_error('Mismatch between numbers of constituent names and their input aliases', subname, __LINE__)
        end if

        if (size(constituent_name) /= size(is_water_species)) then
            call self % model_error('Mismatch between numbers of constituent names and their waterness', subname, __LINE__)
        end if

        if (size(constituent_name) == 0) then
            if (.not. self % les_model) then
                ! If constituent definitions are empty and LES is disabled, `qv` is the only
                ! one constituent as per MPAS requirements.
                ! See `dyn_mpas_init_phase3` for details.
                if (self % number_of_constituents /= 1) then
                    call self % model_error( &
                        'Number of constituents must be 1 if no definitions are provided', subname, __LINE__)
                end if

                self % constituent_name(1) = mpas_scalar_qv_standard_name(1)
                self % input_alias(1) = ''
                self % is_water_species(1) = .true.
            else
                ! If constituent definitions are empty and LES is enabled, `qv`, `qc`, and `tke` are the only
                ! three constituents as per MPAS requirements. However, `tke` is excepted because it is
                ! internal to MPAS only.
                ! See `dyn_mpas_init_phase3` for details.
                if (self % number_of_constituents /= 2) then
                    call self % model_error( &
                        'Number of constituents must be 2 if no definitions are provided', subname, __LINE__)
                end if

                self % constituent_name(1) = mpas_scalar_qv_standard_name(1)
                self % input_alias(1) = ''
                self % is_water_species(1) = .true.

                self % constituent_name(2) = mpas_scalar_qc_standard_name(1)
                self % input_alias(2) = ''
                self % is_water_species(2) = .true.
            end if
        else
            if (size(constituent_name) /= self % number_of_constituents) then
                call self % model_error('Mismatch between numbers of constituents and their names', subname, __LINE__)
            end if

            if (any(len_trim(adjustl(constituent_name)) > len(self % constituent_name))) then
                call self % model_error('Constituent names are too long', subname, __LINE__)
            end if

            if (any(len_trim(adjustl(input_alias)) > len(self % input_alias))) then
                call self % model_error('Input aliases are too long', subname, __LINE__)
            end if

            self % constituent_name(:) = adjustl(constituent_name)
            self % input_alias(:) = adjustl(input_alias)
            self % is_water_species(:) = is_water_species(:)

            if (size(self % constituent_name) /= size(index_unique(self % constituent_name))) then
                call self % model_error('Constituent names must be unique', subname, __LINE__)
            end if

            if (count(self % input_alias /= '') /= size(index_unique(pack(self % input_alias, self % input_alias /= '')))) then
                call self % model_error('Non-empty input aliases must be unique', subname, __LINE__)
            end if
        end if

        ! Lower index in `mpas_scalar_qv_standard_name` has higher precedence, with index 1 being exactly what MPAS wants.
        do i = 1, size(mpas_scalar_qv_standard_name)
            index_qv = findloc(self % constituent_name, mpas_scalar_qv_standard_name(i), dim=1, mask=self % is_water_species)

            ! The best candidate of `qv` has been found. Exit prematurely.
            if (index_qv > 0) then
                exit
            end if
        end do

        ! `qv` must be present in constituents as per MPAS requirements. It is a water species by definition.
        ! See `dyn_mpas_init_phase3` for details.
        if (index_qv == 0) then
            call self % model_error('Constituent names must contain one of: ' // &
                stringify(mpas_scalar_qv_standard_name) // ', and it must be a water species', subname, __LINE__)
        end if

        ! Lower index in `mpas_scalar_qc_standard_name` has higher precedence, with index 1 being exactly what MPAS wants.
        do i = 1, size(mpas_scalar_qc_standard_name)
            index_qc = findloc(self % constituent_name, mpas_scalar_qc_standard_name(i), dim=1, mask=self % is_water_species)

            ! The best candidate of `qc` has been found. Exit prematurely.
            if (index_qc > 0) then
                exit
            end if
        end do

        index_tke = 0

        ! If LES is enabled, `qc` and `tke` must be present in constituents as per MPAS requirements.
        ! Otherwise, it is fine to not have them.
        ! See `dyn_mpas_init_phase3` for details.
        if (self % les_model) then
            if (index_qc == 0) then
                call self % model_error('Constituent names must contain one of: ' // &
                    stringify(mpas_scalar_qc_standard_name) // ', and it must be a water species', subname, __LINE__)
            end if

            index_tke = self % number_of_constituents + 1
        end if

        ! Create index mapping between MPAS scalars and constituent names. For example,
        ! MPAS scalar index `i` corresponds to constituent index `index_mpas_scalar_to_constituent(i)`.

        call self % debug_print(log_level_info, 'Creating index mapping between MPAS scalars and CAM-SIMA constituents')

        allocate(self % index_mpas_scalar_to_constituent(self % number_of_constituents), errmsg=cerr, stat=ierr)

        if (ierr /= 0) then
            call self % model_error('Failed to allocate index_mpas_scalar_to_constituent' // new_line('') // &
                'Allocation returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
                subname, __LINE__)
        end if

        self % index_mpas_scalar_to_constituent(:) = [(i, i = 1, self % number_of_constituents)]

        ! Place water species first, then non-water species last, as per MPAS requirements.
        self % index_mpas_scalar_to_constituent(:) = [ &
            pack(self % index_mpas_scalar_to_constituent, self % is_water_species), &
            pack(self % index_mpas_scalar_to_constituent, .not. self % is_water_species) &
        ]

        index_water_start = 1
        index_water_end = count(self % is_water_species)

        ! Create inverse index mapping between MPAS scalars and constituent names. For example,
        ! Constituent index `i` corresponds to MPAS scalar index `index_constituent_to_mpas_scalar(i)`.

        call self % debug_print(log_level_info, 'Creating inverse index mapping between MPAS scalars and CAM-SIMA constituents')

        allocate(self % index_constituent_to_mpas_scalar(self % number_of_constituents), errmsg=cerr, stat=ierr)

        if (ierr /= 0) then
            call self % model_error('Failed to allocate index_constituent_to_mpas_scalar' // new_line('') // &
                'Allocation returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
                subname, __LINE__)
        end if

        self % index_constituent_to_mpas_scalar(:) = 0

        do i = 1, self % number_of_constituents
            self % index_constituent_to_mpas_scalar(self % index_mpas_scalar_to_constituent(i)) = i
        end do

        ! Set the index of `qv` in terms of MPAS scalars.
        index_qv = self % index_constituent_to_mpas_scalar(index_qv)

        ! Set the index of `qc`, if present, in terms of MPAS scalars.
        if (index_qc > 0) then
            index_qc = self % index_constituent_to_mpas_scalar(index_qc)
        end if

        ! Print information about constituents.
        do i = 1, self % number_of_constituents
            call self % debug_print(log_level_verbose, 'Constituent index ' // stringify([i]))
            call self % debug_print(log_level_verbose, '    Constituent name: ' // &
                trim(self % constituent_name(i)))
            if (self % input_alias(i) /= '') then
                call self % debug_print(log_level_verbose, '    Input alias: ' // &
                    trim(self % input_alias(i)))
            end if
            call self % debug_print(log_level_verbose, '    Is water species: ' // &
                stringify([self % is_water_species(i)]))
            call self % debug_print(log_level_verbose, '    Index mapping from constituent to MPAS scalar: ' // &
                stringify([i]) // ' -> ' // stringify([self % index_constituent_to_mpas_scalar(i)]))
        end do

        ! For consistency, add scalar dimensions to MPAS "dimension" pool.

        call self % get_pool_pointer(mpas_pool, 'dim')

        call mpas_pool_add_dimension(mpas_pool, 'index_qv', index_qv)
        call mpas_pool_add_dimension(mpas_pool, 'index_qc', index_qc)
        call mpas_pool_add_dimension(mpas_pool, 'index_tke', index_tke)
        call mpas_pool_add_dimension(mpas_pool, 'moist_start', index_water_start)
        call mpas_pool_add_dimension(mpas_pool, 'moist_end', index_water_end)

        nullify(mpas_pool)

        call self % debug_print(log_level_debug, 'index_qv = ' // stringify([index_qv]))
        call self % debug_print(log_level_debug, 'index_qc = ' // stringify([index_qc]))
        call self % debug_print(log_level_debug, 'index_tke = ' // stringify([index_tke]))
        call self % debug_print(log_level_debug, 'moist_start = ' // stringify([index_water_start]))
        call self % debug_print(log_level_debug, 'moist_end = ' // stringify([index_water_end]))

        call self % define_scalar_field()

        call self % debug_print(log_level_debug, subname // ' completed')
    end subroutine dyn_mpas_define_scalar