diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index dfd1048b82..1893859fe7 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -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 @@ -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. diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 8c56107a4f..788e922ff2 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -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 @@ -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 @@ -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. @@ -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 @@ -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] @@ -244,7 +247,7 @@ 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 @@ -252,8 +255,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, dt, & 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 @@ -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 + ! 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) @@ -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 @@ -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), & @@ -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 + 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 @@ -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.) @@ -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) @@ -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 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 631a47c259..dae52592e9 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -387,7 +387,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & CS%int_tide_input_CSp) call propagate_int_tide(h, tv, CS%int_tide_input%TKE_itidal_input, CS%int_tide_input%tideamp, & - CS%int_tide_input%Nb, dt, G, GV, US, CS%int_tide) + CS%int_tide_input%Nb, CS%int_tide_input%Rho_bot, dt, G, GV, US, CS%int_tide) if (showCallTree) call callTree_waypoint("done with propagate_int_tide (diabatic)") endif ! end CS%use_int_tides diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index 95e33929df..3da21b48fb 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -14,6 +14,7 @@ module MOM_int_tide_input use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_io, only : slasher, vardesc, MOM_read_data +use MOM_interface_heights, only : thickness_to_dz, find_rho_bottom use MOM_isopycnal_slopes, only : vert_fill_TS use MOM_time_manager, only : time_type, set_time, operator(+), operator(<=) use MOM_unit_scaling, only : unit_scale_type @@ -44,7 +45,8 @@ module MOM_int_tide_input !! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] real, allocatable, dimension(:,:) :: TKE_itidal_coef - !< The time-invariant field that enters the TKE_itidal input calculation [R Z3 T-2 ~> J m-2]. + !< The time-invariant field that enters the TKE_itidal input calculation noting that the + !! stratification and perhaps density are time-varying [R Z4 H-1 T-2 ~> J m-2 or J m kg-1]. character(len=200) :: inputdir !< The directory for input files. logical :: int_tide_source_test !< If true, apply an arbitrary generation site @@ -70,7 +72,8 @@ module MOM_int_tide_input TKE_itidal_input, & !< The internal tide TKE input at the bottom of the ocean [R Z3 T-3 ~> W m-2]. h2, & !< The squared topographic roughness height [Z2 ~> m2]. tideamp, & !< The amplitude of the tidal velocities [Z T-1 ~> m s-1]. - Nb !< The bottom stratification [T-1 ~> s-1]. + Nb, & !< The bottom stratification [T-1 ~> s-1]. + Rho_bot !< The bottom density or the Boussinesq reference density [R ~> kg m-3]. end type int_tide_input_type contains @@ -90,9 +93,12 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) !! to the internal tide sources. real, intent(in) :: dt !< The time increment [T ~> s]. type(int_tide_input_CS), pointer :: CS !< This module's control structure. + ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & N2_bot ! The bottom squared buoyancy frequency [T-2 ~> s-2]. + real, dimension(SZI_(G),SZJ_(G)) :: & + Rho_bot ! The average near-bottom density or the Boussinesq reference density [R ~> kg m-3]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & T_f, S_f ! The temperature and salinity in [C ~> degC] and [S ~> ppt] with the values in @@ -121,15 +127,25 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) call vert_fill_TS(h, tv%T, tv%S, CS%kappa_fill*dt, T_f, S_f, G, GV, US, larger_h_denom=.true.) endif - call find_N2_bottom(h, tv, T_f, S_f, itide%h2, fluxes, G, GV, US, N2_bot) + call find_N2_bottom(h, tv, T_f, S_f, itide%h2, fluxes, G, GV, US, N2_bot, Rho_bot) avg_enabled = query_averaging_enabled(CS%diag, time_end=time_end) - !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie - itide%Nb(i,j) = G%mask2dT(i,j) * sqrt(N2_bot(i,j)) - itide%TKE_itidal_input(i,j) = min(CS%TKE_itidal_coef(i,j)*itide%Nb(i,j), CS%TKE_itide_max) - enddo ; enddo + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + itide%Nb(i,j) = G%mask2dT(i,j) * sqrt(N2_bot(i,j)) + itide%TKE_itidal_input(i,j) = min(GV%Z_to_H*CS%TKE_itidal_coef(i,j)*itide%Nb(i,j), CS%TKE_itide_max) + enddo ; enddo + else + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + itide%Nb(i,j) = G%mask2dT(i,j) * sqrt(N2_bot(i,j)) + itide%Rho_bot(i,j) = G%mask2dT(i,j) * Rho_bot(i,j) + itide%TKE_itidal_input(i,j) = min((GV%RZ_to_H*Rho_bot(i,j)) * CS%TKE_itidal_coef(i,j)*itide%Nb(i,j), & + CS%TKE_itide_max) + enddo ; enddo + endif if (CS%int_tide_source_test) then itide%TKE_itidal_input(:,:) = 0.0 @@ -171,7 +187,7 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) end subroutine set_int_tide_input !> Estimates the near-bottom buoyancy frequency (N^2). -subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) +subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot, 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 @@ -186,55 +202,62 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes real, dimension(SZI_(G),SZJ_(G)), intent(out) :: N2_bot !< The squared buoyancy frequency at the !! ocean bottom [T-2 ~> s-2]. + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: rho_bot !< The average density near the ocean + !! bottom [R ~> kg m-3] ! Local variables real, dimension(SZI_(G),SZK_(GV)+1) :: & + pres, & ! The pressure at each interface [R L2 T-2 ~> Pa]. dRho_int ! The unfiltered density differences across interfaces [R ~> kg m-3]. + real, dimension(SZI_(G),SZK_(GV)) :: dz ! Layer thicknesses in depth units [Z ~> m] real, dimension(SZI_(G)) :: & - pres, & ! The pressure at each interface [R L2 T-2 ~> Pa]. Temp_int, & ! The temperature at each interface [C ~> degC] Salin_int, & ! The salinity at each interface [S ~> ppt] drho_bot, & ! The density difference at the bottom of a layer [R ~> kg m-3] h_amp, & ! The amplitude of topographic roughness [Z ~> m]. - hb, & ! The depth below a layer [Z ~> m]. - z_from_bot, & ! The height of a layer center above the bottom [Z ~> m]. + hb, & ! The thickness of the water column below the midpoint of a layer [H ~> m or kg m-2] + z_from_bot, & ! The distance of a layer center from the bottom [Z ~> m] dRho_dT, & ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] dRho_dS ! The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]. - real :: dz_int ! The thickness associated with an interface [Z ~> m]. - real :: G_Rho0 ! The gravitation acceleration divided by the Boussinesq - ! density [Z T-2 R-1 ~> m4 s-2 kg-1]. + real :: dz_int ! The vertical extent of water associated with an interface [Z ~> m] + real :: G_Rho0 ! The gravitational acceleration, sometimes divided by the Boussinesq + ! density [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2]. logical :: do_i(SZI_(G)), do_any integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - G_Rho0 = (US%L_to_Z**2*GV%g_Earth) / GV%Rho0 + G_Rho0 = (US%L_to_Z**2*GV%g_Earth) / GV%H_to_RZ EOSdom(:) = EOS_domain(G%HI) ! Find the (limited) density jump across each interface. do i=is,ie dRho_int(i,1) = 0.0 ; dRho_int(i,nz+1) = 0.0 enddo -!$OMP parallel do default(none) shared(is,ie,js,je,nz,tv,fluxes,G,GV,US,h,T_f,S_f, & -!$OMP h2,N2_bot,G_Rho0,EOSdom) & -!$OMP private(pres,Temp_Int,Salin_Int,dRho_dT,dRho_dS, & -!$OMP hb,dRho_bot,z_from_bot,do_i,h_amp, & -!$OMP do_any,dz_int) & -!$OMP firstprivate(dRho_int) + + !$OMP parallel do default(none) shared(is,ie,js,je,nz,tv,fluxes,G,GV,US,h,T_f,S_f, & + !$OMP h2,N2_bot,rho_bot,G_Rho0,EOSdom) & + !$OMP private(pres,Temp_Int,Salin_Int,dRho_dT,dRho_dS, & + !$OMP dz,hb,dRho_bot,z_from_bot,do_i,h_amp,do_any,dz_int) & + !$OMP firstprivate(dRho_int) do j=js,je + + ! Find the vertical distances across layers. + call thickness_to_dz(h, tv, dz, j, G, GV) + if (associated(tv%eqn_of_state)) then if (associated(fluxes%p_surf)) then - do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo + do i=is,ie ; pres(i,1) = fluxes%p_surf(i,j) ; enddo else - do i=is,ie ; pres(i) = 0.0 ; enddo + do i=is,ie ; pres(i,1) = 0.0 ; enddo endif do K=2,nz do i=is,ie - pres(i) = pres(i) + (GV%g_Earth*GV%H_to_RZ)*h(i,j,k-1) + pres(i,K) = pres(i,K-1) + (GV%g_Earth*GV%H_to_RZ)*h(i,j,k-1) Temp_Int(i) = 0.5 * (T_f(i,j,k) + T_f(i,j,k-1)) Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) enddo - call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT(:), dRho_dS(:), & + call calculate_density_derivs(Temp_int, Salin_int, pres(:,K), dRho_dT(:), dRho_dS(:), & tv%eqn_of_state, EOSdom) do i=is,ie dRho_int(i,K) = max(dRho_dT(i)*(T_f(i,j,k) - T_f(i,j,k-1)) + & @@ -250,7 +273,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) ! Find the bottom boundary layer stratification. do i=is,ie hb(i) = 0.0 ; dRho_bot(i) = 0.0 - z_from_bot(i) = 0.5*GV%H_to_Z*h(i,j,nz) + z_from_bot(i) = 0.5*dz(i,nz) do_i(i) = (G%mask2dT(i,j) > 0.0) h_amp(i) = sqrt(h2(i,j)) enddo @@ -258,16 +281,16 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) do k=nz,2,-1 do_any = .false. do i=is,ie ; if (do_i(i)) then - dz_int = 0.5*GV%H_to_Z*(h(i,j,k) + h(i,j,k-1)) + dz_int = 0.5*(dz(i,k) + dz(i,k-1)) z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above - hb(i) = hb(i) + dz_int + hb(i) = hb(i) + 0.5*(h(i,j,k) + h(i,j,k-1)) dRho_bot(i) = dRho_bot(i) + dRho_int(i,K) if (z_from_bot(i) > h_amp(i)) then if (k>2) then ! Always include at least one full layer. - hb(i) = hb(i) + 0.5*GV%H_to_Z*(h(i,j,k-1) + h(i,j,k-2)) + hb(i) = hb(i) + 0.5*(h(i,j,k-1) + h(i,j,k-2)) dRho_bot(i) = dRho_bot(i) + dRho_int(i,K-1) endif do_i(i) = .false. @@ -283,6 +306,15 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) N2_bot(i,j) = (G_Rho0 * dRho_bot(i)) / hb(i) else ; N2_bot(i,j) = 0.0 ; endif enddo + + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + do i=is,ie + rho_bot(i,j) = GV%Rho0 + enddo + else + ! Average the density over the envelope of the topography. + call find_rho_bottom(h, dz, pres, h_amp, tv, j, G, GV, US, Rho_bot(:,j)) + endif enddo end subroutine find_N2_bottom @@ -359,6 +391,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) units="m s-1", default=0.0, scale=US%m_s_to_L_T) allocate(itide%Nb(isd:ied,jsd:jed), source=0.0) + allocate(itide%Rho_bot(isd:ied,jsd:jed), source=0.0) allocate(itide%h2(isd:ied,jsd:jed), source=0.0) allocate(itide%TKE_itidal_input(isd:ied,jsd:jed), source=0.0) allocate(itide%tideamp(isd:ied,jsd:jed), source=utide) @@ -452,8 +485,8 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) if (max_frac_rough >= 0.0) & itide%h2(i,j) = min((max_frac_rough*(G%bathyT(i,j)+G%Z_ref))**2, itide%h2(i,j)) - ! Compute the fixed part of internal tidal forcing; units are [R Z3 T-2 ~> J m-2] here. - CS%TKE_itidal_coef(i,j) = 0.5*US%L_to_Z*kappa_h2_factor*GV%Rho0*& + ! Compute the fixed part of internal tidal forcing; units are [R Z4 H-1 T-2 ~> J m-2 or J m kg-1] here. + CS%TKE_itidal_coef(i,j) = 0.5*US%L_to_Z*kappa_h2_factor * GV%H_to_RZ * & kappa_itides * itide%h2(i,j) * itide%tideamp(i,j)**2 enddo ; enddo