Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix interior obcs #814

Merged
merged 2 commits into from
Feb 26, 2025
Merged

Fix interior obcs #814

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
321f9c8
Adding a flag to turn off the bug fix (mostly)
kshedstrom Feb 25, 2025
2171723
Fix OBC check if not associated.
kshedstrom Feb 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 58 additions & 8 deletions src/core/MOM_barotropic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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)) - &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
3 changes: 0 additions & 3 deletions src/core/MOM_boundary_update.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) &
Expand Down
3 changes: 2 additions & 1 deletion src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_dynamics_split_RK2b.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
25 changes: 19 additions & 6 deletions src/core/MOM_open_boundary.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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,
Expand All @@ -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)
Expand Down Expand Up @@ -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.")
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
35 changes: 35 additions & 0 deletions src/tracer/MOM_tracer_advect.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -648,6 +648,21 @@ 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)) 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

! update tracer concentration from i-flux and save some diagnostics
do m=1,ntr

Expand Down Expand Up @@ -1039,6 +1054,26 @@ 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)) 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
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
Expand Down