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

AirSeaFluxCode.py

@@ -304,7 +304,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
             ii = False
         else:
             ii = True
-    itera = np.where(itera > n, -1, itera)
+    itera[ind] = -1
     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

@@ -42,7 +42,8 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
                        (0.60 + 0.070*u10n)*0.001, (0.61+0.567/u10n)*0.001))
     elif (meth == "LY04"):
         cdn = np.where(u10n >= 0.5,
-                       (0.142+(2.7/u10n)+(u10n/13.09))*0.001, np.nan)
+                       (0.142+(2.7/u10n)+(u10n/13.09))*0.001,
+                       (0.142+(2.7/0.5)+(0.5/13.09))*0.001)
     else:
         print("unknown method cdn: "+meth)
     return cdn
@@ -86,7 +87,7 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
             zo = 0.013*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
         elif (meth == "C30"):
             a = 0.011*np.ones(Ta.shape)
-            a = np.where(u10n > 10, 0.011+(u10n-10)/(18-10)*(0.018-0.011),
+            a = np.where(u10n > 10, 0.011+(u10n-10)*(0.018-0.011)/(18-10),
                          np.where(u10n > 18, 0.018, a))
             zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
         elif (meth == "C35"):
@@ -688,8 +689,8 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
     monob : float
         Monin-Obukhov length from previous iteration step (m)
     meth : str
-        bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA",
-        "LY04", "C30", "C35", "C40", "ERA5"
+        bulk parameterisation method option: "S80", "S88", "LP82", "YT96",
+        "UA", "LY04", "C30", "C35", "C40", "ERA5"
 
     Returns
     -------
@@ -704,7 +705,7 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
         tsrv = tsr+0.61*t10n*qsr
         monob = ((tv10n*np.power(usr, 2))/(g*kappa*tsrv))
         monob = np.where(np.fabs(monob) < 1, np.where(monob < 0, -1, 1), monob)
-     elif (L == "ERA5"):
+    elif (L == "ERA5"):
         tsrv = tsr+0.61*t10n*qsr
         Rb = ((g*h_in[0]*((2*dt)/(Ta+sst-g*h_in[0])+0.61*dq)) /
               np.power(wind, 2))

get_init.py

@@ -76,6 +76,10 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
     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")
+    elif ((spd.dtype not in ['float64', 'float32']) or
+          (T.dtype not in ['float64', 'float32']) or
+          (SST.dtype not in ['float64', 'float32'])):
+        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"]:
@@ -111,14 +115,16 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
         gust = [1, 1, 1000]
     elif np.all(gust == None):
         gust = [1, 1.2, 800]
+    elif ((np.size(gust) < 3) and (gust == 0)):
+        gust = [0, 0, 0]
     elif (np.size(gust) < 3):
         sys.exit("gust input must be a 3x1 array")
     if (L not in [None, "S80", "ERA5"]):
         sys.exit("L input must be either None, 0, 1, 2 or 3")
     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 == "YT96" or meth == "LY04" or
+                         meth == "UA" or meth == "C30" or meth == "C35"
+                         or meth == "C40")):
         L = "S80"
     elif ((L == None) and (meth == "ERA5")):
         L = "ERA5"