hvisc_clamp_A Subroutine

private pure subroutine hvisc_clamp_A(ah_t, ah_q, bound_coef, dt, idxT, idyT, idxCu, idyCu, idxCv, idyCv, iareaCu, iareaCv, nx, ny, nz)

Per-cell CFL viscosity limiter (MOM6 BOUND_KH). Clamps the T-cell viscosity to Kh_Max_xx and the corner viscosity to Kh_Max_xy, each derived from the actual discrete stress stencil metrics + dt so the explicit forward-Euler viscous update can never overshoot. Replaces the global ah_max cap. On a uniform square grid Kh_Max = bound_coef·0.25/(dt·(1/dx²+ 1/dy²)).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: ah_t(nx,ny,nz)
real(kind=wp), intent(inout) :: ah_q(nx+1,ny+1,nz)
real(kind=wp), intent(in) :: bound_coef
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: idxT(nx,ny)
real(kind=wp), intent(in) :: idyT(nx,ny)
real(kind=wp), intent(in) :: idxCu(nx+1,ny)
real(kind=wp), intent(in) :: idyCu(nx+1,ny)
real(kind=wp), intent(in) :: idxCv(nx,ny+1)
real(kind=wp), intent(in) :: idyCv(nx,ny+1)
real(kind=wp), intent(in) :: iareaCu(nx+1,ny)
real(kind=wp), intent(in) :: iareaCv(nx,ny+1)
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz

Calls

proc~~hvisc_clamp_a~~CallsGraph proc~hvisc_clamp_a hvisc_clamp_A local local proc~hvisc_clamp_a->local

Called by

proc~~hvisc_clamp_a~~CalledByGraph proc~hvisc_clamp_a hvisc_clamp_A proc~ocean_horizontal_viscosity_compute_tendencies_on ocean_horizontal_viscosity_compute_tendencies_on proc~ocean_horizontal_viscosity_compute_tendencies_on->proc~hvisc_clamp_a proc~ocean_horizontal_viscosity_compute_tendencies ocean_horizontal_viscosity_compute_tendencies proc~ocean_horizontal_viscosity_compute_tendencies->proc~ocean_horizontal_viscosity_compute_tendencies_on proc~run_stage run_stage proc~run_stage->proc~ocean_horizontal_viscosity_compute_tendencies proc~run_stage_split run_stage_split proc~run_stage_split->proc~ocean_horizontal_viscosity_compute_tendencies 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

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: denom
real(kind=wp), private :: dx2
real(kind=wp), private :: dy2
integer, private :: i
real(kind=wp), private :: idt
integer, private :: j
integer, private :: k
real(kind=wp), private :: kh_max
real(kind=wp), private :: tx
real(kind=wp), private :: ty

Source Code

   pure subroutine hvisc_clamp_A(ah_t, ah_q, bound_coef, dt, &
                                 idxT, idyT, idxCu, idyCu, idxCv, idyCv, &
                                 iareaCu, iareaCv, nx, ny, nz)
      !! Per-cell CFL viscosity limiter (MOM6 `BOUND_KH`).  Clamps the
      !! T-cell viscosity to `Kh_Max_xx` and the corner viscosity to
      !! `Kh_Max_xy`, each derived from the actual discrete stress
      !! stencil metrics + `dt` so the explicit forward-Euler viscous
      !! update can never overshoot.  Replaces the global `ah_max` cap.
      !! On a uniform square grid `Kh_Max = bound_coef·0.25/(dt·(1/dx²+
      !! 1/dy²))`.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(inout) :: ah_t(nx, ny, nz), ah_q(nx + 1, ny + 1, nz)
      real(wp), intent(in) :: bound_coef, dt
      real(wp), intent(in) :: idxT(nx, ny), idyT(nx, ny)
      real(wp), intent(in) :: idxCu(nx + 1, ny), idyCu(nx + 1, ny)
      real(wp), intent(in) :: idxCv(nx, ny + 1), idyCv(nx, ny + 1)
      real(wp), intent(in) :: iareaCu(nx + 1, ny), iareaCv(nx, ny + 1)
      integer :: i, j, k
      real(wp) :: idt, kh_max, denom, tx, ty, dx2, dy2
      idt = 1.0_wp/dt
      ! T-cell bound: tension stress acts on the two u-faces (i, i+1)
      ! and two v-faces (j, j+1).  dx2/dy2 per cell from the metric
      ! inverses (dx2 = 1/idxT²).  Full T-cell range (incl. the edge
      ! rows whose stress feeds the interior divergence) — Cartesian /
      ! curvilinear metrics are valid in the ghosts, so the bound is
      ! defined everywhere; guard denom>0 for any zeroed land metric.
      do concurrent(k=1:nz, j=1:ny, i=1:nx) &
         local(kh_max, denom, tx, ty, dx2, dy2)
         dx2 = 1.0_wp/(idxT(i, j)*idxT(i, j))
         dy2 = 1.0_wp/(idyT(i, j)*idyT(i, j))
         tx = dy2*(idyT(i, j)/idxT(i, j))*(idyCu(i + 1, j) + idyCu(i, j))* &
              max(idyCu(i + 1, j)*iareaCu(i + 1, j), idyCu(i, j)*iareaCu(i, j))
         ty = dx2*(idxT(i, j)/idyT(i, j))*(idxCv(i, j + 1) + idxCv(i, j))* &
              max(idxCv(i, j + 1)*iareaCv(i, j + 1), idxCv(i, j)*iareaCv(i, j))
         denom = max(tx, ty)
         if (denom > 0.0_wp) then
            kh_max = bound_coef*0.25_wp*idt/denom
            if (ah_t(i, j, k) > kh_max) ah_t(i, j, k) = kh_max
         end if
      end do
      ! Corner bound: shear stress at Bu(i,j) acts on the u-faces
      ! (i, j-1)/(i, j) and v-faces (i-1, j)/(i, j).  Reuse the
      ! T-stencil metric magnitudes at the SW T-cell — the corner
      ! bound differs from the T-cell bound only by which faces it
      ! sums, and on a square grid both collapse to the same value.
      do concurrent(k=1:nz, j=2:ny, i=2:nx) &
         local(kh_max, denom, tx, ty, dx2, dy2)
         dx2 = 1.0_wp/(idxCu(i, j)*idxCu(i, j))
         dy2 = 1.0_wp/(idyCv(i, j)*idyCv(i, j))
         tx = dx2*(idxCu(i, j) + idxCu(i, j - 1))* &
              max(idxCu(i, j)*iareaCu(i, j), idxCu(i, j - 1)*iareaCu(i, j - 1))
         ty = dy2*(idyCv(i, j) + idyCv(i - 1, j))* &
              max(idyCv(i, j)*iareaCv(i, j), idyCv(i - 1, j)*iareaCv(i - 1, j))
         denom = max(tx, ty)
         if (denom > 0.0_wp) then
            kh_max = bound_coef*0.25_wp*idt/denom
            if (ah_q(i, j, k) > kh_max) ah_q(i, j, k) = kh_max
         end if
      end do
   end subroutine hvisc_clamp_A