Contribution of gravity waves to shear in the extratropical lowermost stratosphere: insights from idealized baroclinic life cycle experiments

The following files contains the code and ICON output processed for analysis and visualization for caeting our publication "Contribution of gravity waves on shear in the extratropical lowermost stratosphere: insights from idealized baroclinic life cycle experiments". 

Follwing things are included here:

All the python codes used for the anlaysis purpose.

The data required and the instructions on how to create the figures.Note: The data provided here is already processed for the corresponding figure. The eulerian ICON model output data is not provided here but available upon request if interested.

Figures in the publication. (figures.zip)

Summary and description on the python code provided in each .py files. (python_codes.zip)

LMS data: The LMS data is processed based on the definition " the region between the dynamical tropopause i.e. 3.5 pvu and the 380 K isentropic level, which corresponds to the maximum potential temperature of the tropical lapse rate tropopause. e.g., Holton et al., 1995; Shepherd, 2007. 

Details on calculations and use of python scripts

Calculate_Spectral_GWMF_S2_DCMIP_ml.py

This script is created to calculate gravity waves momentum flux (GWMF) components and absolute GWMF along with the perturbation vertical velocity (w') and shear pertubations from the model level output of ICON simulations, using a spectral filter to calculate u', v', w'. 

The small-scale horizontal divergence and pertubation vertical velocity calculated using prime quantities on the spectral domain.The vertical wind shear and GW momentum fluxes are calculated using:

S² = (du'/dz)² + (dv'/dz)² (in s⁻²)

|GWMF| = ρ* sqrt(\overline{(u'w')}² + \overline{(v'w')}²)

The calculated quantites will be saved in the netcdf format for further analysis. The ml output is then interpolated on the alitude levels to investigate w', GW momentum fluxes and shear on height focussing in the Lowermost stratosphere. (figure 6). Also, the LMS filter is applied to the calculated quantities for figure 11 and 12. For more details see Umbarkar et al. (2025).

Usage:  python Calculate_Spectral_GWMF_S2_DCMIP_ml.py

Calculate_NSR_ml.py

This script is created to calculate GW energetics such as Eddy kinetic energy, potential energy, total GWs energy along with instbilty variables such as static stability, vertical shear and gradient Richardson number from output ml file of ICON adiabtic idealized simulations. Note the static stability aka squared brunt vaisala frequency N^2 is already available as bvf2 in the ICON output files.

The vertical wind shear on full domain is calculated using:

S² = (du/dz)² + (dv/dz)² (in s⁻²)

The gradient Richardson number is given as:  Ri= N²/S²

              Usage: python Calculate_NSR_ml.py

Calculation of relative occurrence frequencies calculation of N² and S² are followed similarly for all sensitivity experiments as in  figure 7 and 8. The data used for plotting is limited to 5 km above and below part is skipped due to the focus region is UTLS and not the Boundary layer. Note that, this selection is arbitraty and in this case is based on the selection of lower Boundary.

Note that, for occurence frequency calculations and plotting , ICON output processed NetCDF files for each simulation to create individual figure named with *_fig*.nc  format. 

e.g. LMS data of REF wind experiment for day 12 used in fig09 named as REF_wind_LMS_filtered_timestep_12_fig09.nc

 Correlation of instability variables (figure 9 and 10)

Recreated as in ( kaluza2021 figure 14a)

The same LMs definition is applied for the turbulence diagnostic analysis i.e. for N²-S² pair plots (figure 9 and 10)For the occurrence frequency calculation, the bin sizes are selected based on the maxima's observed for each axis variables.For example, the maxima is different of included physics experiments and hence the size of grid volumes are defined accordingly in order to compare it with other sensitvity cases.e.g. the bin/box sizes of N² and S² are distributed into 50 bins each regardless the maxima whether it is 10/40 in case of S².

Usage: python plot_fig09.py

How-TO use:

Sripts to plot individual figure are given respective names. Assuming you want to create figure 09, simply read corresponding data files in the repective code  and run the script using python.

Funding

This project has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – TRR 301 – Project-ID 428312742, TPChange:  The Tropopause Region in a Changing Atmosphere (https://tpchange.de/)