From 0169ebdaf172eec33442bfea89b741cfe9de9659 Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Fri, 28 Jul 2017 10:21:59 +0200
Subject: [PATCH] stability.py added

---
 wetb/wind/stability.py | 118 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 118 insertions(+)
 create mode 100644 wetb/wind/stability.py

diff --git a/wetb/wind/stability.py b/wetb/wind/stability.py
new file mode 100644
index 00000000..2af9ea8e
--- /dev/null
+++ b/wetb/wind/stability.py
@@ -0,0 +1,118 @@
+'''
+Created on 05/10/2015
+
+@author: MMPE
+'''
+
+import numpy as np
+from scipy.integrate import trapz
+from scipy.signal.signaltools import detrend
+
+
+def MoninObukhov_length(u,v,w, T):
+    """Calculate the Monin Obukhov length
+
+    Not validated!!!
+
+    parameters
+    ----------
+    u : array_like
+        Horizontal wind fluctuations in mean wind direction
+    v : array_like
+        Horizontal wind fluctuations in perpendicular to mean wind direction
+    w : array_like
+        Vertical wind fluctuations
+    T : array_like
+        Potential temperature (close to sonic temperature)
+    """
+    K = 0.4
+    g = 9.82
+    u = detrend(u)
+    u = u-np.mean(u)
+    v = v-np.mean(v)
+    w = w-np.mean(w)
+    u_star = (np.mean(u*w)**2+np.mean(v*w)**2)**(1/4)
+    t = T-T.mean()
+    wT = np.mean(w*T)
+    return -u_star ** 3 * (T.mean() + 273.15) / (K * g * wT)
+
+
+def L2category(L, full_category_name=False):
+    """Stability category from Monin-Obukhov length
+    
+    Parameters
+    ----------
+    L : float or int
+        Monin-Obukhov length
+    full_category_name : bool, optional
+        If False, default, category ids are returned, e.g. "n" for neutral
+        If True, full name of category are returned, e.g. "Neutral"
+    
+    Returns
+    -------
+    Stability category : str
+    
+        Examples
+    --------
+    >>> L2category(1000)
+    n 
+    """
+    cat_limits = np.array([-50,-100,-200,-500,500,200,50,10])
+    index = np.searchsorted( 1/cat_limits, 1/np.array(L))-1
+    if full_category_name:
+        return np.array(['Very unstable','Unstable','Near unstable','Neutral','Near stable','Stable','Very stable','undefined'])[index]
+    else:
+        return np.array(['vu','u','nu','n','ns','s','vs','-'])[index]
+    
+def MoninObukhov_length2(u_star, w, T, specific_humidity=None):
+    """Calculate the Monin Obukhov length
+
+    Not validated!!!
+
+    parameters
+    ----------
+    u_star :
+        Refencence velocity at hub height
+    w : array_like
+        Vertical wind fluctuations
+    T : array_like
+        Temperature in celcius
+    """
+    K = 0.4
+    g = 9.82
+    if specific_humidity is not None:
+        potential_temperature = (w * (T + 273.15)).mean() + 0.61 * (T + 273.15).mean() * (w * specific_humidity).mean()
+    else:
+        potential_temperature = (w * (T + 273.15)).mean()
+    return -u_star ** 3 * (T.mean() + 273.15) / (K * g * potential_temperature)
+
+
+def humidity_relative2specific(relative_humidity, T, P):
+    """Not validated
+    parameters
+    ----------
+    relative_humidity : float
+        Relative humidity [%]
+    T : float
+        Temperature [C]
+    P : float
+        Barometric pressure [Pa]
+    """
+    return relative_humidity * np.exp(17.67 * T / (T + 273.15 - 29.65)) / 0.263 / P
+
+def humidity_specific2relative2(specific_humidity, T, P):
+    """Not validated
+    parameters
+    ----------
+    specific_humidity : float
+        specific_humidity [kg/kg]
+    T : float
+        Temperature [C]
+    P : float
+        Barometric pressure [Pa]
+    """
+    return 0.263 * P * specific_humidity / np.exp(17.67 * T / (T + 273.15 - 29.65))
+
+if __name__ == "__main__":
+    print (humidity_relative2specific(85, 8, 101325))
+    print (humidity_specific2relative2(5.61, 8, 101325))
-- 
GitLab