Skip to content

Commit

Permalink
+*Non-Boussinesq revision of lateral_mixing_coeffs
Browse files Browse the repository at this point in the history 
  This commit revises lateral_mixing_coeffs to work in an appropriate mixture of
thickness and vertical extent variables to avoid any dependence on the
Boussinesq reference density in non-Boussinesq mode, while retaining the
previous answers in Boussinesq mode.

  This commit adds the new runtime parameter FULL_DEPTH_EADY_GROWTH_RATE to
indicate that the denominator of an Eady growth rate calculation should be based
on the full depth of the water column, rather than the nominal depth of the
bathymetry.  The new option is only the default for fully non-Boussinesq cases.

  A primordial horizontal indexing bug was corrected in the v-direction slope
calculation.  Because it only applies for very shallow bathymetry, does not
appear to impact any existing test cases and went undetected for at least 12
years, it was corrected directly rather than wrapping in another new runtime
flag.  However, this bug is being retained for now in a comment to help with
review and debugging if the answers should change unexpectedly in some yet-to-be
identified configuration.

  Two debugging checksums were added for the output variables calculated in
calc_resoln_function.

  The case of some indices was corrected to follow the MOM6 soft convention using
case to indicate the staggering position of variables.  The previously incorrect
units of one comment were also fixed.  There is a new logical element in the
VarMix_CS type.  To accommodate these changes there are three new internal
variables in calc_slope_functions_using_just_e.  A total of 9 GV%H_to_Z
conversion factors were eliminated with this commit.

  N2 is no longer calculated separately in calc_slope_functions_using_just_e,
but this code is left in a comment as it may be instructive. This commit
involved changing the units of one internal variable in calc_QG_Leith_viscosity
to use inverse thickness units (as its descriptive comment already indicated).
There are already known problems with calc_QG_Leith_viscosity as documented with
a fatal error; this will be addressed in a subsequent commit.

  All answers are bitwise identical in the existing MOM6-examples test suite,
but they will change when fully non-Boussinesq, and there is a new entry in some
MOM_parameter_doc files.
  • Loading branch information
@Hallberg-NOAA @adcroft
Hallberg-NOAA authored and adcroft committed Oct 5, 2023
1 parent 6d68459 commit 15673a4
Showing 1 changed file with 103 additions and 45 deletions.
148 changes: 103 additions & 45 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ module MOM_lateral_mixing_coeffs
use MOM_domains, only : create_group_pass, do_group_pass
use MOM_domains, only : group_pass_type, pass_var, pass_vector
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_interface_heights, only : find_eta
use MOM_interface_heights, only : find_eta, thickness_to_dz
use MOM_isopycnal_slopes, only : calc_isoneutral_slopes
use MOM_grid, only : ocean_grid_type
use MOM_unit_scaling, only : unit_scale_type
Expand Down Expand Up @@ -59,16 +59,21 @@ module MOM_lateral_mixing_coeffs
!! This parameter is set depending on other parameters.
logical :: calculate_depth_fns !< If true, calculate all the depth factors.
!! This parameter is set depending on other parameters.
logical :: calculate_Eady_growth_rate !< If true, calculate all the Eady growth rate.
logical :: calculate_Eady_growth_rate !< If true, calculate all the Eady growth rates.
!! This parameter is set depending on other parameters.
logical :: use_stanley_iso !< If true, use Stanley parameterization in MOM_isopycnal_slopes
logical :: use_simpler_Eady_growth_rate !< If true, use a simpler method to calculate the
!! Eady growth rate that avoids division by layer thickness.
!! This parameter is set depending on other parameters.
logical :: full_depth_Eady_growth_rate !< If true, calculate the Eady growth rate based on an
!! average that includes contributions from sea-level changes
!! in its denominator, rather than just the nominal depth of
!! the bathymetry. This only applies when using the model
!! interface heights as a proxy for isopycnal slopes.
real :: cropping_distance !< Distance from surface or bottom to filter out outcropped or
!! incropped interfaces for the Eady growth rate calc [Z ~> m]
real :: h_min_N2 !< The minimum vertical distance to use in the denominator of the
!! bouyancy frequency used in the slope calculation [Z ~> m]
!! bouyancy frequency used in the slope calculation [H ~> m or kg m-2]

real, allocatable :: SN_u(:,:) !< S*N at u-points [T-1 ~> s-1]
real, allocatable :: SN_v(:,:) !< S*N at v-points [T-1 ~> s-1]
Expand Down Expand Up @@ -449,6 +454,12 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
if (CS%id_Res_fn > 0) call post_data(CS%id_Res_fn, CS%Res_fn_h, CS%diag)
endif

if (CS%debug) then
call hchksum(CS%cg1, "calc_resoln_fn cg1", G%HI, haloshift=1, scale=US%L_T_to_m_s)
call uvchksum("Res_fn_[uv]", CS%Res_fn_u, CS%Res_fn_v, G%HI, haloshift=0, &
scale=1.0, scalar_pair=.true.)
endif

end subroutine calc_resoln_function

!> Calculates and stores functions of isopycnal slopes, e.g. Sx, Sy, S*N, mostly used in the Visbeck et al.
Expand Down Expand Up @@ -684,7 +695,7 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, h, e, dzu, dzv, dzSxN, dzSyN,
integer :: i, j, k, l_seg
logical :: crop

dz_neglect = GV%H_subroundoff * GV%H_to_Z
dz_neglect = GV%dZ_subroundoff
D_scale = CS%Eady_GR_D_scale
if (D_scale<=0.) D_scale = 64.*GV%max_depth ! 0 means use full depth so choose something big
r_crp_dist = 1. / max( dz_neglect, CS%cropping_distance )
Expand Down Expand Up @@ -818,12 +829,16 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
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)
real :: dz_tot(SZI_(G),SZJ_(G)) ! The total thickness of the water columns [Z ~> m]
! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m]
real :: H_cutoff ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
real :: dZ_cutoff ! A minimum water column depth for masking [H ~> m or kg m-2]
real :: h_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [H ~> m or kg m-2].
real :: S2 ! Interface slope squared [Z2 L-2 ~> nondim]
Expand All @@ -834,6 +849,8 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
! the buoyancy frequency squared at u-points [Z T-2 ~> m s-2]
real :: S2N2_v_local(SZI_(G),SZJB_(G),SZK_(GV)) ! The depth integral of the slope times
! the buoyancy frequency squared at v-points [Z T-2 ~> m s-2]
logical :: use_dztot ! If true, use the total water column thickness rather than the
! bathymetric depth for certain calculations.
integer :: is, ie, js, je, nz
integer :: i, j, k
integer :: l_seg
Expand All @@ -851,6 +868,25 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop

h_neglect = GV%H_subroundoff
H_cutoff = real(2*nz) * (GV%Angstrom_H + h_neglect)
dZ_cutoff = real(2*nz) * (GV%Angstrom_Z + GV%dz_subroundoff)

use_dztot = CS%full_depth_Eady_growth_rate ! .or. .not.(GV%Boussinesq or GV%semi_Boussinesq)

if (use_dztot) then
!$OMP parallel do default(shared)
do j=js-1,je+1 ; do i=is-1,ie+1
dz_tot(i,j) = e(i,j,1) - e(i,j,nz+1)
enddo ; enddo
! The following mathematically equivalent expression is more expensive but is less
! sensitive to roundoff for large Z_ref:
! call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
! do j=js-1,je+1
! do i=is-1,ie+1 ; dz_tot(i,j) = 0.0 ; enddo
! do k=1,nz ; do i=is-1,ie+1
! dz_tot(i,j) = dz_tot(i,j) + dz(i,j,k)
! enddo ; enddo
! enddo
endif

! To set the length scale based on the deformation radius, use wave_speed to
! calculate the first-mode gravity wave speed and then blend the equatorial
Expand All @@ -864,85 +900,98 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
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.
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.
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.
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
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.
if (min(h(i,j,k),h(i,j+1,k)) < H_cutoff) E_y(i,J) = 0.
enddo ; enddo
endif

! Calculate N*S*h from this layer and add to the sum
do j=js,je ; do I=is-1,ie
S2 = ( E_x(I,j)**2 + 0.25*( &
(E_y(I,j)**2+E_y(I+1,j-1)**2) + (E_y(I+1,j)**2+E_y(I,j-1)**2) ) )
(E_y(i,J)**2+E_y(i+1,J-1)**2) + (E_y(i+1,J)**2+E_y(i,J-1)**2) ) )
if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < H_cutoff) S2 = 0.0

Hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
Hup = 2.*h(i+1,j,k)*h(i+1,j,k-1) / (h(i+1,j,k) + h(i+1,j,k-1) + h_neglect)
H_geom = sqrt(Hdn*Hup)
N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
if (min(h(i,j,k-1), h(i+1,j,k-1), h(i,j,k), h(i+1,j,k)) < H_cutoff) &
S2 = 0.0
S2N2_u_local(I,j,k) = (H_geom * GV%H_to_Z) * S2 * N2
! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
S2N2_u_local(I,j,k) = (H_geom * S2) * (GV%g_prime(k) / max(Hdn, Hup, CS%h_min_N2) )
enddo ; enddo
do J=js-1,je ; do i=is,ie
S2 = ( E_y(i,J)**2 + 0.25*( &
(E_x(i,J)**2+E_x(i-1,J+1)**2) + (E_x(i,J+1)**2+E_x(i-1,J)**2) ) )
(E_x(I,j)**2+E_x(I-1,j+1)**2) + (E_x(I,j+1)**2+E_x(I-1,j)**2) ) )
if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < H_cutoff) S2 = 0.0

Hdn = 2.*h(i,j,k)*h(i,j,k-1) / (h(i,j,k) + h(i,j,k-1) + h_neglect)
Hup = 2.*h(i,j+1,k)*h(i,j+1,k-1) / (h(i,j+1,k) + h(i,j+1,k-1) + h_neglect)
H_geom = sqrt(Hdn*Hup)
N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
if (min(h(i,j,k-1), h(i,j+1,k-1), h(i,j,k), h(i,j+1,k)) < H_cutoff) &
S2 = 0.0
S2N2_v_local(i,J,k) = (H_geom * GV%H_to_Z) * S2 * N2
! N2 = GV%g_prime(k) / (GV%H_to_Z * max(Hdn, Hup, CS%h_min_N2))
S2N2_v_local(i,J,k) = (H_geom * S2) * (GV%g_prime(k) / (max(Hdn, Hup, CS%h_min_N2)))
enddo ; enddo

enddo ! k

!$OMP parallel do default(shared)
do j=js,je
do I=is-1,ie ; CS%SN_u(I,j) = 0.0 ; enddo
do k=nz,CS%VarMix_Ktop,-1 ; do I=is-1,ie
CS%SN_u(I,j) = CS%SN_u(I,j) + S2N2_u_local(I,j,k)
enddo ; enddo
! SN above contains S^2*N^2*H, convert to vertical average of S*N
do I=is-1,ie
!### Replace G%bathT+G%Z_ref here with (e(i,j,1) - e(i,j,nz+1)).
!SN_u(I,j) = sqrt( SN_u(I,j) / ( max(G%bathyT(i,j), G%bathyT(i+1,j)) + (G%Z_ref + GV%Angstrom_Z) ) )
!The code below behaves better than the line above. Not sure why? AJA
if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > H_cutoff*GV%H_to_Z ) then

if (use_dztot) then
do I=is-1,ie
CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / &
(max(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref) )
else
CS%SN_u(I,j) = 0.0
endif
enddo
max(dz_tot(i,j), dz_tot(i+1,j), GV%dz_subroundoff) )
enddo
else
do I=is-1,ie
if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then
CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / &
(max(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref) )
else
CS%SN_u(I,j) = 0.0
endif
enddo
endif
enddo
!$OMP parallel do default(shared)
do J=js-1,je
do i=is,ie ; CS%SN_v(i,J) = 0.0 ; enddo
do k=nz,CS%VarMix_Ktop,-1 ; do i=is,ie
CS%SN_v(i,J) = CS%SN_v(i,J) + S2N2_v_local(i,J,k)
enddo ; enddo
do i=is,ie
!### Replace G%bathT+G%Z_ref here with (e(i,j,1) - e(i,j,nz+1)).
!SN_v(i,J) = sqrt( SN_v(i,J) / ( max(G%bathyT(i,J), G%bathyT(i,J+1)) + (G%Z_ref + GV%Angstrom_Z) ) )
!The code below behaves better than the line above. Not sure why? AJA
if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > H_cutoff*GV%H_to_Z ) then
if (use_dztot) then
do i=is,ie
CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / &
(max(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref) )
else
CS%SN_v(i,J) = 0.0
endif
enddo
max(dz_tot(i,j), dz_tot(i,j+1), GV%dz_subroundoff) )
enddo
else
do i=is,ie
! There is a primordial horizontal indexing bug on the following line from the previous
! versions of the code. This comment should be deleted by the end of 2024.
! if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then
if ( min(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref > dZ_cutoff ) then
CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / &
(max(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref) )
else
CS%SN_v(i,J) = 0.0
endif
enddo
endif
enddo

end subroutine calc_slope_functions_using_just_e
Expand Down Expand Up @@ -982,7 +1031,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo
real :: Ih ! The inverse of a combination of thicknesses [H-1 ~> m-1 or m2 kg-1]
real :: f ! A copy of the Coriolis parameter [T-1 ~> s-1]
real :: inv_PI3 ! The inverse of pi cubed [nondim]
integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq,nz
integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
Expand All @@ -1002,8 +1051,8 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo
h_at_slope_below = 2. * ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) * h(i+1,j,k+1) ) / &
( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) &
+ ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff**2 )
Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_Z )
dslopex_dz(I,j) = 2. * ( CS%slope_x(i,j,k) - CS%slope_x(i,j,k+1) ) * Ih
Ih = 1. / ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff )
dslopex_dz(I,j) = 2. * ( CS%slope_x(i,j,k) - CS%slope_x(i,j,k+1) ) * (GV%Z_to_H * Ih)
h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
enddo ; enddo

Expand All @@ -1016,8 +1065,8 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo
h_at_slope_below = 2. * ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) * h(i,j+1,k+1) ) / &
( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) &
+ ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff**2 )
Ih = 1./ ( ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) * GV%H_to_Z )
dslopey_dz(i,J) = 2. * ( CS%slope_y(i,j,k) - CS%slope_y(i,j,k+1) ) * Ih
Ih = 1. / ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff )
dslopey_dz(i,J) = 2. * ( CS%slope_y(i,j,k) - CS%slope_y(i,j,k+1) ) * (GV%Z_to_H * Ih)
h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih
enddo ; enddo

Expand Down Expand Up @@ -1143,7 +1192,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
CS%calculate_cg1 = .false.
CS%calculate_Rd_dx = .false.
CS%calculate_res_fns = .false.
CS%use_simpler_Eady_growth_rate = .false.
CS%use_simpler_Eady_growth_rate = .false.
CS%full_depth_Eady_growth_rate = .false.
CS%calculate_depth_fns = .false.
! Read all relevant parameters and write them to the model log.
call log_version(param_file, mdl, version, "")
Expand Down Expand Up @@ -1298,6 +1348,14 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
"The minimum vertical distance to use in the denominator of the "//&
"bouyancy frequency used in the slope calculation.", &
units="m", default=1.0, scale=GV%m_to_H, do_not_log=CS%use_stored_slopes)

call get_param(param_file, mdl, "FULL_DEPTH_EADY_GROWTH_RATE", CS%full_depth_Eady_growth_rate, &
"If true, calculate the Eady growth rate based on average slope times "//&
"stratification that includes contributions from sea-level changes "//&
"in its denominator, rather than just the nominal depth of the bathymetry. "//&
"This only applies when using the model interface heights as a proxy for "//&
"isopycnal slopes.", default=.not.(GV%Boussinesq.or.GV%semi_Boussinesq), &
do_not_log=CS%use_stored_slopes)
endif
endif

Expand Down

0 comments on commit 15673a4

Please sign in to comment.