diff --git a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 index 84edea237..660423910 100644 --- a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 @@ -51,6 +51,8 @@ module ice_dyn_shared ! shared_mem_1d = 1d without mpi call and refactorization to 1d real (kind=dbl_kind), public :: & + dyn_area_min,& ! minimum ice area concentration to activate dynamics + dyn_mass_min,& ! minimum ice mass to activate dynamics (kg/m^2) elasticDamp ! coefficient for calculating the parameter E, elastic damping parameter ! other EVP parameters @@ -64,9 +66,7 @@ module ice_dyn_shared real (kind=dbl_kind), parameter, public :: & u0 = 5e-5_dbl_kind, & ! residual velocity for seabed stress (m/s) cosw = c1 , & ! cos(ocean turning angle) ! turning angle = 0 - sinw = c0 , & ! sin(ocean turning angle) ! turning angle = 0 - a_min = p001 , & ! minimum ice area - m_min = p01 ! minimum ice mass (kg/m^2) + sinw = c0 ! sin(ocean turning angle) ! turning angle = 0 real (kind=dbl_kind), public :: & revp , & ! 0 for classic EVP, 1 for revised EVP @@ -517,8 +517,8 @@ subroutine dyn_prep1 (nx_block, ny_block, & !----------------------------------------------------------------- ! ice extent mask (T-cells) !----------------------------------------------------------------- - tmphm(i,j) = Tmask(i,j) .and. (aice (i,j) > a_min) & - .and. (Tmass(i,j) > m_min) + tmphm(i,j) = Tmask(i,j) .and. (aice (i,j) > dyn_area_min) & + .and. (Tmass(i,j) > dyn_mass_min) !----------------------------------------------------------------- ! augmented mask (land + open ocean) @@ -721,8 +721,8 @@ subroutine dyn_prep2 (nx_block, ny_block, & do i = ilo, ihi iceXmask_old(i,j) = iceXmask(i,j) ! save ! ice extent mask (U-cells) - iceXmask(i,j) = (Xmask(i,j)) .and. (aiX (i,j) > a_min) & - .and. (Xmass(i,j) > m_min) + iceXmask(i,j) = (Xmask(i,j)) .and. (aiX (i,j) > dyn_area_min) & + .and. (Xmass(i,j) > dyn_mass_min) if (iceXmask(i,j)) then icellX = icellX + 1 diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90 index bc6d224a0..b0ed21553 100644 --- a/cicecore/cicedyn/general/ice_init.F90 +++ b/cicecore/cicedyn/general/ice_init.F90 @@ -14,8 +14,8 @@ module ice_init use ice_kinds_mod use ice_communicate, only: my_task, master_task, ice_barrier - use ice_constants, only: c0, c1, c2, c3, c5, c12, p01, p2, p3, p5, p75, p166, & - cm_to_m + use ice_constants, only: c0, c1, c2, c3, c5, c12, & + p001, p01, p2, p3, p5, p75, p166, cm_to_m use ice_exit, only: abort_ice use ice_fileunits, only: nu_nml, nu_diag, nml_filename, diag_type, & ice_stdout, get_fileunit, release_fileunit, bfbflag, flush_fileunit, & @@ -123,7 +123,7 @@ subroutine input_data e_yieldcurve, e_plasticpot, coriolis, & ssh_stress, kridge, brlx, arlx, & deltaminEVP, deltaminVP, capping, & - elasticDamp + elasticDamp, dyn_area_min, dyn_mass_min use ice_dyn_vp, only: & maxits_nonlin, precond, dim_fgmres, dim_pgmres, maxits_fgmres, & maxits_pgmres, monitor_nonlin, monitor_fgmres, & @@ -255,7 +255,8 @@ subroutine input_data ortho_type, seabed_stress, seabed_stress_method, & k1, k2, alphab, threshold_hw, & deltaminEVP, deltaminVP, capping_method, & - Cf, Pstar, Cstar, Ktens + Cf, Pstar, Cstar, Ktens, & + dyn_area_min, dyn_mass_min namelist /shortwave_nml/ & shortwave, albedo_type, snw_ssp_table, & @@ -413,6 +414,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) 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) @@ -1016,6 +1019,8 @@ subroutine input_data call broadcast_scalar(kstrength, master_task) call broadcast_scalar(Pstar, master_task) call broadcast_scalar(Cstar, master_task) + call broadcast_scalar(dyn_area_min, master_task) + call broadcast_scalar(dyn_mass_min, master_task) call broadcast_scalar(krdg_partic, master_task) call broadcast_scalar(krdg_redist, master_task) call broadcast_scalar(mu_rdg, master_task) @@ -2033,6 +2038,8 @@ subroutine input_data tmpstr2 = ' : unknown value' endif write(nu_diag,1020) ' kdyn = ', kdyn,trim(tmpstr2) + write(nu_diag,1003) ' dyn_area_min = ', dyn_area_min,' : min ice area concentration to activate dynamics' + write(nu_diag,1003) ' dyn_mass_min = ', dyn_mass_min,' : min ice mass to activate dynamics (kg/m2)' if (kdyn >= 1) then if (kdyn == 1 .or. kdyn == 2) then if (revised_evp) then diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in index ed9adc1c5..2e5db8493 100644 --- a/configuration/scripts/ice_in +++ b/configuration/scripts/ice_in @@ -198,6 +198,8 @@ reltol_pgmres = 1e-6 algo_nonlin = 'picard' use_mean_vrel = .true. + dyn_area_min = 0.001d0 + dyn_mass_min = 0.01d0 / &shortwave_nml diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index a8140af3a..645426f33 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -27,7 +27,6 @@ section :ref:`tabnamelist`. "a4Db", "history field accumulations, 4D categories, vertical bio grid", "" "a4Ds", "history field accumulations, 4D categories, vertical snow", "" "a4Df", "history field accumulations, 4D categories, fsd", "" - "a_min", "minimum area concentration for computing velocity", "0.001" "a_rapid_mode", "brine channel diameter", "" "add_mpi_barriers", "turns on MPI barriers for communication throttling", "" "advection", "type of advection algorithm used (‘remap’ or ‘upwind’)", "remap" @@ -204,6 +203,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`" "**E**", "", "" "e11, e12, e22", "strain rate tensor components", "" "earea", "area of E-cell", "m\ :math:`^2`" @@ -419,7 +420,6 @@ section :ref:`tabnamelist`. "ltripole_grid", "flag to signal use of tripole grid", "" "Lvap", "latent heat of vaporization for fresh water", "2.501\ :math:`\times` 10\ :math:`^6` J/kg" "**M**", "", "" - "m_min", "minimum mass for computing velocity", "0.01 kg/m\ :math:`^2`" "m_to_cm", "meters to cm conversion", "100." "m1", "constant for lateral melt rate", "1.6\ :math:`\times`\ 10\ :math:`^{-6}` m/s deg\ :math:`^{-m2}`" "m2", "constant for lateral melt rate", "1.36" diff --git a/doc/source/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index 978da7fcb..ce9c3c4ee 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/doc/source/science_guide/sg_dynamics.rst @@ -96,6 +96,10 @@ 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. + Here we summarize the equations and direct the reader to the above references for details. diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 6b528f7cb..2a11d2325 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -497,6 +497,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" "``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"