Fix for rare instability in (probabilistic) seabed stress (#810)

* Modified doc for Dupont et al ref

* Minor modif to v_i calculation in seabed_stress_factor_prob to prevent rare instability

cicecore/cicedyn/dynamics/ice_dyn_shared.F90

@@ -59,7 +59,7 @@ module ice_dyn_shared
          yield_curve      , & ! 'ellipse' ('teardrop' needs further testing)
          visc_method      , & ! method for viscosity calc at U points (C, CD grids)
          seabed_stress_method ! method for seabed stress calculation
-                              ! LKD: Lemieux et al. 2015, probabilistic: Dupont et al. in prep.
+                              ! LKD: Lemieux et al. 2015, probabilistic: Dupont et al. 2022
 
       real (kind=dbl_kind), parameter, public :: &
          u0    = 5e-5_dbl_kind, & ! residual velocity for seabed stress (m/s)
@@ -1347,8 +1347,9 @@ end subroutine seabed_stress_factor_LKD
 ! a normal distribution with sigma_b = 2.5d0. An improvement would
 ! be to provide the distribution based on high resolution data.
 !
-! Dupont, F. Dumont, D., Lemieux, J.F., Dumas-Lefebvre, E., Caya, A.
-! in prep.
+! Dupont, F., D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, A. Caya (2022).
+! A probabilistic seabed-ice keel interaction model, The Cryosphere, 16,
+! 1963-1977.
 !
 ! authors: D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, F. Dupont
 !
@@ -1481,13 +1482,13 @@ subroutine seabed_stress_factor_prob (nx_block, ny_block,          &
             do n =1, ncat
                v_i = v_i + vcat(n)**2 / (max(acat(n), puny))
             enddo
-            v_i = v_i - m_i**2
+            v_i = max((v_i - m_i**2), puny)
 
             mu_i    = log(m_i/sqrt(c1 + v_i/m_i**2)) ! parameters for the log-normal
             sigma_i = sqrt(log(c1 + v_i/m_i**2))
 
             ! max thickness associated with percentile of log-normal PDF
-            ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al.)
+            ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al. 2022)
 
             x_kmax = exp(mu_i + sqrt(c2*sigma_i)*1.9430d0)
 

cicecore/cicedyn/general/ice_init.F90

@@ -392,7 +392,7 @@ subroutine input_data
       dyscale = 1.0_dbl_kind   ! user defined rectgrid y-grid scale factor (e.g., 1.02)
       close_boundaries = .false.   ! true = set land on edges of grid
       seabed_stress= .false.  ! if true, seabed stress for landfast is on
-      seabed_stress_method  = 'LKD'! LKD = Lemieux et al 2015, probabilistic = Dupont et al. in prep
+      seabed_stress_method  = 'LKD'! LKD = Lemieux et al 2015, probabilistic = Dupont et al. 2022
       k1 = 7.5_dbl_kind       ! 1st free parameter for landfast parameterization
       k2 = 15.0_dbl_kind      ! 2nd free parameter (N/m^3) for landfast parametrization
       alphab = 20.0_dbl_kind  ! alphab=Cb factor in Lemieux et al 2015

doc/source/master_list.bib

@@ -1082,6 +1082,16 @@ @article{Bouchat22
 year = {2022}
 }
 
+@Article{Dupont22,
+  author = {F. Dupont and D. Dumont and J.F. Lemieux and E. Dumas-Lefebvre and A. Caya},
+  title = "{A probabilistic seabed-ice keel interaction model}",
+  journal = TC,
+  year = {2022},
+  volume = {16},
+  pages = {1963-1977},
+  url = {https://doi.org/10.5194/tc-16-1963-2022}
+}
+
 @Article{Tsujino18,
   author = "H. Tsujino and S. Urakawa and R.J. Small and W.M. Kim and S.G. Yeager and et al.",
   title = "{JRA‐55 based surface dataset for driving ocean–sea‐ice models (JRA55‐do)}",

doc/source/science_guide/sg_dynamics.rst

@@ -382,7 +382,7 @@ The value of :math:`k_1` can be changed at runtime using the namelist variable `
 
 This more sophisticated grounding parameterization computes the seabed stress based
 on the probability of contact between the ice thickness distribution
-(ITD) and the seabed. Multi-thickness category models such as CICE typically use a
+(ITD) and the seabed :cite:`Dupont22`. Multi-thickness category models such as CICE typically use a
 few thickness categories (5-10). This crude representation of the ITD
 does not resolve the tail of the ITD, which is crucial for grounding
 events.

doc/source/user_guide/ug_case_settings.rst

@@ -485,7 +485,7 @@ dynamics_nml
    "``revised_evp``", "logical", "use revised EVP formulation", "``.false.``"
    "``seabed_stress``", "logical", "use seabed stress parameterization for landfast ice", "``.false.``"
    "``seabed_stress_method``", "``LKD``", "linear keel draft method :cite:`Lemieux16`", "``LKD``"
-   "", "``probabilistic``", "probability of contact method (Dupont et al., in prep)", ""
+   "", "``probabilistic``", "probability of contact method :cite:`Dupont22`", ""
    "``ssh_stress``", "``coupled``", "computed from coupled sea surface height gradient", "``geostrophic``"
    "", "``geostropic``", "computed from ocean velocity", ""
    "``threshold_hw``", "real", "Max water depth for grounding (see :cite:`Amundrud04`)", "30."