update_lay_enth Subroutine

   pure subroutine update_lay_enth(m_lay, sice, enth, ftop, ht_body, fbot, &
                                   dftop_dt, dfbot_dt, dtt, hf_err_rat, &
                                   extra_heat, new_temp, has_temp_max, temp_max)
      !! Conservative per-layer implicit enthalpy update — SIS2
      !! `update_lay_enth` (SIS2_ice_thm.F90:704-945), closed-form
      !! branches only. Port of prototype `sis2_column.py:58-135`.
      !! Four solution branches (massless layer; pin-to-max with
      !! banked `extra_enth`; fresh `sice==0` linear; salty quadratic),
      !! then the three-way explicit-vs-conservation-inverted flux
      !! bookkeeping (prototype :117-133, incl. the `denom > 0` guard).
      !!
      !! `temp_max` is optional in the SIS2 signature; here it is a
      !! `has_temp_max` logical + `temp_max` value pair (device-routine
      !! `optional` dummies are avoided — same-module call sites only).
      !$acc routine seq
      real(wp), intent(in) :: m_lay
         !! Layer mass (kg/m²).
      real(wp), intent(in) :: sice
         !! Layer bulk salinity (PSU); 0 for snow/fresh.
      real(wp), intent(inout) :: enth
         !! Layer specific enthalpy (J/kg); in = prior step, out = new.
      real(wp), intent(inout) :: ftop
         !! Heat flux at the layer's top interface (W/m²); in = prior
         !! estimate, out = updated (explicit or conservation-inverted).
      real(wp), intent(in) :: ht_body
         !! In-layer heating (solar absorption) (W/m²).
      real(wp), intent(inout) :: fbot
         !! Heat flux at the layer's bottom interface (W/m²); in/out as
         !! `ftop`.
      real(wp), intent(in) :: dftop_dt
         !! d(ftop)/d(new_temp) (W/m²/K).
      real(wp), intent(in) :: dfbot_dt
         !! d(fbot)/d(new_temp) (W/m²/K).
      real(wp), intent(in) :: dtt
         !! Timestep (s).
      real(wp), intent(in) :: hf_err_rat
         !! Precomputed `heat_flux_err_rat` (degC*s/J) deciding explicit
         !! vs conservation-inverted flux bookkeeping.
      real(wp), intent(out) :: extra_heat
         !! Banked excess heat when pinned to `temp_max` (J/m²).
      real(wp), intent(out) :: new_temp
         !! Resulting layer temperature (degC).
      logical, intent(in) :: has_temp_max
         !! True when an explicit `temp_max` clamp applies (snow-branch
         !! call site); false uses the freezing point as the max.
      real(wp), intent(in) :: temp_max
         !! Explicit temperature ceiling (degC), used only when
         !! `has_temp_max`.

      real(wp) :: ftop_in, fbot_in, htg, fb, t_fr, enth_fp
      real(wp) :: max_temp, max_enth, enth_in, extra_enth
      real(wp) :: en_j, aa, bb, cc, disc, dt_denth
      real(wp) :: denom, dflux_dtot_dt

      ftop_in = ftop
      fbot_in = fbot
      htg = (ht_body + ftop_in) - fbot_in
      fb = -(dftop_dt - dfbot_dt)

      extra_heat = 0.0_wp
      extra_enth = 0.0_wp
      if (sice > 0.0_wp) then
         t_fr = ice_t_freeze(sice)
         enth_fp = ice_enthalpy_liquid_freeze(sice)
      else
         t_fr = 0.0_wp
         enth_fp = ice_enth_from_ts(0.0_wp, 0.0_wp)
      end if

      max_temp = t_fr
      max_enth = enth_fp
      if (has_temp_max) then
         if (temp_max < t_fr) then
            max_temp = temp_max
            max_enth = ice_enth_from_ts(temp_max, sice)
         end if
      end if

      enth_in = enth

      if (m_lay == 0.0_wp) then
         new_temp = min(htg/fb, max_temp)
         enth = ice_enth_from_ts(new_temp, sice)
      else if (dtt*(htg - fb*max_temp) >= m_lay*(max_enth - enth_in)) then
         ! Heat applied would push the layer above max_temp -> pin and
         ! bank the excess heat.
         extra_enth = m_lay*(enth_in - max_enth) + dtt*(htg - fb*max_temp)
         extra_heat = extra_enth
         new_temp = max_temp
         enth = max_enth
      else if (sice == 0.0_wp) then
         dt_denth = 1.0_wp/ICE_CP_ICE
         enth = enth_fp + (dtt*htg + m_lay*(enth_in - enth_fp))/ &
                (m_lay + dtt*(fb*dt_denth))
         new_temp = dt_denth*((dtt*htg + m_lay*(enth_in - enth_fp))/ &
                              (m_lay + dtt*(fb*dt_denth)))
      else
         en_j = enth_in - ice_enthalpy_liquid(0.0_wp, 0.0_wp)
         aa = m_lay*ICE_CP_ICE + fb*dtt
         bb = -(m_lay*((en_j - (ICE_CP_WATER - ICE_CP_ICE)*t_fr) + ICE_LAT_FUS) + htg*dtt)
         cc = m_lay*ICE_LAT_FUS*t_fr
         disc = max(bb*bb - 4.0_wp*aa*cc, 0.0_wp)
         if (bb >= 0.0_wp) then
            new_temp = -(bb + sqrt(disc))/(2.0_wp*aa)
         else
            new_temp = (2.0_wp*cc)/(-bb + sqrt(disc))
         end if
         ! Cp_ice == Cp_brine -> "keep this solution" (SIS2:837).
         enth = ice_enth_from_ts(new_temp, sice)
      end if

      ! Decide explicit vs. conservation-inverted flux bookkeeping
      ! (SIS2:913-941; prototype :117-133).
      if (abs(hf_err_rat*dftop_dt) <= m_lay) then
         ftop = ftop_in + dftop_dt*new_temp
         if (hf_err_rat*dfbot_dt <= m_lay) then
            fbot = fbot_in + dfbot_dt*new_temp
         else
            fbot = (ht_body + ftop) - (m_lay*(enth - enth_in) + extra_enth)/dtt
         end if
      else if (hf_err_rat*dfbot_dt <= m_lay) then
         fbot = fbot_in + dfbot_dt*new_temp
         ftop = (fbot - ht_body) + (m_lay*(enth - enth_in) + extra_enth)/dtt
      else
         denom = dfbot_dt - dftop_dt
         if (denom > 0.0_wp) then
            dflux_dtot_dt = (htg - (m_lay*(enth - enth_in) + extra_enth)/dtt)/denom
         else
            dflux_dtot_dt = 0.0_wp
         end if
         ftop = ftop_in + dftop_dt*dflux_dtot_dt
         fbot = fbot_in + dfbot_dt*dflux_dtot_dt
      end if
   end subroutine update_lay_enth