Frequency-dependent internal wave drag

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(SZIB_(G),SZJ_(G)) :: Drag_u
+                  ! The zonal acceleration due to frequency-dependent drag [L T-2 ~> m s-2]
+  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
                   ! anomaly [H ~> m or kg m-2]
@@ -1422,6 +1430,42 @@ 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) 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 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
+  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 +1637,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
@@ -2609,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)
@@ -4710,9 +4758,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,12 +5019,21 @@ 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
+  ! Initialize streaming band-pass filters and frequency-dependent drag
   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)
+      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
 
   dtbt_restart = -1.0
@@ -5302,7 +5363,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.)

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<now) then
+    dt = now - CS%old_time
+    CS%old_time = now
+
+    do k=1,CS%nf
+      c1 = CS%filter_omega(k) * dt
+      c2 = 1.0 - CS%filter_alpha(k) * c1
+
+      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
+    enddo ! k=1,CS%nf
+  endif ! (CS%old_time<now)
 
-  dt = now - CS%old_time
-  CS%old_time = now
-
-  ! Timestepping
-  do k=1,CS%nf
-    c1 = CS%filter_omega(k) * dt
-    c2 = 1.0 - CS%filter_alpha(k) * c1
-
-    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
-  enddo
   u1 => CS%u1
 
 end subroutine Filt_accum

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 <MOM_memory.h>
+
+!> 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
+  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
+
+!> 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
+  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
+
+  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
+