vmix_add_kv_ml_invz2 Subroutine

public pure subroutine vmix_add_kv_ml_invz2(grid, this, ms)

Augment this%kv with an extra near-surface viscosity:

kv_extra(z) = kv_ml_invz2 · (hmix_fixed / max(z, dz_min))²

where z is the depth (m) of the layer interface below the free surface. Active only for interfaces whose z < hmix_fixed. Adds to the existing kv (which kt does not see — momentum-only).

Interface convention: kv(:, :, k) lives at the bottom of layer k. Surface-most interior interface is at k = nz (between layer nz and the wall at nz+1). Surface boundary k = nz+1 stays at zero (closed top), and we don’t touch it.

When kv_ml_invz2 <= 0 the routine is a no-op so existing configs are bit-identical.

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_add_kv_ml_invz2~~CallsGraph proc~vmix_add_kv_ml_invz2 vmix_add_kv_ml_invz2 local local proc~vmix_add_kv_ml_invz2->local

Called by

proc~~vmix_add_kv_ml_invz2~~CalledByGraph proc~vmix_add_kv_ml_invz2 vmix_add_kv_ml_invz2 proc~vmix_apply_in_stage vmix_apply_in_stage proc~vmix_apply_in_stage->proc~vmix_add_kv_ml_invz2 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 :: dz_min_safe
real(kind=wp), private :: hmix
real(kind=wp), private :: hmix_sq
integer, private :: i
integer, private :: j
integer, private :: k
real(kind=wp), private :: kv_extra
integer, private :: nx
integer, private :: ny
integer, private :: nz
real(kind=wp), private :: z

Source Code

   pure subroutine vmix_add_kv_ml_invz2(grid, this, ms)
      !! Augment `this%kv` with an extra near-surface viscosity:
      !!
      !!     kv_extra(z) = kv_ml_invz2 · (hmix_fixed / max(z, dz_min))²
      !!
      !! where `z` is the depth (m) of the layer interface below the
      !! free surface.  Active only for interfaces whose `z <
      !! hmix_fixed`.  Adds to the existing `kv` (which `kt` does not
      !! see — momentum-only).
      !!
      !! Interface convention: `kv(:, :, k)` lives at the bottom of
      !! layer k.  Surface-most interior interface is at `k = nz`
      !! (between layer nz and the wall at nz+1).  Surface boundary
      !! `k = nz+1` stays at zero (closed top), and we don't touch it.
      !!
      !! When `kv_ml_invz2 <= 0` the routine is a no-op so existing
      !! configs are bit-identical.
      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) :: kv_extra, hmix, z, dz_min_safe, hmix_sq

      if (this%kv_ml_invz2 <= 0.0_wp) return
      if (this%hmix_fixed <= 0.0_wp) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      hmix = this%hmix_fixed
      hmix_sq = hmix*hmix

      ! Floor on `z` to prevent the 1/z² profile from diverging right
      ! under the surface.  Half a target layer thickness is a safe,
      ! grid-resolved minimum.
      dz_min_safe = 0.5_wp*hmix/real(max(nz, 1), wp)

      ! Interior interfaces: k = nz (surface-most) down to k = 2.
      ! The depth of interface k below the surface is the cumulative
      ! thickness of layers nz, nz-1, ..., k.
      do concurrent(j=1:ny, i=1:nx) local(k, z, kv_extra)
         z = 0.0_wp
         do k = nz, 2, -1
            z = z + ms%h_layer(i, j, k)
            if (z >= hmix) exit
            kv_extra = this%kv_ml_invz2*hmix_sq/(max(z, dz_min_safe)**2)
            this%kv(i, j, k) = this%kv(i, j, k) + kv_extra
         end do
      end do
   end subroutine vmix_add_kv_ml_invz2