surfstress_distributed_impl Subroutine

private pure subroutine surfstress_distributed_impl(h_layer, wet_mask, tau_x, tau_y, du_stress, dv_stress, rho0, h_min, hmix_stress, nx, ny, nz)

DIRECT_STRESS branch — distribute the wind stress across the top hmix_stress metres of the column. For each face, walk layers from k = nz (surface) down to k = 1, accumulating thickness; the surface boundary layer (SBL) is the set of layers whose top sits within hmix_stress of the free surface. The acceleration per layer is:

du_k/dt = (tau / (rho0 · hmix_stress)) · (h_in_sbl_k / h_face_k)

Sum over k: total impulse per unit area = tau / rho0, matching the bed-only formulation. When hmix_stress is smaller than the top layer thickness the SBL is contained in the surface-most layer and the formula collapses to tau / (rho0 · h_face_nz) — bit-identical to the existing kernel. When hmix_stress straddles multiple layers, each gets its fractional acceleration.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(in) :: wet_mask(nx,ny)
real(kind=wp), intent(in) :: tau_x(nx+1,ny)
real(kind=wp), intent(in) :: tau_y(nx,ny+1)
real(kind=wp), intent(inout) :: du_stress(nx+1,ny,nz)
real(kind=wp), intent(inout) :: dv_stress(nx,ny+1,nz)
real(kind=wp), intent(in) :: rho0
real(kind=wp), intent(in) :: h_min
real(kind=wp), intent(in) :: hmix_stress
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz

Calls

proc~~surfstress_distributed_impl~~CallsGraph proc~surfstress_distributed_impl surfstress_distributed_impl local local proc~surfstress_distributed_impl->local

Called by

proc~~surfstress_distributed_impl~~CalledByGraph proc~surfstress_distributed_impl surfstress_distributed_impl proc~ocean_surface_stress_compute_tendencies ocean_surface_stress_compute_tendencies proc~ocean_surface_stress_compute_tendencies->proc~surfstress_distributed_impl proc~run_stage run_stage proc~run_stage->proc~ocean_surface_stress_compute_tendencies proc~run_stage_split run_stage_split proc~run_stage_split->proc~ocean_surface_stress_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 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 :: cumul_h
real(kind=wp), private :: h_face_k
real(kind=wp), private :: h_in_sbl
integer, private :: i
real(kind=wp), private :: inv_rho_hmix
integer, private :: j
integer, private :: k
real(kind=wp), private :: mask_face

Source Code

   pure subroutine surfstress_distributed_impl(h_layer, wet_mask, tau_x, tau_y, &
                                               du_stress, dv_stress, &
                                               rho0, h_min, hmix_stress, &
                                               nx, ny, nz)
      !! DIRECT_STRESS branch — distribute the wind stress across
      !! the top `hmix_stress` metres of the column.  For each face,
      !! walk layers from k = nz (surface) down to k = 1, accumulating
      !! thickness; the surface boundary layer (SBL) is the set of
      !! layers whose top sits within `hmix_stress` of the free
      !! surface.  The acceleration per layer is:
      !!
      !!     du_k/dt = (tau / (rho0 · hmix_stress)) · (h_in_sbl_k / h_face_k)
      !!
      !! Sum over k: total impulse per unit area = tau / rho0,
      !! matching the bed-only formulation.  When `hmix_stress` is
      !! smaller than the top layer thickness the SBL is contained in
      !! the surface-most layer and the formula collapses to
      !! `tau / (rho0 · h_face_nz)` — bit-identical to the existing
      !! kernel.  When `hmix_stress` straddles multiple layers, each
      !! gets its fractional acceleration.
      integer, intent(in)    :: nx, ny, nz
      real(wp), intent(in)    :: h_layer(nx, ny, nz), wet_mask(nx, ny)
      real(wp), intent(in)    :: tau_x(nx + 1, ny), tau_y(nx, ny + 1)
      real(wp), intent(inout) :: du_stress(nx + 1, ny, nz), dv_stress(nx, ny + 1, nz)
      real(wp), intent(in)    :: rho0, h_min, hmix_stress
      integer :: i, j, k
      real(wp) :: inv_rho_hmix, cumul_h, h_face_k, h_in_sbl, mask_face

      inv_rho_hmix = 1.0_wp/(rho0*hmix_stress)

      do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
         du_stress(i, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, j=1:ny + 1, i=1:nx)
         dv_stress(i, j, k) = 0.0_wp
      end do

      do concurrent(j=1:ny, i=2:nx) &
         local(k, cumul_h, h_face_k, h_in_sbl, mask_face)
         mask_face = min(wet_mask(i - 1, j), wet_mask(i, j))
         cumul_h = 0.0_wp
         do k = nz, 1, -1
            if (cumul_h >= hmix_stress) exit
            h_face_k = 0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k))
            if (h_face_k <= 0.0_wp) exit
            h_in_sbl = max(0.0_wp, min(h_face_k, hmix_stress - cumul_h))
            du_stress(i, j, k) = mask_face* &
                                 tau_x(i, j)*inv_rho_hmix* &
                                 (h_in_sbl/max(h_face_k, h_min))
            cumul_h = cumul_h + h_face_k
         end do
      end do
      do concurrent(j=2:ny, i=1:nx) &
         local(k, cumul_h, h_face_k, h_in_sbl, mask_face)
         mask_face = min(wet_mask(i, j - 1), wet_mask(i, j))
         cumul_h = 0.0_wp
         do k = nz, 1, -1
            if (cumul_h >= hmix_stress) exit
            h_face_k = 0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k))
            if (h_face_k <= 0.0_wp) exit
            h_in_sbl = max(0.0_wp, min(h_face_k, hmix_stress - cumul_h))
            dv_stress(i, j, k) = mask_face* &
                                 tau_y(i, j)*inv_rho_hmix* &
                                 (h_in_sbl/max(h_face_k, h_min))
            cumul_h = cumul_h + h_face_k
         end do
      end do
   end subroutine surfstress_distributed_impl