Skip to content

Commit

Permalink
(*)Correct USE_STANLEY domain extents in EOS calls
Browse files Browse the repository at this point in the history 
  Corrected the domain extents used in EOS calls that are triggered with
USE_STANLEY_GM and USE_STANLEY_ISO.  These bugs were accidentally introduced
when the changes adding the MOM_stoch_eos code to the main branch of MOM6 were
merged with changes on the dev/gfdl branch.  Also added a test for cases when
USE_STANLEY_GM is set to true but STANLEY_COEF is negative to reset the internal
versions of this flag to false with a sensible warning message rather than
encountering segmentation faults.  All solutions are bitwise identical in cases
that worked before.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Dec 8, 2022
1 parent 47cc149 commit be923c7
Show file tree
Hide file tree
Showing 2 changed files with 71 additions and 59 deletions.
92 changes: 44 additions & 48 deletions src/core/MOM_isopycnal_slopes.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,12 @@ module MOM_isopycnal_slopes

! This file is part of MOM6. See LICENSE.md for the license.

use MOM_grid, only : ocean_grid_type
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
use MOM_EOS, only : calculate_density_derivs
use MOM_EOS, only : calculate_density_second_derivs
use MOM_debugging, only : hchksum, uvchksum
use MOM_grid, only : ocean_grid_type
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
use MOM_EOS, only : calculate_density_derivs, calculate_density_second_derivs, EOS_domain
use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE
use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S

Expand All @@ -28,13 +28,12 @@ module MOM_isopycnal_slopes
!> Calculate isopycnal slopes, and optionally return other stratification dependent functions such as N^2
!! and dz*S^2*g-prime used, or calculable from factors used, during the calculation.
subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stanley, &
slope_x, slope_y, N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC) !, eta_to_m)
slope_x, slope_y, N2_u, N2_v, dzu, dzv, dzSxN, dzSyN, halo, OBC)
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
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface heights [Z ~> m] or units
!! given by 1/eta_to_m)
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface heights [Z ~> m]
type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
!! thermodynamic variables
real, intent(in) :: dt_kappa_smooth !< A smoothing vertical diffusivity
Expand All @@ -61,15 +60,12 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
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.
! Local variables
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
T, & ! The temperature [C ~> degC], with the values in
! in massless layers filled vertically by diffusion.
S !, & ! The filled salinity [S ~> ppt], with the values in
S ! The filled salinity [S ~> ppt], with the values in
! in massless layers filled vertically by diffusion.
! Rho ! Density itself, when a nonlinear equation of state is not in use [R ~> kg m-3].
real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: &
pres ! The pressure at an interface [R L2 T-2 ~> Pa].
real, dimension(SZI_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() that will be ingored.
Expand All @@ -96,15 +92,17 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
T_hr, & ! Temperature on the interface at the h (+1) point [C ~> degC].
S_hr, & ! Salinity on the interface at the h (+1) point [S ~> ppt]
pres_hr ! Pressure on the interface at the h (+1) point [R L2 T-2 ~> Pa].
real :: drdiA, drdiB ! Along layer zonal- and meridional- potential density
real :: drdjA, drdjB ! gradients in the layers above (A) and below (B) the
! interface times the grid spacing [R ~> kg m-3].
real :: drdiA, drdiB ! Along layer zonal potential density gradients in the layers above (A)
! and below (B) the interface times the grid spacing [R ~> kg m-3].
real :: drdjA, drdjB ! Along layer meridional potential density gradients in the layers above (A)
! and below (B) the interface times the grid spacing [R ~> kg m-3].
real :: drdkL, drdkR ! Vertical density differences across an interface [R ~> kg m-3].
real :: hg2A, hg2B ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
real :: hg2L, hg2R ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
real :: haA, haB, haL, haR ! Arithmetic mean thicknesses [H ~> m or kg m-2].
real :: dzaL, dzaR ! Temporary thicknesses in eta units [Z ~> m].
real :: wtA, wtB, wtL, wtR ! Unscaled weights, with various units.
real :: wtA, wtB ! Unnormalized weights of the slopes above and below [H3 ~> m3 or kg3 m-6]
real :: wtL, wtR ! Unnormalized weights of the slopes to the left and right [H3 Z ~> m4 or kg3 m-5]
real :: drdx, drdy ! Zonal and meridional density gradients [R L-1 ~> kg m-4].
real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4].
real :: slope ! The slope of density surfaces, calculated in a way
Expand All @@ -117,33 +115,34 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
! in roundoff and can be neglected [Z ~> m].
logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
real :: G_Rho0 ! The gravitational acceleration divided by density [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
real :: Z_to_L ! A conversion factor between from units for e to the
! units for lateral distances [L Z-1 ~> 1]
real :: L_to_Z ! A conversion factor between from units for lateral distances
! to the units for e [Z L-1 ~> 1]
real :: H_to_Z ! A conversion factor from thickness units to the units of e [Z H-1 ~> 1 or m3 kg-1]

logical :: present_N2_u, present_N2_v
integer, dimension(2) :: EOSdom_u, EOSdom_v ! Domains for the equation of state calculations at u and v points
logical :: local_open_u_BC, local_open_v_BC ! True if u- or v-face OBCs exist anywhere in the global domain.
integer, dimension(2) :: EOSdom_u ! The shifted I-computational domain to use for equation of
! state calculations at u-points.
integer, dimension(2) :: EOSdom_v ! The shifted i-computational domain to use for equation of
! state calculations at v-points.
integer, dimension(2) :: EOSdom_h1 ! The shifted i-computational domain to use for equation of
! state calculations at h points with 1 extra halo point
integer :: is, ie, js, je, nz, IsdB
integer :: i, j, k
integer :: l_seg
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
EOSdom_h1(:) = EOS_domain(G%HI, halo=halo+1)
else
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
EOSdom_h1(:) = EOS_domain(G%HI, halo=1)
endif
EOSdom_u(1) = is-1 - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1)
EOSdom_v(:) = EOS_domain(G%HI, halo=halo)

nz = GV%ke ; IsdB = G%IsdB


h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2
Z_to_L = US%Z_to_L ; H_to_Z = GV%H_to_Z
! if (present(eta_to_m)) then
! Z_to_L = eta_to_m*US%m_to_L ; H_to_Z = GV%H_to_m / eta_to_m
! endif
L_to_Z = 1.0 / Z_to_L
dz_neglect = GV%H_subroundoff * H_to_Z
dz_neglect = GV%H_subroundoff * GV%H_to_Z

local_open_u_BC = .false.
local_open_v_BC = .false.
Expand Down Expand Up @@ -221,12 +220,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
enddo ; enddo
enddo

EOSdom_u(1) = is-1 - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1)

!$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,dzSxN,EOSdom_u,local_open_u_BC, &
!$OMP dzu,OBC,use_stanley) &
!$OMP h_neglect,dz_neglect,h_neglect2, &
!$OMP present_N2_u,G_Rho0,N2_u,slope_x,dzSxN,EOSdom_u,EOSdom_h1, &
!$OMP local_open_u_BC,dzu,OBC,use_stanley) &
!$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 drho_dT_dT_h,scrap,pres_h,T_h,S_h, &
Expand Down Expand Up @@ -259,7 +256,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
call calculate_density_second_derivs(T_h, S_h, pres_h, &
scrap, scrap, drho_dT_dT_h, scrap, scrap, &
tv%eqn_of_state, dom=[is-1,ie-is+3])
tv%eqn_of_state, dom=EOSdom_h1)
endif

do I=is-1,ie
Expand Down Expand Up @@ -294,7 +291,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
haR = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
if (GV%Boussinesq) then
dzaL = haL * H_to_Z ; dzaR = haR * H_to_Z
dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z
else
dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect
dzaR = 0.5*(e(i+1,j,K-1) - e(i+1,j,K+1)) + dz_neglect
Expand All @@ -318,7 +315,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan

! This estimate of slope is accurate for small slopes, but bounded
! to be between -1 and 1.
mag_grad2 = (Z_to_L*drdx)**2 + drdz**2
mag_grad2 = (US%Z_to_L*drdx)**2 + drdz**2
if (mag_grad2 > 0.0) then
slope = drdx / sqrt(mag_grad2)
else ! Just in case mag_grad2 = 0 ever.
Expand Down Expand Up @@ -351,11 +348,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
enddo ! I
enddo ; enddo ! end of j-loop

EOSdom_v(1) = is - (G%isd-1) ; EOSdom_v(2) = ie - (G%isd-1)

! 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,h_neglect,e,dz_neglect, &
!$OMP h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,dzSyN,EOSdom_v, &
!$OMP dzv,local_open_v_BC,OBC,use_stanley) &
!$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, &
Expand Down Expand Up @@ -393,10 +388,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
call calculate_density_second_derivs(T_h, S_h, pres_h, &
scrap, scrap, drho_dT_dT_h, scrap, scrap, &
tv%eqn_of_state, dom=[is,ie-is+1])
tv%eqn_of_state, dom=EOSdom_v)
call calculate_density_second_derivs(T_hr, S_hr, pres_hr, &
scrap, scrap, drho_dT_dT_hr, scrap, scrap, &
tv%eqn_of_state, dom=[is,ie-is+1])
tv%eqn_of_state, dom=EOSdom_v)
endif
do i=is,ie
if (use_EOS) then
Expand Down Expand Up @@ -430,7 +425,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan
haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
haR = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
if (GV%Boussinesq) then
dzaL = haL * H_to_Z ; dzaR = haR * H_to_Z
dzaL = haL * GV%H_to_Z ; dzaR = haR * GV%H_to_Z
else
dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect
dzaR = 0.5*(e(i,j+1,K-1) - e(i,j+1,K+1)) + dz_neglect
Expand All @@ -454,7 +449,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan

! This estimate of slope is accurate for small slopes, but bounded
! to be between -1 and 1.
mag_grad2 = (Z_to_L*drdy)**2 + drdz**2
mag_grad2 = (US%Z_to_L*drdy)**2 + drdz**2
if (mag_grad2 > 0.0) then
slope = drdy / sqrt(mag_grad2)
else ! Just in case mag_grad2 = 0 ever.
Expand Down Expand Up @@ -513,8 +508,9 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, lar
! Local variables
real :: ent(SZI_(G),SZK_(GV)+1) ! The diffusive entrainment (kappa*dt)/dz
! between layers in a timestep [H ~> m or kg m-2].
real :: b1(SZI_(G)), d1(SZI_(G)) ! b1, c1, and d1 are variables used by the
real :: c1(SZI_(G),SZK_(GV)) ! tridiagonal solver.
real :: b1(SZI_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]
real :: d1(SZI_(G)) ! A variable used by the tridiagonal solver [nondim], d1 = 1 - c1.
real :: c1(SZI_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim].
real :: kap_dt_x2 ! The 2*kappa_dt converted to H units [H2 ~> m2 or kg2 m-4].
real :: h_neglect ! A negligible thickness [H ~> m or kg m-2], to allow for zero thicknesses.
real :: h0 ! A negligible thickness to allow for zero thickness layers without
Expand All @@ -541,7 +537,7 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, lar
T_f(i,j,k) = T_in(i,j,k) ; S_f(i,j,k) = S_in(i,j,k)
enddo ; enddo ; enddo
else
!$OMP parallel do default(shared) private(ent,b1,d1,c1,h_tr)
!$OMP parallel do default(shared) private(ent,b1,d1,c1,h_tr)
do j=js,je
do i=is,ie
ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
Expand Down
Loading

0 comments on commit be923c7

Please sign in to comment.