compute_gprime_impl Subroutine

private pure subroutine compute_gprime_impl(h_layer, b, dpdx_face, dpdy_face, gfs, gint, idxCu, idyCv, nx, ny, nz)

Reduced-gravity / gprime PGF for NK = 2.

Convention: k = 1 bottom (heavier), k = nz = 2 surface (lighter). b(i, j) is bathymetric depth (positive-down), used to recover ∇η: sum_h = h_1 + h_2 = b + η; η = sum_h - b; ∇η = ∇(h_1+h_2) - ∇b. The ∇b subtraction matters: without it ∇H_bathy dominates ∇η by 10³–10⁴ on shelf-break/spoon configs, over-driving the gyre.

a_top = -g_FS · ∇η a_bot = -g_FS · ∇η - g’_int · ∇h_1

u-face gradient (f(i,j) - f(i-1,j)) · idxCu(i,j), mirror for v-face. Higher k stays zero ⇒ no-op on NK > 2.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(in) :: b(nx,ny)
real(kind=wp), intent(inout) :: dpdx_face(nx+1,ny,nz)
real(kind=wp), intent(inout) :: dpdy_face(nx,ny+1,nz)
real(kind=wp), intent(in) :: gfs
real(kind=wp), intent(in) :: gint
real(kind=wp), intent(in) :: idxCu(nx+1,ny)
real(kind=wp), intent(in) :: idyCv(nx,ny+1)
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz

Calls

proc~~compute_gprime_impl~~CallsGraph proc~compute_gprime_impl compute_gprime_impl local local proc~compute_gprime_impl->local

Called by

proc~~compute_gprime_impl~~CalledByGraph proc~compute_gprime_impl compute_gprime_impl proc~ocean_pressure_force_compute ocean_pressure_force_compute proc~ocean_pressure_force_compute->proc~compute_gprime_impl proc~run_stage run_stage proc~run_stage->proc~ocean_pressure_force_compute proc~run_stage_split run_stage_split proc~run_stage_split->proc~ocean_pressure_force_compute 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 :: db_dx
real(kind=wp), private :: db_dy
real(kind=wp), private :: dhbot_dx
real(kind=wp), private :: dhbot_dy
real(kind=wp), private :: dssh_dx
real(kind=wp), private :: dssh_dy
integer, private :: i
integer, private :: j
integer, private :: k
real(kind=wp), private :: sum_h_above
real(kind=wp), private :: sum_h_below
real(kind=wp), private :: sum_h_left
real(kind=wp), private :: sum_h_right

Source Code

   pure subroutine compute_gprime_impl(h_layer, b, dpdx_face, dpdy_face, &
                                       gfs, gint, idxCu, idyCv, nx, ny, nz)
      !! Reduced-gravity / gprime PGF for NK = 2.
      !!
      !! Convention: k = 1 bottom (heavier), k = nz = 2 surface (lighter).
      !! `b(i, j)` is bathymetric depth (positive-down), used to recover ∇η:
      !!   sum_h = h_1 + h_2 = b + η;  η = sum_h - b;  ∇η = ∇(h_1+h_2) - ∇b.
      !! The ∇b subtraction matters: without it ∇H_bathy dominates ∇η by
      !! 10³–10⁴ on shelf-break/spoon configs, over-driving the gyre.
      !!
      !!   a_top = -g_FS · ∇η
      !!   a_bot = -g_FS · ∇η - g'_int · ∇h_1
      !!
      !! u-face gradient `(f(i,j) - f(i-1,j)) · idxCu(i,j)`, mirror for
      !! v-face. Higher k stays zero ⇒ no-op on NK > 2.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: gfs, gint
      real(wp), intent(in)    :: idxCu(nx + 1, ny), idyCv(nx, ny + 1)
      real(wp), intent(in)    :: h_layer(nx, ny, nz)
      real(wp), intent(in)    :: b(nx, ny)
      real(wp), intent(inout) :: dpdx_face(nx + 1, ny, nz)
      real(wp), intent(inout) :: dpdy_face(nx, ny + 1, nz)
      integer :: i, j, k
      real(wp) :: sum_h_left, sum_h_right, sum_h_below, sum_h_above
      real(wp) :: db_dx, db_dy
      real(wp) :: dssh_dx, dssh_dy, dhbot_dx, dhbot_dy

      ! Zero everything first.  Then fill k = 1 and k = 2 explicitly.
      do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
         dpdx_face(i, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, j=1:ny + 1, i=1:nx)
         dpdy_face(i, j, k) = 0.0_wp
      end do

      if (nz < 2) return

      ! u-face (east): two adjacent cells (i-1, j) and (i, j).
      do concurrent(j=1:ny, i=2:nx) &
         local(sum_h_left, sum_h_right, db_dx, dssh_dx, dhbot_dx)
         sum_h_left = h_layer(i - 1, j, 1) + h_layer(i - 1, j, 2)
         sum_h_right = h_layer(i, j, 1) + h_layer(i, j, 2)
         db_dx = (b(i, j) - b(i - 1, j))*idxCu(i, j)
         dssh_dx = (sum_h_right - sum_h_left)*idxCu(i, j) - db_dx
         dhbot_dx = (h_layer(i, j, 1) - h_layer(i - 1, j, 1))*idxCu(i, j)
         dpdx_face(i, j, 2) = -gfs*dssh_dx                 ! top layer (k=nz=2)
         dpdx_face(i, j, 1) = -gfs*dssh_dx - gint*dhbot_dx  ! bottom layer (k=1)
      end do

      ! v-face (north): two adjacent cells (i, j-1) and (i, j).
      do concurrent(j=2:ny, i=1:nx) &
         local(sum_h_below, sum_h_above, db_dy, dssh_dy, dhbot_dy)
         sum_h_below = h_layer(i, j - 1, 1) + h_layer(i, j - 1, 2)
         sum_h_above = h_layer(i, j, 1) + h_layer(i, j, 2)
         db_dy = (b(i, j) - b(i, j - 1))*idyCv(i, j)
         dssh_dy = (sum_h_above - sum_h_below)*idyCv(i, j) - db_dy
         dhbot_dy = (h_layer(i, j, 1) - h_layer(i, j - 1, 1))*idyCv(i, j)
         dpdy_face(i, j, 2) = -gfs*dssh_dy
         dpdy_face(i, j, 1) = -gfs*dssh_dy - gint*dhbot_dy
      end do
   end subroutine compute_gprime_impl