Merge branch 'development' into 'master'

Development

Closes #43

See merge request sbiri/orchestra!2

Code/AirSeaFluxCode_test.py

+import importlib
+import sys
+sys.path.insert(1, '/Users/ricorne/projects/orchestra/Code')
+import   AirSeaFluxCode
+importlib.reload(AirSeaFluxCode)
+from AirSeaFluxCode import *
+import pickle
+
+inDt = pd.read_csv("~/projects/orchestra/Test_Data/data_all.csv")
+date = np.asarray(inDt["Date"])
+lon = np.asarray(inDt["Longitude"])
+lat = np.asarray(inDt["Latitude"])
+spd = np.asarray(inDt["Wind speed"])
+t = np.asarray(inDt["Air temperature"])
+sst = np.asarray(inDt["SST"])
+rh = np.asarray(inDt["RH"])
+p = np.asarray(inDt["P"])
+sw = np.asarray(inDt["Rs"])
+hu = np.asarray(inDt["zu"])
+ht = np.asarray(inDt["zt"])
+hin = np.array([hu, ht, ht])
+Rs = np.asarray(inDt["Rs"])
+
+del hu, ht, inDt
+# run AirSeaFluxCode
+res1 = AirSeaFluxCode(spd, t, sst, lat=lat, hin=hin, P=p, maxiter=10,hum=None,cskin=0,wl=0,
+                      tol=['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1], L="Rb",meth="S80",gust=None)
+
+# old version
+pickle_off = open ("/Users/ricorne/projects/AirSeaFluxCode_master/orchestra/old_code.txt", "rb")
+res = pickle.load(pickle_off)
+
+for i in res1.columns:
+    try:
+        a=res1[i].round(4)
+        b=res[i].round(4)
+        print(a.equals(b))
+    except:
+        print(res1[i].equals(res[i]))

Code/cs_wl_subs.py

+import numpy as np
+from util_subs import (CtoK, kappa)
+
+def cs_C35(sst, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, g):
+    """
+    Computes cool skin
+    following COARE3.5 (Fairall et al. 1996, Edson et al. 2013)
+
+    Parameters
+    ----------
+    sst : float
+        sea surface temperature      [K]
+    qsea : float
+        specific humidity over sea   [g/kg]
+    rho : float
+        density of air               [kg/m^3]
+    Rs : float
+        downward shortwave radiation [Wm-2]
+    Rnl : float
+        net upwelling IR radiation       [Wm-2]
+    cp : float
+       specific heat of air at constant pressure [J/K/kg]
+    lv : float
+       latent heat of vaporization   [J/kg]
+    delta : float
+       cool skin thickness           [m]
+    usr : float
+       friction velocity             [m/s]
+    tsr : float
+       star temperature              [K]
+    qsr : float
+       star humidity                 [g/kg]
+    g : float
+       gravity                      [ms^-2]
+
+    Returns
+    -------
+    dter : float
+        cool skin correction         [K]
+    dqer : float
+       humidity corrction            [g/kg]
+    delta : float
+       cool skin thickness           [m]
+    """
+    # coded following Saunders (1967) with lambda = 6
+    if (np.nanmin(sst) > 200):  # if sst in Kelvin convert to Celsius
+        sst = sst-CtoK
+    # ************  cool skin constants  *******
+    # density of water, specific heat capacity of water, water viscosity,
+    # thermal conductivity of water
+    rhow, cpw, visw, tcw = 1022, 4000, 1e-6, 0.6
+    for i in range(4):
+        aw = 2.1e-5*np.power(np.maximum(sst+3.2, 0), 0.79)
+        bigc = 16*g*cpw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(
+            rho, 2))
+        Rns = 0.945*Rs  # albedo correction
+        shf = rho*cp*usr*tsr
+        lhf = rho*lv*usr*qsr
+        Qnsol = shf+lhf+Rnl
+        fs = 0.065+11*delta-6.6e-5/delta*(1-np.exp(-delta/8.0e-4))
+        qcol = Qnsol+Rns*fs
+        alq = aw*qcol+0.026*np.minimum(lhf, 0)*cpw/lv
+        xlamx = 6*np.ones(sst.shape)
+        xlamx = np.where(alq > 0, 6/(1+(bigc*alq/usr**4)**0.75)**0.333, 6)
+        delta = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
+        delta = np.where(alq > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), delta)
+    dter = qcol*delta/tcw
+    return dter, delta
+# ---------------------------------------------------------------------
+
+
+def delta(aw, Q, usr, g):
+    """
+    Computes the thickness (m) of the viscous skin layer.
+    Based on Fairall et al., 1996 and cited in IFS Documentation Cy46r1
+    eq. 8.155 p. 164
+
+    Parameters
+    ----------
+    aw : float
+        thermal expansion coefficient of sea-water  [1/K]
+    Q : float
+        part of the net heat flux actually absorbed in the warm layer [W/m^2]
+    usr : float
+        friction velocity in the air (u*) [m/s]
+    g : float
+       gravity                      [ms^-2]
+
+    Returns
+    -------
+    delta : float
+        the thickness (m) of the viscous skin layer
+
+    """
+    rhow, visw, tcw = 1025, 1e-6, 0.6
+    # u* in the water
+    usr_w = np.maximum(usr, 1e-4)*np.sqrt(1.2/rhow)  # rhoa=1.2
+    rcst_cs = 16*g*np.power(visw, 3)/np.power(tcw, 2)
+    lm = 6*(1+np.maximum(Q*aw*rcst_cs/np.power(usr_w, 4), 0)**0.75)**(-1/3)
+    ztmp = visw/usr_w
+    delta = np.where(Q > 0, np.minimum(6*ztmp, 0.007), lm*ztmp)
+    return delta
+# ---------------------------------------------------------------------
+
+# version that doesn't output dqer or import/output delta
+# should be very similar to cs_ecmwf
+# def cs_C35(sst, rho, Rs, Rnl, cp, lv, usr, tsr, qsr, g):
+#     """
+#     Compute cool skin.
+
+#     Based on COARE3.5 (Fairall et al. 1996, Edson et al. 2013)
+
+#     Parameters
+#     ----------
+#     sst : float
+#         sea surface temperature      [K]
+#     rho : float
+#         density of air               [kg/m^3]
+#     Rs : float
+#         downward shortwave radiation [Wm-2]
+#     Rnl : float
+#         net upwelling IR radiation       [Wm-2]
+#     cp : float
+#         specific heat of air at constant pressure [J/K/kg]
+#     lv : float
+#         latent heat of vaporization   [J/kg]
+#     delta : float
+#         cool skin thickness           [m]
+#     usr : float
+#         friction velocity             [ms^-1]
+#     tsr : float
+#         star temperature              [K]
+#     qsr : float
+#         star humidity                 [g/kg]
+#     g : float
+#         gravity                      [ms^-2]
+
+#     Returns
+#     -------
+#     dter : float
+#         cool skin correction         [K]
+#     delta : float
+#         cool skin thickness           [m]
+#     """
+#     # coded following Saunders (1967) with lambda = 6
+#     if (np.nanmin(sst) > 200):  # if sst in Kelvin convert to Celsius
+#         sst = sst-CtoK
+#     # ************  cool skin constants  *******
+#     # density of water, specific heat capacity of water, water viscosity,
+#     # thermal conductivity of water
+#     rhow, cpw, visw, tcw = 1022, 4000, 1e-6, 0.6
+#     aw = 2.1e-5*np.power(np.maximum(sst+3.2, 0), 0.79)
+#     bigc = 16*g*cpw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(
+#         rho, 2))
+#     Rns = 0.945*Rs  # albedo correction
+#     shf = rho*cp*usr*tsr
+#     lhf = rho*lv*usr*qsr
+#     Qnsol = shf+lhf+Rnl
+#     d = delta(aw, Qnsol, usr, g)
+#     for i in range(4):
+#         fs = 0.065+11*d-6.6e-5/d*(1-np.exp(-d/8.0e-4))
+#         qcol = Qnsol+Rns*fs
+#         alq = aw*qcol+0.026*np.minimum(lhf, 0)*cpw/lv
+#         xlamx = 6*np.ones(sst.shape)
+#         xlamx = np.where(alq > 0, 6/(1+(bigc*alq/usr**4)**0.75)**0.333, 6)
+#         d = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
+#         d = np.where(alq > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), d)
+#     dter = qcol*d/tcw
+#     return dter
+# ---------------------------------------------------------------------
+
+
+def cs_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, g):
+    """
+    cool skin adjustment based on IFS Documentation cy46r1
+
+    Parameters
+    ----------
+    rho : float
+        density of air               [kg/m^3]
+    Rs : float
+        downward solar radiation [Wm-2]
+    Rnl : float
+        net thermal radiation     [Wm-2]
+    cp : float
+       specific heat of air at constant pressure [J/K/kg]
+    lv : float
+       latent heat of vaporization   [J/kg]
+    usr : float
+       friction velocity         [m/s]
+    tsr : float
+       star temperature              [K]
+    qsr : float
+       star humidity                 [g/kg]
+    sst : float
+        sea surface temperature  [K]
+    g : float
+       gravity                      [ms^-2]
+
+    Returns
+    -------
+    dtc : float
+        cool skin temperature correction [K]
+
+    """
+    if (np.nanmin(sst) < 200):  # if sst in Celsius convert to Kelvin
+        sst = sst+CtoK
+    aw = np.maximum(1e-5, 1e-5*(sst-CtoK))
+    Rns = 0.945*Rs  # (net solar radiation (albedo correction)
+    shf = rho*cp*usr*tsr
+    lhf = rho*lv*usr*qsr
+    Qnsol = shf+lhf+Rnl  # eq. 8.152
+    d = delta(aw, Qnsol, usr, g)
+    for jc in range(4):  # because implicit in terms of delta...
+        # # fraction of the solar radiation absorbed in layer delta eq. 8.153
+        # and Eq.(5) Zeng & Beljaars, 2005
+        fs = 0.065+11*d-6.6e-5/d*(1-np.exp(-d/8e-4))
+        Q = Qnsol+fs*Rns
+        d = delta(aw, Q, usr, g)
+    dtc = Q*d/0.6  # (rhow*cw*kw)eq. 8.151
+    return dtc
+# ---------------------------------------------------------------------
+
+
+def wl_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, skt, dtc, g):
+    """
+    warm layer correction following IFS Documentation cy46r1
+    and aerobulk (Brodeau et al., 2016)
+
+    Parameters
+    ----------
+    rho : float
+        density of air               [kg/m^3]
+    Rs : float
+        downward solar radiation    [Wm-2]
+    Rnl : float
+        net thermal radiation  [Wm-2]
+    cp : float
+       specific heat of air at constant pressure [J/K/kg]
+    lv : float
+       latent heat of vaporization   [J/kg]
+    usr : float
+        friction velocity           [m/s]
+    tsr : float
+       star temperature              [K]
+    qsr : float
+       star humidity                 [g/kg]
+    sst : float
+        bulk sst                    [K]
+    skt : float
+        skin sst from previous step [K]
+    dtc : float
+        cool skin correction        [K]
+    g : float
+       gravity                      [ms^-2]
+
+    Returns
+    -------
+    dtwl : float
+        warm layer correction       [K]
+
+    """
+    if (np.nanmin(sst) < 200):  # if sst in Celsius convert to Kelvin
+        sst = sst+CtoK
+    rhow, cpw, visw, rd0 = 1025, 4190, 1e-6, 3
+    Rns = 0.945*Rs
+    #  Previous value of dT / warm-layer, adapted to depth:
+    aw = 2.1e-5*np.power(np.maximum(sst-CtoK+3.2, 0), 0.79)  # thermal expansion coefficient of sea-water (SST accurate enough#)
+    # *** Rd = Fraction of solar radiation absorbed in warm layer (-)
+    a1, a2, a3 = 0.28, 0.27, 0.45
+    b1, b2, b3 = -71.5, -2.8, -0.06  # [m-1]
+    Rd = 1-(a1*np.exp(b1*rd0)+a2*np.exp(b2*rd0)+a3*np.exp(b3*rd0))
+    shf = rho*cp*usr*tsr
+    lhf = rho*lv*usr*qsr
+    Qnsol = shf+lhf+Rnl
+    usrw = np.maximum(usr, 1e-4)*np.sqrt(1.2/rhow)   # u* in the water
+    zc3 = rd0*kappa*g/np.power(1.2/rhow, 3/2)
+    zc4 = (0.3+1)*kappa/rd0
+    zc5 = (0.3+1)/(0.3*rd0)
+    for jwl in range(10):  # itteration to solve implicitely equation for warm layer
+        dsst = skt-sst-dtc
+        # Buoyancy flux and stability parameter (zdl = -z/L) in water
+        ZSRD = (Qnsol+Rns*Rd)/(rhow*cpw)
+        ztmp = np.maximum(dsst, 0)
+        zdl = np.where(ZSRD > 0, 2*(np.power(usrw, 2) *
+                                    np.sqrt(ztmp/(5*rd0*g*aw/visw)))+ZSRD,
+                       np.power(usrw, 2)*np.sqrt(ztmp/(5*rd0*g*aw/visw)))
+        usr = np.maximum(usr, 1e-4 )
+        zdL = zc3*aw*zdl/np.power(usr, 3)
+        # Stability function Phi_t(-z/L) (zdL is -z/L) :
+        zphi = np.where(zdL > 0, (1+(5*zdL+4*np.power(zdL, 2)) /
+                                  (1+3*zdL+0.25*np.power(zdL, 2)) +
+                                  2/np.sqrt(1-16*(-np.abs(zdL)))),
+                        1/np.sqrt(1-16*(-np.abs(zdL))))
+        zz = zc4*(usrw)/zphi
+        zz = np.maximum(np.abs(zz), 1e-4)*np.sign(zz)
+        dtwl = np.maximum(0, (zc5*ZSRD)/zz)
+    return dtwl
+# ----------------------------------------------------------------------
+
+
+def cs_Beljaars(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, g, Qs):
+    """
+    cool skin adjustment based on Beljaars (1997)
+    air-sea interaction in the ECMWF model
+
+    Parameters
+    ----------
+    rho : float
+        density of air           [kg/m^3]
+    Rs : float
+        downward solar radiation [Wm-2]
+    Rnl : float
+        net thermal radiaion     [Wm-2]
+    cp : float
+       specific heat of air at constant pressure [J/K/kg]
+    lv : float
+       latent heat of vaporization   [J/kg]
+    usr : float
+       friction velocity         [m/s]
+    tsr : float
+       star temperature              [K]
+    qsr : float
+       star humidity                 [g/kg]
+    g : float
+       gravity                      [ms^-2]
+    Qs : float
+      radiation balance
+
+    Returns
+    -------
+    Qs : float
+      radiation balance
+    dtc : float
+        cool skin temperature correction [K]
+
+    """
+    tcw = 0.6       # thermal conductivity of water (at 20C) [W/m/K]
+    visw = 1e-6     # kinetic viscosity of water [m^2/s]
+    rhow = 1025     # Density of sea-water [kg/m^3]
+    cpw = 4190      # specific heat capacity of water
+    aw = 3e-4       # thermal expansion coefficient [K-1]
+    Rns = 0.945*Rs  # net solar radiation (albedo correction)
+    shf = rho*cp*usr*tsr
+    lhf = rho*lv*usr*qsr
+    Q = Rnl+shf+lhf+Qs
+    xt = (16*Q*g*aw*cpw*np.power(rhow*visw, 3))/(np.power(usr, 4) *
+                                                 np.power(rho*tcw, 2))
+    xt1 = 1+np.power(xt, 3/4)
+    # Saunders const  eq. 22
+    ls = np.where(Q > 0, 6/np.power(xt1, 1/3), 6)
+    delta = np.where(Q > 0, (ls*visw)/(np.sqrt(rho/rhow)*usr),
+                     np.where((ls*visw)/(np.sqrt(rho/rhow)*usr) > 0.01, 0.01,
+                              (ls*visw)/(np.sqrt(rho/rhow)*usr)))  # eq. 21
+    # fraction of the solar radiation absorbed in layer delta
+    fc = 0.065+11*delta-6.6e-5*(1-np.exp(-delta/0.0008))/delta
+    Qs = fc*Rns
+    Q = Rnl+shf+lhf+Qs
+    dtc = Q*delta/tcw
+    return Qs, dtc
+# ----------------------------------------------------------------------
+
+
+def get_dqer(dter, sst, qsea, lv):
+    """
+    Calculate humidity correction.
+
+    Parameters
+    ----------
+    dter : float
+        cool skin correction         [K]
+    sst : float
+        sea surface temperature      [K]
+    qsea : float
+        specific humidity over sea   [g/kg]
+    lv : float
+       latent heat of vaporization   [J/kg]
+
+    Returns
+    -------
+    dqer : float
+       humidity correction            [g/kg]
+
+    """
+    if (np.nanmin(sst) < 200):  # if sst in Celsius convert to Kelvin
+        sst = sst+CtoK
+    wetc = 0.622*lv*qsea/(287.1*np.power(sst, 2))
+    dqer = wetc*dter
+    return dqer
+

Code/hum_subs.py

@@ -54,7 +54,8 @@ def VaporPressure(temp, P, phase, meth):
     if (np.nanmin(temp) > 200):  # if Ta in Kelvin convert to Celsius
         temp = temp-273.16
     T = np.copy(temp)+273.16    #  Most formulas use T in [K]
-                                #  Formulas using [C] use the variable temp
+    
+    #  Formulas using [C] use the variable temp
     #  Calculate saturation pressure over liquid water
     if (phase == 'liquid'):
         if (meth == 'HylandWexler' or meth == ''):