Merge pull request #1457 from herrwang0/fix-topoinit_2

Bugfix: a couple of mask_depth related issues

src/initialization/MOM_grid_initialize.F90

@@ -1187,7 +1187,7 @@ end function Adcroft_reciprocal
 !> Initializes the grid masks and any metrics that come with masks already applied.
 !!
 !!    Initialize_masks sets mask2dT, mask2dCu, mask2dCv, and mask2dBu to mask out
-!! flow over any points which are shallower than Dmin and permit an
+!! flow over any points which are shallower than Dmask and permit an
 !! appropriate treatment of the boundary conditions.  mask2dCu and mask2dCv
 !! are 0.0 at any points adjacent to a land point.  mask2dBu is 0.0 at
 !! any land or boundary point.  For points in the interior, mask2dCu,
@@ -1199,7 +1199,7 @@ subroutine initialize_masks(G, PF, US)
   ! Local variables
   real :: m_to_Z_scale ! A unit conversion factor from m to Z.
   real :: m_to_L  ! A unit conversion factor [L m-1 ~> nondim]
-  real :: Dmin       ! The depth for masking in the same units as G%bathyT [Z ~> m].
+  real :: Dmask      ! The depth for masking in the same units as G%bathyT [Z ~> m].
   real :: min_depth  ! The minimum ocean depth in the same units as G%bathyT [Z ~> m].
   real :: mask_depth ! The depth shallower than which to mask a point as land [Z ~> m].
   character(len=40)  :: mdl = "MOM_grid_init initialize_masks"
@@ -1226,39 +1226,39 @@ subroutine initialize_masks(G, PF, US)
                   'MASKING_DEPTH is larger than MINIMUM_DEPTH and therefore ignored.')
   endif
 
-  Dmin = min_depth
-  if (mask_depth /= -9999.*m_to_Z_scale) Dmin = mask_depth
+  Dmask = mask_depth
+  if (mask_depth == -9999.*m_to_Z_scale) Dmask = min_depth
 
   G%mask2dCu(:,:) = 0.0 ; G%mask2dCv(:,:) = 0.0 ; G%mask2dBu(:,:) = 0.0
 
   ! Construct the h-point or T-point mask
   do j=G%jsd,G%jed ; do i=G%isd,G%ied
-    if (G%bathyT(i,j) <= Dmin) then
+    if (G%bathyT(i,j) <= Dmask) then
       G%mask2dT(i,j) = 0.0
     else
       G%mask2dT(i,j) = 1.0
     endif
   enddo ; enddo
 
   do j=G%jsd,G%jed ; do I=G%isd,G%ied-1
-    if ((G%bathyT(i,j) <= Dmin) .or. (G%bathyT(i+1,j) <= Dmin)) then
+    if ((G%bathyT(i,j) <= Dmask) .or. (G%bathyT(i+1,j) <= Dmask)) then
       G%mask2dCu(I,j) = 0.0
     else
       G%mask2dCu(I,j) = 1.0
     endif
   enddo ; enddo
 
   do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied
-    if ((G%bathyT(i,j) <= Dmin) .or. (G%bathyT(i,j+1) <= Dmin)) then
+    if ((G%bathyT(i,j) <= Dmask) .or. (G%bathyT(i,j+1) <= Dmask)) then
       G%mask2dCv(i,J) = 0.0
     else
       G%mask2dCv(i,J) = 1.0
     endif
   enddo ; enddo
 
   do J=G%jsd,G%jed-1 ; do I=G%isd,G%ied-1
-    if ((G%bathyT(i+1,j) <= Dmin) .or. (G%bathyT(i+1,j+1) <= Dmin) .or. &
-        (G%bathyT(i,j) <= Dmin) .or. (G%bathyT(i,j+1) <= Dmin)) then
+    if ((G%bathyT(i+1,j) <= Dmask) .or. (G%bathyT(i+1,j+1) <= Dmask) .or. &
+        (G%bathyT(i,j) <= Dmask) .or. (G%bathyT(i,j+1) <= Dmask)) then
       G%mask2dBu(I,J) = 0.0
     else
       G%mask2dBu(I,J) = 1.0

src/initialization/MOM_shared_initialization.F90

@@ -195,8 +195,8 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US)
   character(len=200) :: topo_edits_file, inputdir ! Strings for file/path
   character(len=40)  :: mdl = "apply_topography_edits_from_file" ! This subroutine's name.
   integer :: i, j, n, ncid, n_edits, i_file, j_file, ndims, sizes(8)
-  logical :: found
   logical :: topo_edits_change_mask
+  real :: min_depth, mask_depth
 
   call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
 
@@ -210,6 +210,17 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US)
   call get_param(param_file, mdl, "ALLOW_LANDMASK_CHANGES", topo_edits_change_mask, &
                  "If true, allow topography overrides to change land mask.", &
                  default=.false.)
+  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
+                 "If MASKING_DEPTH is unspecified, then anything shallower than "//&
+                 "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
+                 "If MASKING_DEPTH is specified, then all depths shallower than "//&
+                 "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
+                 units="m", default=0.0, scale=m_to_Z)
+  call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, &
+                 "The depth below which to mask points as land points, for which all "//&
+                 "fluxes are zeroed out. MASKING_DEPTH needs to be smaller than MINIMUM_DEPTH", &
+                 units="m", default=-9999.0, scale=m_to_Z)
+  if (mask_depth == -9999.*m_to_Z) mask_depth = min_depth
 
   if (len_trim(topo_edits_file)==0) return
 
@@ -249,7 +260,7 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US)
     i = ig(n) - G%isd_global + 2 ! +1 for python indexing and +1 for ig-isd_global+1
     j = jg(n) - G%jsd_global + 2
     if (i>=G%isc .and. i<=G%iec .and. j>=G%jsc .and. j<=G%jec) then
-      if (new_depth(n)/=0.) then
+      if (new_depth(n)*m_to_Z /= mask_depth) then
         write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') &
           'Ocean topography edit: ', n, ig(n), jg(n), D(i,j)/m_to_Z, '->', abs(new_depth(n)), i, j
         D(i,j) = abs(m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
@@ -434,6 +445,10 @@ subroutine limit_topography(D, G, param_file, max_depth, US)
     do j=G%jsd,G%jed ; do i=G%isd,G%ied
       if (D(i,j) > mask_depth) then
         D(i,j) = min( max( D(i,j), min_depth ), max_depth )
+      else
+        ! This statement is required for cases with masked-out PEs over the land,
+        ! to remove the large initialized values (-9e30) from the halos.
+        D(i,j) = mask_depth
       endif
     enddo ; enddo
   endif