"AirSeaFluxCode is developed to provide an easy and accessible way to calculate turbulent surface fluxes (TSFs) from a small number of bulk variables and for a viariety of bulk algorithms. \n",
"\n",
"By running AirSeaFluxCode you can compare different bulk algorithms and to also investigate the effect choices within the implementation of each parameterisation have on the TSFs estimates. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Getting started\n",
"\n",
"Let's first import the basic python packages we will need for reading in our input data, to perform basic statistics and plotting"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# first import all packages you might need\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import netCDF4 as nc\n",
"import numpy as np\n",
"import pandas as pd"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exchange coefficients\n",
"We can implement a function to calculate the exchange coefficients and momentum roughness length given the method and the 10 metre neutral wind speed (u10n). Note that these are the same functions included in flux_subs.py (cdn_calc, ctcqn_calc, cdn_from roughness) that are called in AirSeaFluxCode, just adjusted to run with dummy input data generated in the next step."
" elif (meth == \"ecmwf\" or meth == \"Beljaars\"):\n",
" # eq. (3.26) p.38 over sea IFS Documentation cy46r1\n",
" usr = np.sqrt(cdn*np.power(u10n, 2))\n",
" zot = 0.40*1.5e-5/usr\n",
" zoq = 0.62*1.5e-5/usr\n",
" cqn = kappa**2/np.log(10/zo)/np.log(10/zoq)\n",
" ctn = kappa**2/np.log(10/zo)/np.log(10/zot)\n",
" else:\n",
" print(\"unknown method ctcqn: \"+meth)\n",
" return cdn, ctn, cqn, zo\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then, generate \"dummy\" values for u10n to use as input to blk_neutral."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"u10n = np.zeros(1201)\n",
"for jw in range(2, 1201):\n",
" dw = 0.25*(60/(1201-1))\n",
" if (u10n[jw-1] >= 5):\n",
" dw = 60/(1201-1)\n",
" if (u10n[jw-1] >= 30):\n",
" dw = 2*(60/(1201-1))\n",
" u10n[jw] = u10n[jw-1]+dw"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, you can run blk_neutral for any method (\"S80\", \"S88\", \"LP82(zol<=0)\", \"LP82(zol>0)\", \"YT96\", \"UA\", \"LY04(zol<=0)\", \"LY04(zol>0)\", \"C30\", \"C35\", \"C40\", \"ecmwf\", \"Beljaars\") and with the \"dummy\" u10n we generated.\n",
"AirSeaFluxCode is set up to run in its default setting with a minimum number of input variables, namely wind speed; air temperature; and sea surface temperature. Let's load the code, import some real data composed for testing it (Research vessel data) and run AirSeaFluxCode with default settings (latitude set to 45°N, relative humidity 80%, atmospheric pressure 1013hPa, sensor height 18m, output height 10m, cool skin/warm layer corrections turned off, bulk algorithm Smith 1988, gustiness on, saturation vapour pressure following Buck, 2012, tolerance limits set for flux estimates and height adjusted variables (['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1]), number of iterations are ten, non converged points are set to missing and Monin-Obukhov stability length is calculated following the ECMWF implementation."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from AirSeaFluxCode import AirSeaFluxCode"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"inDt = pd.read_csv(\"data_all.csv\")\n",
"date = np.asarray(inDt[\"Date\"])\n",
"lon = np.asarray(inDt[\"Longitude\"])\n",
"lat = np.asarray(inDt[\"Latitude\"])\n",
"spd = np.asarray(inDt[\"Wind speed\"])\n",
"t = np.asarray(inDt[\"Air temperature\"])\n",
"sst = np.asarray(inDt[\"SST\"])\n",
"rh = np.asarray(inDt[\"RH\"])\n",
"p = np.asarray(inDt[\"P\"])\n",
"sw = np.asarray(inDt[\"Rs\"])\n",
"hu = np.asarray(inDt[\"zu\"])\n",
"ht = np.asarray(inDt[\"zt\"])\n",
"hin = np.array([hu, ht, ht])\n",
"del hu, ht, inDt\n",
"# run AirSeaFluxCode\n",
"res = AirSeaFluxCode(spd, t, sst)\n",
"flg = res[\"flag\"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"res is the output of AirSeaFluxCode which is a dataFrame with keys: \"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\", \"iteration\", \"dter\", \"dqer\", \"dtwl\", \"qair\", \"qsea\", \"Rl\", \"Rs\", \"Rnl\", \"ug\", \"Rib\", \"rh\". Let's plot the flux estimates."
"Now we can use the real input we have for latitude, sensor heights, relative humidity, air pressure and shortwave radiation. We will also use the cool skin correction option for bulk algorithm Coare 3.5 (C35) and run 30 iterations."
