Negative O3 values in the zonal mean

[plesager]

We found negative values in the monthly zonal mean of O3 on pressure levels in our AerChem cmorized data:

Vertical axis is pressure and is upside down. Clearly at the surface, near the poles, something goes wrong. I guess that the high pressure is not always or never reached there and the interpolated value is incorrect. The calculation is done in tm52cmor.py, where the interpolation to pressure levels is done first https://github.com/EC-Earth/ece2cmor3/blob/54e36a782c6d9471d39a80677234b0113155fbf9/ece2cmor3/tm52cmor.py#L432 and then a simple numpy,mean is applied https://github.com/EC-Earth/ece2cmor3/blob/54e36a782c6d9471d39a80677234b0113155fbf9/ece2cmor3/tm52cmor.py#L443

I think that the pressure interpolation is not working as expected. Shouldn't we set to NaN, if the pressure target is out of the range of the pressure input? @tommibergman @noije

[noije]

Yes, when calculating the zonal mean at a certain pressure level, mixing ratios should be set to NaN at elevated locations where the surface pressure is lower than that pressure.

[tommibergman]

I agree, it is not correct. Unfortunately, I cannot fix this now. Just happens that we found this bug at a very unfortunate time.

[plesager]

This can always be recreated from the 3D model level cmorized output. If I have time, I will write a recipe to replace the data in place. That's the easiest solution, or we just do not provide this zonal mean.

[plesager]

Looking at the code, I wonder is a quick fix is possible. The interpolation is now done like this:

https://github.com/EC-Earth/ece2cmor3/blob/54e36a782c6d9471d39a80677234b0113155fbf9/ece2cmor3/tm52cmor.py#L559

but the documentation for the last flag says:

If False, then no extrapolation is done when the pressure level is outside of the range of psfc (psdata in our call).

To test.

[plesager]

I did a quick test setting the last argument to False in the call to Ngl.vinth2p. This seems to be a step in the right direction, although the NaN are not correctly hanlded (I use the same color scale as the in the first plot above):

The top area in purple has the 1e+30 value everywhere.

[plesager]

Further tests with the corrected SP in EC-Earth output. Here is a reference. This a zonal slice (not mean) of O3 on model levels from the AERmon table. Y-axis goes from 0 (bottom) to 1 (top):

There is no spatial average, only a monthly mean. Data range in the entire atmosphere is [8e-9, 9.7e-6]

[plesager]

Now the zonal mean on 39 pressure levels. This involves vertical interpolation, followed by the zonal mean. I am using the correct surface pressure in the vertical interpolation (as fixed in #603). The result is still wrong. The range of data is now [1.6e-8, 1.8e-4] with the largest values at the surface (top plot), but the top of the atmosphere is ok (bottom plot where I use the same color scale as the reference plot above):

[plesager]

So the only improvement so far is that the negative values disappeared.

[plesager]

For the record, I tried to remove the extrapolation out of the original range during the vertical interpolation (setting the last argument to False in the call to Ngl.vinth2p), and accounting for NaN (white area):

It is not clear to me what False in the call to Ngl.vinth2p exactly does, but this seems too aggressive.

[plesager]

I think there is still something wrong with the vertical interpolation on plevel. Needs to be fixed or clarified before we can address the zonal mean.

[noije]

Further tests with the corrected SP in EC-Earth output. Here is a reference. This a zonal slice (not mean) of O3 on model levels from the AERmon table. Y-axis goes from 0 (bottom) to 1 (top):

There is no spatial average, only a monthly mean. Data range in the entire atmosphere is [8e-9, 9.7e-6]

The range is ok, but I wonder about the vertical axis. Is the surface (y=0) really at the bottom of the graph? Does the vertical distribution agree with that shown in https://dev.ec-earth.org/issues/752#note-41?

[noije]

So, in the plot the y-axis goes from 0 at the bottom (=TOA) to 1 at the top (=surface), and the problem appears only near the surface/in the lower troposphere? Is that correct. I have the feeling this is indeed related to the treatment of the points with high surface elevation, where the surface pressure is lower than the required pressure level. These should be set to NaN in the vertical interpolation, and should be excluded when calculating zonal means at those pressure levels. @plesager I understand you already tried that, right?

[plesager]

Re: vertical axis. That's correct. Re: setting NaN. I tried something different, namely avoiding extrapolation out of the range of data of input data. I will try your proposition.

[plesager]

I rerun everything. Things are looking better. Model levels on the left, zonal mean on 39 pressure levels on the right:

Note that there is still some small negative in the interpolated data, as indicated by the data range. They are found near the surface above the Antartic. But this is a lot better than the first plot.

[noije]

I suspect the small negatives appear where all points at that latitude have surface pressure lower than the pressure level of interest. If that is indeed the case, the zonal mean is undefined and these values should in fact be set to NaN.

[plesager]

I think it work the other way around: if the target pressure is smaller than what's available in the column, then you have to extrapolate. And we are in the opposite case. Or do I miss something?

[noije]

Yes, but extrapolation can produce negatives in case of strong vertical gradients. Just assuming the mixing ratios in the lowest model level are valid throughout the whole cell, i.e. down to the surface, and assuming NaN at higher pressures, is the safest thing to do.

[tommibergman]

I think this also resolved by PR #618

[treerink]

Solved by merging #618.