Skip to content

GitLab

Commit 35957d96 authored by sbiri's avatar sbiri
Browse files

Update data_all_stats.txt, data_all.csv, data_all_out.csv, toy_ASFC.py, Documentation.pdf files

parent 7ced8a70
Hide whitespace changes
Inline Side-by-side
Showing with 7086 additions and 0 deletions
Documentation.pdf 0 → 100644
View file @ 35957d96
File added
data_all.csv 0 → 100644
View file @ 35957d96
This source diff could not be displayed because it is too large. You can view the blob instead.
data_all_out.csv 0 → 100644
View file @ 35957d96
This source diff could not be displayed because it is too large. You can view the blob instead.
data_all_stats.txt 0 → 100644
View file @ 35957d96
Input summary
input file name: data_all.csv,
method: UA,
gustiness: [1, 1, 1000],
cskin: 0,
tolerance: ['flux', 0.001, 0.1, 0.1]
| var | mean | median | min | max | 5% | 95% |
|-------|-----------|-----------|-----------|----------|-----------|----------|
| tau | 7.24e-02 | 4.95e-02 | 4.96e-04 | 5.64e-01 | 7.33e-03 | 2.12e-01 |
| shf | 6.80e+00 | 6.16e+00 | -4.06e+01 | 5.56e+01 | -9.07e+00 | 2.52e+01 |
| lhf | 8.26e+01 | 6.71e+01 | -4.15e+01 | 5.46e+02 | -2.89e-01 | 2.13e+02 |
| L | 4.11e+02 | -4.40e+01 | -1.62e+05 | 1.20e+06 | -5.15e+02 | 3.19e+02 |
| cd | 1.15e-03 | 1.16e-03 | 3.44e-04 | 2.02e-03 | 8.74e-04 | 1.37e-03 |
| cdn | 1.15e-03 | 1.12e-03 | 9.70e-04 | 1.71e-03 | 9.86e-04 | 1.41e-03 |
| ct | 1.18e-03 | 1.18e-03 | 3.60e-04 | 2.87e-03 | 8.90e-04 | 1.46e-03 |
| ctn | 1.12e-03 | 1.12e-03 | 1.04e-03 | 1.25e-03 | 1.05e-03 | 1.21e-03 |
| cq | 1.18e-03 | 1.18e-03 | 3.60e-04 | 2.87e-03 | 8.90e-04 | 1.46e-03 |
| cqn | 1.12e-03 | 1.12e-03 | 1.04e-03 | 1.25e-03 | 1.05e-03 | 1.21e-03 |
| tsrv | -5.18e-02 | -5.57e-02 | -2.34e-01 | 8.52e-02 | -1.27e-01 | 3.28e-02 |
| tsr | -2.70e-02 | -2.80e-02 | -1.67e-01 | 8.30e-02 | -8.42e-02 | 3.53e-02 |
| qsr | -1.37e-04 | -1.24e-04 | -5.57e-04 | 4.67e-05 | -3.24e-04 | 6.52e-07 |
| usr | 2.20e-01 | 2.04e-01 | 2.02e-02 | 6.92e-01 | 7.95e-02 | 4.20e-01 |
| psim | 4.77e-01 | 4.56e-01 | -8.86e+00 | 4.10e+00 | -7.94e-01 | 1.75e+00 |
| psit | 9.03e-01 | 8.44e-01 | -8.86e+00 | 5.74e+00 | -8.02e-01 | 2.75e+00 |
| psiq | 9.03e-01 | 8.44e-01 | -8.86e+00 | 5.74e+00 | -8.02e-01 | 2.75e+00 |
| u10n | 6.37e+00 | 6.10e+00 | 1.17e+00 | 1.67e+01 | 2.56e+00 | 1.12e+01 |
| t10n | 1.79e+01 | 1.82e+01 | -2.76e+00 | 2.98e+01 | 2.82e+00 | 2.79e+01 |
| tv10n | 2.93e+02 | 2.93e+02 | 2.71e+02 | 3.08e+02 | 2.77e+02 | 3.04e+02 |
| q10n | 1.06e-02 | 9.67e-03 | -1.49e-03 | 2.52e-02 | 3.58e-03 | 1.78e-02 |
| zo | 8.96e-05 | 6.30e-05 | 2.65e-05 | 6.37e-04 | 2.94e-05 | 2.39e-04 |
| zot | 6.01e-05 | 6.18e-05 | 1.63e-05 | 1.23e-04 | 4.09e-05 | 7.13e-05 |
| zoq | 6.01e-05 | 6.18e-05 | 1.63e-05 | 1.23e-04 | 4.09e-05 | 7.13e-05 |
| urefs | 6.18e+00 | 5.93e+00 | 6.25e-01 | 1.65e+01 | 2.30e+00 | 1.11e+01 |
| trefs | 1.80e+01 | 1.84e+01 | -2.75e+00 | 2.98e+01 | 2.87e+00 | 2.80e+01 |
| qrefs | 1.10e-02 | 1.00e-02 | 7.66e-04 | 2.53e-02 | 3.66e-03 | 1.83e-02 |
| itera | 3.41e+00 | 3.00e+00 | 2.00e+00 | 6.00e+00 | 3.00e+00 | 4.00e+00 |
| dter | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 |
| dqer | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 |
| dtwl | 3.00e-01 | 3.00e-01 | 3.00e-01 | 3.00e-01 | 3.00e-01 | 3.00e-01 |
| qair | 1.09e-02 | 9.90e-03 | 7.43e-04 | 2.53e-02 | 3.64e-03 | 1.82e-02 |
| qsea | 1.47e-02 | 1.39e-02 | 3.26e-03 | 2.61e-02 | 4.76e-03 | 2.46e-02 |
| Rl | 3.70e+02 | 3.70e+02 | 3.70e+02 | 3.70e+02 | 3.70e+02 | 3.70e+02 |
| Rs | 2.61e+02 | 2.33e+02 | 1.67e-02 | 9.71e+02 | 2.84e+01 | 5.37e+02 |
| Rnl | 4.25e+01 | 4.50e+01 | -6.02e+01 | 1.06e+02 | -3.65e+01 | 9.96e+01 |
-------------------------------------------------------------------------------
toy_ASFC.py 0 → 100644
View file @ 35957d96
"""
example of running AirSeaFluxCode with
1. R/V data (data_all.csv) or
2. one day era5 hourly data (era5_r360x180.nc)
compute fluxes
output NetCDF4
and statistics in stats.txt
@author: sbiri
"""
#%% import packages
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import pandas as pd
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 toy_ASFC(inF, outF, gustIn, cskinIn, tolIn, meth):
"""
Parameters
----------
inF : str
input filename either data_all.csv 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.csv"):
#%% load data_all
inDt = pd.read_csv("data_all.csv")
date = np.asarray(inDt["Date"])
lon = np.asarray(inDt["Longitude"])
lat = np.asarray(inDt["Latitude"])
spd = np.asarray(inDt["Wind speed"])
t = np.asarray(inDt["Air temperature"])
sst = np.asarray(inDt["SST"])
rh = np.asarray(inDt["RH"])
p = np.asarray(inDt["P"])
sw = np.asarray(inDt["Rs"])
hu = np.asarray(inDt["zu"])
ht = np.asarray(inDt["zt"])
hin = np.array([hu, ht, ht])
del hu, ht, inDt
#%% 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, L="ecmwf", n=30)
#%% delete variables
del spd, t, sst, rh, p, sw, 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))
date = np.copy(tim)
#%% run AirSeaFluxCode
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, :],
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="ecmwf")
res[x, :, :] = a.T
del a
if (outF[-3:] == '.nc'):
if (inF == 'era5_r360x180.nc'):
#%% 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'))
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'))
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
dtwl[:] = 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'
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, dtwl
del tim, T, Td, p, lw, sw, lsm, spd, hin, latIn
else:
#%% 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')
dtwl = fid.createVariable('dter', '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]
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'
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
else:
#%% save as .csv
np.savetxt(outF, np.vstack((date, lon, lat, res)).T,
delimiter=',',
header="date, lon, lat, tau, shf, lhf, L, cd, cdn, ct, ctn,"
" cq, cqn, tsrv, tsr, qsr, usr, psim, psit, u10n, t10n,"
" tv10n, q10n, zo, zot, zoq, uref, tref, qref, iter, dter,"
" dqer, dtwl, qair, qsea, Rl, Rs, Rnl")
return res, lon, lat
#%% run function
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"]:
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 = ['flux', 1e-3, 0.1, 0.1]
else:
tolIn = eval(tolIn)
ext = ext+'tol'+tolIn[0]
#------------------------------------------------------------------------------
outF = input("Give path and output file name: \n")
if ((outF == '') and (inF == "data_all.csv")):
outF = "out_"+inF[:-4]+"_"+ext+".csv"
elif ((outF == '') and (inF == "era5_r360x180.nc")):
outF = "out_"+inF[:-3]+"_"+ext+".nc"
elif ((outF[-4:] == '.csv') and (inF == 'era5_r360x180.nc')):
outF = outF[:-4]+".nc"
elif ((outF[-3:] != '.nc') and (outF[-4:] != '.csv')):
if (inF == "data_all.csv"):
outF = outF+".csv"
else:
outF = outF+".nc"
else:
outF = outF
#------------------------------------------------------------------------------
print("\n run_ASFC.py, started for method "+meth)
res, lon, lat = toy_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]+'_'+ext+'.png', dpi=300, bbox_inches='tight')
elif (inF == "data_all.csv"):
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)
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"
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 == "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 == "ecmwf")):
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 ",
"dtwl ", "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'))
print('-'*79+'\n', file=open('./stats.txt', 'a'))
elif (inF == "data_all.csv"):
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'))
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'))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment