apply_zonal_baroclinic Subroutine

   pure subroutine apply_zonal_baroclinic(u_layer, h_layer, ubt_end, &
                                          nx_total, ny_total, nz, &
                                          i_w, i_e, j0, j1, &
                                          bc_w, bc_e, clamp_u_w, clamp_u_e)
      !! Set west and east wall-face per-layer u to
      !! Flather mean + zero-gradient baroclinic anomaly (radiating edges)
      !! or uniform clamped_u (CLAMPED edge).
      !! Explicit-shape dummies; one do concurrent per edge (j outer, k inner).
      integer, intent(in) :: nx_total, ny_total, nz
      real(wp), intent(inout) :: u_layer(nx_total + 1, ny_total, nz)
         !! u_face_x_layer, shape (nx_total+1, ny_total, nz).
      real(wp), intent(in)    :: h_layer(nx_total, ny_total, nz)
      real(wp), intent(in)    :: ubt_end(nx_total + 1, ny_total)
         !! bt_ubt_end, shape (nx_total+1, ny_total).
      integer, intent(in) :: i_w, i_e   !! west / east physical wall-face indices
      integer, intent(in) :: j0, j1     !! first / last physical j-cell
      integer, intent(in) :: bc_w, bc_e
      real(wp), intent(in) :: clamp_u_w, clamp_u_e

      integer  :: j, k
      real(wp) :: h_tot, u_bar, u_int
      integer  :: i_int_w, i_int_e  ! first interior face index (one step inside)

      ! Interior face = one face step inward from the wall face.
      ! For west wall at i_w: interior face is i_w+1 (reads h_layer cell i_w).
      ! For east wall at i_e: interior face is i_e-1 (reads h_layer cell i_e-1).
      i_int_w = i_w + 1
      i_int_e = i_e - 1

      ! ---- West wall ----
      if (bc_w == OBC_CLAMPED) then
         do concurrent(k=1:nz, j=j0:j1)
            u_layer(i_w, j, k) = clamp_u_w
         end do
      else if (is_radiating(bc_w)) then
         do concurrent(j=j0:j1) local(k, h_tot, u_bar, u_int)
            ! Sequential k-loop for depth mean (nz is small, no reduction needed).
            h_tot = 0.0_wp
            u_bar = 0.0_wp
            do k = 1, nz
               ! Face i_int_w is between cell i_w and cell i_w+1.
               ! Use the upstream cell (i_w) thickness for the wall side.
               h_tot = h_tot + h_layer(i_w, j, k)
               u_bar = u_bar + u_layer(i_int_w, j, k)*h_layer(i_w, j, k)
            end do
            if (h_tot > 0.0_wp) then
               u_bar = u_bar/h_tot
            else
               u_bar = 0.0_wp
            end if
            do k = 1, nz
               u_int = u_layer(i_int_w, j, k)
               u_layer(i_w, j, k) = ubt_end(i_w, j) + (u_int - u_bar)
            end do
         end do
      end if

      ! ---- East wall ----
      if (bc_e == OBC_CLAMPED) then
         do concurrent(k=1:nz, j=j0:j1)
            u_layer(i_e, j, k) = clamp_u_e
         end do
      else if (is_radiating(bc_e)) then
         do concurrent(j=j0:j1) local(k, h_tot, u_bar, u_int)
            h_tot = 0.0_wp
            u_bar = 0.0_wp
            do k = 1, nz
               ! Use cell (i_e-1) thickness — the last interior cell.
               h_tot = h_tot + h_layer(i_e - 1, j, k)
               u_bar = u_bar + u_layer(i_int_e, j, k)*h_layer(i_e - 1, j, k)
            end do
            if (h_tot > 0.0_wp) then
               u_bar = u_bar/h_tot
            else
               u_bar = 0.0_wp
            end if
            do k = 1, nz
               u_int = u_layer(i_int_e, j, k)
               u_layer(i_e, j, k) = ubt_end(i_e, j) + (u_int - u_bar)
            end do
         end do
      end if
   end subroutine apply_zonal_baroclinic