From 326bec06c76e755fc50da7c34e09848e321a8f49 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 20 Jun 2023 19:15:28 -0400 Subject: [PATCH] +Refactor internal_tides interface Refactors the internal tide code in MOM_internal_tides and MOM_diabatic_driver to consolidate it in the MOM_internal_tides module and allow the control structure for that module to be made opaque. This includes moving the internal wave speed diagnostics and the call to wave_speeds or other code setting the internal wave speeds into propagate_int_tide. The get_param calls for INTERNAL_WAVE_CG1_THRESH and UNIFORM_TEST_CG were moved from the diabatic module to the MOM_internal_tides module. The wave_speed_CS and uniform_test_cg were removed from diabatic_CS and added to int_tide_CS. The Nb argument to propagate_int_tide has been made intent inout, as it is now usually set via the call to wave_speeds in that routine, but for certain tests it could use the value passed in from diabatic_driver. All answers are bitwise identical, but there are changes to public interfaces and types, and the order of some entries in the MOM_parameter_doc files and the available_diags files is changed for some cases. --- .../lateral/MOM_internal_tides.F90 | 65 ++++++++++++++++--- .../vertical/MOM_diabatic_driver.F90 | 57 +--------------- 2 files changed, 60 insertions(+), 62 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index ec07939ee4..8c56107a4f 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -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 @@ -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 @@ -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 @@ -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. @@ -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] @@ -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). @@ -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) @@ -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] @@ -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 "//& - "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 @@ -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) @@ -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') @@ -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 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index b0d04e434c..1bc29ee16f 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -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 @@ -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. @@ -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 @@ -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] @@ -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 @@ -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) @@ -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 @@ -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 "//& @@ -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", & @@ -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)