Fix for sigma_i going to zero in ITD landfast ice. (#189)

- We can't divide by zero, so don't let sigma_i be zero

- Version coming from the CICE Consortium

- Update Dupont citation

docs/zotero.bib

@@ -1,9 +1,10 @@
 @article{dupont2021,
 	author = "F. Dupont and D. Dumont and J.-F. Lemieux and E. Dumas-Lefebvre and A. Caya",
-	year = 2021,
+	year = 2022,
 	title = "A probabilistic seabed-ice keel interaction model",
 	journal = "The Cryosphere",
-	note = "in review"
+	volume = "16",
+	pages = "1963--1977"
 }
 
 @article{lemieux2015,

src/SIS_dyn_cgrid.F90

@@ -1853,8 +1853,9 @@ end subroutine basal_stress_coeff_C
 !! 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
 !!
@@ -1902,8 +1903,8 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS)
   real :: rho_water  ! water density [R ~> kg m-3]
   real :: pi         ! [nondim]
   integer :: i, ii, j, isc, iec, jsc, jec, k, n, ncat
-  real :: ci_u ! Concentration at u-points [nondim]
-  real :: ci_v ! Concentration at u-points [nondim]
+  real :: ci_u       ! Concentration at u-points [nondim]
+  real :: ci_v       ! Concentration at u-points [nondim]
 
   isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
   ncat = IG%CatIce
@@ -1956,50 +1957,54 @@ subroutine basal_stress_coeff_itd(G, IG, IST, sea_lev, CS)
 
         ! parameters for the log-normal
         mu_i    = log(m_i/(CS%onemeter * sqrt(1.0 + v_i/m_i**2)))
-        sigma_i = max(sqrt(log(1.0 + v_i/m_i**2)), CS%puny)
+        sigma_i = sqrt(log(1.0 + 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 = CS%onemeter * exp(mu_i + sqrt(2.0*sigma_i)*CS%basal_stress_cutoff)
-
-        ! Set x_kmax to hlev of the last category where there is ice
-        ! when there is no ice in the last category
-        cut = x_k(CS%ncat_i)
-        do n = ncat,-1,1
-          if (acat(n) < CS%puny) then
-            cut = hin_max(n-1)
-          else
-            exit
-          endif
-        enddo
-        x_kmax = min(cut, x_kmax)
-
-        g_k(:) = exp(-(log(x_k(:)/CS%onemeter) - mu_i) ** 2 / (2.0 * sigma_i ** 2)) / &
-                 (x_k(:) * sigma_i * sqrt(2.0 * pi))
-
-        b_n(:)  = exp(-(y_n(:) - mu_b) ** 2 / (2.0 * CS%sigma_b(i,j) ** 2)) / (CS%sigma_b(i,j) * sqrt(2.0*pi))
-
-        P_x(:) = g_k(:) * wid_i
-        P_y(:) = b_n(:) * wid_b
-
-        do n =1, CS%ncat_i
-          if (x_k(n) > x_kmax) P_x(n)=0.0
-        enddo
-
-        ! calculate Tb factor at t-location
-        do n=1, CS%ncat_i
-          gt(:) = (y_n(:) <= rho_ice*x_k(n)/rho_water)
-          tmp(:) = merge(1,0,gt(:))
-          ii = sum(tmp)
-          if (ii == 0) then
-            tb_tmp(n) = 0.0
-          else
-            tb_tmp(n) = max(CS%basal_stress_mu_s * G%g_Earth * P_x(n) * &
-                        sum(P_y(1:ii)*(rho_ice*x_k(n) - rho_water*y_n(1:ii))), 0.0)
-          endif
-        enddo
-        Tbt(i,j) = sum(tb_tmp) * exp(-CS%lemieux_alphab * (1.0 - atot))
+        ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al. 2022)
+
+        if (sigma_i > 0) then
+          x_kmax = CS%onemeter * exp(mu_i + sqrt(2.0*sigma_i)*CS%basal_stress_cutoff)
+
+          ! Set x_kmax to hlev of the last category where there is ice
+          ! when there is no ice in the last category
+          cut = x_k(CS%ncat_i)
+          do n = ncat,-1,1
+            if (acat(n) < CS%puny) then
+              cut = hin_max(n-1)
+            else
+              exit
+            endif
+          enddo
+          x_kmax = min(cut, x_kmax)
+
+          g_k(:) = exp(-(log(x_k(:)/CS%onemeter) - mu_i) ** 2 / (2.0 * sigma_i ** 2)) / &
+                   (x_k(:) * sigma_i * sqrt(2.0 * pi))
+
+          b_n(:)  = exp(-(y_n(:) - mu_b) ** 2 / (2.0 * CS%sigma_b(i,j) ** 2)) / (CS%sigma_b(i,j) * sqrt(2.0*pi))
+
+          P_x(:) = g_k(:) * wid_i
+          P_y(:) = b_n(:) * wid_b
+
+          do n =1, CS%ncat_i
+            if (x_k(n) > x_kmax) P_x(n)=0.0
+          enddo
+
+          ! calculate Tb factor at t-location
+          do n=1, CS%ncat_i
+            gt(:) = (y_n(:) <= rho_ice*x_k(n)/rho_water)
+            tmp(:) = merge(1,0,gt(:))
+            ii = sum(tmp)
+            if (ii == 0) then
+              tb_tmp(n) = 0.0
+            else
+              tb_tmp(n) = max(CS%basal_stress_mu_s * G%g_Earth * P_x(n) * &
+                          sum(P_y(1:ii)*(rho_ice*x_k(n) - rho_water*y_n(1:ii))), 0.0)
+            endif
+          enddo
+          Tbt(i,j) = sum(tb_tmp) * exp(-CS%lemieux_alphab * (1.0 - atot))
+        else
+          Tbt(i,j) = 0.0
+        endif
       endif
     enddo
   enddo

src/SIS_open_boundary.F90

@@ -270,7 +270,7 @@ subroutine open_boundary_config(G, US, param_file, OBC)
         (OBC%freeslip_vorticity .and. OBC%computed_vorticity) .or.  &
         (OBC%freeslip_vorticity .and. OBC%specified_vorticity) .or.  &
         (OBC%computed_vorticity .and. OBC%specified_vorticity))  &
-         call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"// &
+         call MOM_error(FATAL, "SIS_open_boundary.F90, open_boundary_config:\n"// &
          "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"// &
          "and OBC_IMPORTED_VORTICITY can be True at once.")
     call get_param(param_file, mdl, "OBC_ZERO_STRAIN", OBC%zero_strain, &
@@ -364,7 +364,7 @@ subroutine open_boundary_config(G, US, param_file, OBC)
       elseif (segment_str(1:2) == 'J=') then
         call setup_v_point_obc(OBC, G, US, segment_str, l, param_file, reentrant_x)
       else
-        call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config: "//&
+        call MOM_error(FATAL, "SIS_open_boundary.F90, open_boundary_config: "//&
              "Unable to interpret "//segment_param_str//" = "//trim(segment_str))
       endif
     enddo
@@ -412,7 +412,7 @@ subroutine open_boundary_config(G, US, param_file, OBC)
   ! Safety check
   if ((OBC%open_u_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) .and. &
        .not.G%symmetric ) call MOM_error(FATAL, &
-       "MOM_open_boundary, open_boundary_config: "//&
+       "SIS_open_boundary, open_boundary_config: "//&
        "Symmetric memory must be used when using Flather OBCs.")
 
   if (.not.(OBC%specified_u_BCs_exist_globally .or. OBC%specified_v_BCs_exist_globally .or. &
@@ -652,7 +652,7 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y)
       OBC%segment(l_seg)%specified_grad = .true.
       OBC%segment(l_seg)%g_values_needed = .true.
     else
-      call MOM_error(FATAL, "MOM_open_boundary.F90, setup_u_point_obc: "//&
+      call MOM_error(FATAL, "SIS_open_boundary.F90, setup_u_point_obc: "//&
                      "String '"//trim(action_str(a_loop))//"' not understood.")
     endif
     if (OBC%segment(l_seg)%nudged .or. OBC%segment(l_seg)%nudged_tan) then
@@ -774,7 +774,7 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x)
       OBC%segment(l_seg)%specified_grad = .true.
       OBC%segment(l_seg)%g_values_needed = .true.
     else
-      call MOM_error(FATAL, "MOM_open_boundary.F90, setup_v_point_obc: "//&
+      call MOM_error(FATAL, "SIS_open_boundary.F90, setup_v_point_obc: "//&
                      "String '"//trim(action_str(a_loop))//"' not understood.")
     endif
     if (OBC%segment(l_seg)%nudged .or. OBC%segment(l_seg)%nudged_tan) then
@@ -1039,14 +1039,14 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC)
   do j=G%jsd,G%jed ; do i=G%isd,G%ied
     if (color(i,j) /= color2(i,j)) then
       fatal_error = .True.
-      write(mesg,'("MOM_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", &
+      write(mesg,'("SIS_open_boundary: problem with OBC segments specification at ",I5,",",I5," during\n", &
           "the masking of the outside grid points.")') i, j
       call MOM_error(WARNING,"MOM register_tracer: "//mesg, all_print=.true.)
     endif
     if (color(i,j) == cout) G%bathyT(i,j) = min_depth
   enddo ; enddo
   if (fatal_error) call MOM_error(FATAL, &
-      "MOM_open_boundary: inconsistent OBC segments.")
+      "SIS_open_boundary: inconsistent OBC segments.")
 
   deallocate(color)
   deallocate(color2)