Skip to content

Commit

Permalink
Merge branch 'fix_obc_maskingdepth' into esmg_work
Browse files Browse the repository at this point in the history
  • Loading branch information
@kshedstrom
kshedstrom committed Nov 22, 2024
2 parents 3913e87 + 677baba commit 000ad9d
Show file tree
Hide file tree
Showing 4 changed files with 117 additions and 38 deletions.
20 changes: 20 additions & 0 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2883,6 +2883,26 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
! reservoirs are used.
call open_boundary_register_restarts(HI, GV, US, CS%OBC, CS%tracer_Reg, &
param_file, restart_CSp, use_temperature)
if (turns /= 0) then
if (CS%OBC%radiation_BCs_exist_globally) then
OBC_in%rx_normal => CS%OBC%rx_normal
OBC_in%ry_normal => CS%OBC%ry_normal
endif
if (CS%OBC%oblique_BCs_exist_globally) then
OBC_in%rx_oblique_u => CS%OBC%rx_oblique_u
OBC_in%ry_oblique_u => CS%OBC%ry_oblique_u
OBC_in%rx_oblique_v => CS%OBC%rx_oblique_v
OBC_in%ry_oblique_v => CS%OBC%ry_oblique_v
OBC_in%cff_normal_u => CS%OBC%cff_normal_u
OBC_in%cff_normal_v => CS%OBC%cff_normal_v
endif
if (any(CS%OBC%tracer_x_reservoirs_used)) then
OBC_in%tres_x => CS%OBC%tres_x
endif
if (any(CS%OBC%tracer_y_reservoirs_used)) then
OBC_in%tres_y => CS%OBC%tres_y
endif
endif
endif

if (present(waves_CSp)) then
Expand Down
40 changes: 38 additions & 2 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -742,6 +742,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
integer :: ioff, joff
integer :: l_seg
real :: factor(SZI_(G),SZJ_(G)) ! If non-zero, work on given points.

if (.not.CS%module_is_initialized) call MOM_error(FATAL, &
"btstep: Module MOM_barotropic must be initialized before it is used.")
Expand Down Expand Up @@ -2457,17 +2458,52 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
haloshift=iev-ie, unscale=US%L_to_m**2*GV%H_to_m)
endif

do j=jsv,jev
do i=isv,iev
factor(i,j) = CS%IareaT(i,j)
enddo
enddo

! Update factor so that nothing changes outside of the OBC (problem for interior OBCs only)
if (associated(OBC)) then ; if (OBC%OBC_pe) then
do j=jsv,jev
if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then
do i=isv,iev-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then
if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then
factor(i+1,j) = 0.0
elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then
factor(i,j) = 0.0
endif
endif ; enddo
endif
if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then
do i=isv,iev
if (OBC%segnum_v(i,J-1) /= OBC_NONE) then
if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then
factor(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
factor(i,j) = 0.0
endif
endif
enddo
endif
enddo
endif ; endif

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)) + factor(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 * factor(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
78 changes: 48 additions & 30 deletions src/core/MOM_open_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -352,24 +352,24 @@ module MOM_open_boundary
type(remapping_CS), pointer :: remap_h_CS=> NULL() !< ALE remapping control structure for
!! thickness-based fields on segments
type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries
real, allocatable :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of
real, pointer :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of
!! grid points per timestep [nondim]
real, allocatable :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of
real, pointer :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of
!! grid points per timestep [nondim]
real, allocatable :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared
real, pointer :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared
!! at u points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared
real, pointer :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared
!! at u points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared
real, pointer :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared
!! at v points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared
real, pointer :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared
!! at v points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation
real, pointer :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation
!! rates at u points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation
real, pointer :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation
!! rates at v points for restarts [L2 T-2 ~> m2 s-2]
real, allocatable :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc]
real, allocatable :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc]
real, pointer :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc]
real, pointer :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc]
logical :: debug !< If true, write verbose checksums for debugging purposes.
real :: silly_h !< A silly value of thickness outside of the domain that can be used to test
!! the independence of the OBCs to this external data [Z ~> m].
Expand Down Expand Up @@ -1948,15 +1948,15 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS)
call create_group_pass(OBC%pass_oblique, OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair)
call do_group_pass(OBC%pass_oblique, G%Domain)
endif
if (allocated(OBC%tres_x) .and. allocated(OBC%tres_y)) then
if (associated(OBC%tres_x) .and. associated(OBC%tres_y)) then
do m=1,OBC%ntr
call pass_vector(OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%Domain, To_All+Scalar_Pair)
enddo
elseif (allocated(OBC%tres_x)) then
elseif (associated(OBC%tres_x)) then
do m=1,OBC%ntr
call pass_var(OBC%tres_x(:,:,:,m), G%Domain, position=EAST_FACE)
enddo
elseif (allocated(OBC%tres_y)) then
elseif (associated(OBC%tres_y)) then
do m=1,OBC%ntr
call pass_var(OBC%tres_y(:,:,:,m), G%Domain, position=NORTH_FACE)
enddo
Expand Down Expand Up @@ -2001,16 +2001,16 @@ subroutine open_boundary_dealloc(OBC)
if (allocated(OBC%segment)) deallocate(OBC%segment)
if (allocated(OBC%segnum_u)) deallocate(OBC%segnum_u)
if (allocated(OBC%segnum_v)) deallocate(OBC%segnum_v)
if (allocated(OBC%rx_normal)) deallocate(OBC%rx_normal)
if (allocated(OBC%ry_normal)) deallocate(OBC%ry_normal)
if (allocated(OBC%rx_oblique_u)) deallocate(OBC%rx_oblique_u)
if (allocated(OBC%ry_oblique_u)) deallocate(OBC%ry_oblique_u)
if (allocated(OBC%rx_oblique_v)) deallocate(OBC%rx_oblique_v)
if (allocated(OBC%ry_oblique_v)) deallocate(OBC%ry_oblique_v)
if (allocated(OBC%cff_normal_u)) deallocate(OBC%cff_normal_u)
if (allocated(OBC%cff_normal_v)) deallocate(OBC%cff_normal_v)
if (allocated(OBC%tres_x)) deallocate(OBC%tres_x)
if (allocated(OBC%tres_y)) deallocate(OBC%tres_y)
if (associated(OBC%rx_normal)) nullify(OBC%rx_normal)
if (associated(OBC%ry_normal)) nullify(OBC%ry_normal)
if (associated(OBC%rx_oblique_u)) nullify(OBC%rx_oblique_u)
if (associated(OBC%ry_oblique_u)) nullify(OBC%ry_oblique_u)
if (associated(OBC%rx_oblique_v)) nullify(OBC%rx_oblique_v)
if (associated(OBC%ry_oblique_v)) nullify(OBC%ry_oblique_v)
if (associated(OBC%cff_normal_u)) nullify(OBC%cff_normal_u)
if (associated(OBC%cff_normal_v)) nullify(OBC%cff_normal_v)
if (associated(OBC%tres_x)) nullify(OBC%tres_x)
if (associated(OBC%tres_y)) nullify(OBC%tres_y)
if (associated(OBC%remap_z_CS)) deallocate(OBC%remap_z_CS)
if (associated(OBC%remap_h_CS)) deallocate(OBC%remap_h_CS)
deallocate(OBC)
Expand Down Expand Up @@ -3369,7 +3369,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US,
haloshift=0, symmetric=sym, unscale=1.0/US%L_T_to_m_s**2)
endif
if (OBC%ntr == 0) return
if (.not. allocated (OBC%tres_x) .or. .not. allocated (OBC%tres_y)) return
if (.not. associated (OBC%tres_x) .or. .not. associated (OBC%tres_y)) return
do m=1,OBC%ntr
write(var_num,'(I3.3)') m
call uvchksum("radiation_OBCs: OBC%tres_[xy]_"//var_num, OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%HI, &
Expand Down Expand Up @@ -5053,7 +5053,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 @@ -5063,6 +5065,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 @@ -5153,7 +5161,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 @@ -5481,7 +5489,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(I,j,k)+ &
((u_L_out+a_out)*Reg%Tr(ntr_id)%t(I+ishift,j,k) - &
(u_L_in+a_in)*segment%tr_Reg%Tr(m)%t(I,j,k)))
if (allocated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k)
if (associated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k)
enddo ; endif
enddo
enddo
Expand Down Expand Up @@ -5521,7 +5529,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(i,J,k) + &
((v_L_out+a_out)*Reg%Tr(ntr_id)%t(i,J+jshift,k) - &
(v_L_in+a_in)*segment%tr_Reg%Tr(m)%t(i,J,k)))
if (allocated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k)
if (associated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k)
enddo ; endif
enddo
enddo
Expand Down Expand Up @@ -5597,7 +5605,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell)

! Update tracer concentrations
segment%tr_Reg%Tr(m)%tres(I,j,:) = tr_column(:)
if (allocated(OBC%tres_x)) then ; do k=1,nz
if (associated(OBC%tres_x)) then ; do k=1,nz
OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k)
enddo ; endif

Expand Down Expand Up @@ -5664,7 +5672,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell)

! Update tracer concentrations
segment%tr_Reg%Tr(m)%tres(i,J,:) = tr_column(:)
if (allocated(OBC%tres_y)) then ; do k=1,nz
if (associated(OBC%tres_y)) then ; do k=1,nz
OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k)
enddo ; endif

Expand Down Expand Up @@ -6109,6 +6117,14 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns)
segment%field(n)%buffer_src)
endif

if (allocated(segment_in%field(n)%buffer_dst)) then
call allocate_rotated_array(segment_in%field(n)%buffer_dst, &
lbound(segment_in%field(n)%buffer_dst), turns, &
segment%field(n)%buffer_dst)
call rotate_array(segment_in%field(n)%buffer_dst, turns, &
segment%field(n)%buffer_dst)
endif

segment%field(n)%nk_src = segment_in%field(n)%nk_src

if (allocated(segment_in%field(n)%dz_src)) then
Expand All @@ -6122,6 +6138,8 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns)
segment%field(n)%value = segment_in%field(n)%value
enddo

call rotate_array(segment_in%SSH, turns, segment%SSH)

segment%temp_segment_data_exists = segment_in%temp_segment_data_exists
segment%salt_segment_data_exists = segment_in%salt_segment_data_exists
end subroutine rotate_OBC_segment_data
Expand Down
17 changes: 11 additions & 6 deletions src/tracer/MOM_tracer_advect.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1070,13 +1070,18 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only)
if (associated(OBC)) 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-1 ; if (OBC%segnum_v(i,J) /= OBC_NONE) then
if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then
do_i(i,j+1) = .false.
elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then
do_i(i,j) = .false.
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 ; enddo
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

Expand Down

0 comments on commit 000ad9d

Please sign in to comment.