fixed issue #33

AirSeaFluxCode.py

@@ -192,9 +192,15 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
     flag = np.empty(spd.shape, dtype="object")
     flag[:] = "n"
     if (hum[0] == 'no'):
-        flag = np.where(np.isnan(spd+T+SST+P+Rs+Rl), "m", flag)
+        if (cskin == 1):
+            flag = np.where(np.isnan(spd+T+SST+P+Rs+Rl), "m", flag)
+        else:
+            flag = np.where(np.isnan(spd+T+SST+P), "m", flag)
     else:
-        flag = np.where(np.isnan(spd+T+SST+hum[1]+P+Rs+Rl), "m", flag)
+        if (cskin == 1):
+            flag = np.where(np.isnan(spd+T+SST+hum[1]+P+Rs+Rl), "m", flag)
+        else:
+            flag = np.where(np.isnan(spd+T+SST+hum[1]+P), "m", flag)
         flag = np.where(rh > 100, "r", flag)
 
     dt = Ta - sst
@@ -391,7 +397,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")): # or meth == "C40"
+        elif (gust[0] == 1 and (meth == "C30" or meth == "C35")):
             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))
@@ -526,11 +532,6 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
                                   (np.char.find(flag.astype(str), 'u') == -1) &
                                   (np.char.find(flag.astype(str), 'q') == -1)),
                                  flag+[","]+["o"], flag))
-    # Do not output radiation parameters if they are not input, unless the
-    # cs/wl is switched on
-    if (((cskin == 0) and (wl == 0)) and
-        (np.all(Rl == 370) and np.all(Rs == 150))):
-        Rl, Rs, Rnl = Rl*np.nan, Rs*np.nan, Rnl*np.nan
     # Do not calculate lhf if a measure of humidity is not input
     if (hum[0] == 'no'):
         latent = np.ones(sst.shape)*np.nan

get_init.py

@@ -115,10 +115,6 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
         P = np.ones(spd.shape)*1013
     elif (((np.all(P != None)) or np.all(~np.isnan(P))) and np.size(P) == 1):
         P = np.ones(spd.shape)*np.copy(P)
-    if (np.all(Rl == None) or np.all(np.isnan(Rl))):
-        Rl = np.ones(spd.shape)*370    # set to default for COARE3.5
-    if (np.all(Rs == None) or np.all(np.isnan(Rs))):
-        Rs = np.ones(spd.shape)*150  # set to default for COARE3.5
     if ((cskin == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
                              or meth == "YT96" or meth == "UA" or
                              meth == "LY04")):
@@ -134,6 +130,15 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
             skin = "Beljaars"
     if (wl == None):
         wl = 0
+    if ((cskin == 1 or wl == 1) and (np.all(Rl == None) or np.all(np.isnan(Rl)))
+        and ((np.all(Rs == None) or np.all(np.isnan(Rs))))):
+        sys.exit("Cool skin/warm layer is switched ON; Radiation input should not be"
+                 " empty")
+    elif (cskin == 0):
+        if (np.all(Rl == None)):
+            Rl = np.ones(spd.shape)*np.nan
+        if (np.all(Rs == None)):
+            Rs = np.ones(spd.shape)*np.nan
     if (np.all(gust == None) and (meth == "C30" or meth == "C35")):
         gust = [1, 1.2, 600]
     elif (np.all(gust == None) and (meth == "UA" or meth == "ecmwf" or