dyn_mpas_define_scalar Subroutine

    subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species)
        ! Module(s) from MPAS.
        use mpas_derived_types, only: field3dreal, mpas_pool_type
        use mpas_pool_routines, only: mpas_pool_add_dimension, mpas_pool_get_field

        class(mpas_dynamical_core_type), intent(inout) :: self
        character(*), intent(in) :: constituent_name(:)
        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' &
        ]

        character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_define_scalar'
        integer :: i, j, ierr
        integer :: index_qv, index_water_start, index_water_end
        integer :: time_level
        type(field3dreal), pointer :: field_3d_real
        type(mpas_pool_type), pointer :: mpas_pool

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

        nullify(field_3d_real)
        nullify(mpas_pool)

        if (self % number_of_constituents == 0) then
            call self % model_error('Constituents must be allocated before being defined', subname, __LINE__)
        end if

        ! Input sanitization.

        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 .and. self % number_of_constituents == 1) then
            ! If constituent definitions are empty, `qv` is the only constituent per MPAS requirements.
            ! See `dyn_mpas_init_phase3` for details.
            allocate(self % constituent_name(1), stat=ierr)

            if (ierr /= 0) then
                call self % model_error('Failed to allocate constituent_name', subname, __LINE__)
            end if

            allocate(self % is_water_species(1), stat=ierr)

            if (ierr /= 0) then
                call self % model_error('Failed to allocate is_water_species', subname, __LINE__)
            end if

            self % constituent_name(1) = mpas_scalar_qv_standard_name(1)
            self % is_water_species(1) = .true.
        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

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

            if (ierr /= 0) then
                call self % model_error('Failed to allocate constituent_name', subname, __LINE__)
            end if

            self % constituent_name(:) = adjustl(constituent_name)

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

            if (ierr /= 0) then
                call self % model_error('Failed to allocate is_water_species', subname, __LINE__)
            end if

            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

            ! `qv` must be present in constituents per MPAS requirements. It is a water species by definition.
            ! See `dyn_mpas_init_phase3` for details.
            index_qv = 0

            ! Lower index in `mpas_scalar_qv_standard_name` has higher precedence, with index 1 being exactly what MPAS wants.
            set_index_qv: do i = 1, size(mpas_scalar_qv_standard_name)
                do j = 1, self % number_of_constituents
                    if (self % constituent_name(j) == mpas_scalar_qv_standard_name(i) .and. self % is_water_species(j)) then
                        index_qv = j

                        ! The best candidate of `qv` has been found. Exit prematurely.
                        exit set_index_qv
                    end if
                end do
            end do set_index_qv

            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
        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), stat=ierr)

        if (ierr /= 0) then
            call self % model_error('Failed to allocate index_mpas_scalar_to_constituent', subname, __LINE__)
        end if

        self % index_mpas_scalar_to_constituent(:) = 0
        j = 1

        ! Place water species first per MPAS requirements.
        do i = 1, self % number_of_constituents
            if (self % is_water_species(i)) then
                self % index_mpas_scalar_to_constituent(j) = i
                j = j + 1
            end if
        end do

        index_water_start = 1
        index_water_end = count(self % is_water_species)

        ! Place non-water species second per MPAS requirements.
        do i = 1, self % number_of_constituents
            if (.not. self % is_water_species(i)) then
                self % index_mpas_scalar_to_constituent(j) = i
                j = j + 1
            end if
        end do

        ! 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), stat=ierr)

        if (ierr /= 0) then
            call self % model_error('Failed to allocate index_constituent_to_mpas_scalar', 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)

        ! 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)))
            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

        ! Define "scalars" for MPAS.

        call self % debug_print(log_level_info, 'Defining MPAS scalars')

        call self % get_pool_pointer(mpas_pool, 'state')

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

        ! MPAS "state" pool has two time levels.
        time_level = 2

        do i = 1, time_level
            call mpas_pool_get_field(mpas_pool, 'scalars', field_3d_real, timelevel=i)

            if (.not. associated(field_3d_real)) then
                call self % model_error('Failed to find variable "scalars"', subname, __LINE__)
            end if

            do j = 1, self % number_of_constituents
                field_3d_real % constituentnames(j) = &
                    trim(adjustl(self % constituent_name(self % index_mpas_scalar_to_constituent(j))))

                ! Print information about MPAS scalars. Only do it once.
                if (i == 1) then
                    call self % debug_print(log_level_verbose, 'MPAS scalar index ' // stringify([j]))
                    call self % debug_print(log_level_verbose, '    MPAS scalar name: ' // &
                        trim(field_3d_real % constituentnames(j)))
                    call self % debug_print(log_level_verbose, '    Is water species: ' // &
                        stringify([self % is_water_species(self % index_mpas_scalar_to_constituent(j))]))
                    call self % debug_print(log_level_verbose, '    Index mapping from MPAS scalar to constituent: ' // &
                        stringify([j]) // ' -> ' // stringify([self % index_mpas_scalar_to_constituent(j)]))
                end if
            end do

            nullify(field_3d_real)
        end do

        nullify(mpas_pool)

        ! Define "scalars_tend" for MPAS.

        call self % debug_print(log_level_info, 'Defining MPAS scalar tendencies')

        call self % get_pool_pointer(mpas_pool, 'tend')

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

        ! MPAS "tend" pool only has one time level.
        time_level = 1

        do i = 1, time_level
            call mpas_pool_get_field(mpas_pool, 'scalars_tend', field_3d_real, timelevel=i)

            if (.not. associated(field_3d_real)) then
                call self % model_error('Failed to find variable "scalars_tend"', subname, __LINE__)
            end if

            do j = 1, self % number_of_constituents
                field_3d_real % constituentnames(j) = &
                    'tendency_of_' // trim(adjustl(self % constituent_name(self % index_mpas_scalar_to_constituent(j))))

                ! Print information about MPAS scalar tendencies. Only do it once.
                if (i == 1) then
                    call self % debug_print(log_level_verbose, 'MPAS scalar tendency index ' // stringify([j]))
                    call self % debug_print(log_level_verbose, '    MPAS scalar tendency name: ' // &
                        trim(field_3d_real % constituentnames(j)))
                    call self % debug_print(log_level_verbose, '    Is water species: ' // &
                        stringify([self % is_water_species(self % index_mpas_scalar_to_constituent(j))]))
                    call self % debug_print(log_level_verbose, '    Index mapping from MPAS scalar tendency to constituent: ' // &
                        stringify([j]) // ' -> ' // stringify([self % index_mpas_scalar_to_constituent(j)]))
                end if
            end do

            nullify(field_3d_real)
        end do

        nullify(mpas_pool)

        ! For consistency, also add dimension variables to MPAS "dimension" pool.

        call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'index_qv', index_qv)
        call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'moist_start', index_water_start)
        call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'moist_end', index_water_end)

        call self % debug_print(log_level_debug, 'index_qv = ' // stringify([index_qv]))
        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 % debug_print(log_level_debug, subname // ' completed')
    end subroutine dyn_mpas_define_scalar