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

+Refactor internal_tides interface #383

Merged
merged 2 commits into from
Jun 28, 2023
Merged

+Refactor internal_tides interface #383

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
326bec0
+Refactor internal_tides interface
Hallberg-NOAA Jun 20, 2023
5299478
Merge branch 'dev/gfdl' into refactor_internal_tide_interface
marshallward Jun 27, 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
65 changes: 57 additions & 8 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 @@ -23,6 +23,7 @@ module MOM_internal_tides
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : surface, thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
use MOM_wave_speed, only : wave_speeds, wave_speed_CS, wave_speed_init

implicit none ; private

Expand All @@ -33,12 +34,14 @@ module MOM_internal_tides
public get_lowmode_loss

!> This control structure has parameters for the MOM_internal_tides module
type, public :: int_tide_CS
type, public :: int_tide_CS ; private
logical :: do_int_tides !< If true, use the internal tide code.
integer :: nFreq = 0 !< The number of internal tide frequency bands
integer :: nMode = 1 !< The number of internal tide vertical modes
integer :: nAngle = 24 !< The number of internal tide angular orientations
integer :: energized_angle = -1 !< If positive, only this angular band is energized for debugging purposes
real :: uniform_test_cg !< Uniform group velocity of internal tide
!! for testing internal tides [L T-1 ~> m s-1]
logical :: corner_adv !< If true, use a corner advection rather than PPM.
logical :: upwind_1st !< If true, use a first-order upwind scheme.
logical :: simple_2nd !< If true, use a simple second order (arithmetic mean) interpolation
Expand Down Expand Up @@ -137,11 +140,14 @@ module MOM_internal_tides
!< The internal wave energy density as a function of (i,j,angle); temporary for restart
real, allocatable, dimension(:) :: frequency !< The frequency of each band [T-1 ~> s-1].

type(wave_speed_CS) :: wave_speed !< Wave speed control structure
type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate the
!! timing of diagnostic output.

!>@{ Diag handles
! Diag handles relevant to all modes, frequencies, and angles
integer :: id_cg1 = -1 ! diagnostic handle for mode-1 speed
integer, allocatable, dimension(:) :: id_cn ! diagnostic handle for all mode speeds
integer :: id_tot_En = -1, id_TKE_itidal_input = -1, id_itide_drag = -1
integer :: id_refl_pref = -1, id_refl_ang = -1, id_land_mask = -1
integer :: id_trans = -1, id_residual = -1
Expand Down Expand Up @@ -181,7 +187,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, cn, TKE_itidal_input, vel_btTide, Nb, dt, &
subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, 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 @@ -194,16 +200,18 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, &
!! internal waves [R Z3 T-3 ~> W m-2].
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: vel_btTide !< Barotropic velocity read
!! from file [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Nb !< Near-bottom buoyancy frequency [T-1 ~> s-1].
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, 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
real, dimension(SZI_(G),SZJ_(G),CS%nMode), &
intent(in) :: cn !< The internal wave speeds of each
!! mode [L T-1 ~> m s-1].

! Local variables
real, dimension(SZI_(G),SZJ_(G),2) :: &
test ! A test unit vector used to determine grid rotation in halos [nondim]
real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: &
cn ! baroclinic internal gravity wave speeds for each mode [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: &
tot_En_mode, & ! energy summed over angles only [R Z3 T-2 ~> J m-2]
Ub, & ! near-bottom horizontal velocity of wave (modal) [L T-1 ~> m s-1]
Expand Down Expand Up @@ -254,6 +262,18 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, &
Ub(:,:,:,:) = 0.
Umax(:,:,:,:) = 0.

cn(:,:,:) = 0.

! Set properties related to the internal tides, such as the wave speeds, storing some
! of them in the control structure for this module.
if (CS%uniform_test_cg > 0.0) then
do m=1,CS%nMode ; cn(:,:,m) = CS%uniform_test_cg ; enddo
else
call wave_speeds(h, tv, G, GV, US, CS%nMode, cn, CS%wave_speed, &
CS%w_struct, CS%u_struct, CS%u_struct_max, CS%u_struct_bot, &
Nb, CS%int_w2, CS%int_U2, CS%int_N2w2, full_halos=.true.)
endif

! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.**********************
! This is wrong, of course, but it works reasonably in some cases.
! Uncomment if wave_speed is not used to calculate the true values (BDM).
Expand Down Expand Up @@ -596,6 +616,10 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, &
call enable_averages(dt, time_end, CS%diag)

if (query_averaging_enabled(CS%diag)) then
! Output internal wave modal wave speeds
if (CS%id_cg1 > 0) call post_data(CS%id_cg1, cn(:,:,1),CS%diag)
do m=1,CS%nMode ; if (CS%id_cn(m) > 0) call post_data(CS%id_cn(m), cn(:,:,m), CS%diag) ; enddo

! Output two-dimensional diagnostics
if (CS%id_tot_En > 0) call post_data(CS%id_tot_En, tot_En, CS%diag)
if (CS%id_itide_drag > 0) call post_data(CS%id_itide_drag, drag_scale, CS%diag)
Expand Down Expand Up @@ -2292,6 +2316,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
real, dimension(:,:), allocatable :: ridge_temp ! array for temporary storage of flags
! of cells with double-reflecting ridges [nondim]
logical :: use_int_tides, use_temperature
real :: IGW_c1_thresh ! A threshold first mode internal wave speed below which all higher
! mode speeds are not calculated but simply assigned a speed of 0 [L T-1 ~> m s-1].
real :: kappa_h2_factor ! A roughness scaling factor [nondim]
real :: RMS_roughness_frac ! The maximum RMS topographic roughness as a fraction of the
! nominal ocean depth, or a negative value for no limit [nondim]
Expand Down Expand Up @@ -2322,8 +2348,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
use_temperature = .true.
call read_param(param_file, "ENABLE_THERMODYNAMICS", use_temperature)
if (.not.use_temperature) call MOM_error(FATAL, &
"register_int_tide_restarts: internal_tides only works with "//&
marshallward marked this conversation as resolved.
Show resolved Hide resolved
"ENABLE_THERMODYNAMICS defined.")
"internal_tides_init: internal_tides only works with ENABLE_THERMODYNAMICS defined.")

! Set number of frequencies, angles, and modes to consider
num_freq = 1 ; num_angle = 24 ; num_mode = 1
Expand Down Expand Up @@ -2447,6 +2472,15 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
"CDRAG is the drag coefficient relating the magnitude of "//&
"the velocity field to the bottom stress.", &
units="nondim", default=0.003)
call get_param(param_file, mdl, "INTERNAL_WAVE_CG1_THRESH", IGW_c1_thresh, &
"A minimal value of the first mode internal wave speed below which all higher "//&
"mode speeds are not calculated but are simply reported as 0. This must be "//&
"non-negative for the wave_speeds routine to be used.", &
units="m s-1", default=0.01, scale=US%m_s_to_L_T)

call get_param(param_file, mdl, "UNIFORM_TEST_CG", CS%uniform_test_cg, &
"If positive, a uniform group velocity of internal tide for test case", &
default=-1., units="m s-1", scale=US%m_s_to_L_T)
call get_param(param_file, mdl, "INTERNAL_TIDE_ENERGIZED_ANGLE", CS%energized_angle, &
"If positive, only one angular band of the internal tides "//&
"gets all of the energy. (This is for debugging.)", default=-1)
Expand Down Expand Up @@ -2610,6 +2644,18 @@ 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


! Register maps of reflection parameters
CS%id_refl_ang = register_diag_field('ocean_model', 'refl_angle', diag%axesT1, &
Time, 'Local angle of coastline/ridge/shelf with respect to equator', 'rad')
Expand Down Expand Up @@ -2777,6 +2823,9 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)

enddo

! Initialize the module that calculates the wave speeds.
call wave_speed_init(CS%wave_speed, c1_thresh=IGW_c1_thresh)

end subroutine internal_tides_init

!> This subroutine deallocates the memory associated with the internal tides control structure
Expand Down
57 changes: 3 additions & 54 deletions src/parameterizations/vertical/MOM_diabatic_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,6 @@ module MOM_diabatic_driver
use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs
use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d
use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units
use MOM_wave_speed, only : wave_speeds, wave_speed_CS, wave_speed_init
use MOM_wave_interface, only : wave_parameters_CS
use MOM_stochastics, only : stochastic_CS

Expand Down Expand Up @@ -123,9 +122,6 @@ module MOM_diabatic_driver
!! shear and ePBL diffusivities are used.
real :: ePBL_Prandtl !< The Prandtl number used by ePBL to convert vertical
!! diffusivities into viscosities [nondim].
integer :: nMode = 1 !< Number of baroclinic modes to consider
real :: uniform_test_cg !< Uniform group velocity of internal tide
!! for testing internal tides [L T-1 ~> m s-1]
logical :: useALEalgorithm !< If true, use the ALE algorithm rather than layered
!! isopycnal/stacked shallow water mode. This logical
!! passed by argument to diabatic_driver_init.
Expand Down Expand Up @@ -239,7 +235,6 @@ module MOM_diabatic_driver
type(int_tide_CS) :: int_tide !< Internal tide control structure
type(opacity_CS) :: opacity !< Opacity control structure
type(regularize_layers_CS) :: regularize_layers !< Regularize layer control structure
type(wave_speed_CS) :: wave_speed !< Wave speed control struct

type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass
type(group_pass_type) :: pass_Kv !< For group halo pass
Expand Down Expand Up @@ -297,8 +292,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
! local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
eta ! Interface heights before diapycnal mixing [Z ~> m]
real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: &
cn_IGW ! baroclinic internal gravity wave speeds [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: temp_diag ! Previous temperature for diagnostics [C ~> degC]
real, dimension(SZI_(G)) :: T_freeze, & ! The freezing potential temperature at the current salinity [C ~> degC].
ps ! Surface pressure [R L2 T-2 ~> Pa]
Expand Down Expand Up @@ -392,17 +385,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
! This block provides an interface for the unresolved low-mode internal tide module.
call set_int_tide_input(u, v, h, tv, fluxes, CS%int_tide_input, dt, G, GV, US, &
CS%int_tide_input_CSp)
cn_IGW(:,:,:) = 0.0
if (CS%uniform_test_cg > 0.0) then
do m=1,CS%nMode ; cn_IGW(:,:,m) = CS%uniform_test_cg ; enddo
else
call wave_speeds(h, tv, G, GV, US, CS%nMode, cn_IGW, CS%wave_speed, CS%int_tide%w_struct, &
CS%int_tide%u_struct, CS%int_tide%u_struct_max, CS%int_tide%u_struct_bot, &
CS%int_tide_input%Nb, CS%int_tide%int_w2, CS%int_tide%int_U2, CS%int_tide%int_N2w2, &
full_halos=.true.)
endif

call propagate_int_tide(h, tv, cn_IGW, CS%int_tide_input%TKE_itidal_input, CS%int_tide_input%tideamp, &
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)
if (showCallTree) call callTree_waypoint("done with propagate_int_tide (diabatic)")
endif ! end CS%use_int_tides
Expand Down Expand Up @@ -505,10 +489,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
call diagnoseMLDbyEnergy((/CS%id_MLD_EN1, CS%id_MLD_EN2, CS%id_MLD_EN3/),&
h, tv, G, GV, US, CS%MLD_En_vals, CS%diag)
endif
if (CS%use_int_tides) then
if (CS%id_cg1 > 0) call post_data(CS%id_cg1, cn_IGW(:,:,1),CS%diag)
do m=1,CS%nMode ; if (CS%id_cn(m) > 0) call post_data(CS%id_cn(m), cn_IGW(:,:,m), CS%diag) ; enddo
endif

if (stoch_CS%do_sppt .and. stoch_CS%id_sppt_wts > 0) &
call post_data(stoch_CS%id_sppt_wts, stoch_CS%sppt_wts, CS%diag)
Expand Down Expand Up @@ -2979,8 +2959,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di

! Local variables
real :: Kd ! A diffusivity used in the default for other tracer diffusivities [Z2 T-1 ~> m2 s-1]
real :: IGW_c1_thresh ! A threshold first mode internal wave speed below which all higher
! mode speeds are not calculated but simply assigned a speed of 0 [L T-1 ~> m s-1].
logical :: use_temperature
character(len=20) :: EN1, EN2, EN3

Expand Down Expand Up @@ -3073,23 +3051,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
call get_param(param_file, mdl, "INTERNAL_TIDES", CS%use_int_tides, &
"If true, use the code that advances a separate set of "//&
"equations for the internal tide energy density.", default=.false.)
CS%nMode = 1
if (CS%use_int_tides) then
call get_param(param_file, mdl, "INTERNAL_TIDE_MODES", CS%nMode, &
"The number of distinct internal tide modes "//&
"that will be calculated.", default=1, do_not_log=.true.)
call get_param(param_file, mdl, "INTERNAL_WAVE_CG1_THRESH", IGW_c1_thresh, &
"A minimal value of the first mode internal wave speed below which all higher "//&
"mode speeds are not calculated but are simply reported as 0. This must be "//&
"non-negative for the wave_speeds routine to be used.", &
units="m s-1", default=0.01, scale=US%m_s_to_L_T)
call get_param(param_file, mdl, "UNIFORM_TEST_CG", CS%uniform_test_cg, &
"If positive, a uniform group velocity of internal tide for test case", &
default=-1., units="m s-1", scale=US%m_s_to_L_T)
endif

call get_param(param_file, mdl, "MASSLESS_MATCH_TARGETS", &
CS%massless_match_targets, &

call get_param(param_file, mdl, "MASSLESS_MATCH_TARGETS", CS%massless_match_targets, &
"If true, the temperature and salinity of massless layers "//&
"are kept consistent with their target densities. "//&
"Otherwise the properties of massless layers evolve "//&
Expand Down Expand Up @@ -3197,19 +3160,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz)
endif

if (CS%use_int_tides) then
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
endif

if (use_temperature) then
CS%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff", diag%axesTi, &
Time, "Diffusive diapycnal temperature flux across interfaces", &
Expand Down Expand Up @@ -3504,7 +3454,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
call int_tide_input_init(Time, G, GV, US, param_file, diag, CS%int_tide_input_CSp, &
CS%int_tide_input)
call internal_tides_init(Time, G, GV, US, param_file, diag, CS%int_tide)
call wave_speed_init(CS%wave_speed, c1_thresh=IGW_c1_thresh)
endif

physical_OBL_scheme = (CS%use_bulkmixedlayer .or. CS%use_KPP .or. CS%use_energetic_PBL)
Expand Down