vmix_compute_pp81 Subroutine

public pure subroutine vmix_compute_pp81(grid, this, ms)

Pacanowski-Philander (1981) Richardson-number closure. Inlined kernel — keeps the do concurrent body adjacent to its derived-type accesses. We tried the outer-shim + _impl pattern but NVHPC’s stdpar codegen produced more descriptor-marshalling memcpys at the shim boundary than the direct-access form generates inside the kernel. Direct ms%foo(i,j,k) access is what other working hot-path kernels (continuity, coriolis_adv, the barotropic substep) use; the shim pattern was an experiment that didn’t help here.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_vmix_t), intent(inout) :: this
type(multilayer_state_t), intent(in) :: ms

Calls

proc~~vmix_compute_pp81~~CallsGraph proc~vmix_compute_pp81 vmix_compute_pp81 local local proc~vmix_compute_pp81->local

Called by

proc~~vmix_compute_pp81~~CalledByGraph proc~vmix_compute_pp81 vmix_compute_pp81 proc~vmix_apply_in_stage vmix_apply_in_stage proc~vmix_apply_in_stage->proc~vmix_compute_pp81 proc~run_stage run_stage proc~run_stage->proc~vmix_apply_in_stage proc~run_stage_split run_stage_split proc~run_stage_split->proc~vmix_apply_in_stage proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~run_stage proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~engine_step->proc~ocean_dyn_step_split proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: denom
real(kind=wp), private :: du_dz
real(kind=wp), private :: dv_dz
real(kind=wp), private :: dz_face
integer, private :: i
integer, private :: j
integer, private :: k
real(kind=wp), private :: n2
integer, private :: nx
integer, private :: ny
integer, private :: nz
real(kind=wp), private :: ri
real(kind=wp), private :: ri_factor
real(kind=wp), private :: shear2
real(kind=wp), private :: u_k
real(kind=wp), private :: u_km1
real(kind=wp), private :: v_k
real(kind=wp), private :: v_km1

Source Code

   pure subroutine vmix_compute_pp81(grid, this, ms)
      !! Pacanowski-Philander (1981) Richardson-number closure.
      !! Inlined kernel — keeps the `do concurrent` body adjacent to
      !! its derived-type accesses.  We tried the outer-shim +
      !! `_impl` pattern but NVHPC's stdpar codegen produced more
      !! descriptor-marshalling memcpys at the shim boundary than the
      !! direct-access form generates inside the kernel.  Direct
      !! `ms%foo(i,j,k)` access is what other working hot-path
      !! kernels (continuity, coriolis_adv, the barotropic substep) use; the
      !! shim pattern was an experiment that didn't help here.
      type(hgrid_t), intent(in) :: grid
      type(ocean_vmix_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms

      integer :: i, j, k, nx, ny, nz
      real(wp) :: u_km1, v_km1, u_k, v_k
      real(wp) :: du_dz, dv_dz, shear2, dz_face
      real(wp) :: n2, ri, denom, ri_factor

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml

      ! Boundary interfaces: keep them at zero (closed top + bottom).
      do concurrent(j=1:ny, i=1:nx)
         this%kv(i, j, 1) = 0.0_wp
         this%kt(i, j, 1) = 0.0_wp
         this%kv(i, j, nz + 1) = 0.0_wp
         this%kt(i, j, nz + 1) = 0.0_wp
      end do

      ! Interior interfaces k = 2..nz.  Richardson stability factor
      ! in one pass per (i, j, k).  Loop spans the full (i, j) extent
      ! including walls — face arrays are (nx+1, ny, *) and
      ! (nx, ny+1, *), so i+1 / j+1 reach valid storage.  At closed
      ! walls the face velocities are zero ⇒ shear² floors ⇒ Ri huge
      ! ⇒ same kv as adjacent interior.  A wall-only fallback that
      ! sets kv = pp81_nu_bg breaks horizontal symmetry by O(nu0/denom²)
      ! and seeds a baroclinic-noise instability in stratified
      ! closed basins.
      do concurrent(k=2:nz, j=1:ny, i=1:nx) &
         local(u_km1, v_km1, u_k, v_k, du_dz, dv_dz, &
               shear2, dz_face, n2, ri, denom, ri_factor)
         u_km1 = 0.5_wp*(ms%u_face_x_layer(i, j, k - 1) + ms%u_face_x_layer(i + 1, j, k - 1))
         v_km1 = 0.5_wp*(ms%v_face_y_layer(i, j, k - 1) + ms%v_face_y_layer(i, j + 1, k - 1))
         u_k = 0.5_wp*(ms%u_face_x_layer(i, j, k) + ms%u_face_x_layer(i + 1, j, k))
         v_k = 0.5_wp*(ms%v_face_y_layer(i, j, k) + ms%v_face_y_layer(i, j + 1, k))

         dz_face = 0.5_wp*(ms%h_layer(i, j, k - 1) + ms%h_layer(i, j, k))
         if (dz_face <= 0.0_wp) then
            this%kv(i, j, k) = this%pp81_nu_bg
            this%kt(i, j, k) = this%pp81_kappa_bg
            cycle
         end if

         du_dz = (u_k - u_km1)/dz_face
         dv_dz = (v_k - v_km1)/dz_face
         shear2 = max(du_dz*du_dz + dv_dz*dv_dz, this%shear2_floor)

         n2 = -GRAVITY*(ms%rho_layer(i, j, k) - ms%rho_layer(i, j, k - 1))/ &
              (this%rho0*dz_face)
         ri = n2/shear2
         denom = 1.0_wp + this%pp81_alpha*max(ri, 0.0_wp)
         ri_factor = 1.0_wp/(denom*denom)
         this%kv(i, j, k) = this%pp81_nu_bg + this%pp81_nu0*ri_factor
         this%kt(i, j, k) = this%pp81_kappa_bg + this%pp81_nu0*ri_factor/denom
      end do
   end subroutine vmix_compute_pp81