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.

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

Adding OBC code to slope_x, slope_y #1153

Merged
merged 8 commits into from
Jul 6, 2020
8 changes: 4 additions & 4 deletions src/core/MOM.F90
Original file line number Diff line number Diff line change
Expand Up @@ -994,7 +994,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
call enable_averages(dt_thermo, Time_local+real_to_time(US%T_to_s*(dt_thermo-dt)), CS%diag)
call cpu_clock_begin(id_clock_thick_diff)
if (associated(CS%VarMix)) &
call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix)
call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC)
call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, &
CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp)
call cpu_clock_end(id_clock_thick_diff)
Expand Down Expand Up @@ -1067,7 +1067,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
if (CS%debug) call hchksum(h,"Pre-thickness_diffuse h", G%HI, haloshift=0, scale=GV%H_to_m)

if (associated(CS%VarMix)) &
call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix)
call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC)
call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt, G, GV, US, &
CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp)

Expand Down Expand Up @@ -1479,7 +1479,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
call pass_var(CS%h, G%Domain)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
call calc_depth_function(G, CS%VarMix)
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix)
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
endif
call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, &
CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv)
Expand All @@ -1505,7 +1505,7 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS
call pass_var(CS%h, G%Domain)
call calc_resoln_function(CS%h, CS%tv, G, GV, US, CS%VarMix)
call calc_depth_function(G, CS%VarMix)
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix)
call calc_slope_functions(CS%h, CS%tv, dt_offline, G, GV, US, CS%VarMix, OBC=CS%OBC)
endif
call tracer_hordiff(CS%h, dt_offline, CS%MEKE, CS%VarMix, G, GV, US, &
CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv)
Expand Down
44 changes: 41 additions & 3 deletions src/core/MOM_isopycnal_slopes.F90
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ module MOM_isopycnal_slopes
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
use MOM_EOS, only : int_specific_vol_dp, calculate_density_derivs
use MOM_open_boundary, only : ocean_OBC_type
use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S

implicit none ; private

Expand All @@ -24,7 +26,7 @@ module MOM_isopycnal_slopes

!> Calculate isopycnal slopes, and optionally return N2 used in calculation.
subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
slope_x, slope_y, N2_u, N2_v, halo) !, eta_to_m)
slope_x, slope_y, N2_u, N2_v, halo, OBC) !, eta_to_m)
type(ocean_grid_type), intent(in) :: 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 @@ -44,6 +46,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
optional, intent(inout) :: N2_v !< Brunt-Vaisala frequency squared at
!! interfaces between u-points [T-2 ~> s-2]
integer, optional, intent(in) :: halo !< Halo width over which to compute
type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.

! real, optional, intent(in) :: eta_to_m !< The conversion factor from the units
! (This argument has been tested but for now serves no purpose.) !! of eta to m; US%Z_to_m by default.
Expand Down Expand Up @@ -102,6 +105,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
integer, dimension(2) :: EOSdom_u, EOSdom_v ! Domains for the equation of state calculations at u and v points
integer :: is, ie, js, je, nz, IsdB
integer :: i, j, k
logical :: local_open_u_BC, local_open_v_BC

if (present(halo)) then
is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo
Expand All @@ -118,6 +122,13 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
L_to_Z = 1.0 / Z_to_L
dz_neglect = GV%H_subroundoff * H_to_Z

local_open_u_BC = .false.
local_open_v_BC = .false.
if (present(OBC)) then ; if (associated(OBC)) then
local_open_u_BC = OBC%open_u_BCs_exist_globally
local_open_v_BC = OBC%open_v_BCs_exist_globally
endif ; endif

use_EOS = associated(tv%eqn_of_state)

present_N2_u = PRESENT(N2_u)
Expand Down Expand Up @@ -167,7 +178,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &

!$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, &
!$OMP h_neglect,dz_neglect,Z_to_L,L_to_Z,H_to_Z,h_neglect2, &
!$OMP present_N2_u,G_Rho0,N2_u,slope_x,EOSdom_u) &
!$OMP present_N2_u,G_Rho0,N2_u,slope_x,EOSdom_u,local_open_u_BC) &
!$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, &
!$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
!$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
Expand Down Expand Up @@ -247,6 +258,19 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
else ! With .not.use_EOS, the layers are constant density.
slope_x(I,j,K) = (Z_to_L*(e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j)
endif
if (local_open_u_BC) then
if (OBC%segment(OBC%segnum_u(I,j))%open) then
slope_x(I,j,K) = 0.
! This and/or the masking code below is to make slopes match inside
! land mask. Might not be necessary except for DEBUG output.
! if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then
! slope_x(I+1,j,K) = 0.
! else
! slope_x(I-1,j,K) = 0.
! endif
endif
slope_x(I,j,K) = slope_x(I,j,k) * max(g%mask2dT(i,j),g%mask2dT(i+1,j))
endif

enddo ! I
enddo ; enddo ! end of j-loop
Expand All @@ -256,7 +280,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
! Calculate the meridional isopycnal slope.
!$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, &
!$OMP h,h_neglect,e,dz_neglect,Z_to_L,L_to_Z,H_to_Z, &
!$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,EOSdom_v) &
!$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,EOSdom_v, &
!$OMP local_open_v_BC,OBC%segnum_u,OBC%segnum_v) &
!$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, &
!$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
!$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
Expand Down Expand Up @@ -333,6 +358,19 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
else ! With .not.use_EOS, the layers are constant density.
slope_y(i,J,K) = (Z_to_L*(e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J)
endif
if (local_open_v_BC) then
if (OBC%segment(OBC%segnum_v(i,J))%open) then
slope_y(i,J,K) = 0.
! This and/or the masking code below is to make slopes match inside
! land mask. Might not be necessary except for DEBUG output.
! if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then
! slope_y(i,J+1,K) = 0.
! else
! slope_y(i,J-1,K) = 0.
! endif
endif
slope_y(i,J,K) = slope_y(i,J,k) * max(g%mask2dT(i,j),g%mask2dT(i,j+1))
endif

enddo ! i
enddo ; enddo ! end of j-loop
Expand Down
7 changes: 6 additions & 1 deletion src/core/MOM_open_boundary.F90
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ module MOM_open_boundary
use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
use MOM_diag_mediator, only : diag_ctrl, time_type
use MOM_domains, only : pass_var, pass_vector
use MOM_domains, only : To_All, SCALAR_PAIR, CGRID_NE
use MOM_domains, only : To_All, SCALAR_PAIR, CGRID_NE, CORNER
use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe
use MOM_error_handler, only : NOTE
use MOM_file_parser, only : get_param, log_version, param_file_type, log_param
Expand Down Expand Up @@ -1598,6 +1598,11 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CSp)
if (.not.associated(OBC)) return

id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE)
if (OBC%radiation_BCs_exist_globally) call pass_vector(OBC%rx_normal, OBC%ry_normal, G%Domain, &
To_All+Scalar_Pair)
if (OBC%oblique_BCs_exist_globally) call pass_vector(OBC%rx_oblique, OBC%ry_oblique, G%Domain, &
To_All+Scalar_Pair)
if (associated(OBC%cff_normal)) call pass_var(OBC%cff_normal, G%Domain, position=CORNER)

! The rx_normal and ry_normal arrays used with radiation OBCs are currently in units of grid
! points per timestep, but if this were to be corrected to [L T-1 ~> m s-1] or [T-1 ~> s-1] to
Expand Down
Loading