From 0e7e470ad25c4a54e8e55d8ec966c7de5dba934b Mon Sep 17 00:00:00 2001 From: Elizabeth Hunke Date: Wed, 7 Aug 2024 12:11:20 -0600 Subject: [PATCH] Alternative congelation formulation following Plante et al. 2024 (#494) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add congel_freeze namelist and add 'one-step' option to the default 'two-step' option. Namelist flag congel_freeze chooses which formulation to use. The original is ‘two-step’, since only the mushy boundary moves in the first step and the phase change happens in the next step. Plante et al. (‘one-step’) moves the boundary and performs the phase change in a single time step. Maintain 'two-step' as default for now. --- columnphysics/icepack_therm_vertical.F90 | 40 +++++++++++++++--------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 7414a060..1cf32ef2 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -27,7 +27,7 @@ module icepack_therm_vertical use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair use icepack_parameters, only: rfracmin, rfracmax, dpscale, frzpnd, snwgrain, snwlvlfac use icepack_parameters, only: phi_i_mushy, floeshape, floediam, use_smliq_pnd, snwredist - use icepack_parameters, only: saltflux_option + use icepack_parameters, only: saltflux_option, congel_freeze use icepack_parameters, only: icepack_chkoptargflag use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso @@ -46,6 +46,7 @@ module icepack_therm_vertical use icepack_mushy_physics, only: icepack_mushy_temperature_mush use icepack_mushy_physics, only: liquidus_temperature_mush use icepack_mushy_physics, only: icepack_enthalpy_mush, enthalpy_of_melting + use icepack_mushy_physics, only: enthalpy_mush_liquid_fraction, enthalpy_brine use icepack_aerosol, only: update_aerosol use icepack_isotope, only: update_isotope @@ -86,7 +87,7 @@ subroutine thermo_vertical (nilyr, nslyr, & fsnow, fpond, & fbot, Tbot, & Tsnice, sss, & - rsnw, & + sst, rsnw, & lhcoef, shcoef, & fswsfc, fswint, & Sswabs, Iswabs, & @@ -168,6 +169,7 @@ subroutine thermo_vertical (nilyr, nslyr, & real (kind=dbl_kind), intent(in) :: & fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) Tbot , & ! ice bottom surface temperature (deg C) + sst , & ! sea surface temperature (C) sss ! ocean salinity ! coupler fluxes to atmosphere @@ -413,6 +415,7 @@ subroutine thermo_vertical (nilyr, nslyr, & congel, snoice, & mlt_onset, frz_onset, & zSin, sss, & + sst, & dsnow, rsnw) if (icepack_warnings_aborted(subname)) return @@ -1094,6 +1097,7 @@ subroutine thickness_changes (nilyr, nslyr, & congel, snoice, & mlt_onset, frz_onset,& zSin, sss, & + sst, & dsnow, rsnw) integer (kind=int_kind), intent(in) :: & @@ -1163,6 +1167,7 @@ subroutine thickness_changes (nilyr, nslyr, & zSin ! ice layer salinity (ppt) real (kind=dbl_kind), intent(in) :: & + sst , & ! sea surface temperature (C) sss ! ocean salinity (PSU) ! local variables @@ -1212,9 +1217,9 @@ subroutine thickness_changes (nilyr, nslyr, & qmlt ! enthalpy of melted ice (J m-3) = zero in BL99 formulation real (kind=dbl_kind) :: & - qbotm , & - qbotp , & - qbot0 , & + qbotm , & ! enthalpy of newly formed congelation ice + qbotp , & ! enthalpy needed to grow new congelation ice (includes ocean enthalpy) + qbotw , & ! enthalpy transferred to ocean during congelation freezing mass , & ! total mass from snow density tracers (kg/m^2) massi , & ! ice mass change factor tmp1 ! temporary scalar @@ -1335,15 +1340,22 @@ subroutine thickness_changes (nilyr, nslyr, & if (ktherm == 2) then - qbotm = icepack_enthalpy_mush(Tbot, sss) - qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy) - qbot0 = qbotm - qbotp - - dhi = ebot_gro / qbotp ! dhi > 0 - + if (congel_freeze == 'one-step') then + ! Plante et al., The Cryosphere, 18, 1685-1708, 2024 + qbotm = enthalpy_mush_liquid_fraction(Tbot, phi_i_mushy) + qbotw = enthalpy_brine(sst) + qbotp = qbotm - qbotw + dhi = ebot_gro / qbotp ! dhi > 0 + hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss*phi_i_mushy + else ! two-step + qbotm = icepack_enthalpy_mush(Tbot, sss) + qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy) + qbotw = qbotm - qbotp + dhi = ebot_gro / qbotp ! dhi > 0 + hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss + endif hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm - hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss - emlt_ocn = emlt_ocn - qbot0 * dhi + emlt_ocn = emlt_ocn - qbotw * dhi else @@ -2760,7 +2772,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, & fsnow=fsnow, fpond=fpond, & fbot=fbot, Tbot=Tbot, & Tsnice=Tsnice, sss=sss, & - rsnw=rsnw, & + sst=sst, rsnw=rsnw, & lhcoef=lhcoef, shcoef=shcoef, & fswsfc=fswsfcn (n), fswint=fswintn (n), & Sswabs=Sswabsn(:,n), Iswabs=Iswabsn (:,n), &