Skip to content

Zap residual ice #543

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Nov 21, 2025
Merged

Zap residual ice #543

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
36e5e00
Add area and mass minimum parameters and zap residual ice
eclare108213 Sep 25, 2025
8cbe6ca
change names of minimum area and mass parameters
eclare108213 Nov 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion columnphysics/icepack_flux.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,10 +75,11 @@ subroutine merge_fluxes (aicen, &
dpnd_initial, dpnd_initialn, &
dpnd_dlid, dpnd_dlidn)

! single category fluxes
! concentration is aicen_init in call to subroutine
real (kind=dbl_kind), intent(in) :: &
aicen ! concentration of ice

! single category fluxes
real (kind=dbl_kind), optional, intent(in) :: &
flw , & ! downward longwave flux (W/m**2)
strairxn, & ! air/ice zonal strss, (N/m**2)
Expand Down
25 changes: 22 additions & 3 deletions columnphysics/icepack_itd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ module icepack_itd
use icepack_tracers, only: tr_pond, tr_pond_lvl, nt_alvl
use icepack_tracers, only: n_iso, tr_iso, nt_smice, nt_rsnw, nt_rhos, nt_sice
use icepack_tracers, only: icepack_compute_tracers
use icepack_parameters, only: skl_bgc, z_tracers, hi_min
use icepack_parameters, only: skl_bgc, z_tracers, hi_min, itd_area_min, itd_mass_min
use icepack_parameters, only: kcatbound, kitd, saltflux_option, snwgrain, snwredist
use icepack_therm_shared, only: Tmin
use icepack_warnings, only: warnstr, icepack_warnings_add
Expand Down Expand Up @@ -1074,13 +1074,33 @@ subroutine zap_small_areas (dt, &
n, k, it, & !counting indices
blevels

logical (kind=log_kind) :: &
zap_residual, &
zap_category(ncat)

real (kind=dbl_kind) :: xtmp, sicen ! temporary variables
real (kind=dbl_kind) :: dvssl, dvint ! temporary variables
real (kind=dbl_kind) , dimension (1):: trcr_skl
real (kind=dbl_kind) , dimension (nblyr+1):: bvol

character(len=*),parameter :: subname='(zap_small_areas)'

!-----------------------------------------------------------------
! Flag categories with very small areas and residual ice
!-----------------------------------------------------------------

zap_category(:) = .false.
do n = 1, ncat
if ( abs(aicen(n)) <= puny .and. &
(abs(aicen(n)) /= c0 .or. abs(vicen(n)) /= c0 .or. abs(vsnon(n)) /= c0)) then
zap_category(n) = .true.
endif
enddo

zap_residual = .false.
if (aice < max(itd_area_min, puny) .or. &
aice*rhoi < max(itd_mass_min, puny)) zap_residual = .true. ! all categories

!-----------------------------------------------------------------
! I. Zap categories with very small areas.
!-----------------------------------------------------------------
Expand All @@ -1095,8 +1115,7 @@ subroutine zap_small_areas (dt, &
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//' Zap ice: negative ice area')
return
elseif (abs(aicen(n)) <= puny .and. &
(abs(aicen(n)) /= c0 .or. abs(vicen(n)) /= c0 .or. abs(vsnon(n)) /= c0)) then
elseif (zap_category(n) .or. zap_residual) then

!-----------------------------------------------------------------
! Account for tracers important for conservation
Expand Down
64 changes: 38 additions & 26 deletions columnphysics/icepack_parameters.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -254,17 +254,19 @@ module icepack_parameters
! 1 for exponential redistribution function

real (kind=dbl_kind), public :: &
Cf = 17._dbl_kind ,&! ratio of ridging work to PE change in ridging
Pstar = 2.75e4_dbl_kind ,&! constant in Hibler strength formula
! (kstrength = 0)
Cstar = 20._dbl_kind ,&! constant in Hibler strength formula
! (kstrength = 0)
dragio = 0.00536_dbl_kind ,&! ice-ocn drag coefficient
itd_area_min = 1.e-11_dbl_kind ,&! zap residual ice below a minimum area
itd_mass_min = 1.e-10_dbl_kind ,&! zap residual ice below a minimum mass
Cf = 17._dbl_kind ,&! ratio of ridging work to PE change in ridging
Pstar = 2.75e4_dbl_kind ,&! constant in Hibler strength formula
! (kstrength = 0)
Cstar = 20._dbl_kind ,&! constant in Hibler strength formula
! (kstrength = 0)
dragio = 0.00536_dbl_kind ,&! ice-ocn drag coefficient
thickness_ocn_layer1 = 2.0_dbl_kind,&! thickness of first ocean level (m)
iceruf_ocn = 0.03_dbl_kind ,&! under-ice roughness (m)
gravit = 9.80616_dbl_kind ,&! gravitational acceleration (m/s^2)
mu_rdg = 3.0_dbl_kind ! e-folding scale of ridged ice, krdg_partic=1 (m^0.5)
! (krdg_redist = 1)
iceruf_ocn = 0.03_dbl_kind ,&! under-ice roughness (m)
gravit = 9.80616_dbl_kind ,&! gravitational acceleration (m/s^2)
mu_rdg = 3.0_dbl_kind ! e-folding scale of ridged ice, krdg_partic=1 (m^0.5)
! (krdg_redist = 1)

logical (kind=log_kind), public :: &
calc_dragio = .false. ! if true, calculate dragio from iceruf_ocn and thickness_ocn_layer1
Expand Down Expand Up @@ -571,7 +573,7 @@ subroutine icepack_init_parameters( &
stefan_boltzmann_in, ice_ref_salinity_in, &
Tffresh_in, Lsub_in, Lvap_in, Timelt_in, Tsmelt_in, &
iceruf_in, Cf_in, Pstar_in, Cstar_in, kappav_in, &
kice_in, ksno_in, &
kice_in, ksno_in, itd_area_min_in, itd_mass_min_in, &
zref_in, hs_min_in, snowpatch_in, rhosi_in, sk_l_in, &
saltmax_in, phi_init_in, min_salin_in, Tliquidus_max_in, &
min_bgc_in, dSin0_frazil_in, hi_ssl_in, hs_ssl_in, hs_ssl_min_in, &
Expand Down Expand Up @@ -801,14 +803,16 @@ subroutine icepack_init_parameters( &
!-----------------------------------------------------------------------

real(kind=dbl_kind), intent(in), optional :: &
Cf_in, & ! ratio of ridging work to PE change in ridging
Pstar_in, & ! constant in Hibler strength formula
Cstar_in, & ! constant in Hibler strength formula
dragio_in, & ! ice-ocn drag coefficient
itd_area_min_in, & ! zap residual ice below this minimum area
itd_mass_min_in, & ! zap residual ice below this minimum mass
Cf_in, & ! ratio of ridging work to PE change in ridging
Pstar_in, & ! constant in Hibler strength formula
Cstar_in, & ! constant in Hibler strength formula
dragio_in, & ! ice-ocn drag coefficient
thickness_ocn_layer1_in, & ! thickness of first ocean level (m)
iceruf_ocn_in, & ! under-ice roughness (m)
gravit_in, & ! gravitational acceleration (m/s^2)
iceruf_in ! ice surface roughness (m)
iceruf_ocn_in, & ! under-ice roughness (m)
gravit_in, & ! gravitational acceleration (m/s^2)
iceruf_in ! ice surface roughness (m)

integer (kind=int_kind), intent(in), optional :: & ! defined in namelist
kstrength_in , & ! 0 for simple Hibler (1979) formulation
Expand Down Expand Up @@ -1152,6 +1156,8 @@ subroutine icepack_init_parameters( &
if (present(Tsmelt_in) ) Tsmelt = Tsmelt_in
if (present(ice_ref_salinity_in) ) ice_ref_salinity = ice_ref_salinity_in
if (present(iceruf_in) ) iceruf = iceruf_in
if (present(itd_area_min_in) ) itd_area_min = itd_area_min_in
if (present(itd_mass_min_in) ) itd_mass_min = itd_mass_min_in
if (present(Cf_in) ) Cf = Cf_in
if (present(Pstar_in) ) Pstar = Pstar_in
if (present(Cstar_in) ) Cstar = Cstar_in
Expand Down Expand Up @@ -1580,7 +1586,7 @@ subroutine icepack_query_parameters( &
stefan_boltzmann_out, ice_ref_salinity_out, &
Tffresh_out, Lsub_out, Lvap_out, Timelt_out, Tsmelt_out, &
iceruf_out, Cf_out, Pstar_out, Cstar_out, kappav_out, &
kice_out, ksno_out, &
kice_out, ksno_out, itd_area_min_out, itd_mass_min_out, &
zref_out, hs_min_out, snowpatch_out, rhosi_out, sk_l_out, &
saltmax_out, phi_init_out, min_salin_out, Tliquidus_max_out, &
min_bgc_out, dSin0_frazil_out, hi_ssl_out, hs_ssl_out, hs_ssl_min_out, &
Expand Down Expand Up @@ -1819,14 +1825,16 @@ subroutine icepack_query_parameters( &
!-----------------------------------------------------------------------

real(kind=dbl_kind), intent(out), optional :: &
Cf_out, & ! ratio of ridging work to PE change in ridging
Pstar_out, & ! constant in Hibler strength formula
Cstar_out, & ! constant in Hibler strength formula
dragio_out, & ! ice-ocn drag coefficient
itd_area_min_out, & ! zap residual ice below this minimum area
itd_mass_min_out, & ! zap residual ice below this minimum mass
Cf_out, & ! ratio of ridging work to PE change in ridging
Pstar_out, & ! constant in Hibler strength formula
Cstar_out, & ! constant in Hibler strength formula
dragio_out, & ! ice-ocn drag coefficient
thickness_ocn_layer1_out, & ! thickness of first ocean level (m)
iceruf_ocn_out, & ! under-ice roughness (m)
gravit_out, & ! gravitational acceleration (m/s^2)
iceruf_out ! ice surface roughness (m)
iceruf_ocn_out, & ! under-ice roughness (m)
gravit_out, & ! gravitational acceleration (m/s^2)
iceruf_out ! ice surface roughness (m)

integer (kind=int_kind), intent(out), optional :: & ! defined in namelist
kstrength_out , & ! 0 for simple Hibler (1979) formulation
Expand Down Expand Up @@ -2203,6 +2211,8 @@ subroutine icepack_query_parameters( &
if (present(ice_ref_salinity_out) ) ice_ref_salinity_out = ice_ref_salinity
if (present(iceruf_out) ) iceruf_out = iceruf
if (present(Cf_out) ) Cf_out = Cf
if (present(itd_area_min_out) ) itd_area_min_out = itd_area_min
if (present(itd_mass_min_out) ) itd_mass_min_out = itd_mass_min
if (present(Pstar_out) ) Pstar_out = Pstar
if (present(Cstar_out) ) Cstar_out = Cstar
if (present(kappav_out) ) kappav_out = kappav
Expand Down Expand Up @@ -2502,6 +2512,8 @@ subroutine icepack_write_parameters(iounit)
write(iounit,*) " Tsmelt = ",Tsmelt
write(iounit,*) " ice_ref_salinity = ",ice_ref_salinity
write(iounit,*) " iceruf = ",iceruf
write(iounit,*) " itd_area_min = ",itd_area_min
write(iounit,*) " itd_mass_min = ",itd_mass_min
write(iounit,*) " Cf = ",Cf
write(iounit,*) " Pstar = ",Pstar
write(iounit,*) " Cstar = ",Cstar
Expand Down
13 changes: 9 additions & 4 deletions configuration/driver/icedrv_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ subroutine input_data
real (kind=dbl_kind) :: ustar_min, albicev, albicei, albsnowv, albsnowi, &
ahmax, R_ice, R_pnd, R_snw, dT_mlt, rsnw_mlt, ksno, hi_min, Tliquidus_max, &
mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, &
apnd_sl, tscale_pnd_drain, &
apnd_sl, tscale_pnd_drain, itd_area_min, itd_mass_min, &
a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, dSdt_slow_mode, &
phi_c_slow_mode, phi_i_mushy, kalg, emissivity, floediam, hfrazilmin, &
rsnw_fall, rsnw_tmax, rhosnew, rhosmin, rhosmax, &
Expand Down Expand Up @@ -150,7 +150,7 @@ subroutine input_data
dumpfreq, diagfreq, diag_file, cpl_bgc, &
conserv_check, history_format, &
hi_init_slab, hsno_init_slab, hbar_init_itd, hsno_init_itd, &
sst_init
sst_init, itd_area_min, itd_mass_min

namelist /grid_nml/ &
kcatbound
Expand All @@ -177,6 +177,7 @@ subroutine input_data
hs0, dpscale, frzpnd, &
rfracmin, rfracmax, pndaspect, hs1, &
hp1, apnd_sl

namelist /snow_nml/ &
snwredist, snwgrain, rsnw_fall, rsnw_tmax, &
rhosnew, rhosmin, rhosmax, snwlvlfac, &
Expand Down Expand Up @@ -257,7 +258,8 @@ subroutine input_data
rhosnew_out=rhosnew, rhosmin_out = rhosmin, rhosmax_out=rhosmax, &
windmin_out=windmin, drhosdwind_out=drhosdwind, snwlvlfac_out=snwlvlfac, &
snw_aging_table_out=snw_aging_table, snw_growth_wet_out=snw_growth_wet, &
drsnw_min_out=drsnw_min, snwliq_max_out=snwliq_max)
drsnw_min_out=drsnw_min, snwliq_max_out=snwliq_max, &
itd_area_min_out=itd_area_min, itd_mass_min_out=itd_mass_min)

call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
Expand Down Expand Up @@ -710,6 +712,8 @@ subroutine input_data
write(nu_diag,1005) ' hsno_init_itd = ', hsno_init_itd
write(nu_diag,1005) ' sst_init = ', sst_init
write(nu_diag,1010) ' conserv_check = ', conserv_check
write(nu_diag,1000) ' itd_area_min = ', itd_area_min
write(nu_diag,1000) ' itd_mass_min = ', itd_mass_min
write(nu_diag,1020) ' kitd = ', kitd
write(nu_diag,1020) ' kcatbound = ', kcatbound
write(nu_diag,1020) ' ndtd = ', ndtd
Expand Down Expand Up @@ -1045,7 +1049,8 @@ subroutine input_data
rhosnew_in=rhosnew, rhosmin_in=rhosmin, rhosmax_in=rhosmax, &
windmin_in=windmin, drhosdwind_in=drhosdwind, snwlvlfac_in=snwlvlfac, &
snw_growth_wet_in=snw_growth_wet, drsnw_min_in=drsnw_min, &
snwliq_max_in=snwliq_max)
snwliq_max_in=snwliq_max, itd_area_min_in=itd_area_min, &
itd_mass_min_in=itd_mass_min)
call icepack_init_tracer_sizes(ntrcr_in=ntrcr, &
ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, &
nfsd_in=nfsd, n_iso_in=n_iso, n_aero_in=n_aero)
Expand Down
2 changes: 2 additions & 0 deletions configuration/scripts/icepack_in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@
history_format = 'none'
cpl_bgc = .false.
conserv_check = .false.
itd_area_min = 1.e-11
itd_mass_min = 1.e-10
/

&grid_nml
Expand Down
2 changes: 2 additions & 0 deletions doc/source/icepack_index.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -253,6 +253,8 @@ section :ref:`tabnamelist`.
"istep0", ":math:`\bullet` number of steps taken in previous run", "0"
"istep1", "total number of steps at current time step", ""
"Iswabs", "shortwave radiation absorbed in ice layers", "W/m\ :math:`^2`"
"itd_area_min", ":math:`\bullet` zap residual ice below a minimum area", "1.e-11"
"itd_mass_min", ":math:`\bullet` zap residual ice below a minimum mass", "1.e-10
"**J**", "", ""
"**K**", "", ""
"kalg", ":math:`\bullet` absorption coefficient for algae", ""
Expand Down
24 changes: 15 additions & 9 deletions doc/source/science_guide/sg_thermo.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1966,6 +1966,16 @@ as the product of sensitivity studies to balance the climatological tendencies o
wave fracture and welding. So that results do not vary as the number or range of
floe size categories varies, we fix this scaling coefficient, c_weld.

If ``tr_fsd=false``, lateral melting is accomplished by multiplying the state variables by
:math:`1-r_{side}`, where :math:`r_{side}` is the fraction of ice melted
laterally :cite:`Maykut87,Steele92`, and adjusting the ice
energy and fluxes as appropriate. We assume a floe diameter of 300 m.

If ``tr_fsd=true``, lateral melting is accomplished using the :cite:`Maykut87`
lateral heat flux, but applied to the ice using the prognostic floe size distribution
as described in :cite:`Horvat15` and :cite:`Roach18`. Lateral melt modifies
the ITD and the FSD.

If the latent heat flux is negative (i.e., latent heat is transferred
from the ice to the atmosphere), snow or snow-free ice sublimates at the
top surface. If the latent heat flux is positive, vapor from the
Expand Down Expand Up @@ -2004,15 +2014,11 @@ old and new layers, respectively. The enthalpies of the new layers are
.. math::
q_k = \frac{1}{\Delta h_i} \sum_{m=1}^{N_i} \eta_{km} q_m.

If ``tr_fsd=false``, lateral melting is accomplished by multiplying the state variables by
:math:`1-r_{side}`, where :math:`r_{side}` is the fraction of ice melted
laterally :cite:`Maykut87,Steele92`, and adjusting the ice
energy and fluxes as appropriate. We assume a floe diameter of 300 m.

If ``tr_fsd=true``, lateral melting is accomplished using the :cite:`Maykut87`
lateral heat flux, but applied to the ice using the prognostic floe size distribution
as described in :cite:`Horvat15` and :cite:`Roach18`. Lateral melt modifies
the ITD and the FSD.
Residual amounts of ice may be conservatively removed following the thermodynamics
and ridging calculations based on minimum area and mass parameters ``itd_area_min`` and
``itd_mass_min``. Initializing these parameters to CICE's ``dyn_area_min`` and ``dyn_mass_min``
namelist values ensures consistency between Icepack's thermodynamic and CICE's
dynamic calculations and avoids tiny amounts of residual ice in the solution.

Snow-ice formation
------------------
Expand Down
7 changes: 4 additions & 3 deletions doc/source/user_guide/lg_overview.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,10 @@ Icepack includes options for simulating sea ice thermodynamics, mechanical redis
tracers, including thickness, enthalpy, ice age, first-year ice area, deformed ice area and
volume, melt ponds, and biogeochemistry.

Icepack is called on a grid point by grid point basis. All data is passed in and out of the model
via subroutine interfaces. Fortran "use" statements are not encouraged for accessing data inside
the Icepack model.
Icepack is called on a grid point by grid point basis. All data is passed in and out of the model
via subroutine interfaces. Fortran "use" statements are not encouraged for accessing data inside
the Icepack model from outside. Parameters may be queried and/or reset from Icepack's default
values by calling the query and init interfaces as described in :ref:`sequence_and_interface`.

Icepack does not generally contain any parallelization or I/O. The driver of Icepack is
expected to support
Expand Down
2 changes: 2 additions & 0 deletions doc/source/user_guide/ug_case_settings.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,8 @@ setup_nml
"", "``none``", "no ice", ""
"", "'path/file'", "restart file name", ""
"``istep0``", "integer", "initial time step number", "0"
"``itd_area_min``", "real", "area below which ice is zapped", "1.e-11"
"``itd_mass_min``", "real", "mass below which ice is zapped", "1.e-10"
"``ndtd``", "integer", "number of dynamics/advection/ridging/steps per thermo timestep", "1"
"``npt``", "integer", "total number of time steps to take", "99999"
"``restart``", "logical", "initialize using restart file", "``.false.``"
Expand Down