cavity_far_field_impl Subroutine

   pure subroutine cavity_far_field_impl(nx, ny, nz, far_depth, cover, wet_mask, &
                                         h_layer, hTr_T, hTr_S, u_face_x, v_face_y, &
                                         active, t_far, s_far, u_far, v_far)
      !! Thickness-weighted mean of `(T, S, u, v)` over `far_depth`
      !! METRES below the ice base, with a PARTIAL last layer.
      !!
      !! Bottom-up stack (`k = nz` is the surface layer, i.e. the one
      !! against the ice base), so the walk runs `k = nz, 1, -1` and
      !! stops when the budget is spent.  `far_depth >= column
      !! thickness` ⇒ the whole column, which is the correct limit and
      !! not an error.
      !!
      !! Vanished layers (`h <= H_VANISHED`) are SKIPPED, not clamped —
      !! the D4 taxonomy's skip/merge marker.  They carry no mass, so
      !! including them would be dividing a zero tracer load by a
      !! near-zero thickness.
      !!
      !! Velocities are centred from the C-grid faces per layer BEFORE
      !! the vertical average (`u_c = 1/2 (u_{i} + u_{i+1})`), not after:
      !! the two commute only for a uniform column, and the order that
      !! matches "the speed the ice base feels" is the one that averages
      !! the centred, per-layer velocity.
      !!
      !! `active` is the composed solve mask handed to the melt kernel:
      !! covered AND wet AND the sample found mass.  A covered column
      !! whose entire sample is vanished gets `active = 0` and zeroed
      !! outputs — no melt, no status noise.
      integer, intent(in) :: nx
         !! First dimension (ghosts included).
      integer, intent(in) :: ny
         !! Second dimension.
      integer, intent(in) :: nz
         !! Layer count; `k = nz` is the surface / ice-base layer.
      real(wp), intent(in) :: far_depth
         !! Sampling thickness (m), > 0.
      real(wp), intent(in) :: cover(nx, ny)
         !! Ice-cover fraction (v1 binary).
      real(wp), intent(in) :: wet_mask(nx, ny)
         !! Static wet (1) / land (0) mask.
      real(wp), intent(in) :: h_layer(nx, ny, nz)
         !! Layer thickness (m).
      real(wp), intent(in) :: hTr_T(nx, ny, nz)
         !! `h*T` (degC m).
      real(wp), intent(in) :: hTr_S(nx, ny, nz)
         !! `h*S` ((g/kg) m).
      real(wp), intent(in) :: u_face_x(nx + 1, ny, nz)
         !! Zonal face velocity (m/s).
      real(wp), intent(in) :: v_face_y(nx, ny + 1, nz)
         !! Meridional face velocity (m/s).
      real(wp), intent(out) :: active(nx, ny)
         !! Composed solve mask, 0 or 1.
      real(wp), intent(out) :: t_far(nx, ny)
         !! Sampled temperature (degC); 0 where inactive.
      real(wp), intent(out) :: s_far(nx, ny)
         !! Sampled salinity (g/kg); 0 where inactive.
      real(wp), intent(out) :: u_far(nx, ny)
         !! Sampled centred x velocity (m/s); 0 where inactive.
      real(wp), intent(out) :: v_far(nx, ny)
         !! Sampled centred y velocity (m/s); 0 where inactive.
      integer :: i, j, k
      real(wp) :: acc_h, acc_t, acc_s, acc_u, acc_v, remain, h, w, inv_h

      do concurrent(j=1:ny, i=1:nx) &
         local(k, acc_h, acc_t, acc_s, acc_u, acc_v, remain, h, w, inv_h)
         acc_h = 0.0_wp
         acc_t = 0.0_wp
         acc_s = 0.0_wp
         acc_u = 0.0_wp
         acc_v = 0.0_wp
         remain = far_depth
         if (cover(i, j) > 0.5_wp .and. wet_mask(i, j) > 0.5_wp) then
            do k = nz, 1, -1
               if (remain <= 0.0_wp) exit
               h = h_layer(i, j, k)
               ! vanished-ok: the far-field sampler SKIPS a vanished layer entirely (it carries
               ! no water to melt against) rather than weighting a substituted
               ! concentration into the average.
               if (h <= H_VANISHED) cycle
               w = min(h, remain)
               inv_h = 1.0_wp/h
               acc_h = acc_h + w
               acc_t = acc_t + w*hTr_T(i, j, k)*inv_h
               acc_s = acc_s + w*hTr_S(i, j, k)*inv_h
               acc_u = acc_u + w*0.5_wp*(u_face_x(i, j, k) + u_face_x(i + 1, j, k))
               acc_v = acc_v + w*0.5_wp*(v_face_y(i, j, k) + v_face_y(i, j + 1, k))
               remain = remain - w
            end do
         end if
         if (acc_h > 0.0_wp) then
            active(i, j) = 1.0_wp
            t_far(i, j) = acc_t/acc_h
            s_far(i, j) = acc_s/acc_h
            u_far(i, j) = acc_u/acc_h
            v_far(i, j) = acc_v/acc_h
         else
            active(i, j) = 0.0_wp
            t_far(i, j) = 0.0_wp
            s_far(i, j) = 0.0_wp
            u_far(i, j) = 0.0_wp
            v_far(i, j) = 0.0_wp
         end if
      end do
   end subroutine cavity_far_field_impl