Streaming Filter

The filters and their target frequencies are no longer hard-coded.
Instead, multiple filters with tidal or arbitrary frequencies as
their target frequencies can be turned on. The filter names are
specified in MOM_input and must consist of two letters/numbers. If
the name of a filter is the same as the name of a tidal constituent,
then the corresponding tidal frequency will be used as its target
frequency. Otherwise, the user must specify the target frequency by
"FILTER_${FILTER_NAME}_OMEGA" in MOM_input.

The restarting capability has also been implemented. Because the
filtering is a point-wise operation, all variables are considered
as fields, even if they are velocity components.

Author: c2xu

src/core/MOM_barotropic.F90

@@ -26,8 +26,7 @@ module MOM_barotropic
 use MOM_restart, only : query_initialized, MOM_restart_CS
 use MOM_self_attr_load, only : scalar_SAL_sensitivity
 use MOM_self_attr_load, only : SAL_CS
-use MOM_streaming_filter, only : Filt_register, Filt_accum, Filter_CS
-use MOM_tidal_forcing, only : tidal_frequency
+use MOM_streaming_filter, only : Filt_register, Filt_init, Filt_accum, Filter_CS
 use MOM_time_manager, only : time_type, real_to_time, operator(+), operator(-)
 use MOM_unit_scaling, only : unit_scale_type
 use MOM_variables, only : BT_cont_type, alloc_bt_cont_type
@@ -250,10 +249,8 @@ module MOM_barotropic
   logical :: linearized_BT_PV  !< If true, the PV and interface thicknesses used
                              !! in the barotropic Coriolis calculation is time
                              !! invariant and linearized.
-  logical :: use_filter_m2   !< If true, apply streaming band-pass filter for detecting
-                             !! instantaneous tidal signals.
-  logical :: use_filter_k1   !< If true, apply streaming band-pass filter for detecting
-                             !! instantaneous tidal signals.
+  logical :: use_filter      !< If true, use streaming band-pass filter to detect the
+                             !! instantaneous tidal signals in the simulation.
   logical :: use_wide_halos  !< If true, use wide halos and march in during the
                              !! barotropic time stepping for efficiency.
   logical :: clip_velocity   !< If true, limit any velocity components that are
@@ -297,10 +294,8 @@ module MOM_barotropic
   type(hor_index_type), pointer :: debug_BT_HI => NULL() !< debugging copy of horizontal index_type
   type(SAL_CS), pointer :: SAL_CSp => NULL() !< Control structure for SAL
   type(harmonic_analysis_CS), pointer :: HA_CSp => NULL() !< Control structure for harmonic analysis
-  type(Filter_CS) :: Filt_CS_um2, & !< Control structures for the M2 streaming filter
-                     Filt_CS_vm2, & !< Control structures for the M2 streaming filter
-                     Filt_CS_uk1, & !< Control structures for the K1 streaming filter
-                     Filt_CS_vk1    !< Control structures for the K1 streaming filter
+  type(Filter_CS) :: Filt_CS_u, & !< Control structures for the streaming band-pass filter of ubt
+                     Filt_CS_v    !< Control structures for the streaming band-pass filter of vbt
   logical :: module_is_initialized = .false.  !< If true, module has been initialized
 
   integer :: isdw !< The lower i-memory limit for the wide halo arrays.
@@ -608,8 +603,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     DCor_v, &     ! An averaged total thickness at v points [H ~> m or kg m-2].
     Datv          ! Basin depth at v-velocity grid points times the x-grid
                   ! spacing [H L ~> m2 or kg m-1].
-  real, dimension(:,:), pointer :: um2, uk1, vm2, vk1
-                  ! M2 and K1 velocities from the output of streaming filters [m s-1]
+  real, dimension(:,:,:), pointer :: ufilt, vfilt
+                  ! Filtered velocities from the output of streaming filters [m s-1]
   real, target, dimension(SZIW_(CS),SZJW_(CS)) :: &
     eta, &        ! The barotropic free surface height anomaly or column mass
                   ! anomaly [H ~> m or kg m-2]
@@ -1598,15 +1593,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     endif ; enddo ; enddo
   endif
 
-  ! Here is an example of how the filter equations are time stepped to determine the M2 and K1 velocities.
-  ! The filters are initialized and registered in subroutine barotropic_init.
-  if (CS%use_filter_m2) then
-    call Filt_accum(ubt, um2, CS%Time, US, CS%Filt_CS_um2)
-    call Filt_accum(vbt, vm2, CS%Time, US, CS%Filt_CS_vm2)
-  endif
-  if (CS%use_filter_k1) then
-    call Filt_accum(ubt, uk1, CS%Time, US, CS%Filt_CS_uk1)
-    call Filt_accum(vbt, vk1, CS%Time, US, CS%Filt_CS_vk1)
+  ! Note that the filtered velocities are only updated during the current predictor step,
+  ! and are calculated using the barotropic velocity from the previous correction step.
+  if (CS%use_filter) then
+    call Filt_accum(ubt, ufilt, CS%Time, US, CS%Filt_CS_u)
+    call Filt_accum(vbt, vfilt, CS%Time, US, CS%Filt_CS_v)
   endif
 
   ! Zero out the arrays for various time-averaged quantities.
@@ -4979,6 +4970,12 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
     endif ! len_trim(wave_drag_file) > 0
   endif ! CS%linear_wave_drag
 
+  ! Initialize streaming band-pass filters
+  if (CS%use_filter) then
+    call Filt_init(param_file, US, CS%Filt_CS_u, restart_CS)
+    call Filt_init(param_file, US, CS%Filt_CS_v, restart_CS)
+  endif
+
   CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input
 
   dtbt_tmp = -1.0
@@ -5263,9 +5260,8 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS)
   ! Local variables
   type(vardesc) :: vd(3)
   character(len=40)  :: mdl = "MOM_barotropic"  ! This module's name.
+  integer :: n_filters       !< Number of streaming band-pass filters to be used in the simulation.
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
-  real :: am2, ak1      !< Bandwidth parameters of the M2 and K1 streaming filters [nondim]
-  real :: om2, ok1      !< Target frequencies of the M2 and K1 streaming filters [rad T-1 ~> rad s-1]
 
   isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed
   IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB
@@ -5278,33 +5274,6 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS)
                  "sum(u dh_dt) while also correcting for truncation errors.", &
                  default=.false., do_not_log=.true.)
 
-  call get_param(param_file, mdl, "STREAMING_FILTER_M2", CS%use_filter_m2, &
-                 "If true, turn on streaming band-pass filter for detecting "//&
-                 "instantaneous tidal signals.", default=.false.)
-  call get_param(param_file, mdl, "STREAMING_FILTER_K1", CS%use_filter_k1, &
-                 "If true, turn on streaming band-pass filter for detecting "//&
-                 "instantaneous tidal signals.", default=.false.)
-  call get_param(param_file, mdl, "FILTER_ALPHA_M2", am2, &
-                 "Bandwidth parameter of the streaming filter targeting the M2 frequency. "//&
-                 "Must be positive. To turn off filtering, set FILTER_ALPHA_M2 <= 0.0.", &
-                 default=0.0, units="nondim", do_not_log=.not.CS%use_filter_m2)
-  call get_param(param_file, mdl, "FILTER_ALPHA_K1", ak1, &
-                 "Bandwidth parameter of the streaming filter targeting the K1 frequency. "//&
-                 "Must be positive. To turn off filtering, set FILTER_ALPHA_K1 <= 0.0.", &
-                 default=0.0, units="nondim", do_not_log=.not.CS%use_filter_k1)
-  call get_param(param_file, mdl, "TIDE_M2_FREQ", om2, &
-                 "Frequency of the M2 tidal constituent. "//&
-                 "This is only used if TIDES and TIDE_M2"// &
-                 " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and M2"// &
-                 " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=tidal_frequency("M2"), &
-                 scale=US%T_to_s, do_not_log=.true.)
-  call get_param(param_file, mdl, "TIDE_K1_FREQ", ok1, &
-                 "Frequency of the K1 tidal constituent. "//&
-                 "This is only used if TIDES and TIDE_K1"// &
-                 " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and K1"// &
-                 " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=tidal_frequency("K1"), &
-                 scale=US%T_to_s, do_not_log=.true.)
-
   ALLOC_(CS%ubtav(IsdB:IedB,jsd:jed))      ; CS%ubtav(:,:) = 0.0
   ALLOC_(CS%vbtav(isd:ied,JsdB:JedB))      ; CS%vbtav(:,:) = 0.0
   if (CS%gradual_BT_ICs) then
@@ -5333,22 +5302,17 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS)
   call register_restart_field(CS%dtbt, "DTBT", .false., restart_CS, &
                               longname="Barotropic timestep", units="seconds", conversion=US%T_to_s)
 
-  ! Initialize and register streaming filters
-  if (CS%use_filter_m2) then
-    if (am2 > 0.0 .and. om2 > 0.0) then
-      call Filt_register(am2, om2, 'u', HI, CS%Filt_CS_um2)
-      call Filt_register(am2, om2, 'v', HI, CS%Filt_CS_vm2)
-    else
-      CS%use_filter_m2 = .false.
-    endif
-  endif
-  if (CS%use_filter_k1) then
-    if (ak1 > 0.0 .and. ok1 > 0.0) then
-      call Filt_register(ak1, ok1, 'u', HI, CS%Filt_CS_uk1)
-      call Filt_register(ak1, ok1, 'v', HI, CS%Filt_CS_vk1)
-    else
-      CS%use_filter_k1 = .false.
-    endif
+  ! Initialize and register streaming band-pass filters
+  call get_param(param_file, mdl, "USE_FILTER", CS%use_filter, &
+                 "If true, use streaming band-pass filters to detect the "//&
+                 "instantaneous tidal signals in the simulation.", default=.false.)
+  call get_param(param_file, mdl, "N_FILTERS", n_filters, &
+                 "Number of streaming band-pass filters to be used in the simulation.", &
+                 default=0, do_not_log=.not.CS%use_filter)
+  if (n_filters<=0) CS%use_filter = .false.
+  if (CS%use_filter) then
+    call Filt_register(n_filters, 'ubt', 'u', HI, CS%Filt_CS_u, restart_CS)
+    call Filt_register(n_filters, 'vbt', 'v', HI, CS%Filt_CS_v, restart_CS)
   endif
 
 end subroutine register_barotropic_restarts