diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 172d2459d5..a8b0d3f813 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -16,11 +16,13 @@ module MOM_internal_tides use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type +use MOM_int_tide_input, only: int_tide_input_CS, get_input_TKE, get_barotropic_tidal_vel use MOM_io, only : slasher, MOM_read_data, file_exists, axis_info use MOM_io, only : set_axis_info, get_axis_info use MOM_restart, only : register_restart_field, MOM_restart_CS, restart_init, save_restart use MOM_restart, only : lock_check, restart_registry_lock use MOM_spatial_means, only : global_area_integral +use MOM_string_functions, only: extract_real use MOM_time_manager, only : time_type, time_type_to_real, operator(+), operator(/), operator(-) use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface, thermo_var_ptrs @@ -54,6 +56,9 @@ module MOM_internal_tides !! the default is false; it is always true with aggress_adjust. logical :: use_PPMang !< If true, use PPM for advection of energy in angular space. + real, allocatable, dimension(:,:) :: fraction_tidal_input + !< how the energy from one tidal component is distributed + !! over the various vertical modes, 2d in frequency and mode [nondim] real, allocatable, dimension(:,:) :: refl_angle !< local coastline/ridge/shelf angles read from file [rad] ! (could be in G control structure) @@ -161,7 +166,7 @@ module MOM_internal_tides ! 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_tot_En = -1 integer :: id_refl_pref = -1, id_refl_ang = -1, id_land_mask = -1 integer :: id_trans = -1, id_residual = -1 integer :: id_dx_Cv = -1, id_dy_Cu = -1 @@ -172,7 +177,12 @@ module MOM_internal_tides integer, allocatable, dimension(:,:) :: & id_En_mode, & id_itidal_loss_mode, & + id_leak_loss_mode, & + id_quad_loss_mode, & + id_Froude_loss_mode, & + id_residual_loss_mode, & id_allprocesses_loss_mode, & + id_itide_drag, & id_Ub_mode, & id_cp_mode ! Diag handles considering: all modes, frequencies, and angles @@ -180,6 +190,7 @@ module MOM_internal_tides id_En_ang_mode, & id_itidal_loss_ang_mode integer, allocatable, dimension(:) :: & + id_TKE_itidal_input, & id_Ustruct_mode, & id_Wstruct_mode, & id_int_w2_mode, & @@ -200,8 +211,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, Rho_bot, dt, & - G, GV, US, CS) +subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_CSp, 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. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -209,10 +219,6 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Pointer to thermodynamic variables !! (needed for wave structure). - real, dimension(SZI_(G),SZJ_(G)), intent(in) :: TKE_itidal_input !< The energy input to the - !! 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(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. @@ -220,9 +226,14 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, !! 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_input_CS), intent(in) :: inttide_input_CSp !< Internal tide input control structure type(int_tide_CS), intent(inout) :: CS !< Internal tide control structure ! Local variables + real, dimension(SZI_(G),SZJ_(G),CS%nFreq) :: & + TKE_itidal_input, & !< The energy input to the internal waves [R Z3 T-3 ~> W m-2]. + vel_btTide !< Barotropic velocity read from file [L T-1 ~> m s-1]. + 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) :: & @@ -231,15 +242,22 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, 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] Umax ! Maximum horizontal velocity of wave (modal) [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: & + drag_scale ! bottom drag scale [T-1 ~> s-1] real, dimension(SZI_(G),SZJ_(G)) :: & + tot_vel_btTide2, & tot_En, & ! energy summed over angles, modes, frequencies [R Z3 T-2 ~> J m-2] tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_Froude_loss, tot_residual_loss, tot_allprocesses_loss, & ! energy loss rates summed over angle, freq, and mode [R Z3 T-3 ~> W m-2] htot, & ! The vertical sum of the layer thicknesses [H ~> m or kg m-2] - drag_scale, & ! bottom drag scale [T-1 ~> s-1] itidal_loss_mode, & ! Energy lost due to small-scale wave drag, summed over angles [R Z3 T-3 ~> W m-2] + leak_loss_mode, & + quad_loss_mode, & + Froude_loss_mode, & + residual_loss_mode, & allprocesses_loss_mode ! Total energy loss rates for a given mode and frequency (summed over ! all angles) [R Z3 T-3 ~> W m-2] + 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] @@ -273,7 +291,10 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, en_subRO = 1e-30*US%W_m2_to_RZ3_T3*US%s_to_T ! initialize local arrays - drag_scale(:,:) = 0. + TKE_itidal_input(:,:,:) = 0. + vel_btTide(:,:,:) = 0. + tot_vel_btTide2(:,:) = 0. + drag_scale(:,:,:,:) = 0. Ub(:,:,:,:) = 0. Umax(:,:,:,:) = 0. @@ -329,24 +350,27 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, !enddo ; enddo ; enddo ! Add the forcing.*************************************************************** + + call get_input_TKE(G, TKE_itidal_input, CS%nFreq, inttide_input_CSp) + if (CS%energized_angle <= 0) then - frac_per_sector = 1.0 / real(CS%nAngle * CS%nMode * CS%nFreq) + frac_per_sector = 1.0 / real(CS%nAngle) do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie f2 = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)) if (CS%frequency(fr)**2 > f2) & CS%En(i,j,a,fr,m) = CS%En(i,j,a,fr,m) + dt*frac_per_sector*(1.0-CS%q_itides) * & - TKE_itidal_input(i,j) + CS%fraction_tidal_input(fr,m) * TKE_itidal_input(i,j,fr) enddo ; enddo ; enddo ; enddo ; enddo elseif (CS%energized_angle <= CS%nAngle) then - frac_per_sector = 1.0 / real(CS%nMode * CS%nFreq) + frac_per_sector = 1.0 a = CS%energized_angle do m=1,CS%nMode ; do fr=1,CS%nFreq ; do j=js,je ; do i=is,ie f2 = 0.25*((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & (G%CoriolisBu(I-1,J)**2 + G%CoriolisBu(I,J-1)**2)) if (CS%frequency(fr)**2 > f2) & CS%En(i,j,a,fr,m) = CS%En(i,j,a,fr,m) + dt*frac_per_sector*(1.0-CS%q_itides) * & - TKE_itidal_input(i,j) + CS%fraction_tidal_input(fr,m) * TKE_itidal_input(i,j,fr) enddo ; enddo ; enddo ; enddo else call MOM_error(WARNING, "Internal tide energy is being put into a angular "//& @@ -397,6 +421,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! Propagate the waves. do m=1,CS%nMode ; do fr=1,CS%Nfreq + ! initialize residual loss, will be computed in propagate CS%TKE_residual_loss(:,:,:,fr,m) = 0. call propagate(CS%En(:,:,:,fr,m), cn(:,:,m), CS%frequency(fr), dt, & @@ -479,29 +504,37 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! Extract the energy for mixing due to bottom drag------------------------------- if (CS%apply_bottom_drag) then - do j=js,je ; do i=is,ie ; htot(i,j) = 0.0 ; enddo ; enddo - do k=1,GV%ke ; do j=js,je ; do i=is,ie + do j=jsd,jed ; do i=isd,ied ; htot(i,j) = 0.0 ; enddo ; enddo + + call get_barotropic_tidal_vel(G, vel_btTide, CS%nFreq, inttide_input_CSp) + + do fr=1,CS%Nfreq ; do j=jsd,jed ; do i=isd,ied + tot_vel_btTide2(i,j) = tot_vel_btTide2(i,j) + vel_btTide(i,j,fr)**2 + enddo ; enddo ; enddo + + 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 if (GV%Boussinesq) then ! This is mathematically equivalent to the form in the option below, but they differ at roundoff. - do j=js,je ; do i=is,ie + do m=1,CS%NMode ; do fr=1,CS%Nfreq ; 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 + drag_scale(i,j,fr,m) = CS%cdrag * sqrt(max(0.0, US%L_to_Z**2*tot_vel_btTide2(i,j)**2 + & + tot_En_mode(i,j,fr,m) * GV%RZ_to_H * I_D_here)) * GV%Z_to_H*I_D_here + enddo ; enddo ; enddo ; enddo else - do j=js,je ; do i=is,ie + do m=1,CS%NMode ; do fr=1,CS%Nfreq ; 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 + drag_scale(i,j,fr,m) = (CS%cdrag * (Rho_bot(i,j)*I_mass)) * & + sqrt(max(0.0, US%L_to_Z**2*tot_vel_btTide2(i,j)**2 + & + tot_En_mode(i,j,fr,m) * I_mass)) + enddo ; enddo ; enddo ; enddo endif do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale ! to each En component (technically not correct; fix later) - CS%TKE_quad_loss(i,j,a,fr,m) = CS%En(i,j,a,fr,m) * drag_scale(i,j) ! loss rate - CS%En(i,j,a,fr,m) = CS%En(i,j,a,fr,m) / (1.0 + dt * drag_scale(i,j)) ! implicit update + CS%TKE_quad_loss(i,j,a,fr,m) = CS%En(i,j,a,fr,m) * drag_scale(i,j,fr,m) ! loss rate + CS%En(i,j,a,fr,m) = CS%En(i,j,a,fr,m) / (1.0 + dt * drag_scale(i,j,fr,m)) ! implicit update enddo ; enddo ; enddo ; enddo ; enddo endif ! Check for En<0 - for debugging, delete later @@ -685,9 +718,14 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! 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) - if (CS%id_TKE_itidal_input > 0) call post_data(CS%id_TKE_itidal_input, & - TKE_itidal_input, CS%diag) + do fr=1,CS%nFreq + if (CS%id_TKE_itidal_input(fr) > 0) call post_data(CS%id_TKE_itidal_input(fr), & + TKE_itidal_input(:,:,fr), CS%diag) + enddo + + do m=1,CS%nMode ; do fr=1,CS%nFreq + if (CS%id_itide_drag(fr,m) > 0) call post_data(CS%id_itide_drag(fr,m), drag_scale(:,:,fr,m), CS%diag) + enddo ; enddo ! Output 2-D energy density (summed over angles) for each frequency and mode do m=1,CS%nMode ; do fr=1,CS%Nfreq ; if (CS%id_En_mode(fr,m) > 0) then @@ -780,15 +818,27 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, do m=1,CS%nMode ; do fr=1,CS%Nfreq if (CS%id_itidal_loss_mode(fr,m) > 0 .or. CS%id_allprocesses_loss_mode(fr,m) > 0) then itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well) + leak_loss_mode(:,:) = 0.0 + quad_loss_mode(:,:) = 0.0 + Froude_loss_mode(:,:) = 0.0 + residual_loss_mode(:,:) = 0.0 allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together do a=1,CS%nAngle ; do j=js,je ; do i=is,ie itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + CS%TKE_itidal_loss(i,j,a,fr,m) + leak_loss_mode(i,j) = leak_loss_mode(i,j) + CS%TKE_leak_loss(i,j,a,fr,m) + quad_loss_mode(i,j) = quad_loss_mode(i,j) + CS%TKE_quad_loss(i,j,a,fr,m) + Froude_loss_mode(i,j) = Froude_loss_mode(i,j) + CS%TKE_Froude_loss(i,j,a,fr,m) + residual_loss_mode(i,j) = residual_loss_mode(i,j) + CS%TKE_residual_loss(i,j,a,fr,m) allprocesses_loss_mode(i,j) = allprocesses_loss_mode(i,j) + & ((((CS%TKE_leak_loss(i,j,a,fr,m) + CS%TKE_quad_loss(i,j,a,fr,m)) + & CS%TKE_itidal_loss(i,j,a,fr,m)) + CS%TKE_Froude_loss(i,j,a,fr,m)) + & CS%TKE_residual_loss(i,j,a,fr,m)) enddo ; enddo ; enddo call post_data(CS%id_itidal_loss_mode(fr,m), itidal_loss_mode, CS%diag) + call post_data(CS%id_leak_loss_mode(fr,m), leak_loss_mode, CS%diag) + call post_data(CS%id_quad_loss_mode(fr,m), quad_loss_mode, CS%diag) + call post_data(CS%id_Froude_loss_mode(fr,m), Froude_loss_mode, CS%diag) + call post_data(CS%id_residual_loss_mode(fr,m), residual_loss_mode, CS%diag) call post_data(CS%id_allprocesses_loss_mode(fr,m), allprocesses_loss_mode, CS%diag) endif ; enddo ; enddo @@ -2501,6 +2551,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) 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] real :: period_1 ! The period of the gravest modeled mode [T ~> s] + real :: period ! A tidal period read from namelist [T ~> s] 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 @@ -2516,6 +2567,9 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) character(len=200) :: h2_file character(len=80) :: rough_var ! Input file variable names + character(len=240), dimension(:), allocatable :: energy_fractions + character(len=240) :: periods + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed nz = GV%ke @@ -2539,17 +2593,29 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) if (.not.((num_freq > 0) .and. (num_angle > 0) .and. (num_mode > 0))) return CS%nFreq = num_freq ; CS%nAngle = num_angle ; CS%nMode = num_mode + allocate(energy_fractions(num_freq)) + allocate(CS%fraction_tidal_input(num_freq,num_mode)) + + call read_param(param_file, "ENERGY_FRACTION_PER_MODE", energy_fractions) + + do fr=1,num_freq ; do m=1,num_mode + CS%fraction_tidal_input(fr,m) = extract_real(energy_fractions(fr), " ,", m, 0.) + enddo ; enddo + ! Allocate phase speed array allocate(CS%cp(isd:ied, jsd:jed, num_freq, num_mode), source=0.0) ! Allocate and populate frequency array (each a multiple of first for now) allocate(CS%frequency(num_freq)) - call get_param(param_file, mdl, "FIRST_MODE_PERIOD", period_1, & - "The period of the first mode for internal tides", default=44567., & - units="s", scale=US%s_to_T) + + + ! The periods of the tidal constituents for internal tides raytracing + call read_param(param_file, "TIDAL_PERIODS", periods) do fr=1,num_freq - CS%frequency(fr) = (8.0*atan(1.0) * (real(fr)) / period_1) ! ADDED BDM + period = extract_real(periods, " ,", fr, 0.) + if (period == 0.) call MOM_error(FATAL, "MOM_internal_tides: invalid tidal period") + CS%frequency(fr) = 8.0*atan(1.0)/period enddo ! Read all relevant parameters and write them to the model log. @@ -2858,14 +2924,18 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) CS%id_tot_En = register_diag_field('ocean_model', 'ITide_tot_En', diag%axesT1, & Time, 'Internal tide total energy density', & 'J m-2', conversion=US%RZ3_T3_to_W_m2*US%T_to_s) - ! Register 2-D drag scale used for quadratic bottom drag - CS%id_itide_drag = register_diag_field('ocean_model', 'ITide_drag', diag%axesT1, & - Time, 'Interior and bottom drag internal tide decay timescale', 's-1', conversion=US%s_to_T) - !Register 2-D energy input into internal tides - CS%id_TKE_itidal_input = register_diag_field('ocean_model', 'TKE_itidal_input', diag%axesT1, & - Time, 'Conversion from barotropic to baroclinic tide, '//& - 'a fraction of which goes into rays', & - 'W m-2', conversion=US%RZ3_T3_to_W_m2) + + allocate(CS%id_itide_drag(CS%nFreq, CS%nMode), source=-1) + allocate(CS%id_TKE_itidal_input(CS%nFreq), source=-1) + do fr=1,CS%nFreq + ! Register 2-D energy input into internal tides for each frequency + write(var_name, '("TKE_itidal_input_freq",i1)') fr + write(var_descript, '("a fraction of which goes into rays in frequency ",i1)') fr + + CS%id_TKE_itidal_input(fr) = register_diag_field('ocean_model', var_name, diag%axesT1, & + Time, 'Conversion from barotropic to baroclinic tide, '//& + var_descript, 'W m-2', conversion=US%RZ3_T3_to_W_m2) + enddo ! Register 2-D energy losses (summed over angles, freq, modes) CS%id_tot_leak_loss = register_diag_field('ocean_model', 'ITide_tot_leak_loss', diag%axesT1, & Time, 'Internal tide energy loss to background drag', & @@ -2889,6 +2959,10 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) allocate(CS%id_En_mode(CS%nFreq,CS%nMode), source=-1) allocate(CS%id_En_ang_mode(CS%nFreq,CS%nMode), source=-1) allocate(CS%id_itidal_loss_mode(CS%nFreq,CS%nMode), source=-1) + allocate(CS%id_leak_loss_mode(CS%nFreq,CS%nMode), source=-1) + allocate(CS%id_quad_loss_mode(CS%nFreq,CS%nMode), source=-1) + allocate(CS%id_Froude_loss_mode(CS%nFreq,CS%nMode), source=-1) + allocate(CS%id_residual_loss_mode(CS%nFreq,CS%nMode), source=-1) allocate(CS%id_allprocesses_loss_mode(CS%nFreq,CS%nMode), source=-1) allocate(CS%id_itidal_loss_ang_mode(CS%nFreq,CS%nMode), source=-1) allocate(CS%id_Ub_mode(CS%nFreq,CS%nMode), source=-1) @@ -2929,6 +3003,30 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) CS%id_itidal_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, & diag%axesT1, Time, var_descript, 'W m-2', conversion=US%RZ3_T3_to_W_m2) call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + ! Leakage loss + write(var_name, '("Itide_leak_loss_freq",i1,"_mode",i1)') fr, m + write(var_descript, '("Internal tide energy loss due to leakage from frequency ",i1," mode ",i1)') fr, m + CS%id_leak_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, & + diag%axesT1, Time, var_descript, 'W m-2', conversion=US%RZ3_T3_to_W_m2) + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + ! Quad loss + write(var_name, '("Itide_quad_loss_freq",i1,"_mode",i1)') fr, m + write(var_descript, '("Internal tide energy quad loss from frequency ",i1," mode ",i1)') fr, m + CS%id_quad_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, & + diag%axesT1, Time, var_descript, 'W m-2', conversion=US%RZ3_T3_to_W_m2) + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + ! Froude loss + write(var_name, '("Itide_froude_loss_freq",i1,"_mode",i1)') fr, m + write(var_descript, '("Internal tide energy Froude loss from frequency ",i1," mode ",i1)') fr, m + CS%id_froude_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, & + diag%axesT1, Time, var_descript, 'W m-2', conversion=US%RZ3_T3_to_W_m2) + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) + ! residual losses + write(var_name, '("Itide_residual_loss_freq",i1,"_mode",i1)') fr, m + write(var_descript, '("Internal tide energy residual loss from frequency ",i1," mode ",i1)') fr, m + CS%id_residual_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, & + diag%axesT1, Time, var_descript, 'W m-2', conversion=US%RZ3_T3_to_W_m2) + call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) ! all loss processes write(var_name, '("Itide_allprocesses_loss_freq",i1,"_mode",i1)') fr, m write(var_descript, '("Internal tide energy loss due to all processes from frequency ",i1," mode ",i1)') fr, m @@ -2958,6 +3056,12 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) 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) + ! Register 2-D drag scale used for quadratic bottom drag for each frequency and mode + write(var_name, '("ITide_drag_freq",i1,"_mode",i1)') fr, m + write(var_descript, '("Interior and bottom drag int tide decay timescale in frequency ",i1, " mode ",i1)') fr, m + + CS%id_itide_drag(fr,m) = register_diag_field('ocean_model', var_name, diag%axesT1, Time, & + 's-1', conversion=US%s_to_T) enddo ; enddo diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 5b89c8c726..097628c032 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -391,8 +391,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call set_int_tide_input(u, v, h, tv, fluxes, CS%int_tide_input, dt, G, GV, US, & 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, CS%int_tide_input%Rho_bot, dt, G, GV, US, CS%int_tide_CSp) + call propagate_int_tide(h, tv, CS%int_tide_input%Nb, CS%int_tide_input%Rho_bot, dt, & + G, GV, US, CS%int_tide_input_CSp, CS%int_tide_CSp) + 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 3da21b48fb..7280106125 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -11,11 +11,13 @@ module MOM_int_tide_input use MOM_debugging, only : hchksum use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_file_parser, only : read_param 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_string_functions, only : extractWord use MOM_time_manager, only : time_type, set_time, operator(+), operator(<=) use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d @@ -27,6 +29,7 @@ module MOM_int_tide_input #include public set_int_tide_input, int_tide_input_init, int_tide_input_end +public get_input_TKE, get_barotropic_tidal_vel ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -44,9 +47,13 @@ module MOM_int_tide_input 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] - real, allocatable, dimension(:,:) :: TKE_itidal_coef + real, allocatable, dimension(:,:,:) :: TKE_itidal_coef !< 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]. + real, allocatable, dimension(:,:,:) :: & + TKE_itidal_input, & !< The internal tide TKE input at the bottom of the ocean [R Z3 T-3 ~> W m-2]. + tideamp !< The amplitude of the tidal velocities [Z T-1 ~> m s-1]. + character(len=200) :: inputdir !< The directory for input files. logical :: int_tide_source_test !< If true, apply an arbitrary generation site @@ -59,19 +66,19 @@ module MOM_int_tide_input integer :: int_tide_source_i !< I Location of generation site integer :: int_tide_source_j !< J Location of generation site logical :: int_tide_use_glob_ij !< Use global indices for generation site + integer :: nFreq = 0 !< The number of internal tide frequency bands !>@{ Diagnostic IDs - integer :: id_TKE_itidal_itide = -1, id_Nb = -1, id_N2_bot = -1 + integer, allocatable, dimension(:) :: id_TKE_itidal_itide + integer :: id_Nb = -1, id_N2_bot = -1 !>@} end type int_tide_input_CS !> This type is used to exchange fields related to the internal tides. type, public :: int_tide_input_type real, allocatable, dimension(:,:) :: & - 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]. Rho_bot !< The bottom density or the Boussinesq reference density [R ~> kg m-3]. end type int_tide_input_type @@ -110,6 +117,7 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) integer :: i, j, is, ie, js, je, nz, isd, ied, jsd, jed integer :: i_global, j_global + integer :: fr is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -133,52 +141,55 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) if (GV%Boussinesq .or. GV%semi_Boussinesq) then !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie + 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%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 + CS%TKE_itidal_input(i,j,fr) = min(GV%Z_to_H*CS%TKE_itidal_coef(i,j,fr)*itide%Nb(i,j), CS%TKE_itide_max) + enddo ; enddo ; enddo else !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie + 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) - 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_itidal_input(i,j,fr) = min((GV%RZ_to_H*Rho_bot(i,j)) * CS%TKE_itidal_coef(i,j,fr)*itide%Nb(i,j), & CS%TKE_itide_max) - enddo ; enddo + enddo ; enddo ; enddo endif if (CS%int_tide_source_test) then - itide%TKE_itidal_input(:,:) = 0.0 + CS%TKE_itidal_input(:,:,:) = 0.0 if (time_end <= CS%time_max_source) then if (CS%int_tide_use_glob_ij) then - do j=js,je ; do i=is,ie + do fr=1,CS%nFreq ; do j=js,je ; do i=is,ie i_global = i + G%idg_offset j_global = j + G%jdg_offset if ((i_global == CS%int_tide_source_i) .and. (j_global == CS%int_tide_source_j)) then - itide%TKE_itidal_input(i,j) = 1.0*US%kg_m3_to_R*US%m_to_Z**3*US%T_to_s**3 + CS%TKE_itidal_input(i,j,fr) = 1.0*US%kg_m3_to_R*US%m_to_Z**3*US%T_to_s**3 endif - enddo ; enddo + enddo ; enddo ; enddo else - do j=js,je ; do i=is,ie + do fr=1,CS%nFreq ; do j=js,je ; do i=is,ie ! Input an arbitrary energy point source.id_ if (((G%geoLonCu(I-1,j)-CS%int_tide_source_x) * (G%geoLonBu(I,j)-CS%int_tide_source_x) <= 0.0) .and. & ((G%geoLatCv(i,J-1)-CS%int_tide_source_y) * (G%geoLatCv(i,j)-CS%int_tide_source_y) <= 0.0)) then - itide%TKE_itidal_input(i,j) = 1.0*US%kg_m3_to_R*US%m_to_Z**3*US%T_to_s**3 + CS%TKE_itidal_input(i,j,fr) = 1.0*US%kg_m3_to_R*US%m_to_Z**3*US%T_to_s**3 endif - enddo ; enddo + enddo ; enddo ; enddo endif endif endif if (CS%debug) then call hchksum(N2_bot,"N2_bot",G%HI,haloshift=0, scale=US%s_to_T**2) - call hchksum(itide%TKE_itidal_input,"TKE_itidal_input",G%HI,haloshift=0, & + call hchksum(CS%TKE_itidal_input,"TKE_itidal_input",G%HI,haloshift=0, & scale=US%RZ3_T3_to_W_m2) endif call enable_averages(dt, time_end, CS%diag) - if (CS%id_TKE_itidal_itide > 0) call post_data(CS%id_TKE_itidal_itide, itide%TKE_itidal_input, CS%diag) + do fr=1,CS%nFreq + if (CS%id_TKE_itidal_itide(fr) > 0) call post_data(CS%id_TKE_itidal_itide(fr), & + CS%TKE_itidal_input(isd:ied,jsd:jed,fr), CS%diag) + enddo if (CS%id_Nb > 0) call post_data(CS%id_Nb, itide%Nb, CS%diag) if (CS%id_N2_bot > 0 ) call post_data(CS%id_N2_bot, N2_bot, CS%diag) @@ -319,6 +330,38 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot, rho_bo end subroutine find_N2_bottom +!> Returns TKE_itidal_input +subroutine get_input_TKE(G, TKE_itidal_input, nFreq, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). + real, dimension(SZI_(G),SZJ_(G),nFreq), & + intent(out) :: TKE_itidal_input !< The energy input to the internal waves [R Z3 T-3 ~> W m-2]. + integer, intent(in) :: nFreq !< number of frequencies + type(int_tide_input_CS), target :: CS !< A pointer that is set to point to the control + !! structure for the internal tide input module. + integer :: i,j,fr + + do fr=1,nFreq ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + TKE_itidal_input(i,j,fr) = CS%TKE_itidal_input(i,j,fr) + enddo ; enddo ; enddo + +end subroutine get_input_TKE + +!> Returns barotropic tidal velocities +subroutine get_barotropic_tidal_vel(G, vel_btTide, nFreq, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). + real, dimension(SZI_(G),SZJ_(G),nFreq), & + intent(out) :: vel_btTide !< Barotropic velocity read from file [L T-1 ~> m s-1]. + integer, intent(in) :: nFreq !< number of frequencies + type(int_tide_input_CS), target :: CS !< A pointer that is set to point to the control + !! structure for the internal tide input module. + integer :: i,j,fr + + do fr=1,nFreq ; do j=G%jsd,G%jed ; do i=G%isd,G%ied + vel_btTide(i,j,fr) = CS%tideamp(i,j,fr) + enddo ; enddo ; enddo + +end subroutine get_barotropic_tidal_vel + !> Initializes the data related to the internal tide input module subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) type(time_type), intent(in) :: Time !< The current model time @@ -337,6 +380,9 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) character(len=40) :: mdl = "MOM_int_tide_input" ! This module's name. character(len=200) :: filename, tideamp_file, h2_file ! Input file names or paths character(len=80) :: tideamp_var, rough_var ! Input file variable names + character(len=80) :: var_name + character(len=200) :: var_descript + character(len=200) :: tidefile_varnames real :: mask_itidal ! A multiplicative land mask, 0 or 1 [nondim] real :: max_frac_rough ! The fraction relating the maximum topographic roughness @@ -349,6 +395,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) integer :: tlen_days !< Time interval from start for adding wave source !! for testing internal tides (BDM) integer :: i, j, is, ie, js, je, isd, ied, jsd, jed + integer :: num_freq, fr if (associated(CS)) then call MOM_error(WARNING, "int_tide_input_init called with an associated "// & @@ -390,12 +437,15 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & units="m s-1", default=0.0, scale=US%m_s_to_L_T) + call read_param(param_file, "INTERNAL_TIDE_FREQS", num_freq) + CS%nFreq= num_freq + 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) - allocate(CS%TKE_itidal_coef(isd:ied,jsd:jed), source=0.0) + allocate(CS%TKE_itidal_input(isd:ied,jsd:jed,num_freq), source=0.0) + allocate(CS%tideamp(isd:ied,jsd:jed,num_freq), source=utide) + allocate(CS%TKE_itidal_coef(isd:ied,jsd:jed, num_freq), source=0.0) call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, & "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//& @@ -419,10 +469,13 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc") filename = trim(CS%inputdir) // trim(tideamp_file) call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename) - call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, & - "The name of the tidal amplitude variable in the input file.", & - default="tideamp") - call MOM_read_data(filename, tideamp_var, itide%tideamp, G%domain, scale=US%m_s_to_L_T) + + call read_param(param_file, "INTTIDE_AMP_VARNAMES", tidefile_varnames) + do fr=1,num_freq + tideamp_var = extractWord(tidefile_varnames,fr) + call MOM_read_data(filename, tideamp_var, CS%tideamp(:,:,fr), G%domain, scale=US%m_s_to_L_T) + enddo + endif call get_param(param_file, mdl, "H2_FILE", h2_file, & @@ -475,25 +528,31 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) endif endif - do j=js,je ; do i=is,ie + do fr=1,num_freq ; do j=js,je ; do i=is,ie mask_itidal = 1.0 if (G%bathyT(i,j) + G%Z_ref < min_zbot_itides) mask_itidal = 0.0 - itide%tideamp(i,j) = itide%tideamp(i,j) * mask_itidal * G%mask2dT(i,j) + CS%tideamp(i,j,fr) = CS%tideamp(i,j,fr) * mask_itidal * G%mask2dT(i,j) ! Restrict rms topo to a fraction (often 10 percent) of the column depth. 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 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 + CS%TKE_itidal_coef(i,j,fr) = 0.5*US%L_to_Z*kappa_h2_factor * GV%H_to_RZ * & + kappa_itides * itide%h2(i,j) * CS%tideamp(i,j,fr)**2 + enddo ; enddo ; enddo - CS%id_TKE_itidal_itide = register_diag_field('ocean_model','TKE_itidal_itide',diag%axesT1,Time, & - 'Internal Tide Driven Turbulent Kinetic Energy', & - 'W m-2', conversion=US%RZ3_T3_to_W_m2) + allocate( CS%id_TKE_itidal_itide(num_freq), source=-1) + + do fr=1,num_freq + write(var_name, '("TKE_itidal_itide_freq",i1)') fr + 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) + enddo CS%id_Nb = register_diag_field('ocean_model','Nb_itide',diag%axesT1,Time, & 'Bottom Buoyancy Frequency', 's-1', conversion=US%s_to_T)