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.

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

Potentially incorrect computation of "lagrangian_tendency_of_air_pressure" #40

Closed
nbren12 opened this issue May 5, 2020 · 13 comments
Closed

Potentially incorrect computation of "lagrangian_tendency_of_air_pressure" #40

nbren12 opened this issue May 5, 2020 · 13 comments

Comments

@nbren12
Copy link

nbren12 commented May 5, 2020 •
edited
Loading

cc @lharris4

Just to follow up on a zoom conversation, I am concerned that the non-hydrostatic model computes "omga" incorrectly. For reference, it is my understanding that CF standard name for "omga" is "lagrangian_tendency_of_air_pressure".

As far as I can tell, “omga” (dp/dt with units Pa/s) is computed two places in the model.

For hydrostatic runs, it computes it as "omga = dp/dt= \partial p/\partial t + (u,v,w).grad(p)" here:

GFDL_atmos_cubed_sphere/model/dyn_core.F90

Line 1102 in eee6555

call adv_pe(ua, va, pem, omga, gridstruct, bd, npx, npy, npz, ng)

I think this formula is correct.

For non-hydrostatic, it is computed with the formula omga=delp/delz * w here:

GFDL_atmos_cubed_sphere/model/fv_dynamics.F90

Line 584 in eee6555

omga(i,j,k) = delp(i,j,k)/delz(i,j,k)*w(i,j,k)

I think this formula is a slightly incorrect application of the change of variables formula from Kasahara (1974):
Screen Shot 2020-05-05 at 1 49 37 PM
I think FV3 is missing the two rightmost terms in this formula.

Is my reading of the code correct here? The missing terms may not be large in magnitude, but could be important for detailed budget analyses using "omga" from this model, especially over steeply sloped topography.
@lharris4
Copy link
Contributor

lharris4 commented May 6, 2020 via email

@nbren12
Copy link
Author

nbren12 commented May 6, 2020

I am sorry but this equation makes no physical sense as a diagnostic.

Hmm. The formula "omega=w * delp/delz" is not the correct change of variables formula, since it assumes that p=p(z), when it is actually p=p(x,y,z,t). In other words, it assume that vertical velocity and pressure velocity point in the same direction. This is not true when pressure surfaces vary horizontally and in time, in which case the pressure velocity will include some component of the horizontal velocity in height coordinates. It does not matter if "w" is prognostic or not. The time and space derivatives in eq 3.12 account for the fact that the unit vector perpendicular to a pressure surface is not parallel to the unit vector perpendicular to a height surface.

One possible source of the disagreement between the hydrostatic and
nonhydrostatic calculations is that the pressure surfaces...

I think we can get around this be defining omega=dp^*/dt where p^* is the hydrostatic pressure. This way, "omega" could be interpreted as a mass flux.

@nbren12
Copy link
Author

nbren12 commented May 6, 2020

Please note that the derivation of Eq. 3.12 makes no assumption about the governing equations or hydro-static balance. It's all chain rule.

@lharris4
Copy link
Contributor

lharris4 commented May 7, 2020 via email

@ofuhrer
Copy link
Contributor

ofuhrer commented May 7, 2020

My $0.02:

I think the "formal" definition of \omega is effectively - as Noah points out - the change of pressure in an air parcel over time, i.e. \omega = dp/dt. This is a natural choice as a prognostic variable for hydrostatic models in pressure coordinates.

But, in non-hydrostatic atmospheric models it is typically only a diagnostic to feed into the physics (and w is the prognostic variable). The reason for this is that physical parametrization were often initially written for hydrostatic models with pressure based coordinates. In non-hydrostatic models it is thus often approximated as \omega = -\rho g w (which is what the formula in FV3 is doing for the non-hydrostatic option, and which is a good approximation on the synoptic scale), and simply serves as an approximation of vertical velocity (or vertical mass flux, as Lucas pointed out).

If you need an an exact computation of dp/dt, \omega is probably not the way to go.

@nbren12
Copy link
Author

nbren12 commented May 7, 2020

Thanks @ofuhrer and @lharris4. At this point, I assume we all agree that (please comment if not)

  1. The quantity of interest is dp^*/dt, which is proportional to the mass flux across surface of p^*, where p^* is the hydrostatic component of pressure.
  2. The formula used by the non-hydrostatic model dp^*/dt = delp/delz w is a convenient approximation, but not strictly accurate. There are two exact formulas to computes this. The first is to compute the total/material derivative of p^* in (x,y,z) coordinates as is done by the hydrostatic model. The second is Eq 3.12 of Kasahara (derived from writing w=dz/dt in (x,y,p^*) coordinates and some simple algebra). delp/delz*w is an approximation of this formula which leaves of some terms.

As I said in the opening comment, I doubt this approximation causes large errors on the synoptic scale, but it could be a problem for detailed budget analysis where we need to e.g. close the mass budget of a control volume between two p^* surfaces.

Moreover, Chris B pointed out to me that w vanishes near the ocean surface, but dp^*/dt does not, so the approximate formula could poorly track near-surface mass fluxes.

If you need an an exact computation of dp/dt, \omega is probably not the way to go.

Sorry, which "\omega" are you referring to? I think we must be agreeing...

@lharris4
Copy link
Contributor

lharris4 commented May 8, 2020 via email

@nbren12
Copy link
Author

nbren12 commented May 10, 2020 •
edited
Loading

Both sort of...we are regridding the “omga” diagnostic of the model from model to pressure levels.

However, the choice of vertical discretization shouldn’t matter. omega = dp*/dt is a continuous quantity defined to be the Lagrangian derivative of hydrostatic pressure surfaces. It does not matter what grid it is discretized on. This material derivative is just tricker to compute in a non-pressure coordinate system. It’s pretty similar to how w=dz/dt is the Lagrangian derivative of height regardless of discretization used by the model.

I hope we can agree that omega is not exactly “-rho g * w”. As confirmed by @ofuhrer and conversations with Chris B, this is an approximation which is only exact in certain conditions (anelastic and hydrostatic).

@nbren12 nbren12 mentioned this issue Jun 5, 2020
Merged
nbren12 added a commit to ai2cm/fv3net that referenced this issue Jun 5, 2020
@nbren12
Redefine fine-res budget for new ShiELD (#383)
4cd497f 
The new simulation had different physics options and therefore different diagnostics outputs. Relevant changes included:
- no gfdl microphysics output. Now microphysics is in the physics tendency, and the fv_sat_adj dynamical core phase changes
- New formulation of pressure velocity "omega" resolving this error: NOAA-GFDL/GFDL_atmos_cubed_sphere#40. The correct pressure velocity and corresponding eddy flux are now named "vulcan_omega" and "eddy_flux_vulcan_omega_sphum" for example.

Changes introduced by this commit:
- Adjust hardcoded budget terms to accomodate the changes above. The team suggested that I not overgeneralize this code since won't have a new dataset anytime soon.
- Update test schema to reflect the new training data
- Add unit and long_name information to the outputs
- Simplify the variable renaming logic. Previously "_coarse" was automatically stripped from all the names of the diagnostic output, which made it tricky to know infer what variables were available at a particular step of the pipeline. Now most of the renaming is concentrated into one function.
- Adjust `FineResolutionSources` to reflect this schema change.

[VCMML-310] 

Closes #347

[VCMML-310]: https://vulcan.atlassian.net/browse/VCMML-310
climbfuji pushed a commit to climbfuji/GFDL_atmos_cubed_sphere that referenced this issue Oct 1, 2020
@DusanJovic-NOAA
Merge pull request NOAA-GFDL#40 from NOAA-GFDL/dev/emc
2ec76f8 
Merge 'regional_boundary_update' from upstream
@lharris4 lharris4 mentioned this issue Jan 15, 2021
Merged
@spencerkclark
Copy link
Member

Through #184 I believe the latest changes I made in SHiELD to address this issue were ported here, so I think this can be closed. The commit history appears to be squashed, but I am posting the description of my relevant merge request from GFDL's internal GitLab. This merge request made the "lagrangian_tendency_of_hydrostatic_pressure" the default "omega" diagnostic, and added the option to pass it to the physics if desired.

Thanks to @nbren12 for identifying this issue and providing a suggestion initially for how to fix it in the code.

Always compute omega diagnostic via the method that is used in hydrostatic mode

This MR refactors the logic in dyn_core.F90, fv_mapz.F90, and fv_dynamics.F90 for computing Atm%omga such that it is computed in the same manner in both hydrostatic and non-hydrostatic modes.

Implications of this MR:

  • The "lagrangian_tendency_of_hydrostatic_pressure" diagnostic is now removed, as it now is identical to the "omega" diagnostic. Previously, the "omega" diagnostic in non-hydrostatic mode represented the "local omega" computed by delp / delz * w.
  • A namelist flag has been added to fv_core_nml which allows you to pass the "full omega" to the physics from the dynamical core in non-hydrostatic mode. This flag is called pass_full_omega_to_physics_in_non_hydrostatic_mode; by default it is set to .false. to ensure bitwise reproduction of previous results.
  • Since we now compute omga in the same way in hydrostatic and non-hydrostatic mode, we can remove the do_omega argument to the Lagrangian_to_Eulerian subroutine, as it is now redundant with the last_step argument.

I have been careful to check that:

  • The "omega" diagnostic with these changes exactly matches the "lagrangian_tendency_of_hydrostatic_pressure" diagnostic from before. Additionally, the coarse diagnostics that depended on "lagrangian_tendency_of_hydrostatic_pressure", e.g. the vertical eddy flux diagnostics, are also identical.
  • With the same namelist parameter settings, the model produces bitwise identical results for all other fields; this was determined by looking at the md5sum of the restart files produced at the end of a one hour simulation. This was ensured by the fact that we still compute and pass the "local omega" to the physics in non-hydrostatic mode. Update 2021-06-21: I have now verified that this is true in the hydrostatic case as well.
  • As expected, the "omega" diagnostic in hydrostatic mode is bitwise identical to what it was before these code changes.

@lharris4
Copy link
Contributor

lharris4 commented Jun 29, 2022 via email

@nbren12
Copy link
Author

nbren12 commented Jun 29, 2022

Thanks for that @lharris4. Sounds like the issues is resolved. Closing.

@nbren12 nbren12 closed this as completed Jun 29, 2022
@nbren12
Copy link
Author

nbren12 commented Oct 11, 2022 via email

@lharris4
Copy link
Contributor

lharris4 commented Oct 11, 2022 via email

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
4 participants
@nbren12 @ofuhrer @spencerkclark @lharris4