changed get_L to always output Rib following ecmwf's eq. 3.24.

changed flag "l" in AirSeaFluxCode to be set when Rib<-0.5 or Rib>0.2

AirSeaFluxCode.py

@@ -122,13 +122,13 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
                        34. downward longwave radiation (Rl)
                        35. downward shortwave radiation (Rs)
                        36. downward net longwave radiation (Rnl)
-                       37. flag ("n": normal, "o": out of nominal range,
+                       37. gust wind speed (ug)
+                       38. Bulk Richardson number (Rib)
+                       39. relative humidity (rh)
+                       40. flag ("n": normal, "o": out of nominal range,
                                  "u": u10n<0, "q":q10n<0
-                                 "m": missing, "l": z/L<0.01,
+                                 "m": missing, "l": Rib<-0.5 or Rib>0.2,
                                  "i": convergence fail at n)
-                       38. gust wind speed (ug)
-                       39. Bulk Richardson number (Rib)
-                       40. relative humidity (rh)
 
     2021 / Author S. Biri
     """
@@ -349,7 +349,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
                                   np.power(get_gust(gust[1], tv[ind], usr[ind],
                                   tsrv[ind], gust[2], lat[ind]), 2)))
                                   # Zeng et al. 1998 (20)
-        elif (gust[0] == 1 and (meth == "C30" or meth == "C35")):
+        elif (gust[0] == 1 and (meth == "C30" or meth == "C35")): # or meth == "C40"
             wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
                                 np.power(get_gust(gust[1], Ta[ind], usr[ind],
                                 tsrv[ind], gust[2], lat[ind]), 2))
@@ -427,8 +427,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
     flag = np.where((q10n < 0) & (flag == "n"), "q",
                     np.where((q10n < 0) & (flag != "n"), flag+[","]+["q"],
                              flag))
-    flag = np.where((np.abs(hin[0]/monob) < 0.01) & (flag == "n"), "l",
-                    np.where((np.abs(hin[0]/monob) < 0.01) & (flag != "n"),
+    flag = np.where((Rb < -0.5 | Rb > 0.2) & (flag == "n"), "l",
+                    np.where((Rb < -0.5 | Rb > 0.2) & (flag != "n"),
                              flag+[","]+["l"], flag))
     flag = np.where((itera == -1) & (flag == "n"), "i",
                     np.where((itera == -1) & (flag != "n"), flag+[","]+["i"],

flux_subs.py

@@ -947,25 +947,23 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
     """
     g = gc(lat)
     Rb = np.empty(sst.shape)
+    # as in NCAR, LY04
+    tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
+    # from eq. 3.24 ifs Cy46r1 pp. 37
+    thvs = sst*(1+0.6077*qsea) # virtual SST
+    dthv = (Ta-sst)*(1+0.6077*qair)+0.6077*Ta*(qair-qsea)
+    tv = 0.5*(thvs+Ta*(1+0.6077*qair)) # estimate tv within surface layer
+    # adjust wind to T sensor's height
+    uz = (wind-usr/kappa*(np.log(hin[0]/hin[1])-psim +
+                            psim_calc(hin[1]/monob, meth)))
+    Rb = g*dthv*hin[1]/(tv*uz*uz)
     if (L == "S80"):
-        # as in NCAR, LY04
-        tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
         temp = (g*kappa*tsrv /
                 np.maximum(np.power(usr, 2)*Ta*(1+0.6077*qair), 1e-9))
         temp = np.minimum(np.abs(temp), 200)*np.sign(temp)
         temp = np.minimum(np.abs(temp), 10/hin[0])*np.sign(temp)
         monob = 1/np.copy(temp)
     elif (L == "ecmwf"):
-        Rb = np.empty(sst.shape)
-        tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
-        # from eq. 3.24 ifs Cy46r1 pp. 37
-        thvs = sst*(1+0.6077*qsea) # virtual SST
-        dthv = (Ta-sst)*(1+0.6077*qair)+0.6077*Ta*(qair-qsea)
-        tv = 0.5*(thvs+Ta*(1+0.6077*qair)) # estimate tv within surface layer
-        # adjust wind to T sensor's height
-        uz = (wind-usr/kappa*(np.log(hin[0]/hin[1])-psim +
-                                psim_calc(hin[1]/monob, meth)))
-        Rb = g*dthv*hin[1]/(tv*uz*uz)
         zo = (0.11*visc_air(Ta)/usr+0.018*np.power(usr, 2)/g)
         zot = 0.40*visc_air(Ta)/usr
         zol = (Rb*(np.power(np.log((hin[1]+zo)/zo)-psim_calc((hin[1]+zo) /