ice_thermo_driver_reduce_multicat_impl Subroutine

   pure subroutine ice_thermo_driver_reduce_multicat_impl(wet_mask, part_size, h2o_ocn_to_ice, &
                                                          h2o_ice_to_ocn, heat_to_ocn, sw_thru, &
                                                          ssurf_seam, fb, fb_part_sum, &
                                                          heat_flux_diag, sw_thru_diag, &
                                                          m_melt_diag, &
                                                          salt_flux_diag, dt_therm, nghost, &
                                                          ncat, nx, ny)
      !! ncat>1 SIS2 ITD-mode reduce (module docstring). Same zero/gate/
      !! ordering contract as `ice_thermo_driver_reduce_impl` — the only
      !! change is PART-WEIGHTING the per-category sums (the column's
      !! `h2o_*`/`heat_to_ocn`/`sw_thru` outputs are per unit ICE-COVERED
      !! area in this mode, so a per-cell total needs the `part_size`
      !! weight) and subtracting `fb_part_sum*fb` instead of the bare `fb`
      !! the ncat==1 path subtracts (`fb` is a per-cell flux; `fb_part_sum`
      !! is the fraction of the cell it was actually charged against —
      !! see the module docstring and `rdb_ice_state%fb_part_sum`).
      !!
      !! PR 31: `sw_thru_diag(i,j) = Σ_cat part_size(i,j,cat)*sw_thru(i,j,cat)`
      !! (W/m^2, +down) — the area-weighted twin of the ncat==1 reduce's
      !! lumped sum, matching `heat_to_ocn`'s own part-weighting here.
      !! Same OWN-field / never-folded-into-`heat_flux_diag` contract as
      !! the ncat==1 path.
      !!
      !! 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) :: part_size(nx, ny, 0:ncat)
      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(in) :: fb_part_sum(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 + part_size(i, j, cat)*h2o_ocn_to_ice(i, j, cat)
               sum_ice = sum_ice + part_size(i, j, cat)*h2o_ice_to_ocn(i, j, cat)
               sum_heat = sum_heat + part_size(i, j, cat)*heat_to_ocn(i, j, cat)
               sum_sw = sum_sw + part_size(i, j, cat)*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_part_sum(i, j)*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_multicat_impl