ValueError: The reference state does not converge for Northern Hemisphere. about hn2016_falwa HOT 17 CLOSED

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?

from hn2016_falwa.

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!

from hn2016_falwa.

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)

from hn2016_falwa.

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.

from hn2016_falwa.

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!

from hn2016_falwa.

@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?

from hn2016_falwa.

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.

from hn2016_falwa.

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.

from hn2016_falwa.

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?

from hn2016_falwa.

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

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

-Clare

from hn2016_falwa.