*+Correct oblique OBC restarts

src/core/MOM_barotropic.F90

@@ -1061,6 +1061,30 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     enddo ; enddo
   enddo
 
+  if (apply_OBCs) then
+    do n=1,OBC%number_of_segments
+      if (.not. OBC%segment(n)%on_pe) cycle
+      I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
+      if (OBC%segment(n)%is_N_or_S .and. (J >= Jsq-1) .and. (J <= Jeq+1)) then
+        do i = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+2,OBC%segment(n)%HI%ied)
+          if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
+            gtot_S(i,j+1) = gtot_S(i,j)
+          else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
+            gtot_N(i,j) = gtot_N(i,j+1)
+          endif
+        enddo
+      elseif (OBC%segment(n)%is_E_or_W .and. (I >= Isq-1) .and. (I <= Ieq+1)) then
+        do j = max(Jsq-1,OBC%segment(n)%HI%jsd), min(Jeq+2,OBC%segment(n)%HI%jed)
+          if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
+            gtot_W(i+1,j) = gtot_W(i,j)
+          else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
+            gtot_E(i,j) = gtot_E(i+1,j)
+          endif
+        enddo
+      endif
+    enddo
+  endif
+
   if (CS%tides) then
     call tidal_forcing_sensitivity(G, CS%tides_CSp, det_de)
     if (CS%tidal_sal_bug) then

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_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)