8000 Dimensionally inconsistent stochastic streamfunction? · Issue #354 · NCAR/MOM6 · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

Dimensionally inconsistent stochastic streamfunction? #354

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.

Sign up for GitHub

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

Jump to bottom
Open
Hallberg-NOAA opened this issue Apr 23, 2025 · 9 comments
Open

Dimensionally inconsistent stochastic streamfunction? #354

Hallberg-NOAA opened this issue Apr 23, 2025 · 9 comments

Comments

@Hallberg-NOAA
Copy link

I believe that the expression for the stochastic streamfunction psi is in [L T-1 ~> m s-1] in the expression

  do k=1,GV%ke
    do J=grid%jscB-1,grid%jecB ; do I=grid%iscB-1,grid%iecB
      psi(I,J,k) = sqrt(0.25 * dt * max((CS%skeb_diss(i  ,j  ,k) + CS%skeb_diss(i+1,j+1,k)) + &
                                        (CS%skeb_diss(i  ,j+1,k) + CS%skeb_diss(i+1,j  ,k)), 0.) ) &
                                  * CS%skeb_wts(I,J)
    enddo ; enddo
  enddo

( https://github.com/NOAA-GFDL/MOM6/blame/484b917568f99d89bfe137b7b925fa265599d0b6/src/parameterizations/stochastic/MOM_stochastics.F90#L344-L350 )

However, when psi is used to set velocity anomalies (in [L T ~> m s-1] a few lines later in

    do j=grid%jsc,grid%jec ; do I=grid%iscB,grid%iecB
      ustar(I,j,k) = - (psi(I,J,k) - psi(I,J-1,k)) * CS%taperCu(I,j) * grid%IdyCu(I,j)
      uc(I,j,k) = uc(I,j,k) + ustar(I,j,k)
    enddo ; enddo
    do J=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec
      vstar(i,J,k) =   (psi(I,J,k) - psi(I-1,J,k)) * CS%taperCv(i,J) * grid%IdxCv(i,J)
      vc(i,J,k) = vc(i,J,k) + vstar(i,J,k)
    enddo ; enddo

psi would have to be in units of [L2 T-1 ~> m2 s-1] to be dimensionally consistent.

@iangrooms, I believe that you originally wrote this block of code. Do you agree with my analysis of this situation?

I do not have a fix for this issue, but I have added potentially relevant scaling factors elsewhere in MOM_stochastics.F90 on my incomplete branch at https://github.com/Hallberg-NOAA/MOM6/tree/stochastics_unit_conversion that might be helpful in sorting this out.
@iangrooms
Copy link

The random field CS%skeb_wts in the first block of code has dimensions [L ~> m]. In the first block of code this makes the dimensions of psi be [L2 T-1 ~> m2 s-1] as necessary. At the moment this is pretty opaque because the declaration of skeb_wts doesn't include its dimensions. It's explained in our manuscript, but that isn't published yet. Should I add the dimensions to the declaration on dev/ncar and send it in with the next PR from NCAR? It may be a few months before the next pull from NCAR to main.

@iangrooms
Copy link

I see that the diagnostic output for skeb_wts says units are 'None' which is incorrect. All the diagnostics in that module have 'None' units, so that needs to be fixed.

@Hallberg-NOAA
Copy link
Author

Thank you very much for that clarification. If you would like, I can update my branch related to this issue with this added information about the correct units about skeb_wts, along with the missing conversion arguments that drew my attention to this code in the first place.

@iangrooms
Copy link

Yes, please. Thanks for your attention on this code!

@Hallberg-NOAA
Copy link
Author

I have updated the revisions to the MOM6 stochastics code on my branch at https://github.com/Hallberg-NOAA/MOM6/tree/stochastics_unit_conversion to reflect this conversation. I actually went a bit further than I had originally intended, as I also added the option to use rotationally symmetric expressions in the stochastics code. This all compiles, and by default it is supposed to be bitwise identical with what was there before (at least when dimensional consistency testing is not being used) but as we do not have any cases at GFDL that utilize the stochastics code, I would appreciate some help in testing this out (and perhaps suggesting further revisions), if you would be willing to do so @iangrooms.

There is still one point here that has me confused, though. I understand how run_stochastic_physics_ocn() could set CS%skeb_wts to be appropriate horizontal length scales (in [m]), but I don't see where US%m_to_L is being passed into run_stochastic_physics_ocn() or update_stochastics() so that CS%skeb_wts could be set in rescaled units of [L ~> m]. Is this somehow being done indirectly via a grid spacing or some other variable length that itself is in units of [L ~> m]'? Does CS%skeb_wtsneed to be rescaled after it is set externally? (All this would probably be clearer if I were not looking at the stub version of the relevant code inMOM6/config_src/external/stochastic_physics/stochastic_physics.F90`.)

@iangrooms
Copy link

Thanks for doing this! It was on my to-do list to make the stochastic code rotationally symmetric, so you've saved me the effort. I'm happy to run tests at NCAR.

Developing the stochastic backscatter has been a joint effort; I focused on the MOM6 side while @pjpegion and @nniraj123 worked on the stochastic physics package, which can be found at github.com/NOAA-PSL/stochastic_physics. So I'm not an expert on that code. But I think that code assumes MKS units for everything. The grid information is passed into stochastic_physics by MOM6, in degrees E and N.1 The Earth's radius is hard-coded in the package here using MKS units.

Given the above, is the correct way to do the dimensional consistency by explicitly converting the output of the stochastic physics package from MKS to L ([m ~> L]) before using it in MOM6?

1MOM_stochastics.F90 passes the lat & lon at h and q points to the stochastic physics package in the 'init' subroutine. I noticed that the comments in the declarations of geoLonT and geoLatT say 'at q points' instead of 'at h points' which is probably a copy-paste error.

@Hallberg-NOAA
Copy link
Author
Hallberg-NOAA commented Apr 29, 2025 •
edited
Loading

Looking at the code, I suspect that the stochastic physics code does use MKS units. Moreover, because it is advancing the field in time to get the autocorrelation, I think that it would be a bad idea to try rescaling the field that is being repeatedly passed to the physics. I think that it is best to modify the MOM6 code to do the rescaling where the skeb_wts are used and leave that array in units of [m]. I have updated my branch accordingly.

Thanks also for pointing out the incorrect descriptions of geoLonT and geoLatT, which are indeed at h points.

@iangrooms
Copy link

The code appears to be failing the dimensional scaling test for time. The stochastic physics package has its own input namelist where we specify the decorrelation time of the stochastic forcing in MKS units. MOM6 passes the time step size dt to the stochastic physics package on this line. I suspect that is the source of the dimensional scaling test failure. I'm not sure the best way to fix this; perhaps multiply dt by the appropriate scaling factor in the call to init_stochastic_physics_ocn?

@iangrooms
Copy link

I added unit scaling for dt, as alluded to in your comment about autocorrelation, and it now passes all the dimensional scaling tests. Would you rather I send this as a PR to your branch, or copy your code and submit to NCAR (or something else)? I also added an option to pre-calculate the weights in the smoothing, which should save computation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@Hallberg-NOAA @iangrooms
0