Streaming filter (NOAA-GFDL#675)

In the MOM_streaming_filter module, a system of two coupled ODEs is
configured as a streaming band-pass filter that takes the broadband
model output as the input and returns the filtered signal at a user
specified tidal frequency as the output. It is capable of detecting
the instantaneous tidal signals while the model is running, and can
be used to impose frequency-dependent wave drag or to de-tide model
output. An example of filtering the M2 and K1 barotropic velocities
is provided in the MOM_barotropic module.

Added runtime flags for the M2 and K1 filters. Simplified the code
by removing the control structure of the linked list, and now each
filter has its own control structure.

Using hor_index_type argument to avoid calculation in halo regions.

Author: c2xu

src/core/MOM_barotropic.F90

@@ -26,6 +26,8 @@ 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_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
@@ -248,6 +250,10 @@ 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_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
@@ -291,6 +297,10 @@ 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
   logical :: module_is_initialized = .false.  !< If true, module has been initialized
 
   integer :: isdw !< The lower i-memory limit for the wide halo arrays.
@@ -598,6 +608,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, target, dimension(SZIW_(CS),SZJW_(CS)) :: &
     eta, &        ! The barotropic free surface height anomaly or column mass
                   ! anomaly [H ~> m or kg m-2]
@@ -1586,6 +1598,17 @@ 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)
+  endif
+
   ! Zero out the arrays for various time-averaged quantities.
   if (find_etaav) then
     !$OMP do
@@ -5247,6 +5270,8 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS)
   type(vardesc) :: vd(3)
   character(len=40)  :: mdl = "MOM_barotropic"  ! This module's name.
   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 [T-1 ~> 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
@@ -5259,6 +5284,33 @@ 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="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="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
@@ -5287,6 +5339,24 @@ 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
+  endif
+
 end subroutine register_barotropic_restarts
 
 !> \namespace mom_barotropic

src/parameterizations/lateral/MOM_streaming_filter.F90

@@ -0,0 +1,119 @@
+!> Streaming band-pass filter for detecting the instantaneous tidal signals in the simulation
+module MOM_streaming_filter
+
+use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL
+use MOM_hor_index,     only : hor_index_type
+use MOM_time_manager,  only : time_type, time_type_to_real
+use MOM_unit_scaling,  only : unit_scale_type
+
+implicit none ; private
+
+public Filt_register, Filt_accum
+
+#include <MOM_memory.h>
+
+!> The control structure for storing the filter infomation of a particular field
+type, public :: Filter_CS ; private
+  real       :: a, &                   !< Parameter that determines the bandwidth [nondim]
+                om, &                  !< Target frequency of the filter [T-1 ~> s-1]
+                old_time = -1.0        !< The time of the previous accumulating step [T ~> s]
+  real, allocatable, dimension(:,:) :: s1, & !< Dummy variable [A]
+                                       u1    !< Filtered data [A]
+  !>@{ Lower and upper bounds of input data
+  integer :: is, ie, js, je
+  !>@}
+end type Filter_CS
+
+contains
+
+!> This subroutine registers each of the fields to be filtered.
+subroutine Filt_register(a, om, grid, HI, CS)
+  real,                      intent(in)    :: a      !< Parameter that determines the bandwidth [nondim]
+  real,                      intent(in)    :: om     !< Target frequency of the filter [T-1 ~> s-1]
+  character(len=*),          intent(in)    :: grid   !< Horizontal grid location: h, u, or v
+  type(hor_index_type),      intent(in)    :: HI     !< Horizontal index type structure
+  type(Filter_CS),           intent(out)   :: CS     !< Control structure for the current field
+
+  ! Local variables
+  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
+
+  if (a  <= 0.0) call MOM_error(FATAL, "MOM_streaming_filter: bandwidth <= 0")
+  if (om <= 0.0) call MOM_error(FATAL, "MOM_streaming_filter: target frequency <= 0")
+
+  CS%a  = a
+  CS%om = om
+
+  isd  = HI%isd  ; ied  = HI%ied  ; jsd  = HI%jsd  ; jed  = HI%jed
+  IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB
+
+  select case (trim(grid))
+    case ('h')
+      allocate(CS%s1(isd:ied,jsd:jed))   ; CS%s1(:,:) = 0.0
+      allocate(CS%u1(isd:ied,jsd:jed))   ; CS%u1(:,:) = 0.0
+      CS%is = isd  ; CS%ie = ied  ; CS%js = jsd  ; CS%je = jed
+    case ('u')
+      allocate(CS%s1(IsdB:IedB,jsd:jed)) ; CS%s1(:,:) = 0.0
+      allocate(CS%u1(IsdB:IedB,jsd:jed)) ; CS%u1(:,:) = 0.0
+      CS%is = IsdB ; CS%ie = IedB ; CS%js = jsd  ; CS%je = jed
+    case ('v')
+      allocate(CS%s1(isd:ied,JsdB:JedB)) ; CS%s1(:,:) = 0.0
+      allocate(CS%u1(isd:ied,JsdB:JedB)) ; CS%u1(:,:) = 0.0
+      CS%is = isd  ; CS%ie = ied  ; CS%js = JsdB ; CS%je = JedB
+    case default
+      call MOM_error(FATAL, "MOM_streaming_filter: horizontal grid not supported")
+  end select
+
+end subroutine Filt_register
+
+!> This subroutine timesteps the filter equations. It takes model output u at the current time step as the input,
+!! and returns tidal signal u1 as the output, which is the solution of a set of two ODEs (the filter equations).
+subroutine Filt_accum(u, u1, Time, US, CS)
+  real, dimension(:,:), pointer, intent(out)   :: u1   !< Output of the filter [A]
+  type(time_type),               intent(in)    :: Time !< The current model time
+  type(unit_scale_type),         intent(in)    :: US   !< A dimensional unit scaling type
+  type(Filter_CS),      target,  intent(inout) :: CS   !< Control structure of the MOM_streaming_filter module
+  real, dimension(CS%is:CS%ie,CS%js:CS%je), intent(in) :: u !< Input into the filter [A]
+
+  ! Local variables
+  real    :: now, &              !< The current model time [T ~> s]
+             dt, &               !< Time step size for the filter equations [T ~> s]
+             c1, c2              !< Coefficients for the filter equations [nondim]
+  integer :: i, j, is, ie, js, je
+
+  now = US%s_to_T * time_type_to_real(Time)
+  is = CS%is ; ie = CS%ie ; js = CS%js ; je = CS%je
+
+  ! Initialize u1
+  if (CS%old_time < 0.0) then
+    CS%old_time = now
+    CS%u1(:,:)  = u(:,:)
+  endif
+
+  dt = now - CS%old_time
+  CS%old_time = now
+
+  ! Timestepping
+  c1 = CS%om * dt
+  c2 = 1.0 - CS%a * c1
+
+  do j=js,je ; do i=is,ie
+    CS%s1(i,j) =  c1 *  CS%u1(i,j) + CS%s1(i,j)
+    CS%u1(i,j) = -c1 * (CS%s1(i,j) - CS%a * u(i,j)) + c2 * CS%u1(i,j)
+  enddo; enddo
+  u1 => CS%u1
+
+end subroutine Filt_accum
+
+!> \namespace streaming_filter
+!!
+!! This module detects instantaneous tidal signals in the model output using a set of coupled ODEs (the filter
+!! equations), given the target frequency (om) and the bandwidth parameter (a) of the filter. At each timestep,
+!! the filter takes model output (u) as the input and returns a time series consisting of sinusoidal motions (u1)
+!! near its target frequency. The filtered tidal signals can be used to parameterize frequency-dependent drag, or
+!! to detide the model output. See Xu & Zaron (2024) for detail.
+!!
+!! Reference: Xu, C., & Zaron, E. D. (2024). Detecting instantaneous tidal signals in ocean models utilizing
+!! streaming band-pass filters. Journal of Advances in Modeling Earth Systems. Under review.
+
+end module MOM_streaming_filter
+