Skip to content

Commit

Permalink
Streaming filter
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@c2xu
c2xu committed Dec 6, 2024
1 parent 2b72682 commit d0bb956
Show file tree
Hide file tree
Showing 2 changed files with 113 additions and 114 deletions.
85 changes: 21 additions & 64 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@ module MOM_barotropic
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
Expand Down Expand Up @@ -250,10 +249,9 @@ 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.
integer :: n_filters !< Number of streaming band-pass filters to be used 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
Expand Down Expand Up @@ -297,10 +295,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 on u-grid
Filt_CS_v !< Control structures for the streaming band-pass filter on v-grid
logical :: module_is_initialized = .false. !< If true, module has been initialized

integer :: isdw !< The lower i-memory limit for the wide halo arrays.
Expand Down Expand Up @@ -608,8 +604,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]
Expand Down Expand Up @@ -1598,15 +1594,9 @@ 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)
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.
Expand Down Expand Up @@ -5283,32 +5273,17 @@ 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.)
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
Expand Down Expand Up @@ -5338,24 +5313,6 @@ 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
Expand Down
142 changes: 92 additions & 50 deletions src/parameterizations/lateral/MOM_streaming_filter.F90
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,21 @@
!> Streaming band-pass filter for detecting the instantaneous tidal signals in the simulation
!!
!! Dec 5, 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.

module MOM_streaming_filter

use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL
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_tidal_forcing, only : tidal_frequency
use MOM_time_manager, only : time_type, time_type_to_real
use MOM_unit_scaling, only : unit_scale_type

Expand All @@ -12,108 +25,137 @@ module MOM_streaming_filter

#include <MOM_memory.h>

!> The control structure for storing the filter infomation of a particular field
!> Control structure for the MOM_streaming_filter module
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]
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=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]
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
!> 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

! 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")
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

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
CS%nf = nf

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
CS%is = HI%isd ; CS%ie = HI%ied ; CS%js = HI%jsd ; CS%je = HI%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
CS%is = HI%IsdB ; CS%ie = HI%IedB ; CS%js = HI%jsd ; CS%je = HI%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
CS%is = HI%isd ; CS%ie = HI%ied ; CS%js = HI%JsdB ; CS%je = HI%JedB
case default
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

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))
read(filter_name_str, *) CS%filter_names

do c=1,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")
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)), &
" filter registered, target frequency = ", CS%filter_omega(c), &
", bandwidth = ", CS%filter_alpha(c)
call MOM_error(NOTE, trim(mesg))
enddo

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(:,:,:), 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
integer :: i, j, k, 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(:,:)
do k=1,CS%nf
CS%u1(:,:,k) = u(:,:)
enddo
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
do k=1,CS%nf
c1 = CS%filter_omega(k) * dt
c2 = 1.0 - CS%filter_alpha(k) * c1

do j=js,je ; do i=is,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
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.
!! 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.
!!
!! 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.
!! streaming band-pass filters. Journal of Advances in Modeling Earth Systems, 16, e2024MS004319.
!! https://doi.org/10.1029/2024MS004319

end module MOM_streaming_filter

0 comments on commit d0bb956

Please sign in to comment.