From b204dc9b157e94c9931124ca100ffc64dd6909f6 Mon Sep 17 00:00:00 2001 From: William Xu Date: Thu, 27 Jun 2024 22:18:56 -0400 Subject: [PATCH 1/9] Streaming filter 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. --- src/core/MOM_barotropic.F90 | 31 +++ .../lateral/MOM_streaming_filter.F90 | 178 ++++++++++++++++++ 2 files changed, 209 insertions(+) create mode 100644 src/parameterizations/lateral/MOM_streaming_filter.F90 diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index e97e794cd6..2c2511e353 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -26,6 +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_init, Filt_register, Filt_accum, streaming_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 @@ -291,6 +292,7 @@ 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(streaming_filter_CS) :: Filt_CS !< Control structure for streaming filters 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 +600,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 +1590,13 @@ 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. + call Filt_accum('um2', ubt, um2, CS%Time, US, CS%Filt_CS) + call Filt_accum('uk1', ubt, uk1, CS%Time, US, CS%Filt_CS) + call Filt_accum('vm2', vbt, vm2, CS%Time, US, CS%Filt_CS) + call Filt_accum('vk1', vbt, vk1, CS%Time, US, CS%Filt_CS) + ! Zero out the arrays for various time-averaged quantities. if (find_etaav) then !$OMP do @@ -4454,6 +4465,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, real :: dtbt_tmp ! A temporary copy of CS%dtbt read from a restart file [T ~> s] real :: wave_drag_scale ! A scaling factor for the barotropic linear wave drag ! piston velocities [nondim]. + 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 [s-1] character(len=200) :: inputdir ! The directory in which to find input files. character(len=200) :: wave_drag_file ! The file from which to read the wave ! drag piston velocity. @@ -4721,6 +4734,17 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, "piston velocities.", default=1.0, units="nondim", & do_not_log=.not.CS%linear_wave_drag) + 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") + 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") + om2 = 1.4051890e-4 + ok1 = 0.7292117e-4 + call get_param(param_file, mdl, "CLIP_BT_VELOCITY", CS%clip_velocity, & "If true, limit any velocity components that exceed "//& "CFL_TRUNCATE. This should only be used as a desperate "//& @@ -4962,6 +4986,13 @@ 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 and register streaming filters + call Filt_init(CS%Filt_CS) + call Filt_register('um2', am2, om2, CS%Filt_CS) + call Filt_register('vm2', am2, om2, CS%Filt_CS) + call Filt_register('uk1', ak1, ok1, CS%Filt_CS) + call Filt_register('vk1', ak1, ok1, CS%Filt_CS) + CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input dtbt_tmp = -1.0 diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 new file mode 100644 index 0000000000..96eb65ecab --- /dev/null +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -0,0 +1,178 @@ +!> 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, NOTE +use MOM_time_manager, only : time_type, time_type_to_real +use MOM_unit_scaling, only : unit_scale_type + +implicit none ; private + +public Filt_init, Filt_register, Filt_accum + +#include + +!> The private control structure for storing the filter infomation of a particular field +type, private :: Filt_type + character(len=16) :: key = 'none' !< Name of the current field + real :: a, & !< Parameter that determines the bandwidth [nondim] + om, & !< Target frequency of the filter [s-1] + old_time = -1.0 !< The time of the previous accumulating step [T ~> s] + integer :: is, ie, js, je !< Lower and upper bounds of input data + real, allocatable :: s1(:,:), & !< Dummy variable [A] + u1(:,:) !< Filtered data [A] +end type Filt_type + +!> A linked list of control structures that store the filter infomation of different fields +type, private :: Filt_node + type(Filt_type) :: this !< Control structure of the current field in the list + type(Filt_node), pointer :: next !< The list of other fields +end type Filt_node + +!> The public control structure of the MOM_streaming_filter module +type, public :: streaming_filter_CS ; private + integer :: length !< Number of fields to be filtered + type(Filt_node), pointer :: list => NULL() !< A linked list for storing filter info of different fields +end type streaming_filter_CS + +contains + +!> This subroutine initializes CS%list. +subroutine Filt_init(CS) + type(streaming_filter_CS), intent(out) :: CS !< Control structure of the MOM_streaming_filter module + + ! Local variable + type(Filt_type) :: ha1 !< A temporary, null field used for initializing CS%list + + allocate(CS%list) + CS%list%this = ha1 + nullify(CS%list%next) + CS%length = 0 + +end subroutine Filt_init + +!> This subroutine registers each of the fields to be filtered. +subroutine Filt_register(key, a, om, CS) + character(len=*), intent(in) :: key !< Name of the current field + real, intent(in) :: a !< Parameter that determines the bandwidth [nondim] + real, intent(in) :: om !< Target frequency of the filter [s-1] + type(streaming_filter_CS), intent(inout) :: CS !< Control structure of the MOM_streaming_filter module + + ! Local variables + type(Filt_type) :: ha1 !< Control structure for the current field + type(Filt_node), pointer :: tmp !< A temporary list to hold the current field + character(len=128) :: mesg + + if (a <= 0.0) then + write(mesg, *) "MOM_streaming_filter: bandwidth <= 0, key ", trim(key), " not registered." + call MOM_error(NOTE, trim(mesg)) + return !< Setting a <= 0.0 turns off the filter + endif + if (om <= 0.0) then + write(mesg, *) "MOM_streaming_filter: target frequency <= 0, key ", trim(key), " not registered." + call MOM_error(NOTE, trim(mesg)) + return !< Setting om <= 0.0 turns off the filter + endif + + ha1%key = trim(key) + ha1%a = a + ha1%om = om + allocate(tmp) + tmp%this = ha1 + tmp%next => CS%list + CS%list => tmp + CS%length = CS%length + 1 + +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(key, u, u1, Time, US, CS) + character(len=*), intent(in) :: key !< Name of the current field + real, dimension(:,:), intent(in) :: u !< Input into the filter [A] + 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(streaming_filter_CS), intent(in) :: CS !< Control structure of the MOM_streaming_filter module + + ! Local variables + type(Filt_type), pointer :: ha1 + type(Filt_node), pointer :: tmp + 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] + real, allocatable, target :: u0(:,:) !< Output (zeros) when a field is not registered [A] + integer :: i, j, k, is, ie, js, je + character(len=128) :: mesg + + ! Exit the accumulator and return zeros if no field is registered + if (CS%length == 0) then + is = LBOUND(u,1) ; ie = UBOUND(u,1) ; js = LBOUND(u,2) ; je = UBOUND(u,2) + allocate(u0(is:ie,js:je), source=0.0) ; u0(:,:) = 0.0 ; u1 => u0 ; return + endif + + ! Loop through the full list to find the current field + tmp => CS%list + do k=1,CS%length + ha1 => tmp%this + if (trim(key) == trim(ha1%key)) exit + tmp => tmp%next + + ! Exit the accumulator and return zeros if the field is not registered + if (k == CS%length) then + is = LBOUND(u,1) ; ie = UBOUND(u,1) ; js = LBOUND(u,2) ; je = UBOUND(u,2) + allocate(u0(is:ie,js:je), source=0.0) ; u0(:,:) = 0.0 ; u1 => u0 ; return + endif + enddo + + now = US%s_to_T * time_type_to_real(Time) + + ! Additional processing at the initial accumulating step + if (ha1%old_time < 0.0) then + ha1%old_time = now + + write(mesg,*) "MOM_streaming_filter: initializing filter equations, key = ", trim(ha1%key) + call MOM_error(NOTE, trim(mesg)) + + ha1%is = LBOUND(u,1) ; is = ha1%is + ha1%ie = UBOUND(u,1) ; ie = ha1%ie + ha1%js = LBOUND(u,2) ; js = ha1%js + ha1%je = UBOUND(u,2) ; je = ha1%je + + allocate(ha1%s1(is:ie,js:je), source=0.0) + allocate(ha1%u1(is:ie,js:je), source=0.0) + do j=js,je ; do i=is,ie + ha1%s1(i,j) = 0.0 + ha1%u1(i,j) = u(i,j) + end do ; end do + endif + + dt = now - ha1%old_time + ha1%old_time = now + + is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je + + ! Timestepping + c1 = ha1%om * dt + c2 = 1.0 - ha1%a * c1 + do j=js,je ; do i=is,ie + ha1%s1(i,j) = c1 * (ha1%u1(i,j) + (ha1%a**2) * u(i,j)) + c2 * ha1%s1(i,j) + ha1%u1(i,j) = -c1 * (ha1%s1(i,j) + (-2*ha1%a) * u(i,j)) + c2 * ha1%u1(i,j) + enddo; enddo + + u1 => ha1%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 + From bbb629e25f98aaef3a54b44ba051f4d1778f11ec Mon Sep 17 00:00:00 2001 From: William Xu Date: Mon, 22 Jul 2024 17:24:48 -0500 Subject: [PATCH 2/9] Streaming Filter Fixed dimensional inconsistency of the om2 and ok1 parameters. --- src/core/MOM_barotropic.F90 | 16 +++++++++++++--- .../lateral/MOM_streaming_filter.F90 | 4 ++-- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 2c2511e353..750d1607d0 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -4466,7 +4466,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, real :: wave_drag_scale ! A scaling factor for the barotropic linear wave drag ! piston velocities [nondim]. 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 [s-1] + real :: om2, ok1 !< Target frequencies of the M2 and K1 streaming filters [T-1 ~> s-1] character(len=200) :: inputdir ! The directory in which to find input files. character(len=200) :: wave_drag_file ! The file from which to read the wave ! drag piston velocity. @@ -4742,8 +4742,18 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, "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") - om2 = 1.4051890e-4 - ok1 = 0.7292117e-4 + 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=1.4051890e-4, & + 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=0.7292117e-4, & + scale=US%T_to_s, do_not_log=.true.) call get_param(param_file, mdl, "CLIP_BT_VELOCITY", CS%clip_velocity, & "If true, limit any velocity components that exceed "//& diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index 96eb65ecab..b6f573b68b 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -15,7 +15,7 @@ module MOM_streaming_filter type, private :: Filt_type character(len=16) :: key = 'none' !< Name of the current field real :: a, & !< Parameter that determines the bandwidth [nondim] - om, & !< Target frequency of the filter [s-1] + om, & !< Target frequency of the filter [T-1 ~> s-1] old_time = -1.0 !< The time of the previous accumulating step [T ~> s] integer :: is, ie, js, je !< Lower and upper bounds of input data real, allocatable :: s1(:,:), & !< Dummy variable [A] @@ -54,7 +54,7 @@ end subroutine Filt_init subroutine Filt_register(key, a, om, CS) character(len=*), intent(in) :: key !< Name of the current field real, intent(in) :: a !< Parameter that determines the bandwidth [nondim] - real, intent(in) :: om !< Target frequency of the filter [s-1] + real, intent(in) :: om !< Target frequency of the filter [T-1 ~> s-1] type(streaming_filter_CS), intent(inout) :: CS !< Control structure of the MOM_streaming_filter module ! Local variables From 060eca83bfadbe23333fed30c07d0a2421b8cab4 Mon Sep 17 00:00:00 2001 From: William Xu Date: Sun, 28 Jul 2024 14:28:13 -0500 Subject: [PATCH 3/9] Streaming Filter 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. --- src/core/MOM_barotropic.F90 | 52 +++++-- .../lateral/MOM_streaming_filter.F90 | 144 +++++------------- 2 files changed, 75 insertions(+), 121 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 750d1607d0..9f9b3d52fc 100644 --- a/src/core/MOM_barotropic.F90 +++ b/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_init, Filt_register, Filt_accum, streaming_filter_CS +use MOM_streaming_filter, only : Filt_register, 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 @@ -249,6 +249,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 @@ -292,7 +296,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(streaming_filter_CS) :: Filt_CS !< Control structure for streaming filters + type(Filter_CS) :: Filt_CS_um2, Filt_CS_vm2, & !< Control structures for the M2 streaming filter + Filt_CS_uk1, 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. @@ -1592,10 +1597,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! 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. - call Filt_accum('um2', ubt, um2, CS%Time, US, CS%Filt_CS) - call Filt_accum('uk1', ubt, uk1, CS%Time, US, CS%Filt_CS) - call Filt_accum('vm2', vbt, vm2, CS%Time, US, CS%Filt_CS) - call Filt_accum('vk1', vbt, vk1, CS%Time, US, CS%Filt_CS) + 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 @@ -4734,14 +4743,20 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, "piston velocities.", default=1.0, units="nondim", & do_not_log=.not.CS%linear_wave_drag) + 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") + 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") + 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"// & @@ -4997,11 +5012,22 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, endif ! CS%linear_wave_drag ! Initialize and register streaming filters - call Filt_init(CS%Filt_CS) - call Filt_register('um2', am2, om2, CS%Filt_CS) - call Filt_register('vm2', am2, om2, CS%Filt_CS) - call Filt_register('uk1', ak1, ok1, CS%Filt_CS) - call Filt_register('vk1', ak1, ok1, CS%Filt_CS) + if (CS%use_filter_m2) then + if (am2>0 .and. om2>0) then + call Filt_register(am2, om2, CS%Filt_CS_um2) + call Filt_register(am2, om2, CS%Filt_CS_vm2) + else + CS%use_filter_m2 = .false. + endif + endif + if (CS%use_filter_k1) then + if (ak1>0 .and. ok1>0) then + call Filt_register(ak1, ok1, CS%Filt_CS_uk1) + call Filt_register(ak1, ok1, CS%Filt_CS_vk1) + else + CS%use_filter_k1 = .false. + endif + endif CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index b6f573b68b..aa67ad778f 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -1,165 +1,93 @@ !> 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, NOTE +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL use MOM_time_manager, only : time_type, time_type_to_real use MOM_unit_scaling, only : unit_scale_type implicit none ; private -public Filt_init, Filt_register, Filt_accum +public Filt_register, Filt_accum #include -!> The private control structure for storing the filter infomation of a particular field -type, private :: Filt_type - character(len=16) :: key = 'none' !< Name of the current field +!> 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] - integer :: is, ie, js, je !< Lower and upper bounds of input data real, allocatable :: s1(:,:), & !< Dummy variable [A] u1(:,:) !< Filtered data [A] -end type Filt_type - -!> A linked list of control structures that store the filter infomation of different fields -type, private :: Filt_node - type(Filt_type) :: this !< Control structure of the current field in the list - type(Filt_node), pointer :: next !< The list of other fields -end type Filt_node - -!> The public control structure of the MOM_streaming_filter module -type, public :: streaming_filter_CS ; private - integer :: length !< Number of fields to be filtered - type(Filt_node), pointer :: list => NULL() !< A linked list for storing filter info of different fields -end type streaming_filter_CS + !>@{ Lower and upper bounds of input data + integer :: is, ie, js, je + !>@} +end type Filter_CS contains -!> This subroutine initializes CS%list. -subroutine Filt_init(CS) - type(streaming_filter_CS), intent(out) :: CS !< Control structure of the MOM_streaming_filter module - - ! Local variable - type(Filt_type) :: ha1 !< A temporary, null field used for initializing CS%list - - allocate(CS%list) - CS%list%this = ha1 - nullify(CS%list%next) - CS%length = 0 - -end subroutine Filt_init - !> This subroutine registers each of the fields to be filtered. -subroutine Filt_register(key, a, om, CS) - character(len=*), intent(in) :: key !< Name of the current field +subroutine Filt_register(a, om, CS) real, intent(in) :: a !< Parameter that determines the bandwidth [nondim] real, intent(in) :: om !< Target frequency of the filter [T-1 ~> s-1] - type(streaming_filter_CS), intent(inout) :: CS !< Control structure of the MOM_streaming_filter module + type(Filter_CS), intent(out) :: CS !< Control structure for the current field - ! Local variables - type(Filt_type) :: ha1 !< Control structure for the current field - type(Filt_node), pointer :: tmp !< A temporary list to hold the current field - character(len=128) :: mesg - - if (a <= 0.0) then - write(mesg, *) "MOM_streaming_filter: bandwidth <= 0, key ", trim(key), " not registered." - call MOM_error(NOTE, trim(mesg)) - return !< Setting a <= 0.0 turns off the filter - endif - if (om <= 0.0) then - write(mesg, *) "MOM_streaming_filter: target frequency <= 0, key ", trim(key), " not registered." - call MOM_error(NOTE, trim(mesg)) - return !< Setting om <= 0.0 turns off the filter - endif + 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") - ha1%key = trim(key) - ha1%a = a - ha1%om = om - allocate(tmp) - tmp%this = ha1 - tmp%next => CS%list - CS%list => tmp - CS%length = CS%length + 1 + CS%a = a + CS%om = om 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(key, u, u1, Time, US, CS) - character(len=*), intent(in) :: key !< Name of the current field +subroutine Filt_accum(u, u1, Time, US, CS) real, dimension(:,:), intent(in) :: u !< Input into the filter [A] 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(streaming_filter_CS), intent(in) :: CS !< Control structure of the MOM_streaming_filter module + type(Filter_CS), target, intent(inout) :: CS !< Control structure of the MOM_streaming_filter module ! Local variables - type(Filt_type), pointer :: ha1 - type(Filt_node), pointer :: tmp 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] - real, allocatable, target :: u0(:,:) !< Output (zeros) when a field is not registered [A] - integer :: i, j, k, is, ie, js, je - character(len=128) :: mesg - - ! Exit the accumulator and return zeros if no field is registered - if (CS%length == 0) then - is = LBOUND(u,1) ; ie = UBOUND(u,1) ; js = LBOUND(u,2) ; je = UBOUND(u,2) - allocate(u0(is:ie,js:je), source=0.0) ; u0(:,:) = 0.0 ; u1 => u0 ; return - endif - - ! Loop through the full list to find the current field - tmp => CS%list - do k=1,CS%length - ha1 => tmp%this - if (trim(key) == trim(ha1%key)) exit - tmp => tmp%next - - ! Exit the accumulator and return zeros if the field is not registered - if (k == CS%length) then - is = LBOUND(u,1) ; ie = UBOUND(u,1) ; js = LBOUND(u,2) ; je = UBOUND(u,2) - allocate(u0(is:ie,js:je), source=0.0) ; u0(:,:) = 0.0 ; u1 => u0 ; return - endif - enddo + integer :: i, j, is, ie, js, je now = US%s_to_T * time_type_to_real(Time) ! Additional processing at the initial accumulating step - if (ha1%old_time < 0.0) then - ha1%old_time = now + if (CS%old_time < 0.0) then + CS%old_time = now - write(mesg,*) "MOM_streaming_filter: initializing filter equations, key = ", trim(ha1%key) - call MOM_error(NOTE, trim(mesg)) + CS%is = LBOUND(u,1) ; is = CS%is + CS%ie = UBOUND(u,1) ; ie = CS%ie + CS%js = LBOUND(u,2) ; js = CS%js + CS%je = UBOUND(u,2) ; je = CS%je - ha1%is = LBOUND(u,1) ; is = ha1%is - ha1%ie = UBOUND(u,1) ; ie = ha1%ie - ha1%js = LBOUND(u,2) ; js = ha1%js - ha1%je = UBOUND(u,2) ; je = ha1%je + allocate(CS%s1(is:ie,js:je), source=0.0) + allocate(CS%u1(is:ie,js:je), source=0.0) - allocate(ha1%s1(is:ie,js:je), source=0.0) - allocate(ha1%u1(is:ie,js:je), source=0.0) do j=js,je ; do i=is,ie - ha1%s1(i,j) = 0.0 - ha1%u1(i,j) = u(i,j) - end do ; end do + CS%s1(i,j) = 0.0 + CS%u1(i,j) = u(i,j) + enddo ; enddo endif - dt = now - ha1%old_time - ha1%old_time = now + dt = now - CS%old_time + CS%old_time = now - is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je + is = CS%is ; ie = CS%ie ; js = CS%js ; je = CS%je ! Timestepping - c1 = ha1%om * dt - c2 = 1.0 - ha1%a * c1 + c1 = CS%om * dt + c2 = 1.0 - CS%a * c1 do j=js,je ; do i=is,ie - ha1%s1(i,j) = c1 * (ha1%u1(i,j) + (ha1%a**2) * u(i,j)) + c2 * ha1%s1(i,j) - ha1%u1(i,j) = -c1 * (ha1%s1(i,j) + (-2*ha1%a) * u(i,j)) + c2 * ha1%u1(i,j) + CS%s1(i,j) = c1 * (CS%u1(i,j) + (CS%a**2) * u(i,j)) + c2 * CS%s1(i,j) + CS%u1(i,j) = -c1 * (CS%s1(i,j) + (-2*CS%a) * u(i,j)) + c2 * CS%u1(i,j) enddo; enddo - u1 => ha1%u1 + u1 => CS%u1 end subroutine Filt_accum From 7df399f03c4ee1b324ed957008d5c898f8fb9eab Mon Sep 17 00:00:00 2001 From: William Xu Date: Sun, 28 Jul 2024 16:57:46 -0400 Subject: [PATCH 4/9] Streaming filter Minor fix on documentation --- src/core/MOM_barotropic.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 9f9b3d52fc..06eb2ec143 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -296,8 +296,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, Filt_CS_vm2, & !< Control structures for the M2 streaming filter - Filt_CS_uk1, Filt_CS_vk1 !< Control structures for the K1 streaming filter + 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. From cbebaaedcb59846491990d53e49380128dbc8b80 Mon Sep 17 00:00:00 2001 From: William Xu Date: Sun, 18 Aug 2024 21:00:52 -0300 Subject: [PATCH 5/9] Streaming filter Filter equations reformulated. --- src/parameterizations/lateral/MOM_streaming_filter.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index aa67ad778f..518461f5f6 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -83,8 +83,8 @@ subroutine Filt_accum(u, u1, Time, US, CS) 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%a**2) * u(i,j)) + c2 * CS%s1(i,j) - CS%u1(i,j) = -c1 * (CS%s1(i,j) + (-2*CS%a) * u(i,j)) + c2 * CS%u1(i,j) + 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 From 9aa335d0d3125f58d3392063789d3b6178f25def Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 23 Aug 2024 12:59:30 -0300 Subject: [PATCH 6/9] Streaming filter Using hor_index_type argument to avoid calculation in halo regions. --- src/core/MOM_barotropic.F90 | 102 ++++++++-------- .../lateral/MOM_streaming_filter.F90 | 113 ++++++++++++------ 2 files changed, 128 insertions(+), 87 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 06eb2ec143..0e6f1fa2cb 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1600,12 +1600,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! 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) + call Filt_accum(ubt, um2, 'u', CS%Time, US, CS%Filt_CS_um2) + call Filt_accum(vbt, vm2, 'v', 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) + call Filt_accum(ubt, uk1, 'u', CS%Time, US, CS%Filt_CS_uk1) + call Filt_accum(vbt, vk1, 'v', CS%Time, US, CS%Filt_CS_vk1) endif ! Zero out the arrays for various time-averaged quantities. @@ -4476,8 +4476,6 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, real :: dtbt_tmp ! A temporary copy of CS%dtbt read from a restart file [T ~> s] real :: wave_drag_scale ! A scaling factor for the barotropic linear wave drag ! piston velocities [nondim]. - 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] character(len=200) :: inputdir ! The directory in which to find input files. character(len=200) :: wave_drag_file ! The file from which to read the wave ! drag piston velocity. @@ -4745,33 +4743,6 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, "piston velocities.", default=1.0, units="nondim", & do_not_log=.not.CS%linear_wave_drag) - 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=1.4051890e-4, & - 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=0.7292117e-4, & - scale=US%T_to_s, do_not_log=.true.) - call get_param(param_file, mdl, "CLIP_BT_VELOCITY", CS%clip_velocity, & "If true, limit any velocity components that exceed "//& "CFL_TRUNCATE. This should only be used as a desperate "//& @@ -5013,24 +4984,6 @@ 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 and register streaming filters - if (CS%use_filter_m2) then - if (am2>0 .and. om2>0) then - call Filt_register(am2, om2, CS%Filt_CS_um2) - call Filt_register(am2, om2, CS%Filt_CS_vm2) - else - CS%use_filter_m2 = .false. - endif - endif - if (CS%use_filter_k1) then - if (ak1>0 .and. ok1>0) then - call Filt_register(ak1, ok1, CS%Filt_CS_uk1) - call Filt_register(ak1, ok1, CS%Filt_CS_vk1) - else - CS%use_filter_k1 = .false. - endif - endif - CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input dtbt_tmp = -1.0 @@ -5316,6 +5269,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 @@ -5328,6 +5283,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=1.4051890e-4, & + 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=0.7292117e-4, & + 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 @@ -5356,6 +5338,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 diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index 518461f5f6..4b7e911d3c 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -2,6 +2,7 @@ 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 @@ -13,11 +14,15 @@ module MOM_streaming_filter !> 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 :: s1(:,:), & !< Dummy variable [A] - u1(:,:) !< Filtered data [A] + 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 ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: s1_h, & !< Dummy variable on h grid [A] + u1_h !< Filtered data on h grid [A] + real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: s1_u, & !< Dummy variable on u grid [A] + u1_u !< Filtered data on u grid [A] + real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: s1_v, & !< Dummy variable on v grid [A] + u1_v !< Filtered data on v grid [A] !>@{ Lower and upper bounds of input data integer :: is, ie, js, je !>@} @@ -26,68 +31,104 @@ module MOM_streaming_filter contains !> This subroutine registers each of the fields to be filtered. -subroutine Filt_register(a, om, CS) +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') + ALLOC_(CS%s1_h(isd:ied,jsd:jed)) ; CS%s1_h(:,:) = 0.0 + ALLOC_(CS%u1_h(isd:ied,jsd:jed)) ; CS%u1_h(:,:) = 0.0 + CS%is = isd ; CS%ie = ied ; CS%js = jsd ; CS%je = jed + case ('u') + ALLOC_(CS%s1_u(IsdB:IedB,jsd:jed)) ; CS%s1_u(:,:) = 0.0 + ALLOC_(CS%u1_u(IsdB:IedB,jsd:jed)) ; CS%u1_u(:,:) = 0.0 + CS%is = IsdB ; CS%ie = IedB ; CS%js = jsd ; CS%je = jed + case ('v') + ALLOC_(CS%s1_v(isd:ied,JsdB:JedB)) ; CS%s1_v(:,:) = 0.0 + ALLOC_(CS%u1_v(isd:ied,JsdB:JedB)) ; CS%u1_v(:,:) = 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(:,:), intent(in) :: u !< Input into the filter [A] - 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 +subroutine Filt_accum(u, u1, grid, Time, US, CS) + real, dimension(:,:), pointer, intent(out) :: u1 !< Output of the filter [A] + character(len=*), intent(in) :: grid !< Horizontal grid location: h, u, or v + 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 + 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 - ! Additional processing at the initial accumulating step + ! Initialize u1 if (CS%old_time < 0.0) then CS%old_time = now - CS%is = LBOUND(u,1) ; is = CS%is - CS%ie = UBOUND(u,1) ; ie = CS%ie - CS%js = LBOUND(u,2) ; js = CS%js - CS%je = UBOUND(u,2) ; je = CS%je - - allocate(CS%s1(is:ie,js:je), source=0.0) - allocate(CS%u1(is:ie,js:je), source=0.0) - - do j=js,je ; do i=is,ie - CS%s1(i,j) = 0.0 - CS%u1(i,j) = u(i,j) - enddo ; enddo + select case (trim(grid)) + case ('h') ; CS%u1_h(:,:) = u(:,:) + case ('u') ; CS%u1_u(:,:) = u(:,:) + case ('v') ; CS%u1_v(:,:) = u(:,:) + case default ; call MOM_error(FATAL, "MOM_streaming_filter: horizontal grid not supported") + end select endif dt = now - CS%old_time CS%old_time = now - is = CS%is ; ie = CS%ie ; js = CS%js ; je = CS%je - ! 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 + select case (trim(grid)) + case ('h') + do j=js,je ; do i=is,ie + CS%s1_h(i,j) = c1 * CS%u1_h(i,j) + CS%s1_h(i,j) + CS%u1_h(i,j) = -c1 * (CS%s1_h(i,j) - CS%a * u(i,j)) + c2 * CS%u1_h(i,j) + enddo; enddo + u1 => CS%u1_h + case ('u') + do j=js,je ; do i=is,ie + CS%s1_u(i,j) = c1 * CS%u1_u(i,j) + CS%s1_u(i,j) + CS%u1_u(i,j) = -c1 * (CS%s1_u(i,j) - CS%a * u(i,j)) + c2 * CS%u1_u(i,j) + enddo; enddo + u1 => CS%u1_u + case ('v') + do j=js,je ; do i=is,ie + CS%s1_v(i,j) = c1 * CS%u1_v(i,j) + CS%s1_v(i,j) + CS%u1_v(i,j) = -c1 * (CS%s1_v(i,j) - CS%a * u(i,j)) + c2 * CS%u1_v(i,j) + enddo; enddo + u1 => CS%u1_h + case default + call MOM_error(FATAL, "MOM_streaming_filter: horizontal grid not supported") + end select end subroutine Filt_accum From f84b9f444bde14a7886f8e575ee50b12647ba8df Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 23 Aug 2024 15:56:38 -0300 Subject: [PATCH 7/9] Streaming filter Bug fix. --- src/parameterizations/lateral/MOM_streaming_filter.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index 4b7e911d3c..0af2994107 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -125,7 +125,7 @@ subroutine Filt_accum(u, u1, grid, Time, US, CS) CS%s1_v(i,j) = c1 * CS%u1_v(i,j) + CS%s1_v(i,j) CS%u1_v(i,j) = -c1 * (CS%s1_v(i,j) - CS%a * u(i,j)) + c2 * CS%u1_v(i,j) enddo; enddo - u1 => CS%u1_h + u1 => CS%u1_v case default call MOM_error(FATAL, "MOM_streaming_filter: horizontal grid not supported") end select From 46e479786e369e2779003768adb7d36beda25b11 Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 30 Aug 2024 20:09:20 -0300 Subject: [PATCH 8/9] Streaming filter Minor update on the M2 and K1 frequencies. --- src/core/MOM_barotropic.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 0e6f1fa2cb..56c6e0fbf9 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -27,6 +27,7 @@ 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 @@ -5301,13 +5302,13 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) "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=1.4051890e-4, & + " 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=0.7292117e-4, & + " 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 From cda34b191e60786bb1ed3eb0cbfbc3dd60a78949 Mon Sep 17 00:00:00 2001 From: William Xu Date: Sat, 7 Sep 2024 09:14:50 -0300 Subject: [PATCH 9/9] Streaming filter Using non-static-memory allocation for s1 and u1. --- src/core/MOM_barotropic.F90 | 8 +-- .../lateral/MOM_streaming_filter.F90 | 58 +++++-------------- 2 files changed, 19 insertions(+), 47 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 56c6e0fbf9..0ac3e52a62 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1601,12 +1601,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! 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, 'u', CS%Time, US, CS%Filt_CS_um2) - call Filt_accum(vbt, vm2, 'v', CS%Time, US, CS%Filt_CS_vm2) + 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, 'u', CS%Time, US, CS%Filt_CS_uk1) - call Filt_accum(vbt, vk1, 'v', CS%Time, US, CS%Filt_CS_vk1) + 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. diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index 0af2994107..a91f6661f2 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -17,12 +17,8 @@ module MOM_streaming_filter 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 ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: s1_h, & !< Dummy variable on h grid [A] - u1_h !< Filtered data on h grid [A] - real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: s1_u, & !< Dummy variable on u grid [A] - u1_u !< Filtered data on u grid [A] - real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: s1_v, & !< Dummy variable on v grid [A] - u1_v !< Filtered data on v grid [A] + real, allocatable, dimension(:,:) :: s1, & !< Dummy variable [A] + u1 !< Filtered data [A] !>@{ Lower and upper bounds of input data integer :: is, ie, js, je !>@} @@ -52,16 +48,16 @@ subroutine Filt_register(a, om, grid, HI, CS) select case (trim(grid)) case ('h') - ALLOC_(CS%s1_h(isd:ied,jsd:jed)) ; CS%s1_h(:,:) = 0.0 - ALLOC_(CS%u1_h(isd:ied,jsd:jed)) ; CS%u1_h(:,:) = 0.0 + 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') - ALLOC_(CS%s1_u(IsdB:IedB,jsd:jed)) ; CS%s1_u(:,:) = 0.0 - ALLOC_(CS%u1_u(IsdB:IedB,jsd:jed)) ; CS%u1_u(:,:) = 0.0 + 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') - ALLOC_(CS%s1_v(isd:ied,JsdB:JedB)) ; CS%s1_v(:,:) = 0.0 - ALLOC_(CS%u1_v(isd:ied,JsdB:JedB)) ; CS%u1_v(:,:) = 0.0 + 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") @@ -71,9 +67,8 @@ 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, grid, Time, US, CS) +subroutine Filt_accum(u, u1, Time, US, CS) real, dimension(:,:), pointer, intent(out) :: u1 !< Output of the filter [A] - character(len=*), intent(in) :: grid !< Horizontal grid location: h, u, or v 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 @@ -91,13 +86,7 @@ subroutine Filt_accum(u, u1, grid, Time, US, CS) ! Initialize u1 if (CS%old_time < 0.0) then CS%old_time = now - - select case (trim(grid)) - case ('h') ; CS%u1_h(:,:) = u(:,:) - case ('u') ; CS%u1_u(:,:) = u(:,:) - case ('v') ; CS%u1_v(:,:) = u(:,:) - case default ; call MOM_error(FATAL, "MOM_streaming_filter: horizontal grid not supported") - end select + CS%u1(:,:) = u(:,:) endif dt = now - CS%old_time @@ -107,28 +96,11 @@ subroutine Filt_accum(u, u1, grid, Time, US, CS) c1 = CS%om * dt c2 = 1.0 - CS%a * c1 - select case (trim(grid)) - case ('h') - do j=js,je ; do i=is,ie - CS%s1_h(i,j) = c1 * CS%u1_h(i,j) + CS%s1_h(i,j) - CS%u1_h(i,j) = -c1 * (CS%s1_h(i,j) - CS%a * u(i,j)) + c2 * CS%u1_h(i,j) - enddo; enddo - u1 => CS%u1_h - case ('u') - do j=js,je ; do i=is,ie - CS%s1_u(i,j) = c1 * CS%u1_u(i,j) + CS%s1_u(i,j) - CS%u1_u(i,j) = -c1 * (CS%s1_u(i,j) - CS%a * u(i,j)) + c2 * CS%u1_u(i,j) - enddo; enddo - u1 => CS%u1_u - case ('v') - do j=js,je ; do i=is,ie - CS%s1_v(i,j) = c1 * CS%u1_v(i,j) + CS%s1_v(i,j) - CS%u1_v(i,j) = -c1 * (CS%s1_v(i,j) - CS%a * u(i,j)) + c2 * CS%u1_v(i,j) - enddo; enddo - u1 => CS%u1_v - case default - call MOM_error(FATAL, "MOM_streaming_filter: horizontal grid not supported") - end select + 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