From 321f9c85dcd8c595677143b744185733938720be Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Mon, 24 Feb 2025 21:49:40 -0900 Subject: [PATCH 1/2] Adding a flag to turn off the bug fix (mostly) --- src/core/MOM_barotropic.F90 | 66 ++++++++++++++++++++++++---- src/core/MOM_boundary_update.F90 | 3 -- src/core/MOM_dynamics_split_RK2.F90 | 3 +- src/core/MOM_dynamics_split_RK2b.F90 | 2 +- src/core/MOM_open_boundary.F90 | 25 ++++++++--- src/tracer/MOM_tracer_advect.F90 | 31 +++++++++++++ 6 files changed, 111 insertions(+), 19 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index ccb49e7036..0948580e9f 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -156,6 +156,7 @@ module MOM_barotropic real, allocatable :: frhatu1(:,:,:) !< Predictor step values of frhatu stored for diagnostics [nondim] real, allocatable :: frhatv1(:,:,:) !< Predictor step values of frhatv stored for diagnostics [nondim] + real, allocatable :: IareaT_OBCmask(:,:) !< If non-zero, work on given points [L-2 ~> m-2]. type(BT_OBC_type) :: BT_OBC !< A structure with all of this modules fields !! for applying open boundary conditions. @@ -290,6 +291,8 @@ module MOM_barotropic !! consistent with tidal self-attraction and loading !! used within the barotropic solver logical :: wt_uv_bug = .true. !< If true, recover a bug that wt_[uv] that is not normalized. + logical :: exterior_OBC_bug = .true. !< If true, recover a bug with boundary conditions + !! inside the domain. type(time_type), pointer :: Time => NULL() !< A pointer to the ocean models clock. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate !! the timing of diagnostic output. @@ -1961,7 +1964,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo !$OMP do do j=jsv-1,jev+1 ; do i=isv-1,iev+1 - eta_pred(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT(i,j) * & + eta_pred(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT_OBCmask(i,j) * & ((uhbt_int(I-1,j) - uhbt_int(I,j)) + (vhbt_int(i,J-1) - vhbt_int(i,J))) enddo ; enddo elseif (use_BT_cont) then @@ -1975,13 +1978,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo !$OMP do do j=jsv-1,jev+1 ; do i=isv-1,iev+1 - eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & + eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT_OBCmask(i,j)) * & ((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J))) enddo ; enddo else !$OMP do do j=jsv-1,jev+1 ; do i=isv-1,iev+1 - eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & + eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT_OBCmask(i,j)) * & (((Datu(I-1,j)*ubt(I-1,j) + uhbt0(I-1,j)) - & (Datu(I,j)*ubt(I,j) + uhbt0(I,j))) + & ((Datv(i,J-1)*vbt(i,J-1) + vhbt0(i,J-1)) - & @@ -2488,14 +2491,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (integral_BT_cont) then !$OMP do do j=jsv,jev ; do i=isv,iev - eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT(i,j) * & + eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT_OBCmask(i,j) * & ((uhbt_int(I-1,j) - uhbt_int(I,j)) + (vhbt_int(i,J-1) - vhbt_int(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo else !$OMP do do j=jsv,jev ; do i=isv,iev - eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & + eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT_OBCmask(i,j)) * & ((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo @@ -3283,7 +3286,7 @@ subroutine set_up_BT_OBC(OBC, eta, SpV_avg, BT_OBC, BT_Domain, G, GV, US, CS, MS type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(barotropic_CS), intent(in) :: CS !< Barotropic control structure + type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure integer, intent(in) :: halo !< The extra halo size to use here. logical, intent(in) :: use_BT_cont !< If true, use the BT_cont_types to calculate !! transports. @@ -3339,6 +3342,7 @@ subroutine set_up_BT_OBC(OBC, eta, SpV_avg, BT_OBC, BT_Domain, G, GV, US, CS, MS allocate(BT_OBC%vhbt(isdw:iedw,jsdw-1:jedw), source=0.0) allocate(BT_OBC%vbt_outer(isdw:iedw,jsdw-1:jedw), source=0.0) allocate(BT_OBC%SSH_outer_v(isdw:iedw,jsdw-1:jedw), source=0.0) + BT_OBC%is_alloced = .true. call create_group_pass(BT_OBC%pass_uv, BT_OBC%ubt_outer, BT_OBC%vbt_outer, BT_Domain) call create_group_pass(BT_OBC%pass_uhvh, BT_OBC%uhbt, BT_OBC%vhbt, BT_Domain) @@ -3481,6 +3485,7 @@ subroutine destroy_BT_OBC(BT_OBC) deallocate(BT_OBC%vhbt) deallocate(BT_OBC%vbt_outer) deallocate(BT_OBC%SSH_outer_v) + BT_OBC%is_alloced = .false. endif end subroutine destroy_BT_OBC @@ -4473,7 +4478,7 @@ end subroutine bt_mass_source !! barotropic calculation and initializes any barotropic fields that have not !! already been initialized. subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & - restart_CS, calc_dtbt, BT_cont, SAL_CSp, HA_CSp) + restart_CS, calc_dtbt, BT_cont, OBC, SAL_CSp, HA_CSp) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -4494,7 +4499,8 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & type(BT_cont_type), pointer :: BT_cont !< A structure with elements that describe the !! effective open face areas as a function of !! barotropic flow. - type(SAL_CS), target, optional :: SAL_CSp !< A pointer to the control structure of the + type(ocean_OBC_type), pointer :: OBC !< The open boundary condition structure. + type(SAL_CS), target, optional :: SAL_CSp !< A pointer to the control structure of the !! SAL module. type(harmonic_analysis_CS), target, optional :: HA_CSp !< A pointer to the control structure of the !! harmonic analysis module @@ -4692,6 +4698,10 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & "If true, recover a bug in barotropic solver that uses an unnormalized weight "//& "function for vertical averages of baroclinic velocity and forcing. Default "//& "of this flag is set by VISC_REM_BUG.", default=visc_rem_bug) + call get_param(param_file, mdl, "EXTERIOR_OBC_BUG", CS%exterior_OBC_bug, & + "If true, recover a bug in barotropic solver and other routines when "//& + "boundary contitions interior to the domain are used.", & + default=.true., do_not_log=.true.) call get_param(param_file, mdl, "TIDES", use_tides, & "If true, apply tidal momentum forcing.", default=.false.) if (use_tides .and. present(HA_CSp)) CS%HA_CSp => HA_CSp @@ -4934,11 +4944,49 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & ALLOC_(CS%IdyCv(CS%isdw:CS%iedw,CS%jsdw-1:CS%jedw)) ; CS%IdyCv(:,:) = 0.0 ALLOC_(CS%dy_Cu(CS%isdw-1:CS%iedw,CS%jsdw:CS%jedw)) ; CS%dy_Cu(:,:) = 0.0 ALLOC_(CS%dx_Cv(CS%isdw:CS%iedw,CS%jsdw-1:CS%jedw)) ; CS%dx_Cv(:,:) = 0.0 + allocate(CS%IareaT_OBCmask(isdw:iedw,jsdw:jedw), source=0.0) do j=G%jsd,G%jed ; do i=G%isd,G%ied CS%IareaT(i,j) = G%IareaT(i,j) CS%bathyT(i,j) = G%bathyT(i,j) + CS%IareaT_OBCmask(i,j) = CS%IareaT(i,j) enddo ; enddo + ! Update IareaT_OBCmask so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC) .and. (CS%exterior_OBC_bug .eqv. .false.)) then + if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then + do i=isd,ied + do j=jsd,jed + if (OBC%segnum_u(I-1,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I-1,j))%direction == OBC_DIRECTION_E) then + CS%IareaT_OBCmask(i,j) = 0.0 + endif + endif + if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + CS%IareaT_OBCmask(i,j) = 0.0 + endif + endif + enddo + enddo + endif + if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then + do j=jsd,jed + do i=isd,ied + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + CS%IareaT_OBCmask(i,j) = 0.0 + endif + endif + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + CS%IareaT_OBCmask(i,j) = 0.0 + endif + endif + enddo + enddo + endif + endif + ! Note: G%IdxCu & G%IdyCv may be valid for a smaller extent than CS%IdxCu & CS%IdyCv, even without ! wide halos. do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB @@ -4949,6 +4997,7 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & enddo ; enddo call create_group_pass(pass_static_data, CS%IareaT, CS%BT_domain, To_All) call create_group_pass(pass_static_data, CS%bathyT, CS%BT_domain, To_All) + call create_group_pass(pass_static_data, CS%IareaT_OBCmask, CS%BT_domain, To_All) call create_group_pass(pass_static_data, CS%IdxCu, CS%IdyCv, CS%BT_domain, To_All+Scalar_Pair) call create_group_pass(pass_static_data, CS%dy_Cu, CS%dx_Cv, CS%BT_domain, To_All+Scalar_Pair) call do_group_pass(pass_static_data, CS%BT_domain) @@ -5302,6 +5351,7 @@ subroutine barotropic_end(CS) if (allocated(CS%frhatu1)) deallocate(CS%frhatu1) if (allocated(CS%frhatv1)) deallocate(CS%frhatv1) + if (allocated(CS%IareaT_OBCmask)) deallocate(CS%IareaT_OBCmask) call deallocate_MOM_domain(CS%BT_domain) ! Allocated in restart registration, prior to timestep initialization diff --git a/src/core/MOM_boundary_update.F90 b/src/core/MOM_boundary_update.F90 index 75f69dc779..31863d10c2 100644 --- a/src/core/MOM_boundary_update.F90 +++ b/src/core/MOM_boundary_update.F90 @@ -145,9 +145,6 @@ subroutine update_OBC_data(OBC, G, GV, US, tv, h, CS, Time) type(update_OBC_CS), pointer :: CS !< Control structure for OBCs type(time_type), intent(in) :: Time !< Model time -! Something here... with CS%file_OBC_CSp? -! if (CS%use_files) & -! call update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (CS%use_tidal_bay) & call tidal_bay_set_OBC_data(OBC, CS%tidal_bay_OBC, G, GV, US, h, Time) if (CS%use_Kelvin) & diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 7169a4ee2c..5ab26f1097 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1554,7 +1554,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo call barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, & - CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp, HA_CSp) + CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, & + CS%OBC, CS%SAL_CSp, HA_CSp) if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. & .not. query_initialized(CS%diffv, "diffv", restart_CS)) then diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index 703fcd4ef7..28c8b6b35f 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -1471,7 +1471,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, call barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, & CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, & - CS%SAL_CSp, HA_CSp) + CS%OBC, CS%SAL_CSp, HA_CSp) flux_units = get_flux_units(GV) thickness_units = get_thickness_units(GV) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 33547590de..a91afdb662 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -393,6 +393,7 @@ module MOM_open_boundary logical :: om4_remap_via_sub_cells !< If true, use the OM4 remapping algorithm character(40) :: remappingScheme !< String selecting the vertical remapping scheme type(group_pass_type) :: pass_oblique !< Structure for group halo pass + logical :: exterior_OBC_bug !< If true, use incorrect form of tracers exterior to OBCs. end type ocean_OBC_type !> Control structure for open boundaries that read from files. @@ -561,6 +562,10 @@ subroutine open_boundary_config(G, US, param_file, OBC) "A silly value of velocities used outside of open boundary "//& "conditions for debugging.", units="m/s", default=0.0, scale=US%m_s_to_L_T, & do_not_log=.not.OBC%debug, debuggingParam=.true.) + call get_param(param_file, mdl, "EXTERIOR_OBC_BUG", OBC%exterior_OBC_bug, & + "If true, recover a bug in barotropic solver and other routines when "//& + "boundary contitions interior to the domain are used.", & + default=.true.) reentrant_x = .false. call get_param(param_file, mdl, "REENTRANT_X", reentrant_x, default=.true.) reentrant_y = .false. @@ -5061,7 +5066,9 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) integer :: i, j integer :: l_seg logical :: fatal_error = .False. - real :: min_depth ! The minimum depth for ocean points [Z ~> m] + real :: min_depth ! The minimum depth for ocean points [Z ~> m] + real :: mask_depth ! The masking depth for ocean points [Z ~> m] + real :: Dmask ! The depth for masking in the same units as G%bathyT [Z ~> m]. integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 character(len=256) :: mesg ! Message for error messages. real, allocatable, dimension(:,:) :: color, color2 ! For sorting inside from outside, @@ -5071,6 +5078,12 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & units="m", default=0.0, scale=US%m_to_Z, do_not_log=.true.) + call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, & + units="m", default=-9999.0, scale=US%m_to_Z, do_not_log=.true.) + + Dmask = mask_depth + if (mask_depth == -9999.0*US%m_to_Z) Dmask = min_depth + ! The reference depth on a dyn_horgrid is 0, otherwise would need: min_depth = min_depth - G%Z_ref allocate(color(G%isd:G%ied, G%jsd:G%jed), source=0.0) @@ -5161,7 +5174,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) &"the masking of the outside grid points.")') i, j call MOM_error(WARNING,"MOM mask_outside_OBCs: "//mesg, all_print=.true.) endif - if (color(i,j) == cout) G%bathyT(i,j) = min_depth + if (color(i,j) == cout) G%bathyT(i,j) = Dmask enddo ; enddo if (fatal_error) call MOM_error(FATAL, & "MOM_open_boundary: inconsistent OBC segments.") @@ -5463,8 +5476,8 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (G%mask2dT(I+ishift,j) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep do m=1,segment%tr_Reg%ntseg - ntr_id = segment%tr_reg%Tr(m)%ntr_index - fd_id = segment%tr_reg%Tr(m)%fd_index + ntr_id = segment%tr_Reg%Tr(m)%ntr_index + fd_id = segment%tr_Reg%Tr(m)%fd_index if (fd_id == -1) then resrv_lfac_out = 1.0 resrv_lfac_in = 1.0 @@ -5507,8 +5520,8 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (G%mask2dT(i,j+jshift) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep do m=1,segment%tr_Reg%ntseg - ntr_id = segment%tr_reg%Tr(m)%ntr_index - fd_id = segment%tr_reg%Tr(m)%fd_index + ntr_id = segment%tr_Reg%Tr(m)%ntr_index + fd_id = segment%tr_Reg%Tr(m)%fd_index if (fd_id == -1) then resrv_lfac_out = 1.0 resrv_lfac_in = 1.0 diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 0e129b2d03..6d20072987 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -648,6 +648,19 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & endif enddo + ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC) .and. (OBC%exterior_OBC_bug .eqv. .false.)) then ; if (OBC%OBC_pe) then + if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then + do i=is,ie-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + do_i(i+1,j) = .false. + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + do_i(i,j) = .false. + endif + endif ; enddo + endif + endif ; endif + ! update tracer concentration from i-flux and save some diagnostics do m=1,ntr @@ -1039,6 +1052,24 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & else ; do_i(i,j) = .false. ; endif enddo + ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC) .and. (OBC%exterior_OBC_bug .eqv. .false.)) then ; if (OBC%OBC_pe) then + if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then + do i=is,ie + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + do_i(i,j) = .false. + endif + endif + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + do_i(i,j) = .false. + endif + endif + enddo + endif + endif ; endif + ! update tracer and save some diagnostics do m=1,ntr do i=is,ie ; if (do_i(i,j)) then From 217172316dfa9a905ccfa326ded8654ecc24c33c Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 25 Feb 2025 12:42:40 -0900 Subject: [PATCH 2/2] Fix OBC check if not associated. --- src/tracer/MOM_tracer_advect.F90 | 50 +++++++++++++++++--------------- 1 file changed, 27 insertions(+), 23 deletions(-) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 6d20072987..dd54b55e9d 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -649,17 +649,19 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & enddo ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) - if (associated(OBC) .and. (OBC%exterior_OBC_bug .eqv. .false.)) then ; if (OBC%OBC_pe) then - if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then - do i=is,ie-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then - if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then - do_i(i+1,j) = .false. - elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then - do_i(i,j) = .false. - endif - endif ; enddo + if (associated(OBC)) then + if ((OBC%exterior_OBC_bug .eqv. .false.) .and. (OBC%OBC_pe)) then + if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then + do i=is,ie-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + do_i(i+1,j) = .false. + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + do_i(i,j) = .false. + endif + endif ; enddo + endif endif - endif ; endif + endif ! update tracer concentration from i-flux and save some diagnostics do m=1,ntr @@ -1053,22 +1055,24 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & enddo ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) - if (associated(OBC) .and. (OBC%exterior_OBC_bug .eqv. .false.)) then ; if (OBC%OBC_pe) then - if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then - do i=is,ie - if (OBC%segnum_v(i,J-1) /= OBC_NONE) then - if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then - do_i(i,j) = .false. + if (associated(OBC)) then + if ((OBC%exterior_OBC_bug .eqv. .false.) .and. (OBC%OBC_pe)) then + if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then + do i=is,ie + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + do_i(i,j) = .false. + endif endif - endif - if (OBC%segnum_v(i,J) /= OBC_NONE) then - if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then - do_i(i,j) = .false. + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + do_i(i,j) = .false. + endif endif - endif - enddo + enddo + endif endif - endif ; endif + endif ! update tracer and save some diagnostics do m=1,ntr