Variable gravity

[WChen-NOAA]

Description

The model uses the shallow atmosphere approximation with constant gravity
This is problematic for a space weather model (WAM) where Ztop/r0 )~10%
Implement a scaling factor wherever gravity is applied to account for variations with height:

g(r0/r)^2 while setting var_grav = .true. in the name list
the new gravity as grav_var(i,j,k)

Fixes # (394)

How Has This Been Tested?

The test results are available vi doc

Checklist:

Please check all whether they apply or not

• My code follows the style guidelines of this project
• I have performed a self-review of my own code
• I have commented my code, particularly in hard-to-understand areas
• I have made corresponding changes to the documentation
• My changes generate no new warnings
• Any dependent changes have been merged and published in downstream modules

[lharris4]

Thank you for submitting this PR. I think this is a nice addition to FV3 and an important first step towards deep atmosphere dynamics. The implementation is done well, mostly self-contained and unobtrusive, and largely follows proper FV3 coding standards. I have a few comments about specific implementation details.

Think carefully about how frequently to update var_grav. The Lagrangian vertical coordinate does deform every acoustic timestep, but the variation in gravity g*(r/r0)**2 even over a dt_atmos dynamical timestep will be weak: for an extreme 100 m/s updraft and a 180 s timestep, the elevation change is only 18 km representing a 0.05% change in gravity. My first suggestion would be to simply re-calculate grav_var at the beginning of each dynamics timestep and hold it constant through the acoustic and remapping timesteps; this would reduce the amount of re-calculation and data movement, but stability and accuracy would need to be tested.

We also need to verify reproducibility for both hydrostatic and nonhydrostatic cases when var_grav == .F..

[lharris4]

This is the best place to put the calculation of variable gravity. We only should recalculate at the beginning of a dynamics timestep dt_atmos and not every remapping or acoustic timestep, as we can expect the variation of grav due to deformation of the Lagrangian vertical coordinate to be small during a dynamical timestep.

[lharris4]

I would strongly suggest defining grav_var(:,:,:) as a variable in fv_atmos_type instead of as a module variable in atmosphere_mod. Also, is there a need to have grav_var_h available outside of fv_dynamics()?

[WChen-NOAA]

The code has been changed as recommended (seen in updated fv_arrys.F90 and atmosphere.F90)

[lharris4]

Is there a good case to be made for re-calculating grav_var after every remapping? See my comment above. Do also note that since FV3 uses a vertically-Lagrangian coordinate, the geometric height of each interface will change on each acoustic timestep.

[WChen-NOAA]

I guess it is better to update grav-var after every remapping. The magnitude of both wind speed (could be hundred m/s to thousand m/s) and temperature (thousand degree) are very high. And species density changes rapidly with height and magnitude of geopotential height deviations increase with increasing altitude, which could reach ~3km at 800km. The WAM needs to coupled with an empirical model of Ionosphere Plasmasphere Electrodynamics. A relative precise scale height would get a better performance of the solar activities. Without this implementation, the scale height needs a correction before coupling to IPE while with current implementation, model turns a good result and allows the correction of the scale height to be removed.

[lharris4]

After some discussion, updating on the remapping timestep is probably the best.

[lharris4]

Why does this recalculation need to be done?

[WChen-NOAA]

To my understanding, this is not recalculation. The author only passed the updated gravity at interface into "dyn-core" and "grav_var" is taken as a local variable seen at line 184 and line 309.

[WChen-NOAA]

by the way, the grav_var_h updated at each remapping loop

[lharris4]

OK. I presume this is just to remap grav_var at layer interfaces onto layer centroids, done once per remapping loop.

[lharris4]

This loop can be refactored for better locality and cache utilization. I would recommend computing the surface (npz+1) values in the first loop, then compute from npz to 1 by -1 using a new k-j-i loop. This would also allow you to define newrad as a 2D array, reducing memory usage and possibly better fitting in cache.

[WChen-NOAA]

done

[lharris4]

Try to avoid divisions whenever possible; changing /2.0 to *0.5 is more efficient.

[WChen-NOAA]

changed

[lharris4]

Also replace with *0.5

[WChen-NOAA]

changed

[lharris4]

Good catch. Do note in the PR message that there are some bugfixes for MULTI_GASES (and for molecular_diffusion too).

[lharris4]

This line can be replaced with rdg(i,j,k) = -rdgas*rgrav(i,j,k) to avoid divisions.

[WChen-NOAA]

Thanks and replaced with suggestion

[yangfanglin]

Lucas,

Thank you for taking the time to carefully review this PR. One thing comes to my mind is that wind speed in the thermosphere is normally a few hundred m/s but can reach thousands during geomagnetic storms. Updating var_grav at each acoustic timestep may become necessary. @XiaqiongZhou-NOAA @WChen-NOAA

Think carefully about how frequently to update var_grav. The Lagrangian vertical coordinate does deform every acoustic timestep, but the variation in gravity g*(r/r0)**2 even over a dt_atmos dynamical timestep will be weak: for an extreme 100 m/s updraft and a 180 s timestep, the elevation change is only 18 km representing a 0.05% change in gravity. My first suggestion would be to simply re-calculate grav_var at the beginning of each dynamics timestep and hold it constant through the acoustic and remapping timesteps; this would reduce the amount of re-calculation and data movement, but stability and accuracy would need to be tested.

[lharris4]

One thing comes to my mind is that wind speed in the thermosphere is normally a few hundred m/s but can reach thousands during geomagnetic storms. Updating var_grav at each acoustic timestep may become necessary.

OK. While the wind speed may necessitate a small timestep, I don't know if this would necessarily translate into rapid variations in altitude that would lead to significant local changes in gravity. However again this would have to be determined through testing. The simplest approach would indeed be to simply update each timestep; the method @WChen-NOAA uses is not particularly expensive.

[WChen-NOAA]

Find a small err in one subroutine, and made a correction (works while using subroutine compute_total_energy). See the differences in attached file.
var_grav_PR.pdf

[lharris4]

Find a small err in one subroutine, and made a correction (works while using subroutine compute_total_energy). See the differences in attached file. var_grav_PR.pdf

Good catch---thanks for finding and fixing this. Thanks also for uploading the demonstration of this improvement as a world-visible GitHub attachment.

[XiaqiongZhou-NOAA]

It seems grav in ~/tools is still only constant. Just wondering if this impacts restart, diagnostics and others?

[WChen-NOAA]

@XiaqiongZhou-NOAA, Good question. I am checking again and think over it.

[WChen-NOAA]

@XiaqiongZhou-NOAA
After many tests and found current implementation of variable gravity does not has impact on the restart for L127 settings.
And found the restart issue for L149 setting is not due to the variable gravity implementation but may caused by early implementations for WAM, but this issue can be resolved with the full WAM implementation.
And I am working on the hydra settings for L127.

It seems grav in ~/tools is still only constant. Just wondering if this impacts restart, diagnostics and others?