Skip to content

Enable MPAS dynamical core to compute relative vorticities at cell points #394

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: development
Choose a base branch
from
Open

Enable MPAS dynamical core to compute relative vorticities at cell points #394

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 101 additions & 0 deletions src/dynamics/mpas/driver/dyn_mpas_subdriver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
!
Expand Down