search

NOAA-GFDL / GFDL_atmos_cubed_sphere

The GFDL atmos_cubed_sphere dynamical core code
Other
55 stars 118 forks source link

Potentially incorrect computation of "lagrangian_tendency_of_air_pressure" #40

Closed nbren12 closed 2 years ago

nbren12 commented 4 years ago

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: https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/blob/eee655598d58836a600c52a849e2ed69ab910653/model/dyn_core.F90#L1102 I think this formula is correct.

For non-hydrostatic, it is computed with the formula omga=delp/delz * w here: https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/blob/eee655598d58836a600c52a849e2ed69ab910653/model/fv_dynamics.F90#L584 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 commented 4 years ago

Hi, Noah. I am sorry but this equation makes no physical sense as a diagnostic. If we have a prognostic equation for vertical velocity in the transformed coordinate, it should certainly include the advective term, and that is indeed what Kasahara is doing. However we are not doing this: we are instead taking our (prognostic) vertical velocity and computing the pressure velocity from that.

One possible source of the disagreement between the hydrostatic and nonhydrostatic calculations is that the pressure surfaces (whether mass as in the hydrostatic system or the full pressure in the nonhydrostatic system) themselves are raising and lowering and are not fixed values, and that may be what Kasahara is trying to include here. Note that using ps and the ak/bk coordinates may not give you the exact hydrostatic pressure on the layer interfaces (p_{K+1/2} = pT + \sum{k=1}^{K} delp), since we use a full-mass coordinate and the physics changes the mass (through precipitation processes, vertical moisture and hydrometeor transport, etc.).

Lucas

On Tue, May 5, 2020 at 4:55 PM Noah D Brenowitz notifications@github.com wrote:

cc @lharris4 https://github.com/lharris4

Just to follow up on a zoom conversation, I am concerned that the way method used to compute "omga" is not correct for the non-hydrostatic model. 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: https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/blob/eee655598d58836a600c52a849e2ed69ab910653/model/dyn_core.F90#L1102 I think this formula is correct.

For non-hydrostatic, it is computed with the formula omga=delp/delz * w here: https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/blob/eee655598d58836a600c52a849e2ed69ab910653/model/fv_dynamics.F90#L584 I think this formula is a slightly incorrect application of the change of variables formula from Kasahara (1974) https://journals.ametsoc.org/doi/abs/10.1175/1520-0493%281974%29102%3C0509%3AVVCSUF%3E2.0.CO%3B2 : [image: Screen Shot 2020-05-05 at 1 49 37 PM] https://user-images.githubusercontent.com/1386642/81114691-61508080-8ed7-11ea-997c-d66ba8a87eaf.png 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.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/issues/40, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVAPOX4CP4SRTW3DIZ3RQB4K3ANCNFSM4MZ42IIA .

nbren12 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

Hi, Noah. Thanks. I was thinking myself that it would make more sense to think of omega as a mass flux rather than as an absolute velocity. For example, if the entire column from the surface to 500 mb expands thermally, at 500 mb w > 0 but \omega = 0.

Lucas

On Wed, May 6, 2020 at 4:22 PM Noah D Brenowitz notifications@github.com wrote:

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.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/issues/40#issuecomment-624869871, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVEDWLWOVBT6Q6ZITODRQHBJRANCNFSM4MZ42IIA .

ofuhrer commented 4 years ago

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 commented 4 years ago

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^, wherep^*` 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 commented 4 years ago

Hi, Noah. I think it depends on what levels we are working on. Are you doing your analysis on model levels or on fixed pressure levels?

Lucas

On Thu, May 7, 2020 at 12:12 PM Noah D Brenowitz notifications@github.com wrote:

Thanks @ofuhrer https://github.com/ofuhrer and @lharris4 https://github.com/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/delzw 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...

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/issues/40#issuecomment-625351344, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVBUKXA4PNYHKRS4EI3RQLMX7ANCNFSM4MZ42IIA .

nbren12 commented 4 years ago

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).

spencerkclark commented 2 years ago

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:

I have been careful to check that:

lharris4 commented 2 years ago

Hi, Spencer. Thanks for taking care of this. The nonhydrostatic omega is a nice addition to the code.

@nbren12 thanks for pushing ahead on this, I now see you are right. I added your discussion (with attribution) to the FV3 Documentation when I wrote it last year.

Lucas

On Wed, Jun 29, 2022 at 11:10 AM Spencer Clark @.***> wrote:

Through #184 https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/pull/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 https://github.com/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.

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/issues/40#issuecomment-1170102114, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVGMADL52J2LH4AHE3TVRRRM3ANCNFSM4MZ42IIA . You are receiving this because you were mentioned.Message ID: @.***>

nbren12 commented 2 years ago

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

nbren12 commented 1 year ago

Feel free to close the issue if you think it has been resolved. You are more familiar with this code than I am. Thanks for all your work on this!

On Wed, Jun 29, 2022 at 8:10 AM Spencer Clark @.***> wrote:

Through #184 https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/pull/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 https://github.com/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.

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/issues/40#issuecomment-1170102114, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAKSREUY7TWP2CHCH4TG6VTVRRRMZANCNFSM4MZ42IIA . You are receiving this because you were mentioned.Message ID: @.***>

lharris4 commented 1 year ago

Hi, all. I believe we now compute omega the same way in both hydrostatic and nonhydrostatic dynamics, which as Noah has explained is the correct definition of omega regardless of hydrostaticity.

Lucas

On Tue, Oct 11, 2022 at 5:02 AM Noah D. Brenowitz @.***> wrote:

Feel free to close the issue if you think it has been resolved. You are more familiar with this code than I am. Thanks for all your work on this!

On Wed, Jun 29, 2022 at 8:10 AM Spencer Clark @.***> wrote:

Through #184 https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/pull/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 https://github.com/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.

— Reply to this email directly, view it on GitHub < https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/issues/40#issuecomment-1170102114 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AAKSREUY7TWP2CHCH4TG6VTVRRRMZANCNFSM4MZ42IIA

. You are receiving this because you were mentioned.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/issues/40#issuecomment-1274354562, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVDKGPWTTKC4MKXXZJLWCUUK5ANCNFSM4MZ42IIA . You are receiving this because you were mentioned.Message ID: @.***>