Bugfix: Minor changes on clipping topography from file (#1428)

* Minor changes to topography initialized from files
This is intended to extend the compatibility with negative topography.
1. A bug that occurs when clipping topography with a negative
MINIMUM_DEPTH is fixed.
2. MASKING_DEPTH can now be negative.
3. A warning will be given if MASKING_DEPTH is set to be smaller than
MINIMUM_DEPTH.

* Adding an exception to avoid answer changes
To keep answer unchnaged in the test cases (at land points),
D is clipped at 0.5*min_depth when min_depth > 0.

* Change the if-statement on whether MASK_DEPTH is default to a simple comparison

* Change the if-statement on whether MASK_DEPTH is default

* Change the if-statement comparison on MASK_DEPTH
Reverted back to a simple comparison. A mistake in a previous commit in
MOM_grid_initialize.F90 is corrected.

Co-authored-by: Robert Hallberg <[email protected]>

src/initialization/MOM_grid_initialize.F90

@@ -9,7 +9,7 @@ module MOM_grid_initialize
 use MOM_domains,       only : To_North, To_South, To_East, To_West
 use MOM_domains,       only : MOM_domain_type, clone_MOM_domain, deallocate_MOM_domain
 use MOM_dyn_horgrid,   only : dyn_horgrid_type, set_derived_dyn_horgrid
-use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe
+use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
 use MOM_error_handler, only : callTree_enter, callTree_leave
 use MOM_file_parser,   only : get_param, log_param, log_version, param_file_type
 use MOM_io,            only : MOM_read_data, slasher, file_exists, stdout
@@ -1217,11 +1217,17 @@ subroutine initialize_masks(G, PF, US)
                  units="m", default=0.0, scale=m_to_Z_scale)
   call get_param(PF, mdl, "MASKING_DEPTH", mask_depth, &
                  "The depth below which to mask points as land points, for which all "//&
-                 "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", &
+                 "fluxes are zeroed out. MASKING_DEPTH needs to be smaller than MINIMUM_DEPTH", &
                  units="m", default=-9999.0, scale=m_to_Z_scale)
 
+  if (mask_depth > min_depth) then
+    mask_depth = -9999.0*m_to_Z_scale
+    call MOM_error(WARNING, "MOM_grid_init: initialize_masks "//&
+                  'MASKING_DEPTH is larger than MINIMUM_DEPTH and therefore ignored.')
+  endif
+
   Dmin = min_depth
-  if (mask_depth>=0.) Dmin = mask_depth
+  if (mask_depth /= -9999.*m_to_Z_scale) Dmin = mask_depth
 
   G%mask2dCu(:,:) = 0.0 ; G%mask2dCv(:,:) = 0.0 ; G%mask2dBu(:,:) = 0.0
 

src/initialization/MOM_shared_initialization.F90

@@ -413,17 +413,27 @@ subroutine limit_topography(D, G, param_file, max_depth, US)
                  "The depth below which to mask the ocean as land.", &
                  units="m", default=-9999.0, scale=m_to_Z, do_not_log=.true.)
 
-! Make sure that min_depth < D(x,y) < max_depth
-  if (mask_depth < -9990.*m_to_Z) then
-    do j=G%jsd,G%jed ; do i=G%isd,G%ied
-      D(i,j) = min( max( D(i,j), 0.5*min_depth ), max_depth )
-    enddo ; enddo
+  if (mask_depth > min_depth) then
+    mask_depth = -9999.0*m_to_Z
+    call MOM_error(WARNING, "MOM_shared_initialization: limit_topography "//&
+                  'MASKING_DEPTH is larger than MINIMUM_DEPTH and therefore ignored.')
+  endif
+
+  ! Make sure that min_depth < D(x,y) < max_depth for ocean points
+  if (mask_depth == -9999.*m_to_Z) then
+    if (min_depth > 0.0) then  ! This is retained to avoid answer changes (over the land points) in the test cases.
+      do j=G%jsd,G%jed ; do i=G%isd,G%ied
+        D(i,j) = min( max( D(i,j), 0.5*min_depth ), max_depth )
+      enddo ; enddo
+    else
+      do j=G%jsd,G%jed ; do i=G%isd,G%ied
+        D(i,j) = min( max( D(i,j), min_depth ), max_depth )
+      enddo ; enddo
+    endif
   else
     do j=G%jsd,G%jed ; do i=G%isd,G%ied
-      if (D(i,j)>0.) then
+      if (D(i,j) > mask_depth) then
         D(i,j) = min( max( D(i,j), min_depth ), max_depth )
-      else
-        D(i,j) = 0.
       endif
     enddo ; enddo
   endif