vdiff_apply_momentum Subroutine

   subroutine vdiff_apply_momentum(grid, this, ms, dt, kv_source, &
                                   tau_u, tau_v, lambda_bot_u, lambda_bot_v, rho0, &
                                   visc_rem_u, visc_rem_v, remnant_only, &
                                   kv_corner_source, kv_corner_prandtl, &
                                   lambda_top_u, lambda_top_v, cover_u, cover_v)
      !! Backward-Euler vertical viscosity applied to
      !! `ms%u_face_x_layer` and `ms%v_face_y_layer`.  Per-face
      !! `h_face` averaged from the two abutting cell columns.
      !!
      !! Diffusivity source: either a 3D cell-centred field from a
      !! closure (PP81 / KPP) passed as `kv_source` (shape (nx, ny,
      !! nz+1) — interface-located, k=1 bed, k=nz+1 surface) OR the
      !! scalar `this%K_v_momentum` broadcast into the
      !! `kv_scalar_buf` workspace.  No-op when neither path
      !! supplies a non-zero diffusivity.
      !!
      !! Corner viscosity add-on (`kv_corner_source`, optional): a
      !! CORNER-staggered interface viscosity `(nx+1, ny+1, nz+1)`
      !! (corner (i,j) = SW corner of cell (i,j); k=1 bed interface,
      !! k=nz+1 surface), scaled by `kv_corner_prandtl` and ADDED to
      !! the face viscosity as the 2-point average of the face's two
      !! END corners — u-face (i,j) reads corners (i,j)/(i,j+1), v-face
      !! (i,j) reads corners (i,j)/(i+1,j).  This is the vertex
      !! kappa-shear `Kv` seam (JHL08 vertex form): the corner
      !! viscosity reaches the momentum solve WITHOUT passing through a
      !! tracer point, so it is NOT smoothed by a corner->centre->face
      !! round trip.  The add lands BEFORE the BBL-glue transform (the
      !! reference formulation folds all viscosity contributions into
      !! `Kv_tot` ahead of the coupling-coefficient bottom-boundary
      !! blend).  Absent ⇒ bit-identical.
      !!
      !! Implicit surface-stress / bottom-drag folding (optional, gated
      !! on `this%implicit_stress` / `this%implicit_drag`):
      !!   * `tau_u(nu, nv)` / `tau_v(nu, nv)` — wind stress (N/m²) on the
      !!     respective faces.  Divided by `rho0` and added to the surface
      !!     (`k = nz`) RHS row.  Only read when `implicit_stress`.
      !!   * `lambda_bot_u` / `lambda_bot_v` — bottom-drag Rayleigh rate
      !!     `λ` (1/s) on the respective faces, added to the bed (`k = 1`)
      !!     diagonal as `dt·λ`.  Only read when `implicit_drag`.
      !!   * `lambda_top_u` / `lambda_top_v` + `cover_u` / `cover_v` —
      !!     the ICE-SHELF TOP-drag twin (`implicit_top_drag`): the rate
      !!     is added to the SURFACE (`k = nz`) diagonal as `dt·λ`, and
      !!     the wind-stress RHS on that same row is scaled by
      !!     `(1 − cover)` so an ice-covered face takes the drag and NOT
      !!     the wind.  All four are read only when `implicit_top_drag`
      !!     AND all four are supplied.
      !! All optional + device-resident; absent ⇒ the matrix is built
      !! exactly as before (bit-identical).
      !!
      !! `visc_rem_u(nu, nv, nz)` / `visc_rem_v(nu, nv, nz)` — optional
      !! viscous-remnant γ_k output (the `bt_work%visc_rem_u/v` seam,
      !! PR-19).  Both present ⇒ the SAME already-factorized tridiagonal
      !! solved a second time with RHS ≡ 1, written into these arrays
      !! (see `diffuse_velocity_columns_impl`).  Absent ⇒ no remnant
      !! work is done and the momentum answer is unperturbed
      !! (bit-identical) — this is a pure add-on, not a new physics
      !! path.
      type(hgrid_t), intent(in) :: grid
      type(ocean_vdiff_t), intent(inout) :: this
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      ! EXPLICIT SHAPE — not `(:, :)` — on every optional array below, and
      ! the same on the two `rdb_ocean_dyn` frames that FORWARD them
      ! (`vmix_apply_in_stage`, `visc_rem_precompute`).  Load-bearing, not
      ! decoration.  Each of these is handed straight on to the
      ! EXPLICIT-SHAPE `(nu, nv)` / `(nx, ny, nz)` dummies of
      ! `diffuse_velocity_columns_impl`, i.e. sequence-associated.  An
      ! ASSUMED-SHAPE actual in that position forces the compiler to decide,
      ! at the call, whether the descriptor is contiguous or has to be
      ! packed into a temporary — and gfortran 15 computes that decision
      ! UNCONDITIONALLY, from descriptor extents/strides it only loaded
      ! under the `present()` guard.  An ABSENT optional therefore leaves
      ! the packing flag uninitialised, and the post-call cleanup
      ! (`if (.not. packed .and. ptr /= NULL) free(ptr)`) branches on it:
      ! 160 valgrind "conditional jump depends on uninitialised value(s)"
      ! from 10 contexts on `test_driver_ocean`, all rooted in this
      ! routine's stack frame.  Benign (the null check still gates the free)
      ! but real, and it is the compiler reading OUR uninitialised stack.
      !
      ! Declaring the shape removes the decision at its source: an
      ! explicit-shape dummy is contiguous by definition, so there is no
      ! flag, no temporary, no free, nothing to read — and, unlike the
      ! CONTIGUOUS attribute (tried and reverted during the cavity stack), it
      ! keeps working on nvfortran.  CONTIGUOUS on an assumed-shape OPTIONAL
      ! makes nvfortran
      ! 26.5 hand the callee an offset-adjusted, NON-NULL address for an
      ! ABSENT actual; the implicit present-or-copyin the compiler emits for
      ! the impl's explicit-shape dummy then loses its null short-circuit
      ! and issues a real `cuMemcpyHtoD` off that garbage pointer
      ! (`__pgi_uacc_dataonb ... 'tau_face(:,:)'` -> SIGSEGV, 44/220 tests
      ! on `-stdpar=gpu -gpu=cc70,mem:separate`).  An explicit-shape dummy
      ! carries no descriptor at all, so an absent actual stays a clean
      ! null and the runtime skips it, exactly as it does today.
      !
      ! The bounds are the ones the body already derives and the impl
      ! already assumes: `multilayer_state_t%init` allocates the u-face
      ! arrays `(nx_total + 1, ny_total, nz_ml)` and the v-face arrays
      ! `(nx_total, ny_total + 1, nz_ml)` off this same `grid`.  Every
      ! actual is a whole allocatable slot component or a whole local array,
      ! never a strided section (the GPU path already requires that: these
      ! arrays are device-mapped whole).  Bit-identical: the packing path
      ! was never taken, only speculated about.
      real(wp), intent(in), optional :: kv_source(grid%nx_total, grid%ny_total, ms%nz_ml + 1)
      real(wp), intent(in), optional :: tau_u(grid%nx_total + 1, grid%ny_total)
      real(wp), intent(in), optional :: tau_v(grid%nx_total, grid%ny_total + 1)
      real(wp), intent(in), optional :: lambda_bot_u(grid%nx_total + 1, grid%ny_total)
      real(wp), intent(in), optional :: lambda_bot_v(grid%nx_total, grid%ny_total + 1)
      real(wp), intent(in), optional :: lambda_top_u(grid%nx_total + 1, grid%ny_total)
      real(wp), intent(in), optional :: lambda_top_v(grid%nx_total, grid%ny_total + 1)
      real(wp), intent(in), optional :: cover_u(grid%nx_total + 1, grid%ny_total)
      real(wp), intent(in), optional :: cover_v(grid%nx_total, grid%ny_total + 1)
      real(wp), intent(in), optional :: rho0
         !! Boussinesq reference density for the implicit surface-stress
         !! fold, `dt*(tau/rho0)/h_nz`.  READ ONLY when that fold is active
         !! (`do_stress` below); every other path never reaches it.  There
         !! is deliberately NO defaulted value: the one rho0 of record is
         !! `&ocean_ic_nml rho_0` -> `eos%rho0`, which the production call
         !! site hands over as `surface_stress%rho0`
         !! (`configure_ocean_reference_density`).  Omitting it WITH the
         !! fold active is a wiring bug and fails loud — it used to fall
         !! back to a silent literal 1035, which is how a run configured at
         !! a different rho_0 could fold its wind stress in on the wrong
         !! one.
      real(wp), intent(inout), optional :: visc_rem_u(grid%nx_total + 1, grid%ny_total, ms%nz_ml)
      real(wp), intent(inout), optional :: visc_rem_v(grid%nx_total, grid%ny_total + 1, ms%nz_ml)
         !! Explicit-shape for the reason spelled out on the argument list
         !! above.
      logical, intent(in), optional :: remnant_only
         !! `.true.` = build the matrices and (re)fill `visc_rem_u/v`
         !! WITHOUT solving for / modifying the velocities — the
         !! pre-substep visc_rem refresh (PGF_BUG.md §9; MOM6 computes
         !! `vertvisc_coef` before `btstep`, so its BT weights never lag).
         !! Requires both `visc_rem_u/v` present.  Default `.false.`.
      real(wp), intent(in), optional :: kv_corner_source(grid%nx_total + 1, grid%ny_total + 1, ms%nz_ml + 1)
         !! Corner-staggered interface viscosity source, `(nx+1, ny+1,
         !! nz+1)` — see the corner add-on note above.  Typically the
         !! vertex kappa-shear `kd_corner` carrier (device-resident).
         !! Explicit-shape for the reason spelled out on the argument list
         !! above.
      real(wp), intent(in), optional :: kv_corner_prandtl
         !! Scale applied to `kv_corner_source` (Kv = Pr·Kd; the vertex
         !! kappa-shear `prandtl_turb`).  Default 1.

      integer :: nx, ny, nx_face, ny_uface, nx_vface, ny_face, nz
      logical :: do_stress, do_drag, do_remnant, solve_mom, do_corner, do_top
      real(wp) :: rho0_l, corner_prandtl_l

      nx = grid%nx_total
      ny = grid%ny_total
      nx_face = size(ms%u_face_x_layer, 1)
      ny_uface = size(ms%u_face_x_layer, 2)
      nx_vface = size(ms%v_face_y_layer, 1)
      ny_face = size(ms%v_face_y_layer, 2)
      nz = ms%nz_ml

      ! Only fold the boundary BCs when the slot knob is on AND the caller
      ! supplied the corresponding face field (driver wires both together).
      do_stress = this%implicit_stress .and. present(tau_u) .and. present(tau_v)
      do_drag = this%implicit_drag .and. present(lambda_bot_u) .and. present(lambda_bot_v)
      do_top = this%implicit_top_drag .and. present(lambda_top_u) .and. &
               present(lambda_top_v) .and. present(cover_u) .and. present(cover_v)
      do_remnant = present(visc_rem_u) .and. present(visc_rem_v)
      solve_mom = .true.
      if (present(remnant_only)) solve_mom = .not. remnant_only
      if (.not. solve_mom .and. .not. do_remnant) return   ! nothing to produce
      ! Reference density for the implicit stress fold.  The placeholder is
      ! inert — `rho0_l` only reaches an answer under `do_stress`, and that
      ! branch is guarded below — so it must never be a plausible-looking
      ! seawater density that a caller could silently inherit.
      rho0_l = 1.0_wp
      if (present(rho0)) then
         rho0_l = rho0
      else if (do_stress) then
         error stop "vdiff_apply_momentum: implicit_stress is on and tau_u/tau_v "// &
            "were supplied, but rho0 was not — pass the configured reference "// &
            "density (surface_stress%rho0, from &ocean_ic_nml rho_0)"
      end if
      do_corner = present(kv_corner_source)
      corner_prandtl_l = 1.0_wp
      if (present(kv_corner_prandtl)) corner_prandtl_l = kv_corner_prandtl
      ! Glue without the per-face MOM6 BBL (`bbl_per_face`, set at
      ! configure): the historical constant piston, `kv_bbl =
      ! bbl_piston·hbbl_visc` over `bbl_thick = hbbl_visc`, refilled every
      ! call so a hand-built slot that sets the knobs after `init` sees them.
      if (this%bbl_glue .and. .not. this%bbl_per_face) then
         call fill_bbl_constants(this%kv_bbl_u, this%bbl_thick_u, &
                                 this%bbl_piston*this%hbbl_visc, this%hbbl_visc, nx_face, ny_uface)
         call fill_bbl_constants(this%kv_bbl_v, this%bbl_thick_v, &
                                 this%bbl_piston*this%hbbl_visc, this%hbbl_visc, nx_vface, ny_face)
      end if

      if (present(kv_source)) then
         call diffuse_velocity_columns_impl( &
            nx_face, ny_uface, nz, dt, &
            ms%u_face_x_layer, ms%h_layer, kv_source, ms%wet_mask, &
            .true., nx, ny, this%use_harmonic, &
            this%a_diag_u%data, this%b_diag_u%data, &
            this%c_diag_u%data, this%rhs_u%data, &
            do_stress, do_drag, rho0_l, tau_u, lambda_bot_u, &
            do_top, lambda_top_u, cover_u, &
            solve_mom, do_remnant, visc_rem_u, &
            this%hvel_mom6, this%hbbl_visc, &
            this%bbl_glue, this%hvel_upwind, this%hvel_harmonic, &
            this%kv_bbl_u, this%bbl_thick_u, this%kv_bbl_bg, &
            do_corner, corner_prandtl_l, kv_corner_source, &
            this%zlevel_faces, ms%k_top_u, ms%k_bot_u)
         call diffuse_velocity_columns_impl( &
            nx_vface, ny_face, nz, dt, &
            ms%v_face_y_layer, ms%h_layer, kv_source, ms%wet_mask, &
            .false., nx, ny, this%use_harmonic, &
            this%a_diag_v%data, this%b_diag_v%data, &
            this%c_diag_v%data, this%rhs_v%data, &
            do_stress, do_drag, rho0_l, tau_v, lambda_bot_v, &
            do_top, lambda_top_v, cover_v, &
            solve_mom, do_remnant, visc_rem_v, &
            this%hvel_mom6, this%hbbl_visc, &
            this%bbl_glue, this%hvel_upwind, this%hvel_harmonic, &
            this%kv_bbl_v, this%bbl_thick_v, this%kv_bbl_bg, &
            do_corner, corner_prandtl_l, kv_corner_source, &
            this%zlevel_faces, ms%k_top_v, ms%k_bot_v)
      else
         ! Pure-vdiff no-op short-circuit ONLY when there is also no
         ! boundary forcing to fold; stress/drag/remnant must still be
         ! applied even at K_v = 0 (the surface Ekman + bed sink don't
         ! need interior viscosity to act, and an implicit_drag-only,
         ! K_v=0 config would otherwise silently leave visc_rem stale).
         ! A corner viscosity source likewise keeps the solve alive.
         if (this%K_v_momentum <= 0.0_wp .and. .not. do_stress .and. .not. do_drag &
             .and. .not. do_top .and. .not. do_remnant .and. .not. do_corner &
             .and. .not. this%bbl_glue) return
         call fill_kv_scalar_buf(this%kv_scalar_buf%data, &
                                 this%K_v_momentum, nx, ny, nz)
         call diffuse_velocity_columns_impl( &
            nx_face, ny_uface, nz, dt, &
            ms%u_face_x_layer, ms%h_layer, this%kv_scalar_buf%data, ms%wet_mask, &
            .true., nx, ny, this%use_harmonic, &
            this%a_diag_u%data, this%b_diag_u%data, &
            this%c_diag_u%data, this%rhs_u%data, &
            do_stress, do_drag, rho0_l, tau_u, lambda_bot_u, &
            do_top, lambda_top_u, cover_u, &
            solve_mom, do_remnant, visc_rem_u, &
            this%hvel_mom6, this%hbbl_visc, &
            this%bbl_glue, this%hvel_upwind, this%hvel_harmonic, &
            this%kv_bbl_u, this%bbl_thick_u, this%kv_bbl_bg, &
            do_corner, corner_prandtl_l, kv_corner_source, &
            this%zlevel_faces, ms%k_top_u, ms%k_bot_u)
         call diffuse_velocity_columns_impl( &
            nx_vface, ny_face, nz, dt, &
            ms%v_face_y_layer, ms%h_layer, this%kv_scalar_buf%data, ms%wet_mask, &
            .false., nx, ny, this%use_harmonic, &
            this%a_diag_v%data, this%b_diag_v%data, &
            this%c_diag_v%data, this%rhs_v%data, &
            do_stress, do_drag, rho0_l, tau_v, lambda_bot_v, &
            do_top, lambda_top_v, cover_v, &
            solve_mom, do_remnant, visc_rem_v, &
            this%hvel_mom6, this%hbbl_visc, &
            this%bbl_glue, this%hvel_upwind, this%hvel_harmonic, &
            this%kv_bbl_v, this%bbl_thick_v, this%kv_bbl_bg, &
            do_corner, corner_prandtl_l, kv_corner_source, &
            this%zlevel_faces, ms%k_top_v, ms%k_bot_v)
      end if
   end subroutine vdiff_apply_momentum