diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index aa78c398d..8a139883f 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/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, & diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index 0c9c7ffe5..2d0cffaab 100644 --- a/configuration/scripts/ice_in +++ b/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 diff --git a/configuration/scripts/options/set_nml.gridc b/configuration/scripts/options/set_nml.gridc index a04fab4fd..63c5af2bb 100644 --- a/configuration/scripts/options/set_nml.gridc +++ b/configuration/scripts/options/set_nml.gridc @@ -1,2 +1,4 @@ grid_ice = 'C' +dyn_area_min = 0.001d0 +dyn_mass_min = 0.01d0 diff --git a/configuration/scripts/options/set_nml.gridcd b/configuration/scripts/options/set_nml.gridcd index 7889e64f4..4cf521a4b 100644 --- a/configuration/scripts/options/set_nml.gridcd +++ b/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' diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index f46fbaa8d..479c0a18c 100644 --- a/doc/source/cice_index.rst +++ b/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", "" diff --git a/doc/source/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index ce9c3c4ee..23d1fbc33 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/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. diff --git a/doc/source/science_guide/sg_fundvars.rst b/doc/source/science_guide/sg_fundvars.rst index 2d6f50328..89b095c70 100644 --- a/doc/source/science_guide/sg_fundvars.rst +++ b/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. diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 98bc7268f..3882073de 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/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" diff --git a/icepack b/icepack index 4954a6f90..eedb51924 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit 4954a6f9033f78e5c32bf33780384cbf2d0843e6 +Subproject commit eedb519247f48bce9fa1d1b275b5cc3dc07643e4