diff --git a/run_ASFC.py b/run_ASFC.py
index 1b8793088d152229d981e03a66e67fef026d0825..193684051ca4bec6955594feae4bc0384ea469fd 100644
--- a/run_ASFC.py
+++ b/run_ASFC.py
@@ -2,8 +2,8 @@
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
+compute fluxes
+output NetCDF4
and statistics in stats.txt
@author: sbiri
@@ -24,6 +24,34 @@ def reject_outliers(data, m=2):
def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
+ """
+
+
+ Parameters
+ ----------
+ inF : str
+ input filename either data_all.nc or era5_r360x180.nc
+ outF : str
+ output filename
+ gustIn : float
+ gustiness option e.g. [1, 1.2, 800]
+ cskinIn : int
+ cool skin option input 0 or 1
+ tolIn : float
+ tolerance input option e.g. ['all', 0.01, 0.01, 5e-05, 0.01, 1, 1]
+ meth : str
+ parametrisation method option
+
+ Returns
+ -------
+ res : float
+ AirSeaFluxCode output
+ lon : float
+ longitude from input netCDF file
+ lat : float
+ latitude from input netCDF file
+
+ """
if (inF == "data_all.nc"):
#%% load data_all
fid=nc.Dataset('data_all.nc')
@@ -33,7 +61,6 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
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"])
@@ -46,7 +73,7 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
#%% 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)
+ cskin=cskinIn, meth=meth, out=1, L="ecmwf", n=30)
#%% save NetCDF4
fid = nc.Dataset(outF,'w', format='NETCDF4')
fid.createDimension('lon', len(lon))
@@ -85,6 +112,7 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
itera = fid.createVariable('iter', 'i4', 'time')
dter = fid.createVariable('dter', 'f4', 'time')
dqer = fid.createVariable('dqer', 'f4', 'time')
+ dtwl = fid.createVariable('dter', 'f4', 'time')
qair = fid.createVariable('qair', 'f4', 'time')
qsea = fid.createVariable('qsea', 'f4', 'time')
Rl = fid.createVariable('Rl', 'f4', 'time')
@@ -124,11 +152,12 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
itera[:] = res[27]
dter[:] = res[28]
dqer[:] = res[29]
- qair[:] = res[30]
- qsea[:] = res[31]
- Rl = res[32]
- Rs = res[33]
- Rnl = res[34]
+ dtwl[:] = res[30]
+ qair[:] = res[31]
+ qsea[:] = res[32]
+ Rl = res[33]
+ Rs = res[34]
+ Rnl = res[35]
longitude.long_name = 'Longitude'
longitude.units = 'degrees East'
latitude.long_name = 'Latitude'
@@ -221,7 +250,7 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
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))
+ res = np.zeros((len(tim),len(lon)*len(lat), 36))
# reshape input and run code
for x in range(len(tim)):
a = AirSeaFluxCode(spd.reshape(len(tim), len(lon)*len(lat))[x, :],
@@ -234,7 +263,7 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
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")
+ meth=meth, n=30, L="ecmwf")
res[x, :, :] = a.T
del a
#%% save NetCDF4
@@ -275,6 +304,7 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
itera = fid.createVariable('iter', 'i4', ('time','lat','lon'))
dter = fid.createVariable('dter', 'f4', ('time','lat','lon'))
dqer = fid.createVariable('dqer', 'f4', ('time','lat','lon'))
+ dtwl = fid.createVariable('dtwl', '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'))
@@ -313,12 +343,13 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
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)))
+ dtwl[:] = res[:, :, 29].reshape((len(tim), len(lat), len(lon)))*lsm
+ dqer[:] = res[:, :, 30].reshape((len(tim), len(lat), len(lon)))*lsm
+ qair[:] = res[:, :, 31].reshape((len(tim), len(lat), len(lon)))*lsm
+ qsea[:] = res[:, :, 32].reshape((len(tim), len(lat), len(lon)))*lsm
+ Rl = res[:, :, 33].reshape((len(tim), len(lat), len(lon)))
+ Rs = res[:, :, 34].reshape((len(tim), len(lat), len(lon)))
+ Rnl = res[:, :, 35].reshape((len(tim), len(lat), len(lon)))
longitude.long_name = 'Longitude'
longitude.units = 'degrees East'
latitude.long_name = 'Latitude'
@@ -385,40 +416,67 @@ def run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
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 qair, qsea, Rl, Rs, Rnl, dtwl
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")
+meth = input("Give prefered method: \n")
+while meth not in ["S80", "S88", "LP82", "YT96", "UA", "LY04", "C30", "C35",
+ "C40", "ecmwf","Beljaars"]:
+ print("method unknown")
+ meth = input("Give prefered method: \n")
+else:
+ meth = meth #[meth]
+ext = meth+"_"
+#------------------------------------------------------------------------------
gustIn = input("Give gustiness option (to use default press enter): \n")
if (gustIn == ''):
gustIn = None
+ ext = ext+'gust_'
else:
gustIn = np.asarray(eval(gustIn), dtype=float)
+ if ((np.all(gustIn) == 0)):
+ ext = ext+'nogust_'
+ else:
+ ext = ext+'gust_'
+#------------------------------------------------------------------------------
cskinIn = input("Give cool skin option (to use default press enter): \n")
if (cskinIn == ''):
cskinIn = None
+ if ((cskinIn == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
+ or meth == "YT96" or meth == "UA" or
+ 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")):
+ cskinIn = 1
+ ext = ext+'skin_'
else:
cskinIn = int(cskinIn)
+ if (cskinIn == 0):
+ ext = ext+'noskin_'
+ elif (cskinIn == 1):
+ ext = ext+'skin_'
+#------------------------------------------------------------------------------
tolIn = input("Give tolerance option (to use default press enter): \n")
if (tolIn == ''):
- tolIn = None
+ tolIn = ['flux', 0.01, 1, 1]
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]
+ext = ext+'tol'+tolIn[0]
+#------------------------------------------------------------------------------
outF = input("Give path and output file name: \n")
if (outF == ''):
- outF = inF[:-3]+"_"+meth+".nc"
+ outF = "out_"+inF[:-3]+"_"+ext+".nc"
+elif (outF[-3:] != '.nc'):
+ outF = outF+".nc"
else:
outF = outF
-
+#------------------------------------------------------------------------------
print("\n run_ASFC.py, started for method "+meth)
res, lon, lat = run_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth)
@@ -440,14 +498,16 @@ if (inF == 'era5_r360x180.nc'):
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title(meth+', '+ttl[i])
- plt.savefig('./'+ttl[i][:3]+'_'+meth+'.png', dpi=300, bbox_inches='tight')
+ plt.savefig('./'+ttl[i][:3]+'_'+ext+'.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')
+ plt.title(meth)
+ plt.xlabel("points")
+ plt.ylabel(ttl[i])
+ plt.savefig('./'+ttl[i][:3]+'_'+ext+'.png', dpi=300, bbox_inches='tight')
#%% generate txt file with statistic
if ((cskinIn == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
@@ -455,11 +515,11 @@ if ((cskinIn == None) and (meth == "S80" or meth == "S88" or meth == "LP82"
meth == "LY04")):
cskinIn = 0
elif ((cskinIn == None) and (meth == "C30" or meth == "C35" or meth == "C40"
- or meth == "ERA5" or meth == "Beljaars")):
+ or meth == "ecmwf" 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")):
+elif (np.all(gustIn == None) and (meth == "UA" or meth == "ecmwf")):
gustIn = [1, 1, 1000]
elif np.all(gustIn == None):
gustIn = [1, 1.2, 800]
@@ -477,7 +537,7 @@ ttl = np.asarray(["tau ", "shf ", "lhf ", "L ", "cd ", "cdn ",
"qsr ", "usr ", "psim ", "psit ", "psiq ", "u10n ",
"t10n ", "tv10n", "q10n ", "zo ", "zot ", "zoq ",
"urefs", "trefs", "qrefs", "itera", "dter ", "dqer ",
- "qair ", "qsea ", "Rl ", "Rs ", "Rnl "])
+ "dtwl ", "qair ", "qsea ", "Rl ", "Rs ", "Rnl "])
header = ["var", "mean", "median", "min", "max", "5%", "95%"]
n = np.shape(res)
stats = np.copy(ttl)
@@ -493,6 +553,7 @@ if (inF == 'era5_r360x180.nc'):
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'))
elif (inF == "data_all.nc"):
stats = np.c_[stats, np.nanmean(res, axis=1)]
stats = np.c_[stats, np.nanmedian(res, axis=1)]
@@ -505,3 +566,8 @@ elif (inF == "data_all.nc"):
"2.2e")), file=open('./stats.txt', 'a'))
print('-'*79+'\n', file=open('./stats.txt', 'a'))
+print('input file name: {}, \n method: {}, \n gustiness: {}, \n cskin: {},'
+ ' \n tolerance: {}, \n output is written in: {}'.format(inF, meth,
+ gustIn, cskinIn,
+ tolIn, outF),
+ file=open('./readme.txt', 'w'))