Last updated on Nov 5, 2023

Demo script for the analyses done in Nakamura and Huang (2018, Science)

This is a complimentary demo script that can be used to implement the local wave activity, fluxes and flux convergence/divergence computation required in the analyses presented in Nakamura and Huang, Atmospheric Blocking as a Traffic Jam in the Jet Stream. Science. (2018)

This notebook demonstrate how to compute local wave activity and all the flux terms in equations (2) and (3) in NH2018 with the updated functionality in the python package falwa. To run the script, please install the package falwa using

python setup.py develop

after cloning the GitHub repo.

Please raise an issue in the GitHub repo or contact Clare S. Y. Huang (csyhuang@uchicago.edu) if you have any questions or suggestions regarding the package.

import datetime as dt
from math import pi
import numpy as np
from numpy import dtype
import xarray as xr
import matplotlib.pyplot as plt
%matplotlib inline
from falwa.oopinterface import QGFieldNH18
import falwa.utilities as utilities

Load ERA-Interim reanalysis data retrieved from ECMWF server

The sample script in this directory download_example.py include the code to retrieve zonal wind field U, meridional wind field V and temperature field T at various pressure levels. Given that you have an account on ECMWF server and have the ecmwfapi package installed, you can run the scripts to download data from there:

python download_example.py

u_file = xr.open_dataset('2005-01-23_to_2005-01-30_u.nc')
v_file = xr.open_dataset('2005-01-23_to_2005-01-30_v.nc')
t_file = xr.open_dataset('2005-01-23_to_2005-01-30_t.nc')
ntimes = u_file.time.size
time_array = u_file.time

u_file

<xarray.Dataset>
Dimensions:    (longitude: 240, latitude: 121, level: 37, time: 32)
Coordinates:
  * longitude  (longitude) float32 0.0 1.5 3.0 4.5 ... 354.0 355.5 357.0 358.5
  * latitude   (latitude) float32 90.0 88.5 87.0 85.5 ... -87.0 -88.5 -90.0
  * level      (level) int32 1 2 3 5 7 10 20 30 ... 850 875 900 925 950 975 1000
  * time       (time) datetime64[ns] 2005-01-23 ... 2005-01-30T18:00:00
Data variables:
    u          (time, level, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2018-07-17 16:50:39 GMT by grib_to_netcdf-2.8.0: grib_to_ne...

xarray.Dataset

• Dimensions:
  • longitude: 240
  • latitude: 121
  • level: 37
  • time: 32
• Coordinates: (4)
  • longitude
    (longitude)
    float32
    0.0 1.5 3.0 ... 355.5 357.0 358.5

    units :

    degrees_east

    long_name :

    longitude

    array([  0. ,   1.5,   3. , ..., 355.5, 357. , 358.5], dtype=float32)

  • latitude
    (latitude)
    float32
    90.0 88.5 87.0 ... -88.5 -90.0

    units :

    degrees_north

    long_name :

    latitude

    array([ 90. ,  88.5,  87. ,  85.5,  84. ,  82.5,  81. ,  79.5,  78. ,  76.5,
            75. ,  73.5,  72. ,  70.5,  69. ,  67.5,  66. ,  64.5,  63. ,  61.5,
            60. ,  58.5,  57. ,  55.5,  54. ,  52.5,  51. ,  49.5,  48. ,  46.5,
            45. ,  43.5,  42. ,  40.5,  39. ,  37.5,  36. ,  34.5,  33. ,  31.5,
            30. ,  28.5,  27. ,  25.5,  24. ,  22.5,  21. ,  19.5,  18. ,  16.5,
            15. ,  13.5,  12. ,  10.5,   9. ,   7.5,   6. ,   4.5,   3. ,   1.5,
             0. ,  -1.5,  -3. ,  -4.5,  -6. ,  -7.5,  -9. , -10.5, -12. , -13.5,
           -15. , -16.5, -18. , -19.5, -21. , -22.5, -24. , -25.5, -27. , -28.5,
           -30. , -31.5, -33. , -34.5, -36. , -37.5, -39. , -40.5, -42. , -43.5,
           -45. , -46.5, -48. , -49.5, -51. , -52.5, -54. , -55.5, -57. , -58.5,
           -60. , -61.5, -63. , -64.5, -66. , -67.5, -69. , -70.5, -72. , -73.5,
           -75. , -76.5, -78. , -79.5, -81. , -82.5, -84. , -85.5, -87. , -88.5,
           -90. ], dtype=float32)

  • level
    (level)
    int32
    1 2 3 5 7 ... 900 925 950 975 1000

    units :

    millibars

    long_name :

    pressure_level

    array([   1,    2,    3,    5,    7,   10,   20,   30,   50,   70,  100,  125,
            150,  175,  200,  225,  250,  300,  350,  400,  450,  500,  550,  600,
            650,  700,  750,  775,  800,  825,  850,  875,  900,  925,  950,  975,
           1000], dtype=int32)

  • time
    (time)
    datetime64[ns]
    2005-01-23 ... 2005-01-30T18:00:00

    long_name :

    time

    array(['2005-01-23T00:00:00.000000000', '2005-01-23T06:00:00.000000000',
           '2005-01-23T12:00:00.000000000', '2005-01-23T18:00:00.000000000',
           '2005-01-24T00:00:00.000000000', '2005-01-24T06:00:00.000000000',
           '2005-01-24T12:00:00.000000000', '2005-01-24T18:00:00.000000000',
           '2005-01-25T00:00:00.000000000', '2005-01-25T06:00:00.000000000',
           '2005-01-25T12:00:00.000000000', '2005-01-25T18:00:00.000000000',
           '2005-01-26T00:00:00.000000000', '2005-01-26T06:00:00.000000000',
           '2005-01-26T12:00:00.000000000', '2005-01-26T18:00:00.000000000',
           '2005-01-27T00:00:00.000000000', '2005-01-27T06:00:00.000000000',
           '2005-01-27T12:00:00.000000000', '2005-01-27T18:00:00.000000000',
           '2005-01-28T00:00:00.000000000', '2005-01-28T06:00:00.000000000',
           '2005-01-28T12:00:00.000000000', '2005-01-28T18:00:00.000000000',
           '2005-01-29T00:00:00.000000000', '2005-01-29T06:00:00.000000000',
           '2005-01-29T12:00:00.000000000', '2005-01-29T18:00:00.000000000',
           '2005-01-30T00:00:00.000000000', '2005-01-30T06:00:00.000000000',
           '2005-01-30T12:00:00.000000000', '2005-01-30T18:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (1)
  • u
    (time, level, latitude, longitude)
    float32
    ...

    units :

    m s**-1

    long_name :

    U component of wind

    standard_name :

    eastward_wind

    [34383360 values with dtype=float32]

• Indexes: (4)
  • longitude
    PandasIndex

    PandasIndex(Float64Index([  0.0,   1.5,   3.0,   4.5,   6.0,   7.5,   9.0,  10.5,  12.0,
                   13.5,
                  ...
                  345.0, 346.5, 348.0, 349.5, 351.0, 352.5, 354.0, 355.5, 357.0,
                  358.5],
                 dtype='float64', name='longitude', length=240))

  • latitude
    PandasIndex

    PandasIndex(Float64Index([ 90.0,  88.5,  87.0,  85.5,  84.0,  82.5,  81.0,  79.5,  78.0,
                   76.5,
                  ...
                  -76.5, -78.0, -79.5, -81.0, -82.5, -84.0, -85.5, -87.0, -88.5,
                  -90.0],
                 dtype='float64', name='latitude', length=121))

  • level
    PandasIndex

    PandasIndex(Int64Index([   1,    2,    3,    5,    7,   10,   20,   30,   50,   70,  100,
                 125,  150,  175,  200,  225,  250,  300,  350,  400,  450,  500,
                 550,  600,  650,  700,  750,  775,  800,  825,  850,  875,  900,
                 925,  950,  975, 1000],
               dtype='int64', name='level'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2005-01-23 00:00:00', '2005-01-23 06:00:00',
                   '2005-01-23 12:00:00', '2005-01-23 18:00:00',
                   '2005-01-24 00:00:00', '2005-01-24 06:00:00',
                   '2005-01-24 12:00:00', '2005-01-24 18:00:00',
                   '2005-01-25 00:00:00', '2005-01-25 06:00:00',
                   '2005-01-25 12:00:00', '2005-01-25 18:00:00',
                   '2005-01-26 00:00:00', '2005-01-26 06:00:00',
                   '2005-01-26 12:00:00', '2005-01-26 18:00:00',
                   '2005-01-27 00:00:00', '2005-01-27 06:00:00',
                   '2005-01-27 12:00:00', '2005-01-27 18:00:00',
                   '2005-01-28 00:00:00', '2005-01-28 06:00:00',
                   '2005-01-28 12:00:00', '2005-01-28 18:00:00',
                   '2005-01-29 00:00:00', '2005-01-29 06:00:00',
                   '2005-01-29 12:00:00', '2005-01-29 18:00:00',
                   '2005-01-30 00:00:00', '2005-01-30 06:00:00',
                   '2005-01-30 12:00:00', '2005-01-30 18:00:00'],
                  dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (2)

  Conventions :

  CF-1.6

  history :

  2018-07-17 16:50:39 GMT by grib_to_netcdf-2.8.0: grib_to_netcdf /data/data03/scratch/a7/d3/_mars-atls14-98f536083ae965b31b0d04811be6f4c6-WcwILE.grib -o /data/data02/scratch/e1/58/_grib2netcdf-atls13-a82bacafb5c306db76464bc7e824bb75-I1oqHm.nc -utime

Load the dimension arrays

In this version, the QGField object takes only: - latitude array in degree ascending order, and - pressure level in hPa in decending order (from ground to aloft).

xlon = u_file.longitude.values

# latitude has to be in ascending order
ylat = u_file.latitude.values
if np.diff(ylat)[0]<0:
    print('Flip ylat.')
    ylat = ylat[::-1]

# pressure level has to be in descending order (ascending height)
plev = u_file.level.values
if np.diff(plev)[0]>0:
    print('Flip plev.')
    plev = plev[::-1]

nlon = xlon.size
nlat = ylat.size
nlev = plev.size

Flip ylat.
Flip plev.

clat = np.cos(np.deg2rad(ylat))     # cosine latitude
p0 = 1000.                          # surface pressure [hPa]
kmax = 49                           # number of grid points for vertical extrapolation (dimension of height)
dz = 1000.                          # differential height element
height = np.arange(0,kmax)*dz       # pseudoheight [m]
dphi = np.diff(ylat)[0]*pi/180.     # differential latitudinal element
dlambda = np.diff(xlon)[0]*pi/180.  # differential latitudinal element
hh = 7000.                          # scale height
cp = 1004.                          # heat capacity of dry air
rr = 287.                           # gas constant
omega = 7.29e-5                     # rotation rate of the earth
aa = 6.378e+6                       # earth radius
prefactor = np.array([np.exp(-z/hh) for z in height[1:]]).sum() # integrated sum of density from the level
                                                                #just above the ground (z=1km) to aloft
npart = nlat                        # number of partitions to construct the equivalent latitude grids
maxits = 100000                     # maximum number of iteration in the SOR solver to solve for reference state
tol = 1.e-5                         # tolerance that define convergence of solution
rjac = 0.95                         # spectral radius of the Jacobi iteration in the SOR solver.
jd = nlat//2+1                      # (one plus) index of latitude grid point with value 0 deg
                                    # This is to be input to fortran code. The index convention is different.

Set the level of pressure and the timestamp to display below

tstamp = [dt.datetime(2005,1,23,0,0) + dt.timedelta(seconds=6*3600) * tt for tt in range(ntimes)]
plev_selected = 10 # selected pressure level to display
tstep_selected = 0

Loop through the input file and plot computed quantities

for tstep in range(32):  # or ntimes

    uu = u_file.u.isel(time=tstep).values[::-1, ::-1, :]
    vv = v_file.v.isel(time=tstep).values[::-1, ::-1, :]
    tt = t_file.t.isel(time=tstep).values[::-1, ::-1, :]

    qgfield_object = QGFieldNH18(xlon, ylat, plev, uu, vv, tt, northern_hemisphere_results_only=False)

    qgfield_object.interpolate_fields(return_named_tuple=False)

    qgfield_object.compute_reference_states(return_named_tuple=False)

    qgfield_object.compute_lwa_and_barotropic_fluxes(return_named_tuple=False)

    if tstep == tstep_selected:
        # === Below demonstrate another way to access the computed variables ===
        # 3D Variables that I would choose one pressure level to display
        variables_3d = [
            (qgfield_object.qgpv, 'Quasigeostrophic potential vorticity (QGPV)'),
            (qgfield_object.lwa, 'Local wave activity (LWA)'),
            (qgfield_object.interpolated_u, 'Interpolated zonal wind (u)'),
            (qgfield_object.interpolated_v, 'Interpolated meridional wind (v)')]

        # Reference states to be displayed on y-z plane
        variables_yz = [
            (qgfield_object.qref, 'Qref'),
            (qgfield_object.uref, 'Uref'),
            (qgfield_object.ptref, 'PTref')]

        # Vertically averaged variables to be displayed on x-y plane
        variables_xy = [
            (qgfield_object.adv_flux_f1, 'Advective flux F1'),
            (qgfield_object.adv_flux_f2, 'Advective flux F2'),
            (qgfield_object.adv_flux_f3, 'Advective flux F3'),
            (qgfield_object.convergence_zonal_advective_flux, 'Advective flux convergence -Div(F1+F2+F3)'),
            (qgfield_object.divergence_eddy_momentum_flux, 'divergence_eddy_momentum_flux'),
            (qgfield_object.meridional_heat_flux, 'meridional_heat_flux')
        ]

        # Plot 240 hPa of 3D-variables
        for variable, name in variables_3d:
            plt.figure(figsize=(12,6))
            plt.contourf(xlon, ylat[1:-1], variable[plev_selected, 1:-1, :], 50, cmap='jet')
            if name=='Local wave activity (LWA)':
                plt.axhline(y=0, c='w', lw=30)
            plt.colorbar()
            plt.ylabel('Latitude (deg)')
            plt.xlabel('Longitude (deg)')
            plt.title(name + ' at 240hPa | ' + str(tstamp[tstep]))
            plt.show()

        # Plot reference states
        for variable, name in variables_yz:
            plt.figure(figsize=(6,4))
            plt.contourf(ylat[1:-1], height, variable[:, 1:-1], 50, cmap='jet')
            plt.axvline(x=0, c='w', lw=2)
            plt.xlabel('Latitude (deg)')
            plt.ylabel('Pseudoheight (m)')
            plt.colorbar()
            plt.title(name + ' | ' + str(tstamp[tstep]))
            plt.show()

        # Plot barotropic (2D-)variables
        for variable, name in variables_xy:
            plt.figure(figsize=(12,6))
            plt.contourf(xlon, ylat[1:-1], variable[1:-1, :], 50, cmap='jet')
            plt.axhline(y=0, c='w', lw=30)
            plt.ylabel('Latitude (deg)')
            plt.xlabel('Longitude (deg)')
            plt.colorbar()
            plt.title(name + ' | ' + str(tstamp[tstep]))
            plt.show()

    print('tstep = {}/{}\n'.format(tstep, ntimes))

  1128425192       32684  1132480416  converged at n =          950
  1128425192       32684  1132480416  converged at n =          721

tstep = 0/32

  1128425192       32684  1132480416  converged at n =          953
  1128425192       32684  1132480416  converged at n =          721
tstep = 1/32

  1128425192       32684  1132480416  converged at n =          949
  1128425192       32684  1132480416  converged at n =          703
tstep = 2/32

   853002984       32684  1132480416  converged at n =          948
   853002984       32684  1132480416  converged at n =          711
tstep = 3/32

  1120484072       32684  1132480416  converged at n =          949
  1120484072       32684  1132480416  converged at n =          718
tstep = 4/32

  1120484072       32684  1132480416  converged at n =          948
  1120484072       32684  1132480416  converged at n =          721
tstep = 5/32

    67727080       32684  1132480416  converged at n =          946
    67727080       32684  1132480416  converged at n =          706
tstep = 6/32

    67727080       32684  1132480416  converged at n =          944
    67727080       32684  1132480416  converged at n =          718
tstep = 7/32

    67727080       32684  1132480416  converged at n =          947
    67727080       32684  1132480416  converged at n =          709
tstep = 8/32

  1128425192       32684  1132480416  converged at n =          946
  1128425192       32684  1132480416  converged at n =          716
tstep = 9/32

   853002984       32684  1132480416  converged at n =          942
   853002984       32684  1132480416  converged at n =          714
tstep = 10/32

    67754216       32684  1132480416  converged at n =          943
    67754216       32684  1132480416  converged at n =          719
tstep = 11/32

    67727080       32684  1132480416  converged at n =          941
    67727080       32684  1132480416  converged at n =          716
tstep = 12/32

    67754216       32684  1132480416  converged at n =          942
    67754216       32684  1132480416  converged at n =          725
tstep = 13/32

    67727080       32684  1132480416  converged at n =          945
    67727080       32684  1132480416  converged at n =          722
tstep = 14/32

  1128425192       32684  1132480416  converged at n =          943
  1128425192       32684  1132480416  converged at n =          734
tstep = 15/32

   588451048       32684  1132480416  converged at n =          946
   588451048       32684  1132480416  converged at n =          735
tstep = 16/32

    67754216       32684  1132480416  converged at n =          944
    67754216       32684  1132480416  converged at n =          742
tstep = 17/32

  1128425192       32684  1132480416  converged at n =          945
  1128425192       32684  1132480416  converged at n =          735
tstep = 18/32

   581618920       32684  1132480416  converged at n =          946
   581618920       32684  1132480416  converged at n =          749
tstep = 19/32

   581618920       32684  1132480416  converged at n =          943
   581618920       32684  1132480416  converged at n =          749
tstep = 20/32

    67727080       32684  1132480416  converged at n =          945
    67727080       32684  1132480416  converged at n =          750
tstep = 21/32

  1120484072       32684  1132480416  converged at n =          943
  1120484072       32684  1132480416  converged at n =          742
tstep = 22/32

   588451048       32684  1132480416  converged at n =          940
   588451048       32684  1132480416  converged at n =          752
tstep = 23/32

  1120484072       32684  1132480416  converged at n =          941
  1120484072       32684  1132480416  converged at n =          749
tstep = 24/32

  1128425192       32684  1132480416  converged at n =          940
  1128425192       32684  1132480416  converged at n =          747
tstep = 25/32

   588451048       32684  1132480416  converged at n =          943
   588451048       32684  1132480416  converged at n =          736
tstep = 26/32

  1128425192       32684  1132480416  converged at n =          944
  1128425192       32684  1132480416  converged at n =          738
tstep = 27/32

  1120484072       32684  1132480416  converged at n =          945
  1120484072       32684  1132480416  converged at n =          736
tstep = 28/32

  1128425192       32684  1132480416  converged at n =          945
  1128425192       32684  1132480416  converged at n =          733
tstep = 29/32

    67727080       32684  1132480416  converged at n =          943
    67727080       32684  1132480416  converged at n =          731
tstep = 30/32

   853002984       32684  1132480416  converged at n =          944
   853002984       32684  1132480416  converged at n =          736
tstep = 31/32