From 66739e9bfdd84eb8359c5fbc8c161685f765d257 Mon Sep 17 00:00:00 2001
From: sbiri <sbiri@noc.ac.uk>
Date: Mon, 8 Feb 2021 10:08:20 +0000
Subject: [PATCH] Update run_ASFC.py
---
run_ASFC.py | 134 +++++++++++++++++++++++++++++++++++++++-------------
1 file changed, 100 insertions(+), 34 deletions(-)
diff --git a/run_ASFC.py b/run_ASFC.py
index 1b87930..1936840 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'))
--
GitLab