ice_column_step Subroutine

   pure subroutine ice_column_step(nk, m_snow, m_ice_tot, enth_snow_pt, &
                                   enth_ice_bu, sal_ice_bu, &
                                   sf_0, dsf_dt, sw_dn, tfw, fb, sst, s_surf, dtt, &
                                   do_snow_ice, &
                                   snow, tsurf, h2o_ocn_to_ice, h2o_ice_to_ocn, &
                                   heat_to_ocn, sw_thru, snow_to_ice)
      !! Per-(cell,category) orchestrator: gather (bottom-up -> top-down
      !! flip, TRAP #2), optics, conduction (`ice_temp_sis2`), resize
      !! (snow add, bottom-freeze, top/bottom melt peel, rebalance),
      !! scatter (flip back).
      !!
      !! PR 26 TRAP: `snow` is added inside step 4 (resize), i.e. AFTER
      !! the optics (step 2) and conduction (step 3) already ran on the
      !! PRE-snowfall `m_snow`. That is SIS2's fast/slow split
      !! (`ice_resize_SIS2` runs after the conduction solve,
      !! SIS2_ice_thm.F90:1122) — new snow IS meltable in this same
      !! window (`ice_top_melt_peel` starts at k=0), but its albedo and
      !! conduction effect are felt only on the NEXT window. Do not
      !! "helpfully" move the add earlier.
      !$acc routine seq
      integer, intent(in) :: nk
         !! Number of ice layers (declared first — decl-order).
      real(wp), intent(inout) :: m_snow
         !! Snow mass per unit area (kg/m²).
      real(wp), intent(inout) :: m_ice_tot
         !! Total ice mass per unit area (kg/m²).
      real(wp), intent(inout) :: enth_snow_pt
         !! Snow specific enthalpy (J/kg).
      real(wp), intent(inout) :: enth_ice_bu(nk)
         !! BOTTOM-UP ice specific enthalpies (J/kg) — state order,
         !! `enth_ice_bu(1)` = ice bottom.
      real(wp), intent(inout) :: sal_ice_bu(nk)
         !! BOTTOM-UP ice bulk salinities (PSU) — state order.
      real(wp), intent(in) :: sf_0
         !! Linearized SEB intercept (W/m²), upward-positive.
      real(wp), intent(in) :: dsf_dt
         !! Linearized SEB slope (W/m²/K), upward-positive.
      real(wp), intent(in) :: sw_dn
         !! Downwelling shortwave at the surface (W/m²).
      real(wp), intent(in) :: tfw
         !! Seawater freezing temperature at the ice base (degC).
      real(wp), intent(in) :: fb
         !! Ocean -> ice-base heat flux (W/m²).
      real(wp), intent(in) :: sst
         !! Sea-surface temperature (degC) — feeds the TRAP-#1 liquid
         !! ocean enthalpy.
      real(wp), intent(in) :: s_surf
         !! Sea-surface salinity (PSU) — feeds the TRAP-#1 liquid
         !! ocean enthalpy (unused by the linear formula but kept for
         !! call-site parity, `ice_enthalpy_liquid`).
      real(wp), intent(in) :: dtt
         !! Timestep (s).
      logical, intent(in) :: do_snow_ice
         !! Archimedes freeboard flood gate (`&ocean_ice_nml snow_ice`,
         !! PR 27) — `.false.` (the default) is a bit-identical no-op:
         !! `snow_to_ice` stays 0 and `ice_snow_ice_flood` is never
         !! called.
      real(wp), intent(in) :: snow
         !! New snow mass this window (kg/m²), `fprec*dtt` — PR 26 source
         !! term, `ice_snow_accumulate`'s `snow` argument. 0 (the
         !! `&ocean_ice_nml snowfall=0` default) is a bit-identical no-op.
      real(wp), intent(out) :: tsurf
         !! Surface skin temperature (degC).
      real(wp), intent(out) :: h2o_ocn_to_ice
         !! Mass flux frozen from the ocean onto the ice base (kg/m²).
      real(wp), intent(out) :: h2o_ice_to_ocn
         !! Meltwater mass flux to the ocean (kg/m²), top + bottom peel.
      real(wp), intent(out) :: heat_to_ocn
         !! Leftover melt heat dumped to the ocean (J/m²), top + bottom.
      real(wp), intent(out) :: sw_thru
         !! Shortwave transmitted through the ice to the ocean (W/m²).
      real(wp), intent(out) :: snow_to_ice
         !! Mass converted from snow to the top ice layer this call
         !! (kg/m²), >= 0 — SIS2 `SN2IC` (PR 27). 0 when `do_snow_ice`
         !! is `.false.` or the column is not flooded.

      real(wp) :: enth_loc(0:ICE_NK_MAX), sal_loc(ICE_NK_MAX)
      real(wp) :: albedo, abs_sfc, abs_snow, abs_ocn, abs_int, pen
      real(wp) :: abs_ice_lay(ICE_NK_MAX), sol(0:ICE_NK_MAX)
      real(wp) :: ts_opt, sw_tot, sf_0_eff
      real(wp) :: tmelt, bmelt
      real(wp) :: col_enth_in, col_enth_out, sum_sol, tflux_sfc, tflux_bot
      real(wp) :: m_lay(0:ICE_NK_MAX), enthalpy(0:ICE_NK_MAX + 1), salin(0:ICE_NK_MAX)
      real(wp) :: enth_ocean, salin_freeze, mtot_ice
      integer :: k

      ! ---- 1. Gather + flip (TRAP #2) ----
      do k = 1, nk
         enth_loc(k) = enth_ice_bu(nk + 1 - k)
         sal_loc(k) = sal_ice_bu(nk + 1 - k)
      end do
      enth_loc(0) = enth_snow_pt
      if (m_snow == 0.0_wp) then
         ! Massless snow slot: re-seed from the top ice layer's
         ! temperature at 0 salinity (SIS_slow_thermo.F90:977).
         enth_loc(0) = ice_enth_from_ts(ice_temp_from_en_s(enth_loc(1), sal_loc(1)), 0.0_wp)
      end if

      ! ---- 2. Optics -> sol ----
      if (m_snow > 0.0_wp) then
         ts_opt = ice_temp_from_en_s(enth_loc(0), 0.0_wp)
      else
         ts_opt = ice_temp_from_en_s(enth_loc(1), sal_loc(1))
         ! TODO(PR-3b): carry a true prognostic Tskin; this reuses the
         ! top-ice-layer temperature as a skin-temp stand-in.
      end if
      call ice_optics_csim4(nk, m_snow/ICE_RHO_SNOW, m_ice_tot/ICE_RHO_ICE, &
                            ts_opt, sal_loc(1), albedo, abs_sfc, abs_snow, &
                            abs_ice_lay(1:nk), abs_ocn, abs_int, pen)

      sw_tot = (1.0_wp - albedo)*sw_dn
      sf_0_eff = sf_0 - abs_sfc*sw_tot
      sol(0) = abs_snow*sw_tot
      do k = 1, nk
         sol(k) = abs_ice_lay(k)*sw_tot
      end do
      sw_thru = abs_ocn*sw_tot

      ! ---- 3. Conduction ----
      tmelt = 0.0_wp
      bmelt = 0.0_wp
      enthalpy(0:nk) = enth_loc(0:nk)
      call ice_temp_sis2(nk, m_snow, m_ice_tot, sal_loc(1:nk), enthalpy(0:nk), &
                         sf_0_eff, dsf_dt, sol(0:nk), tfw, fb, dtt, &
                         tsurf, tmelt, bmelt, &
                         col_enth_in, col_enth_out, sum_sol, tflux_sfc, tflux_bot)
      enth_loc(0:nk) = enthalpy(0:nk)

      ! ---- 4. Resize ----
      m_lay(0) = m_snow
      do k = 1, nk
         m_lay(k) = m_ice_tot/real(nk, wp)
      end do
      enth_ocean = ice_enthalpy_liquid(sst, s_surf)  ! TRAP #1
      salin_freeze = ICE_BULK_SALINITY

      salin(0) = 0.0_wp
      salin(1:nk) = sal_loc(1:nk)
      enthalpy(0:nk) = enth_loc(0:nk)

      ! PR 26: snow source term — SIS2's ice_resize_SIS2 runs this FIRST,
      ! before the melt peels (SIS2_ice_thm.F90:1122), so new snow is
      ! meltable in this same window. `enthalpy(0)` is unchanged by the
      ! add (see ice_snow_accumulate's docstring).
      call ice_snow_accumulate(nk, m_lay(0:nk), snow)

      ! Negative-top-melt fold (do_pond=false path, SIS2:1159-1163) —
      ! unreachable in the PR-3a gates but required for SIS2 parity.
      if (tmelt < 0.0_wp) then
         bmelt = bmelt + tmelt
         tmelt = 0.0_wp
      end if

      call ice_bottom_freeze(nk, m_lay(0:nk), enthalpy(0:nk), salin(0:nk), &
                             bmelt, enth_ocean, salin_freeze, h2o_ocn_to_ice)

      heat_to_ocn = 0.0_wp
      h2o_ice_to_ocn = 0.0_wp
      call ice_top_melt_peel(nk, m_lay(0:nk), enthalpy(0:nk), salin(0:nk), tmelt, &
                             heat_to_ocn, h2o_ice_to_ocn)
      call ice_bottom_melt_peel(nk, m_lay(0:nk), enthalpy(0:nk), salin(0:nk), bmelt, &
                                heat_to_ocn, h2o_ice_to_ocn)

      ! PR 27: Archimedes freeboard flood — AFTER the melt peels (so the
      ! freshly-converted mass is not re-melted this step and m_i
      ! reflects the post-melt column) and BEFORE ice_rebalance_layers
      ! (so the new mass in layer 1 IS redistributed across the nk
      ! layers) — SIS2's exact order (TRAP #7). `do_snow_ice=.false.`
      ! (the default) is a bit-identical no-op by inspection.
      snow_to_ice = 0.0_wp
      if (do_snow_ice) then
         call ice_snow_ice_flood(nk, m_lay(0:nk), enthalpy(0:nk), salin(0:nk), &
                                 ICE_RHO_ICE/ICE_RHO_OCEAN, snow_to_ice)
      end if

      call ice_rebalance_layers(nk, m_lay(0:nk), enthalpy(0:nk), salin(0:nk), mtot_ice)

      ! ---- 5. Scatter + flip back ----
      m_ice_tot = mtot_ice
      m_snow = m_lay(0)
      enth_snow_pt = enthalpy(0)
      do k = 1, nk
         enth_ice_bu(nk + 1 - k) = enthalpy(k)
         sal_ice_bu(nk + 1 - k) = salin(k)
      end do
   end subroutine ice_column_step