Skip to content

Commit

Permalink
+wave_speed arg mono_N2_depth in thickness units
Browse files Browse the repository at this point in the history 
  Changed the units of the optional mono_N2_depth argument to wave_speed,
wave_speed_init and wave_speed_set_param in thickness units instead of height
units.  Accordingly, the units of one element each in the diagnostics_CS and
wave_speed_CS and a local variable in VarMix_init are also changed to thickness
units.  The unit descriptions of some comments describing diagnostics were also
amended to also describe the non-Boussinesq versions.  Because this is
essentially just changing when the unit conversion occurs, all answers are
bitwise identical, but there are changes to the units of an optional argument in
3 publicly visible routines.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed Jul 26, 2023
1 parent 636d610 commit 878fd1e
Show file tree
Hide file tree
Showing 3 changed files with 49 additions and 43 deletions.
24 changes: 12 additions & 12 deletions src/diagnostics/MOM_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ module MOM_diagnostics
!! monotonic for the purposes of calculating the equivalent
!! barotropic wave speed [nondim].
real :: mono_N2_depth = -1. !< The depth below which N2 is limited as monotonic for the purposes of
!! calculating the equivalent barotropic wave speed [Z ~> m].
!! calculating the equivalent barotropic wave speed [H ~> m or kg m-2].

type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
!! regulate the timing of diagnostic output.
Expand Down Expand Up @@ -984,7 +984,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_dKEdt > 0) then
! Calculate the time derivative of the layer KE [H L2 T-3 ~> m3 s-3].
! Calculate the time derivative of the layer KE [H L2 T-3 ~> m3 s-3 or W m-2].
do k=1,nz
do j=js,je ; do I=Isq,Ieq
KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * CS%du_dt(I,j,k)
Expand All @@ -1006,7 +1006,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_PE_to_KE > 0) then
! Calculate the potential energy to KE term [H L2 T-3 ~> m3 s-3].
! Calculate the potential energy to KE term [H L2 T-3 ~> m3 s-3 or W m-2].
do k=1,nz
do j=js,je ; do I=Isq,Ieq
KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%PFu(I,j,k)
Expand All @@ -1025,7 +1025,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_KE_BT > 0) then
! Calculate the barotropic contribution to KE term [H L2 T-3 ~> m3 s-3].
! Calculate the barotropic contribution to KE term [H L2 T-3 ~> m3 s-3 or W m-2].
do k=1,nz
do j=js,je ; do I=Isq,Ieq
KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%u_accel_bt(I,j,k)
Expand All @@ -1044,7 +1044,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_KE_Coradv > 0) then
! Calculate the KE source from the combined Coriolis and advection terms [H L2 T-3 ~> m3 s-3].
! Calculate the KE source from the combined Coriolis and advection terms [H L2 T-3 ~> m3 s-3 or W m-2].
! The Coriolis source should be zero, but is not due to truncation errors. There should be
! near-cancellation of the global integral of this spurious Coriolis source.
do k=1,nz
Expand All @@ -1069,7 +1069,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_KE_adv > 0) then
! Calculate the KE source from along-layer advection [H L2 T-3 ~> m3 s-3].
! Calculate the KE source from along-layer advection [H L2 T-3 ~> m3 s-3 or W m-2].
! NOTE: All terms in KE_adv are multiplied by -1, which can easily produce
! negative zeros and may signal a reproducibility issue over land.
! We resolve this by re-initializing and only evaluating over water points.
Expand Down Expand Up @@ -1098,7 +1098,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_KE_visc > 0) then
! Calculate the KE source from vertical viscosity [H L2 T-3 ~> m3 s-3].
! Calculate the KE source from vertical viscosity [H L2 T-3 ~> m3 s-3 or W m-2].
do k=1,nz
do j=js,je ; do I=Isq,Ieq
KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc(I,j,k)
Expand All @@ -1117,7 +1117,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_KE_visc_gl90 > 0) then
! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3].
! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3 or W m-2].
do k=1,nz
do j=js,je ; do I=Isq,Ieq
KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc_gl90(I,j,k)
Expand All @@ -1136,7 +1136,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_KE_stress > 0) then
! Calculate the KE source from surface stress (included in KE_visc) [H L2 T-3 ~> m3 s-3].
! Calculate the KE source from surface stress (included in KE_visc) [H L2 T-3 ~> m3 s-3 or W m-2].
do k=1,nz
do j=js,je ; do I=Isq,Ieq
KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_str(I,j,k)
Expand All @@ -1155,7 +1155,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_KE_horvisc > 0) then
! Calculate the KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3].
! Calculate the KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3 or W m-2].
do k=1,nz
do j=js,je ; do I=Isq,Ieq
KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%diffu(I,j,k)
Expand All @@ -1174,7 +1174,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS
endif

if (CS%id_KE_dia > 0) then
! Calculate the KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3].
! Calculate the KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3 or W m-2].
do k=1,nz
do j=js,je ; do I=Isq,Ieq
KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_dia(I,j,k)
Expand Down Expand Up @@ -1594,7 +1594,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_DEPTH", CS%mono_N2_depth, &
"The depth below which N2 is limited as monotonic for the "// &
"purposes of calculating the equivalent barotropic wave speed.", &
units='m', scale=US%m_to_Z, default=-1.)
units='m', scale=GV%m_to_H, default=-1.)
call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, &
"The fractional tolerance for finding the wave speeds.", &
units="nondim", default=0.001)
Expand Down
64 changes: 35 additions & 29 deletions src/diagnostics/MOM_wave_speed.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ module MOM_wave_speed
!! wave speed [nondim]. This parameter controls the default behavior of
!! wave_speed() which can be overridden by optional arguments.
real :: mono_N2_depth = -1. !< The depth below which N2 is limited as monotonic for the purposes of
!! calculating the equivalent barotropic wave speed [Z ~> m].
!! calculating the equivalent barotropic wave speed [H ~> m or kg m-2].
!! If this parameter is negative, this limiting does not occur.
!! This parameter controls the default behavior of wave_speed() which
!! can be overridden by optional arguments.
Expand Down Expand Up @@ -81,7 +81,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
!! for the purposes of calculating vertical modal structure [nondim].
real, optional, intent(in) :: mono_N2_depth !< A depth below which N2 is limited as
!! monotonic for the purposes of calculating vertical
!! modal structure [Z ~> m].
!! modal structure [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
optional, intent(out) :: modal_structure !< Normalized model structure [nondim]

Expand Down Expand Up @@ -157,9 +157,11 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
real :: sum_hc ! The sum of the layer thicknesses [Z ~> m]
real :: gp ! A limited local copy of gprime [L2 Z-1 T-2 ~> m s-2]
real :: N2min ! A minimum buoyancy frequency, including a slope rescaling factor [L2 Z-2 T-2 ~> s-2]
logical :: below_mono_N2_frac ! True if an interface is below the fractional depth where N2 should not increase.
logical :: below_mono_N2_depth ! True if an interface is below the absolute depth where N2 should not increase.
logical :: l_use_ebt_mode, calc_modal_structure
real :: l_mono_N2_column_fraction ! A local value of mono_N2_column_fraction [nondim]
real :: l_mono_N2_depth ! A local value of mono_N2_column_depth [Z ~> m]
real :: l_mono_N2_depth ! A local value of mono_N2_column_depth [H ~> m or kg m-2]
real :: mode_struct(SZK_(GV)) ! The mode structure [nondim], but it is also temporarily
! in units of [L2 T-2 ~> m2 s-2] after it is modified inside of tdma6.
real :: ms_min, ms_max ! The minimum and maximum mode structure values returned from tdma6 [L2 T-2 ~> m2 s-2]
Expand Down Expand Up @@ -214,18 +216,11 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
c2_scale = US%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.

min_h_frac = tol_Hfrac / real(nz)
!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,tv,&
!$OMP calc_modal_structure,l_use_ebt_mode,modal_structure, &
!$OMP l_mono_N2_column_fraction,l_mono_N2_depth,CS, &
!$OMP Z_to_pres,cg1,g_Rho0,rescale,I_rescale, &
!$OMP better_est,cg1_min2,tol_merge,tol_solve,c2_scale) &
!$OMP private(htot,hmin,kf,H_here,HxT_here,HxS_here,HxR_here, &
!$OMP Hf,Tf,Sf,Rf,pres,T_int,S_int,drho_dT,drho_dS, &
!$OMP drxh_sum,kc,Hc,Hc_H,Tc,Sc,I_Hnew,gprime,&
!$OMP Rc,speed2_tot,Igl,Igu,lam0,lam,lam_it,dlam, &
!$OMP mode_struct,sum_hc,N2min,gp,hw, &
!$OMP ms_min,ms_max,ms_sq,H_top,H_bot,I_Htot,merge, &
!$OMP det,ddet,det_it,ddet_it)
!$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,tv,use_EOS, &
!$OMP CS,min_h_frac,calc_modal_structure,l_use_ebt_mode, &
!$OMP modal_structure,l_mono_N2_column_fraction,l_mono_N2_depth, &
!$OMP Z_to_pres,cg1,g_Rho0,rescale,I_rescale,cg1_min2, &
!$OMP better_est,tol_solve,tol_merge,c2_scale)
do j=js,je
! First merge very thin layers with the one above (or below if they are
! at the top). This also transposes the row order so that columns can
Expand Down Expand Up @@ -335,7 +330,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * max(0.0,Rf(k,i)-Rf(k-1,i))
enddo
endif
endif
endif ! use_EOS

! Find gprime across each internal interface, taking care of convective instabilities by
! merging layers. If the estimated wave speed is too small, simply return zero.
Expand Down Expand Up @@ -452,24 +447,34 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
Igu(1) = 0. ! Neumann condition for pressure modes
sum_hc = Hc(1)
N2min = gprime(2)/Hc(1)

below_mono_N2_frac = .false.
below_mono_N2_depth = .false.
do k=2,kc
hw = 0.5*(Hc(k-1)+Hc(k))
gp = gprime(K)
if (l_mono_N2_column_fraction>0. .or. l_mono_N2_depth>=0.) then
!### Change to: if ( ((htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i)) .or. & ) )
if ( (((G%bathyT(i,j)+G%Z_ref) - sum_hc < l_mono_N2_column_fraction*(G%bathyT(i,j)+G%Z_ref)) .or. &
((l_mono_N2_depth >= 0.) .and. (sum_hc > l_mono_N2_depth))) .and. &
(gp > N2min*hw) ) then
! Filters out regions where N2 increases with depth but only in a lower fraction
! Determine whether N2 estimates should not be allowed to increase with depth.
if (l_mono_N2_column_fraction>0.) then
!### Change to: (htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i))
below_mono_N2_frac = ((G%bathyT(i,j)+G%Z_ref) - GV%H_to_Z*sum_hc < &
l_mono_N2_column_fraction*(G%bathyT(i,j)+G%Z_ref))
endif
if (l_mono_N2_depth >= 0.) below_mono_N2_depth = (sum_hc > GV%H_to_Z*l_mono_N2_depth)

if ( (gp > N2min*hw) .and. (below_mono_N2_frac .or. below_mono_N2_depth) ) then
! Filters out regions where N2 increases with depth, but only in a lower fraction
! of the water column or below a certain depth.
gp = N2min * hw
else
N2min = gp / hw
endif
endif

Igu(k) = 1.0/(gp*Hc(k))
Igl(k-1) = 1.0/(gp*Hc(k-1))
sum_hc = sum_hc + Hc(k)

if (better_est) then
! Estimate that the ebt_mode is sqrt(2) times the speed of the flat bottom modes.
speed2_tot = speed2_tot + 2.0 * gprime(K)*((H_top(K) * H_bot(K)) * I_Htot)
Expand Down Expand Up @@ -690,7 +695,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
H_top, & ! The distance of each filtered interface from the ocean surface [Z ~> m]
H_bot, & ! The distance of each filtered interface from the bottom [Z ~> m]
gprime, & ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].
N2 ! The Brunt Vaissalla freqency squared [T-2 ~> s-2]
N2 ! The buoyancy freqency squared [T-2 ~> s-2]
real, dimension(SZK_(GV),SZI_(G)) :: &
Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m]
Tf, & ! Layer temperatures after very thin layers are combined [C ~> degC]
Expand All @@ -704,7 +709,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
Sc, & ! A column of layer salinities after convective instabilities are removed [S ~> ppt]
Rc, & ! A column of layer densities after convective instabilities are removed [R ~> kg m-3]
Hc_H ! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2]
real :: I_Htot ! The inverse of the total filtered thicknesses [Z ~> m]
real :: I_Htot ! The inverse of the total filtered thicknesses [Z-1 ~> m-1]
real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant and its
! derivative with lam between rows of the Thomas algorithm solver [L2 s2 T-2 m-2 ~> nondim].
! The exact value should not matter for the final result if it is an even power of 2.
Expand Down Expand Up @@ -797,6 +802,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
! Simplifying the following could change answers at roundoff.
Z_to_pres = GV%Z_to_H * (GV%H_to_RZ * GV%g_Earth)
use_EOS = associated(tv%eqn_of_state)

if (CS%c1_thresh < 0.0) &
call MOM_error(FATAL, "INTERNAL_WAVE_CG1_THRESH must be set to a non-negative "//&
"value via wave_speed_init for wave_speeds to be used.")
Expand Down Expand Up @@ -978,7 +984,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
I_Hnew = 1.0 / (Hc(kc) + Hc(kc-1))
Tc(kc-1) = (Hc(kc)*Tc(kc) + Hc(kc-1)*Tc(kc-1)) * I_Hnew
Sc(kc-1) = (Hc(kc)*Sc(kc) + Hc(kc-1)*Sc(kc-1)) * I_Hnew
Hc(kc-1) = (Hc(kc) + Hc(kc-1))
Hc(kc-1) = Hc(kc) + Hc(kc-1)
kc = kc - 1
else ; exit ; endif
enddo
Expand Down Expand Up @@ -1006,7 +1012,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
if (merge) then
! Merge this layer with the one above and backtrack.
Rc(kc) = (Hc(kc)*Rc(kc) + Hf(k,i)*Rf(k,i)) / (Hc(kc) + Hf(k,i))
Hc(kc) = (Hc(kc) + Hf(k,i))
Hc(kc) = Hc(kc) + Hf(k,i)
! Backtrack to remove any convective instabilities above... Note
! that the tolerance is a factor of two larger, to avoid limit how
! far back we go.
Expand All @@ -1019,7 +1025,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
if (merge) then
! Merge the two bottommost layers. At this point kc = k2.
Rc(kc-1) = (Hc(kc)*Rc(kc) + Hc(kc-1)*Rc(kc-1)) / (Hc(kc) + Hc(kc-1))
Hc(kc-1) = (Hc(kc) + Hc(kc-1))
Hc(kc-1) = Hc(kc) + Hc(kc-1)
kc = kc - 1
else ; exit ; endif
enddo
Expand Down Expand Up @@ -1109,7 +1115,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s
do k=1,kc
w2avg = w2avg + 0.5*(mode_struct(K)**2+mode_struct(K+1)**2)*Hc(k) ![Z L4 T-4]
enddo
renorm = sqrt(htot(i)*a_int/w2avg) ![L-2 T-2]
renorm = sqrt(htot(i)*a_int/w2avg) ! [T2 L-2]
do K=1,kc+1 ; mode_struct(K) = renorm * mode_struct(K) ; enddo
! after renorm, mode_struct is again [nondim]

Expand Down Expand Up @@ -1437,7 +1443,7 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de
!! calculating the vertical modal structure [nondim].
real, optional, intent(in) :: mono_N2_depth !< The depth below which N2 is limited
!! as monotonic for the purposes of calculating the
!! vertical modal structure [Z ~> m].
!! vertical modal structure [H ~> m or kg m-2].
logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
!! that recover the remapping answers from 2018. Otherwise
!! use more robust but mathematically equivalent expressions.
Expand Down Expand Up @@ -1489,7 +1495,7 @@ subroutine wave_speed_set_param(CS, use_ebt_mode, mono_N2_column_fraction, mono_
!! calculating the vertical modal structure [nondim].
real, optional, intent(in) :: mono_N2_depth !< The depth below which N2 is limited
!! as monotonic for the purposes of calculating the
!! vertical modal structure [Z ~> m].
!! vertical modal structure [H ~> m or kg m-2].
logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
!! that recover the remapping answers from 2018. Otherwise
!! use more robust but mathematically equivalent expressions.
Expand Down
4 changes: 2 additions & 2 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1103,7 +1103,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
real :: oneOrTwo ! A variable that may be 1 or 2, depending on which form
! of the equatorial deformation radius us used [nondim]
real :: N2_filter_depth ! A depth below which stratification is treated as monotonic when
! calculating the first-mode wave speed [Z ~> m]
! calculating the first-mode wave speed [H ~> m or kg m-2]
real :: KhTr_passivity_coeff ! Coefficient setting the ratio between along-isopycnal tracer
! mixing and interface height mixing [nondim]
real :: absurdly_small_freq ! A miniscule frequency that is used to avoid division by 0 [T-1 ~> s-1]. The
Expand Down Expand Up @@ -1241,7 +1241,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
"The depth below which N2 is monotonized to avoid stratification "//&
"artifacts from altering the equivalent barotropic mode structure. "//&
"This monotonzization is disabled if this parameter is negative.", &
units="m", default=-1.0, scale=US%m_to_Z)
units="m", default=-1.0, scale=GV%m_to_H)
allocate(CS%ebt_struct(isd:ied,jsd:jed,GV%ke), source=0.0)
endif

Expand Down

0 comments on commit 878fd1e

Please sign in to comment.