Streaming filter

The restarting capability has been implemented.

src/core/MOM_barotropic.F90

@@ -26,7 +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_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
@@ -4974,6 +4974,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
@@ -5259,8 +5265,6 @@ 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
@@ -5273,18 +5277,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, "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", CS%n_filters, &
-                 "Number of streaming band-pass filters to be used in the simulation.", &
-                 default=0, do_not_log=.not.CS%use_filter)
-  if (CS%n_filters==0) CS%use_filter = .false.
-  if (CS%use_filter) then
-    call Filt_register(CS%n_filters, 'u', HI, US, param_file, CS%Filt_CS_u)
-    call Filt_register(CS%n_filters, 'v', HI, US, param_file, CS%Filt_CS_v)
-  endif
-
   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
@@ -5313,6 +5305,19 @@ 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 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", CS%n_filters, &
+                 "Number of streaming band-pass filters to be used in the simulation.", &
+                 default=0, do_not_log=.not.CS%use_filter)
+  if (CS%n_filters==0) CS%use_filter = .false.
+  if (CS%use_filter) then
+    call Filt_register(CS%n_filters, 'ubt', 'u', HI, CS%Filt_CS_u, restart_CS)
+    call Filt_register(CS%n_filters, 'vbt', 'v', HI, CS%Filt_CS_v, restart_CS)
+  endif
+
 end subroutine register_barotropic_restarts
 
 !> \namespace mom_barotropic

src/parameterizations/lateral/MOM_streaming_filter.F90

@@ -1,62 +1,67 @@
 !> Streaming band-pass filter for detecting the instantaneous tidal signals in the simulation
 !!
-!! Dec 5, 2024: Major revision.
+!! Dec 9, 2024: Major revision.
 !!
-!! The filters and their target frequencies are no longer hard-coded. Instead, up to 10 filters
-!! with tidal frequencies as their target frequencies and an unspecified number of filters with
-!! arbitrary target frequencies can be turned on. The filter names are specified in MOM_input
-!! and must consist of two letters/numbers. If a filter name 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 provide the target frequency. In either case, the target frequency
-!! is specified by "TIDE_${FILTER_NAME}_FREQ" in MOM_input.
+!! The filters and their target frequencies are no longer hard-coded. Instead, multiple filters
+!! with tidal frequencies 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. In either case, the target frequency is specified by "TIDE_${FILTER_NAME}_FREQ".
+!!
+!! 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.
 
 module MOM_streaming_filter
 
 use MOM_error_handler, only : MOM_mesg, MOM_error, NOTE, FATAL
 use MOM_file_parser,   only : get_param, param_file_type
 use MOM_hor_index,     only : hor_index_type
+use MOM_io,            only : axis_info, set_axis_info
+use MOM_restart,       only : register_restart_field, query_initialized, MOM_restart_CS
 use MOM_tidal_forcing, only : tidal_frequency
 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
+public Filt_register, Filt_init, Filt_accum
 
 #include <MOM_memory.h>
 
 !> Control structure for the MOM_streaming_filter module
 type, public :: Filter_CS ; private
-  integer :: nf                       !< Number of filters to be used in the simulation
+  integer :: nf                        !< Number of filters to be used in the simulation
   !>@{ Lower and upper bounds of input data
   integer :: is, ie, js, je
   !>@}
+  character(len=8) :: key              !< Identifier of the variable to be filtered
   character(len=2), allocatable, dimension(:) :: filter_names !< Names of filters
-  real, allocatable, dimension(:)     :: filter_omega !< Target frequencies of filters [T-1 ~> s-1]
-  real, allocatable, dimension(:)     :: filter_alpha !< Bandwidth parameters of filters [nondim]
-  real, allocatable, dimension(:,:,:) :: s1, &        !< Dummy variable [A]
-                                         u1           !< Filtered data [A]
-  real :: old_time = -1.0             !< The time of the previous accumulating step [T ~> s]
+  real, allocatable, dimension(:)      :: filter_omega !< Target frequencies of filters [rad T-1 ~> rad s-1]
+  real, allocatable, dimension(:)      :: filter_alpha !< Bandwidth parameters of filters [nondim]
+  real, allocatable, dimension(:,:,:)  :: s1, &        !< Dummy variable [A]
+                                          u1           !< Filtered data [A]
+  real :: old_time = -1.0              !< The time of the previous accumulating step [T ~> s]
 end type Filter_CS
 
 contains
 
-!> This subroutine initializes the filters given the number of filters and the grid
-subroutine Filt_register(nf, grid, HI, US, param_file, CS)
-  integer,               intent(in)   :: nf           !< Number of filters to be used in the simulation
-  character(len=*),      intent(in)   :: grid         !< Horizontal grid location: h, u, or v
-  type(hor_index_type),  intent(in)   :: HI           !< Horizontal index type structure
-  type(unit_scale_type), intent(in)   :: US           !< A dimensional unit scaling type
-  type(param_file_type), intent(in)   :: param_file   !< A structure to parse for run-time parameters
-  type(Filter_CS),       intent(out)  :: CS           !< Control structure for MOM_streaming_filter
+!> This subroutine registers the filter variables given the number of filters and the grid
+subroutine Filt_register(nf, key, grid, HI, CS, restart_CS)
+  integer,               intent(in)    :: nf           !< Number of filters to be used in the simulation
+  character(len=*),      intent(in)    :: key          !< Identifier of the variable to be filtered
+  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 of MOM_streaming_filter
+  type(MOM_restart_CS),  intent(inout) :: restart_CS   !< MOM restart control structure
 
   ! Local variables
-  character(len=40)  :: mdl = "MOM_streaming_filter"  !< This module's name
-  character(len=50)  :: filter_name_str               !< List of filters to be registered
-  character(len=200) :: mesg
+  type(axis_info) :: filter_axis(1)
+  real, dimension(:), allocatable :: n_filters         !< Labels of filters [nondim]
   integer :: c
 
-  CS%nf = nf
+  CS%nf  = nf
+  CS%key = key
 
   select case (trim(grid))
     case ('h')
@@ -69,29 +74,57 @@ subroutine Filt_register(nf, grid, HI, US, param_file, CS)
       call MOM_error(FATAL, "MOM_streaming_filter: horizontal grid not supported")
   end select
 
-  allocate(CS%s1(CS%is:CS%ie,CS%js:CS%je,nf)) ; CS%s1(:,:,:) = 0.0
-  allocate(CS%u1(CS%is:CS%ie,CS%js:CS%je,nf)) ; CS%u1(:,:,:) = 0.0
+  allocate(CS%s1(CS%is:CS%ie, CS%js:CS%je, nf)) ; CS%s1(:,:,:) = 0.0
+  allocate(CS%u1(CS%is:CS%ie, CS%js:CS%je, nf)) ; CS%u1(:,:,:) = 0.0
+
+  ! Register restarts for s1 and u1
+  allocate(n_filters(nf))
+
+  do c=1,nf ; n_filters(c) = c ; enddo
+
+  call set_axis_info(filter_axis(1), "n_filters", "", "number of filters", nf, n_filters, "N", 1)
+
+  call register_restart_field(CS%s1(:,:,:), "Filter_"//trim(key)//"_s1", .false., restart_CS, &
+                              longname="Dummy variable for streaming band-pass filter", &
+                              hor_grid=trim(grid), z_grid="1", t_grid="s", extra_axes=filter_axis)
+  call register_restart_field(CS%u1(:,:,:), "Filter_"//trim(key)//"_u1", .false., restart_CS, &
+                              longname="Output of streaming band-pass filter", &
+                              hor_grid=trim(grid), z_grid="1", t_grid="s", extra_axes=filter_axis)
+
+end subroutine Filt_register
+
+!> This subroutine initializes the filters
+subroutine Filt_init(param_file, US, CS, restart_CS)
+  type(param_file_type), intent(in)    :: param_file   !< A structure to parse for run-time parameters
+  type(unit_scale_type), intent(in)    :: US           !< A dimensional unit scaling type
+  type(Filter_CS),       intent(inout) :: CS           !< Control structure of MOM_streaming_filter
+  type(MOM_restart_CS),  intent(in)    :: restart_CS   !< MOM restart control structure
+
+  ! Local variables
+  character(len=40)  :: mdl = "MOM_streaming_filter"   !< This module's name
+  character(len=50)  :: filter_name_str                !< List of filters to be registered
+  character(len=200) :: mesg
+  integer :: c
 
   call get_param(param_file, mdl, "FILTER_NAMES", filter_name_str, &
                  "Names of streaming band-pass filters to be used in the simulation.", &
                  fail_if_missing=.true.)
-  allocate(CS%filter_names(nf))
-  allocate(CS%filter_omega(nf))
-  allocate(CS%filter_alpha(nf))
+  allocate(CS%filter_names(CS%nf))
+  allocate(CS%filter_omega(CS%nf))
+  allocate(CS%filter_alpha(CS%nf))
   read(filter_name_str, *) CS%filter_names
 
-  do c=1,nf
+  do c=1,CS%nf
     ! If filter_name_str consists of tidal constituents, use tidal frequencies.
-    call get_param(param_file, mdl, "TIDE_"//trim(CS%filter_names(c))//"_FREQ", CS%filter_omega(c), &
-                   "Target frequency of the "//trim(CS%filter_names(c))//" filter. "//&
-                   "This is used if USE_FILTER is true and "//trim(CS%filter_names(c))//&
-                   " is in FILTER_NAMES, even if TIDES and TIDE_"//trim(CS%filter_names(c))//&
-                   " are false.", units="s-1", &
-                   default=tidal_frequency(trim(CS%filter_names(c))), scale=US%T_to_s)
-    call get_param(param_file, mdl, "FILTER_"//trim(CS%filter_names(c))//"_ALPHA", CS%filter_alpha(c), &
-                   "Bandwidth parameter of the "//trim(CS%filter_names(c))//" filter. "//&
-                   "Must be positive.", units="nondim", fail_if_missing=.true.)
-    if (CS%filter_omega(c)<=0.0) call MOM_error(FATAL, "MOM_streaming_filter: target frequency <= 0")
+    call get_param(param_file, mdl, "FILTER_"//trim(CS%filter_names(c))//"_OMEGA", &
+                   CS%filter_omega(c), "Target frequency of the "//trim(CS%filter_names(c))//&
+                   " filter. This is used if USE_FILTER is true and "//trim(CS%filter_names(c))//&
+                   " is in FILTER_NAMES.", units="rad s-1", scale=US%T_to_s, default=0.0)
+    call get_param(param_file, mdl, "FILTER_"//trim(CS%filter_names(c))//"_ALPHA", &
+                   CS%filter_alpha(c), "Bandwidth parameter of the "//trim(CS%filter_names(c))//&
+                   " filter. Must be positive.", units="nondim", fail_if_missing=.true.)
+
+    if (CS%filter_omega(c)<=0.0) CS%filter_omega(c) = tidal_frequency(trim(CS%filter_names(c)))
     if (CS%filter_alpha(c)<=0.0) call MOM_error(FATAL, "MOM_streaming_filter: bandwidth <= 0")
 
     write(mesg,*) "MOM_streaming_filter: ", trim(CS%filter_names(c)), &
@@ -100,34 +133,43 @@ subroutine Filt_register(nf, grid, HI, US, param_file, CS)
     call MOM_error(NOTE, trim(mesg))
   enddo
 
-end subroutine Filt_register
+  if (query_initialized(CS%s1, "Filter_"//trim(CS%key)//"_s1", restart_CS)) then
+    write(mesg,*) "MOM_streaming_filter: Dummy variable for filter ", trim(CS%key), &
+                  " found in restart files."
+  else
+    write(mesg,*) "MOM_streaming_filter: Dummy variable for filter ", trim(CS%key), &
+                  " not found in restart files. The filter will spin up from zeros."
+  endif
+  call MOM_error(NOTE, trim(mesg))
+
+  if (query_initialized(CS%u1, "Filter_"//trim(CS%key)//"_u1", restart_CS)) then
+    write(mesg,*) "MOM_streaming_filter: Output of filter ", trim(CS%key), &
+                  " found in restart files."
+  else
+    write(mesg,*) "MOM_streaming_filter: Output of filter ", trim(CS%key), &
+                  " not found in restart files. The filter will spin up from zeros."
+  endif
+  call MOM_error(NOTE, trim(mesg))
+
+end subroutine Filt_init
 
-!> 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).
+!> This subroutine timesteps the filter equations. Here, u is the broadband input signal from the model,
+!! and u1 is the filtered, narrowband output signal, obtained from the solution of 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
+  type(Filter_CS),        target,  intent(inout) :: CS   !< Control structure of MOM_streaming_filter
   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, k, is, ie, js, je
+  integer :: i, j, k
 
   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
-    do k=1,CS%nf
-      CS%u1(:,:,k) = u(:,:)
-    enddo
-  endif
-
+  if (CS%old_time<0.0) CS%old_time = now
   dt = now - CS%old_time
   CS%old_time = now
 
@@ -136,7 +178,7 @@ subroutine Filt_accum(u, u1, Time, US, CS)
     c1 = CS%filter_omega(k) * dt
     c2 = 1.0 - CS%filter_alpha(k) * c1
 
-    do j=js,je ; do i=is,ie
+    do j=CS%js,CS%je ; do i=CS%is,CS%ie
       CS%s1(i,j,k) =  c1 *  CS%u1(i,j,k) + CS%s1(i,j,k)
       CS%u1(i,j,k) = -c1 * (CS%s1(i,j,k) - CS%filter_alpha(k) * u(i,j)) + c2 * CS%u1(i,j,k)
     enddo; enddo
@@ -147,11 +189,8 @@ 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 (filter_omega) and the bandwidth parameter (filter_alpha) of the filter.
-!! At each timestep, the filter takes model output (u) as the input and returns a time series (u1) consisting of
-!! sinusoidal motions 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.
+!! The algorithm detects the instantaneous, narrowband tidal signals (u1) from the broadband model output
+!! (u) by solving a set of coupled ODEs (the filter equations) at each time step.
 !!
 !! 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, 16, e2024MS004319.