From 4068335ef2a8e1208bab5c28b89dcf00ec49ac08 Mon Sep 17 00:00:00 2001 From: William Xu Date: Wed, 11 Dec 2024 22:00:02 -0400 Subject: [PATCH 1/2] Frequency-dependent internal wave drag Utilizing streaming band-pass filters, the frequency-dependent wave drag can be applied to the narrowband tidal velocities, instead of the broadband barotropic velocity. This allows internal wave drag to be optimized for different frequency bands independently, without affecting the low-frequency, non-tidal motions. Also fixed the bug that Filt_accum() does not return a valid output if it is called at the corrector step. --- src/core/MOM_barotropic.F90 | 85 ++++++++--- .../lateral/MOM_streaming_filter.F90 | 30 ++-- .../lateral/MOM_wave_drag.F90 | 135 ++++++++++++++++++ 3 files changed, 216 insertions(+), 34 deletions(-) create mode 100644 src/parameterizations/lateral/MOM_wave_drag.F90 diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 8ba4b106fd..442d9960a6 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -32,6 +32,7 @@ module MOM_barotropic use MOM_variables, only : BT_cont_type, alloc_bt_cont_type use MOM_verticalGrid, only : verticalGrid_type use MOM_variables, only : accel_diag_ptrs +use MOM_wave_drag, only : wave_drag_init, wave_drag_calc, wave_drag_CS implicit none ; private @@ -251,6 +252,8 @@ module MOM_barotropic !! invariant and linearized. logical :: use_filter !< If true, use streaming band-pass filter to detect the !! instantaneous tidal signals in the simulation. + logical :: linear_freq_drag !< If true, apply a linear frequency-dependent drag to the tidal + !! velocities. The streaming band-pass filter must be turned on. 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 @@ -296,6 +299,7 @@ module MOM_barotropic type(harmonic_analysis_CS), pointer :: HA_CSp => NULL() !< Control structure for harmonic analysis type(Filter_CS) :: Filt_CS_u, & !< Control structures for the streaming band-pass filter of ubt Filt_CS_v !< Control structures for the streaming band-pass filter of vbt + type(wave_drag_CS) :: Drag_CS !< Control structures for the frequency-dependent drag logical :: module_is_initialized = .false. !< If true, module has been initialized integer :: isdw !< The lower i-memory limit for the wide halo arrays. @@ -605,6 +609,10 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! spacing [H L ~> m2 or kg m-1]. real, dimension(:,:,:), pointer :: ufilt, vfilt ! Filtered velocities from the output of streaming filters [L T-1 ~> m s-1] + real, dimension(SZIBW_(CS),SZJW_(CS)) :: Drag_u + ! The zonal acceleration due to frequency-dependent drag [L T-2 ~> m s-2] + real, dimension(SZIW_(CS),SZJBW_(CS)) :: Drag_v + ! The meridional acceleration due to frequency-dependent drag [L T-2 ~> m s-2] real, target, dimension(SZIW_(CS),SZJW_(CS)) :: & eta, & ! The barotropic free surface height anomaly or column mass ! anomaly [H ~> m or kg m-2] @@ -1422,6 +1430,41 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo endif + ! Compute instantaneous tidal velocities and apply frequency-dependent drag. + ! Note that the filtered velocities are only updated during the current predictor step, + ! and are calculated using the barotropic velocity from the previous correction step. + if (CS%use_filter .and. CS%linear_freq_drag) 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) + call wave_drag_calc(ufilt, vfilt, Drag_u, Drag_v, G, CS%Drag_CS) + !$OMP do + do j=js,je ; do I=is-1,ie + Htot = 0.5 * (eta(i,j) + eta(i+1,j)) + if (GV%Boussinesq) & + Htot = Htot + 0.5*GV%Z_to_H * (CS%bathyT(i,j) + CS%bathyT(i+1,j)) + if (Htot > 0.0) then + Drag_u(I,j) = Drag_u(I,j) / Htot + BT_force_u(I,j) = BT_force_u(I,j) - Drag_u(I,j) + else + Drag_u(I,j) = 0.0 + endif + enddo ; enddo + !$OMP do + do J=js-1,je ; do i=is,ie + Htot = 0.5 * (eta(i,j) + eta(i,j+1)) + if (GV%Boussinesq) & + Htot = Htot + 0.5*GV%Z_to_H * (CS%bathyT(i,j) + CS%bathyT(i,j+1)) + if (Htot > 0.0) then + Drag_v(i,J) = Drag_v(i,J) / Htot + BT_force_v(i,J) = BT_force_v(i,J) - Drag_v(i,J) + else + Drag_v(i,J) = 0.0 + endif + enddo ; enddo + else + Drag_u(:,:) = 0.0 ; Drag_v(:,:) = 0.0 + endif + if ((Isq > is-1) .or. (Jsq > js-1)) then ! Non-symmetric memory is being used, so the edge values need to be ! filled in with a halo update of a non-symmetric array. @@ -1593,13 +1636,6 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif ; enddo ; enddo endif - ! Note that the filtered velocities are only updated during the current predictor step, - ! and are calculated using the barotropic velocity from the previous correction step. - 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. if (find_etaav) then !$OMP do @@ -2054,12 +2090,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & - ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) + ((Cor_v(i,J) + PFv(i,J)) - (vbt(i,J)*Rayleigh_v(i,J) + Drag_v(i,J))) enddo ; enddo else !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 - v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J) - Drag_v(i,J)) enddo ; enddo endif @@ -2132,13 +2168,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & - ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) + ((Cor_u(I,j) + PFu(I,j)) - (ubt(I,j)*Rayleigh_u(I,j) + Drag_u(I,j))) enddo ; enddo !$OMP end do nowait else !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev - u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j) - Drag_u(I,j)) enddo ; enddo !$OMP end do nowait endif @@ -2209,12 +2245,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & - ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) + ((Cor_u(I,j) + PFu(I,j)) - (ubt(I,j)*Rayleigh_u(I,j) + Drag_u(I,j))) enddo ; enddo else !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev - u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j) - Drag_u(I,j)) enddo ; enddo endif @@ -2298,13 +2334,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & - ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) + ((Cor_v(i,J) + PFv(i,J)) - (vbt(i,J)*Rayleigh_v(i,J) + Drag_v(i,J))) enddo ; enddo !$OMP end do nowait else !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev - v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J) - Drag_v(i,J)) enddo ; enddo !$OMP end do nowait endif @@ -4710,9 +4746,13 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & "using rates set by lin_drag_u & _v divided by the depth of "//& "the ocean. This was introduced to facilitate tide modeling.", & default=.false.) + call get_param(param_file, mdl, "BT_LINEAR_FREQ_DRAG", CS%linear_freq_drag, & + "If true, apply frequency-dependent drag to the tidal velocities. "//& + "The streaming band-pass filter must be turned on.", default=.false.) call get_param(param_file, mdl, "BT_WAVE_DRAG_FILE", wave_drag_file, & "The name of the file with the barotropic linear wave drag "//& - "piston velocities.", default="", do_not_log=.not.CS%linear_wave_drag) + "piston velocities.", default="", & + do_not_log=.not.CS%linear_wave_drag.and..not.CS%linear_freq_drag) call get_param(param_file, mdl, "BT_WAVE_DRAG_VAR", wave_drag_var, & "The name of the variable in BT_WAVE_DRAG_FILE with the "//& "barotropic linear wave drag piston velocities at h points. "//& @@ -4967,10 +5007,17 @@ subroutine barotropic_init(u, v, h, 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 + ! Initialize streaming band-pass filters and frequency-dependent drag + if (CS%use_filter .and. CS%linear_freq_drag) 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) + + if (.not.CS%linear_wave_drag .and. len_trim(wave_drag_file) > 0) then + inputdir = "." ; call get_param(param_file, mdl, "INPUTDIR", inputdir) + wave_drag_file = trim(slasher(inputdir))//trim(wave_drag_file) + call log_param(param_file, mdl, "INPUTDIR/BT_WAVE_DRAG_FILE", wave_drag_file) + endif + call wave_drag_init(param_file, wave_drag_file, G, GV, US, CS%Drag_CS) endif CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input @@ -5302,7 +5349,7 @@ 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 + ! 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.) diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index a2adaa2e92..7a8bc1b774 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -162,22 +162,22 @@ subroutine Filt_accum(u, u1, Time, US, CS) ! Initialize CS%old_time at the first time step if (CS%old_time<0.0) CS%old_time = now - ! This condition implies Filt_accum has already been called in the predictor step - if (CS%old_time>=now) return + ! Timestep the filter equations only if we are in a new time step + if (CS%old_time CS%u1 end subroutine Filt_accum diff --git a/src/parameterizations/lateral/MOM_wave_drag.F90 b/src/parameterizations/lateral/MOM_wave_drag.F90 new file mode 100644 index 0000000000..7ad340995e --- /dev/null +++ b/src/parameterizations/lateral/MOM_wave_drag.F90 @@ -0,0 +1,135 @@ +!> Frequency-dependent linear wave drag + +module MOM_wave_drag + +use MOM_domains, only : pass_vector, To_All, Scalar_Pair +use MOM_error_handler, only : MOM_error, NOTE +use MOM_file_parser, only : get_param, log_param, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_io, only : MOM_read_data, slasher, EAST_FACE, NORTH_FACE +use MOM_unit_scaling, only : unit_scale_type +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +public wave_drag_init, wave_drag_calc + +#include + +!> Control structure for the MOM_wave_drag module +type, public :: wave_drag_CS ; private + integer :: nf !< Number of filters to be used in the simulation + !>@{ Spatially varying, frequency-dependent drag coefficients [H T-1 ~> m s-1] + real, allocatable, dimension(:,:,:) :: coef_u, coef_v + !>@} +end type wave_drag_CS + +contains + +!> This subroutine reads drag coefficients from file. +subroutine wave_drag_init(param_file, wave_drag_file, G, GV, US, CS) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + character(len=*), intent(in) :: wave_drag_file !< The file from which to read drag coefficients + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(wave_drag_CS), intent(out) :: CS !< Control structure of MOM_wave_drag + + ! Local variables + character(len=40) :: mdl = "MOM_wave_drag" !< This module's name + character(len=50) :: filter_name_str !< List of drag coefficients to be used + character(len=2), allocatable, dimension(:) :: filter_names !< Names of drag coefficients + character(len=80) :: var_names(2) !< Names of variables in wave_drag_file + character(len=200) :: mesg + real :: var_scale !< Scaling factors of drag coefficients [nondim] + integer :: c + + ! The number and names of drag coefficients should match those of the streaming filters. + call get_param(param_file, mdl, "N_FILTERS", CS%nf, & + "Number of streaming band-pass filters to be used in the simulation.", & + default=0, do_not_log=.true.) + call get_param(param_file, mdl, "FILTER_NAMES", filter_name_str, & + "Names of streaming band-pass filters to be used in the simulation.", & + do_not_log=.true.) + + allocate(CS%coef_u(G%IsdB:G%IedB,G%jsd:G%jed,CS%nf)) ; CS%coef_u(:,:,:) = 0.0 + allocate(CS%coef_v(G%isd:G%ied,G%JsdB:G%JedB,CS%nf)) ; CS%coef_v(:,:,:) = 0.0 + allocate(filter_names(CS%nf)) ; read(filter_name_str, *) filter_names + + if (len_trim(wave_drag_file) > 0) then + do c=1,CS%nf + call get_param(param_file, mdl, "BT_"//trim(filter_names(c))//"_DRAG_U", & + var_names(1), "The name of the variable in BT_WAVE_DRAG_FILE "//& + "for the drag coefficient of the "//trim(filter_names(c))//& + " frequency at u points.", default="") + call get_param(param_file, mdl, "BT_"//trim(filter_names(c))//"_DRAG_V", & + var_names(2), "The name of the variable in BT_WAVE_DRAG_FILE "//& + "for the drag coefficient of the "//trim(filter_names(c))//& + " frequency at v points.", default="") + call get_param(param_file, mdl, "BT_"//trim(filter_names(c))//"_DRAG_SCALE", & + var_scale, "A scaling factor for the drag coefficient of the "//& + trim(filter_names(c))//" frequency.", default=1.0, units="nondim") + + if (len_trim(var_names(1))+len_trim(var_names(2))>0 .and. var_scale>0.0) then + call MOM_read_data(wave_drag_file, trim(var_names(1)), CS%coef_u(:,:,c), G%Domain, & + position=EAST_FACE, scale=var_scale*GV%m_to_H*US%T_to_s) + call MOM_read_data(wave_drag_file, trim(var_names(2)), CS%coef_v(:,:,c), G%Domain, & + position=NORTH_FACE, scale=var_scale*GV%m_to_H*US%T_to_s) + call pass_vector(CS%coef_u(:,:,c), CS%coef_v(:,:,c), G%domain, & + direction=To_All+SCALAR_PAIR) + + write(mesg, *) "MOM_wave_drag: ", trim(filter_names(c)), & + " coefficients read from file, scaling factor = ", var_scale + call MOM_error(NOTE, trim(mesg)) + endif ! (len_trim(var_names(1))+len_trim(var_names(2))>0 .and. var_scale>0.0) + enddo ! k=1,CS%nf + endif ! (len_trim(wave_drag_file) > 0) + +end subroutine wave_drag_init + +!> This subroutine calculates the sum of the products of the tidal velocities and the scaled +!! frequency-dependent drag for each tidal constituent specified in MOM_input. +subroutine wave_drag_calc(u, v, drag_u, drag_v, G, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(wave_drag_CS), intent(in) :: CS !< Control structure of MOM_wave_drag + !>@{ Tidal velocities from the output of streaming band-pass filters [L T-1 ~> m s-1] + real, dimension(:,:,:), pointer, intent(in) :: u, v + !>@} + !>@{ Sum of products of tidal velocities and scaled frequency-dependent drag [L2 T-2 ~> m2 s-2] + real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), intent(out) :: drag_u + real, dimension(G%isd:G%ied,G%JsdB:G%JedB), intent(out) :: drag_v + !>@} + + ! Local variables + integer :: is, ie, js, je, i, j, k + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + Drag_u(:,:) = 0.0 ; Drag_v(:,:) = 0.0 + + !$OMP do + do k=1,CS%nf ; do j=js,je ; do I=is-1,ie + Drag_u(I,j) = Drag_u(I,j) + u(I,j,k) * CS%coef_u(I,j,k) + enddo ; enddo ; enddo + + !$OMP do + do k=1,CS%nf ; do J=js-1,je ; do i=is,ie + Drag_v(i,J) = Drag_v(i,J) + v(i,J,k) * CS%coef_v(i,J,k) + enddo ; enddo ; enddo + +end subroutine wave_drag_calc + +!> \namespace mom_wave_drag +!! +!! By Chengzhu Xu (chengzhu.xu@oregonstate.edu) and Edward D. Zaron, December 2024 +!! +!! This module calculates the net effects of the frequency-dependent internal wave drag applied to +!! the tidal velocities, and returns the sum of products of frequency-dependent drag coefficients +!! and tidal velocities for each constituent to the MOM_barotropic module for further calculations. +!! It relies on the use of MOM_streaming_filter for determining the tidal velocities. Furthermore, +!! the number of drag coefficients cannot exceed that of the streaming filters, and the names of +!! drag coefficients should match those of the streaming filters. The frequency-dependent drag +!! coefficients are read from the same file for the linear drag coefficients in MOM_barotropic. + +end module MOM_wave_drag + From ed0707480b1420fbb159afbec9e53fbdfe7bda26 Mon Sep 17 00:00:00 2001 From: William Xu Date: Tue, 4 Feb 2025 14:58:47 -0400 Subject: [PATCH 2/2] Frequency-dependent internal wave drag Minor performance optimization. Also, this commit allows the filters to be activated without the frequency-dependent drag being activated. --- src/core/MOM_barotropic.F90 | 44 ++++++++++++------- .../lateral/MOM_wave_drag.F90 | 22 +++++----- 2 files changed, 40 insertions(+), 26 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 442d9960a6..ccb49e7036 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -609,9 +609,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! spacing [H L ~> m2 or kg m-1]. real, dimension(:,:,:), pointer :: ufilt, vfilt ! Filtered velocities from the output of streaming filters [L T-1 ~> m s-1] - real, dimension(SZIBW_(CS),SZJW_(CS)) :: Drag_u + real, dimension(SZIB_(G),SZJ_(G)) :: Drag_u ! The zonal acceleration due to frequency-dependent drag [L T-2 ~> m s-2] - real, dimension(SZIW_(CS),SZJBW_(CS)) :: Drag_v + real, dimension(SZI_(G),SZJB_(G)) :: Drag_v ! The meridional acceleration due to frequency-dependent drag [L T-2 ~> m s-2] real, target, dimension(SZIW_(CS),SZJW_(CS)) :: & eta, & ! The barotropic free surface height anomaly or column mass @@ -1433,9 +1433,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! Compute instantaneous tidal velocities and apply frequency-dependent drag. ! Note that the filtered velocities are only updated during the current predictor step, ! and are calculated using the barotropic velocity from the previous correction step. + if (CS%use_filter) then + call Filt_accum(ubt(G%IsdB:G%IedB,G%jsd:G%jed), ufilt, CS%Time, US, CS%Filt_CS_u) + call Filt_accum(vbt(G%isd:G%ied,G%JsdB:G%JedB), vfilt, CS%Time, US, CS%Filt_CS_v) + endif + if (CS%use_filter .and. CS%linear_freq_drag) 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) call wave_drag_calc(ufilt, vfilt, Drag_u, Drag_v, G, CS%Drag_CS) !$OMP do do j=js,je ; do I=is-1,ie @@ -1461,8 +1464,6 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, Drag_v(i,J) = 0.0 endif enddo ; enddo - else - Drag_u(:,:) = 0.0 ; Drag_v(:,:) = 0.0 endif if ((Isq > is-1) .or. (Jsq > js-1)) then @@ -2090,12 +2091,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & - ((Cor_v(i,J) + PFv(i,J)) - (vbt(i,J)*Rayleigh_v(i,J) + Drag_v(i,J))) + ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) enddo ; enddo else !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv-1,iev+1 - v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J) - Drag_v(i,J)) + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) enddo ; enddo endif @@ -2168,13 +2169,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & - ((Cor_u(I,j) + PFu(I,j)) - (ubt(I,j)*Rayleigh_u(I,j) + Drag_u(I,j))) + ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) enddo ; enddo !$OMP end do nowait else !$OMP do schedule(static) do j=jsv,jev ; do I=isv-1,iev - u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j) - Drag_u(I,j)) + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) enddo ; enddo !$OMP end do nowait endif @@ -2245,12 +2246,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * & - ((Cor_u(I,j) + PFu(I,j)) - (ubt(I,j)*Rayleigh_u(I,j) + Drag_u(I,j))) + ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j)) enddo ; enddo else !$OMP do schedule(static) do j=jsv-1,jev+1 ; do I=isv-1,iev - u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j) - Drag_u(I,j)) + u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j)) enddo ; enddo endif @@ -2334,13 +2335,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * & - ((Cor_v(i,J) + PFv(i,J)) - (vbt(i,J)*Rayleigh_v(i,J) + Drag_v(i,J))) + ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J)) enddo ; enddo !$OMP end do nowait else !$OMP do schedule(static) do J=jsv-1,jev ; do i=isv,iev - v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J) - Drag_v(i,J)) + v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J)) enddo ; enddo !$OMP end do nowait endif @@ -2645,6 +2646,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo endif + if (CS%use_filter .and. CS%linear_freq_drag) then ! Apply frequency-dependent drag + !$OMP do + do j=js,je ; do I=is-1,ie + u_accel_bt(I,j) = u_accel_bt(I,j) - Drag_u(I,j) + enddo ; enddo + !$OMP do + do J=js-1,je ; do i=is,ie + v_accel_bt(i,J) = v_accel_bt(i,J) - Drag_v(i,J) + enddo ; enddo + endif + if (id_clock_calc_post > 0) call cpu_clock_end(id_clock_calc_post) if (id_clock_pass_post > 0) call cpu_clock_begin(id_clock_pass_post) @@ -5008,10 +5020,12 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & endif ! CS%linear_wave_drag ! Initialize streaming band-pass filters and frequency-dependent drag - if (CS%use_filter .and. CS%linear_freq_drag) then + 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 + if (CS%use_filter .and. CS%linear_freq_drag) then if (.not.CS%linear_wave_drag .and. len_trim(wave_drag_file) > 0) then inputdir = "." ; call get_param(param_file, mdl, "INPUTDIR", inputdir) wave_drag_file = trim(slasher(inputdir))//trim(wave_drag_file) diff --git a/src/parameterizations/lateral/MOM_wave_drag.F90 b/src/parameterizations/lateral/MOM_wave_drag.F90 index 7ad340995e..a507c762c1 100644 --- a/src/parameterizations/lateral/MOM_wave_drag.F90 +++ b/src/parameterizations/lateral/MOM_wave_drag.F90 @@ -18,10 +18,9 @@ module MOM_wave_drag !> Control structure for the MOM_wave_drag module type, public :: wave_drag_CS ; private - integer :: nf !< Number of filters to be used in the simulation - !>@{ Spatially varying, frequency-dependent drag coefficients [H T-1 ~> m s-1] - real, allocatable, dimension(:,:,:) :: coef_u, coef_v - !>@} + integer :: nf !< Number of filters to be used in the simulation + real, allocatable, dimension(:,:,:) :: coef_u !< frequency-dependent drag coefficients [H T-1 ~> m s-1] + real, allocatable, dimension(:,:,:) :: coef_v !< frequency-dependent drag coefficients [H T-1 ~> m s-1] end type wave_drag_CS contains @@ -92,13 +91,14 @@ end subroutine wave_drag_init subroutine wave_drag_calc(u, v, drag_u, drag_v, G, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(wave_drag_CS), intent(in) :: CS !< Control structure of MOM_wave_drag - !>@{ Tidal velocities from the output of streaming band-pass filters [L T-1 ~> m s-1] - real, dimension(:,:,:), pointer, intent(in) :: u, v - !>@} - !>@{ Sum of products of tidal velocities and scaled frequency-dependent drag [L2 T-2 ~> m2 s-2] - real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), intent(out) :: drag_u - real, dimension(G%isd:G%ied,G%JsdB:G%JedB), intent(out) :: drag_v - !>@} + real, dimension(:,:,:), pointer, intent(in) :: u !< Zonal velocity from the output of + !! streaming band-pass filters [L T-1 ~> m s-1] + real, dimension(:,:,:), pointer, intent(in) :: v !< Meridional velocity from the output of + !! streaming band-pass filters [L T-1 ~> m s-1] + real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), intent(out) :: drag_u !< Sum of products of filtered velocities + !! and scaled frequency-dependent drag [L2 T-2 ~> m2 s-2] + real, dimension(G%isd:G%ied,G%JsdB:G%JedB), intent(out) :: drag_v !< Sum of products of filtered velocities + !! and scaled frequency-dependent drag [L2 T-2 ~> m2 s-2] ! Local variables integer :: is, ie, js, je, i, j, k