gm_pe_release Subroutine

private pure subroutine gm_pe_release(nx, ny, nz, slope_max, rho0, h_layer, slope_x, slope_y, n2_u, n2_v, khth_u, khth_v, gm_src)

GM potential-energy release at cell centres for the MEKE seam: gm_src = 1/4 * Sum_k rho0 * (KHSlope^2N^2) * h over the four straddling faces, summed over interior interfaces; slope clamped to slope_max, N^2 floored at 0. gm_src >= 0 for a stable tilted column. An interface value is attributed to the layer below it (kb=K-1); bed + surface carry zero slope/N^2.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: slope_max
real(kind=wp), intent(in) :: rho0
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(in) :: slope_x(nx+1,ny,nz+1)
real(kind=wp), intent(in) :: slope_y(nx,ny+1,nz+1)
real(kind=wp), intent(in) :: n2_u(nx+1,ny,nz+1)
real(kind=wp), intent(in) :: n2_v(nx,ny+1,nz+1)
real(kind=wp), intent(in) :: khth_u(nx+1,ny)
real(kind=wp), intent(in) :: khth_v(nx,ny+1)
real(kind=wp), intent(out) :: gm_src(nx,ny)

Calls

proc~~gm_pe_release~~CallsGraph proc~gm_pe_release gm_pe_release local local proc~gm_pe_release->local proc~gm_clamp_slope gm_clamp_slope proc~gm_pe_release->proc~gm_clamp_slope proc~gm_pos_n2 gm_pos_n2 proc~gm_pe_release->proc~gm_pos_n2

Called by

proc~~gm_pe_release~~CalledByGraph proc~gm_pe_release gm_pe_release proc~gm_compute_impl gm_compute_impl proc~gm_compute_impl->proc~gm_pe_release proc~gm_compute_transports gm_compute_transports proc~gm_compute_transports->proc~gm_compute_impl proc~run_gm_step run_gm_step proc~run_gm_step->proc~gm_compute_transports proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_gm_step proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step_split

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: acc
real(kind=wp), private :: fsum
integer, private :: i
integer, private :: j
integer, private :: k
integer, private :: kb
real(kind=wp), private :: n2e
real(kind=wp), private :: n2n
real(kind=wp), private :: n2s
real(kind=wp), private :: n2w
real(kind=wp), private :: sx_e
real(kind=wp), private :: sx_w
real(kind=wp), private :: sy_n
real(kind=wp), private :: sy_s

Source Code

   pure subroutine gm_pe_release(nx, ny, nz, slope_max, rho0, h_layer, &
                                 slope_x, slope_y, n2_u, n2_v, &
                                 khth_u, khth_v, gm_src)
      !! GM potential-energy release at cell centres for the MEKE seam:
      !!   gm_src = 1/4 * Sum_k rho0 * (KH*Slope^2*N^2) * h
      !! over the four straddling faces, summed over interior interfaces;
      !! slope clamped to `slope_max`, N^2 floored at 0.  `gm_src >= 0` for
      !! a stable tilted column.  An interface value is attributed to the
      !! layer below it (kb=K-1); bed + surface carry zero slope/N^2.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: slope_max, rho0
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: slope_x(nx + 1, ny, nz + 1)
      real(wp), intent(in) :: slope_y(nx, ny + 1, nz + 1)
      real(wp), intent(in) :: n2_u(nx + 1, ny, nz + 1)
      real(wp), intent(in) :: n2_v(nx, ny + 1, nz + 1)
      real(wp), intent(in) :: khth_u(nx + 1, ny)
      real(wp), intent(in) :: khth_v(nx, ny + 1)
      real(wp), intent(out) :: gm_src(nx, ny)

      integer :: i, j, k, kb
      real(wp) :: acc, fsum, sx_w, sx_e, sy_s, sy_n
      real(wp) :: n2w, n2e, n2s, n2n

      do concurrent(j=1:ny, i=1:nx) &
         local(k, kb, acc, fsum, sx_w, sx_e, sy_s, sy_n, n2w, n2e, n2s, n2n)
         acc = 0.0_wp
         do k = 2, nz                 ! interior interface index Kr
            kb = k - 1                 ! layer below the interface
            sx_w = gm_clamp_slope(slope_x(i, j, k), slope_max)
            sx_e = gm_clamp_slope(slope_x(i + 1, j, k), slope_max)
            sy_s = gm_clamp_slope(slope_y(i, j, k), slope_max)
            sy_n = gm_clamp_slope(slope_y(i, j + 1, k), slope_max)
            n2w = gm_pos_n2(n2_u(i, j, k))
            n2e = gm_pos_n2(n2_u(i + 1, j, k))
            n2s = gm_pos_n2(n2_v(i, j, k))
            n2n = gm_pos_n2(n2_v(i, j + 1, k))
            fsum = (khth_u(i, j)*sx_w*sx_w*n2w + khth_u(i + 1, j)*sx_e*sx_e*n2e) + &
                   (khth_v(i, j)*sy_s*sy_s*n2s + khth_v(i, j + 1)*sy_n*sy_n*n2n)
            acc = acc + fsum*h_layer(i, j, kb)
         end do
         gm_src(i, j) = 0.25_wp*rho0*acc
      end do
   end subroutine gm_pe_release