ice_cas_to_ist_impl Subroutine

   pure subroutine ice_cas_to_ist_impl(wet_mask, areaT, mca_ice, mca_snow, part_size, m_ice, &
                                       m_snow, mh_lim, roll_factor, nghost, ncat, nx, ny)
      !! SIS2 `cell_ave_state_to_ice_state` (`:540`).  Per cell, per
      !! category: pre-floor category 1, optional thin-ice rolling
      !! (`roll_factor > 0`), general floor, then re-derive
      !! `part_size(c) = mca_ice(c)/m_ice(c)` and
      !! `m_snow(c) = m_ice(c)*(mca_snow(c)/mca_ice(c))` (per-ICE-area).
      !! `part_size(0) = 1 - Sum_c part_size(c)` — MAY be negative here;
      !! `ice_compress_impl` (Phase 4) is what restores >= 0.
      !!
      !! `roll_factor` rolling test — VERIFIED against the literal SIS2
      !! source (`SIS_transport.F90:572-581`, `L_to_H = US%L_to_Z*Rho_ice`;
      !! Roundabout carries no L_to_Z length rescaling, so `L_to_H =
      !! ICE_RHO_ICE`): a VOLUME-vs-VOLUME comparison, NOT a per-area one
      !! — `areaT` does NOT cancel (only the RHS carries it explicitly in
      !! SIS2's actual, non-simplified code):
      !!   `h_eff = m_ice(c)/ICE_RHO_ICE` (pre-roll per-ICE-area thickness,
      !!   a length);
      !!   `if (roll_factor*h_eff**3 > (mca_ice(c)/ICE_RHO_ICE)*areaT(i,j))`
      !!   — LHS a volume (h_eff cubed), RHS `(mca_ice/ICE_RHO_ICE)` is a
      !!   per-CELL-area ice-volume-per-area (a length) times `areaT` =
      !!   the cell's total ice volume;
      !!   then `m_ice(c) = max(mh_lim(1), ICE_RHO_ICE*sqrt((mca_ice(c)*
      !!   areaT(i,j)/ICE_RHO_ICE) / (roll_factor*h_eff)))` — SIS2's sqrt
      !!   argument is `(m_ice*areaT)/(roll_factor*mH_ice)` in ITS
      !!   pre-division units; converting consistently (dividing the
      !!   `mca_ice*areaT` term by `ICE_RHO_ICE` once, matching the
      !!   `(mca_ice/ICE_RHO_ICE)*areaT` volume on the LHS test) gives the
      !!   expression coded below.
      integer, intent(in) :: nghost, ncat, nx, ny
      real(wp), intent(in) :: wet_mask(nx, ny)
      real(wp), intent(in) :: areaT(nx, ny)
      real(wp), intent(in) :: mca_ice(nx, ny, ncat)
      real(wp), intent(in) :: mca_snow(nx, ny, ncat)
      real(wp), intent(inout) :: part_size(nx, ny, 0:ncat)
      real(wp), intent(inout) :: m_ice(nx, ny, ncat)
      real(wp), intent(inout) :: m_snow(nx, ny, ncat)
      real(wp), intent(in) :: mh_lim(ncat + 1)
      real(wp), intent(in) :: roll_factor
      integer :: i, j, c, i_lo, i_hi, j_lo, j_hi
      real(wp) :: h_eff, ice_vol, part_sum

      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, c=1:ncat) local(h_eff, ice_vol)
         if (wet_mask(i, j) > 0.5_wp) then
            if (c == 1 .and. mca_ice(i, j, 1) > 0.0_wp .and. m_ice(i, j, 1) < mh_lim(1)) then
               m_ice(i, j, 1) = mh_lim(1)
            end if
            if (mca_ice(i, j, c) > 0.0_wp) then
               if (roll_factor > 0.0_wp) then
                  h_eff = m_ice(i, j, c)/ICE_RHO_ICE
                  ice_vol = (mca_ice(i, j, c)/ICE_RHO_ICE)*areaT(i, j)
                  if (roll_factor*h_eff**3 > ice_vol) then
                     m_ice(i, j, c) = max(mh_lim(1), &
                                          ICE_RHO_ICE*sqrt(ice_vol/(roll_factor*h_eff)))
                  end if
               end if
               if (m_ice(i, j, c) < mh_lim(1)) m_ice(i, j, c) = mh_lim(1)
               part_size(i, j, c) = mca_ice(i, j, c)/m_ice(i, j, c)
               m_snow(i, j, c) = m_ice(i, j, c)*(mca_snow(i, j, c)/mca_ice(i, j, c))
            else
               part_size(i, j, c) = 0.0_wp
               m_ice(i, j, c) = 0.0_wp
               m_snow(i, j, c) = 0.0_wp
               ! enth/sal held (inert — massless convention)
            end if
         end if
      end do
      do concurrent(j=j_lo:j_hi, i=i_lo:i_hi) local(part_sum, c)
         if (wet_mask(i, j) > 0.5_wp) then
            part_sum = 0.0_wp
            do c = 1, ncat
               part_sum = part_sum + part_size(i, j, c)
            end do
            part_size(i, j, 0) = 1.0_wp - part_sum
         end if
      end do
   end subroutine ice_cas_to_ist_impl