Option to output spectrally-decomposed TOA flux from RRTMG_LW

[AndrewILWilliams]

This PR will close #137 (when finished).

The interface is the same as before, but if you want the spectrally-decomposed OLR outputted, then you pass return_spectral_olr=True to RRTMG_LW when initialized.

I've implemented a test which asserts that assert np.allclose(simps(rcm.OLR_sr, rcm.RRTMG_LW_bands), rcm.OLR), after a 5-year RCE simulation.

I'm a bit confused by why this test is failing at the moment, it says that AttributeError: 'TimeDependentProcess' object has no attribute 'RRTMG_LW_bands', but as far as I can tell, I'm setting it as an attribute on line 56 of climlab/radiation/rrtm/rrtmg_lw.py ?

I think that the FORTRAN is coded correctly, my only confusion is whether I should include the delwave factor in line 546 of climlab/radiation/rrtm/_rrtmg_lw/sourcemods/rrtmg_lw_rtrnmc.f90. I'm not that familiar with how RRTM works so I'll have to defer to other on this.

Thanks for the help so far, @brian-rose !
A

[brian-rose]

Great start! I'll have some detailed feedback for you a little later.

[brian-rose]

Sorry, the AMS meeting has sucked up all my attention this week but I haven't forgot about this!

[AndrewILWilliams]

No problem! Hope it's going well! :)

[brian-rose]

Some (maybe all) the test errors probably came from a failure of the data server that gets called upon first import of climlab.

Try merging the latest master (which will fix an annoying problem with pytest warnings) and let the tests run again.

[AndrewILWilliams]

Hi @brian-rose , I think I've merged the latest master. The CI tests are still failing though...

[AndrewILWilliams]

When I run the tests locally I'm getting a lot of warnings due to ozone stuff:

  /opt/anaconda3/envs/test_env/lib/python3.7/site-packages/climlab/radiation/radiation.py:164: UserWarning: Interpolation of ozone data failed. Setting O3 to zero instead.
    warnings.warn('Interpolation of ozone data failed. Setting O3 to zero instead.')

-- Docs: https://docs.pytest.org/en/stable/warnings.html
========================================================== short test summary info ==========================================================
FAILED climlab/tests/test_rrtm.py::test_spectral_olr - AttributeError: 'TimeDependentProcess' object has no attribute 'RRTMG_LW_bands'
=========================================== 1 failed, 53 passed, 56 warnings in 66.99s (0:01:06) ============================================
 o3neg!           1   0.0000000000000000     
 o3neg!           2   0.0000000000000000     
 o3neg!           3   0.0000000000000000     
 o3neg!           4   0.0000000000000000     
 o3neg!           5   0.0000000000000000     
 o3neg!           6   0.0000000000000000     
 o3neg!           7   0.0000000000000000     
 o3neg!           8   0.0000000000000000     
 ...
 o3neg!          47   0.0000000000000000     
 o3neg!          48   0.0000000000000000     
 o3neg!          49   0.0000000000000000     
 o3neg!          50   0.0000000000000000  

This seems to be due to nans in the ozone array after interpolating to the xTatm grid, even though fill_value='None' should mean that it extrapolates rather than throwing a nan.

Update:
These O3 warnings go away if in climlab.radiation.radiation.default_absorbers, on line 161, you change

O3 = O3source.interp_like(xTatm, kwargs={'fill_value':None})

to

O3 = O3source.interp_like(xTatm, kwargs={'fill_value':'extrapolate'})

I think this is because you're only interpolating over 1d coordinates, and so scipy.interpolate.interp1d() is called, which requires 'fill_value':'extrapolate' in order to extrapolate rather than throwing a nan. If you're interpolating over multidimensional coordinates, then scipy.interpolate.interpnd() is called, which requires 'fill_value':None.

[brian-rose]

@AndrewWilliams3142 I don't think those O3 issues would be causing the failures in test_spectral_olr. It's still not clear to me what's going on. But good to look at things one at a time.

It would be helpful if you put in a separate PR to fix the O3 issue, and we'll make sure that passes all test.

[brian-rose]

For the record, #140 fixed the O3 issue, and has now been merged into this PR.

[codecov-commenter]

Codecov Report

❗ No coverage uploaded for pull request base (main@3c441f7). Click here to learn what that means.
The diff coverage is n/a.

❗ Current head 06af6fe differs from pull request most recent head 81e2b5a. Consider uploading reports for the commit 81e2b5a to get more accurate results

@@          Coverage Diff           @@
##             main    #138   +/-   ##
======================================
  Coverage        ?   0.00%           
======================================
  Files           ?      66           
  Lines           ?    3562           
  Branches        ?       0           
======================================
  Hits            ?       0           
  Misses          ?    3562           
  Partials        ?       0           

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 3c441f7...81e2b5a. Read the comment docs.

[AndrewILWilliams]

Hey @brian-rose , no problem at all. I think I've managed to do the changes you asked for? Tests are being triggered, we'll see if they pass.. :)

[brian-rose]

Excellent, I'll have a careful look and see what needs doing in terms of documentation.

[brian-rose]

Have you tested this with different grids? I'm guessing that olr_sr will need to be reshaped.

[AndrewILWilliams]

Oooh good point, I've only tested for single-column models so far, do you mean to try with a 2d model?

[brian-rose]

Yes exactly. Climlab tries to support arbitrary sized lat-lon grids. The helper function _rrtm_to_climlab() does the necessary reshaping. It can be a huge pain. In this case since there's an extra dimension in the output field it probably won't work out of the box.

[AndrewILWilliams]

Thanks @brian-rose ! I'm a bit confused by this. It seems like _rrtm_to_climlab() just reverses the ordering of the pressure layers?

Also, maybe a dumb question but I didn't know climlab had any models with a longitude component?

[brian-rose]

Also, maybe a dumb question but I didn't know climlab had any models with a longitude component?

Not dumb at all. Support for lat-lon grids is more of an aspirational goal for climlab, and currently spotty at best, e.g. https://github.com/brian-rose/climlab/issues/86. I wasn't thinking clearly in my previous comment. For testing here, we should focus on 2D latitude-pressure grids

[AndrewILWilliams]

Cool, thanks for the clarification! It works well with test_large_grid(), which is a latitude-pressure grid, so is there more to do here? I don't think there's a need to call _rrtmg_to_climlab() on olr_sr because we don't want to flip the ordering of the last dimension

[brian-rose]

Taking a closer look at your code. This is pretty awesome! I just made one small change, removing the ispec argument from the input arguments for climlab.radiation.RRTMG(). It was redundant with the flag return_spectral_olr so best not to expose it to users.

[brian-rose]

A couple of things I'm wondering about:

• You hard-code central_rrtmg_lw_bands as literal constants in rrmtg_lw.py. Is it possible instead to import these values from the Fortran extension, mirroring what we do for nbndlw and ngptlw? Best to avoid doubly defining constants wherever possible.
• I played around with the tests and found that the spectral OLR sum comes within about 0.02 W/m2 of the total OLR. This is certainly close and probably close enough for most applications, but I guess I would have expected to get closer to machine-accuracy level. Are there some details of the spectral integration that we could implement diagnostically to close that budget better?

[AndrewILWilliams]

Thanks @brian-rose ! Good catch on the ispec argument as well. Regarding your other points:

• I think it would be possible to import them from the Fortran. Currently the upper/lower bounds of each band are here, and so it would be possible to import them both and then calculate central_rrtmg_lw_bands from that. I think this also links into a previous question I had actually! Is it possible to make these central_rrtmg_lw_bands an attribute of the output somehow? In xarray I'd think of this as the coordinate of the spectral dimension, but I'm not 100% sure what the analogue would be with climlab's Field object.

• Another good point! I spoke with Mike Iacono at AER before about these changes to RRTMG and this was how he suggested to implement the spectral OLR diagnostic, I can have another look though!

[brian-rose]

• I think it would be possible to import them from the Fortran. Currently the upper/lower bounds of each band are here, and so it would be possible to import them both and then calculate central_rrtmg_lw_bands from that. I think this also links into a previous question I had actually! Is it possible to make these central_rrtmg_lw_bands an attribute of the output somehow? In xarray I'd think of this as the coordinate of the spectral dimension, but I'm not 100% sure what the analogue would be with climlab's Field object.

It would be possible to implement both the upper/lower bounds and centers of each interval as points and bounds properties of the domain object. I originally set up the climlab Field class to understand these kind of finite volume grid concepts.

I'll see if I can mock something up quickly.

On the other hand, I don't think it's worth putting too much effort into these details because we are planning a complete rewrite of the climlab internals that will get rid of the Field and domain objects entirely and use xarray structures instead (longstanding issue #50)

[brian-rose]

Ok it wasn't that quick, because the climlab Field and domain objects are a ridiculous pain to work with, but the latest commit wraps the output in a labeled Field object that contains the spectral band centers.

Easiest way for the user to access this would be using the to_xarray() interface like:

import climlab
state = climlab.column_state()
rad = climlab.radiation.RRTMG(state=state, return_spectral_olr=True)
rad.step_forward()
rad.OLR_sr.to_xarray()

which produces

<xarray.DataArray (depth: 1, wavenumber: 16)>
array([[3.27982316e+01, 3.97397780e+01, 3.60661874e+01, 7.08628464e+00,
        3.10797673e+01, 4.75361385e+01, 1.55245526e+01, 1.76265583e+01,
        1.56793372e+01, 1.76402992e+00, 2.02886234e+00, 1.99654135e+00,
        7.24616719e-01, 3.61792952e-03, 3.41639617e-01, 2.61125909e-01]])
Coordinates:
  * depth       (depth) float64 0.5
  * wavenumber  (wavenumber) float64 180.0 425.0 565.0 ... 2.49e+03 2.925e+03

Ideally we would also wrap the spectral band boundaries into this object. You can access this information through the native climlab interface like

In [7]: rad.OLR_sr.domain.axes['wavenumber'].bounds
Out[7]:
array([  10.,  350.,  500.,  630.,  700.,  820.,  980., 1080., 1180.,
       1390., 1480., 1800., 2080., 2250., 2380., 2600., 3250.])

but the .to_xarray() wrapper doesn't know how to deal with this.

[brian-rose]

Just merged the latest updates from main branch, shouldn't cause any new failures.

[brian-rose]

@AndrewWilliams3142 I'd like to advocate for changing OLR_sr to OLR_spectral so the user-facing diagnostic is less mysterious. Does that sound ok to you?

[AndrewILWilliams]

Thanks so much for doing this @brian-rose ! Yes I couldn't get my head around how to implement this with the Field object. This implementation with the to_xarray() is great though!

[AndrewILWilliams]

@AndrewWilliams3142 I'd like to advocate for changing OLR_sr to OLR_spectral so the user-facing diagnostic is less mysterious. Does that sound ok to you?

Fine by me! :)

[brian-rose]

@AndrewWilliams3142 any interest in putting together a notebook with a demonstration of how to use this new feature? It would make a good addition to the tutorial collection at https://climlab.readthedocs.io/en/latest/tutorial.html

[AndrewILWilliams]

Yep, I'd be down for that! Would it be easier to do this post-merge as a separate PR? I've done some stuff with the spectral OLR already in this repo (reproducing some of the Seeley+Jeevanjee 2020 paper on ECS) and so it would be easy to move that over :)

[brian-rose]

Great stuff!
Yes, I think we should merge this now and work on a tutorial notebook in a separate PR. But it would be great to have a working example to show off for the next climlab release.

I'll open a new issue for that and we can move the conversation over there.

Thanks so much for your contribution here, really great new feature!