From d0bb95620956a89589c243c10395ad60c96aae93 Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 6 Dec 2024 14:45:21 -0400 Subject: [PATCH] Streaming filter 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. --- src/core/MOM_barotropic.F90 | 85 +++-------- .../lateral/MOM_streaming_filter.F90 | 142 ++++++++++++------ 2 files changed, 113 insertions(+), 114 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index af2beca1fb..9641eb294f 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -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 @@ -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 @@ -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. @@ -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] @@ -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. @@ -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 @@ -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 diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index a91f6661f2..3e8a932ef2 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -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 @@ -12,73 +25,97 @@ module MOM_streaming_filter #include -!> 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 @@ -86,20 +123,24 @@ subroutine Filt_accum(u, u1, Time, US, CS) ! 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 @@ -107,13 +148,14 @@ 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