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

Merged
merged 2 commits into from
Jul 3, 2025
Merged

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

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
17216d6
Enable MPAS subdriver to compute relative vorticities at cell points
kuanchihwang May 21, 2025
390c9c8
Merge branch 'development' into develop/mpas-dycore-cell-rel-vort
kuanchihwang Jun 25, 2025
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 @@ -2777,6 +2778,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)
Copy link
Collaborator

@nusbaume nusbaume Jun 1, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just an FYI that we have started trying to provide the Fortran-provided error message from allocate statements as well, e.g.:

character(256) :: errmsg
...
allocate(cell_relative_vorticity(nvertlevels, ncellssolve), stat=ierr, errmsg=errmsg)

Of course this change certainly isn't required for this PR.

Copy link
Collaborator Author

@kuanchihwang kuanchihwang Jun 2, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I have noticed it in recent PR reviews. I will open a dedicated PR to implement this enhancement.


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