diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 875617b3..5cab4f4a 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -110,6 +110,7 @@ end subroutine model_error_if procedure, pass, public :: exchange_halo => dyn_mpas_exchange_halo procedure, pass, public :: compute_unit_vector => dyn_mpas_compute_unit_vector procedure, pass, public :: compute_edge_wind => dyn_mpas_compute_edge_wind + procedure, pass, public :: compute_cell_relative_vorticity => dyn_mpas_compute_cell_relative_vorticity procedure, pass, public :: init_phase4 => dyn_mpas_init_phase4 procedure, pass, public :: run => dyn_mpas_run procedure, pass, public :: final => dyn_mpas_final @@ -2776,6 +2777,106 @@ subroutine dyn_mpas_compute_edge_wind(self, wind_tendency) call self % debug_print(log_level_debug, subname // ' completed') end subroutine dyn_mpas_compute_edge_wind + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_compute_cell_relative_vorticity + ! + !> summary: Compute the relative vorticities at cell points. + !> author: Kuan-Chih Wang + !> date: 2025-04-12 + !> + !> MPAS uses staggered C-grid for spatial discretization, where relative + !> vorticities are located at vertex points because wind vectors are located at + !> edge points. However, physics schemes that use relative vorticities as input + !> usually want them at cell points instead. + !> This subroutine computes the relative vorticity at each cell point from its + !> surrounding vertex points and returns the results. + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_compute_cell_relative_vorticity(self, cell_relative_vorticity) + class(mpas_dynamical_core_type), intent(in) :: self + real(rkind), allocatable, intent(out) :: cell_relative_vorticity(:, :) + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_cell_relative_vorticity' + integer :: i, k + integer :: ierr + integer, pointer :: ncellssolve, nvertlevels + integer, pointer :: kiteforcell(:, :), nedgesoncell(:), verticesoncell(:, :) + real(rkind), pointer :: areacell(:), kiteareasonvertex(:, :), vorticity(:, :) + + nullify(ncellssolve, nvertlevels) + nullify(kiteforcell, nedgesoncell, verticesoncell) + nullify(areacell, kiteareasonvertex, vorticity) + + ! Input. + call self % get_variable_pointer(ncellssolve, 'dim', 'nCellsSolve') + call self % get_variable_pointer(nvertlevels, 'dim', 'nVertLevels') + + call self % get_variable_pointer(kiteforcell, 'mesh', 'kiteForCell') + call self % get_variable_pointer(nedgesoncell, 'mesh', 'nEdgesOnCell') + call self % get_variable_pointer(verticesoncell, 'mesh', 'verticesOnCell') + + call self % get_variable_pointer(areacell, 'mesh', 'areaCell') + call self % get_variable_pointer(kiteareasonvertex, 'mesh', 'kiteAreasOnVertex') + call self % get_variable_pointer(vorticity, 'diag', 'vorticity') + + ! Output. + allocate(cell_relative_vorticity(nvertlevels, ncellssolve), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate cell_relative_vorticity', subname, __LINE__) + end if + + do i = 1, ncellssolve + do k = 1, nvertlevels + cell_relative_vorticity(k, i) = regrid_from_vertex_to_cell(i, k, & + nedgesoncell, verticesoncell, kiteforcell, kiteareasonvertex, areacell, & + vorticity) + end do + end do + + nullify(ncellssolve, nvertlevels) + nullify(kiteforcell, nedgesoncell, verticesoncell) + nullify(areacell, kiteareasonvertex, vorticity) + end subroutine dyn_mpas_compute_cell_relative_vorticity + + !------------------------------------------------------------------------------- + ! function regrid_from_vertex_to_cell + ! + !> summary: Regrid values from vertex points to the specified cell point. + !> author: Kuan-Chih Wang + !> date: 2025-04-12 + !> + !> This function computes the area weighted average (i.e., `cell_value`) at the + !> specified cell point (i.e., `cell_index` and `cell_level`) from the values + !> at its surrounding vertex points (i.e., `vertex_value`). + !> The formulation used here is adapted and generalized from the + !> `atm_compute_solve_diagnostics` subroutine in MPAS. + ! + !------------------------------------------------------------------------------- + pure function regrid_from_vertex_to_cell(cell_index, cell_level, & + nverticesoncell, verticesoncell, kiteforcell, kiteareasonvertex, areacell, & + vertex_value) result(cell_value) + integer, intent(in) :: cell_index, cell_level + integer, intent(in) :: nverticesoncell(:), verticesoncell(:, :), kiteforcell(:, :) + real(rkind), intent(in) :: kiteareasonvertex(:, :), areacell(:) + real(rkind), intent(in) :: vertex_value(:, :) + real(rkind) :: cell_value + + integer :: i, j, vertex_index + + cell_value = 0.0_rkind + + do i = 1, nverticesoncell(cell_index) + j = kiteforcell(i, cell_index) + vertex_index = verticesoncell(i, cell_index) + + cell_value = cell_value + & + kiteareasonvertex(j, vertex_index) * vertex_value(cell_level, vertex_index) + end do + + cell_value = cell_value / areacell(cell_index) + end function regrid_from_vertex_to_cell + !------------------------------------------------------------------------------- ! subroutine dyn_mpas_init_phase4 !