add option 5 for gustiness following UA and shift C35 code to option 6

Add all function (unified, not unified) in cs_wl_subs.py

AirSeaFluxCode/src/AirSeaFluxCode.py

@@ -10,14 +10,14 @@ from cs_wl_subs import *
 
 class S88:
     def _wind_iterate(self, ind):
-        if self.gust[0] in range(1, 6):
+        if self.gust[0] in range(1, 7):
             self.wind[ind] = np.sqrt(np.power(np.copy(self.spd[ind]), 2) +
                                      np.power(get_gust(
                                          self.gust[1], self.gust[2],
                                          self.gust[3], self.theta[ind],
                                          self.usr[ind], self.tsrv[ind],
                                          self.grav[ind]), 2))
-            if self.gust[0] in [3, 4]:
+            if self.gust[0] in [3, 4, 5]:
                 # self.GustFact[ind] = 1
                 # option to not remove GustFact
                 self.u10n[ind] = self.wind[ind]-self.usr[ind]/kappa*(
@@ -92,11 +92,6 @@ class S88:
                     self.SST[ind]), self.rho[ind], self.Rs[ind], self.Rnl[ind],
                     self.cp[ind], self.lv[ind], np.copy(self.tkt[ind]),
                     self.usr[ind], self.tsr[ind], self.qsr[ind], self.grav[ind])
-
-                # self.dter[ind] = cs_C35(np.copy(
-                #     self.skt[ind]), self.rho[ind], self.Rs[ind], self.Rnl[ind],
-                #     self.cp[ind], self.lv[ind], self.usr[ind], self.tsr[ind],
-                #     self.qsr[ind], self.grav[ind])
             elif self.skin == "ecmwf":
                 self.dter[ind] = cs_ecmwf(
                     self.rho[ind], self.Rs[ind], self.Rnl[ind], self.cp[ind],
@@ -406,6 +401,12 @@ class S88:
             self.uref = self.wind-self.usr/kappa*(
                 np.log(self.h_in[0]/self.h_out[0])-self.psim +
                 psim_calc(self.h_out[0]/self.monob, self.meth))
+        elif self.gust[0] == 5:
+            self.u10n = self.spd-self.usr/kappa*(
+                np.log(self.h_in[0]/self.ref10)-self.psim)
+            self.uref = self.spd-self.usr/kappa*(
+                np.log(self.h_in[0]/self.h_out[0])-self.psim +
+                psim_calc(self.h_out[0]/self.monob, self.meth))
         else:
             self.u10n = self.spd-self.usr/kappa/self.GFo*(
                 np.log(self.h_in[0]/self.ref10)-self.psim) # C.4-7
@@ -527,7 +528,7 @@ class S88:
             gust = [0, 0, 0, 0]
 
         assert np.size(gust) == 4, "gust input must be a 4x1 array"
-        assert gust[0] in range(6), "gust at position 0 must be 0 to 5"
+        assert gust[0] in range(7), "gust at position 0 must be 0 to 6"
         self.gust = gust
 
     def _class_flag(self):

AirSeaFluxCode/src/AirSeaFluxCode_dev.py

@@ -10,7 +10,7 @@ from cs_wl_subs import *
 
 class S88:
     def _wind_iterate(self, ind):
-        if self.gust[0] in range(1, 6):
+        if self.gust[0] in range(1, 7):
             self.wind[ind] = np.sqrt(np.power(np.copy(self.spd[ind]), 2) +
                                      np.power(get_gust(
                                          self.gust[1], self.gust[2],
@@ -361,8 +361,6 @@ class S88:
                 self.tau = self.rho*np.power(self.usr, 2)*self.spd/self.wind
                 self.sensible = self.rho*self.cp*self.usr*self.tsr
                 self.latent = self.rho*self.lv*self.usr*self.qsr
-            
-
             # Set the new variables (for comparison against "old")
             new = np.array([np.copy(getattr(self, i)) for i in new_vars])
 
@@ -582,7 +580,7 @@ class S88:
             gust = [0, 0, 0, 0]
 
         assert np.size(gust) == 4, "gust input must be a 4x1 array"
-        assert gust[0] in range(6), "gust at position 0 must be 0 to 5"
+        assert gust[0] in range(7), "gust at position 0 must be 0 to 6"
         self.gust = gust
 
     def _class_flag(self):

AirSeaFluxCode/src/cs_wl_subs.py

 import numpy as np
 from util_subs import (CtoK, kappa)
 
+# def delta(aw, Q, usr, grav, rho, opt):
+#     """
+#     Compute the thickness (m) of the viscous skin layer.
 
-def delta(aw, Q, usr, grav, rho, opt):
-    """
-    Compute 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
+#     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]
-    grav : float
-       gravity                      [ms^-2]
+#     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]
+#     grav : float
+#        gravity                      [ms^-2]
 
-    Returns
-    -------
-    delta : float
-        the thickness (m) of the viscous skin layer
+#     Returns
+#     -------
+#     delta : float
+#         the thickness (m) of the viscous skin layer
 
-    """
-    if opt == "C35":
-        rhow, cw, visw, tcw = 1022, 4000, 1e-6, 0.6
-        # second term in eq. 14
-        tmp = 16*grav*Q*aw*rhow*cw*np.power(rhow*visw, 3)/(
-            np.power(usr, 4)*np.power(rho, 2)*np.power(tcw, 2))
-        # second term in eq. 14
-        # lm = 6*np.ones(usr.shape)
-        lm = np.where(Q > 0, 6/np.power(1+np.power(tmp, 3/4), 1/3), 6) # eq. 14 F96
-        delta = np.minimum(lm*visw/np.sqrt(rho/rhow)/usr, 0.01) # eq. 12 F96
-        delta = np.where(Q > 0, lm*visw/(np.sqrt(rho/rhow)*usr), delta)
-    elif opt == "ecmwf":
-        rhow, cw, visw, tcw = 1025, 4190, 1e-6, 0.6
-        # u* in the water
-        usr_w = np.maximum(usr, 1e-4)*np.sqrt(rho/rhow)  # rhoa=1.2
-        rcst_cs = 16*grav*np.power(visw, 3)/np.power(tcw, 2)#/(rhow*cw)
-        # eq. 8.155 Cy46r1
-        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
+#     """
+#     if opt == "C35":
+#         rhow, cw, visw, tcw = 1022, 4000, 1e-6, 0.6
+#         # second term in eq. 14
+#         tmp = 16*grav*Q*aw*rhow*cw*np.power(rhow*visw, 3)/(
+#             np.power(usr, 4)*np.power(rho, 2)*np.power(tcw, 2))
+#         # second term in eq. 14
+#         # lm = 6*np.ones(usr.shape)
+#         lm = np.where(Q > 0, 6/np.power(1+np.power(tmp, 3/4), 1/3), 6) # eq. 14 F96
+#         delta = np.minimum(lm*visw/np.sqrt(rho/rhow)/usr, 0.01) # eq. 12 F96
+#         delta = np.where(Q > 0, lm*visw/(np.sqrt(rho/rhow)*usr), delta)
+#         # bigc = 16*grav*cw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(
+#         #         rho, 2))
+#         # lm = np.where(Qb > 0, 6/(1+(bigc*Qb/usr**4)**0.75)**0.333, 6)
+#                 # d = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
+#                 # d = np.where(Qb > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), d)
+#     elif opt == "ecmwf":
+#         rhow, cw, visw, tcw = 1025, 4190, 1e-6, 0.6
+#         # u* in the water
+#         usr_w = np.maximum(usr, 1e-4)*np.sqrt(rho/rhow)  # rhoa=1.2
+#         rcst_cs = 16*grav*np.power(visw, 3)/np.power(tcw, 2)#/(rhow*cw)
+#         # eq. 8.155 Cy46r1
+#         # lm = 6/np.power(1+np.power(Q*aw*rcst_cs/np.power(usr_w, 4), 3/4), 1/3)
+#         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
 # ---------------------------------------------------------------------
 
 
@@ -122,6 +127,162 @@ def cs(sst, d, rho, Rs, Rnl, cp, lv, usr, tsr, qsr, grav, opt):
 # ---------------------------------------------------------------------
 
 
+def cs_C35(sst, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, grav):
+    """
+    Compute cool skin.
+
+    Based on 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]
+    grav : 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*grav*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))
+        Q = Qnsol+Rns*fs
+        Qb = aw*Q+0.026*np.minimum(lhf, 0)*cpw/lv
+        xlamx = 6*np.ones(sst.shape)
+        xlamx = np.where(Qb > 0, 6/(1+(bigc*Qb/usr**4)**0.75)**0.333, 6)
+        delta = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
+        delta = np.where(Qb > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), delta)
+    dter = Q*delta/tcw
+    return dter, delta
+# ---------------------------------------------------------------------
+
+
+def delta(aw, Q, usr, grav):
+    """
+    Compute 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]
+    grav : 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*grav*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
+# ---------------------------------------------------------------------
+
+
+def cs_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, grav):
+    """
+    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]
+    grav : 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, grav)
+    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, grav)
+    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, grav):
     """
     Calculate warm layer correction following IFS Documentation cy46r1.
@@ -200,6 +361,68 @@ def wl_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, skt, dtc, grav):
 # ----------------------------------------------------------------------
 
 
+def cs_Beljaars(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, grav, 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]
+    grav : 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*grav*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.

AirSeaFluxCode/src/flux_subs.py

@@ -633,7 +633,7 @@ def apply_GF(gust, spd, wind, step):
         GustFact = wind*0+1
         if gust[0] in [1, 2]:
             GustFact = np.sqrt(wind/spd)
-        elif gust[0] == 5:
+        elif gust[0] == 6:
             # as in C35 matlab code
             GustFact = wind/spd
     elif step == "TSF":

AirSeaFluxCode/src/flux_subs_dev.py

@@ -704,7 +704,7 @@ def apply_GF(gust, spd, wind, step):
         GustFact = wind*0+1
         if gust[0] in [1, 2]:
             GustFact = np.sqrt(wind/spd)
-        elif gust[0] == 5:
+        elif gust[0] == 6:
             # as in C35 matlab code
             GustFact = wind/spd
     elif step == "TSF":