From 27163baf1156bfa81055742daad080ae23b1b9cd Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Mon, 17 Jul 2017 10:19:42 -0800 Subject: [PATCH 1/5] Revert to ROMSish version of Raymond and Kuo OBC. - Still question about how to do MOM-consistent time levels, etc. --- src/core/MOM_barotropic.F90 | 24 ++++++++---------------- src/core/MOM_open_boundary.F90 | 24 ++++++++---------------- 2 files changed, 16 insertions(+), 32 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index fb63d379b5..066792ff14 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -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) @@ -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) @@ -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) @@ -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_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 1271cd0fdc..63d962d6a3 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1201,14 +1201,12 @@ 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)) - & (max(Cy,0.0)*segment%grad_normal(J-1,2,k) + min(Cy,0.0)*segment%grad_normal(J,2,k))) / (cff + Cx) @@ -1250,14 +1248,12 @@ 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)) - & (max(Cy,0.0)*segment%grad_normal(J-1,2,k) + min(Cy,0.0)*segment%grad_normal(J,2,k))) / (cff + Cx) @@ -1299,14 +1295,12 @@ 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)) - & (max(Cx,0.0)*segment%grad_normal(I-1,2,k) + min(Cx,0.0)*segment%grad_normal(I,2,k))) / (cff + Cy) @@ -1349,14 +1343,12 @@ 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)) - & (max(Cx,0.0)*segment%grad_normal(I-1,2,k) + min(Cx,0.0)*segment%grad_normal(I,2,k))) / (cff + Cy) From 73a5428a41564bda608a855fb28bf3aafbefd3c7 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Mon, 17 Jul 2017 13:35:59 -0800 Subject: [PATCH 2/5] Modification of OBLIQUE OBC. - Use new time in interior value of tangential gradient computation. - Not noticeably better or worse than old way but seems like it should be more consistent. --- src/core/MOM_barotropic.F90 | 16 ++++++++-------- src/core/MOM_open_boundary.F90 | 14 ++++++++------ 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 61f9d9fdd6..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 @@ -2448,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 @@ -2513,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 @@ -2558,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 diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 63d962d6a3..a06a6b85d0 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1176,7 +1176,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, u_old, v_old) if (segment%direction == OBC_DIRECTION_E) then I=segment%HI%IscB do k=1,nz ; do j=segment%HI%jsc,segment%HI%jec @@ -1443,11 +1443,13 @@ subroutine open_boundary_zero_normal_flow(OBC, G, u, v) end subroutine open_boundary_zero_normal_flow !> Calculate the tangential gradient of the normal flow at the boundary q-points. -subroutine gradient_at_q_points(G,segment,uvel,vvel) +subroutine gradient_at_q_points(G, segment, uvel, vvel, uold, vold) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(OBC_segment_type), pointer :: segment !< OBC segment structure real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uvel !< zonal velocity real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vvel !< meridional velocity + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uold !< zonal velocity + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vold !< meridional velocity integer :: i,j,k if (.not. segment%on_pe) return @@ -1463,7 +1465,7 @@ subroutine gradient_at_q_points(G,segment,uvel,vvel) do k=1,G%ke do J=segment%HI%JscB,segment%HI%JecB segment%grad_normal(J,1,k) = (uvel(I-1,j+1,k)-uvel(I-1,j,k)) * G%mask2dBu(I-1,J) - segment%grad_normal(J,2,k) = (uvel(I,j+1,k)-uvel(I,j,k)) * G%mask2dBu(I,J) + segment%grad_normal(J,2,k) = (uold(I,j+1,k)-uold(I,j,k)) * G%mask2dBu(I,J) enddo enddo else ! western segment @@ -1471,7 +1473,7 @@ subroutine gradient_at_q_points(G,segment,uvel,vvel) do k=1,G%ke do J=segment%HI%JscB,segment%HI%JecB segment%grad_normal(J,1,k) = (uvel(I+1,j+1,k)-uvel(I+1,j,k)) * G%mask2dBu(I+1,J) - segment%grad_normal(J,2,k) = (uvel(I,j+1,k)-uvel(I,j,k)) * G%mask2dBu(I,J) + segment%grad_normal(J,2,k) = (uold(I,j+1,k)-uold(I,j,k)) * G%mask2dBu(I,J) enddo enddo endif @@ -1486,7 +1488,7 @@ subroutine gradient_at_q_points(G,segment,uvel,vvel) do k=1,G%ke do I=segment%HI%IscB,segment%HI%IecB segment%grad_normal(I,1,k) = (vvel(i+1,J-1,k)-vvel(i,J-1,k)) * G%mask2dBu(I,J-1) - segment%grad_normal(I,2,k) = (vvel(i+1,J,k)-vvel(i,J,k)) * G%mask2dBu(I,J) + segment%grad_normal(I,2,k) = (vold(i+1,J,k)-vold(i,J,k)) * G%mask2dBu(I,J) enddo enddo else ! south segment @@ -1494,7 +1496,7 @@ subroutine gradient_at_q_points(G,segment,uvel,vvel) do k=1,G%ke do I=segment%HI%IscB,segment%HI%IecB segment%grad_normal(I,1,k) = (vvel(i+1,J+1,k)-vvel(i,J+1,k)) * G%mask2dBu(I,J+1) - segment%grad_normal(I,2,k) = (vvel(i+1,J,k)-vvel(i,J,k)) * G%mask2dBu(I,J) + segment%grad_normal(I,2,k) = (vold(i+1,J,k)-vold(i,J,k)) * G%mask2dBu(I,J) enddo enddo endif From 64787ce87e93b6d78ca345e76cf754091a4d54a5 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 19 Jul 2017 15:57:15 -0800 Subject: [PATCH 3/5] +Still more silly value fixes to OBCs. - Not out of the woods yet, but fewer of my tests are showing trouble. It's down to the 3-d cases plus Kelvin/rotate_BT. --- src/core/MOM_continuity_PPM.F90 | 138 ++++++++++++++++++++++++---- src/core/MOM_dynamics_split_RK2.F90 | 2 + 2 files changed, 122 insertions(+), 18 deletions(-) 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 4d86e7f639..3f82e4b7cc 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 From 5a1ceb94a7bf9c98534e1d43a4e93e8d5c836bad Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 25 Jul 2017 07:04:18 -0800 Subject: [PATCH 4/5] Undoing last commit to MOM_open_boundary. - Changing to new time level for gradient computation. --- src/core/MOM_open_boundary.F90 | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index a06a6b85d0..ed4f7f6186 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1176,7 +1176,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_new, v_new, 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 @@ -1443,13 +1443,11 @@ subroutine open_boundary_zero_normal_flow(OBC, G, u, v) end subroutine open_boundary_zero_normal_flow !> Calculate the tangential gradient of the normal flow at the boundary q-points. -subroutine gradient_at_q_points(G, segment, uvel, vvel, uold, vold) +subroutine gradient_at_q_points(G,segment,uvel,vvel) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(OBC_segment_type), pointer :: segment !< OBC segment structure real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uvel !< zonal velocity real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vvel !< meridional velocity - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uold !< zonal velocity - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vold !< meridional velocity integer :: i,j,k if (.not. segment%on_pe) return @@ -1465,7 +1463,7 @@ subroutine gradient_at_q_points(G, segment, uvel, vvel, uold, vold) do k=1,G%ke do J=segment%HI%JscB,segment%HI%JecB segment%grad_normal(J,1,k) = (uvel(I-1,j+1,k)-uvel(I-1,j,k)) * G%mask2dBu(I-1,J) - segment%grad_normal(J,2,k) = (uold(I,j+1,k)-uold(I,j,k)) * G%mask2dBu(I,J) + segment%grad_normal(J,2,k) = (uvel(I,j+1,k)-uvel(I,j,k)) * G%mask2dBu(I,J) enddo enddo else ! western segment @@ -1473,7 +1471,7 @@ subroutine gradient_at_q_points(G, segment, uvel, vvel, uold, vold) do k=1,G%ke do J=segment%HI%JscB,segment%HI%JecB segment%grad_normal(J,1,k) = (uvel(I+1,j+1,k)-uvel(I+1,j,k)) * G%mask2dBu(I+1,J) - segment%grad_normal(J,2,k) = (uold(I,j+1,k)-uold(I,j,k)) * G%mask2dBu(I,J) + segment%grad_normal(J,2,k) = (uvel(I,j+1,k)-uvel(I,j,k)) * G%mask2dBu(I,J) enddo enddo endif @@ -1488,7 +1486,7 @@ subroutine gradient_at_q_points(G, segment, uvel, vvel, uold, vold) do k=1,G%ke do I=segment%HI%IscB,segment%HI%IecB segment%grad_normal(I,1,k) = (vvel(i+1,J-1,k)-vvel(i,J-1,k)) * G%mask2dBu(I,J-1) - segment%grad_normal(I,2,k) = (vold(i+1,J,k)-vold(i,J,k)) * G%mask2dBu(I,J) + segment%grad_normal(I,2,k) = (vvel(i+1,J,k)-vvel(i,J,k)) * G%mask2dBu(I,J) enddo enddo else ! south segment @@ -1496,7 +1494,7 @@ subroutine gradient_at_q_points(G, segment, uvel, vvel, uold, vold) do k=1,G%ke do I=segment%HI%IscB,segment%HI%IecB segment%grad_normal(I,1,k) = (vvel(i+1,J+1,k)-vvel(i,J+1,k)) * G%mask2dBu(I,J+1) - segment%grad_normal(I,2,k) = (vold(i+1,J,k)-vold(i,J,k)) * G%mask2dBu(I,J) + segment%grad_normal(I,2,k) = (vvel(i+1,J,k)-vvel(i,J,k)) * G%mask2dBu(I,J) enddo enddo endif From 3c073b6aa35d9318f25edb5df86b6b2212723e80 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 25 Jul 2017 14:13:17 -0800 Subject: [PATCH 5/5] +Updating OBC ORLANSKI and OBLIQUE algorithms. - We now use u_new/v_new instead of u_old/v_old at the boundary. These values work because we zeroed out baroclinic accelerations at the boundary and added the barotropic accelerations. - Changes affected OBC answers (for the better). --- src/core/MOM_open_boundary.F90 | 44 +++++++++++++++++++++++----------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index ed4f7f6186..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 @@ -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 @@ -1208,7 +1215,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) cff = max(dhdx*dhdx + dhdy*dhdy, eps) 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) @@ -1219,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 @@ -1234,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 @@ -1255,7 +1265,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) cff = max(dhdx*dhdx + dhdy*dhdy, eps) 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) @@ -1266,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 @@ -1281,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 @@ -1302,7 +1315,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) cff = max(dhdx*dhdx + dhdy*dhdy, eps) 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) @@ -1313,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 @@ -1329,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 @@ -1350,7 +1366,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) cff = max(dhdx*dhdx + dhdy*dhdy, eps) 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) @@ -1361,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