dynamics_to_physics_coupling Subroutine

    subroutine dynamics_to_physics_coupling()
        ! Module(s) from CAM-SIMA.
        use cam_logfile, only: debugout_debug, debugout_info
        use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_states, ncells_solve
        ! Module(s) from CESM Share.
        use shr_kind_mod, only: kind_r8 => shr_kind_r8
        ! Module(s) from MPAS.
        use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind

        character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling'
        integer :: column_index
        integer, allocatable :: is_water_species_index(:)
        integer, pointer :: index_qv
        ! Variable name suffixes have the following meanings:
        ! `*_col`: Variable is of each column.
        ! `*_int`: Variable is at layer interfaces.
        ! `*_mid`: Variable is at layer midpoints.
        real(kind_r8), allocatable :: pd_int_col(:), &       ! Dry hydrostatic air pressure (Pa).
                                      pd_mid_col(:), &       ! Dry hydrostatic air pressure (Pa).
                                      p_int_col(:), &        ! Full hydrostatic air pressure (Pa).
                                      p_mid_col(:), &        ! Full hydrostatic air pressure (Pa).
                                      z_int_col(:)           ! Geometric height (m).
        real(kind_r8), allocatable :: dpd_col(:), &          ! Dry air pressure difference (Pa) between layer interfaces.
                                      dp_col(:), &           ! Full air pressure difference (Pa) between layer interfaces.
                                      dz_col(:)              ! Geometric height difference (m) between layer interfaces.
        real(kind_r8), allocatable :: qv_mid_col(:), &       ! Water vapor mixing ratio (kg kg-1).
                                      sigma_all_q_mid_col(:) ! Summation of all water species mixing ratios (kg kg-1).
        real(kind_r8), allocatable :: rhod_mid_col(:), &     ! Dry air density (kg m-3).
                                      rho_mid_col(:)         ! Full air density (kg m-3).
        real(kind_r8), allocatable :: t_mid_col(:), &        ! Temperature (K).
                                      tm_mid_col(:), &       ! Modified "moist" temperature (K).
                                      tv_mid_col(:)          ! Virtual temperature (K).
        real(kind_r8), allocatable :: u_mid_col(:), &        ! Eastward wind velocity (m s-1).
                                      v_mid_col(:), &        ! Northward wind velocity (m s-1).
                                      omega_mid_col(:)       ! Vertical wind velocity (Pa s-1).
        real(kind_dyn_mpas), pointer :: exner(:, :)
        real(kind_dyn_mpas), pointer :: rho_zz(:, :)
        real(kind_dyn_mpas), pointer :: scalars(:, :, :)
        real(kind_dyn_mpas), pointer :: theta_m(:, :)
        real(kind_dyn_mpas), pointer :: ucellzonal(:, :), ucellmeridional(:, :), w(:, :)
        real(kind_dyn_mpas), pointer :: zgrid(:, :)
        real(kind_dyn_mpas), pointer :: zz(:, :)

        call dyn_debug_print(debugout_debug, subname // ' entered')

        call init_shared_variables()

        call dyn_exchange_constituent_states(direction='i', exchange=.true., conversion=.false.)

        call dyn_debug_print(debugout_info, 'Setting physics state variables column by column')

        ! Set variables in the `physics_state` derived type column by column.
        ! This way, peak memory usage can be reduced.
        do column_index = 1, ncells_solve
            call update_shared_variables(column_index)
            call set_physics_state_column(column_index)
        end do

        call set_physics_state_external()

        call final_shared_variables()

        call dyn_debug_print(debugout_debug, subname // ' completed')
    contains
        !> Initialize variables that are shared and repeatedly used by the `update_shared_variables` and
        !> `set_physics_state_column` internal subroutines.
        !> (KCW, 2024-07-20)
        subroutine init_shared_variables()
            ! Module(s) from CAM-SIMA.
            use cam_abortutils, only: check_allocate
            use cam_constituents, only: const_is_water_species, num_advected
            use dyn_comp, only: mpas_dynamical_core
            use vert_coord, only: pver, pverp

            character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::init_shared_variables'
            integer :: i
            integer :: ierr
            logical, allocatable :: is_water_species(:)

            call dyn_debug_print(debugout_info, 'Preparing for dynamics-physics coupling')

            nullify(index_qv)
            nullify(exner)
            nullify(rho_zz)
            nullify(scalars)
            nullify(theta_m)
            nullify(ucellzonal, ucellmeridional, w)
            nullify(zgrid)
            nullify(zz)

            allocate(is_water_species(num_advected), stat=ierr)
            call check_allocate(ierr, subname, &
                'is_water_species(num_advected)', &
                'dyn_coupling', __LINE__)

            do i = 1, num_advected
                is_water_species(i) = const_is_water_species(i)
            end do

            allocate(is_water_species_index(count(is_water_species)), stat=ierr)
            call check_allocate(ierr, subname, &
                'is_water_species_index(count(is_water_species))', &
                'dyn_coupling', __LINE__)

            is_water_species_index(:) = &
                pack([(mpas_dynamical_core % map_mpas_scalar_index(i), i = 1, num_advected)], is_water_species)

            deallocate(is_water_species)

            allocate(pd_int_col(pverp), pd_mid_col(pver), p_int_col(pverp), p_mid_col(pver), z_int_col(pverp), stat=ierr)
            call check_allocate(ierr, subname, &
                'pd_int_col(pverp), pd_mid_col(pver), p_int_col(pverp), p_mid_col(pver), z_int_col(pverp)', &
                'dyn_coupling', __LINE__)

            allocate(dpd_col(pver), dp_col(pver), dz_col(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                'dpd_col(pver), dp_col(pver), dz_col(pver)', &
                'dyn_coupling', __LINE__)

            allocate(qv_mid_col(pver), sigma_all_q_mid_col(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                'qv_mid_col(pver), sigma_all_q_mid_col(pver)', &
                'dyn_coupling', __LINE__)

            allocate(rhod_mid_col(pver), rho_mid_col(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                'rhod_mid_col(pver), rho_mid_col(pver)', &
                'dyn_coupling', __LINE__)

            allocate(t_mid_col(pver), tm_mid_col(pver), tv_mid_col(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                't_mid_col(pver), tm_mid_col(pver), tv_mid_col(pver)', &
                'dyn_coupling', __LINE__)

            allocate(u_mid_col(pver), v_mid_col(pver), omega_mid_col(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                'u_mid_col(pver), v_mid_col(pver), omega_mid_col(pver)', &
                'dyn_coupling', __LINE__)

            call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv')
            call mpas_dynamical_core % get_variable_pointer(exner, 'diag', 'exner')
            call mpas_dynamical_core % get_variable_pointer(rho_zz, 'state', 'rho_zz', time_level=1)
            call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1)
            call mpas_dynamical_core % get_variable_pointer(theta_m, 'state', 'theta_m', time_level=1)
            call mpas_dynamical_core % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal')
            call mpas_dynamical_core % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional')
            call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1)
            call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid')
            call mpas_dynamical_core % get_variable_pointer(zz, 'mesh', 'zz')
        end subroutine init_shared_variables

        !> Finalize variables that are shared and repeatedly used by the `update_shared_variables` and
        !> `set_physics_state_column` internal subroutines.
        !> (KCW, 2024-07-20)
        subroutine final_shared_variables()
            character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::final_shared_variables'

            deallocate(is_water_species_index)
            deallocate(pd_int_col, pd_mid_col, p_int_col, p_mid_col, z_int_col)
            deallocate(dpd_col, dp_col, dz_col)
            deallocate(qv_mid_col, sigma_all_q_mid_col)
            deallocate(rhod_mid_col, rho_mid_col)
            deallocate(t_mid_col, tm_mid_col, tv_mid_col)
            deallocate(u_mid_col, v_mid_col, omega_mid_col)

            nullify(index_qv)
            nullify(exner)
            nullify(rho_zz)
            nullify(scalars)
            nullify(theta_m)
            nullify(ucellzonal, ucellmeridional, w)
            nullify(zgrid)
            nullify(zz)
        end subroutine final_shared_variables

        !> Update variables for the specific column, indicated by `i`. This subroutine and `set_physics_state_column`
        !> should be called in pairs.
        !> (KCW, 2024-07-30)
        subroutine update_shared_variables(i)
            ! Module(s) from CAM-SIMA.
            use dynconst, only: constant_g => gravit, constant_rd => rair, constant_rv => rh2o
            use vert_coord, only: pver, pverp

            integer, intent(in) :: i

            character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::update_shared_variables'
            integer :: k

            ! The summation term of equation 5 in doi:10.1029/2017MS001257.
            sigma_all_q_mid_col(:) = 1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1)

            ! Compute thermodynamic variables.

            ! By definition.
            z_int_col(:) = real(zgrid(:, i), kind_r8)
            dz_col(:) = z_int_col(2:pverp) - z_int_col(1:pver)
            qv_mid_col(:) = real(scalars(index_qv, :, i), kind_r8)
            rhod_mid_col(:) = real(rho_zz(:, i) * zz(:, i), kind_r8)

            ! Equation 5 in doi:10.1029/2017MS001257.
            rho_mid_col(:) = rhod_mid_col(:) * sigma_all_q_mid_col(:)

            ! Hydrostatic equation.
            dpd_col(:) = -rhod_mid_col(:) * constant_g * dz_col(:)
            dp_col(:) = -rho_mid_col(:) * constant_g * dz_col(:)

            ! By definition of Exner function. Also see below.
            tm_mid_col(:) = real(theta_m(:, i) * exner(:, i), kind_r8)

            ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97.
            ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1.
            t_mid_col(:) = tm_mid_col(:) / &
                (1.0_kind_r8 + constant_rv / constant_rd * qv_mid_col(:))

            ! Equation 16 in doi:10.1029/2017MS001257.
            ! The numerator terms are just `tm_mid_col` here (i.e., modified "moist" temperature).
            tv_mid_col(:) = tm_mid_col(:) / sigma_all_q_mid_col(:)

            ! Equation of state.
            pd_mid_col(:) = rhod_mid_col(:) * constant_rd * t_mid_col(:)
            p_mid_col(:) = rho_mid_col(:) * constant_rd * tv_mid_col(:)

            ! By definition.
            pd_int_col(pverp) = pd_mid_col(pver) + 0.5_kind_r8 * dpd_col(pver)
            p_int_col(pverp) = p_mid_col(pver) + 0.5_kind_r8 * dp_col(pver)

            ! Integrate downward.
            do k = pver, 1, -1
                pd_int_col(k) = pd_int_col(k + 1) - dpd_col(k)
                p_int_col(k) = p_int_col(k + 1) - dp_col(k)
            end do

            ! Compute momentum variables.

            ! By definition.
            u_mid_col(:) = real(ucellzonal(:, i), kind_r8)
            v_mid_col(:) = real(ucellmeridional(:, i), kind_r8)
            omega_mid_col(:) = -rhod_mid_col(:) * constant_g * 0.5_kind_r8 * real(w(1:pver, i) + w(2:pverp, i), kind_r8)
        end subroutine update_shared_variables

        !> Set variables for the specific column, indicated by `i`, in the `physics_state` derived type.
        !> This subroutine and `update_shared_variables` should be called in pairs.
        !> (KCW, 2024-07-30)
        subroutine set_physics_state_column(i)
            ! Module(s) from CAM-SIMA.
            use dyn_comp, only: reverse
            use dynconst, only: constant_g => gravit
            use physics_types, only: phys_state

            integer, intent(in) :: i

            character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_column'

            ! Vertical index order is reversed between CAM-SIMA and MPAS.
            ! Always call `reverse` when assigning anything to/from the `physics_state` derived type.

            phys_state % u(i, :) = reverse(u_mid_col)
            phys_state % v(i, :) = reverse(v_mid_col)
            phys_state % omega(i, :) = reverse(omega_mid_col)

            phys_state % psdry(i) = pd_int_col(1)
            phys_state % pintdry(i, :) = reverse(pd_int_col)
            phys_state % pmiddry(i, :) = reverse(pd_mid_col)
            phys_state % pdeldry(i, :) = reverse(-dpd_col)
            phys_state % lnpintdry(i, :) = log(phys_state % pintdry(i, :))
            phys_state % lnpmiddry(i, :) = log(phys_state % pmiddry(i, :))
            phys_state % rpdeldry(i, :) = 1.0_kind_r8 / phys_state % pdeldry(i, :)

            phys_state % ps(i) = p_int_col(1)
            phys_state % pint(i, :) = reverse(p_int_col)
            phys_state % pmid(i, :) = reverse(p_mid_col)
            phys_state % pdel(i, :) = reverse(-dp_col)
            phys_state % lnpint(i, :) = log(phys_state % pint(i, :))
            phys_state % lnpmid(i, :) = log(phys_state % pmid(i, :))
            phys_state % rpdel(i, :) = 1.0_kind_r8 / phys_state % pdel(i, :)

            phys_state % t(i, :) = reverse(t_mid_col)

            phys_state % phis(i) = constant_g * z_int_col(1)
        end subroutine set_physics_state_column

        !> Set variables in the `physics_state` derived type by calling external procedures.
        !> (KCW, 2024-07-30)
        subroutine set_physics_state_external()
            ! Module(s) from CAM-SIMA.
            use cam_abortutils, only: check_allocate, endrun
            use cam_constituents, only: const_qmin, num_advected
            use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update
            use cam_thermo_formula, only: energy_formula_dycore_mpas
            use dyn_comp, only: mpas_dynamical_core
            use dynconst, only: constant_g => gravit
            use physics_types, only: cappav, cp_or_cv_dycore, cpairv, lagrangian_vertical, phys_state, rairv, zvirv
            use runtime_obj, only: cam_runtime_opts
            use string_utils, only: stringify
            use vert_coord, only: pver, pverp
            ! Module(s) from CCPP.
            use cam_ccpp_cap, only: cam_constituents_array, cam_model_const_properties
            use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
            use ccpp_kinds, only: kind_phys
            use geopotential_temp, only: geopotential_temp_run
            use qneg, only: qneg_run
            use static_energy, only: update_dry_static_energy_run
            ! Module(s) from CESM Share.
            use shr_kind_mod, only: len_cx => shr_kind_cx

            character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_external'
            character(len_cx) :: cerr
            integer :: i
            integer :: ierr
            real(kind_phys), allocatable :: minimum_constituents(:)
            real(kind_phys), pointer :: constituents(:, :, :)
            type(ccpp_constituent_prop_ptr_t), pointer :: constituent_properties(:)

            call dyn_debug_print(debugout_info, 'Setting physics state variables externally')

            nullify(constituents)
            nullify(constituent_properties)

            allocate(minimum_constituents(num_advected), stat=ierr)
            call check_allocate(ierr, subname, &
                'minimum_constituents(num_advected)', &
                'dyn_coupling', __LINE__)

            do i = 1, num_advected
                minimum_constituents(i) = const_qmin(i)
            end do

            constituents => cam_constituents_array()

            if (.not. associated(constituents)) then
                call endrun('Failed to find variable "constituents"', subname, __LINE__)
            end if

            constituent_properties => cam_model_const_properties()

            if (.not. associated(constituent_properties)) then
                call endrun('Failed to find variable "constituent_properties"', subname, __LINE__)
            end if

            ! Update `cappav`, `cpairv`, `rairv`, `zvirv`, etc. as needed by calling `cam_thermo_dry_air_update`.
            ! Note that this subroutine expects constituents to be dry.
            call cam_thermo_dry_air_update( &
                constituents, phys_state % t, ncells_solve, pver, cam_runtime_opts % update_thermodynamic_variables())

            ! Update `cp_or_cv_dycore` by calling `cam_thermo_water_update`.
            ! Note that this subroutine expects constituents to be dry.
            call cam_thermo_water_update( &
                constituents, ncells_solve, pver, energy_formula_dycore_mpas, cp_or_cv_dycore)

            ! This variable name is really misleading. It actually represents the reciprocal of Exner function
            ! with respect to surface pressure. This definition is sometimes used for boundary layer work. See
            ! the paragraph below equation 1.5.1c in doi:10.1007/978-94-009-3027-8.
            ! Also note that `cappav` is updated externally by `cam_thermo_dry_air_update`.
            do i = 1, ncells_solve
                phys_state % exner(i, :) = (phys_state % ps(i) / phys_state % pmid(i, :)) ** cappav(i, :)
            end do

            ! Note that constituents become moist after this.
            call dyn_exchange_constituent_states(direction='i', exchange=.false., conversion=.true.)

            ! Impose minimum limits on constituents.
            call qneg_run(subname, ncells_solve, pver, minimum_constituents, constituents, ierr, cerr)

            if (ierr /= 0) then
                call endrun('Failed to impose minimum limits on constituents externally' // new_line('') // &
                    'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
                    subname, __LINE__)
            end if

            ! Set `zi` (i.e., geopotential height at layer interfaces) and `zm` (i.e., geopotential height at layer midpoints).
            ! Note that `rairv` and `zvirv` are updated externally by `cam_thermo_dry_air_update`.
            call geopotential_temp_run( &
                pver, lagrangian_vertical, pver, 1, pverp, 1, num_advected, &
                phys_state % lnpint, phys_state % pint, phys_state % pmid, phys_state % pdel, phys_state % rpdel, phys_state % t, &
                constituents(:, :, mpas_dynamical_core % map_constituent_index(index_qv)), constituents, &
                constituent_properties, rairv, constant_g, zvirv, phys_state % zi, phys_state % zm, ncells_solve, ierr, cerr)

            if (ierr /= 0) then
                call endrun('Failed to set variable "zi" and "zm" externally' // new_line('') // &
                    'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
                    subname, __LINE__)
            end if

            ! Set `dse` (i.e., dry static energy).
            ! Note that `cpairv` is updated externally by `cam_thermo_dry_air_update`.
            call update_dry_static_energy_run( &
                pver, constant_g, phys_state % t, phys_state % zm, phys_state % phis, phys_state % dse, cpairv, ierr, cerr)

            if (ierr /= 0) then
                call endrun('Failed to set variable "dse" externally' // new_line('') // &
                    'External procedure returned with ' // stringify([ierr]) // ': ' // trim(adjustl(cerr)), &
                    subname, __LINE__)
            end if

            deallocate(minimum_constituents)

            nullify(constituents)
            nullify(constituent_properties)
        end subroutine set_physics_state_external
    end subroutine dynamics_to_physics_coupling