*+Added G%OBCmaskCu and G%OBCmaskCv

  Added two new arrays (OBCmaskCu and OBCmaskCv) to the ocean_grid_type and
dyn_horgrid_type to mask out values over land or at open boundary condition
points.  Without open boundary conditions, these arrays are identical to
mask2dCu and mask2dCv.  These arrays are used in some of the lateral
parameterization modules to zero out certain gradient-dependent fluxes at open
boundary points.  With these changes, the Bering test case solutions no longer
exhibit any dependence on whether DEBUG_OBC is true or false or the value of
OBC_SILLY_THICK, so this commit should help to address some of the issues
discussed in github.com/NOAA-GFDL/issues/208.  All answers are bitwise
identical in the MOM6-examples test cases, but they change for some other tests
that use open boundary conditions more extensively, and there are new arrays in
some transparent types.

src/core/MOM_grid.F90

@@ -90,6 +90,7 @@ module MOM_grid
 
   real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
     mask2dCu, &  !< 0 for boundary points and 1 for ocean points on the u grid [nondim].
+    OBCmaskCu, & !< 0 for boundary or OBC points and 1 for ocean points on the u grid [nondim].
     geoLatCu, &  !< The geographic latitude at u points in degrees of latitude or m.
     geoLonCu, &  !< The geographic longitude at u points in degrees of longitude or m.
     dxCu, &      !< dxCu is delta x at u points [L ~> m].
@@ -102,6 +103,7 @@ module MOM_grid
 
   real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: &
     mask2dCv, &  !< 0 for boundary points and 1 for ocean points on the v grid [nondim].
+    OBCmaskCv, & !< 0 for boundary or OBC points and 1 for ocean points on the v grid [nondim].
     geoLatCv, &  !< The geographic latitude at v points in degrees of latitude or m.
     geoLonCv, &  !< The geographic longitude at v points in degrees of longitude or m.
     dxCv, &      !< dxCv is delta x at v points [L ~> m].
@@ -573,7 +575,9 @@ subroutine allocate_metrics(G)
 
   ALLOC_(G%mask2dT(isd:ied,jsd:jed))      ; G%mask2dT(:,:) = 0.0
   ALLOC_(G%mask2dCu(IsdB:IedB,jsd:jed))   ; G%mask2dCu(:,:) = 0.0
+  ALLOC_(G%OBCmaskCu(IsdB:IedB,jsd:jed))  ; G%OBCmaskCu(:,:) = 0.0
   ALLOC_(G%mask2dCv(isd:ied,JsdB:JedB))   ; G%mask2dCv(:,:) = 0.0
+  ALLOC_(G%OBCmaskCv(isd:ied,JsdB:JedB))  ; G%OBCmaskCv(:,:) = 0.0
   ALLOC_(G%mask2dBu(IsdB:IedB,JsdB:JedB)) ; G%mask2dBu(:,:) = 0.0
   ALLOC_(G%geoLatT(isd:ied,jsd:jed))      ; G%geoLatT(:,:) = 0.0
   ALLOC_(G%geoLatCu(IsdB:IedB,jsd:jed))   ; G%geoLatCu(:,:) = 0.0
@@ -637,8 +641,8 @@ subroutine MOM_grid_end(G)
   DEALLOC_(G%areaCu) ; DEALLOC_(G%IareaCu)
   DEALLOC_(G%areaCv)  ; DEALLOC_(G%IareaCv)
 
-  DEALLOC_(G%mask2dT)  ; DEALLOC_(G%mask2dCu)
-  DEALLOC_(G%mask2dCv) ; DEALLOC_(G%mask2dBu)
+  DEALLOC_(G%mask2dT)  ; DEALLOC_(G%mask2dCu) ; DEALLOC_(G%OBCmaskCu)
+  DEALLOC_(G%mask2dCv) ; DEALLOC_(G%OBCmaskCv) ; DEALLOC_(G%mask2dBu)
 
   DEALLOC_(G%geoLatT)  ; DEALLOC_(G%geoLatCu)
   DEALLOC_(G%geoLatCv) ; DEALLOC_(G%geoLatBu)
@@ -686,6 +690,7 @@ end subroutine MOM_grid_end
 !!
 !! Each location also has a 2D mask indicating whether the entire column is land or ocean.
 !! `mask2dT` is 1 if the column is wet or 0 if the T-cell is land.
-!! `mask2dCu` is 1 if both neighboring column are ocean, and 0 if either is land.
+!! `mask2dCu` is 1 if both neighboring columns are ocean, and 0 if either is land.
+!! `OBCmasku` is 1 if both neighboring columns are ocean, and 0 if either is land of if this is OBC point.
 
 end module MOM_grid

src/core/MOM_open_boundary.F90

@@ -137,7 +137,8 @@ module MOM_open_boundary
   logical :: specified      !< Boundary normal velocity fixed to external value.
   logical :: specified_tan  !< Boundary tangential velocity fixed to external value.
   logical :: specified_grad !< Boundary gradient of tangential velocity fixed to external value.
-  logical :: open           !< Boundary is open for continuity solver.
+  logical :: open           !< Boundary is open for continuity solver, and there are no other
+                            !! parameterized mass fluxes at the open boundary.
   logical :: gradient       !< Zero gradient at boundary.
   logical :: values_needed  !< Whether or not any external OBC fields are needed.
   logical :: u_values_needed      !< Whether or not external u OBC fields are needed.
@@ -1963,16 +1964,16 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US)
       do j=segment%HI%jsd,segment%HI%jed
         if (G%mask2dCu(I,j) == 0) OBC%segnum_u(I,j) = OBC_NONE
         if (segment%direction == OBC_DIRECTION_W) then
-          G%mask2dT(i,j) = 0
+          G%mask2dT(i,j) = 0.0
         else
-          G%mask2dT(i+1,j) = 0
+          G%mask2dT(i+1,j) = 0.0
         endif
       enddo
       do J=segment%HI%JsdB+1,segment%HI%JedB-1
         if (segment%direction == OBC_DIRECTION_W) then
-          G%mask2dCv(i,J) = 0
+          G%mask2dCv(i,J) = 0 ; G%OBCmaskCv(i,J) = 0.0
         else
-          G%mask2dCv(i+1,J) = 0
+          G%mask2dCv(i+1,J) = 0.0 ; G%OBCmaskCv(i+1,J) = 0.0
         endif
       enddo
     else
@@ -1981,21 +1982,38 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US)
       do i=segment%HI%isd,segment%HI%ied
         if (G%mask2dCv(i,J) == 0) OBC%segnum_v(i,J) = OBC_NONE
         if (segment%direction == OBC_DIRECTION_S) then
-          G%mask2dT(i,j) = 0
+          G%mask2dT(i,j) = 0.0
         else
-          G%mask2dT(i,j+1) = 0
+          G%mask2dT(i,j+1) = 0.0
         endif
       enddo
       do I=segment%HI%IsdB+1,segment%HI%IedB-1
         if (segment%direction == OBC_DIRECTION_S) then
-          G%mask2dCu(I,j) = 0
+          G%mask2dCu(I,j) = 0.0 ; G%OBCmaskCu(I,j) = 0.0
         else
-          G%mask2dCu(I,j+1) = 0
+          G%mask2dCu(I,j+1) = 0.0 ; G%OBCmaskCu(I,j+1) = 0.0
         endif
       enddo
     endif
   enddo
 
+  do n=1,OBC%number_of_segments
+    segment=>OBC%segment(n)
+    if (.not. (segment%on_pe .and. segment%open)) cycle
+    ! Set the OBCmask values to help eliminate certain terms at u- or v- OBC points.
+    if (segment%is_E_or_W) then
+      I=segment%HI%IsdB
+      do j=segment%HI%jsd,segment%HI%jed
+        G%OBCmaskCu(I,j) = 0.0
+      enddo
+    else
+      J=segment%HI%JsdB
+      do i=segment%HI%isd,segment%HI%ied
+        G%OBCmaskCv(i,J) = 0.0
+      enddo
+    endif
+  enddo
+
   do n=1,OBC%number_of_segments
     segment=>OBC%segment(n)
     if (.not. segment%on_pe .or. .not. segment%specified) cycle

src/core/MOM_transcribe_grid.F90

@@ -76,6 +76,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
     oG%porous_DavgU(I,j) = dG%porous_DavgU(I+ido,j+jdo) - oG%Z_ref
 
     oG%mask2dCu(I,j) = dG%mask2dCu(I+ido,j+jdo)
+    oG%OBCmaskCu(I,j) = dG%OBCmaskCu(I+ido,j+jdo)
     oG%areaCu(I,j) = dG%areaCu(I+ido,j+jdo)
     oG%IareaCu(I,j) = dG%IareaCu(I+ido,j+jdo)
   enddo ; enddo
@@ -92,6 +93,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
     oG%porous_DavgV(i,J) = dG%porous_DavgV(i+ido,J+jdo) - oG%Z_ref
 
     oG%mask2dCv(i,J) = dG%mask2dCv(i+ido,J+jdo)
+    oG%OBCmaskCv(i,J) = dG%OBCmaskCv(i+ido,J+jdo)
     oG%areaCv(i,J) = dG%areaCv(i+ido,J+jdo)
     oG%IareaCv(i,J) = dG%IareaCv(i+ido,J+jdo)
   enddo ; enddo
@@ -152,6 +154,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
   call pass_vector(oG%dxCu, oG%dyCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(oG%dy_Cu, oG%dx_Cv, oG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(oG%mask2dCu, oG%mask2dCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE)
+  call pass_vector(oG%OBCmaskCu, oG%OBCmaskCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(oG%IareaCu, oG%IareaCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(oG%IareaCu, oG%IareaCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(oG%geoLatCu, oG%geoLatCv, oG%Domain, To_All+Scalar_Pair, CGRID_NE)
@@ -230,6 +233,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
     dG%porous_DavgU(I,j) = oG%porous_DavgU(I+ido,j+jdo) + oG%Z_ref
 
     dG%mask2dCu(I,j) = oG%mask2dCu(I+ido,j+jdo)
+    dG%OBCmaskCu(I,j) = oG%OBCmaskCu(I+ido,j+jdo)
     dG%areaCu(I,j) = oG%areaCu(I+ido,j+jdo)
     dG%IareaCu(I,j) = oG%IareaCu(I+ido,j+jdo)
   enddo ; enddo
@@ -246,6 +250,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
     dG%porous_DavgV(i,J) = oG%porous_DavgU(i+ido,J+jdo) + oG%Z_ref
 
     dG%mask2dCv(i,J) = oG%mask2dCv(i+ido,J+jdo)
+    dG%OBCmaskCv(i,J) = oG%OBCmaskCv(i+ido,J+jdo)
     dG%areaCv(i,J) = oG%areaCv(i+ido,J+jdo)
     dG%IareaCv(i,J) = oG%IareaCv(i+ido,J+jdo)
   enddo ; enddo
@@ -307,6 +312,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
   call pass_vector(dG%dxCu, dG%dyCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(dG%dy_Cu, dG%dx_Cv, dG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(dG%mask2dCu, dG%mask2dCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE)
+  call pass_vector(dG%OBCmaskCu, dG%OBCmaskCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(dG%IareaCu, dG%IareaCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(dG%IareaCu, dG%IareaCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE)
   call pass_vector(dG%geoLatCu, dG%geoLatCv, dG%Domain, To_All+Scalar_Pair, CGRID_NE)

src/framework/MOM_dyn_horgrid.F90

@@ -87,6 +87,7 @@ module MOM_dyn_horgrid
 
   real, allocatable, dimension(:,:) :: &
     mask2dCu, &  !< 0 for boundary points and 1 for ocean points on the u grid [nondim].
+    OBCmaskCu, & !< 0 for boundary or OBC points and 1 for ocean points on the u grid [nondim].
     geoLatCu, &  !< The geographic latitude at u points [degrees of latitude] or [m].
     geoLonCu, &  !< The geographic longitude at u points [degrees of longitude] or [m].
     dxCu, &      !< dxCu is delta x at u points [L ~> m].
@@ -99,6 +100,7 @@ module MOM_dyn_horgrid
 
   real, allocatable, dimension(:,:) :: &
     mask2dCv, &  !< 0 for boundary points and 1 for ocean points on the v grid [nondim].
+    OBCmaskCv, & !< 0 for boundary or OBC points and 1 for ocean points on the v grid [nondim].
     geoLatCv, &  !< The geographic latitude at v points [degrees of latitude] or [m].
     geoLonCv, &  !< The geographic longitude at v points [degrees of longitude] or [m].
     dxCv, &      !< dxCv is delta x at v points [L ~> m].
@@ -250,6 +252,8 @@ subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel)
   allocate(G%mask2dCu(IsdB:IedB,jsd:jed), source=0.0)
   allocate(G%mask2dCv(isd:ied,JsdB:JedB), source=0.0)
   allocate(G%mask2dBu(IsdB:IedB,JsdB:JedB), source=0.0)
+  allocate(G%OBCmaskCu(IsdB:IedB,jsd:jed), source=0.0)
+  allocate(G%OBCmaskCv(isd:ied,JsdB:JedB), source=0.0)
   allocate(G%geoLatT(isd:ied,jsd:jed), source=0.0)
   allocate(G%geoLatCu(IsdB:IedB,jsd:jed), source=0.0)
   allocate(G%geoLatCv(isd:ied,JsdB:JedB), source=0.0)
@@ -331,6 +335,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns)
   call rotate_array_pair(G_in%dx_Cv, G_in%dy_Cu, turns, G%dx_Cv, G%dy_Cu)
 
   call rotate_array_pair(G_in%mask2dCu, G_in%mask2dCv, turns, G%mask2dCu, G%mask2dCv)
+  call rotate_array_pair(G_in%OBCmaskCu, G_in%OBCmaskCv, turns, G%OBCmaskCu, G%OBCmaskCv)
   call rotate_array_pair(G_in%areaCu, G_in%areaCv, turns, G%areaCu, G%areaCv)
   call rotate_array_pair(G_in%IareaCu, G_in%IareaCv, turns, G%IareaCu, G%IareaCv)
 
@@ -501,8 +506,8 @@ subroutine destroy_dyn_horgrid(G)
   deallocate(G%areaCu) ; deallocate(G%IareaCu)
   deallocate(G%areaCv)  ; deallocate(G%IareaCv)
 
-  deallocate(G%mask2dT)  ; deallocate(G%mask2dCu)
-  deallocate(G%mask2dCv) ; deallocate(G%mask2dBu)
+  deallocate(G%mask2dT)  ; deallocate(G%mask2dCu) ; deallocate(G%OBCmaskCu)
+  deallocate(G%mask2dCv) ; deallocate(G%OBCmaskCv) ; deallocate(G%mask2dBu)
 
   deallocate(G%geoLatT)  ; deallocate(G%geoLatCu)
   deallocate(G%geoLatCv) ; deallocate(G%geoLatBu)

src/initialization/MOM_grid_initialize.F90

@@ -1206,6 +1206,8 @@ subroutine initialize_masks(G, PF, US)
     else
       G%mask2dCu(I,j) = 1.0
     endif
+    ! This mask may be revised later after the open boundary positions are specified.
+    G%OBCmaskCu(I,j) = G%mask2dCu(I,j)
   enddo ; enddo
 
   do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied
@@ -1214,6 +1216,8 @@ subroutine initialize_masks(G, PF, US)
     else
       G%mask2dCv(i,J) = 1.0
     endif
+    ! This mask may be revised later after the open boundary positions are specified.
+    G%OBCmaskCv(i,J) = G%mask2dCv(i,J)
   enddo ; enddo
 
   do J=G%jsd,G%jed-1 ; do I=G%isd,G%ied-1
@@ -1229,12 +1233,14 @@ subroutine initialize_masks(G, PF, US)
   call pass_vector(G%mask2dCu, G%mask2dCv, G%Domain, To_All+Scalar_Pair, CGRID_NE)
 
   do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB
+    ! This open face length may be revised later.
     G%dy_Cu(I,j) = G%mask2dCu(I,j) * G%dyCu(I,j)
     G%areaCu(I,j) = G%dxCu(I,j) * G%dy_Cu(I,j)
     G%IareaCu(I,j) = G%mask2dCu(I,j) * Adcroft_reciprocal(G%areaCu(I,j))
   enddo ; enddo
 
   do J=G%JsdB,G%JedB ; do i=G%isd,G%ied
+    ! This open face length may be revised later.
     G%dx_Cv(i,J) = G%mask2dCv(i,J) * G%dxCv(i,J)
     G%areaCv(i,J) = G%dyCv(i,J) * G%dx_Cv(i,J)
     G%IareaCv(i,J) = G%mask2dCv(i,J) * Adcroft_reciprocal(G%areaCv(i,J))

src/parameterizations/lateral/MOM_MEKE.F90

@@ -462,7 +462,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       !$OMP parallel do default(shared)
       do j=js-1,je+1 ; do I=is-2,ie+1
         ! MEKE_uflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].
-        MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * G%mask2dCu(I,j)) * &
+        MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * G%OBCmaskCu(I,j)) * &
             (MEKE%MEKE(i+1,j) - MEKE%MEKE(i,j))
       ! This would have units of [R Z L2 T-2 ~> kg s-2]
       ! MEKE_uflux(I,j) = ((G%dy_Cu(I,j)*G%IdxCu(I,j)) * &
@@ -472,7 +472,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       !$OMP parallel do default(shared)
       do J=js-2,je+1 ; do i=is-1,ie+1
         ! MEKE_vflux is used here as workspace with units of [L2 T-2 ~> m2 s-2].
-        MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * G%mask2dCv(i,J)) * &
+        MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * G%OBCmaskCv(i,J)) * &
             (MEKE%MEKE(i,j+1) - MEKE%MEKE(i,j))
       ! This would have units of [R Z L2 T-2 ~> kg s-2]
       ! MEKE_vflux(i,J) = ((G%dx_Cv(i,J)*G%IdyCv(i,J)) * &

src/parameterizations/lateral/MOM_interface_filter.F90

@@ -284,7 +284,7 @@ subroutine filter_interface(h, e, Lsm2_u, Lsm2_v, uhD, vhD, G, GV, US, halo_size
     do I=is-1,ie ; uhtot(I,j) = 0.0 ; enddo
     do K=nz,2,-1
       do I=is-1,ie
-        Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%mask2dCu(I,j)
+        Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%OBCmaskCu(I,j)
 
         Sfn_est = (Lsm2_u(I,j)*G%dy_Cu(I,j)) * (GV%Z_to_H * Slope)
 
@@ -316,7 +316,7 @@ subroutine filter_interface(h, e, Lsm2_u, Lsm2_v, uhD, vhD, G, GV, US, halo_size
     do i=is,ie ; vhtot(i,J) = 0.0 ; enddo
     do K=nz,2,-1
       do i=is,ie
-        Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%mask2dCv(i,J)
+        Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%OBCmaskCv(i,J)
 
         Sfn_est = (Lsm2_v(i,J)*G%dx_Cv(i,J)) * (GV%Z_to_H * Slope)