Merge pull request #309 from mom-ocean/safer-min_thickness

Safer min thickness

src/mom5/ocean_core/ocean_topog.F90

@@ -32,7 +32,15 @@ module ocean_topog_mod
 !  <DATA NAME="min_thickness" TYPE="real">
 !  min_thickness is only used for Mosaic grid. Since there is no kmt available
 !  in mosaic grid, need to set min_thickness to configure kmt based on ht and zw.
-!  Default min_thickness=1.0e-3 metre.
+!  min_thickness should be set in input namelist (default value -1.0 will cause termination).
+!  min_thickness should be much smaller than minimum vertical grid spacing,
+!  but for floating point representation reasons min_thickness should be no less than 
+!  spacing(max_depth,kind=real32) where max_depth is the maximum depth of the model 
+!  and the real32 parameter can be obtained from iso_fortran_env.
+!  min_thickness=1.0e-3 metre is usually a good choice.
+!  The previous default was 1.0 m; this could be used for consistency when
+!  restarting from runs which did not specify a value,
+!  but may cause problems: https://github.com/COSIMA/access-om2/issues/161
 !  </DATA> 
 !  <DATA NAME="kmt_recompute" TYPE="logical">
 !  To recompute the kmt array based on min_thickness.  This step is not recommended
@@ -72,7 +80,7 @@ module ocean_topog_mod
 logical :: flat_bottom          = .false.
 integer :: flat_bottom_kmt      = 50
 real    :: flat_bottom_ht       = 5500.0
-real    :: min_thickness        = 1.0e-3
+real    :: min_thickness        = -1.0
 integer :: kmt_recompute_offset = 0
 logical :: kmt_recompute        = .false.
 logical :: write_topog          = .false.
@@ -190,6 +198,16 @@ subroutine ocean_topog_init (Domain, Grid, grid_file, vert_coordinate_type)
      call read_data(grd_file, "kmt", Grid%kmt(isc:iec,jsc:jec), Domain%domain2d)
 
   elseif(file_exist(ocean_topog)) then
+
+     if(min_thickness < 0.0) then
+        call mpp_error(FATAL, &
+           'ocean_topog_mod: min_thickness must be explicitly specified in ocean_topog_nml. '// &
+           'min_thickness=1.0e-3 metre is usually a good choice. '// &
+           'The previous default was 1.0 m; this could be used for consistency when '// &
+           'restarting from runs which did not specify a value, '// &
+           'but may cause problems: https://github.com/COSIMA/access-om2/issues/161')
+     endif
+
      call read_data(ocean_topog, 'depth', Grid%ht(isc:iec,jsc:jec), Domain%domain2d)
 
      !--- calculate kmt based on ht and zw.
@@ -214,6 +232,16 @@ subroutine ocean_topog_init (Domain, Grid, grid_file, vert_coordinate_type)
 
   if(kmt_recompute) then 
       call mpp_error(NOTE, 'ocean_topog_mod: recomputing kmt array given ht, zw, and min_thickness.')
+
+      if(min_thickness < 0.0) then
+        call mpp_error(FATAL, &
+           'ocean_topog_mod: min_thickness must be explicitly specified in ocean_topog_nml. '// &
+           'min_thickness=1.0e-3 metre is usually a good choice. '// &
+           'The previous default was 1.0 m; this could be used for consistency when '// &
+           'restarting from runs which did not specify a value, '// &
+           'but may cause problems: https://github.com/COSIMA/access-om2/issues/161')
+      endif
+
       do j=jsc,jec
          do i=isc,iec
             if(Grid%kmt(i,j) > 1) then