using Pkg
Pkg.add("IJulia")
Pkg.add("PyPlot")
Pkg.add("PhysOcean")
# to install the developpement version use this command
#Pkg.clone("https://github.com/gher-ulg/PhysOcean.jl")
using IJulia
notebook()
This opens the web-browser with a Jupyter session or adds a new tab into an existing browser window
Navigate through the directories and open the file Ocean Heat Fluxes.ipynb
Jupyter supports Chrome, Safari and Firefox. Internet Explorer is not recommended
(http://jupyter-notebook.readthedocs.io/en/latest/notebook.html#browser-compatibility)
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.
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.
#uncomment the two first lines if PhysOcean is not installed in the virtual machine
#using Pkg
#Pkg.add("PhysOcean")
using PyPlot
using PhysOcean
using DelimitedFiles
using Statistics
using Dates
#We are now going to download all necessary files from dox
filename = "solrad.csv";
if !isfile(filename)
@info("downloading $filename")
cp(download("https://dox.uliege.be/index.php/s/h5vfLDkOQpMQdnB/download"),filename)
else
@info("$filename is already downloaded")
end
filename = "albedoECMWF.txt";
if !isfile(filename)
@info("downloading $filename")
cp(download("https://dox.uliege.be/index.php/s/Syn1221wZ8eR9KF/download"),filename)
else
@info("$filename is already downloaded")
end
filename = "airtemp.csv";
if !isfile(filename)
@info("downloading $filename")
cp(download("https://dox.uliege.be/index.php/s/JhQCD5lDpD30gqP/download"),filename)
else
@info("$filename is already downloaded")
end
filename = "tccECMWF.txt";
if !isfile(filename)
@info("downloading $filename")
cp(download("https://dox.uliege.be/index.php/s/tIQe4sYftzP3AsZ/download"),filename)
else
@info("$filename is already downloaded")
end
filename = "rhECMWF.txt";
if !isfile(filename)
@info("downloading $filename")
cp(download("https://dox.uliege.be/index.php/s/gXELajKkTagXFJL/download"),filename)
else
@info("$filename is already downloaded")
end
filename = "watertemp_calvi.txt";
if !isfile(filename)
@info("downloading $filename")
cp(download("https://dox.uliege.be/index.php/s/UqUvDoqVFrTXf01/download"),filename)
else
@info("$filename is already downloaded")
end
filename = "windspeed.csv";
if !isfile(filename)
@info("downloading $filename")
cp(download("https://dox.uliege.be/index.php/s/NVm4fabFwSId8vf/download"),filename)
else
@info("$filename is already downloaded")
end
filename = "atmpress.csv";
if !isfile(filename)
@info("downloading $filename")
cp(download("https://dox.uliege.be/index.php/s/028JlDsTYxheOti/download"),filename)
else
@info("$filename is already downloaded")
end
In order to calculate the solar flux, we need to know the incoming radiation (measured at Stareso) and the albedo (from teh model ECMWF).
Q = readdlm("solrad.csv");#(W/m²)
@show size(Q)
Q = Q[:,1];
# set all negative values to zero
Q[Q .< 0] .= 0;
data_albedo = readdlm("albedoECMWF.txt");
al = data_albedo[:,1];
date = datetime_matlab.(data_albedo[:,2]);
Qs = solarflux.(Q,al);
size(Qs)
@show date[1:10]
plot(date,Qs);
ylabel("W m⁻²");
xlabel("time");
title("Solar flux")
$\rightarrow$ What do you observe? How is the solar flux changing with time?
$\rightarrow$ Compare the solar flux in a winter day with one in a summer day (make a new plot below)
#Make your plot in this cell comparing the solar flux magnitude in winter and in summer
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.
Ta = readdlm("airtemp.csv");
Ta = Ta[:,1];
tcc = readdlm("tccECMWF.txt");
tcc = tcc[:,1];
#tcc = (tcc*0)+1;
r = readdlm("rhECMWF.txt");
r = r[:,1];
ep = vaporpressure.(Ta);#Ta in C
ep = ep.*r;
Ts = readdlm("watertemp_calvi.txt");
#Ts = Ts[1:30915,1];
Ts = Ts[:,1];
Ts[Ts.==9999].=NaN;
@show size(Ts) size(ep) size(r) size(tcc) size(Ta)
Qb = longwaveflux.(Ts,Ta,ep,tcc);#Ta in C
plot(date,Qb);
ylabel("W m⁻²");
xlabel("time");
title("Longwave flux")
What do you think that the straight line at the end of 2016 is? Chech the input variables to see where this comes from
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)
w = readdlm("windspeed.csv");
w = w[:,1];
@show size(w)
#Ta = readdlm("airtemp.csv");
#Ta = Ta[:,1];
#@show size(Ts)
#Ts = readdlm("watertemp_calvi.csv");
Qc = sensibleflux.(Ts,Ta,w);
plot(date,Qc);
ylabel("W m⁻²");
xlabel("time");
title("Sensible heat flux")
pa = readdlm("atmpress.csv");
pa = pa[:,1];
Qe = latentflux.(Ts,Ta,r,w,pa);
plot(date,Qe);
ylabel("W m⁻²");
xlabel("time");
title("Latent heat flux")
#Here a little exercise: we visualise wind speed, and a filtered version of the wind speed to see long-term changes
# The function gaussfilter (xd = gaussfilter(x,WW) ) performs a Gaussian filtering with a window width of WW
#Parameter WW in the units of the input dataset
wd = gaussfilter(w,3*30*24); #20-minute time step * 3 = 1hour; *24 = 1 day * 30 = 1 month filter length
plot(date,w,label ="Wind")
#hold(true)
plot(date,wd,"r",label ="1-month filter");
ylabel("m s⁻¹");
xlabel("time");
title("Wind")
We are going to average the heat fluxes time series. In this example we do daily averages.
date = readdlm("rhECMWF.txt");#(W/m²)
date = datetime_matlab.(date[:,2]);
@show date[1:10]
#We look for all indices of time in a given day, and average all values at those indices
s=1;
floor_date = Dates.Date.(date);
#floor_date = [floor(d,Dates.Day) for d in date]
date_daily = minimum(floor_date):Dates.Day(1):maximum(floor_date)
ndates = length(date_daily)
Qsd = zeros(ndates)
Qbd = zeros(ndates)
Qcd = zeros(ndates)
Qed = zeros(ndates)
for i = date_daily
Qsd[s] = mean(Qs[floor_date .== i]);
Qbd[s] = mean(Qb[floor_date .== i]);
Qcd[s] = mean(Qc[floor_date .== i]);
Qed[s] = mean(Qe[floor_date .== i]);
s=s+1;
end
@show date_daily[1:10]
plot(date,Qb,label ="Longwave flux")
#hold(true)
plot(date_daily,Qbd,"r",label ="Daily average");
ylabel("W m⁻²");
xlabel("time");
legend();
$\rightarrow$ Try visualizing other components of the heat fluxes and their filtered version
Same exercice, but with monthly averages
m = Dates.month.(date)
dates_monthly = DateTime(2016,1,15):Dates.Month(1):DateTime(2016,12,15);
Qsm = zeros(12)
Qbm = zeros(12)
Qcm = zeros(12)
Qem = zeros(12)
for i = 1:12
Qsm[i] = nanmean(Qs[m .== i]);
Qbm[i] = nanmean(Qb[m .== i]);
Qcm[i] = nanmean(Qc[m .== i]);
Qem[i] = nanmean(Qe[m .== i]);
end
plot(date,Qb,label ="Longwave flux")
#hold(true)
plot(dates_monthly,Qbm,"r",label ="Montly average");
ylabel("W m⁻²");
xlabel("time");
legend();
$\rightarrow$ Try this monthly average for other components of the heat fluxes
Now that we have the 4 components of the air-sea heat fluxes, we can calculate the total heat flux
Qi = Qs -(Qb+Qe+Qc);
low = 1*30*24; #what is the length of this filter in days?
Qif = gaussfilter(Qi,low);
Qsf = gaussfilter(Qs,low);
Qbf = gaussfilter(Qb,low);
Qcf = gaussfilter(Qc,low);
Qef = gaussfilter(Qe,low);
@show low
# plot of the raw fluxes (with no filtering)
figure()
plot(date,Qif,label="Qi")
# datetick("x")
grid("on")
#hold(true)
plot(date,Qs,"r")
plot(date,Qb,"g")
plot(date,Qc,"m")
plot(date,Qe,"k")
legend(["Qi","Qs","Qb","Qc","Qe"]);
ylabel("W m⁻²");
xlabel("time");
title("Total fluxes")
$\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