Repository: brian-rose/climlab

Issue:
Interpolation of ozone fails with surface pressures higher than 1056 hPa
Trying to create a radiation model for a massive atmosphere, following method laid out in #39, creates a new error at model initialization time:

First here is the hack to create a domain with altered surface pressure:

```
import numpy as np
import climlab
def make_column(num_lev=20, ps=const.ps, water_depth=5.):
    state = climlab.column_state(num_lev=num_lev, num_lat=1, water_depth=water_depth)
    lev = state.Tatm.domain.lev
    lev.bounds = np.linspace(0., ps, num_lev+1)
    lev.points = lev.bounds[:-1] + np.diff(lev.bounds)/2.
    lev.delta = np.abs(np.diff(lev.bounds))
    return state
```

Now I try to use this to create a radiation model with `climlab.radiation.RRTMG`:

```
state = make_column(num_lev=10, ps=1056.)
r = climlab.radiation.RRTMG(state=state)
```

This works fine. But any larger surface pressure throws an error:

```
state = make_column(num_lev=10, ps=1057.)
r = climlab.radiation.RRTMG(state=state)
```
produces
```
ValueError: A value in x_new is above the interpolation range.
```

This occurs because of code in `climlab.radiation.radiation.py` that attempts to initialize the ozone distribution by interpolation from a climatology file with maximum surface pressure of 1056 hPa.

It's possible to sidestep this problem by passing the `absorber_vmr` dictionary as input to `climlab.radiation.RRTMG`. But we need a more graceful handling of this case. Perhaps just assuming that ozone is zero at higher pressures, and alerting the user about this.

Discussion:

Patch:
diff --git a/climlab/__init__.py b/climlab/__init__.py
index ab3f7fe..e531af9 100644
--- a/climlab/__init__.py
+++ b/climlab/__init__.py
@@ -7,7 +7,7 @@ Nevertheless also the underlying code of the ``climlab`` architecture
 has been documented for a comprehensive understanding and traceability.
 '''
 
-__version__ = '0.7.1.dev0'
+__version__ = '0.7.1.dev1'
 
 # this should ensure that we can still import constants.py as climlab.constants
 from .utils import constants
diff --git a/climlab/radiation/radiation.py b/climlab/radiation/radiation.py
index 7b6dd33..f9a7134 100644
--- a/climlab/radiation/radiation.py
+++ b/climlab/radiation/radiation.py
@@ -16,13 +16,17 @@ State:
 Input arguments (both LW and SW):
 
     - ``specific_humidity`` (kg/kg)
-    - ``absorber_vmr`` (dict of volumetric mixing ratios)
+    - ``absorber_vmr = None`` (dictionary of volumetric mixing ratios. Default values supplied if ``None``)
     - ``cldfrac``      (layer cloud fraction
     - ``clwp``         (in-cloud liquid water path (g/m2))
     - ``ciwp = 0.``,     # in-cloud ice water path (g/m2)
     - ``r_liq = 0.``,    # Cloud water drop effective radius (microns)
     - ``r_ice = 0.``,    # Cloud ice particle effective size (microns)
-    - ``ozone_file = 'apeozone_cam3_5_54.nc'``  (file with ozone distribution)
+    - ``ozone_file = 'apeozone_cam3_5_54.nc'``  (file with ozone distribution --
+    ignored if ``absorber_vmr`` is given)
+
+If ``absorber_vmr = None`` then ozone will be interpolated to the model grid
+from a climatology file, or set to zero if ``ozone_file = None``.
 
 Additional input arguments for SW:
 
@@ -128,34 +132,38 @@ def default_absorbers(Tatm,
     absorber_vmr['CFC22'] = 0.
     absorber_vmr['CCL4']  = 0.
 
-    datadir = os.path.join(os.path.dirname(__file__), 'data', 'ozone')
-    ozonefilepath = os.path.join(datadir, ozone_file)
-    #  Open the ozone data file
-    if verbose:
-        print('Getting ozone data from', ozonefilepath)
-    ozonedata = nc.Dataset(ozonefilepath)
-    ozone_lev = ozonedata.variables['lev'][:]
-    ozone_lat = ozonedata.variables['lat'][:]
-    #  zonal and time average
-    ozone_zon = np.mean(ozonedata.variables['OZONE'], axis=(0,3))
-    ozone_global = np.average(ozone_zon, weights=np.cos(np.deg2rad(ozone_lat)), axis=1)
-    lev = Tatm.domain.axes['lev'].points
-    if Tatm.shape == lev.shape:
-        # 1D interpolation on pressure levels using global average data
-        f = interp1d(ozone_lev, ozone_global)
-        #  interpolate data to model levels
-        absorber_vmr['O3'] = f(lev)
-    else:
-        #  Attempt 2D interpolation in pressure and latitude
-        f2d = interp2d(ozone_lat, ozone_lev, ozone_zon)
-        try:
-            lat = Tatm.domain.axes['lat'].points
+    # Ozone: start with all zeros, interpolate to data if we can
+    absorber_vmr['O3'] = np.zeros_like(Tatm)
+    if ozone_file is not None:
+        datadir = os.path.join(os.path.dirname(__file__), 'data', 'ozone')
+        ozonefilepath = os.path.join(datadir, ozone_file)
+        #  Open the ozone data file
+        if verbose:
+            print('Getting ozone data from', ozonefilepath)
+        ozonedata = nc.Dataset(ozonefilepath)
+        ozone_lev = ozonedata.variables['lev'][:]
+        ozone_lat = ozonedata.variables['lat'][:]
+        #  zonal and time average
+        ozone_zon = np.mean(ozonedata.variables['OZONE'], axis=(0,3))
+        ozone_global = np.average(ozone_zon,
+                    weights=np.cos(np.deg2rad(ozone_lat)), axis=1)
+        lev = Tatm.domain.axes['lev'].points
+        if Tatm.shape == lev.shape:
+            # 1D interpolation on pressure levels using global average data
+            f = interp1d(ozone_lev, ozone_global)
+            #  interpolate data to model levels
+            absorber_vmr['O3'] = f(lev)
+        else:
+            #  Attempt 2D interpolation in pressure and latitude
             f2d = interp2d(ozone_lat, ozone_lev, ozone_zon)
-            absorber_vmr['O3'] = f2d(lat, lev).transpose()
-        except:
-            print('Interpolation of ozone data failed.')
-            print('Reverting to default O3.')
-            absorber_vmr['O3'] = np.zeros_like(Tatm)
+            try:
+                lat = Tatm.domain.axes['lat'].points
+                f2d = interp2d(ozone_lat, ozone_lev, ozone_zon)
+                absorber_vmr['O3'] = f2d(lat, lev).transpose()
+            except:
+                print('Interpolation of ozone data failed.')
+                print('Reverting to default O3.')
+                absorber_vmr['O3'] = np.zeros_like(Tatm)
     return absorber_vmr
 
 def init_interface(field):
diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml
index 3302166..71f308f 100644
--- a/conda-recipe/meta.yaml
+++ b/conda-recipe/meta.yaml
@@ -1,5 +1,5 @@
 {% set name = "climlab" %}
-{% set version = "0.7.1.dev0" %}
+{% set version = "0.7.1.dev1" %}
 
 package:
   name: {{ name|lower }}
diff --git a/setup.py b/setup.py
index 03f544d..a632935 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 import os, sys
 import textwrap
 
-VERSION = '0.7.1.dev0'
+VERSION = '0.7.1.dev1'
 
 # BEFORE importing setuptools, remove MANIFEST. Otherwise it may not be
 # properly updated when the contents of directories change (true for distutils,


Test patch:
diff --git a/climlab/tests/test_rrtm.py b/climlab/tests/test_rrtm.py
index 9613b19..e3737c9 100644
--- a/climlab/tests/test_rrtm.py
+++ b/climlab/tests/test_rrtm.py
@@ -150,3 +150,19 @@ def test_latitude():
     grad = np.diff(model.Ts, axis=0)
     assert np.all(grad[0:(int(num_lat/2)-1)] > 0.)
     assert np.all(grad[int(num_lat/2):] < 0.)
+
+@pytest.mark.compiled
+@pytest.mark.fast
+def test_no_ozone():
+    '''When user gives None as the ozone_file, the model is initialized
+    with zero ozone. This should work on arbitrary grids.'''
+    ps = 1060.
+    num_lev=4000
+    state = climlab.column_state(num_lev=num_lev, num_lat=1, water_depth=5.)
+    lev = state.Tatm.domain.lev
+    lev.bounds = np.linspace(0., ps, num_lev+1)
+    lev.points = lev.bounds[:-1] + np.diff(lev.bounds)/2.
+    lev.delta = np.abs(np.diff(lev.bounds))
+    #  Create a RRTM radiation model
+    rad = climlab.radiation.RRTMG(state=state, ozone_file=None)
+    assert np.all(rad.absorber_vmr['O3']==0.)