*Merge branch 'user/ksh/open_bc' of https://github.com/ESMG/MOM6 into…

… ESMG-user/ksh/open_bc

- Changes answers for circle_obcs

src/core/MOM_barotropic.F90

@@ -2404,8 +2404,8 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans,
         elseif (OBC%segment(OBC%segnum_u(I,j))%oblique) then
           grad(I,J) = (ubt_old(I,j+1) - ubt_old(I,j)) * G%mask2dBu(I,J)
           grad(I,J-1) = (ubt_old(I,j) - ubt_old(I,j-1)) * G%mask2dBu(I,J-1)
-          grad(I-1,J) = (ubt_old(I-1,j+1) - ubt_old(I-1,j)) * G%mask2dBu(I-1,J)
-          grad(I-1,J-1) = (ubt_old(I-1,j) - ubt_old(I-1,j-1)) * G%mask2dBu(I-1,J-1)
+          grad(I-1,J) = (ubt(I-1,j+1) - ubt(I-1,j)) * G%mask2dBu(I-1,J)
+          grad(I-1,J-1) = (ubt(I-1,j) - ubt(I-1,j-1)) * G%mask2dBu(I-1,J-1)
           dhdt = ubt_old(I-1,j)-ubt(I-1,j) !old-new
           dhdx = ubt(I-1,j)-ubt(I-2,j) !in new time backward sasha for I-1
 !         if (OBC%segment(OBC%segnum_u(I,j))%oblique) then
@@ -2418,14 +2418,12 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans,
             endif
 !         endif
           if (dhdt*dhdx < 0.0) dhdt = 0.0
-          if (dhdx == 0.0) dhdx=eps  ! avoid segv
-          Cx = min(dhdt/dhdx,rx_max) ! default to normal flow only
+          Cx = min(dhdt*dhdx,rx_max) ! default to normal flow only
 !         Cy = 0
           cff = max(dhdx*dhdx, eps)
 !         if (OBC%segment(OBC%segnum_u(I,j))%oblique) then
             cff = max(dhdx*dhdx + dhdy*dhdy, eps)
-            if (dhdy==0.) dhdy=eps ! avoid segv
-            Cy = min(cff, max(dhdt/dhdy, -cff))
+            Cy = min(cff, max(dhdt*dhdy, -cff))
 !         endif
           ubt(I,j) = ((cff*ubt_old(I,j) + Cx*ubt(I-1,j)) - &
               (max(Cy,0.0)*grad(I,J-1) + min(Cy,0.0)*grad(I,J))) / (cff + Cx)
@@ -2450,8 +2448,8 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans,
         elseif (OBC%segment(OBC%segnum_u(I,j))%oblique) then
           grad(I,J) = (ubt_old(I,j+1) - ubt_old(I,j)) * G%mask2dBu(I,J)
           grad(I,J-1) = (ubt_old(I,j) - ubt_old(I,j-1)) * G%mask2dBu(I,J-1)
-          grad(I+1,J) = (ubt_old(I+1,j+1) - ubt_old(I+1,j)) * G%mask2dBu(I+1,J)
-          grad(I+1,J-1) = (ubt_old(I+1,j) - ubt_old(I+1,j-1)) * G%mask2dBu(I+1,J-1)
+          grad(I+1,J) = (ubt(I+1,j+1) - ubt(I+1,j)) * G%mask2dBu(I+1,J)
+          grad(I+1,J-1) = (ubt(I+1,j) - ubt(I+1,j-1)) * G%mask2dBu(I+1,J-1)
           dhdt = ubt_old(I+1,j)-ubt(I+1,j) !old-new
           dhdx = ubt(I+1,j)-ubt(I+2,j) !in new time backward sasha for I+1
 !         if (OBC%segment(OBC%segnum_u(I,j))%oblique) then
@@ -2464,14 +2462,12 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans,
             endif
 !         endif
           if (dhdt*dhdx < 0.0) dhdt = 0.0
-          if (dhdx == 0.0) dhdx=eps  ! avoid segv
-          Cx = min(dhdt/dhdx,rx_max) ! default to normal flow only
+          Cx = min(dhdt*dhdx,rx_max) ! default to normal flow only
 !         Cy = 0
           cff = max(dhdx*dhdx, eps)
 !         if (OBC%segment(OBC%segnum_u(I,j))%oblique) then
             cff = max(dhdx*dhdx + dhdy*dhdy, eps)
-            if (dhdy==0.) dhdy=eps ! avoid segv
-            Cy = min(cff,max(dhdt/dhdy,-cff))
+            Cy = min(cff,max(dhdt*dhdy,-cff))
 !         endif
           ubt(I,j) = ((cff*ubt_old(I,j) + Cx*ubt(I+1,j)) - &
               (max(Cy,0.0)*grad(I,J-1) + min(Cy,0.0)*grad(I,J))) / (cff + Cx)
@@ -2517,8 +2513,8 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans,
         elseif (OBC%segment(OBC%segnum_v(i,J))%oblique) then
           grad(I,J) = (vbt_old(i+1,J) - vbt_old(i,J)) * G%mask2dBu(I,J)
           grad(I-1,J) = (vbt_old(i,J) - vbt_old(i-1,J)) * G%mask2dBu(I-1,J)
-          grad(I,J-1) = (vbt_old(i+1,J-1) - vbt_old(i,J-1)) * G%mask2dBu(I,J-1)
-          grad(I-1,J-1) = (vbt_old(i,J-1) - vbt_old(i-1,J-1)) * G%mask2dBu(I-1,J-1)
+          grad(I,J-1) = (vbt(i+1,J-1) - vbt(i,J-1)) * G%mask2dBu(I,J-1)
+          grad(I-1,J-1) = (vbt(i,J-1) - vbt(i-1,J-1)) * G%mask2dBu(I-1,J-1)
           dhdt = vbt_old(i,J-1)-vbt(i,J-1) !old-new
           dhdy = vbt(i,J-1)-vbt(i,J-2) !in new time backward sasha for J-1
 !         if (OBC%segment(OBC%segnum_v(i,J))%oblique) then
@@ -2531,14 +2527,12 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans,
             endif
 !         endif
           if (dhdt*dhdy < 0.0) dhdt = 0.0
-          if (dhdy == 0.0) dhdy=eps  ! avoid segv
-          Cy = min(dhdt/dhdy,rx_max) ! default to normal flow only
+          Cy = min(dhdt*dhdy,rx_max) ! default to normal flow only
 !         Cx = 0
           cff = max(dhdy*dhdy, eps)
 !         if (OBC%segment(OBC%segnum_v(i,J))%oblique) then
             cff = max(dhdx*dhdx + dhdy*dhdy, eps)
-            if (dhdx==0.) dhdx=eps ! avoid segv
-            Cx = min(cff,max(dhdt/dhdx,-cff))
+            Cx = min(cff,max(dhdt*dhdx,-cff))
 !         endif
           vbt(i,J) = ((cff*vbt_old(i,J) + Cy*vbt(i,J-1)) - &
             (max(Cx,0.0)*grad(I-1,J) + min(Cx,0.0)*grad(I,J))) / (cff + Cy)
@@ -2564,8 +2558,8 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans,
         elseif (OBC%segment(OBC%segnum_v(i,J))%oblique) then
           grad(I,J) = (vbt_old(i+1,J) - vbt_old(i,J)) * G%mask2dBu(I,J)
           grad(I-1,J) = (vbt_old(i,J) - vbt_old(i-1,J)) * G%mask2dBu(I-1,J)
-          grad(I,J+1) = (vbt_old(i+1,J+1) - vbt_old(i,J+1)) * G%mask2dBu(I,J+1)
-          grad(I-1,J+1) = (vbt_old(i,J+1) - vbt_old(i-1,J+1)) * G%mask2dBu(I-1,J+1)
+          grad(I,J+1) = (vbt(i+1,J+1) - vbt(i,J+1)) * G%mask2dBu(I,J+1)
+          grad(I-1,J+1) = (vbt(i,J+1) - vbt(i-1,J+1)) * G%mask2dBu(I-1,J+1)
           dhdt = vbt_old(i,J+1)-vbt(i,J+1) !old-new
           dhdy = vbt(i,J+1)-vbt(i,J+2) !in new time backward sasha for J+1
 !         if (OBC%segment(OBC%segnum_v(i,J))%oblique) then
@@ -2578,14 +2572,12 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans,
             endif
 !         endif
           if (dhdt*dhdy < 0.0) dhdt = 0.0
-          if (dhdy == 0.0) dhdy=eps  ! avoid segv
-          Cy = min(dhdt/dhdy,rx_max) ! default to normal flow only
+          Cy = min(dhdt*dhdy,rx_max) ! default to normal flow only
 !         Cx = 0
           cff = max(dhdy*dhdy, eps)
 !         if (OBC%segment(OBC%segnum_v(i,J))%oblique) then
             cff = max(dhdx*dhdx + dhdy*dhdy, eps)
-            if (dhdx==0.) dhdx=eps ! avoid segv
-            Cx = min(cff,max(dhdt/dhdx,-cff))
+            Cx = min(cff,max(dhdt*dhdx,-cff))
 !         endif
           vbt(i,J) = ((cff*vbt_old(i,J) + Cy*vbt(i,J+1)) - &
             (max(Cx,0.0)*grad(I-1,J) + min(Cx,0.0)*grad(I,J))) / (cff + Cy)

src/core/MOM_continuity_PPM.F90

@@ -363,7 +363,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, &
       enddo ; enddo
     else
       call PPM_reconstruction_x(h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), G, LB, &
-                                2.0*GV%Angstrom, CS%monotonic, simple_2nd=CS%simple_2nd)
+                                2.0*GV%Angstrom, CS%monotonic, simple_2nd=CS%simple_2nd, OBC=OBC)
     endif
     do I=ish-1,ieh ; visc_rem(I,k) = 1.0 ; enddo
   enddo
@@ -375,17 +375,19 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, &
         I=segment%HI%IsdB
         do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
           h_in(i+1,j,k) = h_in(i,j,k)
-          h_L(i+1,j,k) = h_in(i,j,k)
-          h_R(i+1,j,k) = h_in(i,j,k)
-          h_R(i,j,k) = h_in(i,j,k)
+!         h_L(i+1,j,k) = h_in(i,j,k)
+!         h_R(i+1,j,k) = h_in(i,j,k)
+!         h_L(i,j,k) = h_in(i,j,k)
+!         h_R(i,j,k) = h_in(i,j,k)
         enddo ; enddo
       elseif (segment%direction == OBC_DIRECTION_W) then
         I=segment%HI%IsdB
         do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
             h_in(i,j,k) = h_in(i+1,j,k)
-            h_L(i,j,k) = h_in(i+1,j,k)
-            h_R(i,j,k) = h_in(i+1,j,k)
-            h_L(i+1,j,k) = h_in(i+1,j,k)
+!           h_L(i,j,k) = h_in(i+1,j,k)
+!           h_R(i,j,k) = h_in(i+1,j,k)
+!           h_L(i+1,j,k) = h_in(i+1,j,k)
+!           h_R(i+1,j,k) = h_in(i+1,j,k)
         enddo ; enddo
       endif
     enddo
@@ -1138,7 +1140,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, &
       enddo ; enddo
     else
       call PPM_reconstruction_y(h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), G, LB, &
-                                2.0*GV%Angstrom, CS%monotonic, simple_2nd=CS%simple_2nd)
+                                2.0*GV%Angstrom, CS%monotonic, simple_2nd=CS%simple_2nd, OBC=OBC)
     endif
     do i=ish,ieh ; visc_rem(i,k) = 1.0 ; enddo
   enddo
@@ -1150,17 +1152,19 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, &
         J=segment%HI%JsdB
         do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
           h_in(i,j+1,k) = h_in(i,j,k)
-          h_L(i,j+1,k) = h_in(i,j,k)
-          h_R(i,j+1,k) = h_in(i,j,k)
-          h_R(i,j,k) = h_in(i,j,k)
+!         h_L(i,j+1,k) = h_in(i,j,k)
+!         h_R(i,j+1,k) = h_in(i,j,k)
+!         h_L(i,j,k) = h_in(i,j,k)
+!         h_R(i,j,k) = h_in(i,j,k)
         enddo ; enddo
       elseif (segment%direction == OBC_DIRECTION_S) then
         J=segment%HI%JsdB
         do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
           h_in(i,j,k) = h_in(i,j+1,k)
-          h_L(i,j,k) = h_in(i,j+1,k)
-          h_R(i,j,k) = h_in(i,j+1,k)
-          h_L(i,j+1,k) = h_in(i,j+1,k)
+!         h_L(i,j,k) = h_in(i,j+1,k)
+!         h_R(i,j,k) = h_in(i,j+1,k)
+!         h_L(i,j+1,k) = h_in(i,j+1,k)
+!         h_R(i,j+1,k) = h_in(i,j+1,k)
         enddo ; enddo
       endif
     enddo
@@ -1824,7 +1828,7 @@ subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0,
 end subroutine set_merid_BT_cont
 
 !> Calculates left/right edge values for PPM reconstruction.
-subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd)
+subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd, OBC)
   type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
   real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness, in H.
   real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_L  !< Left thickness in the
@@ -1840,6 +1844,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
   logical, optional,                 intent(in)  :: simple_2nd !< If true, use the
                     !! arithmetic mean thicknesses as the default edge values
                     !! for a simple 2nd order scheme.
+  type(ocean_OBC_type), pointer,     optional    :: OBC !< Open boundaries control structure.
 
   ! Local variables with useful mnemonic names.
   real, dimension(SZI_(G),SZJ_(G))  :: slp ! The slopes.
@@ -1848,10 +1853,18 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
   real :: dMx, dMn
   logical :: use_CW84, use_2nd
   character(len=256) :: mesg
-  integer :: i, j, isl, iel, jsl, jel, stencil
+  integer :: i, j, isl, iel, jsl, jel, n, stencil
+  logical :: local_open_BC
+  type(OBC_segment_type), pointer :: segment
 
   use_CW84 = .false. ; if (present(monotonic)) use_CW84 = monotonic
   use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
+
+  local_open_BC = .false.
+  if (present(OBC)) then ; if (associated(OBC)) then
+    local_open_BC = OBC%open_u_BCs_exist_globally
+  endif ; endif
+
   isl = LB%ish-1 ; iel = LB%ieh+1 ; jsl = LB%jsh ; jel = LB%jeh
 
   ! This is the stencil of the reconstruction, not the scheme overall.
@@ -1892,6 +1905,21 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
       endif
     enddo; enddo
 
+    if (local_open_BC) then
+      do n=1, OBC%number_of_segments
+        segment => OBC%segment(n)
+        if (.not. segment%on_pe .or. segment%specified) cycle
+        if (segment%direction == OBC_DIRECTION_E .or. &
+            segment%direction == OBC_DIRECTION_W) then
+          I=segment%HI%IsdB
+          do j=segment%HI%jsd,segment%HI%jed
+            slp(i+1,j) = 0.0
+            slp(i,j) = 0.0
+          enddo
+        endif
+      enddo
+    endif
+
     do j=jsl,jel ; do i=isl,iel
       ! Neighboring values should take into account any boundaries.  The 3
       ! following sets of expressions are equivalent.
@@ -1905,6 +1933,32 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
     enddo; enddo
   endif
 
+  if (local_open_BC) then
+    do n=1, OBC%number_of_segments
+      segment => OBC%segment(n)
+      if (.not. segment%on_pe .or. segment%specified) cycle
+      if (segment%direction == OBC_DIRECTION_E) then
+        I=segment%HI%IsdB
+        do j=segment%HI%jsd,segment%HI%jed
+!         h_in(i+1,j) = h_in(i,j)
+          h_L(i+1,j) = h_in(i,j)
+          h_R(i+1,j) = h_in(i,j)
+          h_L(i,j) = h_in(i,j)
+          h_R(i,j) = h_in(i,j)
+        enddo
+      elseif (segment%direction == OBC_DIRECTION_W) then
+        I=segment%HI%IsdB
+        do j=segment%HI%jsd,segment%HI%jed
+!         h_in(i,j) = h_in(i+1,j)
+          h_L(i,j) = h_in(i+1,j)
+          h_R(i,j) = h_in(i+1,j)
+          h_L(i+1,j) = h_in(i+1,j)
+          h_R(i+1,j) = h_in(i+1,j)
+        enddo
+      endif
+    enddo
+  endif
+
   if (use_CW84) then
     call PPM_limit_CW84(h_in, h_L, h_R, G, isl, iel, jsl, jel)
   else
@@ -1915,7 +1969,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
 end subroutine PPM_reconstruction_x
 
 !> Calculates left/right edge values for PPM reconstruction.
-subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd)
+subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_2nd, OBC)
   type(ocean_grid_type),             intent(in)  :: G    !< Ocean's grid structure.
   real, dimension(SZI_(G),SZJ_(G)),  intent(in)  :: h_in !< Layer thickness, in H.
   real, dimension(SZI_(G),SZJ_(G)),  intent(out) :: h_L  !< Left thickness in the
@@ -1931,6 +1985,7 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
   logical, optional,                 intent(in)  :: simple_2nd !< If true, use the
                     !! arithmetic mean thicknesses as the default edge values
                     !! for a simple 2nd order scheme.
+  type(ocean_OBC_type), pointer,     optional    :: OBC !< Open boundaries control structure.
 
   ! Local variables with useful mnemonic names.
   real, dimension(SZI_(G),SZJ_(G))  :: slp ! The slopes.
@@ -1939,10 +1994,18 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
   real :: dMx, dMn
   logical :: use_CW84, use_2nd
   character(len=256) :: mesg
-  integer :: i, j, isl, iel, jsl, jel, stencil
+  integer :: i, j, isl, iel, jsl, jel, n, stencil
+  logical :: local_open_BC
+  type(OBC_segment_type), pointer :: segment
 
   use_CW84 = .false. ; if (present(monotonic)) use_CW84 = monotonic
   use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
+
+  local_open_BC = .false.
+  if (present(OBC)) then ; if (associated(OBC)) then
+    local_open_BC = OBC%open_u_BCs_exist_globally
+  endif ; endif
+
   isl = LB%ish ; iel = LB%ieh ; jsl = LB%jsh-1 ; jel = LB%jeh+1
 
   ! This is the stencil of the reconstruction, not the scheme overall.
@@ -1983,6 +2046,21 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
       endif
     enddo ; enddo
 
+    if (local_open_BC) then
+      do n=1, OBC%number_of_segments
+        segment => OBC%segment(n)
+        if (.not. segment%on_pe .or. segment%specified) cycle
+        if (segment%direction == OBC_DIRECTION_S .or. &
+            segment%direction == OBC_DIRECTION_N) then
+          J=segment%HI%JsdB
+          do i=segment%HI%isd,segment%HI%ied
+            slp(i,j+1) = 0.0
+            slp(i,j) = 0.0
+          enddo
+        endif
+      enddo
+    endif
+
     do j=jsl,jel ; do i=isl,iel
       ! Neighboring values should take into account any boundaries.  The 3
       ! following sets of expressions are equivalent.
@@ -1994,6 +2072,30 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_
     enddo ; enddo
   endif
 
+  if (local_open_BC) then
+    do n=1, OBC%number_of_segments
+      segment => OBC%segment(n)
+      if (.not. segment%on_pe .or. segment%specified) cycle
+      if (segment%direction == OBC_DIRECTION_N) then
+        J=segment%HI%JsdB
+        do i=segment%HI%isd,segment%HI%ied
+          h_L(i,j+1) = h_in(i,j)
+          h_R(i,j+1) = h_in(i,j)
+          h_L(i,j) = h_in(i,j)
+          h_R(i,j) = h_in(i,j)
+        enddo
+      elseif (segment%direction == OBC_DIRECTION_S) then
+        J=segment%HI%JsdB
+        do i=segment%HI%isd,segment%HI%ied
+          h_L(i,j) = h_in(i,j+1)
+          h_R(i,j) = h_in(i,j+1)
+          h_L(i,j+1) = h_in(i,j+1)
+          h_R(i,j+1) = h_in(i,j+1)
+        enddo
+      endif
+    enddo
+  endif
+
   if (use_CW84) then
     call PPM_limit_CW84(h_in, h_L, h_R, G, isl, iel, jsl, jel)
   else

src/core/MOM_dynamics_split_RK2.F90

@@ -50,6 +50,7 @@ module MOM_dynamics_split_RK2
 use MOM_MEKE_types,            only : MEKE_type
 use MOM_open_boundary,         only : ocean_OBC_type, radiation_open_bdry_conds
 use MOM_open_boundary,         only : open_boundary_zero_normal_flow
+use MOM_open_boundary,         only : open_boundary_test_extern_h
 use MOM_PressureForce,         only : PressureForce, PressureForce_init, PressureForce_CS
 use MOM_set_visc,              only : set_viscous_ML, set_visc_CS
 use MOM_tidal_forcing,         only : tidal_forcing_init, tidal_forcing_CS
@@ -321,6 +322,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
   endif
 
   if (associated(CS%OBC)) then
+!   call open_boundary_test_extern_h(G, CS%OBC, h)
     do k=1,nz ; do j=js-1,je+1 ; do I=is-2,ie+1
       u_old_rad_OBC(I,j,k) = u_av(I,j,k)
     enddo ; enddo ; enddo