ValueError: The reference state does not converge for Northern Hemisphere.

[hs2112]

Hello again! I was doing a trial run with some data I had and ran into the given error. The interpolated fields look good to me but qgfield_object.compute_reference_states(northern_hemisphere_results_only=False) fails. I will take a closer look at the values tomorrow but please let me know if I should look for something specific for debugging this. Thanks in advance.

[csyhuang]

Hi @hs2112 , this behavior is observed sometimes with fine resolution data and where the QGPV field around equator changes sign. May you possibly share a snapshot of u, v and temperature field at the instance it fails (maybe in netCDF file - you may upload it on https://gofile.io/ and send me the link to the files here) such that I can check for you? Please also share the python code you use here. Thanks!

[hs2112]

Hi @csyhuang. Is it not possible to compute LWA for a subset of NH instead of the entire globe?

Here is the untreated data (u, v, t combined) on the filtered timestep - https://gofile.io/d/Wqurnq. But I tried on some other timesteps and it throws the same error every time. Note that my pressure levels are spaced by as much as 100hPa in mid troposphere and kmax=17 is done based on the vertical extent available.

I am running the following lines in python -

tstep=0
uu = u_file.variables['u'][tstep, ::-1, ::-1, :].data
vv = v_file.variables['v'][tstep, ::-1, ::-1, :].data
tt = t_file.variables['t'][tstep, ::-1, ::-1, :].data
qgfield_object = QGField(xlon, ylat, plev, uu, vv, tt, kmax=17)
qgpv[tstep, :, :, :], interpolated_u[tstep, :, :, :], interpolated_v[tstep, :, :, :], \
    interpolated_theta[tstep, :, :, :], static_stability = qgfield_object.interpolate_fields()
qref[tstep, :, :], uref[tstep, :, :], ptref[tstep, :, :] = \
    qgfield_object.compute_reference_states(northern_hemisphere_results_only=False)

[chpolste]

Hello @hs2112,

I had a quick look at the dataset and can reproduce your issue:

>>> import xarray as xr
>>> from hn2016_falwa.xarrayinterface import QGDataset
>>> data = xr.open_dataset("export_tstep_0.nc")
>>> qgd = QGDataset(data, qgfield_kwargs={
...     "kmax": 17
... })
>>> qgd.interpolate_fields()
<xarray.Dataset>
>>> qgd.compute_reference_states()
 Maxits exceeded
           0       32764  1196400352  converged at n =       100001
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File ".../python3.9/site-packages/hn2016_falwa/xarrayinterface.py", line 245, in compute_reference_states
    out_fields = _map_collect(
  File ".../python3.9/site-packages/hn2016_falwa/xarrayinterface.py", line 45, in _map_collect
    for name, y in zip(names, f(x)):
  File ".../python3.9/site-packages/hn2016_falwa/xarrayinterface.py", line 246, in <lambda>
    lambda field: field.compute_reference_states(northern_hemisphere_results_only),
  File ".../python3.9/site-packages/hn2016_falwa/oopinterface.py", line 534, in compute_reference_states
    raise ValueError("The reference state does not converge for Northern Hemisphere.")
ValueError: The reference state does not converge for Northern Hemisphere.

I noticed that the data has 0.25° horizontal spacing, which is a very high resolution and might be an issue for the numerical solver. If I coarsen the data to 1°, the reference state solver converges:

>>> coarsened = data.isel({
...     "latitude":  slice(None, None, 4),
...     "longitude": slice(None, None, 4)
... })
>>> qgd = QGDataset(coarsened, qgfield_kwargs={
...     "kmax": 17
... })
>>> qgd.interpolate_fields()
<xarray.Dataset>
>>> qgd.compute_reference_states()
  1196400384       32764  1196400088  converged at n =          716
           0           0           0  converged at n =          721
<xarray.Dataset>
Dimensions:  (height: 17, ylat: 181)
Coordinates:
  * height   (height) float64 0.0 1e+03 2e+03 3e+03 ... 1.4e+04 1.5e+04 1.6e+04
  * ylat     (ylat) float64 -90.0 -89.0 -88.0 -87.0 ... 87.0 88.0 89.0 90.0
Data variables:
    qref     (height, ylat) float64 -0.0 -0.0 -0.0 -0.0 -0.0 ... 0.0 0.0 0.0 0.0
    uref     (height, ylat) float32 0.0 0.0 0.0 0.0 ... 11.65 8.511 4.507 0.503
    ptref    (height, ylat) float32 208.3 209.4 211.6 ... 443.1 443.4 443.5
Attributes: (12/13)
    ...

I'm only using every 4th gridpoint here, which is generally not a good way to coarsen data. It's just for testing.

Is moving to a coarser horizontal resolution, e.g. 1°, an option for you?

[hs2112]

I see! The coarser resolution might not be a problem given that I'm trying to find high amplitude Rossby wave propagation/ breaking events. I will implement this tomorrow and let you know how it goes. Thank you so much for the timely responses.

[csyhuang]

Thanks @chpolste for the quick solution! Someone emailed me about the same issue with 0.25 deg resolution dataset before, and ended up with the solution of coarse sampling (1 deg) too (it worked). @hs2112 please keep us posted if that works on your side. Thanks!

[hs2112]

@chpolste @csyhuang Confirming that the reference state converged for 1 degree resolution. Could you also confirm if it is possible to compute LWA for a subset of NH, say beyond subtropics, instead of the entire globe?

[csyhuang]

Hi @hs2112 , yes it is possible to compute LWA for a subset of NH. In Neal et al 2022, we compute LWA using a reference state computed using data points from 5 deg N to the pole. You can refer to scripts/nhn_grl2022/sample_run_script.py for the usage. Here I use the data you provided to compute Qref and Uref (At the moment, this is not available in the xarray interface yet):

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from hn2016_falwa.oopinterface import QGField
data = xr.open_dataset("export_tstep_0.nc")
coarsened = data.isel({
    "latitude":  slice(None, None, 4),
    "longitude": slice(None, None, 4)})

uu = coarsened.variables['u'][:, ::-1, :]
vv = coarsened.variables['v'][:, ::-1, :]
tt = coarsened.variables['t'][:, ::-1, :]
xlon = np.arange(360)
ylat = np.arange(-90, 91, 1)
plev = coarsened.variables['level']
kmax = 17
dz = 1000.
height = np.arange(0, kmax)*dz

eq_boundary_index = 5  # use domain from 5N to North pole
qgfield_object = QGField(xlon, ylat, plev, uu, vv, tt, kmax=kmax, dz=dz, eq_boundary_index=eq_boundary_index)

qgpv_temp, interpolated_u_temp, interpolated_v_temp, interpolated_avort_temp, interpolated_theta_temp, \
static_stability_n, static_stability_s, tn0, ts0 = qgfield_object._interpolate_field_dirinv()

qref, uref, tref, fawa, ubar, tbar = qgfield_object._compute_qref_fawa_and_bc()

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
cs1 = ax1.contourf(ylat[90:], height, np.swapaxes(qref, 0, 1), 30, cmap='rainbow')
ax1.set_title('Qref')
cbar1 = fig.colorbar(cs1)
cs2 = ax2.contourf(ylat[90+eq_boundary_index:], height, np.swapaxes(uref, 0, 1), np.arange(-20, 30, 5), cmap='rainbow')
ax2.set_title('Uref')
cbar2 = fig.colorbar(cs2)
plt.show()

The resultant Qref and Uref (only available from 5N to the North Pole) looks like:

(x-label: latitude. y-label: pseudoheight) Is that something you want to do?

Update: Please find the most updated release on the master branch (v0.6.4) which I used to produce the figure above.

[csyhuang]

Please let me know if you have any more questions, or else I will close this ticket in a week (by 2023/4/4). Thanks.

[hs2112]

Great @csyhuang! Sorry I go caught up with some things; was not able to reply. Are the reference states are a function of the horizontal extent of data provided?

Had another quick query - Currently I'm inputting the entire vertical extent upto a certain level, so that the library can interpolate to peudoheight levels. But I'm interested in vertical levels of 10km to 12km and the pressure levels corresponding to these vertical levels are between 200 and 300hPa. So can I just upload these levels instead of the entire vertical extent upto 200hPa?

[csyhuang]

Hi @hs2112 ,

The current version of package only supports LWA budget analysis over the whole vertical column(so the entire vertical extent of data is necessary, since the vertical integral function is written for this configuration).

What analysis do you want to with levels between 200 and 300hPa? If you just want to LWA for those levels, you can use those levels of QGPV to compute the wave activity via the eqvlat_fawa function to get reference state for a pressure level surface and then the lwa function to compute LWA (a bugfix version is available on the release0.7.0 branch):

https://github.com/csyhuang/hn2016_falwa/blob/c6df73fd142f1a18beccff8137410046857ab9c5/hn2016_falwa/basis.py#L11

Is this something you're interested in computing? In any case, let me know.

-Clare

[hs2112]

Hi @csyhuang ,

So, compute_lwa_and_barotropic_fluxes outputs LWA at various pressure levels; does it give vertically integrated values up to each level?

I am looking to compute LWA corresponding to any given pressure level(s) only. I will eventually look at the vertical structure of wave activity itself.

Can I use QGPV for inputting to eqvlat_lwa, and is there a function to calculate QGPV using u, v, z, t fields with minimum hard-coding involved? Can't I use qgpv_eqlat_lwa directly, as given in Example_qgpv.ipynb?

[csyhuang]

@hs2112 Sure. You can use qgpv_eqlat_lwa as shown in Example_qgpv.ipynb to compute LWA. qgpv_eqlat_lwa is just a wrapper calling eqlat_lwa. You can look at the code base to see what it does: https://github.com/csyhuang/hn2016_falwa/blob/9acd8106cf8a359979d4fb6c763ced1cd91d666b/hn2016_falwa/wrapper.py#L197

Can you state clearly what you want to calculate? If you only want to look at LWA without looking at the fluxes (which requires computation of Uref, and that requires data for all pressure level), this is what you can do:

1. Put in the data among pressure levels just below and above 200-300hPa into QGField to make sure static stability is properly calculated for layers between 200-300hPa.

2. Still do the 3 steps:

   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)

   and access the QGPV(z, lat, lon) and LWA(z, lat, lon) calculated: qgfield.qgpv and qgfield.lwa. The values are on regular psudoheight grid. You can get the values for your desired pressure levels my interpolation back to your original grid.

(Note that if you do what I suggested above, Uref, Tref and wave activity flux quantities will not be correct and shall not be analyzed.)

Would that work? Let me know if you have questions.

[hs2112]

@csyhuang "If you only want to look at LWA without looking at the fluxes (which requires computation of Uref, and that requires data for all pressure level)," - How do you calculate LWA without Uref? Since it is a function of (y,z) I thought it should be readily calculable for each pressure level separately (NH16).

Sorry if I was not clear before. This is exactly what I want to do - I want to compute an array of hemisphere-wide LWA corresponding to an array of pressure levels between 900hPa and 200hPa, instead of vertically averaged. I am mainly interested in the LWA field at higher levels for establishing stationary wave activity during a certain period, but I would like to generate it at all possible levels above surface. Hence, I was asking about the LWA field that is an output of compute_lwa_and_barotropic_fluxes, which is available at pressure levels.

So, basically, I just want to use the QGPV-based LWA formulation per pressure level, as established in NH16. Of course, I would like to analyze the contributing fluxes too, if possible.

[hs2112]

I do currently have daily average 1x1 outputs from compute_lwa_and_barotropic_fluxes, and have been looking at the lwa field between 5N and 40N over and around South Asia (largely between Africa and west Pacific), where the signals are low amplitude and the field is rather noisy. Would the 276 hPa field be accurately represented, or largely bear the low level signature due to the vertical averaging? It would be helpful if you could comment on this.

My concern is that the vertical structure of the wave field might not be barotropic, especially due to the Himalayan topography and secondary gradients, so wanted to implement LWA ideas from NH16, to each pressure level separately. Does that make sense?

[csyhuang]

@hs2112 : Thanks for explaining your use case. Now I get your concern.

First of all, the functionalities available in our python package were only the ones we have used so far in our publications. My PhD thesis focuses on blocking over the ocean (blocking is the main subject of HN16, NH18 and NHN22 as well), and the structure of blocking is largely barotropic (an example is shown in HN16 Fig. 10), so it suffices to look at the barotropic component of fluxes for those studies. It does not mean that other studies of different focus shall only look at barotropic components.

Your case involves grid points over the continents, so the lower level fields are likely noisy or unreliable. (Are you using reanalysis data? It also depends on the treatment of topography in data assimilation process.) Computation of LWA per pressure level only requires Qref but not Uref (so what I suggested above can do the job):

Since layerwise LWA (not vertical average) is also computed within the method compute_lwa_and_barotropic_fluxes, you still have to run it but you can simply take qgfield.qgpv and qgfield.lwa from qgfield object, which are 3D variables in (z,y,x), and discard the rest that you don't need. you can then analyze QGPV and LWA over only the upper levels.

Would that work?

[hs2112]

@csyhuang Yes, I am using ERA5 reanalysis that uses a seemingly advanced Land Data Assimilation system, but I agree and am happy looking at upper levels.

I am already using qgpv and lwa generated by calling interpolate_fields and compute_lwa_and_barotropic_fluxes on the qgfield object, and following the 3 steps you suggested as pasted below. So, guess there's no need to make any changes (not using qgpv_eqlat_lwa separately) for layerwise LWA.

1. Did you say that the fluxes and some other fields may not be reliable to analyze layerwise?
2. Would you recommend trying out any other changes like using actual geopotential heights instead of the pseudoheight hardcoding?

qgfield_object = QGField(xlon, ylat, plev, uu, vv, tt, kmax=kmax)

qgpv[tstep, :, :, :], interpolated_u[tstep, :, :, :], interpolated_v[tstep, :, :, :], \
    interpolated_theta[tstep, :, :, :], static_stability = qgfield_object.interpolate_fields()

qref[tstep, :, :], uref[tstep, :, :], ptref[tstep, :, :] = \
    qgfield_object.compute_reference_states(northern_hemisphere_results_only=False)

adv_flux_f1[tstep, :, :], \
adv_flux_f2[tstep, :, :], \
adv_flux_f3[tstep, :, :], \
adv_flux_conv[tstep, :, :], \
divergence_eddy_momentum_flux[tstep, :, :], \
meridional_heat_flux[tstep, :, :], \
lwa_baro[tstep, :, :], \
u_baro[tstep, :, :], \
lwa[tstep, :, :, :] \
    = qgfield_object.compute_lwa_and_barotropic_fluxes(northern_hemisphere_results_only=False)

Thanks!

[hs2112]

I think my main problem is that the LWA field looks noisy. When I dig to see why, it roughly reflects the noisy QGPV field computed in compute_lwa_and_barotropic_fluxes, and the meridional extent also looks larger than indicated in QGPV. Pasting the QGPV (*10^4) and LWA plots (ignore the contours) along with ERA5 PV which is much more coherent.