Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

+wave_speed arg mono_N2_depth in thickness units #424

Merged
merged 3 commits into from
Jul 26, 2023
Merged

+wave_speed arg mono_N2_depth in thickness units #424

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
60089b3
+wave_speed arg mono_N2_depth in thickness units
Hallberg-NOAA Jul 24, 2023
7012d26
Merge branch 'dev/gfdl' into mono_N2_depth_units
marshallward Jul 25, 2023
1943a9f
Merge branch 'dev/gfdl' into mono_N2_depth_units
marshallward Jul 26, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 12 additions & 12 deletions src/diagnostics/MOM_diagnostics.F90
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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