Update flux_subs.py, AirSeaFluxCode.py, get_init.py, toy_ASFC.py, hum_subs.py files

Deleted Documentation.pdf, data_all_out.csv, data_all_stats.txt, readme.txt files

AirSeaFluxCode.py

@@ -122,9 +122,9 @@ 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, "ul": spd<2m/s,
+                       37. flag ("n": normal, "o": out of nominal range,
                                  "u": u10n<0, "q":q10n<0
-                                 "t": DT>10, "l": z/L<0.01,
+                                 "m": missing, "l": z/L<0.01,
                                  "i": convergence fail at n)
 
     2021 / Author S. Biri
@@ -138,7 +138,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
                                                               skin, wl, gust, L,
                                                               tol, meth, qmeth)
     flag = np.ones(spd.shape, dtype="object")*"n"
-    flag = np.where(spd < 2, "ul", flag)
+    flag = np.where(np.isnan(spd+T+SST+lat+hum[1]+P+Rs) & (flag == "n"),
+                    "m", np.where(np.isnan(spd+T+SST+lat+hum[1]+P+Rs) &
+                                  (flag != "n"), flag+[","]+["m"], flag))
     ref_ht = 10        # reference height
     h_in = get_heights(hin, len(spd))  # heights of input measurements/fields
     h_out = get_heights(hout, 1)       # desired height of output variables
@@ -152,6 +154,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
                   np.copy(T)*np.power(1000/P,287.1/1004.67))  # potential T
     sst = np.where(SST < 200, np.copy(SST)+CtoK, np.copy(SST))
     qair, qsea = get_hum(hum, T, sst, P, qmeth)
+    Rb = np.empty(sst.shape)
     #lapse rate
     tlapse = gamma_moist(SST, T, qair/1000)
     Ta = np.where(T < 200, np.copy(T)+CtoK+tlapse*h_in[1],
@@ -163,26 +166,24 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
 
     dt = Ta - sst
     dq = qair - qsea
-    flag = np.where((dt > 10) & (flag == "n"), "t",
-                    np.where((dt > 10) & (flag != "n"), flag+[","]+["t"],
-                             flag))
+
     #  first guesses
     t10n, q10n = np.copy(Ta), np.copy(qair)
-    tv10n = t10n*(1+0.61*q10n)
+    tv10n = t10n*(1+0.6077*q10n)
     #  Zeng et al. 1998
-    tv=th*(1.+0.61*qair)   # virtual potential T
-    dtv=dt*(1.+0.61*qair)+0.61*th*dq
+    tv=th*(1+0.6077*qair)   # virtual potential T
+    dtv=dt*(1+0.6077*qair)+0.6077*th*dq
     # ------------
     rho = P*100/(287.1*tv10n)
     lv = (2.501-0.00237*(sst-CtoK))*1e6
     cp = 1004.67*(1 + 0.00084*qsea)
     u10n = np.copy(spd)
-    cdn = cdn_calc(u10n, Ta, None, lat, meth)
-    ctn, ct, cqn, cq = (np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan,
+    cd10n = cdn_calc(u10n, Ta, None, lat, meth)
+    ct10n, ct, cq10n, cq = (np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan,
                         np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan)
     psim, psit, psiq = (np.zeros(spd.shape), np.zeros(spd.shape),
                         np.zeros(spd.shape))
-    cd = cd_calc(cdn, h_in[0], ref_ht, psim)
+    cd = cd_calc(cd10n, h_in[0], ref_ht, psim)
     tsr, tsrv = np.zeros(spd.shape), np.zeros(spd.shape)
     qsr = np.zeros(spd.shape)
     # cskin parameters
@@ -228,22 +229,23 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
         elif (tol[0] == 'all'):
             old = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
                             np.copy(tau), np.copy(sensible), np.copy(latent)])
-        cdn[ind] = cdn_calc(u10n[ind], Ta[ind], None, lat[ind], meth)
-        if (np.all(np.isnan(cdn))):
+        cd10n[ind] = cdn_calc(u10n[ind], Ta[ind], None, lat[ind], meth)
+        if (np.all(np.isnan(cd10n))):
             break
-            logging.info('break %s at iteration %s cdn<0', meth, it)
-        zo[ind] = ref_ht/np.exp(kappa/np.sqrt(cdn[ind]))
+            logging.info('break %s at iteration %s cd10n<0', meth, it)
+        zo[ind] = ref_ht/np.exp(kappa/np.sqrt(cd10n[ind]))
         psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth)
-        cd[ind] = cd_calc(cdn[ind], h_in[0, ind], ref_ht, psim[ind])
-        ctn[ind], cqn[ind] = ctcqn_calc(h_in[1, ind]/monob[ind], cdn[ind],
-                                        u10n[ind], zo[ind], Ta[ind], meth)
+        cd[ind] = cd_calc(cd10n[ind], h_in[0, ind], ref_ht, psim[ind])
+        ct10n[ind], cq10n[ind] = ctcqn_calc(h_in[1, ind]/monob[ind],
+                                            cd10n[ind], u10n[ind], zo[ind],
+                                            Ta[ind], meth)
         zot[ind] = ref_ht/(np.exp(np.power(kappa, 2) /
-                           (ctn[ind]*np.log(ref_ht/zo[ind]))))
+                           (ct10n[ind]*np.log(ref_ht/zo[ind]))))
         zoq[ind] = ref_ht/(np.exp(np.power(kappa, 2) /
-                           (cqn[ind]*np.log(ref_ht/zo[ind]))))
+                           (cq10n[ind]*np.log(ref_ht/zo[ind]))))
         psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
         psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
-        ct[ind], cq[ind] = ctcq_calc(cdn[ind], cd[ind], ctn[ind], cqn[ind],
+        ct[ind], cq[ind] = ctcq_calc(cd10n[ind], cd[ind], ct10n[ind], cq10n[ind],
                                       h_in[1, ind], h_in[2, ind], ref_ht,
                                       psit[ind], psiq[ind])
         usr[ind], tsr[ind], qsr[ind] = get_strs(h_in[:, ind], monob[ind],
@@ -325,12 +327,14 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
                      tsr[ind]/kappa*(np.log(h_in[1, ind]/ref_ht)-psit[ind]))
         q10n[ind] = (qair[ind] -
                      qsr[ind]/kappa*(np.log(h_in[2, ind]/ref_ht)-psiq[ind]))
-        tv10n[ind] = t10n[ind]*(1+0.61*q10n[ind])
-        tsrv[ind], monob[ind] = get_L(L, lat[ind], usr[ind], tsr[ind],
-                                      qsr[ind], t10n[ind], h_in[:, ind],
-                                      Ta[ind], sst[ind],
-                                      qair[ind], qsea[ind], q10n[ind],
-                                      wind[ind], np.copy(monob[ind]), meth)
+        tv10n[ind] = t10n[ind]*(1+0.6077*q10n[ind])
+        tsrv[ind], monob[ind], Rb[ind] = get_L(L, lat[ind], usr[ind], tsr[ind],
+                                               qsr[ind], h_in[:, ind], Ta[ind],
+                                               sst[ind]-dter[ind]*cskin+dtwl[ind]*wl,
+                                               qair[ind], qsea[ind], wind[ind],
+                                               np.copy(monob[ind]), psim[ind],
+                                               meth)
+        # sst[ind]-dter[ind]*cskin+dtwl[ind]*wl
         psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth)
         psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
         psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
@@ -395,14 +399,24 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
     # calculate output parameters
     rho = (0.34838*P)/(tv10n)
     t10n = t10n-(273.16+tlapse*ref_ht)
-    zo = ref_ht/np.exp(kappa/cdn**0.5)
-    zot = ref_ht/(np.exp(kappa**2/(ctn*np.log(ref_ht/zo))))
-    zoq = ref_ht/(np.exp(kappa**2/(cqn*np.log(ref_ht/zo))))
+    # solve for zo from cd10n
+    zo = ref_ht/np.exp(kappa/np.sqrt(cd10n))
+    # adjust neutral cdn at any output height
+    cdn = np.power(kappa/np.log(hout/zo), 2)
+    cd = cd_calc(cdn, h_in[0], h_out[0], psim)
+     # solve for zot, zoq from ct10n, cq10n
+    zot = ref_ht/(np.exp(kappa**2/(ct10n*np.log(ref_ht/zo))))
+    zoq = ref_ht/(np.exp(kappa**2/(cq10n*np.log(ref_ht/zo))))
+    # adjust neutral ctn, cqn at any output height
+    ctn =np.power(kappa,2)/(np.log(hout/zo)*np.log(hout/zot))
+    cqn =np.power(kappa,2)/(np.log(hout/zo)*np.log(hout/zoq))
+    ct, cq = ctcq_calc(cdn, cd, ctn, cqn, h_in[1], h_in[2], h_out[1],
+                       psit, psiq)
     uref = (spd-usr/kappa*(np.log(h_in[0]/h_out[0])-psim +
             psim_calc(h_out[0]/monob, meth)))
     tref = (Ta-tsr/kappa*(np.log(h_in[1]/h_out[1])-psit +
             psit_calc(h_out[0]/monob, meth)))
-    tref = tref-(273.16+tlapse*h_out[1])
+    tref = tref-(CtoK+tlapse*h_out[1])
     qref = (qair-qsr/kappa*(np.log(h_in[2]/h_out[2]) -
             psit+psit_calc(h_out[2]/monob, meth)))
     flag = np.where((q10n < 0) & (flag == "n"), "q",
@@ -414,7 +428,41 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
     flag = np.where((itera == -1) & (flag == "n"), "i",
                     np.where((itera == -1) & (flag != "n"), flag+[","]+["i"],
                              flag))
-    res = np.zeros((36, len(spd)))
+    if (meth == "S80"):
+        flag = np.where(((u10n < 6) | (u10n > 22)) & (flag == "n"), "o",
+                        np.where(((u10n < 6) | (u10n > 22)) & (flag != "n"),
+                                 flag+[","]+["o"], flag))
+    elif (meth == "LP82"):
+        flag = np.where(((u10n < 3) | (u10n > 25)) & (flag == "n"), "o",
+                        np.where(((u10n < 3) | (u10n > 25)) & (flag != "n"),
+                                 flag+[","]+["o"], flag))
+    elif (meth == "YT96"):
+        flag = np.where(((u10n < 3) | (u10n > 26)) & (flag == "n"), "o",
+                        np.where(((u10n < 3) | (u10n > 26)) & (flag != "n"),
+                                 flag+[","]+["o"], flag))
+    elif (meth == "UA"):
+        flag = np.where(((u10n < 0.5) | (u10n > 18)) & (flag == "n"), "o",
+                        np.where(((u10n < 0.5) | (u10n > 18)) & (flag != "n"),
+                                 flag+[","]+["o"], flag))
+    elif (meth == "LY04"):
+        flag = np.where((u10n < 0.5) & (flag == "n"), "o",
+                        np.where((u10n < 0.5) & (flag != "n"),
+                                 flag+[","]+["o"], flag))
+    if (hum == None):
+        rh = np.ones(sst.shape)*80
+    elif (hum[0] == 'rh'):
+        rh = hum[1]
+        rh = np.where(rh > 100, np.nan, rh)
+    elif (hum[0] == 'Td'):
+        Td = hum[1] # dew point temperature (K)
+        Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td))
+        T = np.where(T < 200, np.copy(T)+CtoK, np.copy(T))
+        esd = 611.21*np.exp(17.502*((Td-273.16)/(Td-32.19)))
+        es = 611.21*np.exp(17.502*((T-273.16)/(T-32.19)))
+        rh = 100*esd/es
+        rh = np.where(rh > 100, np.nan, rh)
+
+    res = np.zeros((39, len(spd)))
     res[0][:] = tau
     res[1][:] = sensible
     res[2][:] = latent
@@ -451,21 +499,24 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
     res[33][:] = Rl
     res[34][:] = Rs
     res[35][:] = Rnl
+    res[36][:] = np.sqrt(np.power(wind, 2)-np.power(spd, 2))
+    res[37][:] = Rb
+    res[38][:] = rh
 
     if (out == 0):
         res[:, ind] = np.nan
     # set missing values where data have non acceptable values
-    res = [np.where(spd < 0, np.nan, res[i][:]) for i in range(36)]
-    res = [np.where(q10n < 0, np.nan, res[i][:]) for i in range(36)]
-    res = np.asarray(res)
+    res = np.asarray([np.where((spd < 0) | (q10n < 0), np.nan,
+                               res[i][:]) for i in range(39)])
     # output with pandas
     resAll = pd.DataFrame(data=res.T, index=range(len(spd)),
                         columns=["tau", "shf", "lhf", "L", "cd", "cdn", "ct",
-                                "ctn", "cq", "cqn", "tsrv", "tsr", "qsr",
-                                "usr", "psim", "psit","psiq", "u10n", "t10n",
-                                "tv10n", "q10n", "zo", "zot", "zoq", "uref",
-                                "tref", "qref", "iteration", "dter", "dqer",
-                                "dtwl", "qair", "qsea", "Rl", "Rs", "Rnl"])
+                                 "ctn", "cq", "cqn", "tsrv", "tsr", "qsr",
+                                 "usr", "psim", "psit","psiq", "u10n", "t10n",
+                                 "tv10n", "q10n", "zo", "zot", "zoq", "uref",
+                                 "tref", "qref", "iteration", "dter", "dqer",
+                                 "dtwl", "qair", "qsea", "Rl", "Rs", "Rnl",
+                                 "ug", "Rib", "rh"])
     resAll["flag"] = flag
     return resAll
     

data_all_stats.txt

-Input summary
-input file name: data_all.csv, 
- method: UA, 
- gustiness: [1, 1, 1000], 
- cskin: 0, 
- tolerance: ['flux', 0.001, 0.1, 0.1]
-| var   | mean      | median    | min       | max      | 5%        | 95%      |
-|-------|-----------|-----------|-----------|----------|-----------|----------|
-| tau   | 7.24e-02  | 4.95e-02  | 4.96e-04  | 5.64e-01 | 7.33e-03  | 2.12e-01 |
-| shf   | 6.80e+00  | 6.16e+00  | -4.06e+01 | 5.56e+01 | -9.07e+00 | 2.52e+01 |
-| lhf   | 8.26e+01  | 6.71e+01  | -4.15e+01 | 5.46e+02 | -2.89e-01 | 2.13e+02 |
-| L     | 4.11e+02  | -4.40e+01 | -1.62e+05 | 1.20e+06 | -5.15e+02 | 3.19e+02 |
-| cd    | 1.15e-03  | 1.16e-03  | 3.44e-04  | 2.02e-03 | 8.74e-04  | 1.37e-03 |
-| cdn   | 1.15e-03  | 1.12e-03  | 9.70e-04  | 1.71e-03 | 9.86e-04  | 1.41e-03 |
-| ct    | 1.18e-03  | 1.18e-03  | 3.60e-04  | 2.87e-03 | 8.90e-04  | 1.46e-03 |
-| ctn   | 1.12e-03  | 1.12e-03  | 1.04e-03  | 1.25e-03 | 1.05e-03  | 1.21e-03 |
-| cq    | 1.18e-03  | 1.18e-03  | 3.60e-04  | 2.87e-03 | 8.90e-04  | 1.46e-03 |
-| cqn   | 1.12e-03  | 1.12e-03  | 1.04e-03  | 1.25e-03 | 1.05e-03  | 1.21e-03 |
-| tsrv  | -5.18e-02 | -5.57e-02 | -2.34e-01 | 8.52e-02 | -1.27e-01 | 3.28e-02 |
-| tsr   | -2.70e-02 | -2.80e-02 | -1.67e-01 | 8.30e-02 | -8.42e-02 | 3.53e-02 |
-| qsr   | -1.37e-04 | -1.24e-04 | -5.57e-04 | 4.67e-05 | -3.24e-04 | 6.52e-07 |
-| usr   | 2.20e-01  | 2.04e-01  | 2.02e-02  | 6.92e-01 | 7.95e-02  | 4.20e-01 |
-| psim  | 4.77e-01  | 4.56e-01  | -8.86e+00 | 4.10e+00 | -7.94e-01 | 1.75e+00 |
-| psit  | 9.03e-01  | 8.44e-01  | -8.86e+00 | 5.74e+00 | -8.02e-01 | 2.75e+00 |
-| psiq  | 9.03e-01  | 8.44e-01  | -8.86e+00 | 5.74e+00 | -8.02e-01 | 2.75e+00 |
-| u10n  | 6.37e+00  | 6.10e+00  | 1.17e+00  | 1.67e+01 | 2.56e+00  | 1.12e+01 |
-| t10n  | 1.79e+01  | 1.82e+01  | -2.76e+00 | 2.98e+01 | 2.82e+00  | 2.79e+01 |
-| tv10n | 2.93e+02  | 2.93e+02  | 2.71e+02  | 3.08e+02 | 2.77e+02  | 3.04e+02 |
-| q10n  | 1.06e-02  | 9.67e-03  | -1.49e-03 | 2.52e-02 | 3.58e-03  | 1.78e-02 |
-| zo    | 8.96e-05  | 6.30e-05  | 2.65e-05  | 6.37e-04 | 2.94e-05  | 2.39e-04 |
-| zot   | 6.01e-05  | 6.18e-05  | 1.63e-05  | 1.23e-04 | 4.09e-05  | 7.13e-05 |
-| zoq   | 6.01e-05  | 6.18e-05  | 1.63e-05  | 1.23e-04 | 4.09e-05  | 7.13e-05 |
-| urefs | 6.18e+00  | 5.93e+00  | 6.25e-01  | 1.65e+01 | 2.30e+00  | 1.11e+01 |
-| trefs | 1.80e+01  | 1.84e+01  | -2.75e+00 | 2.98e+01 | 2.87e+00  | 2.80e+01 |
-| qrefs | 1.10e-02  | 1.00e-02  | 7.66e-04  | 2.53e-02 | 3.66e-03  | 1.83e-02 |
-| itera | 3.41e+00  | 3.00e+00  | 2.00e+00  | 6.00e+00 | 3.00e+00  | 4.00e+00 |
-| dter  | 0.00e+00  | 0.00e+00  | 0.00e+00  | 0.00e+00 | 0.00e+00  | 0.00e+00 |
-| dqer  | 0.00e+00  | 0.00e+00  | 0.00e+00  | 0.00e+00 | 0.00e+00  | 0.00e+00 |
-| dtwl  | 3.00e-01  | 3.00e-01  | 3.00e-01  | 3.00e-01 | 3.00e-01  | 3.00e-01 |
-| qair  | 1.09e-02  | 9.90e-03  | 7.43e-04  | 2.53e-02 | 3.64e-03  | 1.82e-02 |
-| qsea  | 1.47e-02  | 1.39e-02  | 3.26e-03  | 2.61e-02 | 4.76e-03  | 2.46e-02 |
-| Rl    | 3.70e+02  | 3.70e+02  | 3.70e+02  | 3.70e+02 | 3.70e+02  | 3.70e+02 |
-| Rs    | 2.61e+02  | 2.33e+02  | 1.67e-02  | 9.71e+02 | 2.84e+01  | 5.37e+02 |
-| Rnl   | 4.25e+01  | 4.50e+01  | -6.02e+01 | 1.06e+02 | -3.65e+01 | 9.96e+01 |
--------------------------------------------------------------------------------
-

flux_subs.py

@@ -6,7 +6,7 @@ from util_subs import (CtoK, kappa, gc, visc_air)
 
 def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
     """
-    Calculates neutral drag coefficient
+    Calculates 10m neutral drag coefficient
 
     Parameters
     ----------
@@ -56,7 +56,7 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
 
 def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
     """
-    Calculates neutral drag coefficient from roughness length
+    Calculates 10m neutral drag coefficient from roughness length
 
     Parameters
     ----------
@@ -97,9 +97,6 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
             zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
         elif (meth == "C35"):
             a = 0.011*np.ones(Ta.shape)
-            # a = np.where(u10n > 19, 0.0017*19-0.0050,
-            #             np.where((u10n > 7) & (u10n <= 18),
-            #                       0.0017*u10n-0.0050, a))
             a = np.where(u10n > 19, 0.0017*19-0.0050, 0.0017*u10n-0.0050)
             zo = 0.11*visc_air(Ta)/usr+a*np.power(usr, 2)/g
         elif (meth == "C40"):
@@ -117,7 +114,7 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
 # ---------------------------------------------------------------------
 
 
-def cd_calc(cdn, height, ref_ht, psim):
+def cd_calc(cdn, hin, hout, psim):
     """
     Calculates drag coefficient at reference height
 
@@ -125,9 +122,9 @@ def cd_calc(cdn, height, ref_ht, psim):
     ----------
     cdn : float
         neutral drag coefficient
-    height : float
+    hin : float
         original sensor height  [m]
-    ref_ht : float
+    hout : float
         reference height        [m]
     psim : float
         momentum stability function
@@ -136,14 +133,14 @@ def cd_calc(cdn, height, ref_ht, psim):
     -------
     cd : float
     """
-    cd = (cdn/np.power(1+(np.sqrt(cdn)*(np.log(height/ref_ht)-psim))/kappa, 2))
+    cd = (cdn/np.power(1+(np.sqrt(cdn)*(np.log(hin/hout)-psim))/kappa, 2))
     return cd
 # ---------------------------------------------------------------------
 
 
 def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
     """
-    Calculates neutral heat and moisture exchange coefficients
+    Calculates 10m neutral heat and moisture exchange coefficients
 
     Parameters
     ----------
@@ -209,9 +206,6 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
                        1.0e-4/np.power(rr, 0.55)) # temperature roughness
         zoq = np.where(2.0e-5/np.power(rr,0.22) > 1.1e-4/np.power(rr,0.9),
                        1.1e-4/np.power(rr,0.9), 2.0e-5/np.power(rr,0.22))
-        # moisture roughness determined by the CLIMODE, GASEX and CBLAST data
-#        zoq = np.where(5e-5/np.power(rr, 0.6) > 1.15e-4, 1.15e-4,
-#                       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 == "ecmwf" or meth == "Beljaars"):
@@ -227,7 +221,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
 # ---------------------------------------------------------------------
 
 
-def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, ref_ht, psit, psiq):
+def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, hout, psit, psiq):
     """
     Calculates heat and moisture exchange coefficients at reference height
 
@@ -245,8 +239,8 @@ def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, ref_ht, psit, psiq):
         original temperature sensor height [m]
     hq : float
         original moisture sensor height    [m]
-    ref_ht : float
-        reference height                   [m]
+    hout : float
+        output height                   [m]
     psit : float
         heat stability function
     psiq : float
@@ -260,9 +254,9 @@ def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, ref_ht, psit, psiq):
        moisture exchange coefficient
     """
     ct = (ctn*np.sqrt(cd/cdn) /
-          (1+ctn*((np.log(ht/ref_ht)-psit)/(kappa*np.sqrt(cdn)))))
+          (1+ctn*((np.log(ht/hout)-psit)/(kappa*np.sqrt(cdn)))))
     cq = (cqn*np.sqrt(cd/cdn) /
-          (1+cqn*((np.log(hq/ref_ht)-psiq)/(kappa*np.sqrt(cdn)))))
+          (1+cqn*((np.log(hq/hout)-psiq)/(kappa*np.sqrt(cdn)))))
     return ct, cq
 # ---------------------------------------------------------------------
 
@@ -736,8 +730,6 @@ def cs_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, lat):
         # # 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))
-        # fs = np.maximum(0.065+11*delta-(6.6e-5/delta)*(1-np.exp(-delta/8e-4)),
-        #                 0.01)
         Q = Qnsol-fs*Rns
         d = delta(aw, Q, usr, lat)
     dtc = Q*d/0.6  # (rhow*cw*kw)eq. 8.151
@@ -918,8 +910,8 @@ def get_gust(beta, Ta, usr, tsrv, zi, lat):
 # ---------------------------------------------------------------------
 
 
-def get_L(L, lat, usr, tsr, qsr, t10n, hin, Ta, sst, qair, qsea, q10n,
-          wind, monob, meth):
+def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
+          meth):
     """
     calculates Monin-Obukhov length and virtual star temperature
 
@@ -968,6 +960,7 @@ def get_L(L, lat, usr, tsr, qsr, t10n, hin, Ta, sst, qair, qsea, q10n,
 
     """
     g = gc(lat)
+    Rb = np.empty(sst.shape)
     if (L == "S80"):
         # as in NCAR, LY04
         tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
@@ -977,20 +970,26 @@ def get_L(L, lat, usr, tsr, qsr, t10n, hin, Ta, sst, qair, qsea, q10n,
         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
-        Rb = g*hin[1]/(wind*wind)*((Ta-sst)/(Ta-0.5*(Ta-sst+g*hin[1] /
-                                                     (1005+1860*qair))) +
-                                   0.6077*(qair-qsea))
+        # 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[0]+zo)/zo)-psim_calc((hin[0]+zo) /
+        zol = (Rb*(np.power(np.log((hin[1]+zo)/zo)-psim_calc((hin[1]+zo) /
                                                               monob, meth) +
                             psim_calc(zo/monob, meth), 2) /
-                   (np.log((hin[0]+zo)/zot) -
-                    psit_calc((hin[0]+zo)/monob, meth) +
+                   (np.log((hin[1]+zo)/zot) -
+                    psit_calc((hin[1]+zo)/monob, meth) +
                     psit_calc(zot/monob, meth))))
-        monob = hin[0]/zol
-    return tsrv, monob
+        monob = hin[1]/zol 
+    return tsrv, monob, Rb
 #------------------------------------------------------------------------------
 
 

get_init.py

@@ -137,7 +137,7 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, skin, wl, gust, L, tol, meth,
     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, 0, 1, 2 or 3")
+        sys.exit("L input must be either None, S80 or ecmwf")
     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"
@@ -149,4 +149,4 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, skin, wl, gust, L, tol, meth,
         tol = ['flux', 1e-3, 0.1, 0.1]
     elif (tol[0] not in ['flux', 'ref', 'all']):
         sys.exit("unknown tolerance input")
-    return lat, P, Rl, Rs, cskin, skin, wl, gust, tol, L
+    return lat, P, Rl, Rs, cskin, skin, wl, gust, tol, L
\ No newline at end of file

hum_subs.py

@@ -384,6 +384,7 @@ def get_hum(hum, T, sst, P, qmeth):
         if (np.all(RH < 1)):
             sys.exit("input relative humidity units should be \%")
             qair, qsea = np.nan, np.nan
+        RH = np.where(RH > 100, np.nan, RH)  # ensure RH <=100
         qsea = qsat_sea(sst, P, qmeth)/1000    # surface water q (kg/kg)
         qair = qsat_air(T, P, RH, qmeth)/1000  # q of air (kg/kg)
     elif (hum[0] == 'q'):
@@ -396,6 +397,7 @@ def get_hum(hum, T, sst, P, qmeth):
         esd = 611.21*np.exp(17.502*((Td-273.16)/(Td-32.19)))
         es = 611.21*np.exp(17.502*((T-273.16)/(T-32.19)))
         RH = 100*esd/es
+        RH = np.where(RH > 100, np.nan, RH) # ensure RH <=100
         qair = qsat_air(T, P, RH, qmeth)/1000  # q of air (kg/kg)
         qsea = qsat_sea(sst, P, qmeth)/1000    # surface water q (kg/kg)
     return qair, qsea

readme.txt

-input file name: data_all.csv, 
- method: UA, 
- gustiness: [1, 1, 1000], 
- cskin: 0, 
- tolerance: ['all', 0.01, 0.01, 1e-05, 0.001, 0.1, 0.1], 
- output is written in: data_all_out.csv