Skip to content

Commit

Permalink
+(*)Fix sign of neutral_slope_x diag with no EOS
Browse files Browse the repository at this point in the history 
  Corrected the sign convention used for the neutral_slope_x and neutral_slope_y
diagnostics in cases when there is not an equation of state being used to match
the usual sign convention (interfaces sloping up to the east or north is a
positive slope) and to match what is done when an equation of state is being
used to calculate neutral density slopes rather than just basing the slopes of
the interfaces alone.  The same correction was made to the slope_x and slope_y
arguments returned from calc_isoneutral_slopes, and self-consistent changes were
made to the overturning streamfunction calculations in thickness_diffuse_full.

  In addition, there were some corrections made to the documentation at the end
of MOM_thickness_diffuse.F90, and the calculate_slopes argument to
calc_slope_functions_using_just_e() that was always hard-coded to .true. was
eliminated to avoid any confusion about whether these changes might propagate
beyond the interface height diffusion calculations.  This commit addresses most
aspects of the issue discussed at github.com/NOAA-GFDL/issues/359, although
it does not change the sign convention for the overturning streamfunction.

  All solutions are bitwise identical, but there is a change in the sign
convention for two diagnostics and two arguments to a publicly visible interface
in pure layer mode when no equation of state is used.
  • Loading branch information
@Hallberg-NOAA @c2xu
Hallberg-NOAA authored and c2xu committed Jun 11, 2024
1 parent 3b8ad83 commit a7b82cc
Show file tree
Hide file tree
Showing 3 changed files with 27 additions and 41 deletions.
12 changes: 6 additions & 6 deletions src/core/MOM_isopycnal_slopes.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -356,8 +356,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
slope = 0.0
endif
else ! With .not.use_EOS, the layers are constant density.
slope = (e(i,j,K)-e(i+1,j,K)) * G%IdxCu(I,j)
slope = (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)
endif

if (local_open_u_BC) then
l_seg = OBC%segnum_u(I,j)
if (l_seg /= OBC_NONE) then
Expand All @@ -372,7 +373,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
! endif
endif
endif
slope = slope * max(g%mask2dT(i,j),g%mask2dT(i+1,j))
slope = slope * max(G%mask2dT(i,j), G%mask2dT(i+1,j))
endif
slope_x(I,j,K) = slope
if (present(dzSxN)) &
Expand Down Expand Up @@ -504,11 +505,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
else ! Just in case mag_grad2 = 0 ever.
slope = 0.0
endif


else ! With .not.use_EOS, the layers are constant density.
slope = (e(i,j,K)-e(i,j+1,K)) * G%IdyCv(i,J)
slope = (e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)
endif

if (local_open_v_BC) then
l_seg = OBC%segnum_v(i,J)
if (l_seg /= OBC_NONE) then
Expand All @@ -523,7 +523,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
! endif
endif
endif
slope = slope * max(g%mask2dT(i,j),g%mask2dT(i,j+1))
slope = slope * max(G%mask2dT(i,j), G%mask2dT(i,j+1))
endif
slope_y(i,J,K) = slope
if (present(dzSyN)) &
Expand Down
40 changes: 13 additions & 27 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -497,8 +497,7 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, halo=1, OBC=OBC)
call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS)
else
!call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, CS%slope_x, CS%slope_y)
call calc_slope_functions_using_just_e(h, G, GV, US, CS, e, .true.)
call calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
endif
endif

Expand Down Expand Up @@ -821,16 +820,14 @@ end subroutine calc_Eady_growth_rate_2D

!> The original calc_slope_function() that calculated slopes using
!! interface positions only, not accounting for density variations.
subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slopes)
subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface position [Z ~> m]
! type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
logical, intent(in) :: calculate_slopes !< If true, calculate slopes
!! internally otherwise use slopes stored in CS
! Local variables
real :: E_x(SZIB_(G),SZJ_(G)) ! X-slope of interface at u points [Z L-1 ~> nondim] (for diagnostics)
real :: E_y(SZI_(G),SZJB_(G)) ! Y-slope of interface at v points [Z L-1 ~> nondim] (for diagnostics)
Expand Down Expand Up @@ -894,28 +891,17 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
!$OMP parallel do default(shared) private(E_x,E_y,S2,Hdn,Hup,H_geom,N2)
do k=nz,CS%VarMix_Ktop,-1

if (calculate_slopes) then
! Calculate the interface slopes E_x and E_y and u- and v- points respectively
do j=js-1,je+1 ; do I=is-1,ie
E_x(I,j) = (e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j)
! Mask slopes where interface intersects topography
if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) E_x(I,j) = 0.
enddo ; enddo
do J=js-1,je ; do i=is-1,ie+1
E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J)
! Mask slopes where interface intersects topography
if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(i,J) = 0.
enddo ; enddo
else ! This branch is not used.
do j=js-1,je+1 ; do I=is-1,ie
E_x(I,j) = CS%slope_x(I,j,k)
if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) E_x(I,j) = 0.
enddo ; enddo
do J=js-1,je ; do i=is-1,ie+1
E_y(i,J) = CS%slope_y(i,J,k)
if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(i,J) = 0.
enddo ; enddo
endif
! Calculate the interface slopes E_x and E_y and u- and v- points respectively
do j=js-1,je+1 ; do I=is-1,ie
E_x(I,j) = (e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j)
! Mask slopes where interface intersects topography
if (min(h(i,j,k),h(i+1,j,k)) < H_cutoff) E_x(I,j) = 0.
enddo ; enddo
do J=js-1,je ; do i=is-1,ie+1
E_y(i,J) = (e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J)
! Mask slopes where interface intersects topography
if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(i,J) = 0.
enddo ; enddo

! Calculate N*S*h from this layer and add to the sum
do j=js,je ; do I=is-1,ie
Expand Down
16 changes: 8 additions & 8 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1041,10 +1041,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
if (present_slope_x) then
Slope = slope_x(I,j,k)
else
Slope = ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) * G%OBCmaskCu(I,j)
Slope = ((e(i+1,j,K)-e(i,j,K))*G%IdxCu(I,j)) * G%OBCmaskCu(I,j)
endif
if (CS%id_slope_x > 0) CS%diagSlopeX(I,j,k) = Slope
Sfn_unlim_u(I,K) = ((KH_u(I,j,K)*G%dy_Cu(I,j))*Slope)
Sfn_unlim_u(I,K) = -(KH_u(I,j,K)*G%dy_Cu(I,j))*Slope
dzN2_u(I,K) = GV%g_prime(K)
endif ! if (use_EOS)
else ! if (k > nk_linear)
Expand Down Expand Up @@ -1361,10 +1361,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
if (present_slope_y) then
Slope = slope_y(i,J,k)
else
Slope = ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) * G%OBCmaskCv(i,J)
Slope = ((e(i,j+1,K)-e(i,j,K))*G%IdyCv(i,J)) * G%OBCmaskCv(i,J)
endif
if (CS%id_slope_y > 0) CS%diagSlopeY(I,j,k) = Slope
Sfn_unlim_v(i,K) = ((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope)
Sfn_unlim_v(i,K) = -((KH_v(i,J,K)*G%dx_Cv(i,J))*Slope)
dzN2_v(i,K) = GV%g_prime(K)
endif ! if (use_EOS)
else ! if (k > nk_linear)
Expand Down Expand Up @@ -2451,19 +2451,19 @@ end subroutine thickness_diffuse_end
!! to the potential density slope
!! \f[
!! \vec{\psi} = - \kappa_h \frac{\nabla_z \rho}{\partial_z \rho}
!! = \frac{g\kappa_h}{\rho_o} \frac{\nabla \rho}{N^2} = \kappa_h \frac{M^2}{N^2}
!! = \frac{g\kappa_h}{\rho_o} \frac{\nabla \rho}{N^2} = -\kappa_h \frac{M^2}{N^2}
!! \f]
!! but for robustness the scheme is implemented as
!! \f[
!! \vec{\psi} = \kappa_h \frac{M^2}{\sqrt{N^4 + M^4}}
!! \vec{\psi} = -\kappa_h \frac{M^2}{\sqrt{N^4 + M^4}}
!! \f]
!! since the quantity \f$\frac{M^2}{\sqrt{N^2 + M^2}}\f$ is bounded between $-1$ and $1$ and does not change sign
!! since the quantity \f$\frac{M^2}{\sqrt{N^4 + M^4}}\f$ is bounded between $-1$ and $1$ and does not change sign
!! if \f$N^2<0\f$.
!!
!! Optionally, the method of Ferrari et al, 2010, can be used to obtain the streamfunction which solves the
!! vertically elliptic equation:
!! \f[
!! \gamma_F \partial_z c^2 \partial_z \psi - N_*^2 \psi = ( 1 + \gamma_F ) \kappa_h N_*^2 \frac{M^2}{\sqrt{N^4+M^4}}
!! \gamma_F \partial_z c^2 \partial_z \psi - N_*^2 \psi = -( 1 + \gamma_F ) \kappa_h N_*^2 \frac{M^2}{\sqrt{N^4+M^4}}
!! \f]
!! which recovers the previous streamfunction relation in the limit that \f$ c \rightarrow 0 \f$.
!! Here, \f$c=\max(c_{min},c_g)\f$ is the maximum of either \f$c_{min}\f$ and either the first baroclinic mode
Expand Down

0 comments on commit a7b82cc

Please sign in to comment.