continuity_apply_fluxes Subroutine

   pure subroutine continuity_apply_fluxes(ms, dt, h_min)
      !! **Test-only** (no production caller): unsplit apply, paired with
      !! `continuity_compute_fluxes` as the split path's reference oracle.
      !! Per-layer forward-Euler thickness update.
      !!
      !! Optional `h_min` (m): when > 0 and the active vcoord is
      !! VCOORD_LAGRANGIAN, applies max(h_new, h_min) on the h-update so
      !! grounding layers cannot go below the floor.  0.0 (default absent) ⇒
      !! original update verbatim (bit-identical).  The `mass_budget_continuity`
      !! accumulator records the RAW divergence tendency regardless — the floor
      !! injection shows up in the mass Error diagnostic (R2).
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      real(wp), intent(in), optional :: h_min
         !! Minimum-thickness floor (m). 0 or absent ⇒ off ⇒ bit-identical.
      integer :: i, j, k, nx, ny, nz
      real(wp) :: h_min_use

      h_min_use = 0.0_wp
      if (present(h_min)) h_min_use = h_min

      nx = size(ms%h_layer, 1)
      ny = size(ms%h_layer, 2)
      nz = size(ms%h_layer, 3)

      if (h_min_use > 0.0_wp) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            ms%h_layer(i, j, k) = max(ms%h_layer(i, j, k) - dt*ms%flux_h_layer(i, j, k), h_min_use)
            ms%mass_budget_continuity(i, j, k) = &
               ms%mass_budget_continuity(i, j, k) - dt*ms%flux_h_layer(i, j, k)
         end do
      else
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            ms%h_layer(i, j, k) = ms%h_layer(i, j, k) - dt*ms%flux_h_layer(i, j, k)
            ms%mass_budget_continuity(i, j, k) = &
               ms%mass_budget_continuity(i, j, k) - dt*ms%flux_h_layer(i, j, k)
         end do
      end if
   end subroutine continuity_apply_fluxes