diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 9641eb294f..b424f8a464 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_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 diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index 3e8a932ef2..eac8d13c91 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/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 !> 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.