ice_temp_sis2 Subroutine

   pure subroutine ice_temp_sis2(nk, m_snow, m_ice_tot, sice, enthalpy, &
                                 sf_0, dsf_dt, sol, tfw, fb, dtt, &
                                 tsurf, tmelt, bmelt, &
                                 col_enth_in, col_enth_out, sum_sol, &
                                 tflux_sfc, tflux_bot)
      !! SEB + vertical-conduction column solve — SIS2 `ice_temp_SIS2`
      !! (SIS2_ice_thm.F90:169-540). Port of prototype
      !! `sis2_column.py:151-412`. TOP-DOWN column (index 0 = snow,
      !! 1..nk = ice top->bottom) — see module docstring TRAP #2.
      !!
      !! The five diag outputs (`col_enth_in`, `col_enth_out`,
      !! `sum_sol`, `tflux_sfc`, `tflux_bot`) are ALWAYS computed (a
      !! handful of flops on a small column) and feed the
      !! `column_energy_closure` test's identity: `col_enth_out -
      !! col_enth_in == sum_sol + tflux_sfc + tflux_bot` (`tflux_bot`
      !! ADDED — already the signed contribution, prototype :198-227).
      !! `col_enth_out` is measured AFTER the conservative update but
      !! BEFORE the liq-lim clamp (prototype `col_enth2b`) — the clamp
      !! moves energy into tmelt/bmelt, outside this identity.
      !$acc routine seq
      integer, intent(in) :: nk
         !! Number of ice layers (declared first — decl-order).
      real(wp), intent(in) :: m_snow
         !! Snow mass per unit area (kg/m²).
      real(wp), intent(in) :: m_ice_tot
         !! Total ice mass per unit area (kg/m²).
      real(wp), intent(in) :: sice(nk)
         !! TOP-DOWN ice bulk salinities (PSU).
      real(wp), intent(inout) :: enthalpy(0:nk)
         !! TOP-DOWN specific enthalpies (J/kg): 0 = snow, 1..nk = ice.
      real(wp), intent(in) :: sf_0
         !! Linearized SEB intercept (W/m²), upward-positive: `SF(T) =
         !! sf_0 + dsf_dt*T`.
      real(wp), intent(in) :: dsf_dt
         !! Linearized SEB slope (W/m²/K), upward-positive.
      real(wp), intent(in) :: sol(0:nk)
         !! Absorbed solar per layer (W/m²), TOP-DOWN.
      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²); post-hoc residual only
         !! — TRAP #3.
      real(wp), intent(in) :: dtt
         !! Timestep (s).
      real(wp), intent(out) :: tsurf
         !! Surface skin temperature (degC).
      real(wp), intent(inout) :: tmelt
         !! Accumulated top melting energy (J/m²); caller zeroes per step.
      real(wp), intent(inout) :: bmelt
         !! Accumulated bottom melting/freezing energy (J/m²); caller
         !! zeroes per step.
      real(wp), intent(out) :: col_enth_in
         !! Column enthalpy Σ m_lay*enth BEFORE anything (diag).
      real(wp), intent(out) :: col_enth_out
         !! Column enthalpy Σ m_lay*enth AFTER the conservative update,
         !! BEFORE the liq-lim clamp (diag).
      real(wp), intent(out) :: sum_sol
         !! Σ sol*dtt over the column (diag, J/m²).
      real(wp), intent(out) :: tflux_sfc
         !! Time-integrated surface heat flux into the column (diag,
         !! J/m²).
      real(wp), intent(out) :: tflux_bot
         !! Time-integrated basal heat flux into the column (diag,
         !! J/m²).

      real(wp) :: temp_ic(0:ICE_NK_MAX), tfi(ICE_NK_MAX)
      real(wp) :: ml_ice, ml_snow, hl_ice_eff, hsnow_eff, tsf
      real(wp) :: kk, k10, k0a, k0skin, k0a_x_ta
      real(wp) :: m_lay(0:ICE_NK_MAX)
      real(wp) :: bb(0:ICE_NK_MAX), cc(0:ICE_NK_MAX + 1), cc_bb(0:ICE_NK_MAX)
      real(wp) :: temp_est(0:ICE_NK_MAX)
      real(wp) :: heat_flux_int(-1:ICE_NK_MAX)
      real(wp) :: b_denom_1, i_bb, comp_rat, tsurf_est, m_pond
      real(wp) :: heat_flux_err_rat
      real(wp) :: e_extra, e_extra_sum, ftop_new, fbot_new, snow_temp_max, snow_temp_new
      real(wp) :: enth_liq_lim, i_liq_lim
      integer :: k

      ! ---- T<->E inversion of the incoming state (top-down) ----
      temp_ic(0) = ice_temp_from_en_s(enthalpy(0), 0.0_wp)
      do k = 1, nk
         temp_ic(k) = ice_temp_from_en_s(enthalpy(k), sice(k))
      end do

      ml_ice = m_ice_tot/real(nk, wp)
      ml_snow = m_snow
      do k = 1, nk
         tfi(k) = ice_t_freeze(sice(k))
      end do

      hl_ice_eff = max(ml_ice/ICE_RHO_ICE, ICE_H_LO_LIM)
      hsnow_eff = ml_snow/ICE_RHO_SNOW + max(1.0e-35_wp, 1.0e-20_wp*ICE_H_LO_LIM)

      tsf = tfi(1)
      if (ml_snow > 0.0_wp) tsf = 0.0_wp

      kk = ICE_K_ICE/hl_ice_eff
      k10 = 2.0_wp*(ICE_K_SNOW*ICE_K_ICE)/(hl_ice_eff*ICE_K_SNOW + hsnow_eff*ICE_K_ICE)
      k0a = (ICE_K_SNOW*dsf_dt)/(0.5_wp*dsf_dt*hsnow_eff + ICE_K_SNOW)
      k0skin = 2.0_wp*ICE_K_SNOW/hsnow_eff
      k0a_x_ta = (ICE_K_SNOW*sf_0)/(0.5_wp*dsf_dt*hsnow_eff + ICE_K_SNOW)

      m_lay(0) = ml_snow
      do k = 1, nk
         m_lay(k) = ml_ice
      end do

      col_enth_in = 0.0_wp
      do k = 0, nk
         col_enth_in = col_enth_in + m_lay(k)*enthalpy(k)
      end do

      ! ---- Effective layer heat capacities bb(k) — TRAP #4 ----
      bb(0) = ml_snow*ICE_CP_ICE
      do k = 1, nk
         if (tfi(k) >= 0.0_wp) then
            bb(k) = ml_ice*ICE_CP_ICE
         else if (temp_ic(k) < tfi(k)) then
            bb(k) = ml_ice*(ICE_CP_ICE - (tfi(k)/temp_ic(k)**2)* &
                            (ICE_LAT_FUS - (ICE_CP_BRINE - ICE_CP_ICE)*temp_ic(k)))
         else
            bb(k) = ml_ice*(ICE_CP_BRINE - ICE_LAT_FUS/tfi(k))
         end if
      end do

      ! ---- Coupling coefficients cc — TRAP #3 (bottom couples to tfw) ----
      cc(0) = k0a*dtt
      cc(1) = k10*dtt
      do k = 2, nk
         cc(k) = kk*dtt
      end do
      cc(nk + 1) = 2.0_wp*kk*dtt

      ! ---- UP sweep ----
      b_denom_1 = bb(nk) + cc(nk + 1)
      i_bb = 1.0_wp/(b_denom_1 + cc(nk))
      temp_est(nk) = ((sol(nk)*dtt + bb(nk)*temp_ic(nk)) + cc(nk + 1)*tfw)*i_bb
      comp_rat = b_denom_1*i_bb
      cc_bb(nk) = cc(nk)*i_bb

      do k = nk - 1, 1, -1
         b_denom_1 = bb(k) + comp_rat*cc(k + 1)
         i_bb = 1.0_wp/(b_denom_1 + cc(k))
         temp_est(k) = ((sol(k)*dtt + bb(k)*temp_ic(k)) + cc(k + 1)*temp_est(k + 1))*i_bb
         comp_rat = b_denom_1*i_bb
         cc_bb(k) = cc(k)*i_bb
      end do

      b_denom_1 = bb(0) + comp_rat*cc(1)
      i_bb = 1.0_wp/(b_denom_1 + cc(0))
      temp_est(0) = (((sol(0)*dtt + bb(0)*temp_ic(0)) - k0a_x_ta*dtt) + cc(1)*temp_est(1))*i_bb

      tsurf_est = (k0skin*temp_est(0) - sf_0)/(dsf_dt + k0skin)

      m_pond = 0.0_wp
      if (tsurf_est > tsf .or. m_pond > 0.0_wp) then
         tsurf_est = tsf
         i_bb = 1.0_wp/(b_denom_1 + k0skin*dtt)
         temp_est(0) = min(tsf, &
                           (((sol(0)*dtt + bb(0)*temp_ic(0)) + k0skin*dtt*tsf) + &
                            cc(1)*temp_est(1))*i_bb)
      end if

      ! ---- DOWN sweep ----
      do k = 1, nk
         temp_est(k) = min(temp_est(k) + cc_bb(k)*temp_est(k - 1), tfi(k))
      end do

      ! ---- Quasi-conservative re-solve via laytemp_sis2 — TRAP #5 ----
      if (nk == 1) then
         temp_est(1) = laytemp_sis2(ml_ice, tfi(1), &
                                    sol(1) + (2.0_wp*kk*tfw + k10*temp_est(0)), &
                                    2.0_wp*kk + k10, temp_ic(1), dtt)
      else
         temp_est(nk) = laytemp_sis2(ml_ice, tfi(nk), &
                                     sol(nk) + kk*(2.0_wp*tfw + temp_est(nk - 1)), &
                                     3.0_wp*kk, temp_ic(nk), dtt)
         do k = nk - 1, 2, -1
            temp_est(k) = laytemp_sis2(ml_ice, tfi(k), &
                                       sol(k) + kk*(temp_est(k - 1) + temp_est(k + 1)), &
                                       2.0_wp*kk, temp_ic(k), dtt)
         end do
         temp_est(1) = laytemp_sis2(ml_ice, tfi(1), &
                                    sol(1) + (kk*temp_est(2) + k10*temp_est(0)), &
                                    kk + k10, temp_ic(1), dtt)
      end if

      temp_est(0) = laytemp_sis2(ml_snow, 0.0_wp, &
                                 sol(0) + (k10*temp_est(1) - k0a_x_ta), &
                                 k10 + k0a, temp_ic(0), dtt)
      tsurf = (k0skin*temp_est(0) - sf_0)/(dsf_dt + k0skin)

      ! ---- Conservative DOWN pass: actually update enthalpies ----
      heat_flux_err_rat = 0.7071_wp*dtt*ICE_TEMP_RANGE_EST/ &
                          (ICE_TEMP_RANGE_EST*ICE_CP_ICE + ICE_LAT_FUS)

      e_extra_sum = 0.0_wp
      sum_sol = 0.0_wp
      do k = 0, nk
         sum_sol = sum_sol + sol(k)
      end do
      sum_sol = sum_sol*dtt

      if (tsurf > tsf .or. m_pond > 0.0_wp) then
         tsurf = tsf
         if (ml_snow > 0.0_wp) then
            heat_flux_int(-1) = k0skin*tsf
            heat_flux_int(0) = -k10*temp_est(1)
            call update_lay_enth(ml_snow, 0.0_wp, enthalpy(0), heat_flux_int(-1), &
                                 sol(0), heat_flux_int(0), -k0skin, k10, dtt, &
                                 heat_flux_err_rat, e_extra, snow_temp_new, &
                                 .false., 0.0_wp)
            tmelt = tmelt + e_extra - dtt*((sf_0 + dsf_dt*tsf) + heat_flux_int(-1))
            e_extra_sum = e_extra_sum + e_extra
            tflux_sfc = dtt*heat_flux_int(-1)
         else
            enthalpy(0) = ice_enth_from_ts(tsf, 0.0_wp)
            heat_flux_int(0) = k10*(tsf - temp_est(1))
            heat_flux_int(-1) = heat_flux_int(0)
            tmelt = tmelt + dtt*((sol(0) - (sf_0 + dsf_dt*tsf)) - heat_flux_int(0))
            tflux_sfc = dtt*heat_flux_int(0)
         end if
      else
         heat_flux_int(-1) = -k0a_x_ta
         heat_flux_int(0) = -k10*temp_est(1)
         snow_temp_max = (tsf*(dsf_dt + k0skin) + sf_0)/k0skin
         call update_lay_enth(ml_snow, 0.0_wp, enthalpy(0), heat_flux_int(-1), &
                              sol(0), heat_flux_int(0), -k0a, k10, dtt, &
                              heat_flux_err_rat, e_extra, snow_temp_new, &
                              .true., snow_temp_max)
         tsurf = (k0skin*snow_temp_new - sf_0)/(dsf_dt + k0skin)
         e_extra_sum = e_extra_sum + e_extra
         tmelt = tmelt + e_extra
         tflux_sfc = dtt*heat_flux_int(-1)
      end if

      do k = 1, nk - 1
         heat_flux_int(k) = -kk*temp_est(k + 1)
         ftop_new = heat_flux_int(k - 1)
         fbot_new = heat_flux_int(k)
         call update_lay_enth(ml_ice, sice(k), enthalpy(k), ftop_new, &
                              sol(k), fbot_new, 0.0_wp, kk, dtt, &
                              heat_flux_err_rat, e_extra, snow_temp_new, &
                              .false., 0.0_wp)
         heat_flux_int(k - 1) = ftop_new
         heat_flux_int(k) = fbot_new
         e_extra_sum = e_extra_sum + e_extra
         if (k <= nk/2) then
            tmelt = tmelt + e_extra
         else
            bmelt = bmelt + e_extra
         end if
      end do

      heat_flux_int(nk) = -2.0_wp*kk*tfw
      ftop_new = heat_flux_int(nk - 1)
      fbot_new = heat_flux_int(nk)
      call update_lay_enth(ml_ice, sice(nk), enthalpy(nk), ftop_new, &
                           sol(nk), fbot_new, 0.0_wp, 2.0_wp*kk, dtt, &
                           heat_flux_err_rat, e_extra, snow_temp_new, &
                           .false., 0.0_wp)
      heat_flux_int(nk - 1) = ftop_new
      heat_flux_int(nk) = fbot_new
      e_extra_sum = e_extra_sum + e_extra
      bmelt = bmelt + e_extra
      ! ---- END conservative update of enthalpy ----

      col_enth_out = 0.0_wp
      do k = 0, nk
         col_enth_out = col_enth_out + m_lay(k)*enthalpy(k)
      end do

      tflux_bot = -heat_flux_int(nk)*dtt

      ! TRAP #3: fb enters ONLY here, as a post-hoc bmelt residual.
      bmelt = bmelt + (dtt*fb - tflux_bot)

      ! ---- Excess-heat clamp to liq_lim (SIS2:500-524) ----
      enth_liq_lim = ice_enth_from_ts(0.0_wp, 0.0_wp)
      if (enthalpy(0) > enth_liq_lim) then
         e_extra = (enthalpy(0) - enth_liq_lim)*ml_snow
         tmelt = tmelt + e_extra
         enthalpy(0) = enth_liq_lim
      end if

      i_liq_lim = 1.0_wp/ICE_LIQ_LIM
      do k = 1, nk
         enth_liq_lim = ice_enth_from_ts(tfi(k)*i_liq_lim, sice(k))
         if (enthalpy(k) > enth_liq_lim) then
            e_extra = (enthalpy(k) - enth_liq_lim)*ml_ice
            enthalpy(k) = enth_liq_lim
            if (k <= nk/2) then
               tmelt = tmelt + e_extra
            else
               bmelt = bmelt + e_extra
            end if
         end if
      end do
   end subroutine ice_temp_sis2