ice_thermo_driver_reduce_impl Subroutine

   pure subroutine ice_thermo_driver_reduce_impl(wet_mask, h2o_ocn_to_ice, h2o_ice_to_ocn, &
                                                 heat_to_ocn, sw_thru, ssurf_seam, fb, &
                                                 heat_flux_diag, sw_thru_diag, m_melt_diag, &
                                                 salt_flux_diag, &
                                                 dt_therm, nghost, ncat, nx, ny)
      !! Device kernel over PHYSICAL cells (ghosts excluded). Reduces the
      !! per-category column outputs into the per-cell coupling diags.
      !! `heat_flux_diag`/`sw_thru_diag`/`m_melt_diag` are zeroed
      !! unconditionally at loop top (this module owns them outright);
      !! `salt_flux_diag` is NOT zeroed — `ice_frazil_uptake` (which runs
      !! first this window, see the module docstring) already zeroed +
      !! wrote it, and this kernel ADDS its net-melt contribution on top,
      !! gated on wet to avoid land noise (an unwetted cell adds exactly 0).
      !!
      !! PR 31: `sw_thru_diag(i,j) = Σ_cat sw_thru(i,j,cat)` (W/m^2, +down)
      !! — the lumped ncat==1 reduction (no part-weight, matching
      !! `heat_flux_diag`). It is written to its OWN field, NEVER folded
      !! into `heat_flux_diag`: `heat_flux_diag` stays the non-shortwave
      !! heat share, and `ice_ocean_sw_flux` delivers `sw_thru_diag` to
      !! `Q_heat` exactly once (see that routine for the no-double-count
      !! argument across both component modes).
      !!
      !! ICE-FREE CELLS ARE A CLEAN NO-OP: on an ice-free wet cell the column
      !! zeroes all its per-cat outputs (h2o_*/heat_to_ocn = 0), and
      !! `ice_compute_basal_flux` gates `fb = 0` there (no ice base ⇒ no basal
      !! flux — the same ice-presence threshold the column uses). So
      !! `heat_flux_diag = 0/dt - 0 = 0`, `m_melt_diag = 0`, and the net-melt
      !! salt add is `0*(...) = 0`. The whole open-ocean surface therefore sees
      !! Q_heat = Q_salt = 0 from the ice path, as physics requires.
      !!
      !! Decl-order: all integer dims declared before the explicit-shape
      !! arrays that use them. Inner `if`/serial `do cat` — never a
      !! masked `do concurrent` header.
      integer, intent(in) :: nghost, ncat, nx, ny
      real(wp), intent(in) :: wet_mask(nx, ny)
      real(wp), intent(in) :: h2o_ocn_to_ice(nx, ny, ncat)
      real(wp), intent(in) :: h2o_ice_to_ocn(nx, ny, ncat)
      real(wp), intent(in) :: heat_to_ocn(nx, ny, ncat)
      real(wp), intent(in) :: sw_thru(nx, ny, ncat)
      real(wp), intent(in) :: ssurf_seam(nx, ny)
      real(wp), intent(in) :: fb(nx, ny)
      real(wp), intent(inout) :: heat_flux_diag(nx, ny)
      real(wp), intent(inout) :: sw_thru_diag(nx, ny)
      real(wp), intent(inout) :: m_melt_diag(nx, ny)
      real(wp), intent(inout) :: salt_flux_diag(nx, ny)
      real(wp), intent(in) :: dt_therm

      integer :: i, j, cat, i_lo, i_hi, j_lo, j_hi
      real(wp) :: sum_ocn, sum_ice, sum_heat, sum_sw, m_net

      i_lo = nghost + 1
      i_hi = nx - nghost
      j_lo = nghost + 1
      j_hi = ny - nghost

      do concurrent(j=j_lo:j_hi, i=i_lo:i_hi) &
         local(cat, sum_ocn, sum_ice, sum_heat, sum_sw, m_net)
         heat_flux_diag(i, j) = 0.0_wp
         sw_thru_diag(i, j) = 0.0_wp
         m_melt_diag(i, j) = 0.0_wp
         if (wet_mask(i, j) > 0.5_wp) then
            sum_ocn = 0.0_wp
            sum_ice = 0.0_wp
            sum_heat = 0.0_wp
            sum_sw = 0.0_wp
            do cat = 1, ncat
               sum_ocn = sum_ocn + h2o_ocn_to_ice(i, j, cat)
               sum_ice = sum_ice + h2o_ice_to_ocn(i, j, cat)
               sum_heat = sum_heat + heat_to_ocn(i, j, cat)
               sum_sw = sum_sw + sw_thru(i, j, cat)
            end do
            m_net = sum_ocn - sum_ice
            m_melt_diag(i, j) = sum_ice - sum_ocn
            heat_flux_diag(i, j) = sum_heat/dt_therm - fb(i, j)
            sw_thru_diag(i, j) = sum_sw
            salt_flux_diag(i, j) = salt_flux_diag(i, j) &
                                   + m_net*(ssurf_seam(i, j) - ICE_BULK_SALINITY)/dt_therm
         end if
      end do
   end subroutine ice_thermo_driver_reduce_impl