C40 option removed

AirSeaFluxCode.py

@@ -58,7 +58,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
             default for UA, ecmwf [1, 1, 1000]
             default else [1, 1.2, 800]
         meth : str
-            "S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35", "C40",
+            "S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
             "ecmwf", "Beljaars"
         qmeth : str
             is the saturation evaporation method to use amongst
@@ -86,17 +86,17 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
     Returns
     -------
         res : array that contains
-                       1. momentum flux (N/m^2)
-                       2. sensible heat (W/m^2)
-                       3. latent heat (W/m^2)
-                       4. Monin-Obhukov length (mb)
+                       1. momentum flux       (N/m^2)
+                       2. sensible heat       (W/m^2)
+                       3. latent heat         (W/m^2)
+                       4. Monin-Obhukov length (m)
                        5. drag coefficient (cd)
                        6. neutral drag coefficient (cdn)
-                       7. heat exhange coefficient (ct)
-                       8. neutral heat exhange coefficient (ctn)
+                       7. heat exchange coefficient (ct)
+                       8. neutral heat exchange coefficient (ctn)
                        9. moisture exhange coefficient (cq)
-                       10. neutral moisture exhange coefficient (cqn)
-                       11. star virtual temperature (tsrv)
+                       10. neutral moisture exchange coefficient (cqn)
+                       11. star virtual temperatcure (tsrv)
                        12. star temperature (tsr)
                        13. star specific humidity (qsr)
                        14. star wind speed (usr)
@@ -110,7 +110,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
                        22. surface roughness length (zo)
                        23. heat roughness length (zot)
                        24. moisture roughness length (zoq)
-                       25. velocity at reference height (uref)
+                       25. wind speed at reference height (uref)
                        26. temperature at reference height (tref)
                        27. specific humidity at reference height (qref)
                        28. number of iterations until convergence
@@ -349,8 +349,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))

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 10m neutral drag coefficient
+    Calculates neutral drag coefficient
 
     Parameters
     ----------
@@ -33,7 +33,7 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
                        np.where((u10n < 11) & (u10n >= 4), 1.2*0.001,
                                 (0.49+0.065*u10n)*0.001))
     elif (meth == "S88" or meth == "UA" or meth == "ecmwf" or meth == "C30" or
-          meth == "C35" or meth == "C40" or meth == "Beljaars"):
+          meth == "C35" or meth == "Beljaars"): 
         cdn = cdn_from_roughness(u10n, Ta, None, lat, meth)
     elif (meth == "YT96"):
         # for u<3 YT96 convert usr in eq. 21 to cdn
@@ -56,7 +56,7 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
 
 def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
     """
-    Calculates 10m neutral drag coefficient from roughness length
+    Calculates neutral drag coefficient from roughness length
 
     Parameters
     ----------
@@ -99,10 +99,6 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
             a = 0.011*np.ones(Ta.shape)
             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"):
-            a = 0.011*np.ones(Ta.shape)
-            a = np.where(u10n > 22, 0.0016*22-0.0035, 0.0016*u10n-0.0035)
-            zo = a*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr # surface roughness
         elif ((meth == "ecmwf" or meth == "Beljaars")):
             # eq. (3.26) p.38 over sea IFS Documentation cy46r1
             zo = 0.018*np.power(usr, 2)/g+0.11*visc_air(Ta)/usr
@@ -140,7 +136,7 @@ def cd_calc(cdn, hin, hout, psim):
 
 def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
     """
-    Calculates 10m neutral heat and moisture exchange coefficients
+    Calculates neutral heat and moisture exchange coefficients
 
     Parameters
     ----------
@@ -198,16 +194,6 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
         zot=zoq  # temperature roughness
         cqn = kappa**2/np.log(10/zo)/np.log(10/zoq)
         ctn = kappa**2/np.log(10/zo)/np.log(10/zot)
-    elif (meth == "C40"):
-        usr = np.sqrt(cdn*np.power(u10n, 2))
-        rr = zo*usr/visc_air(Ta)
-        zot = np.where(1.0e-4/np.power(rr, 0.55) > 2.4e-4/np.power(rr, 1.2),
-                       2.4e-4/np.power(rr, 1.2),
-                       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))
-        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"):
         # eq. (3.26) p.38 over sea IFS Documentation cy46r1
         usr = np.sqrt(cdn*np.power(u10n, 2))
@@ -236,11 +222,11 @@ def ctcq_calc(cdn, cd, ctn, cqn, ht, hq, hout, psit, psiq):
     cqn : float
         neutral moisture exchange coefficient
     ht : float
-        original temperature height [m]
+        original temperature sensor height [m]
     hq : float
-        original moisture height    [m]
+        original moisture sensor height    [m]
     hout : float
-        reference height            [m]
+        reference height                   [m]
     psit : float
         heat stability function
     psiq : float
@@ -276,7 +262,7 @@ def get_stabco(meth="S80"):
     alpha, beta, gamma = 0, 0, 0
     if (meth == "S80" or meth == "S88" or meth == "LY04" or
         meth == "UA" or meth == "ecmwf" or meth == "C30" or
-        meth == "C35" or meth == "C40" or meth == "Beljaars"):
+        meth == "C35" or meth == "Beljaars"):
         alpha, beta, gamma = 16, 0.25, 5  # Smith 1980, from Dyer (1974)
     elif (meth == "LP82"):
         alpha, beta, gamma = 16, 0.25, 7
@@ -308,7 +294,7 @@ def psim_calc(zol, meth="S80"):
     """
     if (meth == "ecmwf"):
         psim = psim_ecmwf(zol)
-    elif (meth == "C30" or meth == "C35" or meth == "C40"):
+    elif (meth == "C30" or meth == "C35"):  # or meth == "C40"
         psim = psiu_26(zol, meth)
     elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
         psim = np.where(zol < 0, psim_conv(zol, meth), psi_Bel(zol))
@@ -337,7 +323,7 @@ def psit_calc(zol, meth="S80"):
     if (meth == "ecmwf"):
         psit = np.where(zol < 0, psi_conv(zol, meth),
                         psi_ecmwf(zol))
-    elif (meth == "C30" or meth == "C35" or meth == "C40"):
+    elif (meth == "C30" or meth == "C35"):
         psit = psit_26(zol)
     elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
         psit = np.where(zol < 0, psi_conv(zol, meth), psi_Bel(zol))
@@ -514,7 +500,7 @@ def psiu_26(zol, meth):
                         4*np.arctan(1)/np.sqrt(3), np.nan)
         f = np.power(zol, 2)/(1+np.power(zol, 2))
         psi = np.where(zol < 0, (1-f)*psik+f*psic, psi)
-    elif (meth == "C35" or meth == "C40"):
+    elif (meth == "C35"):
         dzol = np.where(0.35*zol > 50, 50, 0.35*zol)  # stable
         a, b, c, d = 0.7, 3/4, 5, 0.35
         psi = np.where(zol > 0, -(a*zol+b*(zol-c/d)*np.exp(-dzol)+b*c/d),
@@ -949,7 +935,7 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
         Monin-Obukhov length from previous iteration step (m)
     meth : str
         bulk parameterisation method option: "S80", "S88", "LP82", "YT96",
-        "UA", "LY04", "C30", "C35", "C40", "ecmwf", "Beljaars"
+        "UA", "LY04", "C30", "C35", "ecmwf", "Beljaars"
 
     Returns
     -------
@@ -988,7 +974,7 @@ def get_L(L, lat, usr, tsr, qsr, hin, Ta, sst, qair, qsea, wind, monob, psim,
                    (np.log((hin[1]+zo)/zot) -
                     psit_calc((hin[1]+zo)/monob, meth) +
                     psit_calc(zot/monob, meth))))
-        monob = hin[1]/zol 
+        monob = hin[1]/zol
     return tsrv, monob, Rb
 #------------------------------------------------------------------------------
 
@@ -1032,7 +1018,7 @@ def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, dtwl, ct, cq,
         warm layer correction switch
     meth : str
         bulk parameterisation method option: "S80", "S88", "LP82", "YT96", "UA",
-        "LY04", "C30", "C35", "C40", "ecmwf", "Beljaars"
+        "LY04", "C30", "C35", "ecmwf", "Beljaars"
 
     Returns
     -------
@@ -1096,7 +1082,7 @@ def get_strs(hin, monob, wind, zo, zot, zoq, dt, dq, dter, dqer, dtwl, ct, cq,
                                          (np.log(monob/zoq)+5-5*zoq/monob +
                                           5*np.log(hin[2]/monob) +
                                           hin[2]/monob-1))))
-    elif (meth == "C30" or meth == "C35" or meth == "C40"):
+    elif (meth == "C30" or meth == "C35"):
         usr = (wind*kappa/(np.log(hin[0]/zo)-psiu_26(hin[0]/monob, meth)))
         tsr = ((dt+dter*cskin-dtwl*wl)*(kappa/(np.log(hin[1]/zot) -
                                        psit_26(hin[1]/monob))))

get_init.py

 import numpy as np
 import sys
 
-def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, meth,
-             qmeth):
+def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol,
+             meth, qmeth):
     """
     Checks initial input values and sets defaults if needed
 
@@ -17,7 +17,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
         sea surface temperature in K
     lat : float
         latitude (deg), default 45deg
-    hum : float
+    hum : array
         relative humidity, if None is set to 80%
     P : float
         air pressure (hPa), default 1013hPa
@@ -40,14 +40,14 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
         default else [1, 1.2, 800]
     L : int
         Monin-Obukhov length definition options
-    tol : float
+    tol : array
         4x1 or 7x1 [option, lim1-3 or lim1-6]
         option : 'flux' to set tolerance limits for fluxes only lim1-3
         option : 'ref' to set tolerance limits for height adjustment lim-1-3
         option : 'all' to set tolerance limits for both fluxes and height
                  adjustment lim1-6 ['all', 0.01, 0.01, 5e-05, 1e-3, 0.1, 0.1]
     meth : str
-        "S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ecmwf",
+        "S80","S88","LP82","YT96","UA","LY04","C30","C35","ecmwf",
         "Beljaars"
     qmeth : str
         is the saturation evaporation method to use amongst
@@ -90,7 +90,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
         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", "ecmwf", "Beljaars"]:
+                    "ecmwf", "Beljaars"]:
         sys.exit("unknown method")
     if qmeth not in ["HylandWexler", "Hardy", "Preining", "Wexler",
                      "GoffGratch", "WMO", "MagnusTetens", "Buck", "Buck2",
@@ -119,11 +119,10 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
                              or meth == "YT96" or meth == "UA" or
                              meth == "LY04")):
         cskin = 0
-    elif ((cskin == None) and (meth == "C30" or meth == "C35" or meth == "C40"
+    elif ((cskin == None) and (meth == "C30" or meth == "C35"
                                or meth == "ecmwf" or meth == "Beljaars")):
         cskin = 1
-        if ((skin == None) and (meth == "C30" or meth == "C35"
-                                or meth == "C40")):
+        if ((skin == None) and (meth == "C30" or meth == "C35")):
             skin = "C35"
         elif ((skin == None) and (meth == "ecmwf")):
             skin = "ecmwf"
@@ -131,8 +130,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
             skin = "Beljaars"
     if (wl == None):
         wl = 0
-    if (np.all(gust == None) and (meth == "C30" or meth == "C35" or
-                                  meth == "C40")):
+    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
                                     meth == "Beljaars")):
@@ -148,7 +146,7 @@ def get_init(spd, T, SST, lat, hum, P, Rl, Rs, cskin, skin, wl, gust, L, tol, me
     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 == "Beljaars")):
+                         or meth == "Beljaars")):
         L = "S80"
     elif ((L == None) and (meth == "ecmwf")):
         L = "ecmwf"

toy_ASFC.py

@@ -493,7 +493,7 @@ start_time = time.perf_counter()
 inF = input("Give input file name (data_all.csv or era5_r360x180.nc): \n")
 meth = input("Give prefered method: \n")
 while meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
-                   "C40", "ecmwf","Beljaars"]:
+                   "ecmwf","Beljaars"]:
     print("method unknown")
     meth = input("Give prefered method: \n")
 else:
@@ -519,8 +519,8 @@ if (cskinIn == ''):
                              meth == "LY04")):
         cskinIn = 0
         ext = ext+'noskin_'
-    elif ((cskinIn == None) and (meth == "C30" or meth == "C35" or meth == "C40"
-                               or meth == "ecmwf" or meth == "Beljaars")):
+    elif ((cskinIn == None) and (meth == "C30" or meth == "C35"
+                                 or meth == "ecmwf" or meth == "Beljaars")):
         cskinIn = 1
         ext = ext+'skin_'
 else:
@@ -589,10 +589,10 @@ if ((cskinIn == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
                            or meth == "YT96" or meth == "UA" or
                            meth == "LY04")):
    cskinIn = 0
-elif ((cskinIn == None) and (meth == "C30" or meth == "C35" or meth == "C40"
+elif ((cskinIn == None) and (meth == "C30" or meth == "C35"
                              or meth == "ecmwf" or meth == "Beljaars")):
    cskinIn = 1
-if (np.all(gustIn == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
+if (np.all(gustIn == None) and (meth == "C30" or meth == "C35")):
     gustIn = [1, 1.2, 600]
 elif (np.all(gustIn == None) and (meth == "UA" or meth == "ecmwf")):
     gustIn = [1, 1, 1000]