diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 7657a7d703..88dc987a6b 100644 --- a/src/core/MOM_barotropic.F90 +++ b/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) diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index b3635ca0b6..d67868a67f 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/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 diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 512774522f..d1d49723d3 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/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 diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 1271cd0fdc..c6d9d7bc2b 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1150,9 +1150,13 @@ end subroutine open_boundary_impose_land_mask subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new !< New u values on open boundaries + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new !< On exit, new u values on open boundaries + !! On entry, the old time-level v but + !! including barotropic accelerations. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u_old !< Original unadjusted u - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new !< New v values on open boundaries + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new !< On exit, new v values on open boundaries. + !! On entry, the old time-level v but + !! including barotropic accelerations. real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v_old !< Original unadjusted v real, intent(in) :: dt !< Appropriate timestep ! Local variables @@ -1176,7 +1180,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) do n=1,OBC%number_of_segments segment=>OBC%segment(n) if (.not. segment%on_pe) cycle - if (segment%oblique) call gradient_at_q_points(G,segment,u_old,v_old) + if (segment%oblique) call gradient_at_q_points(G,segment,u_new,v_new) if (segment%direction == OBC_DIRECTION_E) then I=segment%HI%IscB do k=1,nz ; do j=segment%HI%jsc,segment%HI%jec @@ -1187,7 +1191,10 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max) ! outward phase speed rx_avg = (1.0-gamma_u)*segment%rx_normal(I,j,k) + gamma_u*rx_new segment%rx_normal(I,j,k) = rx_avg - segment%normal_vel(I,j,k) = (u_old(I,j,k) + rx_avg*u_new(I-1,j,k)) / (1.0+rx_avg) + ! The new boundary value is interpolated between future interior + ! value, u_new(I-1) and past boundary value but with barotropic + ! accelerations, u_new(I). + segment%normal_vel(I,j,k) = (u_new(I,j,k) + rx_avg*u_new(I-1,j,k)) / (1.0+rx_avg) elseif (segment%oblique) then dhdt = u_old(I-1,j,k)-u_new(I-1,j,k) !old-new dhdx = u_new(I-1,j,k)-u_new(I-2,j,k) !in new time backward sasha for I-1 @@ -1201,16 +1208,14 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) 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 radiation + Cx = min(dhdt*dhdx,rx_max) ! default to normal radiation ! Cy = 0.0 cff = max(dhdx*dhdx,eps) ! if (segment%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 - segment%normal_vel(I,j,k) = ((cff*u_old(I,j,k) + Cx*u_new(I-1,j,k)) - & + segment%normal_vel(I,j,k) = ((cff*u_new(I,j,k) + Cx*u_new(I-1,j,k)) - & (max(Cy,0.0)*segment%grad_normal(J-1,2,k) + min(Cy,0.0)*segment%grad_normal(J,2,k))) / (cff + Cx) elseif (segment%gradient) then segment%normal_vel(I,j,k) = u_new(I-1,j,k) @@ -1221,7 +1226,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) else tau = segment%Tnudge_out endif - segment%normal_vel(I,j,k) = u_new(I,j,k) + dt*tau*(segment%nudged_normal_vel(I,j,k) - u_old(I,j,k)) + segment%normal_vel(I,j,k) = u_new(I,j,k) + dt*tau*(segment%nudged_normal_vel(I,j,k) - u_new(I,j,k)) endif enddo; enddo endif @@ -1236,7 +1241,10 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max) rx_avg = (1.0-gamma_u)*segment%rx_normal(I,j,k) + gamma_u*rx_new segment%rx_normal(I,j,k) = rx_avg - segment%normal_vel(I,j,k) = (u_old(I,j,k) + rx_avg*u_new(I+1,j,k)) / (1.0+rx_avg) + ! The new boundary value is interpolated between future interior + ! value, u_new(I+1) and past boundary value but with barotropic + ! accelerations, u_new(I). + segment%normal_vel(I,j,k) = (u_new(I,j,k) + rx_avg*u_new(I+1,j,k)) / (1.0+rx_avg) elseif (segment%oblique) then dhdt = u_old(I+1,j,k)-u_new(I+1,j,k) !old-new dhdx = u_new(I+1,j,k)-u_new(I+2,j,k) !in new time forward sasha for I+1 @@ -1250,16 +1258,14 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) 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 (segment%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 - segment%normal_vel(I,j,k) = ((cff*u_old(I,j,k) + Cx*u_new(I+1,j,k)) - & + segment%normal_vel(I,j,k) = ((cff*u_new(I,j,k) + Cx*u_new(I+1,j,k)) - & (max(Cy,0.0)*segment%grad_normal(J-1,2,k) + min(Cy,0.0)*segment%grad_normal(J,2,k))) / (cff + Cx) elseif (segment%gradient) then segment%normal_vel(I,j,k) = u_new(I+1,j,k) @@ -1270,7 +1276,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) else tau = segment%Tnudge_out endif - segment%normal_vel(I,j,k) = u_new(I,j,k) + dt*tau*(segment%nudged_normal_vel(I,j,k) - u_old(I,j,k)) + segment%normal_vel(I,j,k) = u_new(I,j,k) + dt*tau*(segment%nudged_normal_vel(I,j,k) - u_new(I,j,k)) endif enddo; enddo endif @@ -1285,7 +1291,10 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max) ry_avg = (1.0-gamma_v)*segment%rx_normal(I,j,k) + gamma_v*ry_new segment%rx_normal(i,J,k) = ry_avg - segment%normal_vel(i,J,k) = (v_old(i,J,k) + ry_avg*v_new(i,J-1,k)) / (1.0+ry_avg) + ! The new boundary value is interpolated between future interior + ! value, v_new(J-1) and past boundary value but with barotropic + ! accelerations, v_new(J). + segment%normal_vel(i,J,k) = (v_new(i,J,k) + ry_avg*v_new(i,J-1,k)) / (1.0+ry_avg) elseif (segment%oblique) then dhdt = v_old(i,J-1,k)-v_new(i,J-1,k) !old-new dhdy = v_new(i,J-1,k)-v_new(i,J-2,k) !in new time backward sasha for J-1 @@ -1299,16 +1308,14 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) 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 (segment%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 - segment%normal_vel(i,J,k) = ((cff*v_old(i,J,k) + Cy*v_new(i,J-1,k)) - & + segment%normal_vel(i,J,k) = ((cff*v_new(i,J,k) + Cy*v_new(i,J-1,k)) - & (max(Cx,0.0)*segment%grad_normal(I-1,2,k) + min(Cx,0.0)*segment%grad_normal(I,2,k))) / (cff + Cy) elseif (segment%gradient) then segment%normal_vel(i,J,k) = v_new(i,J-1,k) @@ -1319,7 +1326,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) else tau = segment%Tnudge_out endif - segment%normal_vel(i,J,k) = v_new(i,J,k) + dt*tau*(segment%nudged_normal_vel(i,J,k) - v_old(i,J,k)) + segment%normal_vel(i,J,k) = v_new(i,J,k) + dt*tau*(segment%nudged_normal_vel(i,J,k) - v_new(i,J,k)) endif enddo; enddo endif @@ -1335,7 +1342,10 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max) ry_avg = (1.0-gamma_v)*segment%rx_normal(I,j,k) + gamma_v*ry_new segment%rx_normal(i,J,k) = ry_avg - segment%normal_vel(i,J,k) = (v_old(i,J,k) + ry_avg*v_new(i,J+1,k)) / (1.0+ry_avg) + ! The new boundary value is interpolated between future interior + ! value, v_new(J+1) and past boundary value but with barotropic + ! accelerations, v_new(J). + segment%normal_vel(i,J,k) = (v_new(i,J,k) + ry_avg*v_new(i,J+1,k)) / (1.0+ry_avg) elseif (segment%oblique) then dhdt = v_old(i,J+1,k)-v_new(i,J+1,k) !old-new dhdy = v_new(i,J+1,k)-v_new(i,J+2,k) !in new time backward sasha for J-1 @@ -1349,16 +1359,14 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) 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 (segment%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 - segment%normal_vel(i,J,k) = ((cff*v_old(i,J,k) + Cy*v_new(i,J+1,k)) - & + segment%normal_vel(i,J,k) = ((cff*v_new(i,J,k) + Cy*v_new(i,J+1,k)) - & (max(Cx,0.0)*segment%grad_normal(I-1,2,k) + min(Cx,0.0)*segment%grad_normal(I,2,k))) / (cff + Cy) elseif (segment%gradient) then segment%normal_vel(i,J,k) = v_new(i,J+1,k) @@ -1369,7 +1377,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) else tau = segment%Tnudge_out endif - segment%normal_vel(i,J,k) = v_new(i,J,k) + dt*tau*(segment%nudged_normal_vel(i,J,k) - v_old(i,J,k)) + segment%normal_vel(i,J,k) = v_new(i,J,k) + dt*tau*(segment%nudged_normal_vel(i,J,k) - v_new(i,J,k)) endif enddo; enddo end if