From 5234503d09974bcf1328f47dd134e57f448a0dc4 Mon Sep 17 00:00:00 2001
From: "J.R. Leeman" <kd5wxb@gmail.com>
Date: Thu, 26 Oct 2017 09:19:11 -0600
Subject: [PATCH 1/4] Cleanup precipitable water function.

* Correct error in formula documentation.
* Add bottom pressure argument.
* Make default of top pressure None.
---
 metpy/calc/indices.py | 25 +++++++++++++++++--------
 metpy/calc/tools.py   |  6 +++++-
 2 files changed, 22 insertions(+), 9 deletions(-)

diff --git a/metpy/calc/indices.py b/metpy/calc/indices.py
index 526d90b362e..b5401a6f648 100644
--- a/metpy/calc/indices.py
+++ b/metpy/calc/indices.py
@@ -15,12 +15,12 @@
 
 @exporter.export
 @check_units('[temperature]', '[pressure]', '[pressure]')
-def precipitable_water(dewpt, pressure, top=400 * units('hPa')):
+def precipitable_water(dewpt, pressure, bottom=None, top=None):
     r"""Calculate precipitable water through the depth of a sounding.
 
-    Default layer depth is sfc-400 hPa. Formula used is:
+    Formula used is:
 
-    .. math:: \frac{1}{pg} \int\limits_0^d x \,dp
+    .. math::  -\frac{1}{\rho_l g} \int\limits_{p_\text{bottom}}^{p_\text{top}} r dp
 
     from [Tsonis2008]_, p. 170.
 
@@ -30,8 +30,10 @@ def precipitable_water(dewpt, pressure, top=400 * units('hPa')):
         Atmospheric dewpoint profile
     pressure : `pint.Quantity`
         Atmospheric pressure profile
+    bottom: `pint.Quantity`, optional
+        Bottom of the layer, specified in pressure. Defaults to None (highest pressure).
     top: `pint.Quantity`, optional
-        The top of the layer, specified in pressure. Defaults to 400 hPa.
+        The top of the layer, specified in pressure. Defaults to None (lowest pressure).
 
     Returns
     -------
@@ -39,15 +41,22 @@ def precipitable_water(dewpt, pressure, top=400 * units('hPa')):
         The precipitable water in the layer
 
     """
-    sort_inds = np.argsort(pressure[::-1])
+    # Sort pressure and dewpoint to be in decreasing pressure order (increasing height)
+    sort_inds = np.argsort(pressure)[::-1]
     pressure = pressure[sort_inds]
     dewpt = dewpt[sort_inds]
 
-    pres_layer, dewpt_layer = get_layer(pressure, dewpt, depth=pressure[0] - top)
+    if top is None:
+        top = np.nanmin(pressure) * pressure.units
+
+    if bottom is None:
+        bottom = np.nanmax(pressure) * pressure.units
+
+    pres_layer, dewpt_layer = get_layer(pressure, dewpt, bottom=bottom, depth=bottom - top)
 
     w = mixing_ratio(saturation_vapor_pressure(dewpt_layer), pres_layer)
-    # Since pressure is in decreasing order, pw will be the negative of what we want.
-    # Thus the *-1
+
+    # Since pressure is in decreasing order, pw will be the opposite sign of that expected.
     pw = -1. * (np.trapz(w.magnitude, pres_layer.magnitude) * (w.units * pres_layer.units) /
                 (g * rho_l))
     return pw.to('millimeters')
diff --git a/metpy/calc/tools.py b/metpy/calc/tools.py
index 99adc72ef18..d9b0af22f37 100644
--- a/metpy/calc/tools.py
+++ b/metpy/calc/tools.py
@@ -323,6 +323,10 @@ def _get_bound_pressure_height(pressure, bound, heights=None, interpolate=True):
         The bound pressure and height.
 
     """
+    sort_inds = np.argsort(pressure)[::-1]
+    pressure = pressure[sort_inds]
+    if heights is not None:
+        heights = heights[sort_inds]
     # Bound is given in pressure
     if bound.dimensionality == {'[length]': -1.0, '[mass]': 1.0, '[time]': -2.0}:
         # If the bound is in the pressure data, we know the pressure bound exactly
@@ -538,7 +542,7 @@ def get_layer(pressure, *args, **kwargs):
 
     # If the bottom is not specified, make it the surface pressure
     if bottom is None:
-        bottom = pressure[0]
+        bottom = np.nanmax(pressure) * pressure.units
 
     bottom_pressure, bottom_height = _get_bound_pressure_height(pressure, bottom,
                                                                 heights=heights,

From 49049b321aba95456ae9ecb73f13c53da2d10eef Mon Sep 17 00:00:00 2001
From: "J.R. Leeman" <kd5wxb@gmail.com>
Date: Thu, 26 Oct 2017 09:21:16 -0600
Subject: [PATCH 2/4] Update test to have top of 400 hPa (old default value).

---
 metpy/calc/tests/test_indices.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/metpy/calc/tests/test_indices.py b/metpy/calc/tests/test_indices.py
index 1183f5b30bc..0c1d34ef109 100644
--- a/metpy/calc/tests/test_indices.py
+++ b/metpy/calc/tests/test_indices.py
@@ -21,7 +21,8 @@ def test_precipitable_water():
     """Test precipitable water with observed sounding."""
     with UseSampleData():
         data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC', source='wyoming')
-    pw = precipitable_water(data.variables['dewpoint'][:], data.variables['pressure'][:])
+    pw = precipitable_water(data.variables['dewpoint'][:], data.variables['pressure'][:],
+                            top=400 * units.hPa)
     truth = (0.8899441949243486 * units('inches')).to('millimeters')
     assert_array_equal(pw, truth)
 

From b522ed926efb4c81e46ee0211026e675a9c25f40 Mon Sep 17 00:00:00 2001
From: "J.R. Leeman" <kd5wxb@gmail.com>
Date: Thu, 26 Oct 2017 13:45:35 -0600
Subject: [PATCH 3/4] Update reference to Salby 1996.

---
 docs/references.rst   | 3 ---
 metpy/calc/indices.py | 2 +-
 2 files changed, 1 insertion(+), 4 deletions(-)

diff --git a/docs/references.rst b/docs/references.rst
index 772bcaac47c..3dcab24cf9a 100644
--- a/docs/references.rst
+++ b/docs/references.rst
@@ -94,9 +94,6 @@ References
            composite and significant tornado parameters. Preprints, 22nd Conf. on Severe Local
            Storms, Hyannis, MA, Amer. Meteor. Soc.
 
-.. [Tsonis2008] Tsonis, A. A., 2008: An introduction to atmospheric thermodynamics.
-           Cambridge Univ. Press, Cambridge.
-
 .. [Ziv1994] Ziv, B., and Alpert, P., 1994: Isobaric to Isentropic Interpolation Errors
            and Implication to Potential Vorticity Analysis. *J. Appl. Meteor.*, **33**,
            694-703, doi:`10.1175/1520-0450(1994)033%3C0694:ITIIEA%3E2.0.CO;2
diff --git a/metpy/calc/indices.py b/metpy/calc/indices.py
index b5401a6f648..c5d84323517 100644
--- a/metpy/calc/indices.py
+++ b/metpy/calc/indices.py
@@ -22,7 +22,7 @@ def precipitable_water(dewpt, pressure, bottom=None, top=None):
 
     .. math::  -\frac{1}{\rho_l g} \int\limits_{p_\text{bottom}}^{p_\text{top}} r dp
 
-    from [Tsonis2008]_, p. 170.
+    from [Salby1996]_, p. 28.
 
     Parameters
     ----------

From 3c2f10d003535f24796aa47dc99bbdec463948a8 Mon Sep 17 00:00:00 2001
From: "J.R. Leeman" <kd5wxb@gmail.com>
Date: Thu, 26 Oct 2017 15:00:16 -0600
Subject: [PATCH 4/4] Add test for precipitable water with no bounds given.

---
 metpy/calc/tests/test_indices.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/metpy/calc/tests/test_indices.py b/metpy/calc/tests/test_indices.py
index 0c1d34ef109..5067fb70cb8 100644
--- a/metpy/calc/tests/test_indices.py
+++ b/metpy/calc/tests/test_indices.py
@@ -27,6 +27,18 @@ def test_precipitable_water():
     assert_array_equal(pw, truth)
 
 
+def test_precipitable_water_no_bounds():
+    """Test precipitable water with observed sounding and no bounds given."""
+    with UseSampleData():
+        data = get_upper_air_data(datetime(2016, 5, 22, 0), 'DDC', source='wyoming')
+    dewpoint = data.variables['dewpoint'][:]
+    pressure = data.variables['pressure'][:]
+    inds = pressure >= 400 * units.hPa
+    pw = precipitable_water(dewpoint[inds], pressure[inds])
+    truth = (0.8899441949243486 * units('inches')).to('millimeters')
+    assert_array_equal(pw, truth)
+
+
 def test_mean_pressure_weighted():
     """Test pressure-weighted mean wind function with vertical interpolation."""
     with UseSampleData():