Icepack area underflow patch (#184)

New parameter RIDGE_AREA_UNDERFLOW is introduced in order to avoid the occurrence of undetectably small fractional area and negative ice volume following ridging adjustments from Icepack. The default is 0.0, which does not introduce changes in existing runs. A reasonable value would be 10^-26 which would correspond to an Angstrom-scale ice patch within a km-scale grid cell. This addresses #178 and #183

src/ice_ridge.F90

@@ -11,7 +11,9 @@ module ice_ridging_mod
 ! 1) implement new snow_to_ocean diagnostic to record this flux.               !
 ! 2) implement ridging_rate diagnostics: ridging_shear, ridging_conv           !
 ! 3) implement "do_j" style optimization as in "compress_ice" or               !
-!    "adjust_ice_categories" (SIS_transport.F90) if deemed necessary           !
+!    "adjust_ice_categories" (SIS_transport.F90) if deemed necessary
+!
+!
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
 
 use SIS_diag_mediator, only : post_SIS_data, query_SIS_averaging_enabled, SIS_diag_ctrl
@@ -43,9 +45,13 @@ module ice_ridging_mod
 
 type, public :: ice_ridging_CS ; private
   logical :: &
-       new_rdg_partic = .false., & !< .true. = new participation, .false. = Thorndike et al 75
-       new_rdg_redist = .false.    !< .true. = new redistribution, .false. = Hibler 80
+  new_rdg_partic = .false., & !< .true. = new participation, .false. = Thorndike et al 75
+  new_rdg_redist = .false.    !< .true. = new redistribution, .false. = Hibler 80
   real :: mu_rdg = 3.0 !< e-folding scale of ridged ice, new_rdg_partic (m^0.5)
+  real :: area_underflow = 0.0 !< a non-dimesional fractional area underflow limit for the sea-ice
+                               !! ridging scheme. This is defaulted to zero, but a reasonable
+                               !! value might be 10^-26 which for a km square grid cell
+                               !! would equate to an Angstrom scale ice patch.
 end type ice_ridging_CS
 
 contains
@@ -72,6 +78,10 @@ subroutine ice_ridging_init(G, IG, PF, CS, US)
     call get_param(PF, mdl, "RIDGE_MU", CS%mu_rdg, &
                    "E-folding scale of ridge ice from Lipscomb et al. 2007", &
                    units="m^0.5", default=3.0)
+    call get_param(PF, mdl, "RIDGE_AREA_UNDERFLOW", CS%area_underflow, &
+                   "A fractional area limit below which ice fraction is set to zero "//&
+                   "A reasonable default value for a km scale grid cell is 10^-24.",&
+                   units="none", default=0.0)
   endif
 
   ncat=IG%CatIce ! The number of sea-ice thickness categories
@@ -239,7 +249,11 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, CS, US, d
 
   call icepack_query_tracer_sizes(ncat_out=ncat_out,ntrcr_out=ntrcr_out, nilyr_out=nilyr_out, nslyr_out=nslyr_out)
 
-  if (nIlyr .ne. nilyr_out .or. nSlyr .ne. nslyr_out ) call SIS_error(FATAL,'nilyr or nslyr mismatch with Icepack')
+  if (nIlyr .ne. nilyr_out .or. nSlyr .ne. nslyr_out ) &
+    call SIS_error(FATAL,"Oops!! It looks like you are trying to use sea-ice ridging "//&
+                         "but did not include the Icepack (https://github.com/CICE-Consortium/Icepack)"//&
+                         "source code repository in your compilation procedure, and are instead using the default "//&
+                         "stub routine contained in config_src/external. Adjust your compilation accordingly." )
 
   ! copy strain calculation code from SIS_C_dynamics; might be a more elegant way ...
   !
@@ -475,7 +489,10 @@ subroutine ice_ridging(IST, G, IG, mca_ice, mca_snow, mca_pond, TrReg, CS, US, d
 
       ! ! output: snow/ice masses/thicknesses
       do k=1,nCat
-
+        if (aicen(k) < CS%area_underflow) then
+           aicen(k)=0.0
+           vicen(k)=0.0
+        endif
         if (aicen(k) > 0.0) then
           IST%part_size(i,j,k)  = aicen(k)
           mca_ice(i,j,k)  = vicen(k)*Rho_ice * US%m_to_Z