removed strict constraint in get_L

renamed L options "S80" -> "tsrv" and "ecmwf" -> "Rb"
changed default L option relative to method used in get_init

AirSeaFluxCode.py

@@ -80,9 +80,10 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
             set 1 to keep points
         L : str
            Monin-Obukhov length definition options
-           "S80"  : default for S80, S88, LP82, YT96 and LY04
-           "ecmwf" : following ecmwf (IFS Documentation cy46r1), default for
-           ecmwf
+           "tsrv"  : default for "S80", "S88", "LP82", "YT96", "UA", "LY04",
+                     "C30", "C35"
+           "Rb" : following ecmwf (IFS Documentation cy46r1), default for
+                  "ecmwf", "Beljaars"
     Returns
     -------
         res : array that contains

flux_subs.py

@@ -935,6 +935,7 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, zo,
     Rb = np.empty(sst.shape)
     # as in aerobulk One_on_L in mod_phymbl.f90
     tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
+    # tsrv = tsr+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)
@@ -943,13 +944,12 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, zo,
     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"):
+    if (L == "tsrv"):
         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"):
+    elif (L == "Rb"):
         zol = (Rb*(np.power(np.log((hin[1]+zo)/zo)-psim_calc((hin[1]+zo) /
                                                               monob, meth) +
                             psim_calc(zo/monob, meth), 2) /

get_init.py

@@ -146,10 +146,12 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
         gust = [0, 0, 0]
     elif (np.size(gust) < 3):
         sys.exit("gust input must be a 3x1 array")
-    if (L not in [None, "S80", "ecmwf"]):
-        sys.exit("L input must be either None, S80 or ecmwf")
-    if (L == None):
-        L = "ecmwf"
+    if (L not in [None, "tsrv", "Rb"]):
+        sys.exit("L input must be either None, tsrv or Rb")
+    if ((L == None) and ((meth != "ecmwf") and (meth != "Beljaars"))):
+        L = "tsrv"
+    elif ((L == None) and ((meth == "ecmwf") or (meth == "Beljaars"))):
+        L = "Rb"
     if (tol == None):
         tol = ['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1]
     elif (tol[0] not in ['flux', 'ref', 'all']):