Update flux_subs.py

flux_subs.py

@@ -978,8 +978,9 @@ def get_L(L, lat, usr, tsr, qsr, t10n, hin, Ta, sst, qair, qsea, q10n,
         monob = 1/np.copy(temp)
     elif (L == "ecmwf"):
         tsrv = tsr*(1+0.6077*qair)+0.6077*Ta*qsr
-        Rb = ((g*hin[1])/np.power(wind, 2))*(2*(Ta-sst)/(Ta+sst-g*hin[1]) +
-                                             0.6077*(qair-qsea))
+        Rb = g*hin[1]/(wind*wind)*((Ta-sst)/(Ta-0.5*(Ta-sst+g*hin[1] /
+                                                     (1005+1860*qair))) +
+                                   0.6077*(qair-qsea))
         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) /