New function open_boundary_impose_normal_slope()

- initialize_masks() in MOM_grid_initialization.F90 used to
  read the open-boundary parameters and then adjust the topography
  on the edge of the domain only.
- There is no mention of "Flather" in MOM_grid_initialization.F90
  any more.
- MOM_fixed_initialization.F90 now calls this new routine prior
  to calling initialize_masks().
- I had to move a pass_var(G%bathyT) out of initialize_masks() and
  up to MOM_initialize_fixed() to do this.
- This is part of a larger re-factor of OBC code.
- No answer changes.

src/core/MOM_open_boundary.F90

@@ -21,6 +21,7 @@ module MOM_open_boundary
 public open_boundary_init
 public open_boundary_query
 public open_boundary_end
+public open_boundary_impose_normal_slope
 public Radiation_Open_Bdry_Conds
 public set_Flather_positions
 public set_Flather_data
@@ -132,6 +133,14 @@ subroutine open_boundary_config(G, param_file, OBC)
                  "Apply a Flather open boundary condition on the southern\n"//&
                  "side of the global domain", &
                  default=.false.)
+
+  ! Safety check
+  if ((OBC%apply_OBC_u_flather_west .or. OBC%apply_OBC_v_flather_south) .and. &
+      .not.G%symmetric ) call MOM_error(FATAL, &
+                 "MOM_open_boundary, open_boundary_config: "//&
+                 "Symmetric memory must be used when APPLY_OBC_U_FLATHER_WEST "//&
+                 "or APPLY_OBC_V_FLATHER_SOUTH is true.")
+
   if (.not.(OBC%apply_OBC_u .or. OBC%apply_OBC_v .or. &
             OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south .or. &
             OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west)) then
@@ -199,6 +208,32 @@ subroutine open_boundary_end(OBC)
   deallocate(OBC)
 end subroutine open_boundary_end
 
+!> Sets the slope of bathymetry normal to an open bounndary to zero.
+subroutine open_boundary_impose_normal_slope(OBC, G, depth)
+  type(ocean_OBC_type),              pointer       :: OBC !< Open boundary control structure
+  type(ocean_grid_type),             intent(inout) :: G !< Ocean grid structure
+  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: depth !< Bathymetry at h-points
+  ! Local variables
+  integer :: i, j
+
+  if (.not.associated(OBC)) return
+
+  if (associated(OBC%OBC_kind_u)) then
+    do j=G%jsd,G%jed ; do I=G%isd,G%ied-1
+      if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) depth(i+1,j) = depth(i,j)
+      if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) depth(i,j) = depth(i+1,j)
+    enddo ; enddo
+  endif
+
+  if (associated(OBC%OBC_kind_v)) then
+    do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied
+      if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) depth(i,j+1) = depth(i,j)
+      if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) depth(i,j) = depth(i,j+1)
+    enddo ; enddo
+  endif
+
+end subroutine open_boundary_impose_normal_slope
+
 !> Diagnose radiation conditions at open boundaries
 subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
                                      h_new, h_old, G)

src/initialization/MOM_fixed_initialization.F90

@@ -20,7 +20,7 @@ module MOM_fixed_initialization
 use MOM_grid_initialize, only : initialize_masks, set_grid_metrics
 use MOM_open_boundary, only : ocean_OBC_type
 use MOM_open_boundary, only : open_boundary_config, open_boundary_query
-use MOM_open_boundary, only : set_Flather_positions
+use MOM_open_boundary, only : set_Flather_positions, open_boundary_impose_normal_slope
 use MOM_string_functions, only : uppercase
 use user_initialization, only : user_initialize_topography, USER_set_OBC_positions
 use DOME_initialization, only : DOME_initialize_topography, DOME_set_OBC_positions
@@ -104,8 +104,13 @@ subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir)
     call set_Flather_positions(G, OBC)
   endif
 
-!    This call sets seamasks that prohibit flow over any point with  !
-!  a bottom that is shallower than min_depth from PF.                !
+  ! To initialize masks, the bathymetry in halo regions must be filled in
+  call pass_var(G%bathyT, G%Domain)
+
+  ! Make bathymetry consistent with open boundaries
+  call open_boundary_impose_normal_slope(OBC, G, G%bathyT)
+
+  ! This call sets masks that prohibit flow over any point interpreted as land
   call initialize_masks(G, PF)
   if (debug) then
     call hchksum(G%bathyT, 'MOM_initialize_fixed: depth ', G%HI, haloshift=1)

src/initialization/MOM_grid_initialize.F90

@@ -1343,8 +1343,6 @@ subroutine initialize_masks(G, PF)
 
   real :: Dmin, min_depth, mask_depth
   integer :: i, j
-  logical :: apply_OBC_u_flather_east, apply_OBC_u_flather_west
-  logical :: apply_OBC_v_flather_north, apply_OBC_v_flather_south
   character(len=40)  :: mod = "MOM_grid_init initialize_masks"
 
   call callTree_enter("initialize_masks(), MOM_grid_initialize.F90")
@@ -1358,65 +1356,13 @@ subroutine initialize_masks(G, PF)
                  "The depth below which to mask points as land points, for which all\n"//&
                  "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", &
                  units="m", default=-9999.0)
-  call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", apply_OBC_u_flather_east,&
-                 "Apply a Flather open boundary condition on the eastern \n"//&
-                 "side of the global domain", default=.false.)
-  call get_param(PF, mod, "APPLY_OBC_U_FLATHER_WEST", apply_OBC_u_flather_west,&
-                 "Apply a Flather open boundary condition on the western \n"//&
-                 "side of the global domain", default=.false.)
-  call get_param(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", apply_OBC_v_flather_north,&
-                 "Apply a Flather open boundary condition on the northern \n"//&
-                 "side of the global domain", default=.false.)
-  call get_param(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,&
-                 "Apply a Flather open boundary condition on the southern \n"//&
-                 "side of the global domain", default=.false.)
-
-  if ((apply_OBC_u_flather_west .or. apply_OBC_v_flather_south) .and. &
-      .not.G%symmetric ) &
-    call MOM_error(FATAL, "Symmetric memory must be used when "//&
-      "APPLY_OBC_U_FLATHER_WEST or APPLY_OBC_V_FLATHER_SOUTH is true.")
 
   Dmin = min_depth
   if (mask_depth>=0.) Dmin = mask_depth
 
-  call pass_var(G%bathyT, G%Domain)
   G%mask2dCu(:,:) = 0.0 ; G%mask2dCv(:,:) = 0.0 ; G%mask2dBu(:,:) = 0.0
 
-  ! Extrapolate the bottom depths at any points that are subject to Flather
-  ! open boundary conditions.  This should be generalized for Flather OBCs
-  ! that are not necessarily at the edges of the domain.
-  if (apply_OBC_u_flather_west) then
-    do j=G%jsd,G%jed ; do I=G%isd+1,G%ied
-      if ((I+G%idg_offset) == G%isg) then
-        G%bathyT(i-1,j) = G%bathyT(i,j)
-      endif
-    enddo; enddo       
-  endif
-
-  if (apply_OBC_u_flather_east) then
-    do j=G%jsd,G%jed ; do I=G%isd,G%ied-1
-      if ((i+G%idg_offset) == G%ieg) then
-        G%bathyT(i+1,j) = G%bathyT(i,j)
-      endif
-    enddo; enddo    
-  endif
-
-  if (apply_OBC_v_flather_north) then
-    do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied
-      if ((j+G%jdg_offset) == G%jeg) then
-        G%bathyT(i,j+1) = G%bathyT(i,j)
-      endif
-    enddo; enddo    
-  endif
-
-  if (apply_OBC_v_flather_south) then
-    do J=G%jsd+1,G%jed ; do i=G%isd,G%ied
-      if ((J+G%jdg_offset) == G%jsg) then
-        G%bathyT(i,j-1) = G%bathyT(i,j)
-      endif
-    enddo; enddo
-  endif
-
+  ! 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
       G%mask2dT(i,j) = 0.0