Skip to content

Commit

Permalink
add underflow management, hchksum reproduce
Browse files Browse the repository at this point in the history
  • Loading branch information
Raphael Dussin authored and Raphael Dussin committed Aug 6, 2024
1 parent 45e29e7 commit ceb92ee
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 27 deletions.
10 changes: 5 additions & 5 deletions src/core/MOM_interface_heights.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -394,23 +394,23 @@ end subroutine find_col_avg_SpV

!> 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, h_bot, k_bot)
subroutine find_rho_bottom(G, GV, US, tv, h, dz, pres_int, dz_avg, j, Rho_bot, h_bot, k_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
type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
!! thermodynamic fields.
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].
real, dimension(SZI_(G)), optional, intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2]
integer, dimension(SZI_(G)), optional, intent(out) :: k_bot !< Bottom boundary layer top layer index
real, dimension(SZI_(G)), intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2]
integer, dimension(SZI_(G)), intent(out) :: k_bot !< Bottom boundary layer top layer index

! Local variables
real :: hb(SZI_(G)) ! Running sum of the thickness in the bottom boundary layer [H ~> m or kg m-2]
Expand Down
31 changes: 27 additions & 4 deletions src/parameterizations/lateral/MOM_internal_tides.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,8 @@ module MOM_internal_tides
real, allocatable, dimension(:,:,:) :: int_N2w2 !< Depth-integrated Brunt Vaissalla freqency times
!! vertical profile squared, for each mode [H T-2 ~> m s-2 or kg m-2 s-2]
real :: q_itides !< fraction of local dissipation [nondim]
real :: En_sum !< global sum of energy for use in debugging, in MKS units [J]
real :: En_sum !< global sum of energy for use in debugging, in MKS units [H Z2 T-2 L2 ~> m5 s-2 or J]
real :: En_underflow !< A minuscule amount of energy [H Z2 T-2 ~> m3 s-2 or J m-2]
type(time_type), pointer :: Time => NULL() !< A pointer to the model's clock.
character(len=200) :: inputdir !< directory to look for coastline angle file
real :: decay_rate !< A constant rate at which internal tide energy is
Expand Down Expand Up @@ -1047,6 +1048,11 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C
call post_data(CS%id_En_mode(fr,m), tot_En, CS%diag)
endif ; enddo ; enddo

! fix underflows
do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
if (abs(CS%En(i,j,a,fr,m)) < CS%En_underflow) CS%En(i,j,a,fr,m) = 0.0
enddo ; enddo ; enddo ; enddo ; enddo

! split energy array into multiple restarts
do fr=1,CS%Nfreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied
CS%En_restart_mode1(i,j,a,fr) = CS%En(i,j,a,fr,1)
Expand Down Expand Up @@ -1367,8 +1373,8 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G)), intent(in) :: h_bot !< Bottom boundary layer thickness [H ~> m]
integer, dimension(SZI_(G)), intent(in) :: k_bot !< Bottom boundary layer top layer index [nondim]
real, dimension(SZI_(G)), intent(in) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2]
integer, dimension(SZI_(G)), intent(in) :: k_bot !< Bottom boundary layer top layer index
integer, intent(in) :: j !< The j-index to work on
real, dimension(SZI_(G),SZK_(GV)), intent(in) :: N2_lay !< The squared buoyancy frequency of the
!! layers [T-2 ~> s-2].
Expand Down Expand Up @@ -2162,6 +2168,11 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss)
LB%jsh = jsh ; LB%jeh = jeh ; LB%ish = ish ; LB%ieh = ieh
call propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, CS%nAngle, CS, LB, residual_loss)

! fix underflows
do a=1,na ; do j=jsh,jeh ; do i=ish,ieh
if (abs(En(i,j,a)) < CS%En_underflow) En(i,j,a) = 0.0
enddo ; enddo ; enddo

if (CS%debug) then
do m=1,CS%nMode ; do fr=1,CS%Nfreq
call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: after propagate_x')
Expand All @@ -2184,6 +2195,11 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss)
LB%jsh = jsh ; LB%jeh = jeh ; LB%ish = ish ; LB%ieh = ieh
call propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, CS%nAngle, CS, LB, residual_loss)

! fix underflows
do a=1,na ; do j=jsh,jeh ; do i=ish,ieh
if (abs(En(i,j,a)) < CS%En_underflow) En(i,j,a) = 0.0
enddo ; enddo ; enddo

call pass_var(En, G%domain)
call pass_var(residual_loss, G%domain)

Expand All @@ -2202,6 +2218,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss)
if (is_root_pe()) print *, 'propagate: bottom of routine', CS%En_sum
enddo ; enddo
endif

end subroutine propagate

!> This subroutine does first-order corner advection. It was written with the hopes
Expand Down Expand Up @@ -3360,6 +3377,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
real :: HZ2_T2_to_J_m2 ! unit conversion factor for Energy from internal to mks
real :: HZ2_T3_to_W_m2 ! unit conversion factor for TKE from internal to mks
real :: W_m2_to_HZ2_T3 ! unit conversion factor for TKE from mks to internal
real :: J_m2_to_HZ2_T2 ! unit conversion factor for Energy from mks to internal
integer :: num_angle, num_freq, num_mode, m, fr
integer :: isd, ied, jsd, jed, a, id_ang, i, j, nz
type(axes_grp) :: axes_ang
Expand All @@ -3384,6 +3402,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
HZ2_T2_to_J_m2 = GV%H_to_kg_m2*(US%Z_to_m**2)*(US%s_to_T**2)
HZ2_T3_to_W_m2 = GV%H_to_kg_m2*(US%Z_to_m**2)*(US%s_to_T**3)
W_m2_to_HZ2_T3 = GV%kg_m2_to_H*(US%m_to_Z**2)*(US%T_to_s**3)
J_m2_to_HZ2_T2 = GV%kg_m2_to_H*(US%m_to_Z**2)*(US%T_to_s**2)

CS%initialized = .true.

Expand Down Expand Up @@ -3550,7 +3569,11 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
call get_param(param_file, mdl, "EN_CHECK_TOLERANCE", CS%En_check_tol, &
"An energy density tolerance for flagging points with small negative "//&
"internal tide energy.", &
units="J m-2", default=1.0, scale=W_m2_to_HZ2_T3, &
units="J m-2", default=1.0, scale=J_m2_to_HZ2_T2, &
do_not_log=.not.CS%apply_Froude_drag)
call get_param(param_file, mdl, "EN_UNDERFLOW", CS%En_underflow, &
"A small energy density below which Energy is set to zero.", &
units="J m-2", default=1.0e-100, scale=J_m2_to_HZ2_T2, &
do_not_log=.not.CS%apply_Froude_drag)
call get_param(param_file, mdl, "CDRAG", CS%cdrag, &
"CDRAG is the drag coefficient relating the magnitude of "//&
Expand Down
39 changes: 22 additions & 17 deletions src/parameterizations/vertical/MOM_internal_tide_input.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ module MOM_int_tide_input
logical :: debug !< If true, write verbose checksums for debugging.
type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
!! regulate the timing of diagnostic output.
real :: TKE_itide_max !< Maximum Internal tide conversion
real :: TKE_itide_maxi !< Maximum Internal tide conversion
!! available to mix above the BBL [H Z2 T-3 ~> m3 s-3 or W m-2]
real :: kappa_fill !< Vertical diffusivity used to interpolate sensible values
!! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
Expand Down Expand Up @@ -105,7 +105,9 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS)
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].
Rho_bot, & ! The average near-bottom density or the Boussinesq reference density [R ~> kg m-3].
h_bot ! Bottom boundary layer thickness [H ~> m or kg m-2].
integer, dimension(SZI_(G),SZJ_(G)) :: k_bot ! Bottom boundary layer top layer index.

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
Expand Down Expand Up @@ -140,7 +142,7 @@ 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, Rho_bot)
call find_N2_bottom(G, GV, US, tv, fluxes, h, T_f, S_f, itide%h2, N2_bot, Rho_bot, h_bot, k_bot)

avg_enabled = query_averaging_enabled(CS%diag, time_end=time_end)

Expand All @@ -149,15 +151,15 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS)
do fr=1,CS%nFreq ; do j=js,je ; do i=is,ie
itide%Nb(i,j) = G%mask2dT(i,j) * sqrt(N2_bot(i,j))
CS%TKE_itidal_input(i,j,fr) = min(GV%RZ_to_H*GV%Z_to_H*CS%TKE_itidal_coef(i,j,fr)*itide%Nb(i,j), &
CS%TKE_itide_max)
CS%TKE_itide_maxi)
enddo ; enddo ; enddo
else
!$OMP parallel do default(shared)
do fr=1,CS%nFreq ; 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)
CS%TKE_itidal_input(i,j,fr) = min((GV%RZ_to_H*GV%RZ_to_H*Rho_bot(i,j))*CS%TKE_itidal_coef(i,j,fr)*itide%Nb(i,j), &
CS%TKE_itide_max)
CS%TKE_itide_maxi)
enddo ; enddo ; enddo
endif

Expand Down Expand Up @@ -204,23 +206,26 @@ 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, rho_bot)
subroutine find_N2_bottom(G, GV, US, tv, fluxes, h, T_f, S_f, h2, N2_bot, Rho_bot, h_bot, k_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]
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to the
!! thermodynamic fields
type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: T_f !< Temperature after vertical filtering to
!! smooth out the values in thin layers [C ~> degC].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: S_f !< Salinity after vertical filtering to
!! smooth out the values in thin layers [S ~> ppt].
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h2 !< Bottom topographic roughness [Z2 ~> m2].
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
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Rho_bot !< The average density near the ocean
!! bottom [R ~> kg m-3]
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2]
integer, dimension(SZI_(G),SZJ_(G)), intent(out) :: k_bot !< Bottom boundary layer top layer index

! Local variables
real, dimension(SZI_(G),SZK_(GV)+1) :: &
pres, & ! The pressure at each interface [R L2 T-2 ~> Pa].
Expand All @@ -234,9 +239,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot, rho_bo
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].
h_bot ! The bottom boundary layer thickness [H ~> m or kg m-2]
integer, dimension(SZI_(G)) :: k_bot ! Bottom boundary layer top layer index
dRho_dS ! The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-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
Expand All @@ -255,7 +258,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot, rho_bo
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,rho_bot,h_bot,k_bot,G_Rho0,EOSdom) &
!$OMP h2,N2_bot,Rho_bot,h_bot,k_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)
Expand Down Expand Up @@ -332,7 +335,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot, rho_bo
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), h_bot, k_bot)
call find_rho_bottom(G, GV, US, tv, h, dz, pres, h_amp, j, Rho_bot(:,j), h_bot(:,j), k_bot(:,j))
endif
enddo

Expand Down Expand Up @@ -401,6 +404,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)
real :: kappa_h2_factor ! factor for the product of wavenumber * rms sgs height [nondim].
real :: kappa_itides ! topographic wavenumber and non-dimensional scaling [L-1 ~> m-1]
real :: min_zbot_itides ! Minimum ocean depth for internal tide conversion [Z ~> m].
real :: HZ2_T3_to_W_m2 ! unit conversion factor for TKE from internal to mks
real :: W_m2_to_HZ2_T3 ! unit conversion factor for TKE from mks to internal
integer :: tlen_days !< Time interval from start for adding wave source
!! for testing internal tides (BDM)
Expand All @@ -427,7 +431,8 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)

CS%diag => diag

W_m2_to_HZ2_T3 = GV%m_to_H*(US%m_to_Z**2)*(US%T_to_s**3)
HZ2_T3_to_W_m2 = GV%H_to_kg_m2*(US%Z_to_m**2)*(US%s_to_T**3)
W_m2_to_HZ2_T3 = GV%kg_m2_to_H*(US%m_to_Z**2)*(US%T_to_s**3)

! Read all relevant parameters and write them to the model log.
call log_version(param_file, mdl, version, "")
Expand Down Expand Up @@ -467,7 +472,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)
call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, &
"A scaling factor for the roughness amplitude with "//&
"INT_TIDE_DISSIPATION.", units="nondim", default=1.0)
call get_param(param_file, mdl, "TKE_ITIDE_MAX", CS%TKE_itide_max, &
call get_param(param_file, mdl, "TKE_ITIDE_MAX", CS%TKE_itide_maxi, &
"The maximum internal tide energy source available to mix "//&
"above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
units="W m-2", default=1.0e3, scale=W_m2_to_HZ2_T3)
Expand Down Expand Up @@ -563,7 +568,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)
write(var_descript, '("Internal Tide Driven Turbulent Kinetic Energy in frequency ",i1)') fr

CS%id_TKE_itidal_itide(fr) = register_diag_field('ocean_model',var_name,diag%axesT1,Time, &
var_descript, 'W m-2', conversion=US%RZ3_T3_to_W_m2)
var_descript, 'W m-2', conversion=HZ2_T3_to_W_m2)
enddo

CS%id_Nb = register_diag_field('ocean_model','Nb_itide',diag%axesT1,Time, &
Expand Down
2 changes: 1 addition & 1 deletion src/parameterizations/vertical/MOM_set_diffusivity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1197,7 +1197,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &

! Average over the larger of the envelope of the topography or a minimal distance.
do i=is,ie ; dz_BBL_avg(i) = max(h_amp(i), CS%dz_BBL_avg_min) ; enddo
call find_rho_bottom(h, dz, pres, dz_BBL_avg, tv, j, G, GV, US, Rho_bot, h_bot, k_bot)
call find_rho_bottom(G, GV, US, tv, h, dz, pres, dz_BBL_avg, j, Rho_bot, h_bot, k_bot)

end subroutine find_N2

Expand Down

0 comments on commit ceb92ee

Please sign in to comment.