%****************************************************** Sonic anemometer ********************************************************
% This program reads sonic anemometer raw data and computes 1-hour averaged turbulent fluxes and other statistics
% The program is adopted for the data collected at Tiksi site (Russia, Siberia) during 2012 (NOAA SEARCH program)
% This sub-program is launched by 'read_Tiksi_hr_2012.m' or 'read_Tiksi_day_2012.m'
% Written by Andrey Grachev <Andrey.Grachev@noaa.gov>, April-May 2013
% The program is based on the NOAA ETL/PSD3 programs 'read_sonic_hr_10.m', 'read_Eureka_sonic_0_hr_2009_egu.m', read_sonic_20Hz_05hr_Materhorn2012_es2_level_1_v2.m etc
% Original program 'read_sonic_hr' was written by Chris Fairall. Help from Ludovic Bariteau is appreciated 
%********************************************************************************************************************************
% New for version 2: 
% Computation of epsilom, Nt, Monin-Obukhov universal functions (phi_w, phi_epsilon etc.) were added
% These results were saved for quick analysis
%********************************************************************************************************************************

disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>');
disp('read_sonic_10Hz_1hr_Tiksi_2012_9m_v2.m')
disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>');

warning off all;

%*************************************************** The Tiksi Observatory ******************************************************

% The Tiksi Observatory (Russia) located at the mouth of the Lena River Delta on the Northern edge of the Eurasian continent (71.6 N, 128.9 E)
% Time difference between Tiksi/Yakutsk and UTC/GMT is +10 hr (May 2013)
% The Tiksi weather station was established at the Polyarka settlement on August 12, 1932 by the chief management of the northern sea route. 
% The centerpiece of the Tiksi Observatory currently is a 20-m tower. Three levels of the tower (3.5, 9.5, 15.5 m) are instrumented with three-axis sonic anemometer/thermometers
% A fast-response infrared open-path gas analyzer Licor-7500 (H2O and CO2 measurements) is also mounted on the second level. 
% Mean temperature and humidity are measured at three levels and can used to evaluate the vertical T/RH gradients. Shortwave and longwave radiation are also measured.

%************************************************ Sonic Anemometer Information **************************************************
 
% Model: Applied Technologies, Inc. (ATI), three-axis sonic anemometer/thermometers (K-probe)
% Sampling frequency: 10 Hertz. Raw data file: 1-hr long
z_level_n = 3;     % Height of measurements = 9.5 m (only one level of the flux measurements in 2012)
sonic_9m_dir = 240;  % Orientation (azimuth) of the sonic anemoneter relatively true North
% Here orientation of the sonic anemoneter (sonic_9m_dir = 197 deg = 180 + 17 deg) is provided by Irina Repina
% NOAA ATI sonic anemometer is approximately pointed out to South-South-West (opposite to the Arctic Ocean)
% Axis 'x' (u-wind component) is aligned with the ATI sonic anemometer array
% Positive u-component represents wind into the array, parallel to the support boom (toward the tower), i.e. wind from the un-obstructed direction. 
% If you're looking into the array, along the boom, toward the tower, then positive v-wind component is to your left (from right to left). 
% These are right handed with an upward (positive) w-wind component

% Sonic anemometer orientation was again measured by Matt Okraszewski and Vasily Kustov in June 2013 during visit in Tiksi 
% Orientation was 180 deg (true South) upon arrival in Tiksi (supposedly because earlier adjusment by Arsenii)
% Orientation is set up 190 deg (45 deg from corner) on June 17, 2013
% This value is close to 197 deg found previously by Irina Repina

samp_freq = 20;        % sonic sampling rate in Hertz (10 Hz)
npos = 3600*samp_freq; % most number of possible data points (36000)
sfreq = 1/samp_freq;
nx = samp_freq*(3600-1);
  
% Sonic output data have the following format:
% 1 column: time, 1 hr length (maximum), sampling frequency is 10 Hz.
% 2 column: u, divide by 100 to put in units of m/s
% 3 column: v, divide by 100 to put in units of m/s
% 4 column: w, divide by 100 to put in units of m/s
% 5 column: T, divide by 100 to put in units of deg C
% Here u, v , and w are the longitudinal, lateral, and vertical velocity components, T is the air temperature (or "sonic temperature").

T    = zeros(1,npos); 
ur   = zeros(1,npos);
vr   = zeros(1,npos);
wr   = zeros(1,npos);
Tr   = zeros(1,npos);
stux = zeros(5,npos);

%************************************************** Open a Sonic Data File ******************************************************

% Read raw data from a sonic marked as 'tik1' (NOAA server):
e = ['X:\tiksi\FLUXTOWER\raw\sonic\9meter\2012\tik112' ddd hh '_raw.txt']; % data path for Year 2012  
% Read raw data from a sonic marked as 'tik1' (C:\ 'Hard drive'):
% e = ['C:\Data\AG-Raw_Data\Arctic\Tiksi\Fluxtower\Raw\Sonic\9meter\2012\tik112' ddd hh '_raw.txt']; % data path for Year 2012  
flist = fopen(e, 'r');

% If the data file is not exist --> Replace with NaNs
if flist > 0  % if the file exists
    % Read the header:
    stuff   = fgets(flist);          % read the header
    disp(stuff);
    rmon    = str2num(stuff(1:2));   % month
    rdat    = str2num(stuff(4:5));   % day
    ryear   = str2num(stuff(7:10));  % year (2012)
    rhr     = str2num(stuff(12:13)); % hour (from 00 to 23)
    rmin    = str2num(stuff(15:16)); % min (00)
    sonicid = str2num(stuff(49:54));
    stux = fscanf(flist,'%7g %6g %5g %5g %5g',[5,36000]);      % Tiksi sonic data format
else
    rmon=NaN;
    rdat=NaN;
    ryear=NaN;
    rhr=NaN;
    rmin=NaN;
    sonicid=NaN;
    rmon=NaN;
    stux=ones(5,36000)+NaN;
end

[a,np] = size(stux);

minr = floor(stux(1,1:np)/100000); 
secr = (stux(1,1:np)-minr*100000)/1000;
tmrw = minr*60+secr;       % put in units of seconds

uraw = stux(2,1:np)/100;   % put in units of m/s
vraw = stux(3,1:np)/100;   % put in units of m/s
wraw = stux(4,1:np)/100;   % put in units of m/s
Traw = stux(5,1:np)/100;   % put in units of deg C
nsqst = 1; nsqen = np;     % the unbroken sequence to do the FFT on

% Replace '-99.99' values by 'NaN':
uraw(uraw==-99.99)=NaN;
vraw(vraw==-99.99)=NaN;
wraw(wraw==-99.99)=NaN;
Traw(Traw==-99.99)=NaN;

% Replace '-Inf' values by 'NaN':
% uraw(uraw==-Inf)=NaN;
% vraw(vraw==-Inf)=NaN;
% wraw(wraw==-Inf)=NaN;
% Traw(Traw==-Inf)=NaN;

%*************************************** Despike procedure. Run function 'despik.m' *****************************************
on=1:1:np;
uz(1:np,1)=uraw';
uz(1:np,2)=on';
u=despik(uz);
uz(1:np,1)=vraw';
v=despik(uz);
uz(1:np,1)=wraw';
w=despik(uz);
uz(1:np,1)=Traw';
T=despik(uz);

% Replace _all_ data by NaN if length of a raw data file is less than 2^13 = 8192 points
% 8192 points = 8192/10 sec  = 13.65333 min (for samp_freq = 10 Hz)
if np < 8192
   u=ones(1,np)*NaN;
   v=ones(1,np)*NaN;
   w=ones(1,np)*NaN;
   T=ones(1,np)*NaN;
end;

% Replace _all_ data by NaN if more than 50% of the raw data of u- or v- or w-wind speed component are 'bad' (-9999=NaN)
if length(u(isnan(u))) > np/2 || length(v(isnan(v))) > np/2 || length(w(isnan(w))) > np/2
   u=ones(1,np)*NaN;
   v=ones(1,np)*NaN;
   w=ones(1,np)*NaN;
   T=ones(1,np)*NaN;
end;

% Replace missing values (NaN) using mean values:
% This method is preferable in the case if there are too many NaNs or long gaps
u(isnan(u)) = mean(u(~isnan(u)));
v(isnan(v)) = mean(v(~isnan(v)));
w(isnan(w)) = mean(w(~isnan(w)));
T(isnan(T)) = mean(T(~isnan(T)));

% Replace missing values (NaN) using interpolation based on the time reference base:
%  (Another possible solution)
% if flist > 0, % if the file exists
%    xi = 1:np; % creates "time" number
%    u = interp1(xi(~isnan(u(xi))),u(~isnan(u(xi))),xi); 
%    v = interp1(xi(~isnan(v(xi))),v(~isnan(v(xi))),xi);
%    w = interp1(xi(~isnan(w(xi))),w(~isnan(w(xi))),xi);
%    T = interp1(xi(~isnan(T(xi))),T(~isnan(T(xi))),xi);
% else
%   stux = ones(5,36000) + NaN;
% end

% Computes means for entire time block (1 hr or less)
um = mean(u);
vm = mean(v);
wm = mean(w);
Tm = mean(T);

%************************* Rotate the anemometer coordinate system to a streamline coordinate system ************************
 
% We used the most common method, which is a double rotation of the anemometer coordinate system, 
% to compute the longitudinal, lateral, and vertical velocity components in real time (Kaimal and Finnigan 1994, Sect. 6.6).
s=sqrt(um^2+vm^2);
thet=atan2(vm,um);
phi=atan2(wm,s);
rot(1,1)=cos(phi)*cos(thet);
rot(1,2)=cos(phi)*sin(thet);
rot(1,3)=sin(phi);
rot(2,1)=-sin(thet);
rot(2,2)=cos(thet);
rot(2,3)=0;
rot(3,1)=-sin(phi)*cos(thet);
rot(3,2)=-sin(phi)*sin(thet);
rot(3,3)=cos(phi);
for j=1:np;
   x=[u(j) v(j) w(j)];
   xr=rot*x';
   ur(j)=xr(1);
   vr(j)=xr(2);
   wr(j)=xr(3);
end;
% !!! However, in situations characterized by a complex or sloping topography or strong canopy dishomogeneity, the planar fit method may be the preferred choice.
% See papers by Wiczak et al. (2001, BLM, 99, 127-150) and van Dijk et al. (2004, Internal Report, Wageningen University) for discussion

% Mean rotated wind speed components:
% vrm and wrm can be used for QC purposes (theoretically vrm = 0 and wrm = 0)
urm = mean(ur);
vrm = mean(vr);
wrm = mean(wr);

% *** Figure 1 ++++++++++++++++++++++++++++++++++++++++++
% Plot full 1-hr time series of the rotated wind components and sonic temperature
t = 0.1:0.1:np/samp_freq;
nt = np;

figure (1); clf
plot(t(1:nt)./60,ur(1:nt), t(1:nt)./60,vr(1:nt), t(1:nt)./60,wr(1:nt), t(1:nt)./60,T(1:nt)-Tm+2, [0 4000],[0 0],'k-');
set(gca, 'linewidth',2, 'fontsize',11, 'XLim',[0 60], 'XTick',[0:10:60])
xlabel('Time, min','fontsize',14); 
ylabel('u,  v,  w,  T - T_m+ 2', 'fontsize',14, 'VerticalAlignment','bottom');
title(['Rotated sonic components: u-blue v-green w-red T-cyan; YD = ' ddd ', ' hh ':00 (UTC)'])  

% Instantaneous relative direction of the wind speed vector (relative to sonic anemometer x-axis):
idir = 1:nt;
wind_dir(:,idir) = atan2(-v(:,idir),u(:,idir))*180/pi;          % Instantaneous relative direction of the wind speed vector (deg)

% Convert wind direction range to [0 360] segment if relative wind direction has only negative values in the 1/2 hr time series:
if (wind_dir(:,idir) < 0)
    wind_dir(:,idir) = wind_dir(:,idir) + 360;
end

% Convert wind direction range to [0 360] segment if relative wind direction has both possitive and negative values in the 1/2 hr time series:
% jj=find(wind_dir(:,ii) < 0);
% wind_dir(jj) = wind_dir(jj) + 360;

% *** Figure 2 ++++++++++++++++++++++++++++++++++++++++++
% Plot full 1-hr time series of the instantaneous relative direction of the wind speed vector (relative to sonic anemometer x-axis)
figure (2); clf
plot(t(1:nt)./60,wind_dir(1:nt), [0 4000],[0 0],'k-');
% axis([0 3600 0 360]);
set(gca, 'linewidth',2, 'fontsize',11, 'XLim',[0 60], 'XTick',[0:10:60])
% set(gca,'yTick',-180:30:360);
xlabel('Time, min','fontsize',14); 
ylabel('Relative wind direction, deg', 'fontsize',14, 'VerticalAlignment','bottom');
title(['Relative wind direction; YD = ' ddd ', ' hh ':00 (UTC)'])  

%******************************************** Compute Sonic Spectra & Cospectra *********************************************

% Define 'nf', the length in gridpoints of the FFT windowing segment
% Here np = length(T) = length(Traw) = length(w) = length(wraw) = ... - length of a raw data file
if np < 8192,
   nf=np;            % condition to avoid different matrix dimensions for short raw data series; 
%   nf=2^12;          % 2^12 = 4096 points  (or = 4096/10 sec  =  6.82667 min for samp_freq = 10 Hz)
elseif np < 16384,
   nf=2^13;          % 2^13 = 8192 points  (or = 8192/10 sec  = 13.65333 min for samp_freq = 10 Hz)
elseif np < 32768,
   nf=2^14;          % 2^14 = 16384 points (or = 16384/10 sec = 27.30667 min for samp_freq = 10 Hz)
else
   nf=2^15;          % 2^15 = 32768 points (or = 32768/10 sec = 54.61333 min for samp_freq = 10 Hz)
end;

nf_min = nf/600;     % nf in minutes

disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>');
disp(['The length of the raw data for a Sonic Anemometer, Number of lines, np = ' num2str(length(T)) ' or ' num2str(length(T)/600) ' min'])
disp(['The length of the FFT segment for a Sonic Anemometer, nf = ' num2str(nf) ' or ' num2str(nf_min) ' min'])
disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>');

% Instantaneous relative direction of the wind speed vector in the rotated coordinate system:
wdirr = atan2(-vr,ur)*180/pi;          % Instantaneous relative direction of the wind speed vector (deg)

us=ur(1:nf);
vs=vr(1:nf);
ws=wr(1:nf);
Ts=T(1:nf);
wdirs=wdirr(1:nf);

% Means for FFT segment (Note that they differ from hourly means, see above):
usm = mean(us);
vsm = mean(vs);
wsm = mean(ws);
Tsm = mean(Ts);
wdirsm = mean(wdirs);

% >>> Perform Fast Fourier Transform (FFT) to compute power spectra and cospectra: 

% In the following version linear detrend is used (no overlaping)
[su,F]=psd(us-usm,nf,10,hamming(nf),'linear');           % (psd = Power Spectral Density)
[sv,F]=psd(vs-vsm,nf,10,hamming(nf),'linear');
[sw,F]=psd(ws-wsm,nf,10,hamming(nf),'linear');
[sT,F]=psd(Ts-Tsm,nf,10,hamming(nf),'linear');
[swu,F]=csd(ws-wsm,us-usm,nf,10,hamming(nf),'linear');   % (csd = Cross Spectral Density)
[swv,F]=csd(ws-wsm,vs-vsm,nf,10,hamming(nf),'linear');
[swT,F]=csd(ws-wsm,Ts-Tsm,nf,10,hamming(nf),'linear');

% In addition to Chris' original code, computation of CuT, CvT, and Cuv are added by AG
% Cospectra CuT & CvT are associated with the horizontal heat flux
[suT,F]=csd(us-usm,Ts-Tsm,nf,10,hamming(nf),'linear');
[svT,F]=csd(vs-vsm,Ts-Tsm,nf,10,hamming(nf),'linear');
[suv,F]=csd(us-usm,vs-vsm,nf,10,hamming(nf),'linear');
% Also spectrum of wind speed direction is added (AG)
[swdir,F]=psd(wdirs-wdirsm,nf,10,hamming(nf),'linear');

% In the following version no overlaping and no linear detrend is used
%[su,F]=psd(us-usm,nf,10,hamming(nf),'none');           % (psd = Power Spectral Density)
%[sv,F]=psd(vs-vsm,nf,10,hamming(nf),'none');
%[sw,F]=psd(ws-wsm,nf,10,hamming(nf),'none');
%[sT,F]=psd(Ts-Tsm,nf,10,hamming(nf),'none');
%[swu,F]=csd(ws-wsm,us-usm,nf,10,hamming(nf),'none');   % (csd = Cross Spectral Density)
%[swv,F]=csd(ws-wsm,vs-vsm,nf,10,hamming(nf),'none');
%[swT,F]=csd(ws-wsm,Ts-Tsm,nf,10,hamming(nf),'none');

% In the SHEBA program the following procedure was used (see 'read_sonics_day_mldsk4.m')
%[su,F]=psd(us,nfw,10,hamming(nfw),novrlp,'linear');
%[sv,F]=psd(vs,nfw,10,hamming(nfw),novrlp,'linear');
%[sw,F]=psd(ws,nfw,10,hamming(nfw),novrlp,'linear');
%[sT,F]=psd(Ts,nfw,10,hamming(nfw),novrlp,'linear');
%[swu,F]=csd(ws,us,nfw,10,hamming(nfw),novrlp,'linear');
%[swv,F]=csd(ws,vs,nfw,10,hamming(nfw),novrlp,'linear');
%[swT,F]=csd(ws,Ts,nfw,10,hamming(nfw),novrlp,'linear');

% >>> Normalize spectra properly (only for psd & csd functions):
su=su*2/10; sv=sv*2/10; sw=sw*2/10; sT=sT*2/10; swu=swu*2/10; swv=swv*2/10; swT=swT*2/10;
suT=suT*2/10; svT=svT*2/10; suv=suv*2/10;
swdir=swdir*2/10;

% ===============================================================================================================
% !!! MatLab Warning: PSD has been replaced by SPECTRUM objects. PSD still works but may be removed in the future. Use SPECTRUM (or its functional form PWELCH) instead.
% !!! MatLab Warning: CSD has been replaced by CPSD. CSD still works but may be removed in future. Use CPSD instead.
% ===============================================================================================================

% >>> Spectra smoothing procedure
n=nf/2;
c=0.1;
jx=1;
dx=1;
ix=1;
inx=1;
while jx<n,
   	dx=(dx*exp(c));
   	d=floor(dx);
   	acu=0;acv=0;acw=0;acT=0;acwu=0;acwv=0;acwT=0;acuT=0;acvT=0;acuv=0;acwdir=0;ac2=0; 
   	k=0;
   	for jx=ix:ix+d-1;
   	   if jx==n,
      	   break;
      	end;
      	acu=acu+su(jx);
      	acv=acv+sv(jx);
      	acw=acw+sw(jx);
      	acT=acT+sT(jx);
      	acwu=acwu+swu(jx);
      	acwv=acwv+swv(jx);
      	acwT=acwT+swT(jx);
        acuT=acuT+suT(jx);
        acvT=acvT+svT(jx);
        acuv=acuv+suv(jx);
        acwdir=acwdir+swdir(jx);
      	ac2=ac2+F(jx);
      	k=k+1;
   	end;
   	sus(inx)=acu/k;
   	svs(inx)=acv/k;
   	sws(inx)=acw/k;
   	sTs(inx)=acT/k;
   	cwus(inx)=acwu/k;
   	cwvs(inx)=acwv/k;
   	cwTs(inx)=acwT/k;
    cuTs(inx)=acuT/k;
    cvTs(inx)=acvT/k;
    cuvs(inx)=acuv/k;
    swdirs(inx)=acwdir/k;
   	fs(inx)=ac2/k;
   	dfs(inx)=F(jx)-F(ix)+F(2);
   	ix=jx+1;
   	inx=inx+1; 
end;

% The total frequency domain 'fs' in spectral or cospectral plots contains:
% 74 spectral points (for 2^15 = 32768 data points = 54.61333 min), size(fs) = [1 74]
% 67 spectral points (for 2^14 = 16384 data points = 27.30667 min), size(fs) = [1 67]
% 60 spectral points (for 2^13 =  8192 data points = 13.65333 min), size(fs) = [1 60]
% 53 spectral points (for 2^12 =  4096 data points =  6.82667 min), size(fs) = [1 53]

%********************************************** Plot Sonic Spectra & Cospectra **********************************************

% *** Figure 3 ++++++++++++++++++++++++++++++++++++++++++
% plot smoothed variance spectra
figure (3); clf
loglog(fs,sus.*fs, fs,svs.*fs, fs,sws.*fs, fs,sTs.*fs);
set(gca, 'linewidth',2, 'fontsize',11, 'XLim',[0.0001 10], 'XTick',[0.0001 0.001 0.01 0.1 1 10])
xlabel('Frequency,  f (Hz)','fontsize',14);
ylabel('f*S_{x}(f)', 'fontsize',14, 'VerticalAlignment','bottom');
title(['Variance spectra: u-blue v-green w-red T-cyan; YD = ' ddd ', hr = ' hh ':00 (UTC)'])

% *** Figure 4 ++++++++++++++++++++++++++++++++++++++++++
% plot smoothed cospectra for Cwu, Cwv, CwT
figure (4); clf
semilogx(fs,cwus.*fs, fs,cwvs.*fs, fs,cwTs.*fs, [1e-4 10],[0 0],'k-');
set(gca, 'linewidth',2, 'fontsize',11, 'XLim',[0.0001 10], 'XTick',[0.0001 0.001 0.01 0.1 1 10])
xlabel('Frequency,  f (Hz)','fontsize',14);
ylabel('f*C_{wx}(f)', 'fontsize',14, 'VerticalAlignment','bottom');
title(['Cospectra: <wu>-blue <wv>-green <wT>-red; YD = ' ddd ', hr = ' hh ':00 (UTC)'])

% *** Figure 5 ++++++++++++++++++++++++++++++++++++++++++
% plot smoothed cospectra for CuT, CvT, Cuv 
figure (5); clf
semilogx(fs,cuTs.*fs, fs,cvTs.*fs, fs,cuvs.*fs, [1e-4 10],[0 0],'k-');
set(gca, 'linewidth',2, 'fontsize',11, 'XLim',[0.0001 10], 'XTick',[0.0001 0.001 0.01 0.1 1 10])
xlabel('Frequency,  f (Hz)','fontsize',14);
ylabel('f*C_{xy}(f)', 'fontsize',14, 'VerticalAlignment','bottom');
title(['Cospectra: <uT>-blue <vT>-green <uv>-red; YD = ' ddd ', hr = ' hh ':00 (UTC)']) 

% *** Figure 6 ++++++++++++++++++++++++++++++++++++++++++
% plot smoothed variance spectrum of wind speed direction
% figure (6); clf
% loglog(fs,swdirs.*fs);
% set(gca, 'linewidth',2, 'fontsize',11, 'XLim',[0.0001 10], 'XTick',[0.0001 0.001 0.01 0.1 1 10])
% xlabel('Frequency,  f (Hz)','fontsize',14);
% ylabel('f*S_{wdir}(f)', 'fontsize',14, 'VerticalAlignment','bottom');
% title(['Variance spectrum of wind direction; YD = ' ddd ', hr = ' hh ':00 (UTC)']) 

% *** Figure 6 ++++++++++++++++++++++++++++++++++++++++++
% plot smoothed cospectra for CwT, CuT, CvT 
figure (6); clf
semilogx(fs,cwTs.*fs,'-', fs,cuTs.*fs,'-', fs,cvTs.*fs,'-.', [1e-4 10],[0 0],'k-');
set(gca, 'linewidth',2, 'fontsize',11, 'XLim',[0.0001 10], 'XTick',[0.0001 0.001 0.01 0.1 1 10])
xlabel('Frequency,  f (Hz)','fontsize',14);
ylabel('f*C_{xy}(f)', 'fontsize',14, 'VerticalAlignment','bottom');
title(['Cospectra: <wT>-blue <uT>-green <vT>-red; YD = ' ddd ', hr = ' hh ':00 (UTC)']) 

%++++++++++++++++++++++++ Important constants and equations ++++++++++++++++++++++++++++++++
% !!! Tsm - temperature in deg C
tdk = 273.15;                  % Absolute Zero
Rgas = 287.1;                  % Universal gas constant, J/(kg K)
sigma = 5.67e-8;               % Stefan-Boltzman constant, W/(m^2 K^4)
p = 1013;                      % air pressure (mb) when no pressure measurements
qo = 0;                        % air humidity when no humidity measurements   
rho = p*100/(Rgas*(Tsm+tdk)*(1+0.61e-3*qo)); % density of air derived from ideal gas equation, kg/m^3
cp = 1005.6+0.017211*Tsm+0.000392*Tsm*Tsm; % isobaric specific heat of air (median), J/(kg K), (from Ed Andreas)
visa = 1.326e-5*(1+6.542e-3*Tsm+8.301e-6*Tsm*Tsm-4.84e-9*Tsm*Tsm*Tsm); % m^2/s, molecular viscosity from the TOGA CORE bulk algorithm
kt = (0.02411*(1+0.003309*Tsm-1.441e-6*Tsm*Tsm))/(rho*cp); % m^2/s, koefficient of molecular thermal diffusivity (from Ed Andreas)
Le = (2.501-.00237*Tsm)*1e6;   % Latent heat of vaporization, J/kg

%+++++++++++++++++++++++++++++ Wind speed and direction ++++++++++++++++++++++++++++++++++++

wsp = sqrt(um^2 + vm^2 + wm^2);       % 1/2 hour averaged wind speed (m/s)
% wsp = sqrt(usm^2 + vsm^2 + wsm^2);    % 27.3067 min averaged wind speed (m/s)

wdir = atan2(-vm,um)*180/pi;          % 1/2 hour averaged relative wind direction (deg)
% Convert wind direction range to [0 360] segment:
if (wdir < 0)
    wdir = wdir + 360;
end

% All directions for wind and stress here are calculated using the meteorological convention ('from').
% All directions are relatively true (Geographic) North, not magnetic North!

% For true wind direction add a true orientation of a sonic anemometer:
truewdir = wdir + sonic_9m_dir;                 % 1-hr averaged true wind direction (deg)
% Here orientation of the sonic anemoneter (sonic_9m_dir = 180 + 17 deg) is provided by Irina Repina
% Convert wind direction range to [0 360] segment:
if (truewdir > 360)
    truewdir = truewdir - 360;
end

disp('***************************************************************');
disp(['Sonic sonic anemometer at measurement height = ' num2str(z_level_n) ' m'])
disp(['Wind speed = ' num2str(wsp) ' m/s']);
disp(['Relative wind direction = ' num2str(wdir) ' deg']);
disp(['True wind direction = ' num2str(truewdir) ' deg']);
disp(['Air (sonic) temperature = ' num2str(Tsm) ' deg C']);

%++++++++++++++++++++++++++++ Fluxes Based on the Cospectra +++++++++++++++++++++++++++++++

% >>> Compute covariances ("fluxes") based on the appropriate cospectra integration (the total flux is the sum of cospectra frequency components)
% The turbulent fluxes and variances are derived through frequency integration of the appropriate cospectra and spectra
cwux = sum(cwus.*dfs);
cwvx = sum(cwvs.*dfs);
cuvx = sum(cuvs.*dfs);
cwTx = sum(cwTs.*dfs);
cuTx = sum(cuTs.*dfs);
cvTx = sum(cvTs.*dfs);
% >>> Compute standard deviations based on the appropriate spectra integration:
sigu_spc = sqrt(sum(sus.*dfs));       % based on the Su-spectra integration
sigv_spc = sqrt(sum(svs.*dfs));       % based on the Sv-spectra integration
sigw_spc = sqrt(sum(sws.*dfs));       % based on the Sw-spectra integration
sigT_spc = sqrt(sum(sTs.*dfs));       % based on the St-spectra integration
sigwdir_spc = sqrt(sum(swdirs.*dfs)); % based on the Sdir-spectra integration
ssm = sqrt(usm*usm+vsm*vsm);

% !! Note that sigwdir_spc may be used for quality control (or QC for short) of the data
% For example sigwdir_spc > 15 deg may be considered as "bad" data (non-stationary)
% However this QC is not applicable for free convection limit (light winds)

% >>> Wind stress components and the friction velocity:
wu_csp = real(cwux);   % wu-covariance based on the wu-cospectra integration
wv_csp = real(cwvx);   % wv-covariance based on the wv-cospectra integration
uv_csp = real(cuvx);   % uv-covariance based on the uv-cospectra integration
ustar = - sign(wu_csp)*sqrt(abs(wu_csp)); % the friction velocity based only on the downstream, uw,  stress components (m/s)
% ustar = - sign(wu_csp)*((abs(wu_csp))^2 + (abs(wv_csp))^2)^(1/4);  % ustar is based on both stress components (m/s)

% >>> Sensible heat flux:
% Note that a sonic anemometer/thermometer measures the so-called 'sonic' temperature which is close to the virtual temperature
% Thus <w't'> here is a buoyancy flux
% Moisture correction is necessary to obtain the correct sensible heat flux (e.g., Kaimal and Finnigan, 1994)
% However, in Arctic and Antarctic this value very close to the temperature flux (i.e. the sensible heat flux)
wT_csp = real(cwTx);   % wt-covariance, vertical flux of the sonic temperature  [deg m/s]
uT_csp = real(cuTx);   % ut-covariance, horizontal flux of the sonic temperature (along-wind component) [deg m/s]
vT_csp = real(cvTx);   % vt-covariance, horizontal flux of the sonic temperature (cross-wind component) [deg m/s]
Hs = wT_csp*rho*cp;    % sensible heat flux (W/m2) - Based on the sonic temperature!!!
Tstar = - wT_csp./abs(ustar); % the temperature scale 

% The loop below is added for cases when number of points < 32768 (54.61333 min) or < 16384 (27.30667 min)
% It defines frequencies for fluxes in the high-frequency subrange and defines the inertial subrange
if np < 8192,
   fsl = 4;  fsn = 53;         % 2^12 = 4096 points  (or = 4096/10 sec  =  6.82667 min for samp_freq = 10 Hz)
   fsi01=34; fsi02=35; fsi03=36; fsi04=37; fsi05=38; fsi06=39; fsi07=40; fsi08=41; fsi09=42; fsi10=43; fsi11=44; fsi12=45; 
elseif np < 16384,
   fsl = 6;  fsn = 60;         % 2^13 = 8192 points  (or =  8192/10 sec = 13.65333 min for samp_freq = 10 Hz)
   fsi01=40; fsi02=41; fsi03=42; fsi04=43; fsi05=44; fsi06=45; fsi07=46; fsi08=47; fsi09=48; fsi10=49; fsi11=50; fsi12=51; 
elseif np < 32768,
   fsl = 9;  fsn = 67;         % 2^14 = 16384 points (or = 16384/10 sec = 27.30667 min for samp_freq = 10 Hz)
   fsi01=47; fsi02=48; fsi03=49; fsi04=50; fsi05=51; fsi06=52; fsi07=53; fsi08=54; fsi09=55; fsi10=56; fsi11=57; fsi12=58; 
else
   fsl = 13; fsn = 74;         % 2^15 = 32768 points (or = 32768/10 sec = 54.61333 min for samp_freq = 10 Hz)
   fsi01=54; fsi02=55; fsi03=56; fsi04=57; fsi05=58; fsi06=59; fsi07=60; fsi08=61; fsi09=62; fsi10=63; fsi11=64; fsi12=65; 
end;

% >>> Compute covariances ("fluxes") in the high-frequency part of the appropriate cospectra
% To separate the small-scale turbulence from mesoscale motions (mostly in SBL), a low-frequency cut-off was applied on the spectra/cospectra as a lower limit of integration.
% Low-frequency cut-off fsn = fs(13) = 0.0064 Hz (the 13-th spectral value) for 2^15 = 32768 data points = 54.61333 min, 
% Low-frequency cut-off fsn = fs(9)  = 0.0064 Hz (the  9-th spectral value) for 2^14 = 16384 data points = 27.30667 min, 
% Low-frequency cut-off fsn = fs(6)  = 0.0061 Hz (the  6-th spectral value) for 2^13 =  8192 data points = 13.65333 min
% These cut-off values are consistent with the low-frequency cut-off at 0.0061 Hz for SHEBA data (Grachev et al. 2013, BLM, 147(1), 51-82, doi 10.1007/s10546-012-9771-0) 
% The upper limit of integration is 5 Hz (the Nyquist frequency)
wu_csp_hi = sum(real(cwus(fsl:fsn)).*dfs(fsl:fsn));
wv_csp_hi = sum(real(cwvs(fsl:fsn)).*dfs(fsl:fsn));
uv_csp_hi = sum(real(cuvs(fsl:fsn)).*dfs(fsl:fsn));
wT_csp_hi = sum(real(cwTs(fsl:fsn)).*dfs(fsl:fsn));
uT_csp_hi = sum(real(cuTs(fsl:fsn)).*dfs(fsl:fsn));
vT_csp_hi = sum(real(cvTs(fsl:fsn)).*dfs(fsl:fsn));
ustar_hi = - sign(wu_csp_hi)*sqrt(abs(wu_csp_hi)); % friction velocity based on the high-frequency part of downstream, uw_hi,  stress components (m/s)
Tstar_hi = - wT_csp_hi./abs(ustar_hi); % the temperature scale based on the high-frequency part of the turbulent fluxes
Hs_hi = wT_csp_hi*rho*cp;    % sensible heat flux (sonic temperature!) based on the high-frequency part of the cospectrum(W/m2)
% >>> Compute standard deviations in the high-frequency part of the appropriate spectra
sigu_spc_hi = sqrt(sum(sus(fsl:fsn).*dfs(fsl:fsn)));
sigv_spc_hi = sqrt(sum(svs(fsl:fsn).*dfs(fsl:fsn)));
sigw_spc_hi = sqrt(sum(sws(fsl:fsn).*dfs(fsl:fsn)));
sigT_spc_hi = sqrt(sum(sTs(fsl:fsn).*dfs(fsl:fsn)));
sigwdir_spc_hi = sqrt(sum(swdirs(fsl:fsn).*dfs(fsl:fsn))); 

% >>> Some Monin-Obukhov (MO) parameters based on the local turbulent measurements:
% Monin-Obukhov stability parameter, z/L:
zeta_level_n = - ((0.4*9.81)./(Tsm+tdk)).*(z_level_n.*wT_csp./(ustar.^3));
% Drag coefficient, Cd:
Cd = - wu_csp./(wsp.^2);
% MO universal functions for the standard deviations:
phi_u = sigu_spc./ustar; phi_v = sigv_spc./ustar; phi_w = sigw_spc./ustar; phi_T = sigT_spc./abs(Tstar);
% MO universal function for the horizontal heat flux:
phi_uT = uT_csp./wT_csp;

% >>> Some Monin-Obukhov (MO) parameters based on the local turbulent measurements (same as above but based on the high-frequency part of the turbulent fluxes):
% Monin-Obukhov stability parameter based on the high-frequency part of the turbulent fluxes, z/L_hi:
zeta_level_n_hi = - ((0.4*9.81)./(Tsm+tdk)).*(z_level_n.*wT_csp_hi./(ustar_hi.^3));
% Drag coefficient based on the high-frequency part of the momentum flux, Cd_hi:
Cd_hi = - wu_csp_hi./(wsp.^2);
% MO universal functions for the standard deviations based on the high-frequency part of the turbulent fluxes:
phi_u_hi = sigu_spc_hi./ustar_hi; phi_v_hi = sigv_spc_hi./ustar_hi; phi_w_hi = sigw_spc_hi./ustar_hi; phi_T_hi = sigT_spc_hi./abs(Tstar_hi);
% MO universal function for the horizontal heat flux based on the high-frequency part of the turbulent fluxes:
phi_uT_hi = uT_csp_hi./wT_csp_hi;

disp('Turbulent fluxes based on the cospectra:');
disp(['Friction velocity = ' num2str(ustar) ' m/s']);
disp(['Sensible heat flux = ' num2str(Hs) ' W/m2 (based on the sonic temperature!)']);
disp(['Drag coefficient, Cd = ' num2str(Cd) ' ']);
disp(['Monin-Obukhov (MO) stability parameter, z/L = ' num2str(zeta_level_n) ' ']);
disp(['MO universal functions for standard deviations (hi-freq), phi_u = ' num2str(phi_u_hi) '; phi_v = ' num2str(phi_v_hi) '; phi_w = ' num2str(phi_w_hi) '; phi_T = ' num2str(phi_T_hi) '']);
disp(['MO universal functions for horizontal heat flux (hi-freq), phi_uT = <uT>/<wT> = ' num2str(phi_uT_hi) '']);

%+++++++++++++++++++ Compute structure parameters in the inertial subrange +++++++++++++++++
%+++++++++++++++++++++++++++++ –5/3 Kolmogorov power law +++++++++++++++++++++++++++++++++++

% Structure parameters are used to derive dissipation rate of the turbulent kinetic energy ('epsilon') and half the temperature variance (Nt)
% The inertial subrange is considered in the frequency subrange ~ 0.6-2 Hz
% Generally this subrange depends on a height of measurements and may be different for w-component
% For 2^15 = 32768 data points = 54.61333 min, the structure parameters are computed in the frequency domain between the 54th and 65th spectral values:
%    fs(fsi01:fsi12) = fs(54:65) = [0.6662 0.7372 0.8156 0.9023 0.9981 1.1041 1.2213 1.3507 1.4937 1.6518 1.8265 2.0197]
% For 2^14 = 16384 data points = 27.30667 min, the structure parameters are computed in the frequency domain between the 47th and 58th spectral values: 
%    fs(fsi01:fsi12) = fs(47:58) = [0.6531 0.7233 0.8011 0.8871 0.9824 1.0876 1.2039 1.3324 1.4743 1.6312 1.8045 1.9962]
% For 2^13 =  8192 data points = 13.65333 min, the structure parameters are computed in the frequency domain between the 40th and 51st spectral values:
%    fs(fsi01:fsi12) = fs(40:51) = [0.6342 0.7037 0.7806 0.8655 0.9595 1.0638 1.1792 1.3062 1.4465 1.6022 1.7743 1.9647]
% For 2^12 =  4096 data points =  6.82667 min, the structure parameters are computed in the frequency domain between the 34th and 45th spectral values:
%    fs(fsi01:fsi12) = fs(34:45) = [0.6738 0.7495 0.8337 0.9265 1.0291 1.1426 1.2683 1.4075 1.5613 1.7310 1.9189 2.1277]

% Structure parameters definition: 
% Cu^2 = 4*alpha*epsilon^{2/3} (alpha=0.55); Ct^2 = 4*beta*Nt*epsilon^{-1/3} (beta=0.8)
% Relationships between the structure parameters and the frequency spectra:
% Cu^2 = 4*(2pi/U)^{2/3}*Su(f)*f^{-5/3}; Ct^2 = 4*(2pi/U)^{2/3}*St(f)*f^{-5/3};
% Dimensions: [Cu^2] = [Ustar^2/z^{2/3}]; [Ct^2] = [Tstar^2/z^{2/3}];
% Local isotropy in the inertial subrange assumes: 
% Fv(k)=Fw(k)=(4/3)Fu(k), or Sv(f)=Sw(f)=(4/3)Su(f), or Cv^2=Cw^2=(4/3)Cu^2.
% Conversion from wavenumber to frequency scales: k*Fu(k)=f*Su(f).
% See details in the textbook by J.C. Kaimal & J.J. Finnigan "Atmospheric Boundary Layer Flows" (1994)

gfac = 4*(2*pi/wsp)^.667;
cu2 = gfac*median(sus(fsi01:fsi12).*fs(fsi01:fsi12).^1.667);              % Cu^2, U structure parameter [x^2/m^2/3 = m^1.33/s^2]
cv2 = gfac*median(svs(fsi01:fsi12).*fs(fsi01:fsi12).^1.667);              % Cv^2, V structure parameter [x^2/m^2/3 = m^1.33/s^2]
cw2 = gfac*median(sws(fsi01:fsi12).*fs(fsi01:fsi12).^1.667);              % Cw^2, W structure parameter [x^2/m^2/3 = m^1.33/s^2]
cT2 = gfac*median((sTs(fsi01:fsi12)-min(sTs)).*fs(fsi01:fsi12).^1.667);   % Ct^2, T structure parameter [x^2/m^2/3 = deg^2/m^2/3]; (term -min(sTs) reduces some noise - proposed by Chris F.)
% cT2 = gfac*median((sTs(fsi01:fsi12)).*fs(fsi01:fsi12).^1.667);            % Ct^2, T structure parameter [x^2/m^2/3 = deg^2/m^2/3]

%+++++++++++++++++ Compute dissipation rate of the turbulent kinetic energy ++++++++++++++++

% Non-dimensional dissipation rate of the turbulent kinetic energy (TKE):
% Phi_epsilon = kappa*z*epsilon/Ustar^3 - see Eq. (1.30) in Kaimal & Finnigan (1994)
% Substitution of epsilon = (Cu^2/4*alpha)^(3/2) leads to
% Phi_epsilon = kappa*[(z^2/3)*Cu^2/4*alpha*Ustar^2] = kappa*(Cu^2n/4*alpha)^3/2, where

alphaK = 0.55; % - the Kolmogorov constant (alpha=0.55 - see above)
% Note that Cu^2=(3/4)Cv^2=(3/4)Cw^2

% >>> For epsilon_u (Based on Su) ++++++++++++++++++++++++++++++++++++++ 
% Dissipation rate of the turbulent kinetic energy based on the energy spectra of the longitudinal velocity component in the inertial subrange [m^2/s^3]:
epsilon_u = (cu2/(4*alphaK)).^(3/2); 
% >>> For epsilon_v (Based on Sv) ++++++++++++++++++++++++++++++++++++++ 
% Dissipation rate of the turbulent kinetic energy based on the energy spectra of the lateral velocity component in the inertial subrange [m^2/s^3]:
epsilon_v = (3/4)*(cv2/(4*alphaK)).^(3/2); 
% >>> For epsilon_w (Based on Sw) ++++++++++++++++++++++++++++++++++++++ 
% Dissipation rate of the turbulent kinetic energy based on the energy spectra of the vertical velocity component in the inertial subrange [m^2/s^3]:
epsilon_w = (3/4)*(cw2/(4*alphaK)).^(3/2); 
% >>> Median epsilon (median of epsilon_u, epsilon_v, & epsilon_w) ++++++++++++++++++++++++++++++++++++++ 
% Dissipation rate of the turbulent kinetic energy = median of the values derived from Su, Sv, & Sw [m^2/s^3]:
epsilon = median([epsilon_u epsilon_v epsilon_w],2);
% Monin-Obukhov universal function Phi_epsilon based on the median epsilon:
Phi_epsilon = (0.4*z_level_n.*epsilon)./(ustar.^3);
% Monin-Obukhov universal function Phi_epsilon based on the median epsilon based on the high-frequency part of the momentum flux:
Phi_epsilon_hi = (0.4*z_level_n.*epsilon)./(ustar_hi.^3);

disp(['MO universal functions for dissipation rate of the (TKE) (hi-freq), phi_epsilon = 0.4*z*epsilon/Ustar^3 = ' num2str(Phi_epsilon_hi) '']);

%+++++ Compute the dissipation (destruction) rate for half the temperature variance +++++++
%++++++++++++++++++ –5/3 Obukhov-Corrsin power law for the passive scalar ++++++++++++++++++

%  Temperature Structure Parameter: Ct^2 = 4*beta*Nt*epsilon^{-1/3} (beta=0.8)
betaK = 0.8; % - the Kolmogorov (Obukhov-Corrsin) constant for the passive scalar (beta=0.8 - see above)

% >>> Nt (Based on Ct^2 derived from St) ++++++++++++++++++++++++++++++++++++++ 
% The dissipation (destruction) rate for half the temperature variance [deg^2/s]:
Nt = (cT2.*(epsilon.^(1/3)))/(4*betaK); 
% Monin-Obukhov universal function Phi_Nt:
Phi_Nt = (0.4*z_level_n.*Nt)./(ustar.*Tstar^2);
% Monin-Obukhov universal function Phi_Nt based on the high-frequency part of the turbulent fluxes:
Phi_Nt_hi = (0.4*z_level_n.*Nt)./(ustar_hi.*Tstar_hi^2);

disp(['MO universal functions for destruction rate for half the temperature variance (hi-freq), phi_Nt = 0.4*z*Nt/(Ustar*Tstar^2) = ' num2str(Phi_Nt_hi) '']);
disp('***************************************************************');

%++++++++++++++++++++++ Compute spectral slopes in the inertial subrange +++++++++++++++++++

% >>> Spectrla slopes are computed in the frequency domain defined above
disp('Test of the –5/3 Kolmogorov and Obukhov-Corrsin power laws (in the frequency domain 0.67–2.02 Hz):');

% Compute a spectral slope in the inertial subrange for Su. 
% Individual slopes:
nSu_1 = log(sus(:,fsi07)./sus(:,fsi01))./log(fs(fsi07)./fs(fsi01)); 
nSu_2 = log(sus(:,fsi08)./sus(:,fsi02))./log(fs(fsi08)./fs(fsi02)); 
nSu_3 = log(sus(:,fsi09)./sus(:,fsi03))./log(fs(fsi09)./fs(fsi03)); 
nSu_4 = log(sus(:,fsi10)./sus(:,fsi04))./log(fs(fsi10)./fs(fsi04)); 
nSu_5 = log(sus(:,fsi11)./sus(:,fsi05))./log(fs(fsi11)./fs(fsi05)); 
nSu_6 = log(sus(:,fsi12)./sus(:,fsi06))./log(fs(fsi12)./fs(fsi06)); 
% Median spectral slope in the inertial subrange:
nSu = median([nSu_1 nSu_2 nSu_3 nSu_4 nSu_5 nSu_6],2);
disp(['Median spectral slope in the inertial subrange for Su : = ' num2str(nSu)]) 

% Compute a spectral slope in the inertial subrange for Sv. 
% Individual slopes:
nSv_1 = log(svs(:,fsi07)./svs(:,fsi01))./log(fs(fsi07)./fs(fsi01)); 
nSv_2 = log(svs(:,fsi08)./svs(:,fsi02))./log(fs(fsi08)./fs(fsi02)); 
nSv_3 = log(svs(:,fsi09)./svs(:,fsi03))./log(fs(fsi09)./fs(fsi03)); 
nSv_4 = log(svs(:,fsi10)./svs(:,fsi04))./log(fs(fsi10)./fs(fsi04)); 
nSv_5 = log(svs(:,fsi11)./svs(:,fsi05))./log(fs(fsi11)./fs(fsi05)); 
nSv_6 = log(svs(:,fsi12)./svs(:,fsi06))./log(fs(fsi12)./fs(fsi06)); 
% Median spectral slope in the inertial subrange:
nSv = median([nSv_1 nSv_2 nSv_3 nSv_4 nSv_5 nSv_6],2);
disp(['Median spectral slope in the inertial subrange for Sv : = ' num2str(nSv)]) 

% Compute a spectral slope in the inertial subrange for Sw. 
% Individual slopes:
nSw_1 = log(sws(:,fsi07)./sws(:,fsi01))./log(fs(fsi07)./fs(fsi01)); 
nSw_2 = log(sws(:,fsi08)./sws(:,fsi02))./log(fs(fsi08)./fs(fsi02)); 
nSw_3 = log(sws(:,fsi09)./sws(:,fsi03))./log(fs(fsi09)./fs(fsi03)); 
nSw_4 = log(sws(:,fsi10)./sws(:,fsi04))./log(fs(fsi10)./fs(fsi04)); 
nSw_5 = log(sws(:,fsi11)./sws(:,fsi05))./log(fs(fsi11)./fs(fsi05)); 
nSw_6 = log(sws(:,fsi12)./sws(:,fsi06))./log(fs(fsi12)./fs(fsi06)); 
% Median spectral slope in the inertial subrange:
nSw = median([nSw_1 nSw_2 nSw_3 nSw_4 nSw_5 nSw_6],2);
disp(['Median spectral slope in the inertial subrange for Sw : = ' num2str(nSw)]) 

% Compute a spectral slope in the inertial subrange for St. 
% Individual slopes:
nSt_1 = log(sTs(:,fsi07)./sTs(:,fsi01))./log(fs(fsi07)./fs(fsi01)); 
nSt_2 = log(sTs(:,fsi08)./sTs(:,fsi02))./log(fs(fsi08)./fs(fsi02)); 
nSt_3 = log(sTs(:,fsi09)./sTs(:,fsi03))./log(fs(fsi09)./fs(fsi03)); 
nSt_4 = log(sTs(:,fsi10)./sTs(:,fsi04))./log(fs(fsi10)./fs(fsi04)); 
nSt_5 = log(sTs(:,fsi11)./sTs(:,fsi05))./log(fs(fsi11)./fs(fsi05)); 
nSt_6 = log(sTs(:,fsi12)./sTs(:,fsi06))./log(fs(fsi12)./fs(fsi06)); 
% Median spectral slope in the inertial subrange:
nSt = median([nSt_1 nSt_2 nSt_3 nSt_4 nSt_5 nSt_6],2);
disp(['Median spectral slope in the inertial subrange for St : = ' num2str(nSt)]) 

% Compute a spectral slope in the inertial subrange for Swdirs (wind direction spectrum). 
% Individual slopes:
nSwdir_1 = log(swdirs(:,fsi07)./swdirs(:,fsi01))./log(fs(fsi07)./fs(fsi01)); 
nSwdir_2 = log(swdirs(:,fsi08)./swdirs(:,fsi02))./log(fs(fsi08)./fs(fsi02)); 
nSwdir_3 = log(swdirs(:,fsi09)./swdirs(:,fsi03))./log(fs(fsi09)./fs(fsi03)); 
nSwdir_4 = log(swdirs(:,fsi10)./swdirs(:,fsi04))./log(fs(fsi10)./fs(fsi04)); 
nSwdir_5 = log(swdirs(:,fsi11)./swdirs(:,fsi05))./log(fs(fsi11)./fs(fsi05)); 
nSwdir_6 = log(swdirs(:,fsi12)./swdirs(:,fsi06))./log(fs(fsi12)./fs(fsi06)); 
% Median spectral slope in the inertial subrange:
nSwdir = median([nSwdir_1 nSwdir_2 nSwdir_3 nSwdir_4 nSwdir_5 nSwdir_6],2);
disp(['Median spectral slope in the inertial subrange for Swdir : = ' num2str(nSwdir)]) 

% Note that the spectral slopes can be used as QC parameters, e.g. as indicator of resolution limit of a sonic anemometer
% This resolution leads to a step ladder appearance in the data time series, and the turbulent fluxes cannot be reliably calculated
% See Fig. 1b in Vickers & Mahrt (1997)J. Atmos. Oc. Tech. 14(3): 512–526
% Eventually, in the very stable regime, the spectral slope levels off asymptotically at zero which corresponds to 'white noise' in the sensor 
% See spectra at level 5 in Fig. 3 in Grachev et al. (2008) Acta Geophysica. 56(1): 142–166.
% For example, conditions nSu, nSv, nSu, and nSt > - 0.5 (or even > - 1) can be used as QC thresholds

%+++++++++++++++Second-order moments of atmospheric turbulence (including fluxes)+++++++++++
%+++++++++++++++ (Direct variances and covariances based on Reynolds averaging) ++++++++++++

% Removing linear trends from the data (detrending the data):
us_dtr = detrend(us);
vs_dtr = detrend(vs);
ws_dtr = detrend(ws);
Ts_dtr = detrend(Ts);

% Trend lines:
us_trend = us - us_dtr; 
vs_trend = vs - vs_dtr; 
ws_trend = ws - ws_dtr; 
Ts_trend = Ts - Ts_dtr; 

% Mean of the detrended data (should be very close to zero)
us_dtr_m = mean(us_dtr);
vs_dtr_m = mean(vs_dtr);
ws_dtr_m = mean(ws_dtr);
Ts_dtr_m = mean(Ts_dtr);

% >>> Covariances and Variances of the detrending data:
cov_wu = cov(ws_dtr,us_dtr);
cov_wv = cov(ws_dtr,vs_dtr);
cov_uv = cov(us_dtr,vs_dtr);
cov_wT = cov(ws_dtr,Ts_dtr);
cov_uT = cov(us_dtr,Ts_dtr);
cov_vT = cov(vs_dtr,Ts_dtr);
wu_cov = cov_wu(1,2); 
wv_cov = cov_wv(1,2); 
uv_cov = cov_uv(1,2); 
wT_cov = cov_wT(1,2); 
uT_cov = cov_uT(1,2); 
vT_cov = cov_vT(1,2); 
sigu_cov = sqrt(cov_wu(2,2)); % Standard deviation of the u-wind component
sigv_cov = sqrt(cov_wv(2,2)); % Standard deviation of the v-wind component
sigw_cov = sqrt(cov_wu(1,1)); % Standard deviation of the w-wind component
sigT_cov = sqrt(cov_wT(2,2)); % Standard deviation of the the air temperature (sonic temperature)
sigu_std = std(us_dtr); % Standard deviation of the u-wind component
sigv_std = std(vs_dtr); % Standard deviation of the v-wind component
sigw_std = std(ws_dtr); % Standard deviation of the w-wind component
sigT_std = std(Ts_dtr); % Standard deviation of the the air temperature (sonic temperature)

% >>> Second-order moments (directly or by definition)
% Linear detrending but no block-averaging of the data!
% >>> Note that this is the same as 'cov' MatLab command above 
% nf = 32768 (54.6133 min averaged data)
uu_dtr = (sum(us_dtr.*us_dtr))./nf;
uv_dtr = (sum(us_dtr.*vs_dtr))./nf;
uw_dtr = (sum(us_dtr.*ws_dtr))./nf;
vv_dtr = (sum(vs_dtr.*vs_dtr))./nf;
vw_dtr = (sum(vs_dtr.*ws_dtr))./nf;
ww_dtr = (sum(ws_dtr.*ws_dtr))./nf;
uT_dtr = (sum(us_dtr.*Ts_dtr))./nf;
vT_dtr = (sum(vs_dtr.*Ts_dtr))./nf;
wT_dtr = (sum(ws_dtr.*Ts_dtr))./nf;
TT_dtr = (sum(Ts_dtr.*Ts_dtr))./nf;
sigu_dtr = sqrt(uu_dtr); % Standard deviation of the u-wind component
sigv_dtr = sqrt(vv_dtr); % Standard deviation of the v-wind component
sigw_dtr = sqrt(ww_dtr); % Standard deviation of the w-wind component
sigT_dtr = sqrt(TT_dtr); % Standard deviation of the the air temperature (sonic temperature)

% Note that sigu_cov = sigu_std = sigu_dtr (the same for 'v', 'w', and 'T') - different methods of computation

% >>> Standard deviations of the detrending second-order moments (fluxes and variances):
sig_uv_dtr = std(us_dtr.*vs_dtr);
sig_uw_dtr = std(us_dtr.*ws_dtr);
sig_vw_dtr = std(vs_dtr.*ws_dtr);
sig_uT_dtr = std(us_dtr.*Ts_dtr);
sig_vT_dtr = std(vs_dtr.*Ts_dtr);
sig_wT_dtr = std(ws_dtr.*Ts_dtr);
sig_uu_dtr = std(us_dtr.*us_dtr);
sig_vv_dtr = std(vs_dtr.*vs_dtr);
sig_ww_dtr = std(ws_dtr.*ws_dtr);
sig_TT_dtr = std(Ts_dtr.*Ts_dtr);

% >>> Second-order moments (directly or by definition)
% Same as above, but no detrending and no block-averaging of the data!
% nf = 32768 (54.6133 min averaged data)
uu_direct = (sum((us-usm).*(us-usm)))./nf;
uv_direct = (sum((us-usm).*(vs-vsm)))./nf;
uw_direct = (sum((us-usm).*(ws-wsm)))./nf;
vv_direct = (sum((vs-vsm).*(vs-vsm)))./nf;
vw_direct = (sum((vs-vsm).*(ws-wsm)))./nf;
ww_direct = (sum((ws-wsm).*(ws-wsm)))./nf;
uT_direct = (sum((us-usm).*(Ts-Tsm)))./nf;
vT_direct = (sum((vs-vsm).*(Ts-Tsm)))./nf;
wT_direct = (sum((ws-wsm).*(Ts-Tsm)))./nf;
TT_direct = (sum((Ts-Tsm).*(Ts-Tsm)))./nf;
sigu_direct = sqrt(uu_direct);
sigv_direct = sqrt(vv_direct);
sigw_direct = sqrt(ww_direct);
sigT_direct = sqrt(TT_direct);

% Ratio of the fluxes (standard deviations) derived from covariances to the appropriate values derived from cospectra (spectra): 
wu_ratio = wu_cov/wu_csp;
wv_ratio = wv_cov/wv_csp;
wT_ratio = wT_cov/wT_csp;
sigu_ratio = sigu_cov/sigu_spc;
sigv_ratio = sigv_cov/sigv_spc;
sigw_ratio = sigw_cov/sigw_spc;
sigT_ratio = sigT_cov/sigT_spc;

% See details of the linear detrending in:
% Gash, J. H. C. and A. D. Culf. 1996. Applying linear de-trend to eddy correlation data in real time. Boundary-Layer Meteorology, 79: 301-306.

disp('***************************************************************');
disp('Ratio of the fluxes (standard deviations) derived from covariances to the appropriate values derived from cospectra (spectra):');
disp(['wu_cov/wu_csp = ' num2str(wu_ratio) '; wT_cov/wT_csp = ' num2str(wT_ratio) '; wv_cov/wv_csp = ' num2str(wv_ratio)]);
% disp(['wu_cov / wu_csp = ' num2str(wu_ratio)]);
% disp(['wv_cov / wv_csp = ' num2str(wv_ratio)]);
% disp(['wT_cov / wT_csp = ' num2str(wT_ratio)]);
disp(['sigu_cov/sigu_spc = ' num2str(sigu_ratio) '; sigv_cov/sigv_spc = ' num2str(sigv_ratio) '; sigw_cov/sigw_spc = ' num2str(sigw_ratio) '; sigT_cov/sigT_spc = ' num2str(sigT_ratio)]);
% disp(['sigu_cov / sigu_spc = ' num2str(sigu_ratio)]);
% disp(['sigv_cov / sigv_spc = ' num2str(sigv_ratio)]);
% disp(['sigw_cov / sigw_spc = ' num2str(sigw_ratio)]);
% disp(['sigT_cov / sigT_spc = ' num2str(sigT_ratio)]);
disp('***************************************************************');

%++++++++++++++++++++ Third-order moments of atmospheric turbulence ++++++++++++++++++++++++
% Linear detrending but no block-averaging of the data!
% nf = 32768 (54.6133 min averaged data)

uuu_dtr = (sum(us_dtr.*us_dtr.*us_dtr))./nf;
uuv_dtr = (sum(us_dtr.*us_dtr.*vs_dtr))./nf;
uuw_dtr = (sum(us_dtr.*us_dtr.*ws_dtr))./nf;
uvv_dtr = (sum(us_dtr.*vs_dtr.*vs_dtr))./nf;
uvw_dtr = (sum(us_dtr.*vs_dtr.*ws_dtr))./nf;
uww_dtr = (sum(us_dtr.*ws_dtr.*ws_dtr))./nf;
vvv_dtr = (sum(vs_dtr.*vs_dtr.*vs_dtr))./nf;
vvw_dtr = (sum(vs_dtr.*vs_dtr.*ws_dtr))./nf;
vww_dtr = (sum(vs_dtr.*ws_dtr.*ws_dtr))./nf;
www_dtr = (sum(ws_dtr.*ws_dtr.*ws_dtr))./nf;

uuT_dtr = (sum(us_dtr.*us_dtr.*Ts_dtr))./nf;
uvT_dtr = (sum(us_dtr.*vs_dtr.*Ts_dtr))./nf;
uwT_dtr = (sum(us_dtr.*ws_dtr.*Ts_dtr))./nf;
vvT_dtr = (sum(vs_dtr.*vs_dtr.*Ts_dtr))./nf;
vwT_dtr = (sum(vs_dtr.*ws_dtr.*Ts_dtr))./nf;
wwT_dtr = (sum(ws_dtr.*ws_dtr.*Ts_dtr))./nf;
uTT_dtr = (sum(us_dtr.*Ts_dtr.*Ts_dtr))./nf;
vTT_dtr = (sum(vs_dtr.*Ts_dtr.*Ts_dtr))./nf;
wTT_dtr = (sum(ws_dtr.*Ts_dtr.*Ts_dtr))./nf;
TTT_dtr = (sum(Ts_dtr.*Ts_dtr.*Ts_dtr))./nf;

%+++++++++++++++++++++++++++++++ Skewness & Kurtosis +++++++++++++++++++++++++++++++++++++++

% >>> Skewness is a measure of symmetry, or more precisely, the lack of symmetry.
% The skewness for a normal distribution is zero, and any symmetric data should have a skewness near zero.
% Negative values for the skewness indicate data that are skewed left and positive values for the skewness indicate data that are skewed right.

% Skewness of the linear detrended data:
Skew_u = skewness(us_dtr);
Skew_v = skewness(vs_dtr);
Skew_w = skewness(ws_dtr);
Skew_T = skewness(Ts_dtr);
Skew_uw = skewness(us_dtr.*ws_dtr);
Skew_vw = skewness(vs_dtr.*ws_dtr);
Skew_wT = skewness(ws_dtr.*Ts_dtr);
Skew_uT = skewness(us_dtr.*Ts_dtr);

% Direct definition of skewness (note that this is the same as 'skewness' MatLab command above):
Skew_u0=uuu_dtr./sigu_dtr.^3; Skew_v0=vvv_dtr./sigv_dtr.^3; Skew_w0=www_dtr./sigw_dtr.^3; Skew_T0=TTT_dtr./sigT_dtr.^3;
% Use this definition in the case if MatLab Statistics Toolbox is not available

% >>> Kurtosis is a measure of whether the data are peaked or flat relative to a normal distribution.
% Data sets with low kurtosis tend to have a flat top near the mean rather than a sharp peak.
% The kurtosis for a standard normal distribution is three.

% Kurtosis of the linear detrended data:
Kurt_u = kurtosis(us_dtr);
Kurt_v = kurtosis(vs_dtr);
Kurt_w = kurtosis(ws_dtr);
Kurt_T = kurtosis(Ts_dtr);
Kurt_uw = kurtosis(us_dtr.*ws_dtr);
Kurt_vw = kurtosis(vs_dtr.*ws_dtr);
Kurt_wT = kurtosis(ws_dtr.*Ts_dtr);
Kurt_uT = kurtosis(us_dtr.*Ts_dtr);

% Some forth-order moments (for kurtosis calculations):
uuuu_dtr = (sum(us_dtr.*us_dtr.*us_dtr.*us_dtr))./nf; vvvv_dtr = (sum(vs_dtr.*vs_dtr.*vs_dtr.*vs_dtr))./nf; 
wwww_dtr = (sum(ws_dtr.*ws_dtr.*ws_dtr.*ws_dtr))./nf; TTTT_dtr = (sum(Ts_dtr.*Ts_dtr.*Ts_dtr.*Ts_dtr))./nf;
% Direct definition of kurtosis (note that this is the same as 'kurtosis' MatLab command above):
Kurt_u0=uuuu_dtr./sigu_dtr.^4; Kurt_v0=vvvv_dtr./sigv_dtr.^4; Kurt_w0=wwww_dtr./sigw_dtr.^4; Kurt_T0=TTTT_dtr./sigT_dtr.^4;

% Different examples of skewness/kurtosis and PDFs in CBL and SBL can be found in:
% Chu et al (1996) WRR v32(6) 1681–1688; Graf et al (2010) BLM v134(3) 459-486; Mahrt et al (2012) JPO v42(7) 1134-1142

disp(['Skewness(w) = ' num2str(Skew_w) '; Skewness(T) = ' num2str(Skew_T) '; Skewness(u) = ' num2str(Skew_u) '; Skewness(v) = ' num2str(Skew_v)]);
disp(['Kurtosis(w) = ' num2str(Kurt_w) '; Kurtosis(T) = ' num2str(Kurt_T) '; Kurtosis(u) = ' num2str(Kurt_u) '; Kurtosis(v) = ' num2str(Kurt_v)]);
disp('***************************************************************');

disp(['Skewness(uw) = ' num2str(Skew_uw) '; Skewness(wT) = ' num2str(Skew_wT) '; Skewness(uT) = ' num2str(Skew_uT)]);
disp(['Kurtosis(uw) = ' num2str(Kurt_uw) '; Kurtosis(wT) = ' num2str(Kurt_wT) '; Kurtosis(uT) = ' num2str(Kurt_uT)]);
disp('***************************************************************');

%+++++++ Quality Control (QC) parameters derived from sonic anemometer for the 1-hr time series +++++++++++

disp('>>> Quality Control (QC) parameters for the 1-hr time series:');
disp('+++++');
% Flagged if the values of different variables exceed user-defined (user-selected) thresholds
% Angle of attack:
Phix = phi*180/pi;
disp(['Angle of attack = ' num2str(Phix) ' deg']);
% If angle of attack is large (say > 15 deg) data should be filtered out or a correction to compensate for the angle of attack error should be applied, e.g. see
% van der Molen, M. K. J. H. C. Gash, and J. A. Elbers (2004) Sonic anemometer (co)sine response and flux measurement: II. The effect of introducing an angle of attack dependent calibration. Agricultural and Forest Meterology, 122: 95-109.
% Nakai, Taro; van der Molen, M. K.; Gash, J. H. C.; Kodama, Yuji (2006) Correction of sonic anemometer angle of attack errors. Agricultural and Forest Meteorology, 136. 19-30.
disp('+++++');
% Non-Stationarity of the data in the 1-hr time series
disp(['Standard deviation of the wind speed direction = ' num2str(sigwdir_spc) ' deg']);
% If sigwdir_spc is large (e.g. > 15 deg) data may be considered as "bad" data (non-stationary)
% However this QC is not applicable for free convection limit (light winds)
% Another QC parameters: Steadiness of horizontal wind and sonic temperature (non-stationary data)
disp('+++++');
DeltaU = us_trend(:,nf) - us_trend(:,1);
DeltaV = vs_trend(:,nf) - vs_trend(:,1);
DeltaT = Ts_trend(:,nf) - Ts_trend(:,1);
disp(['Steadiness of the along-wind component (trend), delta U = ' num2str(DeltaU) ' m/s']);
disp(['Steadiness of the cross-wind component (trend), delta V = ' num2str(DeltaV) ' m/s']);
disp(['Steadiness of the sonic temperature (trend), delta T = ' num2str(DeltaT) ' deg C']);
% For example, |DeltaU| > 2 m/s, |DeltaV| > 2 m/s, |DeltaT| > 2 deg C can be used as QC thresholds

% Recall that spectral slopes and rotated mean wind components, vrm & wrm, can also be used as QC parameters (see earlier)

% Different issues of QC are discussed in:
% Foken, T. and B. Wichura. (1996) Tools for quality assessment of surface-based flux measurements. Agricultural and Forest Meteorology, 78: 83-105.
% Vickers D., Mahrt L. (1997) Quality control and flux sampling problems for tower and aircraft data. J. Atmos. Oc. Tech. 14(3): 512–526
% Foken et al. (2004) Edited by X. Lee, et al. Post-field quality control, in Handbook of micrometeorology: A guide for surface flux measurements, 81-108.
% Burba, G., and D. Anderson, (2010) A Brief Practical Guide to Eddy Covariance Flux Measurements: Principles and Workflow Examples for Scientific and Industrial Applications. LI-COR, Lincoln, USA, Hard- and Softbound, 211 pp. 

% *****************************************************************************************************************************
                % Save fluxes, spectra, cospectra, third moments etc. for the sonic anemometer as separate files
% *****************************************************************************************************************************

jday = str2num(ddd);
hour = str2num(hh);
jdate = jday+(hour/24);         % Decimal 2012 Year Day ('Julian Day')

%++++++++++++++++++++++++++++++++++ Output flux file +++++++++++++++++++++++++++++++++++++++

outputflux = [jday hour jdate sonicid wsp wdir truewdir Tsm...
    wu_csp wv_csp uv_csp wT_csp uT_csp vT_csp sigu_spc sigv_spc sigw_spc sigT_spc sigwdir_spc ...
    wu_csp_hi wv_csp_hi uv_csp_hi wT_csp_hi uT_csp_hi vT_csp_hi sigu_spc_hi sigv_spc_hi sigw_spc_hi sigT_spc_hi sigwdir_spc_hi...
    wu_cov wv_cov uv_cov wT_cov uT_cov vT_cov sigu_cov sigv_cov sigw_cov sigT_cov...
    cu2 cv2 cw2 cT2 epsilon Nt nSu nSv nSw nSt Phix DeltaU DeltaV DeltaT vrm wrm];
save(['C:\Data\AG-Raw_Data\Arctic\Tiksi\Fluxtower\Processed\Sonic\9meter\2012\Fluxes\Flux_Tiksi_sonic9m_2012_' ddd '_' hh '_raw.txt'], 'outputflux', '-ascii', '-tabs') 

%+++++++++++++++++++++++++++++++++++ Output spectra file +++++++++++++++++++++++++++++++++++

outputspectra = [fs' sus' svs' sws' sTs' swdirs'];
save(['C:\Data\AG-Raw_Data\Arctic\Tiksi\Fluxtower\Processed\Sonic\9meter\2012\Spectra\Spectra_Tiksi_sonic9m_2012_' ddd '_' hh '_raw.txt'], 'outputspectra', '-ascii', '-tabs') 

%++++++++++++++++++++++++++++++++++ Output cospectra file ++++++++++++++++++++++++++++++++++

outputcospectra = [fs' cwus' cwvs' cwTs' cuvs' cuTs' cvTs'];
save(['C:\Data\AG-Raw_Data\Arctic\Tiksi\Fluxtower\Processed\Sonic\9meter\2012\Cospectra\Cospectra_Tiksi_sonic9m_2012_' ddd '_' hh '_raw.txt'], 'outputcospectra', '-ascii', '-tabs') 

%+++++++++++++++++++++ Output Third-order moments file (+ kurtosis) ++++++++++++++++++++++++

output3rd_moments = [jday hour jdate sonicid uuu_dtr uuv_dtr uuw_dtr uvv_dtr uvw_dtr uww_dtr vvv_dtr vvw_dtr vww_dtr www_dtr ...
    uuT_dtr uvT_dtr uwT_dtr vvT_dtr vwT_dtr wwT_dtr uTT_dtr vTT_dtr wTT_dtr TTT_dtr ...
    Skew_u Skew_v Skew_w Skew_T Skew_uw Skew_vw Skew_wT Skew_uT Kurt_u Kurt_v Kurt_w Kurt_T Kurt_uw Kurt_vw Kurt_wT Kurt_uT];
save(['C:\Data\AG-Raw_Data\Arctic\Tiksi\Fluxtower\Processed\Sonic\9meter\2012\Third Moments\3rd_moments_Tiksi_sonic9m_2012_' ddd '_' hh '_raw.txt'], 'output3rd_moments', '-ascii', '-tabs') 

%+++++++++++++++++++++++++++++ Output of the 'quick products' +++++++++++++++++++++++++++++

output_quick_products = [jday hour jdate sonicid wsp truewdir Tsm ustar Hs Tstar zeta_level_n Cd phi_u phi_v phi_w phi_T phi_uT Phi_epsilon Phi_Nt...
    ustar_hi Hs_hi Tstar_hi zeta_level_n_hi Cd_hi phi_u_hi phi_v_hi phi_w_hi phi_T_hi phi_uT_hi Phi_epsilon_hi Phi_Nt_hi];
save(['C:\Data\AG-Raw_Data\Arctic\Tiksi\Fluxtower\Processed\Sonic\9meter\2012\Quick Products\quick_products_Tiksi_sonic9m_2012_' ddd '_' hh '_raw.txt'], 'output_quick_products', '-ascii', '-tabs') 

%***************************************************** The End ****************************************************************

disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>');
disp('***************************************************************');
disp(['Sonic anemometer, the Year Day: ' ddd ' , Hour: ' hh ' is done ...'])
disp('***************************************************************');
disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>');
disp('***************************************************************');
disp('The end of the "read_sonic_10Hz_1hr_Tiksi_2012_9m_v1.m" file')
disp('***************************************************************');
disp('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>');
disp(' ');
disp(' ');
% =========================================================================

