module subroutine dyn_exchange_constituent_states(direction, exchange, conversion)
! Module(s) from CAM-SIMA.
use cam_abortutils, only: check_allocate, endrun
use cam_constituents, only: const_is_dry, const_is_water_species, num_advected
use cam_logfile, only: debugout_debug, debugout_info
use dyn_comp, only: dyn_debug_print, kind_dyn_mpas, mpas_dynamical_core
use dyn_grid, only: ncells_solve
use dyn_procedures, only: reverse
use physics_types, only: phys_state
use vert_coord, only: pver
! Module(s) from CCPP.
use cam_ccpp_cap, only: cam_constituents_array
use ccpp_kinds, only: kind_phys
! Module(s) from CESM Share.
use shr_kind_mod, only: kind_r8 => shr_kind_r8, &
len_cx => shr_kind_cx
character(*), intent(in) :: direction
logical, intent(in) :: exchange
logical, intent(in) :: conversion
character(*), parameter :: subname = 'dyn_coupling::dyn_exchange_constituent_states'
character(len_cx) :: cerr
integer :: i, j
integer :: ierr
integer, allocatable :: is_water_species_index(:)
logical, allocatable :: is_conversion_needed(:)
logical, allocatable :: is_water_species(:)
real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory.
real(kind_r8), allocatable :: sigma_all_q(:) ! Summation of all water species mixing ratios.
real(kind_dyn_mpas), pointer :: scalars(:, :, :) ! This points to MPAS memory.
call dyn_debug_print(debugout_debug, subname // ' entered')
select case (trim(adjustl(direction)))
case ('e', 'export')
if (exchange) then
call dyn_debug_print(debugout_info, 'Setting MPAS state "scalars" from physics state "constituents"')
end if
if (conversion) then
call dyn_debug_print(debugout_info, 'Converting MPAS state "scalars"')
end if
case ('i', 'import')
if (exchange) then
call dyn_debug_print(debugout_info, 'Setting physics state "constituents" from MPAS state "scalars"')
end if
if (conversion) then
call dyn_debug_print(debugout_info, 'Converting physics state "constituents"')
end if
case default
call endrun('Unsupported exchange direction "' // trim(adjustl(direction)) // '"', subname, __LINE__)
end select
nullify(constituents)
nullify(scalars)
allocate(is_conversion_needed(num_advected), errmsg=cerr, stat=ierr)
call check_allocate(ierr, subname, &
'is_conversion_needed(num_advected)', &
file='dyn_coupling', line=__LINE__, errmsg=trim(adjustl(cerr)))
allocate(is_water_species(num_advected), errmsg=cerr, stat=ierr)
call check_allocate(ierr, subname, &
'is_water_species(num_advected)', &
file='dyn_coupling', line=__LINE__, errmsg=trim(adjustl(cerr)))
do j = 1, num_advected
! All constituent mixing ratios in MPAS are dry.
! Therefore, conversion in between is needed for any constituent mixing ratios that are not dry in CAM-SIMA.
is_conversion_needed(j) = .not. const_is_dry(j)
is_water_species(j) = const_is_water_species(j)
end do
allocate(is_water_species_index(count(is_water_species)), errmsg=cerr, stat=ierr)
call check_allocate(ierr, subname, &
'is_water_species_index(count(is_water_species))', &
file='dyn_coupling', line=__LINE__, errmsg=trim(adjustl(cerr)))
allocate(sigma_all_q(pver), errmsg=cerr, stat=ierr)
call check_allocate(ierr, subname, &
'sigma_all_q(pver)', &
file='dyn_coupling', line=__LINE__, errmsg=trim(adjustl(cerr)))
constituents => cam_constituents_array()
if (.not. associated(constituents)) then
call endrun('Failed to find variable "constituents"', subname, __LINE__)
end if
call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1)
if (trim(adjustl(direction)) == 'e' .or. trim(adjustl(direction)) == 'export') then
do i = 1, ncells_solve
if (conversion .and. any(is_conversion_needed)) then
! The summation term of equation 8 in doi:10.1029/2017MS001257.
! Using equation 7 here is not possible because it requires all constituent mixing ratio to be moist
! on the RHS of it. There is no such guarantee in CAM-SIMA.
sigma_all_q(:) = reverse(phys_state % pdel(i, :) / phys_state % pdeldry(i, :))
end if
! `j` is indexing into `scalars`, so it is regarded as MPAS scalar index.
do j = 1, num_advected
if (exchange) then
! Vertical index order is reversed between CAM-SIMA and MPAS.
scalars(j, :, i) = &
real(reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j))), kind_dyn_mpas)
end if
if (conversion .and. is_conversion_needed(mpas_dynamical_core % map_constituent_index(j))) then
! Equation 8 in doi:10.1029/2017MS001257.
scalars(j, :, i) = &
real(real(scalars(j, :, i), kind_r8) * sigma_all_q(:), kind_dyn_mpas)
end if
end do
end do
else
is_water_species_index(:) = &
pack([(mpas_dynamical_core % map_mpas_scalar_index(i), i = 1, num_advected)], is_water_species)
do i = 1, ncells_solve
if (conversion .and. any(is_conversion_needed)) then
! The summation term of equation 8 in doi:10.1029/2017MS001257.
sigma_all_q(:) = reverse(1.0_kind_r8 + sum(real(scalars(is_water_species_index, :, i), kind_r8), 1))
end if
! `j` is indexing into `constituents`, so it is regarded as constituent index.
do j = 1, num_advected
if (exchange) then
! Vertical index order is reversed between CAM-SIMA and MPAS.
constituents(i, :, j) = &
reverse(real(scalars(mpas_dynamical_core % map_mpas_scalar_index(j), :, i), kind_r8))
end if
if (conversion .and. is_conversion_needed(j)) then
! Equation 8 in doi:10.1029/2017MS001257.
constituents(i, :, j) = &
constituents(i, :, j) / sigma_all_q(:)
end if
end do
end do
end if
deallocate(is_conversion_needed)
deallocate(is_water_species)
deallocate(is_water_species_index)
deallocate(sigma_all_q)
nullify(constituents)
nullify(scalars)
if (trim(adjustl(direction)) == 'e' .or. trim(adjustl(direction)) == 'export') then
! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually.
call mpas_dynamical_core % exchange_halo('scalars')
end if
call dyn_debug_print(debugout_debug, subname // ' completed')
end subroutine dyn_exchange_constituent_states