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
apcraig merged 2 commits into from
Nov 21, 2025
Merged

Zap residual ice #543

Show file tree
Hide file tree
Changes from 1 commit
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, min_area, min_mass
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(min_area, puny) .or. &
aice*rhoi < max(min_mass, puny)) zap_residual = .true. ! all categories
Copy link
Contributor

@anton-seaice anton-seaice Oct 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about aborting if min_area or min_mass is less than `puny? Those setings mean the user has tried to do something which the model doesn't support.

e.g. could implementing as max(min_area, puny) leave people confused about why there is no ice with concentration less than puny, greater than min_area

Copy link
Contributor Author

@eclare108213 eclare108213 Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess someone might try to set min_area and min_mass both equal to zero in an attempt to turn off zap_residual, which wouldn't work. We can certainly implement a flag to turn it off, as @NickSzapiro-NOAA suggests. Then under the covers we could do min_area = min(min_area, puny) at initialization to enforce the puny minimum, which I do think we should keep, although I'd rather not obscure what's happening. What do you suggest, @anton-seaice ?

Copy link
Contributor

@anton-seaice anton-seaice Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes the namelist flag to enable/disable zap_residual is ok - (I think it should be enabled by default though, but we could revisit the default in a year or two)

Does setting min_area < puny means no ice would be zapped ? (as thermodynamics wouldn't be active?)

I think abort if the user tries to set it less than puny, we could write in the error message to use the "zap_residual" flag to disable zapping if needed ?

Copy link
Contributor Author

@eclare108213 eclare108213 Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Setting min_area < puny means that ice areas less than puny would still be zapped, which (if I remember correctly) is the current behavior. In that case, we could write a warning that areas less than puny will still be zapped and offer the namelist flag as an alternative to turn off zapping completely, with no need to abort. I agree that having zap_residual on as the default is what we want.


!-----------------------------------------------------------------
! 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
16 changes: 14 additions & 2 deletions columnphysics/icepack_parameters.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -254,6 +254,8 @@ module icepack_parameters
! 1 for exponential redistribution function

real (kind=dbl_kind), public :: &
min_area = 1.e-11_dbl_kind ,&! zap residual ice below a minimum area
min_mass = 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)
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, min_area_in, min_mass_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,6 +803,8 @@ subroutine icepack_init_parameters( &
!-----------------------------------------------------------------------

real(kind=dbl_kind), intent(in), optional :: &
min_area_in, & ! zap residual ice below this minimum area
min_mass_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
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(min_area_in) ) min_area = min_area_in
if (present(min_mass_in) ) min_mass = min_mass_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, min_area_out, min_mass_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,6 +1825,8 @@ subroutine icepack_query_parameters( &
!-----------------------------------------------------------------------

real(kind=dbl_kind), intent(out), optional :: &
min_area_out, & ! zap residual ice below this minimum area
min_mass_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
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(min_area_out) ) min_area_out = min_area
if (present(min_mass_out) ) min_mass_out = min_mass
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,*) " min_area = ",min_area
write(iounit,*) " min_mass = ",min_mass
write(iounit,*) " Cf = ",Cf
write(iounit,*) " Pstar = ",Pstar
write(iounit,*) " Cstar = ",Cstar
Expand Down
12 changes: 8 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, min_area, min_mass, &
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, min_area, min_mass

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, &
min_area_out=min_area, min_mass_out=min_mass)

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) ' min_area = ', min_area
write(nu_diag,1000) ' min_mass = ', min_mass
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,7 @@ 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, min_area_in=min_area, min_mass_in=min_mass)
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.
min_area = 1.e-11
min_mass = 1.e-10
/

&grid_nml
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 ``min_area`` and
``min_mass``. 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"
"``min_area``", "real", "area below which ice is zapped", ""
"``min_mass``", "real", "mass below which ice is zapped", ""
"``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