Skip to content

Commit

Permalink
Alternative congelation formulation following Plante et al. 2024 (CIC…
Browse files Browse the repository at this point in the history 
…E-Consortium#494)

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.
  • Loading branch information
@eclare108213 @dabail10
eclare108213 authored and dabail10 committed Aug 19, 2024
1 parent 064efdd commit 6d902e5
Show file tree
Hide file tree
Showing 10 changed files with 102 additions and 41 deletions.
52 changes: 35 additions & 17 deletions columnphysics/icepack_parameters.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,6 +161,11 @@ module icepack_parameters
! only for use with tr_aero or tr_zaero
conserv_check = .false. ! if true, do conservations checks and abort

character(len=char_len), public :: &
congel_freeze = 'two-step' ! congelation computation
! 'two-step' = original formulation
! 'one-step' = Plante et al, The Cryosphere, 2024

character(len=char_len), public :: &
tfrz_option = 'mushy' ! form of ocean freezing temperature
! 'minus1p8' = -1.8 C
Expand Down Expand Up @@ -476,7 +481,7 @@ subroutine icepack_init_parameters( &
atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
atmiter_conv_in, calc_dragio_in, &
tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, &
apnd_sl_in, saltflux_option_in, &
saltflux_option_in, congel_freeze_in, apnd_sl_in, &
floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, &
dpscale_in, rfracmin_in, rfracmax_in, pndaspect_in, hs1_in, hp1_in, &
bgc_flux_type_in, z_tracers_in, scale_bgc_in, solve_zbgc_in, &
Expand Down Expand Up @@ -586,15 +591,20 @@ subroutine icepack_init_parameters( &
phi_i_mushy_in ! liquid fraction of congelation ice

character(len=*), intent(in), optional :: &
tfrz_option_in ! form of ocean freezing temperature
! 'minus1p8' = -1.8 C
! 'linear_salt' = -depressT * sss
! 'mushy' conforms with ktherm=2
congel_freeze_in ! congelation computation
! 'two-step' = original formulation
! 'one-step' = Plante et al, The Cryosphere, 2024

character(len=*), intent(in), optional :: &
saltflux_option_in ! Salt flux computation
! 'constant' reference value of ice_ref_salinity
! 'prognostic' prognostic salt flux
tfrz_option_in ! form of ocean freezing temperature
! 'minus1p8' = -1.8 C
! 'linear_salt' = -depressT * sss
! 'mushy' conforms with ktherm=2

character(len=*), intent(in), optional :: &
saltflux_option_in ! Salt flux computation
! 'constant' reference value of ice_ref_salinity
! 'prognostic' prognostic salt flux

!-----------------------------------------------------------------------
! Parameters for radiation
Expand Down Expand Up @@ -964,6 +974,7 @@ subroutine icepack_init_parameters( &
if (present(highfreq_in) ) highfreq = highfreq_in
if (present(natmiter_in) ) natmiter = natmiter_in
if (present(atmiter_conv_in) ) atmiter_conv = atmiter_conv_in
if (present(congel_freeze_in) ) congel_freeze = congel_freeze_in
if (present(tfrz_option_in) ) tfrz_option = tfrz_option_in
if (present(saltflux_option_in) ) saltflux_option = saltflux_option_in
if (present(kitd_in) ) kitd = kitd_in
Expand Down Expand Up @@ -1208,7 +1219,7 @@ subroutine icepack_query_parameters( &
atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
atmiter_conv_out, calc_dragio_out, &
tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, &
apnd_sl_out, saltflux_option_out, &
saltflux_option_out, congel_freeze_out, apnds_sl_out, &
floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, &
dpscale_out, rfracmin_out, rfracmax_out, pndaspect_out, hs1_out, hp1_out, &
bgc_flux_type_out, z_tracers_out, scale_bgc_out, solve_zbgc_out, &
Expand Down Expand Up @@ -1327,16 +1338,21 @@ subroutine icepack_query_parameters( &
phi_i_mushy_out ! liquid fraction of congelation ice

character(len=*), intent(out), optional :: &
tfrz_option_out ! form of ocean freezing temperature
! 'minus1p8' = -1.8 C
! 'constant' = Tocnfrz
! 'linear_salt' = -depressT * sss
! 'mushy' conforms with ktherm=2
congel_freeze_out ! congelation computation
! 'two-step' = original formulation
! 'one-step' = Plante et al, The Cryosphere, 2024

character(len=*), intent(out), optional :: &
tfrz_option_out ! form of ocean freezing temperature
! 'minus1p8' = -1.8 C
! 'constant' = Tocnfrz
! 'linear_salt' = -depressT * sss
! 'mushy' conforms with ktherm=2

character(len=*), intent(out), optional :: &
saltflux_option_out ! Salt flux computation
! 'constant' reference value of ice_ref_salinity
! 'prognostic' prognostic salt flux
saltflux_option_out ! Salt flux computation
! 'constant' reference value of ice_ref_salinity
! 'prognostic' prognostic salt flux


!-----------------------------------------------------------------------
Expand Down Expand Up @@ -1739,6 +1755,7 @@ subroutine icepack_query_parameters( &
if (present(highfreq_out) ) highfreq_out = highfreq
if (present(natmiter_out) ) natmiter_out = natmiter
if (present(atmiter_conv_out) ) atmiter_conv_out = atmiter_conv
if (present(congel_freeze_out) ) congel_freeze_out = congel_freeze
if (present(tfrz_option_out) ) tfrz_option_out = tfrz_option
if (present(saltflux_option_out) ) saltflux_option_out = saltflux_option
if (present(kitd_out) ) kitd_out = kitd
Expand Down Expand Up @@ -1947,6 +1964,7 @@ subroutine icepack_write_parameters(iounit)
write(iounit,*) " highfreq = ", highfreq
write(iounit,*) " natmiter = ", natmiter
write(iounit,*) " atmiter_conv = ", atmiter_conv
write(iounit,*) " congel_freeze = ", trim(congel_freeze)
write(iounit,*) " tfrz_option= ", trim(tfrz_option)
write(iounit,*) " saltflux_option = ", trim(saltflux_option)
write(iounit,*) " kitd = ", kitd
Expand Down
40 changes: 26 additions & 14 deletions columnphysics/icepack_therm_vertical.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -86,7 +87,7 @@ subroutine thermo_vertical (nilyr, nslyr, &
fsnow, fpond, &
fbot, Tbot, &
Tsnice, sss, &
rsnw, &
sst, rsnw, &
lhcoef, shcoef, &
fswsfc, fswint, &
Sswabs, Iswabs, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -1023,6 +1026,7 @@ subroutine thickness_changes (nilyr, nslyr, &
congel, snoice, &
mlt_onset, frz_onset,&
zSin, sss, &
sst, &
dsnow, rsnw)

integer (kind=int_kind), intent(in) :: &
Expand Down Expand Up @@ -1092,6 +1096,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
Expand Down Expand Up @@ -1141,9 +1146,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
Expand Down Expand Up @@ -1264,15 +1269,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

Expand Down Expand Up @@ -2683,7 +2695,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), &
Expand Down
9 changes: 6 additions & 3 deletions configuration/driver/icedrv_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ subroutine input_data
natmiter, kitd, kcatbound

character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
cpl_frazil, tfrz_option, saltflux_option, &
cpl_frazil, congel_freeze, tfrz_option, saltflux_option, &
frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table

logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain, update_ocn_f
Expand Down Expand Up @@ -176,12 +176,12 @@ subroutine input_data
update_ocn_f, l_mpond_fresh, ustar_min, &
fbot_xfer_type, oceanmixed_ice, emissivity, &
formdrag, highfreq, natmiter, &
atmiter_conv, calc_dragio, &
atmiter_conv, calc_dragio, congel_freeze, &
tfrz_option, saltflux_option, ice_ref_salinity, &
default_season, wave_spec_type, cpl_frazil, &
precip_units, fyear_init, ycycle, &
atm_data_type, ocn_data_type, bgc_data_type, &
lateral_flux_type, &
lateral_flux_type, &
atm_data_file, ocn_data_file, bgc_data_file, &
ice_data_file, &
atm_data_format, ocn_data_format, bgc_data_format, &
Expand Down Expand Up @@ -229,6 +229,7 @@ subroutine input_data
dSdt_slow_mode_out=dSdt_slow_mode, &
phi_c_slow_mode_out=phi_c_slow_mode, Tliquidus_max_out=Tliquidus_max, &
phi_i_mushy_out=phi_i_mushy, conserv_check_out=conserv_check, &
congel_freeze_out=congel_freeze, &
tfrz_option_out=tfrz_option, saltflux_option_out=saltflux_option, &
ice_ref_salinity_out=ice_ref_salinity, kalg_out=kalg, &
fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, &
Expand Down Expand Up @@ -801,6 +802,7 @@ subroutine input_data
write(nu_diag,1005) ' hi_min = ', hi_min
write(nu_diag,1030) ' fbot_xfer_type = ', trim(fbot_xfer_type)
write(nu_diag,1010) ' oceanmixed_ice = ', oceanmixed_ice
write(nu_diag,1030) ' congel_freeze = ', trim(congel_freeze)
write(nu_diag,1030) ' tfrz_option = ', trim(tfrz_option)
write(nu_diag,*) ' saltflux_option = ', &
trim(saltflux_option)
Expand Down Expand Up @@ -987,6 +989,7 @@ subroutine input_data
dSdt_slow_mode_in=dSdt_slow_mode, &
phi_c_slow_mode_in=phi_c_slow_mode, Tliquidus_max_in=Tliquidus_max, &
phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, &
congel_freeze_in=congel_freeze, &
tfrz_option_in=tfrz_option, saltflux_option_in=saltflux_option, &
ice_ref_salinity_in=ice_ref_salinity, kalg_in=kalg, &
fbot_xfer_type_in=fbot_xfer_type, &
Expand Down
1 change: 1 addition & 0 deletions configuration/scripts/icepack_in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,7 @@
cpl_frazil = 'fresh_ice_correction'
update_ocn_f = .false.
l_mpond_fresh = .false.
congel_freeze = 'two-step'
tfrz_option = 'mushy'
ice_ref_salinity = 4.0
saltflux_option = 'constant'
Expand Down
2 changes: 1 addition & 1 deletion configuration/scripts/machines/env.chicoma_intel
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ setenv ICE_MACHINE_WKDIR /lustre/scratch5/$user/ICEPACK_RUNS
setenv ICE_MACHINE_INPUTDATA /usr/projects/climate/eclare/DATA/Consortium
setenv ICE_MACHINE_BASELINE /lustre/scratch5/$user/ICEPACK_BASELINE
setenv ICE_MACHINE_SUBMIT "sbatch "
setenv ICE_MACHINE_ACCT t23_cice
setenv ICE_MACHINE_ACCT t24_cice
setenv ICE_MACHINE_QUEUE "debug"
setenv ICE_MACHINE_TPNODE 128 # tasks per node
setenv ICE_MACHINE_BLDTHRDS 12
Expand Down
1 change: 1 addition & 0 deletions configuration/scripts/options/set_nml.congel
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
congel_freeze = 'one-step'
2 changes: 2 additions & 0 deletions configuration/scripts/tests/base_suite.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ smoke col 1x1 debug,run1year,fluxopenw
smoke col 1x1 debug,run1year,dyn,fluxopenw
smoke col 1x1 debug,run1year,debug,snicartest
smoke col 1x1 debug,run1year,debug,snicar
smoke col 1x1 debug,run1year,debug,congel
restart col 1x1 debug
restart col 1x1 diag1
restart col 1x1 pondlvl
Expand All @@ -43,4 +44,5 @@ restart col 1x1 modal
restart col 1x1 fluxopenw
restart col 1x1 snicartest
restart col 1x1 snicar
restart col 1x1 congel

16 changes: 13 additions & 3 deletions doc/source/master_list.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ @string{JGR
@string{JGRO = {J. Geophys. Res. Oceans}}
@string{JGRB = {J. Geophys. Res. Biogeo.}}
@string{JGRA = {J. Geophys. Res. Atmos.}}
@string{JCT = {J. Comput. Phys.}}
@string{TC = {The Cryosphere}}
@string{QJRMS = {Quart. J. Royal Met. Soc.}}
@string{GRL = {Geophys. Res. Lett.}}
@string{JAS = {J. Atmos. Sci.}}
Expand Down Expand Up @@ -244,7 +244,7 @@ @Article{Maykut95
@Article{Murray96
author = "R.J. Murray",
title = "{Explicit generation of orthogonal grids for ocean models}",
journal = JCT,
journal = JCP,
year = {1996},
volume = {126},
pages = {251-273},
Expand Down Expand Up @@ -805,13 +805,23 @@ @article{Roberts19
@article{Dang19,
author = {Dang, C. and Zender, C. S. and Flanner, M. G.},
title = {Intercomparison and improvement of two-stream shortwave radiative transfer schemes in Earth system models for a unified treatment of cryospheric surfaces},
journal = {The Cryosphere},
journal = {TC},
volume = 13,
pages = {2325--2343},
doi = {10.5194/tc-13-2325-2019},
year = {2019}
}

@article{Plante24,
author = {Plante, M. and J.-F. Lemieux and L. B. Tremblay and A. Tivy and J. Angnatok and F. Roy and G. Smith and F. Dupont and A. K. Turner},
title = {Using Icepack to reproduce ice mass balance buoy observations in landfast ice: improvements from the mushy-layer thermodynamcis},
journal = {TC},
volume = 18,
pages = {1685--1708},
doi = {10.5194/tc-18-1685-2024},
year = {20124}
}

% **********************************************
% For new entries, see example entry in BIB_TEMPLATE.txt
% **********************************************
18 changes: 15 additions & 3 deletions doc/source/science_guide/sg_thermo.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1737,9 +1737,21 @@ conductive heat flux at the bottom surface:
If ice is melting at the bottom surface, :math:`q`
in Equation :eq:`bottom-melting` is the enthalpy of the bottom ice layer. If
ice is growing, :math:`q` is the enthalpy of new ice with temperature
:math:`T_f` and salinity :math:`S_{max}` (``ktherm`` = 1) or ocean surface
salinity (``ktherm`` = 2). This ice is added to the bottom layer.
ice is growing, :math:`q` is the enthalpy of new ice added to the bottom layer with
temperature :math:`T_f` and salinity :math:`S_{max}` (``ktherm`` = 1) for the Bitz99
:cite:`Bitz99` formulation. The original mushy thermodynamics
(``ktherm`` = 2, ``congel_freeze`` = 'two-step') formed new ice in two steps, first
moving the lower ice boundary into the ocean to form a mushy layer with an initial
liquid fraction :math:`phi_{init}`, then freezing the new ice in the next step. In
this case, the freezing temperature is calculated using the ocean surface salinity
:math:`SSS`. A second mushy thermo method
(``ktherm`` = 2, ``congel_freeze`` = 'one-step') freezes the ice immediately using
the brine salinity :math:`phi_{init} * SSS` :cite:`Plante24`.
In the two-step method, enthalpy not
yet used to freeze the ice is returned to the ocean, causing frazil ice to form
instead. Together, the congelation and frazil ice add up to a similar total amount
of new ice, but the differing freezing mechanism complicates comparisons with observational
data and there is an unnecessary lag in ice-ocean coupling.

In general, frazil ice formed in the ocean is added to the thinnest ice
category. The new ice is grown in the open water area of the grid cell
Expand Down
Loading

0 comments on commit 6d902e5

Please sign in to comment.