Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

updating #5

Merged
merged 56 commits into from
Sep 23, 2024
Merged
Changes from 1 commit
Commits
Show all changes
56 commits
Select commit Hold shift + click to select a range
6628dbf
(*)Parenthesize squares of wind stresses for FMAs
Hallberg-NOAA Mar 1, 2024
f0e61f3
(*)Parenthesize continuity_PPM curv_3 expressions
Hallberg-NOAA Feb 29, 2024
4f710ef
(*)Add parentheses for oblique OBCs with FMAs
Hallberg-NOAA Mar 1, 2024
9172cd5
(*)Add parentheses for density_integrals with FMAs
Hallberg-NOAA Mar 1, 2024
24091cc
(*)Add parentheses to 4 EOS int_density routines
Hallberg-NOAA Mar 1, 2024
8066a3d
(*)Simplify density integral parentheses
Hallberg-NOAA Mar 4, 2024
99fd957
(*)Parenthesize PressureForce_Montgomery for FMAs
Hallberg-NOAA Mar 1, 2024
307a4e2
(*)Parenthesize calc_isoneutral_slopes for FMAs
Hallberg-NOAA Mar 1, 2024
ce559ce
(*)Parenthesize MOM_calc_varT for FMAs
Hallberg-NOAA Mar 1, 2024
c344a11
(*)Parenthesize tracer_hordiff for FMAs
Hallberg-NOAA Mar 1, 2024
b2beab2
(*)Parenthesize iceberg_forces for FMAs
Hallberg-NOAA Mar 1, 2024
5398e6f
(*)Parenthesize CoriolisStokes and LA_Stk for FMAs
Hallberg-NOAA Mar 1, 2024
56d053a
+(*)Add and use G%Coriolis2Bu
Hallberg-NOAA Mar 1, 2024
49419f7
(*)Parenthesize thickness_diffuse for FMAs
Hallberg-NOAA Mar 1, 2024
03dc6f9
(*)Parenthesize Zanna_Bolton for FMAs
Hallberg-NOAA Mar 1, 2024
654cd4a
(*)Parenthesize MOM_internal_tides for FMAs
Hallberg-NOAA Mar 1, 2024
ebf02a9
(*)Parenthesize find_uv_at_h for FMAs
Hallberg-NOAA Mar 1, 2024
c0bef18
(*)Parenthesize set_viscous_ML for FMAs
Hallberg-NOAA Mar 1, 2024
0b50a15
(*)Rearrange calc_kappa_shear_vertex for FMAs
Hallberg-NOAA Mar 1, 2024
f0c52dd
(*)Parenthesize MOM_set_diffusivity for FMAs
Hallberg-NOAA Mar 1, 2024
64b851c
(*)Parenthesize CorAdCalc for FMAs
Hallberg-NOAA Mar 1, 2024
46e8b66
(*)Parenthesize MOM_barotropic for FMAs
Hallberg-NOAA Mar 1, 2024
6216fa1
(*)Parenthesize MOM_lateral_mixing_coeffs for FMAs
Hallberg-NOAA Mar 2, 2024
ffef92f
(*)Parenthesize MOM_hor_visc for FMAs
Hallberg-NOAA Apr 18, 2024
fc2af28
(*)Parenthesize initialization squares for FMAs
Hallberg-NOAA Mar 3, 2024
44f1130
(*)Parenthesize parameterization squares for FMAs
Hallberg-NOAA Mar 3, 2024
182223c
(*)Parenthesize diagnostics for FMAs
Hallberg-NOAA Apr 30, 2024
e810ac5
(*)Parenthesize tracer_advect PPM edge values
Hallberg-NOAA May 5, 2024
ffa766b
(*)More parentheses in density_integrals for FMAs
Hallberg-NOAA Jul 31, 2024
fd82861
(*)Add parentheses in end_value_h4 for FMAs
Hallberg-NOAA Aug 2, 2024
ce58a32
Merge pull request #1634 from Hallberg-NOAA/FMA_rotational_symmetry_main
Hallberg-NOAA Aug 24, 2024
a6dd0fd
Makedep output cleanup and PEP8 fixes
marshallward Aug 6, 2024
1eccd28
Merge branch 'main' into dev/gfdl
Hallberg-NOAA Sep 4, 2024
91eee52
Inline harmonic analysis
c2xu Sep 6, 2024
2316ae5
diffusivities from internal tides ray tracing algo (#677)
raphaeldussin Sep 9, 2024
95744a7
Streaming filter (#675)
c2xu Sep 10, 2024
ffff6f3
+Optionally use SSH in calculate density for PGF
Hallberg-NOAA Jul 29, 2024
9b9c165
(*)Refactor p_ave calculation
Hallberg-NOAA Sep 7, 2024
5fc90eb
Rotate ice shelf forcing and initialization
Hallberg-NOAA Aug 17, 2024
70a48e3
+Add MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP
claireyung Apr 30, 2024
8520c9f
+Add top mass_weight_in_PGF option to 13 integrals
Hallberg-NOAA Jul 24, 2024
e172fe8
+Add CORRECTION_INTXPA
claireyung Apr 30, 2024
15ea628
+Add CORRECTION_INTXPA_5PT
claireyung Apr 30, 2024
1b9bf67
+Add RESET_INTXPA_INTEGRAL
claireyung Apr 30, 2024
7a9545a
Revisions of sub-ice pressure gradient fixes
Hallberg-NOAA Jul 19, 2024
4cf1590
*Refactor CORRECTION_INTX_PA
Hallberg-NOAA Jul 29, 2024
15fd31c
*Non-Boussinesq code for RESET_INTXPA_INTEGRAL
Hallberg-NOAA Aug 5, 2024
5fceecf
+(*)Add 5-point quadrature in RESET_INTXPA_INTEGRAL
Hallberg-NOAA Aug 20, 2024
bdf4b9e
+(*)Eliminate CORRECTION_INTXPA_5PT
Hallberg-NOAA Sep 16, 2024
0363d2b
*Set MASS_WEIGHT_IN_PRESSURE_GRADIENT in .testing
Hallberg-NOAA Sep 16, 2024
1830b8e
Dummy code to suppress errors in posix.F90
marshallward Sep 6, 2024
e05cc01
F2023: Fix argument orders and IO statements
marshallward Sep 9, 2024
b67e93a
Reorder arguments in FMS_cap functions
marshallward Sep 12, 2024
b2db6bf
CI: Fortran 2018 testing
marshallward Sep 12, 2024
b3d7348
Change the default of VISC_REM_CONT_HVEL_FIX
herrwang0 Sep 12, 2024
ba59078
Separate scalar diagnostics for each ice sheet + parameters to contro…
alex-huth Sep 18, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Streaming filter (NOAA-GFDL#675)
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.

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.

Using hor_index_type argument to avoid calculation in halo regions.
c2xu authored Sep 10, 2024

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
commit 95744a71bf558c9802d02ee8562d6c1d1c147c83
70 changes: 70 additions & 0 deletions src/core/MOM_barotropic.F90
Original file line number Diff line number Diff line change
@@ -26,6 +26,8 @@ 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_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
@@ -248,6 +250,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
@@ -291,6 +297,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, & !< 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.
@@ -598,6 +608,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 +1598,17 @@ 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)
endif

! Zero out the arrays for various time-averaged quantities.
if (find_etaav) then
!$OMP do
@@ -5247,6 +5270,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
@@ -5259,6 +5284,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=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.)

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
@@ -5287,6 +5339,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
119 changes: 119 additions & 0 deletions src/parameterizations/lateral/MOM_streaming_filter.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,119 @@
!> 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, 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

implicit none ; private

public Filt_register, Filt_accum

#include <MOM_memory.h>

!> 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, dimension(:,:) :: s1, & !< Dummy variable [A]
u1 !< Filtered data [A]
!>@{ Lower and upper bounds of input data
integer :: is, ie, js, je
!>@}
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

! 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')
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')
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')
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")
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(:,:), 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

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(:,:)
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
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.
!!
!! 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