Skip to content

Commit ae6b26f

Browse files
Hallberg-NOAAHallberg-NOAA
committed
+Add calc_isoneutral_slopes optional argument OBC_N2
 Added the new optional argument OBC_N2 to calc_isoneutral_slopes that can be use to specify that only interior points are used to calculate the buoyancy frequencies at velocity points when open boundary conditions are in use. When the new OBC_N2 argument is absent or false, the solutions are identical to what they had been previously, but this gives a dependence on the arbitrary outer thicknesses when open boundary conditions are used. By default, all answers are bitwise identical, but there is a new optional argument to a publicly visible interface.
1 parent cbb6172 commit ae6b26f

File tree

1 file changed

+76
-8
lines changed

1 file changed

+76
-8
lines changed

‎src/core/MOM_isopycnal_slopes.F90‎

Lines changed: 76 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -28,8 +28,8 @@ module MOM_isopycnal_slopes
2828

2929
!> Calculate isopycnal slopes, and optionally return other stratification dependent functions such as N^2
3030
!! and dz*S^2*g-prime used, or calculable from factors used, during the calculation.
31-
subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, &
32-
slope_x, slope_y, N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC)
31+
subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, slope_x, slope_y, &
32+
N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC, OBC_N2)
3333
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
3434
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
3535
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
@@ -61,6 +61,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
6161
!! Eady growth rate at v-points. [Z T-1 ~> m s-1]
6262
integer, optional, intent(in) :: halo !< Halo width over which to compute
6363
type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
64+
logical, optional, intent(in) :: OBC_N2 !< If present and true, use interior data
65+
!! to calculate stratification at open boundary
66+
!! condition faces.
6467

6568
! Local variables
6669
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
@@ -127,6 +130,8 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
127130

128131
logical :: present_N2_u, present_N2_v
129132
logical :: local_open_u_BC, local_open_v_BC ! True if u- or v-face OBCs exist anywhere in the global domain.
133+
logical :: OBC_friendly ! If true, open boundary conditions are in use and only interior data should
134+
! be used to calculate N2 at OBC faces.
130135
integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of
131136
! state calculations at u-points.
132137
integer, dimension(2) :: EOSdom_v ! The shifted i-computational domain to use for equation of
@@ -155,9 +160,11 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
155160

156161
local_open_u_BC = .false.
157162
local_open_v_BC = .false.
163+
OBC_friendly = .false.
158164
if (present(OBC)) then ; if (associated(OBC)) then
159165
local_open_u_BC = OBC%open_u_BCs_exist_globally
160166
local_open_v_BC = OBC%open_v_BCs_exist_globally
167+
if (present(OBC_N2)) OBC_friendly = OBC_N2
161168
endif ; endif
162169

163170
use_EOS = associated(tv%eqn_of_state)
@@ -241,12 +248,12 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
241248
enddo
242249

243250
do I=is-1,ie
244-
GxSpV_u(I) = G_Rho0 !This will be changed if both use_EOS and allocated(tv%SpV_avg) are true
251+
GxSpV_u(I) = G_Rho0 ! This will be changed if both use_EOS and allocated(tv%SpV_avg) are true
245252
enddo
246253
!$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, &
247254
!$OMP h_neglect,dz_neglect,h_neglect2, &
248255
!$OMP present_N2_u,G_Rho0,N2_u,slope_x,dzSxN,EOSdom_u,EOSdom_h1, &
249-
!$OMP local_open_u_BC,dzu,OBC,use_stanley) &
256+
!$OMP local_open_u_BC,dzu,OBC,use_stanley,OBC_friendly) &
250257
!$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, &
251258
!$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
252259
!$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, &
@@ -266,6 +273,22 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
266273
T_u(I) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1)))
267274
S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1)))
268275
enddo
276+
if (OBC_friendly) then
277+
do I=is-1,ie
278+
l_seg = OBC%segnum_u(I,j)
279+
if (l_seg /= OBC_NONE) then
280+
if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then
281+
pres_u(I) = pres(i,j,K)
282+
T_u(I) = 0.5*(T(i,j,k) + T(i,j,k-1))
283+
S_u(I) = 0.5*(S(i,j,k) + S(i,j,k-1))
284+
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then
285+
pres_u(I) = pres(i+1,j,K)
286+
T_u(I) = 0.5*(T(i+1,j,k) + T(i+1,j,k-1))
287+
S_u(I) = 0.5*(S(i+1,j,k) + S(i+1,j,k-1))
288+
endif
289+
endif
290+
enddo
291+
endif
269292
call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, &
270293
tv%eqn_of_state, EOSdom_u)
271294
if (present_N2_u .or. (present(dzSxN))) then
@@ -338,8 +361,21 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
338361
! The expression for drdz above is mathematically equivalent to:
339362
! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
340363
! ((hg2L/haL) + (hg2R/haR))
341-
! This is the gradient of density along geopotentials.
364+
! which is an estimate of the gradient of density across geopotentials.
342365
if (present_N2_u) then
366+
if (OBC_friendly) then ; if (OBC%segnum_u(I,j) /= OBC_NONE) then
367+
l_seg = OBC%segnum_u(I,j)
368+
if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then
369+
drdz = drdkL / dzaL ! Note that drdz is not used for slopes at OBC faces.
370+
if (use_EOS .and. allocated(tv%SpV_avg)) &
371+
GxSpV_u(I) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
372+
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then
373+
drdz = drdkR / dzaR
374+
if (use_EOS .and. allocated(tv%SpV_avg)) &
375+
GxSpV_u(I) = GV%g_Earth * 0.5 * (tv%SpV_avg(i+1,j,k) + tv%SpV_avg(i+1,j,k-1))
376+
endif
377+
endif ; endif
378+
343379
N2_u(I,j,K) = GxSpV_u(I) * drdz * G%mask2dCu(I,j) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
344380
endif
345381

@@ -375,6 +411,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
375411
endif
376412
slope = slope * max(G%mask2dT(i,j), G%mask2dT(i+1,j))
377413
endif
414+
378415
slope_x(I,j,K) = slope
379416
if (present(dzSxN)) &
380417
dzSxN(I,j,K) = sqrt( GxSpV_u(I) * max(0., (wtL * ( dzaL * drdkL )) &
@@ -391,7 +428,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
391428
!$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, &
392429
!$OMP h,h_neglect,e,dz_neglect, &
393430
!$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,dzSyN,EOSdom_v, &
394-
!$OMP dzv,local_open_v_BC,OBC,use_stanley) &
431+
!$OMP dzv,local_open_v_BC,OBC,use_stanley,OBC_friendly) &
395432
!$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, &
396433
!$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
397434
!$OMP drho_dT_dT_h,scrap,pres_h,T_h,S_h, &
@@ -411,6 +448,22 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
411448
T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1)))
412449
S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1)))
413450
enddo
451+
if (OBC_friendly) then
452+
do i=is,ie
453+
l_seg = OBC%segnum_v(i,J)
454+
if (l_seg /= OBC_NONE) then
455+
if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then
456+
pres_v(i) = pres(i,j,K)
457+
T_v(i) = 0.5*(T(i,j,k) + T(i,j,k-1))
458+
S_v(i) = 0.5*(S(i,j,k) + S(i,j,k-1))
459+
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then
460+
pres_v(i) = pres(i,j+1,K)
461+
T_v(i) = 0.5*(T(i,j+1,k) + T(i,j+1,k-1))
462+
S_v(i) = 0.5*(S(i,j+1,k) + S(i,j+1,k-1))
463+
endif
464+
endif
465+
enddo
466+
endif
414467
call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, &
415468
tv%eqn_of_state, EOSdom_v)
416469

@@ -490,8 +543,23 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
490543
! The expression for drdz above is mathematically equivalent to:
491544
! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
492545
! ((hg2L/haL) + (hg2R/haR))
493-
! This is the gradient of density along geopotentials.
494-
if (present_N2_v) N2_v(i,J,K) = GxSpV_v(i) * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
546+
! which is an estimate of the gradient of density across geopotentials.
547+
if (present_N2_v) then
548+
if (OBC_friendly) then ; if (OBC%segnum_v(i,J) /= OBC_NONE) then
549+
l_seg = OBC%segnum_v(i,J)
550+
if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then
551+
drdz = drdkL / dzaL ! Note that drdz is not used for slopes at OBC faces.
552+
if (use_EOS .and. allocated(tv%SpV_avg)) &
553+
GxSpV_v(i) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j,k-1))
554+
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then
555+
drdz = drdkL / dzaL
556+
if (use_EOS .and. allocated(tv%SpV_avg)) &
557+
GxSpV_v(i) = GV%g_Earth * 0.5 * (tv%SpV_avg(i,j+1,k) + tv%SpV_avg(i,j+1,k-1))
558+
endif
559+
endif ; endif
560+
561+
N2_v(i,J,K) = GxSpV_v(i) * drdz * G%mask2dCv(i,J) ! Square of buoyancy freq. [L2 Z-2 T-2 ~> s-2]
562+
endif
495563

496564
if (use_EOS) then
497565
drdy = ((wtA * drdjA + wtB * drdjB) / (wtA + wtB) - &

0 commit comments

Comments
 (0)