diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index bcc7e39a6030cfd4881044b951464906b02f4d1f..2a232856ad0419b2099eb8733719e3037e237cc2 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -53,8 +53,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
"S80","S88","LP82","YT96","UA","LY04","C30","C35","C40","ERA5"
qmeth : str
is the saturation evaporation method to use amongst
- "HylandWexler","Hardy","Preining","Wexler","GoffGratch","CIMO",
- "MagnusTetens","Buck","Buck2","WMO","WMO2000","Sonntag","Bolton",
+ "HylandWexler","Hardy","Preining","Wexler","GoffGratch","WMO",
+ "MagnusTetens","Buck","Buck2","WMO2018","Sonntag","Bolton",
"IAPWS","MurphyKoop"]
default is Buck2
tol : float
@@ -76,10 +76,10 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
Returns
-------
res : array that contains
- 1. momentum flux (W/m^2)
+ 1. momentum flux (N/m^2)
2. sensible heat (W/m^2)
3. latent heat (W/m^2)
- 4. Monin-Obhukov length (mb)
+ 4. Monin-Obhukov length (m)
5. drag coefficient (cd)
6. neutral drag coefficient (cdn)
7. heat exhange coefficient (ct)
@@ -104,10 +104,15 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
26. temperature at reference height (tref)
27. specific humidity at reference height (qref)
28. number of iterations until convergence
- ind : int
- the indices in the matrix for the points that did not converge
- after the maximum number of iterations
- The code is based on bform.f and flux_calc.R modified by S. Biri
+ 29. cool-skin temperature depression (dter)
+ 30. cool-skin humidity depression (dqer)
+ 31. specific humidity of air (qair)
+ 32. specific humidity at sea surface (qsea)
+ 33. downward longwave radiation (Rl)
+ 34. downward shortwave radiation (Rs)
+ 35. downward net longwave radiation (Rnl)
+
+ 2020 / Author S. Biri
"""
logging.basicConfig(filename='flux_calc.log',
format='%(asctime)s %(message)s',level=logging.INFO)
@@ -135,6 +140,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
np.nanmedian(qsea), np.nanmedian(qair))
if (np.all(np.isnan(qsea)) or np.all(np.isnan(qair))):
print("qsea and qair cannot be nan")
+ # tlapse = gamma_moist(SST, T, qair/1000) lapse rate can be computed like this
dt = Ta - sst
dq = qair - qsea
# first guesses
@@ -289,7 +295,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
elif (tol[0] == 'all'):
new = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
np.copy(tau), np.copy(sensible), np.copy(latent)])
- d = np.fabs(new-old)
+ d = np.abs(new-old)
if (tol[0] == 'flux'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
(d[2, :] > tol[3]))
@@ -305,6 +311,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None,
else:
ii = True
itera[ind] = -1
+ # itera = np.where(itera > n, -1, itera)
logging.info('method %s | # of iterations:%s', meth, it)
logging.info('method %s | # of points that did not converge :%s', meth,
ind[0].size)
diff --git a/Documentation.pdf b/Documentation.pdf
index 581a3fee7b30335ac5992acf745a39b327c8d9b3..35ffac4308bd8857772f8a4709e8d2b25e049479 100644
Binary files a/Documentation.pdf and b/Documentation.pdf differ
diff --git a/data_Mxtr.nc b/data_Mxtr.nc
deleted file mode 100644
index f30b901aaf12ae257c5bc247f81df37e63270198..0000000000000000000000000000000000000000
Binary files a/data_Mxtr.nc and /dev/null differ
diff --git a/data_all.nc b/data_all.nc
new file mode 100644
index 0000000000000000000000000000000000000000..b503e7e8c2c7089cf46f55f4d75266a153a9572f
Binary files /dev/null and b/data_all.nc differ
diff --git a/data_mod.nc b/data_mod.nc
deleted file mode 100644
index 90dc029b78c0f9eb958ded8ba8daa43a5d0cf759..0000000000000000000000000000000000000000
Binary files a/data_mod.nc and /dev/null differ
diff --git a/data_trp.nc b/data_trp.nc
deleted file mode 100644
index b6deea86dacab0d4af5d8a13ddfaee37b8f99ef7..0000000000000000000000000000000000000000
Binary files a/data_trp.nc and /dev/null differ
diff --git a/data_xtr.nc b/data_xtr.nc
deleted file mode 100644
index 2d3f3f93d6eedccf81f749680847c6c9218c98a9..0000000000000000000000000000000000000000
Binary files a/data_xtr.nc and /dev/null differ
diff --git a/era5_r360x180.nc b/era5_r360x180.nc
new file mode 100644
index 0000000000000000000000000000000000000000..f9954e6936de308c8ca63f430083af45af090ecb
Binary files /dev/null and b/era5_r360x180.nc differ
diff --git a/flux_subs.py b/flux_subs.py
index d69ec4bf0871e8c84f6de007e576246575499a14..cb31715b4305fb75bf763ad4068c2b95f83a4624 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -32,7 +32,7 @@ def cdn_calc(u10n, Ta, Tp, lat, meth="S80"):
np.where((u10n <= 25) & (u10n >= 11),
(0.49+0.065*u10n)*0.001, 1.14*0.001))
elif (meth == "S88" or meth == "UA" or meth == "ERA5" or meth == "C30" or
- meth == "C35" or meth == "C40"):
+ meth == "C35" or meth == "C40" or meth == "Beljaars"):
cdn = cdn_from_roughness(u10n, Ta, None, lat, meth)
elif (meth == "YT96"):
# for u<3 same as S80
@@ -101,7 +101,7 @@ def cdn_from_roughness(u10n, Ta, Tp, lat, meth="S88"):
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 == "ERA5"):
+ elif ((meth == "ERA5" 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
else:
@@ -209,7 +209,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, meth="S80"):
# 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 == "ERA5"):
+ elif (meth == "ERA5" or meth == "Beljaars"):
# eq. (3.26) p.38 over sea IFS Documentation cy46r1
usr = np.sqrt(cdn*np.power(u10n, 2))
zot = 0.40*visc_air(Ta)/usr
@@ -274,8 +274,8 @@ 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 == "ERA5" or meth == "C30" or meth == "C35" or
- meth == "C40"):
+ meth == "UA" or meth == "ERA5" or meth == "C30" or
+ meth == "C35" or meth == "C40" 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,6 +308,8 @@ def psim_calc(zol, meth="S80"):
psim = psim_era5(zol)
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))
else:
psim = np.where(zol < 0, psim_conv(zol, meth),
psim_stab(zol, meth))
@@ -334,6 +336,8 @@ def psit_calc(zol, meth="S80"):
psi_era5(zol))
elif (meth == "C30" or meth == "C35" or meth == "C40"):
psit = psit_26(zol)
+ elif (meth == "Beljaars"): # Beljaars (1997) eq. 16, 17
+ psit = np.where(zol < 0, psi_conv(zol, meth), psi_Bel(zol))
else:
psit = np.where(zol < 0, psi_conv(zol, meth),
psi_stab(zol, meth))
@@ -341,6 +345,26 @@ def psit_calc(zol, meth="S80"):
# ---------------------------------------------------------------------
+def psi_Bel(zol):
+ """ Calculates momentum/heat stability function
+
+ Parameters
+ ----------
+ zol : float
+ height over MO length
+ meth : str
+ parameterisation method
+
+ Returns
+ -------
+ psit : float
+ """
+ a, b, c, d = 0.7, 0.75, 5, 0.35
+ psi = -(a*zol+b*(zol-c/d)*np.exp(-d*zol)+b*c/d)
+ return psi
+# ---------------------------------------------------------------------
+
+
def psi_era5(zol):
""" Calculates heat stability function for stable conditions
for method ERA5
@@ -578,8 +602,11 @@ def get_skin(sst, qsea, rho, Rl, Rs, Rnl, cp, lv, tkt, usr, tsr, qsr, lat):
Returns
-------
dter : float
+ cool-skin temperature depression
dqer : float
-
+ cool-skin humidity depression
+ tkt : float
+ cool skin thickness
"""
# coded following Saunders (1967) with lambda = 6
g = gc(lat, None)
@@ -717,7 +744,7 @@ def get_L(L, lat, usr, tsr, qsr, t10n, tv10n, qair, h_in, T, Ta, th, tv, sst,
(np.log((h_in[0]+zo)/zot) -
psit_calc((h_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
- monob = h_in[0]/zol
+ monob = h_in[0]/zol
return tsrv, monob
#------------------------------------------------------------------------------
diff --git a/get_init.py b/get_init.py
index c5a346d8969749d1382dbb566f7547926c49e9c8..8b6ce7b4796af8a30c16ee62ed64b45341ceb43a 100644
--- a/get_init.py
+++ b/get_init.py
@@ -73,6 +73,7 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
MO length switch
"""
+ # check if input is correct (type, size, value) and set defaults
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")
@@ -82,11 +83,11 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
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"]:
+ "C40", "ERA5", "Beljaars"]:
sys.exit("unknown method")
- if qmeth not in ["HylandWexler", "Hardy", "Preining", "Wexler", "CIMO",
- "GoffGratch", "MagnusTetens", "Buck", "Buck2", "WMO",
- "WMO2000", "Sonntag", "Bolton", "IAPWS", "MurphyKoop"]:
+ if qmeth not in ["HylandWexler", "Hardy", "Preining", "Wexler",
+ "GoffGratch", "WMO", "MagnusTetens", "Buck", "Buck2",
+ "WMO2018", "Sonntag", "Bolton", "IAPWS", "MurphyKoop"]:
sys.exit("unknown q-method")
if (np.all(np.isnan(spd)) or np.all(np.isnan(T)) or np.all(np.isnan(SST))):
sys.exit("input wind, T or SST is empty")
@@ -107,11 +108,13 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
meth == "LY04")):
cskin = 0
elif ((cskin == None) and (meth == "C30" or meth == "C35" or meth == "C40"
- or meth == "ERA5")):
+ or meth == "ERA5" or meth == "Beljaars")):
cskin = 1
- 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" or
+ meth == "C40")):
gust = [1, 1.2, 600]
- elif (np.all(gust == None) and (meth == "UA" or meth == "ERA5")):
+ elif (np.all(gust == None) and (meth == "UA" or meth == "ERA5" or
+ meth == "Beljaars")):
gust = [1, 1, 1000]
elif np.all(gust == None):
gust = [1, 1.2, 800]
@@ -124,7 +127,7 @@ def get_init(spd, T, SST, lat, P, Rl, Rs, cskin, gust, L, tol, meth, qmeth):
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 == "C40" or meth == "Beljaars")):
L = "S80"
elif ((L == None) and (meth == "ERA5")):
L = "ERA5"
diff --git a/hum_subs.py b/hum_subs.py
index 8528b80c1a8eda0384a7bc910fe8e921c44313df..b244b44b137a81a1cf26fa200a8f1d99d5c60170 100644
--- a/hum_subs.py
+++ b/hum_subs.py
@@ -24,19 +24,16 @@ def VaporPressure(temp, P, phase, meth):
'ice' : Calculate vapor pressure over ice
meth : str
formula to be used
- MartiMauersberger : vaporpressure formula from Marti Mauersberger
+ Hardy : vaporpressure formula from Hardy (1998)
MagnusTetens : vaporpressure formula from Magnus Tetens
GoffGratch : vaporpressure formula from Goff Gratch
- Buck_original : vaporpressure formula from Buck (original)
- Buck_manual : vaporpressure formula from the Buck manual
- CIMO : vaporpressure formula recommended by CIMO
- Vaisala : vaporpressure formula from Vaisala
- WMO_Goff : vaporpressure formula from WMO, which should have been Goff
- WMO2000 : vaporpressure formula from WMO (2000) containing a typo
+ Buck : vaporpressure formula from Buck (1981)
+ Buck2 : vaporpressure formula from the Buck (2012)
+ WMO : vaporpressure formula from WMO (1988)
+ WMO2018 : vaporpressure formula from WMO (2018)
Wexler : vaporpressure formula from Wexler (1976)
Sonntag : vaporpressure formula from Sonntag (1994)
Bolton : vaporpressure formula from Bolton (1980)
- Fukuta : vaporpressure formula from Fukuta (2003)
HylandWexler : vaporpressure formula from Hyland and Wexler (1983)
IAPWS : vaporpressure formula from IAPWS (2002)
Preining : vaporpressure formula from Preining (2002)
@@ -111,12 +108,6 @@ def VaporPressure(temp, P, phase, meth):
1.3816e-7*(np.power(10, 11.344*(1-T/Ts))-1) +
8.1328e-3*(np.power(10, -3.49149*(Ts/T-1))-1) +
np.log10(ews))
- if (meth == 'CIMO'):
- """Source: Annex 4B, Guide to Meteorological Instruments and
- Methods of Observation, WMO Publication No 8, 7th edition, Geneva,
- 2008. (CIMO Guide)"""
- Psat = (6.112*np.exp(17.62*temp/(243.12+temp)) *
- (1.0016+3.15e-6*P-0.074/P))
if (meth == 'MagnusTetens'):
"""Source: Murray, F. W., On the computation of \
saturation vapor pressure, J. Appl. Meteorol., \
@@ -134,7 +125,7 @@ def VaporPressure(temp, P, phase, meth):
if (meth == 'Buck2'):
"""Bucks vapor pressure formulation based on Tetens formula.
Source: Buck Research, Model CR-1A Hygrometer Operating Manual,
- Sep 2001"""
+ May 2012"""
Psat = (6.1121*np.exp((18.678-(temp)/234.5)*(temp)/(257.14+temp)) *
(1+1e-4*(7.2+P*(0.0320)+5.9e-6*np.power(T, 2))))
if (meth == 'WMO'):
@@ -147,18 +138,12 @@ def VaporPressure(temp, P, phase, meth):
Ts = 273.16 # triple point temperature in K
Psat = np.power(10, 10.79574*(1-Ts/T)-5.028*np.log10(T/Ts) +
1.50475e-4*(1-np.power(10, -8.2969*(T/Ts-1))) +
- 0.42873e-3*(np.power(10, -4.76955*(1-Ts/T))-1) +
- 0.78614)
- if (meth == 'WMO2000'):
- """WMO formulation, which is very similar to Goff & Gratch.
- Source : WMO technical regulations, WMO-NO 49, Vol I, General
- Meteorological Standards and Recommended Practices, App. A,
- Corrigendum Aug 2000."""
- Ts = 273.16 # triple point temperature in K
- Psat = np.power(10, 10.79574*(1-Ts/T)-5.028*np.log10(T/Ts) +
- 1.50475e-4*(1-np.power(10, -8.2969*(T/Ts-1))) +
- 0.42873e-3*(np.power(10, -4.76955*(1.-Ts/T))-1) +
+ 0.42873e-3*(np.power(10, 4.76955*(1-Ts/T))-1) + # in eq. 13 is -4.76955; in aerobulk is like this
0.78614)
+ if (meth == 'WMO2018'):
+ """WMO 2018 edition. Annex 4.B, eq. 4.B.1, 4.B.2, 4.B.5 """
+ Psat = 6.112*np.exp(17.62*temp/(243.12+temp))*(1.0016+3.15e-6*P -
+ 0.074/P)
if (meth == 'Sonntag'):
"""Source: Sonntag, D., Advancements in the field of hygrometry,
Meteorol. Z., N. F., 3, 51-66, 1994."""
@@ -413,3 +398,38 @@ def get_hum(hum, T, sst, P, qmeth):
qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (g/kg)
return qair, qsea
#------------------------------------------------------------------------------
+
+
+def gamma_moist(sst, t, q):
+ """
+ Computes the moist adiabatic lapse-rate
+
+ Parameters
+ ----------
+ sst : float
+ sea surface temperature [K]
+ t : float
+ air temperature [K]
+ q : float
+ specific humidity of air [kg/kg]
+
+ Returns
+ -------
+ gamma : float
+ lapse rate [K/m]
+
+ """
+ if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
+ sst = sst+CtoK
+ if (np.nanmin(t) < 200): # if sst in Celsius convert to Kelvin
+ t = t+CtoK
+
+ t = np.maximum(t, 180)
+ q = np.maximum(q, 1e-6)
+ w = q/(1-q) # mixing ratio w = q/(1-q)
+ iRT = 1/(287.05*t)
+ # latent heat of vaporization of water as a function of temperature
+ lv = (2.501-0.00237*(sst-CtoK))*1e6
+ gamma = 9.8*(1+lv*w*iRT)/(1005+np.power(lv, 2)*w*(287.05/461.495)*iRT/t)
+ return gamma
+#------------------------------------------------------------------------------
diff --git a/run_ASFC.py b/run_ASFC.py
new file mode 100644
index 0000000000000000000000000000000000000000..1b8793088d152229d981e03a66e67fef026d0825
--- /dev/null
+++ b/run_ASFC.py
@@ -0,0 +1,507 @@
+"""
+example of running AirSeaFluxCode with
+1. R/V data (data_all.nc) or
+2. one day era5 hourly data (era5_r360x180.nc)
+computes fluxes
+outputs NetCDF4
+and statistics in stats.txt
+
+@author: sbiri
+"""
+#%% import packages
+import matplotlib.pyplot as plt
+import netCDF4 as nc
+import numpy as np
+from AirSeaFluxCode import AirSeaFluxCode
+import time
+from tabulate import tabulate
+#%%
+def reject_outliers(data, m=2):
+ x = np.copy(data)
+ x = np.where(np.abs(x - np.nanmean(x)) < m*np.nanstd(x),
+ x, np.nan)
+ return x
+
+
+def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
+ if (inF == "data_all.nc"):
+ #%% load data_all
+ fid=nc.Dataset('data_all.nc')
+ date = np.array(fid.variables["date"])
+ lon = np.array(fid.variables["lon"])
+ lat = np.array(fid.variables["lat"])
+ spd = np.array(fid.variables["spd"])
+ t = np.array(fid.variables["t"])
+ sst = np.array(fid.variables["sst"])
+ lat = np.array(fid.variables["lat"])
+ rh = np.array(fid.variables["rh"])
+ p = np.array(fid.variables["p"])
+ sw = np.array(fid.variables["sw"])
+ hu = np.array(fid.variables["spd"].getncattr("height"))
+ ht = np.array(fid.variables["t"].getncattr("height"))
+ hin = np.array([hu, ht, ht])
+ del hu, ht
+ fid.close()
+ del fid
+ #%% run AirSeaFluxCode
+ res = AirSeaFluxCode(spd, t, sst, lat=lat, hum=['rh', rh], P=p,
+ hin=hin, Rs=sw, tol=tolIn, gust=gustIn,
+ cskin=cskinIn, meth=meth, out=1, L="ERA5", n=30)
+ #%% save NetCDF4
+ fid = nc.Dataset(outF,'w', format='NETCDF4')
+ fid.createDimension('lon', len(lon))
+ fid.createDimension('lat', len(lat))
+ fid.createDimension('time', None)
+ longitude = fid.createVariable('lon', 'f4', 'lon')
+ latitude = fid.createVariable('lat', 'f4', 'lat')
+ Date = fid.createVariable('Date', 'i4', 'time')
+ tau = fid.createVariable('tau', 'f4', 'time')
+ sensible = fid.createVariable('shf', 'f4', 'time')
+ latent = fid.createVariable('lhf', 'f4', 'time')
+ monob = fid.createVariable('MO', 'f4', 'time')
+ cd = fid.createVariable('cd', 'f4', 'time')
+ cdn = fid.createVariable('cdn', 'f4', 'time')
+ ct = fid.createVariable('ct', 'f4', 'time')
+ ctn = fid.createVariable('ctn', 'f4', 'time')
+ cq = fid.createVariable('cq', 'f4', 'time')
+ cqn = fid.createVariable('cqn', 'f4', 'time')
+ tsrv = fid.createVariable('tsrv', 'f4', 'time')
+ tsr = fid.createVariable('tsr', 'f4', 'time')
+ qsr = fid.createVariable('qsr', 'f4', 'time')
+ usr = fid.createVariable('usr', 'f4', 'time')
+ psim = fid.createVariable('psim', 'f4', 'time')
+ psit = fid.createVariable('psit', 'f4', 'time')
+ psiq = fid.createVariable('psiq', 'f4', 'time')
+ u10n = fid.createVariable('u10n', 'f4', 'time')
+ t10n = fid.createVariable('t10n', 'f4', 'time')
+ tv10n = fid.createVariable('tv10n', 'f4', 'time')
+ q10n = fid.createVariable('q10n', 'f4', 'time')
+ zo = fid.createVariable('zo', 'f4', 'time')
+ zot = fid.createVariable('zot', 'f4', 'time')
+ zoq = fid.createVariable('zoq', 'f4', 'time')
+ urefs = fid.createVariable('uref', 'f4', 'time')
+ trefs = fid.createVariable('tref', 'f4', 'time')
+ qrefs = fid.createVariable('qref', 'f4', 'time')
+ itera = fid.createVariable('iter', 'i4', 'time')
+ dter = fid.createVariable('dter', 'f4', 'time')
+ dqer = fid.createVariable('dqer', 'f4', 'time')
+ qair = fid.createVariable('qair', 'f4', 'time')
+ qsea = fid.createVariable('qsea', 'f4', 'time')
+ Rl = fid.createVariable('Rl', 'f4', 'time')
+ Rs = fid.createVariable('Rs', 'f4', 'time')
+ Rnl = fid.createVariable('Rnl', 'f4', 'time')
+
+ longitude[:] = lon
+ latitude[:] = lat
+ Date[:] = date
+ tau[:] = res[0]
+ sensible[:] = res[1]
+ latent[:] = res[2]
+ monob[:] = res[3]
+ cd[:] = res[4]
+ cdn[:] = res[5]
+ ct[:] = res[6]
+ ctn[:] = res[7]
+ cq[:] = res[8]
+ cqn[:] = res[9]
+ tsrv[:] = res[10]
+ tsr[:] = res[11]
+ qsr[:] = res[12]
+ usr[:] = res[13]
+ psim[:] = res[14]
+ psit[:] = res[15]
+ psiq[:] = res[16]
+ u10n[:] = res[17]
+ t10n[:] = res[18]
+ tv10n[:] = res[19]
+ q10n[:] = res[20]
+ zo[:] = res[21]
+ zot[:] = res[22]
+ zoq[:] = res[23]
+ urefs[:] = res[24]
+ trefs[:] = res[25]
+ qrefs[:] = res[26]
+ itera[:] = res[27]
+ dter[:] = res[28]
+ dqer[:] = res[29]
+ qair[:] = res[30]
+ qsea[:] = res[31]
+ Rl = res[32]
+ Rs = res[33]
+ Rnl = res[34]
+ longitude.long_name = 'Longitude'
+ longitude.units = 'degrees East'
+ latitude.long_name = 'Latitude'
+ latitude.units = 'degrees North'
+ Date.long_name = "calendar date"
+ Date.units = "YYYYMMDD UTC"
+ tau.long_name = 'Wind stress'
+ tau.units = 'N/m^2'
+ sensible.long_name = 'Sensible heat fluxe'
+ sensible.units = 'W/m^2'
+ latent.long_name = 'Latent heat flux'
+ latent.units = 'W/m^2'
+ monob.long_name = 'Monin-Obukhov length'
+ monob.units = 'm'
+ cd.long_name = 'Drag coefficient'
+ cd.units = ''
+ cdn.long_name = 'Neutral Drag coefficient'
+ cdn.units = ''
+ ct.long_name = 'Heat exchange coefficient'
+ ct.units = ''
+ ctn.long_name = 'Neutral Heat exchange coefficient'
+ ctn.units = ''
+ cq.long_name = 'Moisture exchange coefficient'
+ cq.units = ''
+ cqn.long_name = 'Neutral Moisture exchange coefficient'
+ cqn.units = ''
+ tsrv.long_name = 'star virtual temperature'
+ tsrv.units = 'degrees Celsius'
+ tsr.long_name = 'star temperature'
+ tsr.units = 'degrees Celsius'
+ qsr.long_name = 'star specific humidity'
+ qsr.units = 'gr/kgr'
+ usr.long_name = 'friction velocity'
+ usr.units = 'm/s'
+ psim.long_name = 'Momentum stability function'
+ psit.long_name = 'Heat stability function'
+ u10n.long_name = '10m neutral wind speed'
+ u10n.units = 'm/s'
+ t10n.long_name = '10m neutral temperature'
+ t10n.units = 'degrees Celsius'
+ tv10n.long_name = '10m neutral virtual temperature'
+ tv10n.units = 'degrees Celsius'
+ q10n.long_name = '10m neutral specific humidity'
+ q10n.units = 'gr/kgr'
+ zo.long_name = 'momentum roughness length'
+ zo.units = 'm'
+ zot.long_name = 'temperature roughness length'
+ zot.units = 'm'
+ zoq.long_name = 'moisture roughness length'
+ zoq.units = 'm'
+ urefs.long_name = 'wind speed at ref height'
+ urefs.units = 'm/s'
+ trefs.long_name = 'temperature at ref height'
+ trefs.units = 'degrees Celsius'
+ qrefs.long_name = 'specific humidity at ref height'
+ qrefs.units = 'gr/kgr'
+ qair.long_name = 'specific humidity of air'
+ qair.units = 'gr/kgr'
+ qsea.long_name = 'specific humidity over water'
+ qsea.units = 'gr/kgr'
+ itera.long_name = 'number of iterations'
+ fid.close()
+ #%% delete variables
+ del longitude, latitude, Date, tau, sensible, latent, monob, cd, cdn
+ del ct, ctn, cq, cqn, tsrv, tsr, qsr, usr, psim, psit, psiq, u10n, t10n
+ del tv10n, q10n, zo, zot, zoq, urefs, trefs, qrefs, itera, dter, dqer
+ del qair, qsea, Rl, Rs, Rnl
+ del t, rh, date, p, sw, spd, hin
+
+ elif (inF == 'era5_r360x180.nc'):
+ #%% load era5_r360x180.nc
+ fid = nc.Dataset(inF)
+ lon = np.array(fid.variables["lon"])
+ lat = np.array(fid.variables["lat"])
+ T = np.array(fid.variables["t2m"])
+ tim = np.array(fid.variables["time"])
+ Td = np.array(fid.variables["d2m"])
+ sst = np.array(fid.variables["sst"])
+ sst = np.where(sst < -100, np.nan, sst)
+ p = np.array(fid.variables["msl"])/100 # to set hPa
+ lw = np.array(fid.variables["strd"])/60/60
+ sw = np.array(fid.variables["ssrd"])/60/60
+ u = np.array(fid.variables["u10"])
+ v = np.array(fid.variables["v10"])
+ lsm = np.array(fid.variables["lsm"])
+ fid.close()
+ spd = np.sqrt(np.power(u, 2)+np.power(v, 2))
+ del u, v, fid
+ lsm = np.where(lsm > 0, np.nan, 1) # reverse 0 on land 1 over ocean
+ hin = np.array([10, 2, 2])
+ latIn = np.tile(lat, (len(lon), 1)).T.reshape(len(lon)*len(lat))
+ #%% run AirSeaFluxCode
+ res = np.zeros((len(tim),len(lon)*len(lat), 35))
+ # reshape input and run code
+ for x in range(len(tim)):
+ a = AirSeaFluxCode(spd.reshape(len(tim), len(lon)*len(lat))[x, :],
+ T.reshape(len(tim), len(lon)*len(lat))[x, :],
+ sst.reshape(len(tim), len(lon)*len(lat))[x, :],
+ lat=latIn,
+ hum=['Td', Td.reshape(len(tim), len(lon)*len(lat))[x, :]],
+ P=p.reshape(len(tim), len(lon)*len(lat))[x, :],
+ hin=hin,
+ Rs=sw.reshape(len(tim), len(lon)*len(lat))[x, :],
+ Rl=lw.reshape(len(tim), len(lon)*len(lat))[x, :],
+ gust=gustIn, cskin=cskinIn, tol=tolIn, qmeth='WMO',
+ meth=meth, n=30, L="ERA5")
+ res[x, :, :] = a.T
+ del a
+ #%% save NetCDF4
+ fid = nc.Dataset(outF,'w', format='NETCDF4')
+ fid.createDimension('lon', len(lon))
+ fid.createDimension('lat', len(lat))
+ fid.createDimension('time', None)
+ longitude = fid.createVariable('lon', 'f4', 'lon')
+ latitude = fid.createVariable('lat', 'f4', 'lat')
+ Date = fid.createVariable('Date', 'i4', 'time')
+ tau = fid.createVariable('tau', 'f4', ('time','lat','lon'))
+ sensible = fid.createVariable('shf', 'f4', ('time','lat','lon'))
+ latent = fid.createVariable('lhf', 'f4', ('time','lat','lon'))
+ monob = fid.createVariable('MO', 'f4', ('time','lat','lon'))
+ cd = fid.createVariable('cd', 'f4', ('time','lat','lon'))
+ cdn = fid.createVariable('cdn', 'f4', ('time','lat','lon'))
+ ct = fid.createVariable('ct', 'f4', ('time','lat','lon'))
+ ctn = fid.createVariable('ctn', 'f4', ('time','lat','lon'))
+ cq = fid.createVariable('cq', 'f4', ('time','lat','lon'))
+ cqn = fid.createVariable('cqn', 'f4', ('time','lat','lon'))
+ tsrv = fid.createVariable('tsrv', 'f4', ('time','lat','lon'))
+ tsr = fid.createVariable('tsr', 'f4', ('time','lat','lon'))
+ qsr = fid.createVariable('qsr', 'f4', ('time','lat','lon'))
+ usr = fid.createVariable('usr', 'f4', ('time','lat','lon'))
+ psim = fid.createVariable('psim', 'f4', ('time','lat','lon'))
+ psit = fid.createVariable('psit', 'f4', ('time','lat','lon'))
+ psiq = fid.createVariable('psiq', 'f4', ('time','lat','lon'))
+ u10n = fid.createVariable('u10n', 'f4', ('time','lat','lon'))
+ t10n = fid.createVariable('t10n', 'f4', ('time','lat','lon'))
+ tv10n = fid.createVariable('tv10n', 'f4', ('time','lat','lon'))
+ q10n = fid.createVariable('q10n', 'f4', ('time','lat','lon'))
+ zo = fid.createVariable('zo', 'f4', ('time','lat','lon'))
+ zot = fid.createVariable('zot', 'f4', ('time','lat','lon'))
+ zoq = fid.createVariable('zoq', 'f4', ('time','lat','lon'))
+ urefs = fid.createVariable('uref', 'f4', ('time','lat','lon'))
+ trefs = fid.createVariable('tref', 'f4', ('time','lat','lon'))
+ qrefs = fid.createVariable('qref', 'f4', ('time','lat','lon'))
+ itera = fid.createVariable('iter', 'i4', ('time','lat','lon'))
+ dter = fid.createVariable('dter', 'f4', ('time','lat','lon'))
+ dqer = fid.createVariable('dqer', 'f4', ('time','lat','lon'))
+ qair = fid.createVariable('qair', 'f4', ('time','lat','lon'))
+ qsea = fid.createVariable('qsea', 'f4', ('time','lat','lon'))
+ Rl = fid.createVariable('Rl', 'f4', ('time','lat','lon'))
+ Rs = fid.createVariable('Rs', 'f4', ('time','lat','lon'))
+ Rnl = fid.createVariable('Rnl', 'f4', ('time','lat','lon'))
+
+ longitude[:] = lon
+ latitude[:] = lat
+ Date[:] = tim
+ tau[:] = res[:, :, 0].reshape((len(tim), len(lat), len(lon)))*lsm
+ sensible[:] = res[:, :, 1].reshape((len(tim), len(lat), len(lon)))*lsm
+ latent[:] = res[:, :, 2].reshape((len(tim), len(lat), len(lon)))*lsm
+ monob[:] = res[:, :, 3].reshape((len(tim), len(lat), len(lon)))*lsm
+ cd[:] = res[:, :, 4].reshape((len(tim), len(lat), len(lon)))*lsm
+ cdn[:] = res[:, :, 5].reshape((len(tim), len(lat), len(lon)))*lsm
+ ct[:] = res[:, :, 6].reshape((len(tim), len(lat), len(lon)))*lsm
+ ctn[:] = res[:, :, 7].reshape((len(tim), len(lat), len(lon)))*lsm
+ cq[:] = res[:, :, 8].reshape((len(tim), len(lat), len(lon)))*lsm
+ cqn[:] = res[:, :, 9].reshape((len(tim), len(lat), len(lon)))*lsm
+ tsrv[:] = res[:, :, 10].reshape((len(tim), len(lat), len(lon)))*lsm
+ tsr[:] = res[:, :, 11].reshape((len(tim), len(lat), len(lon)))*lsm
+ qsr[:] = res[:, :, 12].reshape((len(tim), len(lat), len(lon)))*lsm
+ usr[:] = res[:, :, 13].reshape((len(tim), len(lat), len(lon)))*lsm
+ psim[:] = res[:, :, 14].reshape((len(tim), len(lat), len(lon)))*lsm
+ psit[:] = res[:, :, 15].reshape((len(tim), len(lat), len(lon)))*lsm
+ psiq[:] = res[:, :, 16].reshape((len(tim), len(lat), len(lon)))*lsm
+ u10n[:] = res[:, :, 17].reshape((len(tim), len(lat), len(lon)))*lsm
+ t10n[:] = res[:, :, 18].reshape((len(tim), len(lat), len(lon)))*lsm
+ tv10n[:] = res[:, :, 19].reshape((len(tim), len(lat), len(lon)))*lsm
+ q10n[:] = res[:, :, 20].reshape((len(tim), len(lat), len(lon)))*lsm
+ zo[:] = res[:, :, 21].reshape((len(tim), len(lat), len(lon)))*lsm
+ zot[:] = res[:, :, 22].reshape((len(tim), len(lat), len(lon)))*lsm
+ zoq[:] = res[:, :, 23].reshape((len(tim), len(lat), len(lon)))*lsm
+ urefs[:] = res[:, :, 24].reshape((len(tim), len(lat), len(lon)))*lsm
+ trefs[:] = res[:, :, 25].reshape((len(tim), len(lat), len(lon)))*lsm
+ qrefs[:] = res[:, :, 26].reshape((len(tim), len(lat), len(lon)))*lsm
+ itera[:] = res[:, :, 27].reshape((len(tim), len(lat), len(lon)))*lsm
+ dter[:] = res[:, :, 28].reshape((len(tim), len(lat), len(lon)))*lsm
+ dqer[:] = res[:, :, 29].reshape((len(tim), len(lat), len(lon)))*lsm
+ qair[:] = res[:, :, 30].reshape((len(tim), len(lat), len(lon)))*lsm
+ qsea[:] = res[:, :, 31].reshape((len(tim), len(lat), len(lon)))*lsm
+ Rl = res[:, :, 32].reshape((len(tim), len(lat), len(lon)))
+ Rs = res[:, :, 33].reshape((len(tim), len(lat), len(lon)))
+ Rnl = res[:, :, 34].reshape((len(tim), len(lat), len(lon)))
+ longitude.long_name = 'Longitude'
+ longitude.units = 'degrees East'
+ latitude.long_name = 'Latitude'
+ latitude.units = 'degrees North'
+ Date.long_name = "calendar date"
+ Date.units = "YYYYMMDD UTC"
+ tau.long_name = 'Wind stress'
+ tau.units = 'N/m^2'
+ sensible.long_name = 'Sensible heat fluxe'
+ sensible.units = 'W/m^2'
+ latent.long_name = 'Latent heat flux'
+ latent.units = 'W/m^2'
+ monob.long_name = 'Monin-Obukhov length'
+ monob.units = 'm'
+ cd.long_name = 'Drag coefficient'
+ cd.units = ''
+ cdn.long_name = 'Neutral Drag coefficient'
+ cdn.units = ''
+ ct.long_name = 'Heat exchange coefficient'
+ ct.units = ''
+ ctn.long_name = 'Neutral Heat exchange coefficient'
+ ctn.units = ''
+ cq.long_name = 'Moisture exchange coefficient'
+ cq.units = ''
+ cqn.long_name = 'Neutral Moisture exchange coefficient'
+ cqn.units = ''
+ tsrv.long_name = 'star virtual temperature'
+ tsrv.units = 'degrees Celsius'
+ tsr.long_name = 'star temperature'
+ tsr.units = 'degrees Celsius'
+ qsr.long_name = 'star specific humidity'
+ qsr.units = 'gr/kgr'
+ usr.long_name = 'friction velocity'
+ usr.units = 'm/s'
+ psim.long_name = 'Momentum stability function'
+ psit.long_name = 'Heat stability function'
+ u10n.long_name = '10m neutral wind speed'
+ u10n.units = 'm/s'
+ t10n.long_name = '10m neutral temperature'
+ t10n.units = 'degrees Celsius'
+ tv10n.long_name = '10m neutral virtual temperature'
+ tv10n.units = 'degrees Celsius'
+ q10n.long_name = '10m neutral specific humidity'
+ q10n.units = 'gr/kgr'
+ zo.long_name = 'momentum roughness length'
+ zo.units = 'm'
+ zot.long_name = 'temperature roughness length'
+ zot.units = 'm'
+ zoq.long_name = 'moisture roughness length'
+ zoq.units = 'm'
+ urefs.long_name = 'wind speed at ref height'
+ urefs.units = 'm/s'
+ trefs.long_name = 'temperature at ref height'
+ trefs.units = 'degrees Celsius'
+ qrefs.long_name = 'specific humidity at ref height'
+ qrefs.units = 'gr/kgr'
+ qair.long_name = 'specific humidity of air'
+ qair.units = 'gr/kgr'
+ qsea.long_name = 'specific humidity over water'
+ qsea.units = 'gr/kgr'
+ itera.long_name = 'number of iterations'
+ fid.close()
+ #%% delete variables
+ del longitude, latitude, Date, tau, sensible, latent, monob, cd, cdn
+ del ct, ctn, cq, cqn, tsrv, tsr, qsr, usr, psim, psit, psiq, u10n, t10n
+ del tv10n, q10n, zo, zot, zoq, urefs, trefs, qrefs, itera, dter, dqer
+ del qair, qsea, Rl, Rs, Rnl
+ del T, Td, tim, p, lw, sw, lsm, spd, hin, latIn
+ return res, lon, lat
+#%% run function
+start_time = time.perf_counter()
+inF = input("Give input file name (data_all.nc or era5_r360x180.nc): \n")
+gustIn = input("Give gustiness option (to use default press enter): \n")
+if (gustIn == ''):
+ gustIn = None
+else:
+ gustIn = np.asarray(eval(gustIn), dtype=float)
+cskinIn = input("Give cool skin option (to use default press enter): \n")
+if (cskinIn == ''):
+ cskinIn = None
+else:
+ cskinIn = int(cskinIn)
+tolIn = input("Give tolerance option (to use default press enter): \n")
+if (tolIn == ''):
+ tolIn = None
+else:
+ tolIn = eval(tolIn)
+meth = input("Give prefered method: \n")
+while meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
+ "C40", "ERA5","Beljaars"]:
+ print("method unknown")
+ meth = input("Give prefered method: \n")
+else:
+ meth = meth #[meth]
+outF = input("Give path and output file name: \n")
+if (outF == ''):
+ outF = inF[:-3]+"_"+meth+".nc"
+else:
+ outF = outF
+
+print("\n run_ASFC.py, started for method "+meth)
+
+res, lon, lat = run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth)
+print("run_ASFC.py took ", np.round((time.perf_counter()-start_time)/60, 2),
+ "minutes to run")
+
+#%% generate flux plots
+if (inF == 'era5_r360x180.nc'):
+ cm = plt.cm.get_cmap('RdYlBu')
+ ttl = ["tau (Nm$^{-2}$)", "shf (Wm$^{-2}$)", "lhf (Wm$^{-2}$)"]
+ for i in range(3):
+ plt.figure()
+ plt.contourf(lon, lat,
+ np.nanmean(res[:, :, i], axis=0).reshape(len(lat),
+ len(lon)),
+ 100, cmap=cm)
+ plt.colorbar()
+ plt.tight_layout()
+ plt.xlabel("Longitude")
+ plt.ylabel("Latitude")
+ plt.title(meth+', '+ttl[i])
+ plt.savefig('./'+ttl[i][:3]+'_'+meth+'.png', dpi=300, bbox_inches='tight')
+elif (inF == "data_all.nc"):
+ ttl = ["tau (Nm$^{-2}$)", "shf (Wm$^{-2}$)", "lhf (Wm$^{-2}$)"]
+ for i in range(3):
+ plt.figure()
+ plt.plot(res[i],'.c', markersize=1)
+ plt.title(meth+', '+ttl[i])
+ plt.savefig('./'+ttl[i][:3]+'_'+meth+'.png', dpi=300, bbox_inches='tight')
+
+#%% generate txt file with statistic
+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"
+ or meth == "ERA5" or meth == "Beljaars")):
+ cskinIn = 1
+if (np.all(gustIn == None) and (meth == "C30" or meth == "C35" or meth == "C40")):
+ gustIn = [1, 1.2, 600]
+elif (np.all(gustIn == None) and (meth == "UA" or meth == "ERA5")):
+ gustIn = [1, 1, 1000]
+elif np.all(gustIn == None):
+ gustIn = [1, 1.2, 800]
+elif ((np.size(gustIn) < 3) and (gustIn == 0)):
+ gust = [0, 0, 0]
+if (tolIn == None):
+ tolIn = ['flux', 0.01, 1, 1]
+
+print("Input summary", file=open('./stats.txt', 'a'))
+print('input file name: {}, \n method: {}, \n gustiness: {}, \n cskin: {},'
+ ' \n tolerance: {}'.format(inF, meth, gustIn, cskinIn, tolIn),
+ file=open('./stats.txt', 'a'))
+ttl = np.asarray(["tau ", "shf ", "lhf ", "L ", "cd ", "cdn ",
+ "ct ", "ctn ", "cq ", "cqn ", "tsrv ", "tsr ",
+ "qsr ", "usr ", "psim ", "psit ", "psiq ", "u10n ",
+ "t10n ", "tv10n", "q10n ", "zo ", "zot ", "zoq ",
+ "urefs", "trefs", "qrefs", "itera", "dter ", "dqer ",
+ "qair ", "qsea ", "Rl ", "Rs ", "Rnl "])
+header = ["var", "mean", "median", "min", "max", "5%", "95%"]
+n = np.shape(res)
+stats = np.copy(ttl)
+if (inF == 'era5_r360x180.nc'):
+ stats = np.c_[stats, np.nanmean(res.reshape(n[0]*n[1], n[2]), axis=0)]
+ stats = np.c_[stats, np.nanmedian(res.reshape(n[0]*n[1], n[2]), axis=0)]
+ stats = np.c_[stats, np.nanmin(res.reshape(n[0]*n[1], n[2]), axis=0)]
+ stats = np.c_[stats, np.nanmax(res.reshape(n[0]*n[1], n[2]), axis=0)]
+ stats = np.c_[stats, np.nanpercentile(res.reshape(n[0]*n[1], n[2]), 5,
+ axis=0)]
+ stats = np.c_[stats, np.nanpercentile(res.reshape(n[0]*n[1], n[2]), 95,
+ axis=0)]
+ print(tabulate(stats, headers=header, tablefmt="github", numalign="left",
+ floatfmt=("s", "2.2e", "2.2e", "2.2e", "2.2e", "2.2e",
+ "2.2e")), file=open('./stats.txt', 'a'))
+elif (inF == "data_all.nc"):
+ stats = np.c_[stats, np.nanmean(res, axis=1)]
+ stats = np.c_[stats, np.nanmedian(res, axis=1)]
+ stats = np.c_[stats, np.nanmin(res, axis=1)]
+ stats = np.c_[stats, np.nanmax(res, axis=1)]
+ stats = np.c_[stats, np.nanpercentile(res, 5, axis=1)]
+ stats = np.c_[stats, np.nanpercentile(res, 95, axis=1)]
+ print(tabulate(stats, headers=header, tablefmt="github", numalign="left",
+ floatfmt=("s", "2.2e", "2.2e", "2.2e", "2.2e", "2.2e",
+ "2.2e")), file=open('./stats.txt', 'a'))
+ print('-'*79+'\n', file=open('./stats.txt', 'a'))
+