Update AirSeaFluxCode.py, flux_subs.py, get_init.py, hum_subs.py, data_all.nc,...

Update AirSeaFluxCode.py, flux_subs.py, get_init.py, hum_subs.py, data_all.nc, Documentation.pdf, era5_r360x180.nc, run_ASFC.py files
Deleted data_Mxtr.nc, data_mod.nc, data_trp.nc, data_xtr.nc files

AirSeaFluxCode.py

@@ -53,8 +53,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
             "S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ERA5"
         qmeth : str
             is the saturation evaporation method to use amongst
-            "HylandWexler","Hardy","Preining","Wexler","GoffGratch","CIMO",
-            "MagnusTetens","Buck","Buck2","WMO","WMO2000","Sonntag","Bolton",
+            "HylandWexler","Hardy","Preining","Wexler","GoffGratch","WMO",
+            "MagnusTetens","Buck","Buck2","WMO2018","Sonntag","Bolton",
             "IAPWS","MurphyKoop"]
             default is Buck2
         tol : float
@@ -76,10 +76,10 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
     Returns
     -------
         res : array that contains
-                       1. momentum flux (W/m^2)
+                       1. momentum flux (N/m^2)
                        2. sensible heat (W/m^2)
                        3. latent heat (W/m^2)
-                       4. Monin-Obhukov length (mb)
+                       4. Monin-Obhukov length (m)
                        5. drag coefficient (cd)
                        6. neutral drag coefficient (cdn)
                        7. heat exhange coefficient (ct)
@@ -104,10 +104,15 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
                        26. temperature at reference height (tref)
                        27. specific humidity at reference height (qref)
                        28. number of iterations until convergence
-        ind : int
-            the indices in the matrix for the points that did not converge
-            after the maximum number of iterations
-    The code is based on bform.f and flux_calc.R modified by S. Biri
+                       29. cool-skin temperature depression (dter)
+                       30. cool-skin humidity depression (dqer)
+                       31. specific humidity of air (qair)
+                       32. specific humidity at sea surface (qsea)
+                       33. downward longwave radiation (Rl)
+                       34. downward shortwave radiation (Rs)
+                       35. downward net longwave radiation (Rnl)
+
+    2020 / Author S. Biri
     """
     logging.basicConfig(filename='flux_calc.log',
                         format='%(asctime)s %(message)s',level=logging.INFO)
@@ -135,6 +140,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
                   np.nanmedian(qsea), np.nanmedian(qair))
     if (np.all(np.isnan(qsea)) or np.all(np.isnan(qair))):
         print("qsea and qair cannot be nan")
+    # tlapse = gamma_moist(SST, T, qair/1000) lapse rate can be computed like this
     dt = Ta - sst
     dq = qair - qsea
     #  first guesses
@@ -289,7 +295,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
         elif (tol[0] == 'all'):
             new = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
                             np.copy(tau), np.copy(sensible), np.copy(latent)])
-        d = np.fabs(new-old)
+        d = np.abs(new-old)
         if (tol[0] == 'flux'):
             ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
                             (d[2, :] > tol[3]))
@@ -305,6 +311,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
         else:
             ii = True
     itera[ind] = -1
+    # itera = np.where(itera > n, -1, itera)
     logging.info('method %s | # of iterations:%s', meth, it)
     logging.info('method %s | # of points that did not converge :%s', meth,
                   ind[0].size)

flux_subs.py

@@ -32,7 +32,7 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
                        np.where((u10n <= 25) & (u10n >= 11),
                        (0.49+0.065*u10n)*0.001, 1.14*0.001))
     elif (meth == "S88" or meth == "UA" or meth == "ERA5" or meth == "C30" or
-          meth == "C35" or meth == "C40"):
+          meth == "C35" or meth == "C40" or meth == "Beljaars"):
         cdn = cdn_from_roughness(u10n, Ta, None, lat, meth)
     elif (meth == "YT96"):
         # for u<3 same as S80
@@ -101,7 +101,7 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
             a = 0.011*np.ones(Ta.shape)
             a = np.where(u10n > 22, 0.0016*22-0.0035, 0.0016*u10n-0.0035)
             zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr # surface roughness
-        elif (meth == "ERA5"):
+        elif ((meth == "ERA5" or meth == "Beljaars")):
             # eq. (3.26) p.38 over sea IFS Documentation cy46r1
             zo = 0.018*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
         else:
@@ -209,7 +209,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
 #                       5e-5/np.power(rr, 0.6))  # moisture roughness as in C30
         cqn = kappa**2/np.log(10/zo)/np.log(10/zoq)
         ctn = kappa**2/np.log(10/zo)/np.log(10/zot)
-    elif (meth == "ERA5"):
+    elif (meth == "ERA5" or meth == "Beljaars"):
         # eq. (3.26) p.38 over sea IFS Documentation cy46r1
         usr = np.sqrt(cdn*np.power(u10n, 2))
         zot = 0.40*visc_air(Ta)/usr
@@ -274,8 +274,8 @@ def get_stabco(meth="S80"):
     """
     alpha, beta, gamma = 0, 0, 0
     if (meth == "S80" or meth == "S88" or meth == "LY04" or
-        meth == "UA" or meth == "ERA5" or meth == "C30" or meth == "C35" or
-        meth == "C40"):
+        meth == "UA" or meth == "ERA5" or meth == "C30" or
+        meth == "C35" or meth == "C40" or meth == "Beljaars"):
         alpha, beta, gamma = 16, 0.25, 5  # Smith 1980, from Dyer (1974)
     elif (meth == "LP82"):
         alpha, beta, gamma = 16, 0.25, 7
@@ -308,6 +308,8 @@ def psim_calc(zol, meth="S80"):
         psim = psim_era5(zol)
     elif (meth == "C30" or meth == "C35" or meth == "C40"):
         psim = psiu_26(zol, meth)
+    elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
+        psim = np.where(zol < 0, psim_conv(zol, meth), psi_Bel(zol))
     else:
         psim = np.where(zol < 0, psim_conv(zol, meth),
                         psim_stab(zol, meth))
@@ -334,6 +336,8 @@ def psit_calc(zol, meth="S80"):
                         psi_era5(zol))
     elif (meth == "C30" or meth == "C35" or meth == "C40"):
         psit = psit_26(zol)
+    elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
+        psit = np.where(zol < 0, psi_conv(zol, meth), psi_Bel(zol))
     else:
         psit = np.where(zol < 0, psi_conv(zol, meth),
                         psi_stab(zol, meth))
@@ -341,6 +345,26 @@ def psit_calc(zol, meth="S80"):
 # ---------------------------------------------------------------------
 
 
+def psi_Bel(zol):
+    """ Calculates momentum/heat stability function
+
+    Parameters
+    ----------
+    zol : float
+        height over MO length
+    meth : str
+        parameterisation method
+
+    Returns
+    -------
+    psit : float
+    """
+    a, b, c, d = 0.7, 0.75, 5, 0.35
+    psi = -(a*zol+b*(zol-c/d)*np.exp(-d*zol)+b*c/d)
+    return psi
+# ---------------------------------------------------------------------
+
+
 def psi_era5(zol):
     """ Calculates heat stability function for stable conditions
         for method ERA5
@@ -578,8 +602,11 @@ def get_skin(sst, qsea, rho, Rl, Rs, Rnl, cp, lv, tkt, usr, tsr, qsr, lat):
     Returns
     -------
     dter : float
+       cool-skin temperature depression
     dqer : float
-
+       cool-skin humidity depression
+    tkt  : float
+       cool skin thickness
     """
     # coded following Saunders (1967) with lambda = 6
     g = gc(lat, None)
@@ -717,7 +744,7 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
                    (np.log((h_in[0]+zo)/zot) -
                     psit_calc((h_in[0]+zo)/monob, meth) +
                     psit_calc(zot/monob, meth))))
-        monob = h_in[0]/zol
+        monob = h_in[0]/zol        
     return tsrv, monob
 #------------------------------------------------------------------------------
 

get_init.py

@@ -73,6 +73,7 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
         MO length switch
 
     """
+    # check if input is correct (type, size, value) and set defaults
     if ((type(spd) != np.ndarray) or (type(T) != np.ndarray) or
          (type(SST) != np.ndarray)):
         sys.exit("input type of spd, T and SST should be numpy.ndarray")
@@ -82,11 +83,11 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
         sys.exit("input dtype of spd, T and SST should be float")
     # if input values are nan break
     if meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
-                    "C40","ERA5"]:
+                    "C40", "ERA5", "Beljaars"]:
         sys.exit("unknown method")
-    if qmeth not in ["HylandWexler", "Hardy", "Preining", "Wexler", "CIMO",
-                      "GoffGratch", "MagnusTetens", "Buck", "Buck2", "WMO",
-                      "WMO2000", "Sonntag", "Bolton", "IAPWS", "MurphyKoop"]:
+    if qmeth not in ["HylandWexler", "Hardy", "Preining", "Wexler",
+                     "GoffGratch", "WMO", "MagnusTetens", "Buck", "Buck2",
+                     "WMO2018", "Sonntag", "Bolton", "IAPWS", "MurphyKoop"]:
         sys.exit("unknown q-method")
     if (np.all(np.isnan(spd)) or np.all(np.isnan(T)) or np.all(np.isnan(SST))):
         sys.exit("input wind, T or SST is empty")
@@ -107,11 +108,13 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
                              meth == "LY04")):
         cskin = 0
     elif ((cskin == None) and (meth == "C30" or meth == "C35" or meth == "C40"
-                               or meth == "ERA5")):
+                               or meth == "ERA5" or meth == "Beljaars")):
         cskin = 1
-    if (np.all(gust == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
+    if (np.all(gust == None) and (meth == "C30" or meth == "C35" or
+                                  meth == "C40")):
         gust = [1, 1.2, 600]
-    elif (np.all(gust == None) and (meth == "UA" or meth == "ERA5")):
+    elif (np.all(gust == None) and (meth == "UA" or meth == "ERA5" or
+                                    meth == "Beljaars")):
         gust = [1, 1, 1000]
     elif np.all(gust == None):
         gust = [1, 1.2, 800]
@@ -124,7 +127,7 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
     if ((L == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
                          or meth == "YT96" or meth == "LY04" or
                          meth == "UA" or meth == "C30" or meth == "C35"
-                         or meth == "C40")):
+                         or meth == "C40" or meth == "Beljaars")):
         L = "S80"
     elif ((L == None) and (meth == "ERA5")):
         L = "ERA5"

hum_subs.py

@@ -24,19 +24,16 @@ def VaporPressure(temp, P, phase, meth):
         'ice' : Calculate vapor pressure over ice
     meth : str
         formula to be used
-        MartiMauersberger   : vaporpressure formula from Marti Mauersberger
+        Hardy               : vaporpressure formula from Hardy (1998)
         MagnusTetens        : vaporpressure formula from Magnus Tetens
         GoffGratch          : vaporpressure formula from Goff Gratch
-        Buck_original       : vaporpressure formula from Buck (original)
-        Buck_manual         : vaporpressure formula from the Buck manual
-        CIMO                : vaporpressure formula recommended by CIMO
-        Vaisala             : vaporpressure formula from Vaisala
-        WMO_Goff            : vaporpressure formula from WMO, which should have been Goff
-        WMO2000             : vaporpressure formula from WMO (2000) containing a typo
+        Buck                : vaporpressure formula from Buck (1981)
+        Buck2               : vaporpressure formula from the Buck (2012)
+        WMO                 : vaporpressure formula from WMO (1988)
+        WMO2018             : vaporpressure formula from WMO (2018)
         Wexler              : vaporpressure formula from Wexler (1976)
         Sonntag             : vaporpressure formula from Sonntag (1994)
         Bolton              : vaporpressure formula from Bolton (1980)
-        Fukuta              : vaporpressure formula from Fukuta (2003)
         HylandWexler        : vaporpressure formula from Hyland and Wexler (1983)
         IAPWS               : vaporpressure formula from IAPWS (2002)
         Preining            : vaporpressure formula from Preining (2002)
@@ -111,12 +108,6 @@ def VaporPressure(temp, P, phase, meth):
                             1.3816e-7*(np.power(10, 11.344*(1-T/Ts))-1) +
                             8.1328e-3*(np.power(10, -3.49149*(Ts/T-1))-1) +
                             np.log10(ews))
-        if (meth == 'CIMO'):
-            """Source: Annex 4B, Guide to Meteorological Instruments and
-            Methods of Observation, WMO Publication No 8, 7th edition, Geneva,
-            2008. (CIMO Guide)"""
-            Psat = (6.112*np.exp(17.62*temp/(243.12+temp)) *
-                    (1.0016+3.15e-6*P-0.074/P))
         if (meth == 'MagnusTetens'):
             """Source: Murray, F. W., On the computation of \
                          saturation vapor pressure, J. Appl. Meteorol., \
@@ -134,7 +125,7 @@ def VaporPressure(temp, P, phase, meth):
         if (meth == 'Buck2'):
             """Bucks vapor pressure formulation based on Tetens formula.
             Source: Buck Research, Model CR-1A Hygrometer Operating Manual,
-            Sep 2001"""
+            May 2012"""
             Psat = (6.1121*np.exp((18.678-(temp)/234.5)*(temp)/(257.14+temp)) *
                     (1+1e-4*(7.2+P*(0.0320)+5.9e-6*np.power(T, 2))))
         if (meth == 'WMO'):
@@ -147,18 +138,12 @@ def VaporPressure(temp, P, phase, meth):
             Ts = 273.16  # triple point temperature in K
             Psat = np.power(10, 10.79574*(1-Ts/T)-5.028*np.log10(T/Ts) +
                             1.50475e-4*(1-np.power(10, -8.2969*(T/Ts-1))) +
-                            0.42873e-3*(np.power(10, -4.76955*(1-Ts/T))-1) +
-                            0.78614)
-        if (meth == 'WMO2000'):
-            """WMO formulation, which is very similar to Goff & Gratch.
-            Source : WMO technical regulations, WMO-NO 49, Vol I, General
-            Meteorological Standards and Recommended Practices, App. A,
-            Corrigendum Aug 2000."""
-            Ts = 273.16  # triple point temperature in K
-            Psat = np.power(10, 10.79574*(1-Ts/T)-5.028*np.log10(T/Ts) +
-                            1.50475e-4*(1-np.power(10, -8.2969*(T/Ts-1))) +
-                            0.42873e-3*(np.power(10, -4.76955*(1.-Ts/T))-1) +
+                            0.42873e-3*(np.power(10, 4.76955*(1-Ts/T))-1) +  # in eq. 13 is -4.76955; in aerobulk is like this
                             0.78614)
+        if (meth == 'WMO2018'):
+            """WMO 2018 edition. Annex 4.B, eq. 4.B.1, 4.B.2, 4.B.5 """
+            Psat = 6.112*np.exp(17.62*temp/(243.12+temp))*(1.0016+3.15e-6*P -
+                                                           0.074/P)
         if (meth == 'Sonntag'):
             """Source: Sonntag, D., Advancements in the field of hygrometry,
             Meteorol. Z., N. F., 3, 51-66, 1994."""
@@ -413,3 +398,38 @@ def get_hum(hum, T, sst, P, qmeth):
         qsea = qsat_sea(sst, P, qmeth)/1000    # surface water q (g/kg)
     return qair, qsea
 #------------------------------------------------------------------------------
+
+
+def gamma_moist(sst, t, q):
+    """
+    Computes the moist adiabatic lapse-rate
+
+    Parameters
+    ----------
+    sst : float
+        sea surface temperature [K]
+    t : float
+        air temperature [K]
+    q : float
+        specific humidity of air [kg/kg]
+
+    Returns
+    -------
+    gamma : float
+        lapse rate [K/m]
+
+    """
+    if (np.nanmin(sst) < 200):  # if sst in Celsius convert to Kelvin
+        sst = sst+CtoK
+    if (np.nanmin(t) < 200):  # if sst in Celsius convert to Kelvin
+        t = t+CtoK
+
+    t = np.maximum(t, 180)
+    q = np.maximum(q,  1e-6)
+    w = q/(1-q) # mixing ratio w = q/(1-q)
+    iRT = 1/(287.05*t)
+    # latent heat of vaporization of water as a function of temperature
+    lv = (2.501-0.00237*(sst-CtoK))*1e6
+    gamma = 9.8*(1+lv*w*iRT)/(1005+np.power(lv, 2)*w*(287.05/461.495)*iRT/t)
+    return gamma
+#------------------------------------------------------------------------------