change dynamics minimum area and mass values and add Icepack parameters to remove residual ice

cicecore/cicedyn/general/ice_init.F90

@@ -157,7 +157,7 @@ subroutine input_data
         phi_c_slow_mode, phi_i_mushy, kalg, atmiter_conv, Pstar, Cstar, &
         sw_frac, sw_dtemp, floediam, hfrazilmin, iceruf, iceruf_ocn, &
         rsnw_fall, rsnw_tmax, rhosnew, rhosmin, rhosmax, Tliquidus_max, &
-        windmin, drhosdwind, snwlvlfac, tscale_pnd_drain
+        windmin, drhosdwind, snwlvlfac, tscale_pnd_drain, itd_area_min, itd_mass_min
 
       integer (kind=int_kind) :: ktherm, kstrength, krdg_partic, krdg_redist, natmiter, &
         kitd, kcatbound, ktransport
@@ -178,6 +178,7 @@ subroutine input_data
 
       integer (kind=int_kind) :: rplvl, rptopo, rpsealvl
       real (kind=dbl_kind)    :: Cf, ksno, puny, ice_ref_salinity, Tocnfrz
+      real (kind=dbl_kind), parameter :: ice_init_spval = -999._dbl_kind
 
       character (len=char_len) :: abort_list
       character (len=char_len)      :: nml_name ! namelist name
@@ -203,7 +204,7 @@ subroutine input_data
         hist_time_axis,                                                 &
         print_global,   print_points,   latpnt,          lonpnt,        &
         debug_forcing,  histfreq,       histfreq_n,      hist_avg,      &
-        hist_suffix, history_deflate, history_chunksize,                &
+        hist_suffix,    history_deflate, history_chunksize,             &
         history_dir,    history_file,   history_precision, cpl_bgc,     &
         histfreq_base,  dumpfreq_base,  timer_stats,     memory_stats,  &
         conserv_check,  debug_model,    debug_model_step,               &
@@ -240,8 +241,7 @@ subroutine input_data
         a_rapid_mode,   Rac_rapid_mode,  aspect_rapid_mode,             &
         dSdt_slow_mode, phi_c_slow_mode, phi_i_mushy,                   &
         floediam,       hfrazilmin,      Tliquidus_max,   hi_min,       &
-        tscale_pnd_drain
-
+        itd_area_min,   itd_mass_min,    tscale_pnd_drain
 
       namelist /dynamics_nml/ &
         kdyn,           ndte,           revised_evp,    yield_curve,    &
@@ -419,8 +419,8 @@ subroutine input_data
       kstrength = 1           ! 1 = Rothrock 75 strength, 0 = Hibler 79
       Pstar = 2.75e4_dbl_kind ! constant in Hibler strength formula (kstrength = 0)
       Cstar = 20._dbl_kind    ! constant in Hibler strength formula (kstrength = 0)
-      dyn_area_min = p001     ! minimum ice area concentration to activate dynamics
-      dyn_mass_min = p01      ! minimum ice mass to activate dynamics (kg/m^2)
+      dyn_area_min = 1.e-11_dbl_kind ! minimum ice area concentration to activate dynamics
+      dyn_mass_min = 1.e-10_dbl_kind ! minimum ice mass to activate dynamics (kg/m^2)
       krdg_partic = 1         ! 1 = new participation, 0 = Thorndike et al 75
       krdg_redist = 1         ! 1 = new redistribution, 0 = Hibler 80
       mu_rdg = 3              ! e-folding scale of ridged ice, krdg_partic=1 (m^0.5)
@@ -487,6 +487,8 @@ subroutine input_data
       cpl_frazil = 'fresh_ice_correction' ! type of coupling for frazil ice
       ustar_min = 0.005        ! minimum friction velocity for ocean heat flux (m/s)
       hi_min = p01             ! minimum ice thickness allowed (m)
+      itd_area_min = ice_init_spval ! zap residual ice below a minimum area
+      itd_mass_min = ice_init_spval ! zap residual ice below a minimum mass
       iceruf = 0.0005_dbl_kind ! ice surface roughness at atmosphere interface (m)
       iceruf_ocn = 0.03_dbl_kind ! under-ice roughness (m)
       calc_dragio = .false.    ! compute dragio from iceruf_ocn and thickness of first ocean level
@@ -922,6 +924,24 @@ subroutine input_data
       if (trim(diag_type) == 'file') call get_fileunit(nu_diag)
 #endif
 
+      ! To remove small amounts of residual ice that are not handled by either dynamics or
+      ! column physics, the minimum area and mass parameters should be set to the same values
+      ! in both places. The default sets the column physics (itd) parameters to the dynamics
+      ! values (available in namelist). itd_area_min and itd_mass_min can be added to the
+      ! namelist file ice_in and set to different values, if desired. Setting them to
+      ! zero turns off residual zapping completely.
+      if (itd_area_min /= ice_init_spval .or. itd_mass_min /= ice_init_spval) then
+         ! allow itd and dyn parameters to be different
+         write(nu_diag,*) subname//' WARNING: zap_residual parameters are reset in namelist'
+      elseif (itd_area_min == c0 .or. itd_mass_min == c0) then
+         ! turn off residual zapping in Icepack
+         write(nu_diag,*) subname//' WARNING: zap_residual is turned off'
+      else
+         ! itd and dyn parameters are the same by default
+         itd_area_min = dyn_area_min ! zap residual ice below dynamics minimum area
+         itd_mass_min = dyn_mass_min ! zap residual ice below dynamics minimum mass
+      endif
+
       !-----------------------------------------------------------------
       ! broadcast namelist settings
       !-----------------------------------------------------------------
@@ -1145,6 +1165,8 @@ subroutine input_data
       call broadcast_scalar(l_mpond_fresh,        master_task)
       call broadcast_scalar(ustar_min,            master_task)
       call broadcast_scalar(hi_min,               master_task)
+      call broadcast_scalar(itd_area_min,         master_task)
+      call broadcast_scalar(itd_mass_min,         master_task)
       call broadcast_scalar(iceruf,               master_task)
       call broadcast_scalar(iceruf_ocn,           master_task)
       call broadcast_scalar(calc_dragio,          master_task)
@@ -2400,6 +2422,10 @@ subroutine input_data
          write(nu_diag,1030) ' fbot_xfer_type   = ', trim(fbot_xfer_type),trim(tmpstr2)
          write(nu_diag,1000) ' ustar_min        = ', ustar_min,' : minimum value of ocean friction velocity'
          write(nu_diag,1000) ' hi_min           = ', hi_min,' : minimum ice thickness allowed (m)'
+         write(nu_diag,1000) ' puny             = ', puny,' : general-use minimum value'
+         write(nu_diag,*) ' Ice thickness category areas smaller than puny are always removed.'
+         write(nu_diag,1003) ' itd_area_min     = ', itd_area_min,' : zap residual ice below a minimum area'
+         write(nu_diag,1003) ' itd_mass_min     = ', itd_mass_min,' : zap residual ice below a minimum mass'
          if (calc_dragio) then
             tmpstr2 = ' : dragio computed from iceruf_ocn'
          else
@@ -2410,8 +2436,11 @@ subroutine input_data
             write(nu_diag,1002) ' iceruf_ocn       = ', iceruf_ocn,' : under-ice roughness length'
          endif
 
+         write(nu_diag,*) ' '
+         write(nu_diag,*) ' Floe size distribution and waves'
+         write(nu_diag,*) '---------------------------------'
+         write(nu_diag,1002) ' floediam         = ', floediam, ' : constant floe diameter'
          if (tr_fsd) then
-            write(nu_diag,1002) ' floediam         = ', floediam, ' constant floe diameter'
             if (wave_spec) then
                tmpstr2 = ' : use wave spectrum for floe size distribution'
             else
@@ -2799,7 +2828,7 @@ subroutine input_data
          atmbndy_in=atmbndy, calc_strair_in=calc_strair, formdrag_in=formdrag, highfreq_in=highfreq, &
          kitd_in=kitd, kcatbound_in=kcatbound, hs0_in=hs0, dpscale_in=dpscale, frzpnd_in=frzpnd, &
          rfracmin_in=rfracmin, rfracmax_in=rfracmax, pndaspect_in=pndaspect, hs1_in=hs1, hp1_in=hp1, &
-         apnd_sl_in=apnd_sl, &
+         apnd_sl_in=apnd_sl, itd_area_min_in=itd_area_min, itd_mass_min_in=itd_mass_min, &
          ktherm_in=ktherm, calc_Tsfc_in=calc_Tsfc, conduct_in=conduct, semi_implicit_Tsfc_in=semi_implicit_Tsfc, &
          a_rapid_mode_in=a_rapid_mode, Rac_rapid_mode_in=Rac_rapid_mode, vapor_flux_correction_in=vapor_flux_correction, &
          floediam_in=floediam, hfrazilmin_in=hfrazilmin, Tliquidus_max_in=Tliquidus_max, &

configuration/scripts/ice_in

@@ -201,8 +201,8 @@
     reltol_pgmres   = 1e-6
     algo_nonlin     = 'picard'
     use_mean_vrel   = .true.
-    dyn_area_min    = 0.001d0
-    dyn_mass_min    = 0.01d0
+    dyn_area_min    = 1.e-11
+    dyn_mass_min    = 1.e-10
 /
 
 &shortwave_nml

configuration/scripts/options/set_nml.gridc

@@ -1,2 +1,4 @@
 grid_ice = 'C'
+dyn_area_min    = 0.001d0
+dyn_mass_min    = 0.01d0
 

configuration/scripts/options/set_nml.gridcd

@@ -1,4 +1,6 @@
 grid_ice = 'C_override_D'
+dyn_area_min    = 0.001d0
+dyn_mass_min    = 0.01d0
 # visc_method=avg_zeta causes some gridcd tests to abort, use avg_strength for now
 visc_method = 'avg_strength'
 

doc/source/cice_index.rst

@@ -120,7 +120,7 @@ section :ref:`tabnamelist`.
    "cosw", "cosine of the turning angle in water", "1."
    "coszen", "cosine of the zenith angle", ""
    "Cp", "proportionality constant for potential energy", "kg/m\ :math:`^2`/s\ :math:`^2`"
-   "cpl_frazil", ":math:`\bullet` type of frazil ice coupling", ""
+   "cpl_frazil", "type of frazil ice coupling", ""
    "cp_air", "specific heat of air", "1005.0 J/kg/K"
    "cp_ice", "specific heat of fresh ice", "2106. J/kg/K"
    "cp_ocn", "specific heat of sea water", "4218. J/kg/K"
@@ -208,8 +208,8 @@ section :ref:`tabnamelist`.
    "dvidtd", "ice volume tendency due to dynamics/transport", "m/s"
    "dvidtt", "ice volume tendency due to thermodynamics", "m/s"
    "dvirdg(n)dt", "ice volume ridging rate (category n)", "m/s"
-   "dyn_area_min", "minimum area concentration for computing velocity", "0.001"
-   "dyn_mass_min", "minimum mass for computing velocity", "0.01 kg/m\ :math:`^2`"
+   "dyn_area_min", "minimum area concentration for computing velocity", "1.e-11"
+   "dyn_mass_min", "minimum mass for computing velocity", "1.e-10 kg/m\ :math:`^2`"
    "**E**", "", ""                       
    "e11, e12, e22", "strain rate tensor components", ""
    "earea", "area of E-cell", "m\ :math:`^2`"
@@ -390,6 +390,8 @@ section :ref:`tabnamelist`.
    "istep0", "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", "zap residual ice below a minimum area", "dyn_area_min"
+   "itd_mass_min", "zap residual ice below a minimum mass", "dyn_mass_min"
    "**J**", "", ""
    "**K**", "", ""
    "kalg", "absorption coefficient for algae", ""

doc/source/science_guide/sg_dynamics.rst

@@ -96,9 +96,12 @@ Note that the VP solver has not yet been tested on the ``tx1`` grid.
 
 The EVP, rEVP, EAP and VP approaches are all available with the B grid. However, at the moment, only the EVP and rEVP schemes are possible with the C grid.
 
-The dynamics are solved for all gridcells with area concentration greater than ``dyn_area_min`` and mass
-greater than ``dyn_mass_min``.  These parameters are respectively 0.001 and 0.01 by default but can be set in 
-namelist.  Lower values can improve the solution but also lead to instabilities.
+The dynamics are solved for all grid cells with area concentration greater than ``dyn_area_min``
+and mass greater than ``dyn_mass_min``.  These parameters can be set in namelist.  Lower
+values can improve the solution with increased computational expense due to additional
+calculations in grid cells with small amounts of ice, but can also lead to instabilities.
+For this reason, default values in the code and base namelist file are set for the B-grid,
+with different values provided for C-grid tests.
 
 Here we summarize the equations and
 direct the reader to the above references for details.

doc/source/science_guide/sg_fundvars.rst

@@ -13,7 +13,6 @@ modeling is to describe the evolution of the ice thickness distribution
 (ITD) in time and space.
 
 In addition to an ice thickness distribution, CICE includes an optional capability for a floe size distribution.
-
 Ice floe horizontal size may change through vertical and lateral growth and melting of existing floes, freezing of new ice, wave breaking, and welding of floes in freezing conditions.  The floe size distribution (FSD) is a probability function that characterizes this variability. The scheme is based on the theoretical framework described in :cite:`Horvat15` for a joint floe size and thickness distribution (FSTD), and was implemented by :cite:`Roach18` and :cite:`Roach19`.  The joint floe size distribution is carried as an area-weighted tracer, defined as the fraction of ice belonging to a given thickness category with lateral floe size belong to a given floe size class. This development includes interactions between sea ice and ocean surface waves. Input data on ocean surface wave spectra at a single time is provided for testing, but as with the other CICE datasets, it should not be used for production runs or publications. It is not recommended to use the FSD without ocean surface waves.
 
 Additional information about the ITD and joint FSTD for CICE can be found in the
@@ -113,3 +112,15 @@ the beginning of the timestep. Rather than recompute the albedo and
 shortwave components at the beginning of the next timestep using new
 values of the downwelling shortwave forcing, the shortwave components
 computed at the end of the last timestep are scaled for the new forcing.
+
+In Icepack, residual amounts of ice may be conservatively removed 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 column physics and
+CICE's dynamic calculations by avoiding residual ice not handled in either
+place. However, ``dyn_area_min`` and ``dyn_mass_min`` should be relatively
+small to avoid removing too much ice. The default behavior sets the column
+physics (itd) parameters to the dynamics values (from namelist).
+``itd_area_min`` and ``itd_mass_min`` can be added to the namelist file
+**ice_in** and set to different values, if desired. Setting them to zero
+turns off residual zapping completely.

doc/source/user_guide/ug_case_settings.rst

@@ -466,6 +466,8 @@ thermo_nml
    "``ktherm``", "``-1``", "thermodynamic model disabled", "1"
    "", "``1``", "Bitz and Lipscomb thermodynamic model", ""
    "", "``2``", "mushy-layer thermodynamic model", ""
+   "``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"
    "``phi_c_slow_mode``", ":math:`0<\phi_c < 1`", "critical liquid fraction", "0.05"
    "``phi_i_mushy``", ":math:`0<\phi_i < 1`", "solid fraction at lower boundary", "0.85"
    "``Rac_rapid_mode``", "real", "critical Rayleigh number", "10.0"
@@ -512,8 +514,8 @@ dynamics_nml
    "``deltaminVP``", "real", "minimum delta for viscosities", "2e-9"
    "``dim_fgmres``", "integer", "maximum number of Arnoldi iterations for FGMRES solver", "50"
    "``dim_pgmres``", "integer", "maximum number of Arnoldi iterations for PGMRES preconditioner", "5"
-   "``dyn_area_min``", "real", "min ice area concentration to activate dynamics", "0.001"
-   "``dyn_mass_min``", "real", "min ice mass to activate dynamics (kg/m\ :math:`^2`)", "0.01"
+   "``dyn_area_min``", "real", "min ice area concentration to activate dynamics", "1.e-11"
+   "``dyn_mass_min``", "real", "min ice mass to activate dynamics (kg/m\ :math:`^2`)", "1.e-10"
    "``e_plasticpot``", "real", "aspect ratio of elliptical plastic potential", "2.0"
    "``e_yieldcurve``", "real", "aspect ratio of elliptical yield curve", "2.0"
    "``elasticDamp``", "real", "elastic damping parameter", "0.36"