OceanHeatFluxes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Installing Julia and Jupyter\n",
    "\n",
    "\n",
    "* Copy the data files and the ipynb file on your computer (and remember where you put them ☺)\n",
    "* Download and install Julia from http://julialang.org/downloads/\n",
    "* Open Julia and enter the following commands:\n",
    "\n",
    "```julia\n",
    "using Pkg\n",
    "Pkg.add(\"IJulia\")\n",
    "Pkg.add(\"PyPlot\")\n",
    "Pkg.add(\"PhysOcean\")\n",
    "# to install the developpement version use this command\n",
    "#Pkg.clone(\"https://github.com/gher-ulg/PhysOcean.jl\")\n",
    "```\n",
    "\n",
    "* Note that some commands may take some time to execute.\n",
    "* Load the IJulia package and start the Jupyter notebook interface\n",
    "\n",
    "```julia\n",
    "using IJulia\n",
    "notebook()\n",
    "```\n",
    "\n",
    "* This opens the web-browser with a Jupyter session or adds a new tab into an existing browser window\n",
    "\n",
    "* Navigate through the directories and open the file `Ocean Heat Fluxes.ipynb`\n",
    "\n",
    "* Jupyter supports Chrome, Safari and Firefox. Internet Explorer is not recommended\n",
    "(http://jupyter-notebook.readthedocs.io/en/latest/notebook.html#browser-compatibility)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise on Air-Sea fluxes\n",
    "\n",
    "We are going to calculate the air-sea fluxes on the Bay of Calvi during one year (2016) using a combination of in situ data (from the meteorological station in Stareso) and modeled data (when no in situ data are available). The names of the files read indicate if they come from a model (ECMWF) or not.\n",
    "\n",
    "If working on the Virtual Macine, you need to install one package (the Pkg.add(\"PhysOcean\") line in the cell below) before you can call it.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#uncomment the two first lines if PhysOcean is not installed in the virtual machine\n",
    "#using Pkg\n",
    "#Pkg.add(\"PhysOcean\")\n",
    "using PyPlot\n",
    "using PhysOcean\n",
    "using DelimitedFiles\n",
    "using Statistics\n",
    "using Dates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We are now going to download all necessary files from dox\n",
    "\n",
    "filename = \"solrad.csv\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/h5vfLDkOQpMQdnB/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded\")\n",
    "end\n",
    "\n",
    "filename = \"albedoECMWF.txt\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/Syn1221wZ8eR9KF/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded\")\n",
    "end\n",
    "\n",
    "filename = \"airtemp.csv\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/JhQCD5lDpD30gqP/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded\")\n",
    "end\n",
    "\n",
    "filename = \"tccECMWF.txt\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/tIQe4sYftzP3AsZ/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded\")\n",
    "end\n",
    "\n",
    "filename = \"rhECMWF.txt\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/gXELajKkTagXFJL/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded\")\n",
    "end\n",
    "\n",
    "filename = \"watertemp_calvi.txt\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/UqUvDoqVFrTXf01/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded\")\n",
    "end\n",
    "\n",
    "filename = \"windspeed.csv\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/NVm4fabFwSId8vf/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded\")\n",
    "end\n",
    "\n",
    "filename = \"atmpress.csv\";\n",
    "if !isfile(filename)\n",
    "    @info(\"downloading $filename\")\n",
    "    cp(download(\"https://dox.uliege.be/index.php/s/028JlDsTYxheOti/download\"),filename)\n",
    "else\n",
    "    @info(\"$filename is already downloaded\")\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solar flux\n",
    "\n",
    "In order to calculate the solar flux, we need to know the incoming radiation (measured at Stareso) and the albedo (from teh model ECMWF)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = readdlm(\"solrad.csv\");#(W/m²)\n",
    "@show size(Q)\n",
    "Q = Q[:,1];\n",
    "\n",
    "# set all negative values to zero \n",
    "Q[Q .< 0] .= 0;\n",
    "\n",
    "data_albedo = readdlm(\"albedoECMWF.txt\");\n",
    "al = data_albedo[:,1];\n",
    "date = datetime_matlab.(data_albedo[:,2]);\n",
    "\n",
    "Qs = solarflux.(Q,al);\n",
    "size(Qs)\n",
    "@show date[1:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(date,Qs);\n",
    "ylabel(\"W m⁻²\");\n",
    "xlabel(\"time\");\n",
    "title(\"Solar flux\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\rightarrow$ What do you observe? How is the solar flux changing with time?\n",
    "\n",
    "$\\rightarrow$ Compare the solar flux in a winter day with one in a summer day (make a new plot below)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Make your plot in this cell comparing the solar flux magnitude in winter and in summer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Longwave flux\n",
    "\n",
    "We now calculate the longwave flux. If you recall from the equation given in class, we will need the air and water temperature, the total cloud cover and the relative humidity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ta = readdlm(\"airtemp.csv\");\n",
    "Ta = Ta[:,1];\n",
    "\n",
    "tcc = readdlm(\"tccECMWF.txt\");\n",
    "tcc = tcc[:,1];\n",
    "#tcc = (tcc*0)+1;\n",
    "\n",
    "r = readdlm(\"rhECMWF.txt\");\n",
    "r = r[:,1];\n",
    "\n",
    "ep = vaporpressure.(Ta);#Ta in C\n",
    "ep = ep.*r;\n",
    "\n",
    "Ts = readdlm(\"watertemp_calvi.txt\");\n",
    "#Ts = Ts[1:30915,1];\n",
    "Ts = Ts[:,1];\n",
    "Ts[Ts.==9999].=NaN;\n",
    "@show size(Ts) size(ep) size(r) size(tcc) size(Ta)\n",
    "Qb = longwaveflux.(Ts,Ta,ep,tcc);#Ta in C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot(date,Qb);\n",
    "ylabel(\"W m⁻²\");\n",
    "xlabel(\"time\");\n",
    "title(\"Longwave flux\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do you think that the straight line at the end of 2016 is? Chech the input variables to see where this comes from"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sensible heat flux\n",
    "\n",
    "Which variables do we need to calculate the sensible heat flux? (some have already been read above we do not need to read them again)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = readdlm(\"windspeed.csv\");\n",
    "w = w[:,1];\n",
    "@show size(w)\n",
    "#Ta = readdlm(\"airtemp.csv\");\n",
    "#Ta = Ta[:,1];\n",
    "#@show size(Ts)\n",
    "#Ts = readdlm(\"watertemp_calvi.csv\");\n",
    "\n",
    "Qc = sensibleflux.(Ts,Ta,w);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot(date,Qc);\n",
    "ylabel(\"W m⁻²\");\n",
    "xlabel(\"time\");\n",
    "title(\"Sensible heat flux\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Latent heat flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pa = readdlm(\"atmpress.csv\");\n",
    "pa = pa[:,1];\n",
    "\n",
    "Qe = latentflux.(Ts,Ta,r,w,pa); "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(date,Qe);\n",
    "ylabel(\"W m⁻²\");\n",
    "xlabel(\"time\");\n",
    "title(\"Latent heat flux\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Here a little exercise: we visualise wind speed, and a filtered version of the wind speed to see long-term changes\n",
    "# The function gaussfilter (xd = gaussfilter(x,WW) ) performs a Gaussian filtering with a window width of WW\n",
    "#Parameter WW in the units of the input dataset\n",
    "\n",
    "wd = gaussfilter(w,3*30*24); #20-minute time step * 3 = 1hour; *24 = 1 day * 30 = 1 month filter length\n",
    "plot(date,w,label =\"Wind\")\n",
    "#hold(true)\n",
    "plot(date,wd,\"r\",label =\"1-month filter\");\n",
    "ylabel(\"m s⁻¹\");\n",
    "xlabel(\"time\");\n",
    "title(\"Wind\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Daily averages\n",
    "\n",
    "We are going to average the heat fluxes time series. In this example we do daily averages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "date = readdlm(\"rhECMWF.txt\");#(W/m²)\n",
    "date = datetime_matlab.(date[:,2]);\n",
    "@show date[1:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We look for all indices of time in a given day, and average all values at those indices \n",
    "s=1;\n",
    "\n",
    "floor_date = Dates.Date.(date);\n",
    "#floor_date = [floor(d,Dates.Day) for d in date]\n",
    "date_daily = minimum(floor_date):Dates.Day(1):maximum(floor_date)\n",
    "\n",
    "ndates = length(date_daily)\n",
    "Qsd = zeros(ndates)\n",
    "Qbd = zeros(ndates)\n",
    "Qcd = zeros(ndates)\n",
    "Qed = zeros(ndates)\n",
    "\n",
    "for i = date_daily\n",
    "    Qsd[s] = mean(Qs[floor_date .== i]);\n",
    "    Qbd[s] = mean(Qb[floor_date .== i]);\n",
    "    Qcd[s] = mean(Qc[floor_date .== i]);\n",
    "    Qed[s] = mean(Qe[floor_date .== i]);\n",
    "    s=s+1;\n",
    "end\n",
    "@show date_daily[1:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(date,Qb,label =\"Longwave flux\")\n",
    "#hold(true)\n",
    "plot(date_daily,Qbd,\"r\",label =\"Daily average\");\n",
    "ylabel(\"W m⁻²\");\n",
    "xlabel(\"time\");\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\rightarrow$ Try visualizing other components of the heat fluxes and their filtered version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monthly averages\n",
    "\n",
    "Same exercice, but with monthly averages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Dates.month.(date)\n",
    "dates_monthly = DateTime(2016,1,15):Dates.Month(1):DateTime(2016,12,15);\n",
    "\n",
    "Qsm = zeros(12)\n",
    "Qbm = zeros(12)\n",
    "Qcm = zeros(12)\n",
    "Qem = zeros(12)\n",
    "\n",
    "for i = 1:12\n",
    "    Qsm[i] = nanmean(Qs[m .== i]);\n",
    "    Qbm[i] = nanmean(Qb[m .== i]);\n",
    "    Qcm[i] = nanmean(Qc[m .== i]);\n",
    "    Qem[i] = nanmean(Qe[m .== i]);\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(date,Qb,label =\"Longwave flux\")\n",
    "#hold(true)\n",
    "plot(dates_monthly,Qbm,\"r\",label =\"Montly average\");\n",
    "ylabel(\"W m⁻²\");\n",
    "xlabel(\"time\");\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\rightarrow$ Try this monthly average for other components of the heat fluxes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total flux\n",
    "\n",
    "Now that we have the 4 components of the air-sea heat fluxes, we can calculate the total heat flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Qi = Qs -(Qb+Qe+Qc);\n",
    "low = 1*30*24; #what is the length of this filter in days? \n",
    "\n",
    "Qif = gaussfilter(Qi,low);\n",
    "Qsf = gaussfilter(Qs,low);\n",
    "Qbf = gaussfilter(Qb,low);\n",
    "Qcf = gaussfilter(Qc,low);\n",
    "Qef = gaussfilter(Qe,low);\n",
    "@show low"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# plot of the raw fluxes (with no filtering)\n",
    "figure()\n",
    "plot(date,Qif,label=\"Qi\")\n",
    "# datetick(\"x\")\n",
    "grid(\"on\")\n",
    "#hold(true)\n",
    "plot(date,Qs,\"r\")\n",
    "plot(date,Qb,\"g\")\n",
    "plot(date,Qc,\"m\")\n",
    "plot(date,Qe,\"k\")\n",
    "legend([\"Qi\",\"Qs\",\"Qb\",\"Qc\",\"Qe\"]);\n",
    "ylabel(\"W m⁻²\");\n",
    "xlabel(\"time\");\n",
    "title(\"Total fluxes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\rightarrow$ Now make the same plot but with different filter lengths, to see the effect of the filtering in the data, and also to better see the different curves. Do a plot with 3-month filtering and 6-month filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.5.2",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}