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

Jump to bottom

+*Use actual density for internal tide bottom drag #415

Merged
merged 3 commits into from
Jul 31, 2023
Merged

+*Use actual density for internal tide bottom drag #415

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e62a930
+Add find_rho_bottom
Hallberg-NOAA Jul 20, 2023
c0e9b26
+*Non-Boussinesq internal tide drag uses density
Hallberg-NOAA Jul 20, 2023
4666feb
Merge branch 'dev/gfdl' into nonBous_internal_tide_input
marshallward Jul 31, 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
132 changes: 132 additions & 0 deletions src/core/MOM_interface_heights.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ module MOM_interface_heights

public find_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple
public calc_derived_thermo
public find_rho_bottom

!> Calculates the heights of the free surface or all interfaces from layer thicknesses.
interface find_eta
Expand Down Expand Up @@ -313,6 +314,137 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)

end subroutine calc_derived_thermo


!> Determine the in situ density averaged over a specified distance from the bottom,
!! calculating it as the inverse of the mass-weighted average specific volume.
subroutine find_rho_bottom(h, dz, pres_int, dz_avg, tv, j, G, GV, US, Rho_bot)
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),SZK_(GV)), &
intent(in) :: dz !< Height change across layers [Z ~> m]
real, dimension(SZI_(G),SZK_(GV)+1), &
intent(in) :: pres_int !< Pressure at each interface [R L2 T-2 ~> Pa]
real, dimension(SZI_(G)), intent(in) :: dz_avg !< The vertical distance over which to average [Z ~> m]
type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
!! thermodynamic fields.
integer, intent(in) :: j !< j-index of row to work on
real, dimension(SZI_(G)), intent(out) :: Rho_bot !< Near-bottom density [R ~> kg m-3].

! Local variables
real :: hb(SZI_(G)) ! Running sum of the thickness in the bottom boundary layer [H ~> m or kg m-2]
real :: SpV_h_bot(SZI_(G)) ! Running sum of the specific volume times thickness in the bottom
! boundary layer [R-1 H ~> m4 kg-1 or m]
real :: dz_bbl_rem(SZI_(G)) ! Vertical extent of the boundary layer that has yet to be accounted
! for [Z ~> m]
real :: h_bbl_frac(SZI_(G)) ! Thickness of the fractional layer that makes up the top of the
! boundary layer [H ~> m or kg m-2]
real :: T_bbl(SZI_(G)) ! Temperature of the fractional layer that makes up the top of the
! boundary layer [C ~> degC]
real :: S_bbl(SZI_(G)) ! Salinity of the fractional layer that makes up the top of the
! boundary layer [S ~> ppt]
real :: P_bbl(SZI_(G)) ! Pressure the top of the boundary layer [R L2 T-2 ~> Pa]
real :: dp(SZI_(G)) ! Pressure change across the fractional layer that makes up the top
! of the boundary layer [R L2 T-2 ~> Pa]
real :: SpV_bbl(SZI_(G)) ! In situ specific volume of the fractional layer that makes up the
! top of the boundary layer [R-1 ~> m3 kg-1]
real :: frac_in ! The fraction of a layer that is within the bottom boundary layer [nondim]
logical :: do_i(SZI_(G)), do_any
logical :: use_EOS
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
integer :: i, k, is, ie, nz

is = G%isc ; ie = G%iec ; nz = GV%ke

use_EOS = associated(tv%T) .and. associated(tv%S) .and. associated(tv%eqn_of_state)

if (GV%Boussinesq .or. GV%semi_Boussinesq .or. .not.allocated(tv%SpV_avg)) then
do i=is,ie
rho_bot(i) = GV%Rho0
enddo
else
! Check that SpV_avg has been set.
if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, &
"find_rho_bottom called in fully non-Boussinesq mode with invalid values of SpV_avg.")

! Set the bottom density to the inverse of the in situ specific volume averaged over the
! specified distance, with care taken to avoid having compressibility lead to an imprint
! of the layer thicknesses on this density.
do i=is,ie
hb(i) = 0.0 ; SpV_h_bot(i) = 0.0
dz_bbl_rem(i) = G%mask2dT(i,j) * max(0.0, dz_avg(i))
do_i(i) = .true.
if (G%mask2dT(i,j) <= 0.0) then
! Set acceptable values for calling the equation of state over land.
T_bbl(i) = 0.0 ; S_bbl(i) = 0.0 ; dp(i) = 0.0 ; P_bbl(i) = 0.0
SpV_bbl(i) = 1.0 ! This value is arbitrary, provided it is non-zero.
h_bbl_frac(i) = 0.0
do_i(i) = .false.
endif
enddo

do k=nz,1,-1
do_any = .false.
do i=is,ie ; if (do_i(i)) then
if (dz(i,k) < dz_bbl_rem(i)) then
! This layer is fully within the averaging depth.
SpV_h_bot(i) = SpV_h_bot(i) + h(i,j,k) * tv%SpV_avg(i,j,k)
dz_bbl_rem(i) = dz_bbl_rem(i) - dz(i,k)
hb(i) = hb(i) + h(i,j,k)
do_any = .true.
else
if (dz(i,k) > 0.0) then
frac_in = dz_bbl_rem(i) / dz(i,k)
else
frac_in = 0.0
endif
if (use_EOS) then
! Store the properties of this layer to determine the average
! specific volume of the portion that is within the BBL.
T_bbl(i) = tv%T(i,j,k) ; S_bbl(i) = tv%S(i,j,k)
dp(i) = frac_in * (GV%g_Earth*GV%H_to_RZ * h(i,j,k))
P_bbl(i) = pres_int(i,K) + (1.0-frac_in) * (GV%g_Earth*GV%H_to_RZ * h(i,j,k))
else
SpV_bbl(i) = tv%SpV_avg(i,j,k)
endif
h_bbl_frac(i) = frac_in * h(i,j,k)
dz_bbl_rem(i) = 0.0
do_i(i) = .false.
endif
endif ; enddo
if (.not.do_any) exit
enddo
do i=is,ie ; if (do_i(i)) then
! The nominal bottom boundary layer is thicker than the water column, but layer 1 is
! already included in the averages. These values are set so that the call to find
! the layer-average specific volume will behave sensibly.
if (use_EOS) then
T_bbl(i) = tv%T(i,j,1) ; S_bbl(i) = tv%S(i,j,1)
dp(i) = 0.0
P_bbl(i) = pres_int(i,1)
else
SpV_bbl(i) = tv%SpV_avg(i,j,1)
endif
h_bbl_frac(i) = 0.0
endif ; enddo

if (use_EOS) then
! Find the average specific volume of the fractional layer atop the BBL.
EOSdom(:) = EOS_domain(G%HI)
call average_specific_vol(T_bbl, S_bbl, P_bbl, dp, SpV_bbl, tv%eqn_of_state, EOSdom)
endif

do i=is,ie
if (hb(i) + h_bbl_frac(i) < GV%H_subroundoff) h_bbl_frac(i) = GV%H_subroundoff
rho_bot(i) = G%mask2dT(i,j) * (hb(i) + h_bbl_frac(i)) / (SpV_h_bot(i) + h_bbl_frac(i)*SpV_bbl(i))
enddo
endif

end subroutine find_rho_bottom


!> Converts thickness from geometric height units to thickness units, perhaps via an
!! inversion of the integral of the density in pressure using variables stored in
!! the thermo_var_ptrs type when in non-Boussinesq mode.
Expand Down
91 changes: 55 additions & 36 deletions src/parameterizations/lateral/MOM_internal_tides.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,9 +78,10 @@ module MOM_internal_tides
real, allocatable, dimension(:,:,:,:,:) :: TKE_Froude_loss
!< energy lost due to wave breaking [R Z3 T-3 ~> W m-2]
real, allocatable, dimension(:,:) :: TKE_itidal_loss_fixed
!< Fixed part of the energy lost due to small-scale drag [R L-2 Z3 ~> kg m-2] here;
!! This will be multiplied by N and the squared near-bottom velocity to get
!! the energy losses in [R Z3 T-3 ~> W m-2]
!< Fixed part of the energy lost due to small-scale drag [R Z3 L-2 ~> kg m-2] here;
!! This will be multiplied by N and the squared near-bottom velocity (and by
!! the near-bottom density in non-Boussinesq mode) to get the energy losses
!! in [R Z4 H-1 L-2 ~> kg m-2 or m]
real, allocatable, dimension(:,:,:,:,:) :: TKE_itidal_loss
!< energy lost due to small-scale wave drag [R Z3 T-3 ~> W m-2]
real, allocatable, dimension(:,:,:,:,:) :: TKE_residual_loss
Expand Down Expand Up @@ -120,7 +121,7 @@ module MOM_internal_tides
real :: cdrag !< The bottom drag coefficient [nondim].
real :: drag_min_depth !< The minimum total ocean thickness that will be used in the denominator
!! of the quadratic drag terms for internal tides when
!! INTERNAL_TIDE_QUAD_DRAG is true [Z ~> m]
!! INTERNAL_TIDE_QUAD_DRAG is true [H ~> m or kg m-2]
logical :: apply_background_drag
!< If true, apply a drag due to background processes as a sink.
logical :: apply_bottom_drag
Expand Down Expand Up @@ -187,7 +188,7 @@ module MOM_internal_tides

!> Calls subroutines in this file that are needed to refract, propagate,
!! and dissipate energy density of the internal tide.
subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, &
subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, dt, &
G, GV, US, CS)
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
Expand All @@ -203,6 +204,8 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, &
real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: Nb !< Near-bottom buoyancy frequency [T-1 ~> s-1].
!! In some cases the input values are used, but in
!! others this is set along with the wave speeds.
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Rho_bot !< Near-bottom density or the Boussinesq
!! reference density [R ~> kg m-3].
real, intent(in) :: dt !< Length of time over which to advance
!! the internal tides [T ~> s].
type(int_tide_CS), intent(inout) :: CS !< Internal tide control structure
Expand All @@ -228,8 +231,8 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, &
real :: frac_per_sector ! The inverse of the number of angular, modal and frequency bins [nondim]
real :: f2 ! The squared Coriolis parameter interpolated to a tracer point [T-2 ~> s-2]
real :: Kmag2 ! A squared horizontal wavenumber [L-2 ~> m-2]
real :: I_D_here ! The inverse of the local depth [Z-1 ~> m-1]
real :: I_rho0 ! The inverse fo the Boussinesq density [R-1 ~> m3 kg-1]
real :: I_D_here ! The inverse of the local water column thickness [H-1 ~> m-1 or m2 kg-1]
real :: I_mass ! The inverse of the local water mass [R-1 Z-1 ~> m2 kg-1]
real :: freq2 ! The frequency squared [T-2 ~> s-2]
real :: PE_term ! total potential energy of profile [R Z ~> kg m-2]
real :: KE_term ! total kinetic energy of profile [R Z ~> kg m-2]
Expand All @@ -244,16 +247,15 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, &
real :: En_initial, Delta_E_check ! Energies for debugging [R Z3 T-2 ~> J m-2]
real :: TKE_Froude_loss_check, TKE_Froude_loss_tot ! Energy losses for debugging [R Z3 T-3 ~> W m-2]
character(len=160) :: mesg ! The text of an error message
integer :: a, m, fr, i, j, k, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm
integer :: a, m, fr, i, j, k, is, ie, js, je, isd, ied, jsd, jed, nAngle
integer :: id_g, jd_g ! global (decomp-invar) indices (for debugging)
type(group_pass_type), save :: pass_test, pass_En
type(time_type) :: time_end
logical:: avg_enabled

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nAngle = CS%NAngle
nzm = GV%ke
I_rho0 = 1.0 / GV%Rho0

cn_subRO = 1e-30*US%m_s_to_L_T
en_subRO = 1e-30*US%W_m2_to_RZ3_T3*US%s_to_T

Expand Down Expand Up @@ -429,11 +431,20 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, &
do k=1,GV%ke ; do j=jsd,jed ; do i=isd,ied
htot(i,j) = htot(i,j) + h(i,j,k)
enddo ; enddo ; enddo
do j=jsd,jed ; do i=isd,ied
I_D_here = 1.0 / (max(GV%H_to_Z*htot(i,j), CS%drag_min_depth))
drag_scale(i,j) = CS%cdrag * sqrt(max(0.0, US%L_to_Z**2*vel_btTide(i,j)**2 + &
tot_En(i,j) * I_rho0 * I_D_here)) * I_D_here
enddo ; enddo
if (GV%Boussinesq) then
marshallward marked this conversation as resolved.
Show resolved Hide resolved
! This is mathematically equivalent to the form in the option below, but they differ at roundoff.
do j=jsd,jed ; do i=isd,ied
I_D_here = 1.0 / (max(htot(i,j), CS%drag_min_depth))
drag_scale(i,j) = CS%cdrag * sqrt(max(0.0, US%L_to_Z**2*vel_btTide(i,j)**2 + &
tot_En(i,j) * GV%RZ_to_H * I_D_here)) * GV%Z_to_H*I_D_here
enddo ; enddo
else
do j=jsd,jed ; do i=isd,ied
I_mass = GV%RZ_to_H / (max(htot(i,j), CS%drag_min_depth))
drag_scale(i,j) = (CS%cdrag * (Rho_bot(i,j)*I_mass)) * &
sqrt(max(0.0, US%L_to_Z**2*vel_btTide(i,j)**2 + tot_En(i,j) * I_mass))
enddo ; enddo
endif
do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
! to each En component (technically not correct; fix later)
Expand Down Expand Up @@ -504,7 +515,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, &
! Finally, apply loss
if (CS%apply_wave_drag) then
! Calculate loss rate and apply loss over the time step
call itidal_lowmode_loss(G, US, CS, Nb, Ub, CS%En, CS%TKE_itidal_loss_fixed, &
call itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, CS%En, CS%TKE_itidal_loss_fixed, &
CS%TKE_itidal_loss, dt, full_halos=.false.)
endif
! Check for En<0 - for debugging, delete later
Expand Down Expand Up @@ -782,18 +793,21 @@ end subroutine sum_En

!> Calculates the energy lost from the propagating internal tide due to
!! scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001).
subroutine itidal_lowmode_loss(G, US, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos)
subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos)
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
type(int_tide_CS), intent(in) :: CS !< Internal tide control structure
real, dimension(G%isd:G%ied,G%jsd:G%jed), &
intent(in) :: Nb !< Near-bottom stratification [T-1 ~> s-1].
real, dimension(G%isd:G%ied,G%jsd:G%jed), &
intent(in) :: Rho_bot !< Near-bottom density [R ~> kg m-3].
real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%nFreq,CS%nMode), &
intent(inout) :: Ub !< RMS (over one period) near-bottom horizontal
!! mode velocity [L T-1 ~> m s-1].
real, dimension(G%isd:G%ied,G%jsd:G%jed), &
intent(in) :: TKE_loss_fixed !< Fixed part of energy loss [R L-2 Z3 ~> kg m-2]
!! (rho*kappa*h^2).
intent(in) :: TKE_loss_fixed !< Fixed part of energy loss [R Z4 H-1 L-2 ~> kg m-2 or m]
!! (rho*kappa*h^2) or (kappa*h^2).
real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
intent(inout) :: En !< Energy density of the internal waves [R Z3 T-2 ~> J m-2].
real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
Expand Down Expand Up @@ -830,14 +844,18 @@ subroutine itidal_lowmode_loss(G, US, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss,
enddo

! Calculate TKE loss rate; units of [R Z3 T-3 ~> W m-2] here.
TKE_loss_tot = q_itides * TKE_loss_fixed(i,j) * Nb(i,j) * Ub(i,j,fr,m)**2
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
marshallward marked this conversation as resolved.
Show resolved Hide resolved
TKE_loss_tot = q_itides * GV%Z_to_H * TKE_loss_fixed(i,j) * Nb(i,j) * Ub(i,j,fr,m)**2
else
TKE_loss_tot = q_itides * (GV%RZ_to_H * Rho_bot(i,j)) * TKE_loss_fixed(i,j) * Nb(i,j) * Ub(i,j,fr,m)**2
endif

! Update energy remaining (this is a pseudo implicit calc)
! (E(t+1)-E(t))/dt = -TKE_loss(E(t+1)/E(t)), which goes to zero as E(t+1) goes to zero
if (En_tot > 0.0) then
do a=1,CS%nAngle
frac_per_sector = En(i,j,a,fr,m)/En_tot
TKE_loss(i,j,a,fr,m) = frac_per_sector*TKE_loss_tot ! Wm-2
TKE_loss(i,j,a,fr,m) = frac_per_sector*TKE_loss_tot ! [R Z3 T-3 ~> W m-2]
loss_rate = TKE_loss(i,j,a,fr,m) / (En(i,j,a,fr,m) + En_negl) ! [T-1 ~> s-1]
En(i,j,a,fr,m) = En(i,j,a,fr,m) / (1.0 + dt*loss_rate)
enddo
Expand Down Expand Up @@ -2458,8 +2476,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
call get_param(param_file, mdl, "INTERNAL_TIDE_DRAG_MIN_DEPTH", CS%drag_min_depth, &
"The minimum total ocean thickness that will be used in the denominator "//&
"of the quadratic drag terms for internal tides.", &
units="m", default=1.0, scale=US%m_to_Z, do_not_log=.not.CS%apply_bottom_drag)
CS%drag_min_depth = MAX(CS%drag_min_depth, GV%H_subroundoff * GV%H_to_Z)
units="m", default=1.0, scale=GV%m_to_H, do_not_log=.not.CS%apply_bottom_drag)
CS%drag_min_depth = MAX(CS%drag_min_depth, GV%H_subroundoff)
call get_param(param_file, mdl, "INTERNAL_TIDE_FROUDE_DRAG", CS%apply_Froude_drag, &
"If true, apply wave breaking as a sink.", &
default=.false.)
Expand Down Expand Up @@ -2543,9 +2561,10 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
else
h2(i,j) = max(h2(i,j), 0.0)
endif
! Compute the fixed part; units are [R L-2 Z3 ~> kg m-2] here
! will be multiplied by N and the squared near-bottom velocity to get into [R Z3 T-3 ~> W m-2]
CS%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*GV%Rho0 * US%L_to_Z*kappa_itides * h2(i,j)
! Compute the fixed part; units are [R Z4 H-1 L-2 ~> kg m-2 or m] here
! will be multiplied by N and the squared near-bottom velocity (and by the
! near-bottom density in non-Boussinesq mode) to get into [R Z3 T-3 ~> W m-2]
CS%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor* GV%H_to_RZ * US%L_to_Z*kappa_itides * h2(i,j)
enddo ; enddo

deallocate(h2)
Expand Down Expand Up @@ -2644,16 +2663,16 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
enddo
call pass_var(CS%residual,G%domain)

CS%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, &
Time, 'First baroclinic mode (eigen) speed', 'm s-1', conversion=US%L_T_to_m_s)
allocate(CS%id_cn(CS%nMode), source=-1)
do m=1,CS%nMode
write(var_name, '("cn_mode",i1)') m
write(var_descript, '("Baroclinic (eigen) speed of mode ",i1)') m
CS%id_cn(m) = register_diag_field('ocean_model',var_name, diag%axesT1, &
Time, var_descript, 'm s-1', conversion=US%L_T_to_m_s)
call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
enddo
CS%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, &
Time, 'First baroclinic mode (eigen) speed', 'm s-1', conversion=US%L_T_to_m_s)
allocate(CS%id_cn(CS%nMode), source=-1)
do m=1,CS%nMode
write(var_name, '("cn_mode",i1)') m
write(var_descript, '("Baroclinic (eigen) speed of mode ",i1)') m
CS%id_cn(m) = register_diag_field('ocean_model',var_name, diag%axesT1, &
Time, var_descript, 'm s-1', conversion=US%L_T_to_m_s)
call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
enddo


! Register maps of reflection parameters
Expand Down
Loading