laytemp_sis2 Function

public pure function laytemp_sis2(m, t_fr, qf, bf, tp, dtt) result(new_temp)

Per-layer implicit heat-budget solve for the new layer temperature — SIS2 laytemp_SIS2 (SIS2_ice_thm.F90:544-700), ICE_CP_BRINE == ICE_CP_ICE closed-form branches only (the Newton/false-position refinement at :625-692 is dead code under the simplification and is deliberately NOT ported). Port of prototype sis2_column.py:21-55.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: m

Layer mass (kg/m²).
real(kind=wp), intent(in) :: t_fr

Layer freezing temperature (degC); 0 for snow/fresh water.
real(kind=wp), intent(in) :: qf

Forcing heat flux into the layer (W/m²).
real(kind=wp), intent(in) :: bf

Implicit coupling coefficient to the neighbour temperature (W/m²/K).
real(kind=wp), intent(in) :: tp

Previous-step layer temperature (degC).
real(kind=wp), intent(in) :: dtt

Timestep (s).

Return Value real(kind=wp)


Called by

proc~~laytemp_sis2~~CalledByGraph proc~laytemp_sis2 laytemp_sis2 proc~ice_temp_sis2 ice_temp_sis2 proc~ice_temp_sis2->proc~laytemp_sis2 proc~ice_column_step ice_column_step proc~ice_column_step->proc~ice_temp_sis2 proc~ice_thermo_columns ice_thermo_columns proc~ice_thermo_columns->proc~ice_column_step proc~ice_thermo_driver_step ice_thermo_driver_step proc~ice_thermo_driver_step->proc~ice_thermo_columns proc~engine_step_ice engine_step_ice proc~engine_step_ice->proc~ice_thermo_driver_step

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: aa
real(kind=wp), private :: bb
real(kind=wp), private :: cc
real(kind=wp), private :: disc
real(kind=wp), private :: e0

Source Code

   pure function laytemp_sis2(m, t_fr, qf, bf, tp, dtt) result(new_temp)
      !! Per-layer implicit heat-budget solve for the new layer
      !! temperature — SIS2 `laytemp_SIS2` (SIS2_ice_thm.F90:544-700),
      !! `ICE_CP_BRINE == ICE_CP_ICE` closed-form branches only (the
      !! Newton/false-position refinement at :625-692 is dead code
      !! under the simplification and is deliberately NOT ported).
      !! Port of prototype `sis2_column.py:21-55`.
      !$acc routine seq
      real(wp), intent(in) :: m
         !! Layer mass (kg/m²).
      real(wp), intent(in) :: t_fr
         !! Layer freezing temperature (degC); 0 for snow/fresh water.
      real(wp), intent(in) :: qf
         !! Forcing heat flux into the layer (W/m²).
      real(wp), intent(in) :: bf
         !! Implicit coupling coefficient to the neighbour temperature
         !! (W/m²/K).
      real(wp), intent(in) :: tp
         !! Previous-step layer temperature (degC).
      real(wp), intent(in) :: dtt
         !! Timestep (s).
      real(wp) :: new_temp

      real(wp) :: e0, aa, bb, cc, disc

      if (t_fr == 0.0_wp) then
         ! Fresh water / snow linear branch (SIS2:585-592).
         new_temp = (m*ICE_CP_ICE*tp + qf*dtt)/(m*ICE_CP_ICE + bf*dtt)
      else
         if (tp >= t_fr) then
            e0 = ICE_CP_WATER*(tp - t_fr)
         else
            ! (Cp_brine - Cp_ice) term vanishes under the simplification.
            e0 = ICE_CP_ICE*(tp - t_fr) - ICE_LAT_FUS*(1.0_wp - t_fr/tp)
         end if

         if (m*e0 + dtt*(qf - bf*t_fr) >= 0.0_wp) then
            ! Layer would be fully melted -> pin to freezing (SIS2:606).
            new_temp = t_fr
         else
            aa = m*ICE_CP_ICE + bf*dtt
            bb = -(m*((e0 + ICE_LAT_FUS) + ICE_CP_ICE*t_fr) + qf*dtt)
            cc = m*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 -> the quadratic root IS the final
            ! answer; the Newton/false-position loop is not ported.
         end if
      end if

      new_temp = min(new_temp, t_fr)
   end function laytemp_sis2