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

Frequency-dependent internal wave drag #722

Open
wants to merge 4 commits into
base: dev/gfdl
Choose a base branch
from

Conversation

c2xu
Copy link

@c2xu c2xu commented Sep 12, 2024

Utilizing streaming bandpass filters, the linear wave drag can now be applied to the M2 and K1 velocities independently, without affecting the broadband barotropic velocity fields. This implementation would be particularly useful (and perhaps necessary) for correctly representing tides in climate simulations, because if the linear drag is applied to the broadband barotropic velocity, the low-frequency, non-tidal motions would also be damped.

Utilizing streaming bandpass filters, the linear wave drag can now be
applied to the M2 and K1 velocities independently, without affecting
the broadband barotropic velocity fields.
Simplified the implementation by allowing both linear wave drag and
frequency-dependent drag to be turned on at the same time.
call MOM_read_data(wave_drag_file, m2_drag_u, CS%lin_drag_um2, G%Domain, &
position=EAST_FACE, scale=GV%m_to_H*US%T_to_s)
call pass_var(CS%lin_drag_um2, G%Domain)
CS%lin_drag_um2(:,:) = m2_drag_scale * CS%lin_drag_um2(:,:)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that it would be clearer and more efficient if you were to include m2_drag_scale in the scale= argument of the MOM_read_data call two lines earlier, so that it became scale=m2_drag_scale*GV%m_to_H*US%T_to_s.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I followed what has been done with lin_drag_h and wave_drag_scale if CS%linear_wave_drag is true. I made changes for both CS%linear_wave_drag and CS%linear_freq_drag.

if (CS%linear_freq_drag) then
!$OMP do
do j=js,je ; do I=is-1,ie
if (CS%lin_drag_um2(I,j) > 0.0 .or. CS%lin_drag_uk1(I,j) > 0.0) then
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It has been our experience that moving logical tests outside of do loops is almost always more efficient, even if it means duplicating a little code or (in this case) adding 0 to some expressions that did not need it. Compliers would do optimizations like this if they have enough information, but a compiler would not know, for instance, that there should never be negative values of lin_drag_um2 and lin_drag_uk1.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again I was following what has been done in the case of CS%linear_wave_drag. Here, it seems not possible to completely eliminate logical tests inside do loops, because the calculation should only be performed where Htot is non-zero. I did remove the two individual logical tests for CS%lin_drag_um2(I,j) > 0.0 and CS%lin_drag_uk1(I,j) > 0.0 in both loops though, and this should add some efficiency.

@@ -2054,7 +2103,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
do J=jsv-1,jev ; do i=isv-1,iev+1
vel_prev = vbt(i,J)
vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + &
dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J)))
dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J) - Drag_v(i,J)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As written, Drag_u and Drag_v do not appear to be changing within the barotropic time-stepping loops, so is there a reason why they could not be included in BT_force_u and BT_force_v?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point. The drag is essentially an external forcing.

@@ -2063,12 +2112,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))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here and elsewhere , Drag_u and Drag_v have been added to sums without any parentheses to specify the order of operations. Because floating point math is not commutative we very strongly encourage sums of 3 or more numbers to include enough parentheses to specific the order of arithmetic, as is discussed at https://github.com/NOAA-GFDL/MOM6/wiki/Code-style-guide. ((a + b) + c is does not always give the same value as a + (b + c); for example, (-1. + 1.) + 1e-20 = 1e-20 but -1. + (1. + 1e-20) = 0.0 when evaluated using IEEE 64-bit floating point numbers.)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did not realize this but it has been updated now.

@Hallberg-NOAA Hallberg-NOAA added enhancement New feature or request Parameter change Input parameter changes (addition, removal, or description) labels Sep 16, 2024
@c2xu c2xu marked this pull request as draft September 16, 2024 17:10
@c2xu c2xu force-pushed the c2xu/linear_freq_drag branch 2 times, most recently from 8fbeaf3 to 398b79c Compare September 16, 2024 17:26
Minor update for performance optimization.
@c2xu
Copy link
Author

c2xu commented Sep 20, 2024

I tested the new implementation and after a small bug fix, the results are identical as before.

@c2xu c2xu marked this pull request as ready for review September 20, 2024 23:55
@herrwang0
Copy link

It looks to me the newly-added Drag_[uv] forcing term does not change for each baroclinic time step. And as @Hallberg-NOAA suggested and the changes adopted since, it makes more sense to add it to BT_force_[uv] (vertical sum of baroclinic forcings).

Following the same logic, I kind of wanted to open an discussion on whether BT solver would be the best place to house the calculation of Drag_[uv]. In my opinion, here are the pros and cons for having it in BT solver.

  • Pros
  1. The filter as it stands is based on barotropic velocities and we can avoid repeating calculations of barotropic velocities.
  2. Drag_[uv] is essentially a 2D force on the whole column, so it is a barotropic force.
  • Cons
  1. The term is essentially updated every baroclinic time step. So there is no difference between doing the calculation inside and outside of BT solver. For BT solver, the term is no different from any other vertically-summed baroclinic forcing terms. I always think (and I could be wrong) BT solver should be reserved for calculating forcings that are time-varying for each BT step.
  2. If you look from a momentum budget perspective, the journey through BT solver makes Drag_[uv] embedded in [uv]_accel_bt, along with time-varying barotropic forcings, which will be diagnostics "[uv]_BT_accel". This seems counterintuitive and confusing to me.
  3. I wonder if it makes it difficult for potentially expansion of the functionality in the future. What if we want to apply it to non-barotropic velocity (say, bottom velocity)?

@c2xu
Copy link
Author

c2xu commented Sep 22, 2024

It looks to me the newly-added Drag_[uv] forcing term does not change for each baroclinic time step. And as @Hallberg-NOAA suggested and the changes adopted since, it makes more sense to add it to BT_force_[uv] (vertical sum of baroclinic forcings).

Following the same logic, I kind of wanted to open an discussion on whether BT solver would be the best place to house the calculation of Drag_[uv]. In my opinion, here are the pros and cons for having it in BT solver.

  • Pros
  1. The filter as it stands is based on barotropic velocities and we can avoid repeating calculations of barotropic velocities.
  2. Drag_[uv] is essentially a 2D force on the whole column, so it is a barotropic force.
  • Cons
  1. The term is essentially updated every baroclinic time step. So there is no difference between doing the calculation inside and outside of BT solver. For BT solver, the term is no different from any other vertically-summed baroclinic forcing terms. I always think (and I could be wrong) BT solver should be reserved for calculating forcings that are time-varying for each BT step.
  2. If you look from a momentum budget perspective, the journey through BT solver makes Drag_[uv] embedded in [uv]_accel_bt, along with time-varying barotropic forcings, which will be diagnostics "[uv]_BT_accel". This seems counterintuitive and confusing to me.
  3. I wonder if it makes it difficult for potentially expansion of the functionality in the future. What if we want to apply it to non-barotropic velocity (say, bottom velocity)?

Thanks for the comments, but I think ideally, for better accuracy, the tidal velocities [uv]m2 and [uv]k1, and hence Drag_[uv], should be calculated for each barotropic time step instead of baroclinic time step, though this seems to significantly increase the computational cost (about 20%-30% increase of computation time for one-layer, tide only simulations, though this could change depending on the configuration). There is also the complexity related to the prediction steps versus correction steps when calculating [uv]m2 and [uv]k1 using the current streaming filter module.

For my own configuration (one-layer, tide only model), the current implementation seems to work fine. But we should definitely keep this conversation open, as we might want to modify the implementation depending on the application. For example, in 3D models, we might want to apply the drag to certain layers instead of the barotropic velocity. And in simulations where the baroclinic time step is large, then the drag should be time-varying for each BT step.

Bug fix for memory allocation of uninitialized tidal velocities.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request Parameter change Input parameter changes (addition, removal, or description)
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants