mask_layer_velocities Subroutine

   pure subroutine mask_layer_velocities(grid, metrics, ms, bt_work)
      !! Zero the per-layer face velocities at land faces (spec §14 C4 /
      !! R4(b.2) / MOM6 `up = mask2dCu·(u+dt·accel)`).  Runs once per
      !! RK2 stage AFTER
      !! `apply_bt_correction` and inside the RK2 average so the
      !! additive layer tendencies + the BT correction cannot leave a
      !! re-ingested land-face velocity (a slow conservation leak).  The
      !! transports themselves already ride zeroed metrics; this resets
      !! the prognostic velocity so `derive_bt_from_layers` next step
      !! sees zero there.  All-wet ⇒ `wet_u/v≡1` ⇒ literal no-op.
      !!
      !! When `bt_work` is passed AND wet/dry is enabled
      !! (docs/ocean_wetdry_plan.md §4), the DYNAMIC bed-blocking face
      !! mask (`wd_open_u/v`, filled by the last BT substep's Pass 2)
      !! composes multiplicatively with the static one — layer momentum
      !! at a dynamically blocked (drying-front) face resets to zero, so
      !! rewetting starts from rest.  Knob off / absent ⇒ the ORIGINAL
      !! loops run untouched (byte-identical).
      !!
      !! When `&vcoord_nml zfixed_closed_faces` is on, the PER-LAYER
      !! z-level mask `metrics%open_u/open_v` composes multiplicatively
      !! too — this is the "no normal velocity" half of the z-level wall
      !! (the "no mass or tracer flux" half is continuity's, and
      !! "free-slip" is the horizontal-viscosity kernels').  It runs
      !! AFTER `apply_bt_correction` at every call site, so a barotropic
      !! increment can never be left behind at a closed face.  Knob off
      !! ⇒ the ORIGINAL loops, byte-identical.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(multilayer_state_t), intent(inout) :: ms
      type(barotropic_workstate_t), intent(in), optional :: bt_work
      integer :: i, j, k, nz, nx_face, ny_uface, nx_vface, ny_face
      logical :: wd_on

      nz = ms%nz_ml
      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)
      wd_on = .false.
      if (present(bt_work)) wd_on = bt_work%wetdry_enable

      if (wd_on) then
         do concurrent(k=1:nz, j=1:ny_uface, i=1:nx_face)
            ms%u_face_x_layer(i, j, k) = metrics%wet_u(i, j)* &
                                         bt_work%wd_open_u(i, j)* &
                                         ms%u_face_x_layer(i, j, k)
         end do
         do concurrent(k=1:nz, j=1:ny_face, i=1:nx_vface)
            ms%v_face_y_layer(i, j, k) = metrics%wet_v(i, j)* &
                                         bt_work%wd_open_v(i, j)* &
                                         ms%v_face_y_layer(i, j, k)
         end do
      else
         do concurrent(k=1:nz, j=1:ny_uface, i=1:nx_face)
            ms%u_face_x_layer(i, j, k) = metrics%wet_u(i, j)*ms%u_face_x_layer(i, j, k)
         end do
         do concurrent(k=1:nz, j=1:ny_face, i=1:nx_vface)
            ms%v_face_y_layer(i, j, k) = metrics%wet_v(i, j)*ms%v_face_y_layer(i, j, k)
         end do
      end if

      ! z-level closed faces: a SEPARATE host-gated pass, not folded into
      ! the branches above, so neither the all-wet nor the wet/dry loop
      ! changes textually when the knob is off.
      if (metrics%use_closed_faces) then
         do concurrent(k=1:nz, j=1:ny_uface, i=1:nx_face)
            ms%u_face_x_layer(i, j, k) = metrics%open_u(i, j, k)* &
                                         ms%u_face_x_layer(i, j, k)
         end do
         do concurrent(k=1:nz, j=1:ny_face, i=1:nx_vface)
            ms%v_face_y_layer(i, j, k) = metrics%open_v(i, j, k)* &
                                         ms%v_face_y_layer(i, j, k)
         end do
      end if
   end subroutine mask_layer_velocities