sima0_13_002: Refactor vertical coordinate computation in MPAS dynamical core

src/dynamics/mpas/driver/dyn_mpas_procedures.F90

@@ -15,7 +15,7 @@ module dyn_mpas_procedures
 
     private
     ! Computational procedures.
-    ! None in this group for now.
+    public :: dzw_of_rdzw, dzu_of_dzw, zu_of_dzw, zw_of_dzw
     ! Utility procedures.
     public :: almost_divisible
     public :: almost_equal
@@ -25,6 +25,26 @@ module dyn_mpas_procedures
     public :: stringify
     public :: tokenize
 
+    interface dzw_of_rdzw
+        module procedure dzw_of_rdzw_real32
+        module procedure dzw_of_rdzw_real64
+    end interface dzw_of_rdzw
+
+    interface dzu_of_dzw
+        module procedure dzu_of_dzw_real32
+        module procedure dzu_of_dzw_real64
+    end interface dzu_of_dzw
+
+    interface zu_of_dzw
+        module procedure zu_of_dzw_real32
+        module procedure zu_of_dzw_real64
+    end interface zu_of_dzw
+
+    interface zw_of_dzw
+        module procedure zw_of_dzw_real32
+        module procedure zw_of_dzw_real64
+    end interface zw_of_dzw
+
     interface almost_divisible
         module procedure almost_divisible_real32
         module procedure almost_divisible_real64
@@ -47,6 +67,138 @@ module dyn_mpas_procedures
         module procedure tokenize_into_tokens_separator
     end interface tokenize
 contains
+    !> Compute the differences in \( \zeta \) between w-wind levels `dzw` from the reciprocal differences in \( \zeta \) between
+    !> w-wind levels `rdzw`, where \( \zeta \) is the vertical coordinate, and w-wind levels are synonymous with layer interfaces
+    !> in MPAS.
+    !> (KCW, 2025-10-20)
+    pure elemental function dzw_of_rdzw_real32(rdzw) result(dzw)
+        use, intrinsic :: iso_fortran_env, only: real32
+
+        real(real32), intent(in) :: rdzw
+        real(real32) :: dzw
+
+        dzw = 1.0_real32 / rdzw
+    end function dzw_of_rdzw_real32
+
+    !> Compute the differences in \( \zeta \) between w-wind levels `dzw` from the reciprocal differences in \( \zeta \) between
+    !> w-wind levels `rdzw`, where \( \zeta \) is the vertical coordinate, and w-wind levels are synonymous with layer interfaces
+    !> in MPAS.
+    !> (KCW, 2025-10-20)
+    pure elemental function dzw_of_rdzw_real64(rdzw) result(dzw)
+        use, intrinsic :: iso_fortran_env, only: real64
+
+        real(real64), intent(in) :: rdzw
+        real(real64) :: dzw
+
+        dzw = 1.0_real64 / rdzw
+    end function dzw_of_rdzw_real64
+
+    !> Compute the differences in \( \zeta \) between u-wind levels `dzu` from the differences in \( \zeta \) between
+    !> w-wind levels `dzw`, where \( \zeta \) is the vertical coordinate, u-wind and w-wind levels are synonymous with
+    !> layer midpoints and interfaces in MPAS, respectively.
+    !> (KCW, 2025-10-20)
+    pure function dzu_of_dzw_real32(dzw) result(dzu)
+        use, intrinsic :: iso_fortran_env, only: real32
+
+        real(real32), intent(in) :: dzw(:)
+        real(real32) :: dzu(size(dzw) - 1)
+
+        integer :: k
+
+        do k = 1, size(dzu)
+            dzu(k) = 0.5_real32 * (dzw(k) + dzw(k + 1))
+        end do
+    end function dzu_of_dzw_real32
+
+    !> Compute the differences in \( \zeta \) between u-wind levels `dzu` from the differences in \( \zeta \) between
+    !> w-wind levels `dzw`, where \( \zeta \) is the vertical coordinate, u-wind and w-wind levels are synonymous with
+    !> layer midpoints and interfaces in MPAS, respectively.
+    !> (KCW, 2025-10-20)
+    pure function dzu_of_dzw_real64(dzw) result(dzu)
+        use, intrinsic :: iso_fortran_env, only: real64
+
+        real(real64), intent(in) :: dzw(:)
+        real(real64) :: dzu(size(dzw) - 1)
+
+        integer :: k
+
+        do k = 1, size(dzu)
+            dzu(k) = 0.5_real64 * (dzw(k) + dzw(k + 1))
+        end do
+    end function dzu_of_dzw_real64
+
+    !> Compute the \( \zeta \) coordinates at u-wind levels `zu` from the differences in \( \zeta \) between w-wind levels `dzw`,
+    !> where \( \zeta \) is the vertical coordinate, u-wind and w-wind levels are synonymous with layer midpoints and interfaces
+    !> in MPAS, respectively.
+    !> (KCW, 2025-10-20)
+    pure function zu_of_dzw_real32(dzw) result(zu)
+        use, intrinsic :: iso_fortran_env, only: real32
+
+        real(real32), intent(in) :: dzw(:)
+        real(real32) :: zu(size(dzw))
+
+        integer :: k
+        real(real32) :: zw(size(dzw) + 1)
+
+        zw(:) = zw_of_dzw(dzw)
+
+        do k = 1, size(zu)
+            zu(k) = 0.5_real32 * (zw(k) + zw(k + 1))
+        end do
+    end function zu_of_dzw_real32
+
+    !> Compute the \( \zeta \) coordinates at u-wind levels `zu` from the differences in \( \zeta \) between w-wind levels `dzw`,
+    !> where \( \zeta \) is the vertical coordinate, u-wind and w-wind levels are synonymous with layer midpoints and interfaces
+    !> in MPAS, respectively.
+    !> (KCW, 2025-10-20)
+    pure function zu_of_dzw_real64(dzw) result(zu)
+        use, intrinsic :: iso_fortran_env, only: real64
+
+        real(real64), intent(in) :: dzw(:)
+        real(real64) :: zu(size(dzw))
+
+        integer :: k
+        real(real64) :: zw(size(dzw) + 1)
+
+        zw(:) = zw_of_dzw(dzw)
+
+        do k = 1, size(zu)
+            zu(k) = 0.5_real64 * (zw(k) + zw(k + 1))
+        end do
+    end function zu_of_dzw_real64
+
+    !> Compute the \( \zeta \) coordinates at w-wind levels `zw` from the differences in \( \zeta \) between w-wind levels `dzw`,
+    !> where \( \zeta \) is the vertical coordinate, and w-wind levels are synonymous with layer interfaces in MPAS.
+    !> (KCW, 2025-10-20)
+    pure function zw_of_dzw_real32(dzw) result(zw)
+        use, intrinsic :: iso_fortran_env, only: real32
+
+        real(real32), intent(in) :: dzw(:)
+        real(real32) :: zw(size(dzw) + 1)
+
+        integer :: k
+
+        do k = 1, size(zw)
+            zw(k) = sum(dzw(1:k - 1))
+        end do
+    end function zw_of_dzw_real32
+
+    !> Compute the \( \zeta \) coordinates at w-wind levels `zw` from the differences in \( \zeta \) between w-wind levels `dzw`,
+    !> where \( \zeta \) is the vertical coordinate, and w-wind levels are synonymous with layer interfaces in MPAS.
+    !> (KCW, 2025-10-20)
+    pure function zw_of_dzw_real64(dzw) result(zw)
+        use, intrinsic :: iso_fortran_env, only: real64
+
+        real(real64), intent(in) :: dzw(:)
+        real(real64) :: zw(size(dzw) + 1)
+
+        integer :: k
+
+        do k = 1, size(zw)
+            zw(k) = sum(dzw(1:k - 1))
+        end do
+    end function zw_of_dzw_real64
+
     !> Test if `a` is divisible by `b`, where `a` and `b` are both reals.
     !> (KCW, 2024-05-25)
     pure elemental function almost_divisible_real32(a, b) result(almost_divisible)

src/dynamics/mpas/driver/dyn_mpas_subdriver.F90

@@ -2441,6 +2441,7 @@ subroutine dyn_mpas_compute_unit_vector(self)
         class(mpas_dynamical_core_type), intent(in) :: self
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_unit_vector'
+        real(rkind), parameter :: pi = acos(-1.0_rkind) ! Pi in native precision.
         integer :: i
         integer, pointer :: ncells
         real(rkind), pointer :: latcell(:), loncell(:)
@@ -2467,6 +2468,12 @@ subroutine dyn_mpas_compute_unit_vector(self)
 
         call self % debug_print(log_level_info, 'Computing unit vectors')
 
+        ! Make sure longitude values are within the range of [0, 2 * pi) according to MPAS mesh specifications.
+        ! Although stand-alone MPAS can still accept longitude values outside this range without issues, CAM-SIMA is less
+        ! tolerant about it. Note that there is a distinction between the `mod` and `modulo` intrinsic procedures.
+        ! The former uses truncated division, while the latter uses floored division.
+        loncell(:) = modulo(loncell, 2.0_rkind * pi)
+
         do i = 1, ncells
             east(1, i) = -sin(loncell(i))
             east(2, i) =  cos(loncell(i))
@@ -3370,7 +3377,7 @@ end function dyn_mpas_map_constituent_index
     subroutine dyn_mpas_get_local_mesh_dimension(self, &
             ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels)
         class(mpas_dynamical_core_type), intent(in) :: self
-        integer, intent(out) :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels
+        integer, optional, intent(out) :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_local_mesh_dimension'
         integer, pointer :: ncells_pointer
@@ -3401,16 +3408,35 @@ subroutine dyn_mpas_get_local_mesh_dimension(self, &
         call self % get_variable_pointer(nverticessolve_pointer, 'dim', 'nVerticesSolve')
         call self % get_variable_pointer(nvertlevels_pointer, 'dim', 'nVertLevels')
 
-        ncells = ncells_pointer                  ! Number of cells, including halo cells.
-        ncells_solve = ncellssolve_pointer       ! Number of cells, excluding halo cells.
-        nedges = nedges_pointer                  ! Number of edges, including halo edges.
-        nedges_solve = nedgessolve_pointer       ! Number of edges, excluding halo edges.
-        nvertices = nvertices_pointer            ! Number of vertices, including halo vertices.
-        nvertices_solve = nverticessolve_pointer ! Number of vertices, excluding halo vertices.
+        if (present(ncells)) then
+            ncells = ncells_pointer                  ! Number of cells, including halo cells.
+        end if
+
+        if (present(ncells_solve)) then
+            ncells_solve = ncellssolve_pointer       ! Number of cells, excluding halo cells.
+        end if
+
+        if (present(nedges)) then
+            nedges = nedges_pointer                  ! Number of edges, including halo edges.
+        end if
+
+        if (present(nedges_solve)) then
+            nedges_solve = nedgessolve_pointer       ! Number of edges, excluding halo edges.
+        end if
+
+        if (present(nvertices)) then
+            nvertices = nvertices_pointer            ! Number of vertices, including halo vertices.
+        end if
+
+        if (present(nvertices_solve)) then
+            nvertices_solve = nverticessolve_pointer ! Number of vertices, excluding halo vertices.
+        end if
 
         ! Vertical levels are not decomposed.
         ! All tasks have the same number of vertical levels.
-        nvertlevels = nvertlevels_pointer
+        if (present(nvertlevels)) then
+            nvertlevels = nvertlevels_pointer
+        end if
 
         nullify(ncells_pointer)
         nullify(ncellssolve_pointer)
@@ -3445,8 +3471,8 @@ subroutine dyn_mpas_get_global_mesh_dimension(self, &
         use mpas_dmpar, only: mpas_dmpar_max_int, mpas_dmpar_sum_int
 
         class(mpas_dynamical_core_type), intent(in) :: self
-        integer, intent(out) :: ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max
-        real(rkind), intent(out) :: sphere_radius
+        integer, optional, intent(out) :: ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max
+        real(rkind), optional, intent(out) :: sphere_radius
 
         character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_global_mesh_dimension'
         integer, pointer :: maxedges_pointer
@@ -3471,18 +3497,35 @@ subroutine dyn_mpas_get_global_mesh_dimension(self, &
         call self % get_variable_pointer(nverticessolve_pointer, 'dim', 'nVerticesSolve')
         call self % get_variable_pointer(nvertlevels_pointer, 'dim', 'nVertLevels')
 
-        call mpas_dmpar_sum_int(self % domain_ptr % dminfo, ncellssolve_pointer, ncells_global)
-        call mpas_dmpar_sum_int(self % domain_ptr % dminfo, nedgessolve_pointer, nedges_global)
-        call mpas_dmpar_sum_int(self % domain_ptr % dminfo, nverticessolve_pointer, nvertices_global)
+        if (present(ncells_global)) then
+            call mpas_dmpar_sum_int(self % domain_ptr % dminfo, ncellssolve_pointer, ncells_global)
+        end if
+
+        if (present(nedges_global)) then
+            call mpas_dmpar_sum_int(self % domain_ptr % dminfo, nedgessolve_pointer, nedges_global)
+        end if
+
+        if (present(nvertices_global)) then
+            call mpas_dmpar_sum_int(self % domain_ptr % dminfo, nverticessolve_pointer, nvertices_global)
+        end if
 
         ! Vertical levels are not decomposed.
         ! All tasks have the same number of vertical levels.
-        nvertlevels = nvertlevels_pointer
+        if (present(nvertlevels)) then
+            nvertlevels = nvertlevels_pointer
+        end if
 
-        call mpas_dmpar_max_int(self % domain_ptr % dminfo, ncellssolve_pointer, ncells_max)
+        if (present(ncells_max)) then
+            call mpas_dmpar_max_int(self % domain_ptr % dminfo, ncellssolve_pointer, ncells_max)
+        end if
 
-        nedges_max = maxedges_pointer
-        sphere_radius = self % domain_ptr % sphere_radius
+        if (present(nedges_max)) then
+            nedges_max = maxedges_pointer
+        end if
+
+        if (present(sphere_radius)) then
+            sphere_radius = self % domain_ptr % sphere_radius
+        end if
 
         nullify(maxedges_pointer)
         nullify(ncellssolve_pointer)

src/dynamics/mpas/dyn_grid_impl.F90

@@ -135,6 +135,7 @@ subroutine init_reference_pressure()
         use cam_history_support, only: add_vert_coord
         use cam_logfile, only: debugout_debug, debugout_verbose
         use dyn_comp, only: dyn_debug_print, mpas_dynamical_core
+        use dyn_procedures, only: reverse
         use dynconst, only: constant_p0 => pref
         use ref_pres, only: ref_pres_init
         use std_atm_profile, only: std_atm_pres
@@ -143,6 +144,8 @@ subroutine init_reference_pressure()
         ! Module(s) from CESM Share.
         use shr_kind_mod, only: kind_r8 => shr_kind_r8, &
                                 len_cx => shr_kind_cx
+        ! Module(s) from MPAS.
+        use dyn_mpas_procedures, only: dzw_of_rdzw, zu_of_dzw, zw_of_dzw
 
         character(*), parameter :: subname = 'dyn_grid::init_reference_pressure'
         character(len_cx) :: cerr
@@ -178,7 +181,7 @@ subroutine init_reference_pressure()
         call check_allocate(ierr, subname, 'dzw(pver)', &
             file='dyn_grid', line=__LINE__, errmsg=trim(adjustl(cerr)))
 
-        dzw(:) = 1.0_kind_r8 / real(rdzw(:), kind_r8)
+        dzw(:) = dzw_of_rdzw(real(rdzw, kind_r8))
 
         nullify(rdzw)
 
@@ -191,13 +194,8 @@ subroutine init_reference_pressure()
 
         ! In MPAS, zeta coordinates are stored in increasing order (i.e., bottom to top of atmosphere).
         ! In CAM-SIMA, however, index order is reversed (i.e., top to bottom of atmosphere).
-        ! Compute in reverse below.
-        zw(pverp) = 0.0_kind_r8
-
-        do k = pver, 1, -1
-            zw(k) = zw(k + 1) + dzw(pver - k + 1)
-            zu(k) = 0.5_kind_r8 * (zw(k + 1) + zw(k))
-        end do
+        zw(:) = reverse(zw_of_dzw(dzw))
+        zu(:) = reverse(zu_of_dzw(dzw))
 
         deallocate(dzw)