physics_to_dynamics_coupling Subroutine

    subroutine physics_to_dynamics_coupling()
        ! Module(s) from CAM-SIMA.
        use cam_logfile, only: debugout_debug
        use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_states
        ! 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::physics_to_dynamics_coupling'
        integer, pointer :: index_qv
        real(kind_r8), allocatable :: qv_prev(:, :) ! Water vapor mixing ratio (kg kg-1)
                                                    ! before being updated by physics.
        real(kind_dyn_mpas), pointer :: rho_zz(:, :)
        real(kind_dyn_mpas), pointer :: scalars(:, :, :)
        real(kind_dyn_mpas), pointer :: zz(:, :)

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

        call init_shared_variables()

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

        call set_mpas_physics_tendency_ru()
        call set_mpas_physics_tendency_rho()
        call set_mpas_physics_tendency_rtheta()

        call final_shared_variables()

        call dyn_debug_print(debugout_debug, subname // ' completed')
    contains
        !> Initialize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines.
        !> (KCW, 2024-09-13)
        subroutine init_shared_variables()
            ! Module(s) from CAM-SIMA.
            use cam_abortutils, only: check_allocate
            use cam_logfile, only: debugout_info
            use dyn_comp, only: mpas_dynamical_core, ncells_solve
            use vert_coord, only: pver

            character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::init_shared_variables'
            integer :: ierr

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

            nullify(index_qv)
            nullify(rho_zz)
            nullify(scalars)
            nullify(zz)

            call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv')
            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(zz, 'mesh', 'zz')

            allocate(qv_prev(pver, ncells_solve), stat=ierr)
            call check_allocate(ierr, subname, &
                'qv_prev(pver, ncells_solve)', &
                'dyn_coupling', __LINE__)

            ! Save water vapor mixing ratio before being updated by physics because `set_mpas_physics_tendency_rtheta`
            ! needs it. This must be done before calling `dyn_exchange_constituent_states`.
            qv_prev(:, :) = real(scalars(index_qv, :, 1:ncells_solve), kind_r8)
        end subroutine init_shared_variables

        !> Finalize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines.
        !> (KCW, 2024-09-13)
        subroutine final_shared_variables()
            character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::final_shared_variables'

            deallocate(qv_prev)

            nullify(index_qv)
            nullify(rho_zz)
            nullify(scalars)
            nullify(zz)
        end subroutine final_shared_variables

        !> Set MPAS physics tendency `tend_ru_physics` (i.e., "coupled" tendency of horizontal velocity at edge interfaces
        !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`.
        !> (KCW, 2024-09-11)
        subroutine set_mpas_physics_tendency_ru()
            ! Module(s) from CAM-SIMA.
            use cam_logfile, only: debugout_info
            use dyn_comp, only: reverse, mpas_dynamical_core, ncells_solve
            use physics_types, only: phys_tend

            character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_ru'
            integer :: i
            real(kind_dyn_mpas), pointer :: u_tendency(:, :), v_tendency(:, :)

            call dyn_debug_print(debugout_info, 'Setting MPAS physics tendency "tend_ru_physics"')

            nullify(u_tendency, v_tendency)

            call mpas_dynamical_core % get_variable_pointer(u_tendency, 'tend_physics', 'tend_uzonal')
            call mpas_dynamical_core % get_variable_pointer(v_tendency, 'tend_physics', 'tend_umerid')

            ! Vertical index order is reversed between CAM-SIMA and MPAS.
            ! Always call `reverse` when assigning anything to/from the `physics_tend` derived type.
            do i = 1, ncells_solve
                u_tendency(:, i) = real(reverse(phys_tend % dudt_total(i, :)) * real(rho_zz(:, i), kind_r8), kind_dyn_mpas)
                v_tendency(:, i) = real(reverse(phys_tend % dvdt_total(i, :)) * real(rho_zz(:, i), kind_r8), kind_dyn_mpas)
            end do

            nullify(u_tendency, v_tendency)

            call mpas_dynamical_core % compute_edge_wind(wind_tendency=.true.)
        end subroutine set_mpas_physics_tendency_ru

        !> Set MPAS physics tendency `tend_rho_physics` (i.e., "coupled" tendency of dry air density due to physics).
        !> In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`.
        !> (KCW, 2024-09-11)
        subroutine set_mpas_physics_tendency_rho()
            ! Module(s) from CAM-SIMA.
            use cam_logfile, only: debugout_info
            use dyn_comp, only: mpas_dynamical_core, ncells_solve

            character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rho'
            real(kind_dyn_mpas), pointer :: rho_tendency(:, :)

            call dyn_debug_print(debugout_info, 'Setting MPAS physics tendency "tend_rho_physics"')

            nullify(rho_tendency)

            call mpas_dynamical_core % get_variable_pointer(rho_tendency, 'tend_physics', 'tend_rho_physics')

            ! The material derivative of `rho` (i.e., dry air density) is zero for incompressible fluid.
            rho_tendency(:, 1:ncells_solve) = 0.0_kind_dyn_mpas

            nullify(rho_tendency)

            ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
            call mpas_dynamical_core % exchange_halo('tend_rho_physics')
        end subroutine set_mpas_physics_tendency_rho

        !> Set MPAS physics tendency `tend_rtheta_physics` (i.e., "coupled" tendency of modified "moist" potential temperature
        !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`.
        !> (KCW, 2024-09-19)
        subroutine set_mpas_physics_tendency_rtheta()
            ! Module(s) from CAM-SIMA.
            use cam_abortutils, only: check_allocate
            use cam_logfile, only: debugout_info
            use dyn_comp, only: reverse, mpas_dynamical_core, ncells_solve
            use dynconst, only: constant_rd => rair, constant_rv => rh2o
            use physics_types, only: dtime_phys, phys_tend
            use vert_coord, only: pver

            character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rtheta'
            integer :: i
            integer :: ierr
            ! Variable name suffixes have the following meanings:
            ! `*_col`: Variable is of each column.
            ! `*_prev`: Variable is before being updated by physics.
            ! `*_curr`: Variable is after being updated by physics.
            real(kind_r8), allocatable :: qv_col_prev(:), qv_col_curr(:)         ! Water vapor mixing ratio (kg kg-1).
            real(kind_r8), allocatable :: rhod_col(:)                            ! Dry air density (kg m-3).
            real(kind_r8), allocatable :: t_col_prev(:), t_col_curr(:)           ! Temperature (K).
            real(kind_r8), allocatable :: theta_col_prev(:), theta_col_curr(:)   ! Potential temperature (K).
            real(kind_r8), allocatable :: thetam_col_prev(:), thetam_col_curr(:) ! Modified "moist" potential temperature (K).
            real(kind_dyn_mpas), pointer :: theta_m(:, :)
            real(kind_dyn_mpas), pointer :: theta_m_tendency(:, :)

            call dyn_debug_print(debugout_info, 'Setting MPAS physics tendency "tend_rtheta_physics"')

            nullify(theta_m)
            nullify(theta_m_tendency)

            allocate(qv_col_prev(pver), qv_col_curr(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                'qv_col_prev(pver), qv_col_curr(pver)', &
                'dyn_coupling', __LINE__)

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

            allocate(t_col_prev(pver), t_col_curr(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                't_col_prev(pver), t_col_curr(pver)', &
                'dyn_coupling', __LINE__)

            allocate(theta_col_prev(pver), theta_col_curr(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                'theta_col_prev(pver), theta_col_curr(pver)', &
                'dyn_coupling', __LINE__)

            allocate(thetam_col_prev(pver), thetam_col_curr(pver), stat=ierr)
            call check_allocate(ierr, subname, &
                'thetam_col_prev(pver), thetam_col_curr(pver)', &
                'dyn_coupling', __LINE__)

            call mpas_dynamical_core % get_variable_pointer(theta_m, 'state', 'theta_m', time_level=1)
            call mpas_dynamical_core % get_variable_pointer(theta_m_tendency, 'tend_physics', 'tend_rtheta_physics')

            ! Set `theta_m_tendency` column by column. This way, peak memory usage can be reduced.
            do i = 1, ncells_solve
                qv_col_curr(:) = real(scalars(index_qv, :, i), kind_r8)
                qv_col_prev(:) = qv_prev(:, i)
                rhod_col(:) = real(rho_zz(:, i) * zz(:, i), kind_r8)

                thetam_col_prev(:) = real(theta_m(:, i), kind_r8)
                theta_col_prev(:) = thetam_col_prev(:) / (1.0_kind_r8 + constant_rv / constant_rd * qv_col_prev(:))
                t_col_prev(:) = t_of_theta_rhod_qv(theta_col_prev, rhod_col, qv_col_prev)

                ! Vertical index order is reversed between CAM-SIMA and MPAS.
                ! Always call `reverse` when assigning anything to/from the `physics_tend` derived type.
                t_col_curr(:) = t_col_prev(:) + reverse(phys_tend % dtdt_total(i, :)) * dtime_phys
                theta_col_curr(:) = theta_of_t_rhod_qv(t_col_curr, rhod_col, qv_col_curr)
                thetam_col_curr(:) = theta_col_curr(:) * (1.0_kind_r8 + constant_rv / constant_rd * qv_col_curr(:))

                theta_m_tendency(:, i) = &
                    real((thetam_col_curr(:) - thetam_col_prev(:)) * real(rho_zz(:, i), kind_r8) / dtime_phys, kind_dyn_mpas)
            end do

            deallocate(qv_col_prev, qv_col_curr)
            deallocate(rhod_col)
            deallocate(t_col_prev, t_col_curr)
            deallocate(theta_col_prev, theta_col_curr)
            deallocate(thetam_col_prev, thetam_col_curr)

            nullify(theta_m)
            nullify(theta_m_tendency)

            ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
            call mpas_dynamical_core % exchange_halo('tend_rtheta_physics')
        end subroutine set_mpas_physics_tendency_rtheta

        !> Compute temperature `t` as a function of potential temperature `theta`, dry air density `rhod` and water vapor
        !> mixing ratio `qv`. Essentially,
        !> \( T = \theta^{\frac{C_p}{C_v}} [\frac{\rho_d R_d (1 + \frac{R_v}{R_d} q_v)}{P_0}]^{\frac{R_d}{C_v}} \).
        !> The formulation comes from Poisson equation with equation of state plugged in and arranging
        !> for temperature. This function is the exact inverse of `theta_of_t_rhod_qv`, which means that:
        !> `t == t_of_theta_rhod_qv(theta_of_t_rhod_qv(t, rhod, qv), rhod, qv)`.
        !> (KCW, 2024-09-13)
        pure elemental function t_of_theta_rhod_qv(theta, rhod, qv) result(t)
            ! Module(s) from CAM-SIMA.
            use dynconst, only: constant_cpd => cpair, constant_p0 => pref, &
                                constant_rd => rair, constant_rv => rh2o

            real(kind_r8), intent(in) :: theta, rhod, qv
            real(kind_r8) :: t

            real(kind_r8) :: constant_cvd ! Specific heat of dry air at constant volume.

            ! Mayer's relation.
            constant_cvd = constant_cpd - constant_rd

            ! Poisson equation with equation of state plugged in and arranging for temperature. For equation of state,
            ! it can be shown that the effect of water vapor can be passed on to the temperature term entirely such that
            ! dry air density and dry air gas constant can be used at all times. This modified "moist" temperature is
            ! described herein:
            ! 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.
            !
            ! In all, solve the below equation set for $T$ in terms of $\theta$, $\rho_d$ and $q_v$:
            ! \begin{equation*}
            !     \begin{cases}
            !         \theta &= T (\frac{P_0}{P})^{\frac{R_d}{C_p}} \\
            !         P &= \rho_d R_d T_m \\
            !         T_m &= T (1 + \frac{R_v}{R_d} q_v)
            !     \end{cases}
            ! \end{equation*}
            t = (theta ** (constant_cpd / constant_cvd)) * &
                (((rhod * constant_rd * (1.0_kind_r8 + constant_rv / constant_rd * qv)) / constant_p0) ** &
                (constant_rd / constant_cvd))
        end function t_of_theta_rhod_qv

        !> Compute potential temperature `theta` as a function of temperature `t`, dry air density `rhod` and water vapor
        !> mixing ratio `qv`. Essentially,
        !> \( \theta = T^{\frac{C_v}{C_p}} [\frac{P_0}{\rho_d R_d (1 + \frac{R_v}{R_d} q_v)}]^{\frac{R_d}{C_p}} \).
        !> The formulation comes from Poisson equation with equation of state plugged in and arranging
        !> for potential temperature. This function is the exact inverse of `t_of_theta_rhod_qv`, which means that:
        !> `theta == theta_of_t_rhod_qv(t_of_theta_rhod_qv(theta, rhod, qv), rhod, qv)`.
        !> (KCW, 2024-09-13)
        pure elemental function theta_of_t_rhod_qv(t, rhod, qv) result(theta)
            ! Module(s) from CAM-SIMA.
            use dynconst, only: constant_cpd => cpair, constant_p0 => pref, &
                                constant_rd => rair, constant_rv => rh2o

            real(kind_r8), intent(in) :: t, rhod, qv
            real(kind_r8) :: theta

            real(kind_r8) :: constant_cvd ! Specific heat of dry air at constant volume.

            ! Mayer's relation.
            constant_cvd = constant_cpd - constant_rd

            ! Poisson equation with equation of state plugged in and arranging for potential temperature. For equation of state,
            ! it can be shown that the effect of water vapor can be passed on to the temperature term entirely such that
            ! dry air density and dry air gas constant can be used at all times. This modified "moist" temperature is
            ! described herein:
            ! 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.
            !
            ! In all, solve the below equation set for $\theta$ in terms of $T$, $\rho_d$ and $q_v$:
            ! \begin{equation*}
            !     \begin{cases}
            !         \theta &= T (\frac{P_0}{P})^{\frac{R_d}{C_p}} \\
            !         P &= \rho_d R_d T_m \\
            !         T_m &= T (1 + \frac{R_v}{R_d} q_v)
            !     \end{cases}
            ! \end{equation*}
            theta = (t ** (constant_cvd / constant_cpd)) * &
                ((constant_p0 / (rhod * constant_rd * (1.0_kind_r8 + constant_rv / constant_rd * qv))) ** &
                (constant_rd / constant_cpd))
        end function theta_of_t_rhod_qv
    end subroutine physics_to_dynamics_coupling