toy_ASFC.py

"""
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
        temp = 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)
        res = temp.loc[:,"tau":"Rnl"]
        res = res.to_numpy().T
        del temp
        #%% 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)):
            temp = 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")
            a = temp.loc[:,"tau":"Rnl"]
            a = a.to_numpy()
            del temp
            res[x, :, :] = a
            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, psiq, 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'))