Streaming filter

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
+