% This routine loads equidistant 5-Hz data of w (vertical velocity) and T
% (temperature) and calculates the sensible heat flux.
%
% Author:   Andreas Muschinski
% Date:     1 November 2009
%
%-------------------------------------------------------------------------

% Load 5-Hz data from 16 September:
load wT_5Hz_LUKE_20090916.mat;

% Pressure: 
p = 100000;

% Specific heat capacity of dry air at constant pressure:
c_p = 1004;

% Specific gas constant for dry air:
R_d = 287;

%------------------------------------------------------------------------
% 30-min estimates:
N_avg = 48;
N_block = 30 * 60* 5;

for i=1:N_avg
    t_30min(i) = mean(t_5Hz( (i-1)*N_block+1:i*N_block ) );
    T_mean  = mean(T_5Hz( (i-1)*N_block+1:i*N_block ) );
    w_mean  = mean(w_5Hz( (i-1)*N_block+1:i*N_block ) );
    H_30min(i) = c_p*p/(R_d*T_mean) * mean ( (T_5Hz((i-1)*N_block+1:i*N_block)-T_mean).* (w_5Hz((i-1)*N_block+1:i*N_block)-w_mean) );
end

%------------------------------------------------------------------------
% 10-min estimates:
N_avg = 144;
N_block = 10 * 60* 5;

for i=1:N_avg
    t_10min(i) = mean(t_5Hz( (i-1)*N_block+1:i*N_block ) );
    T_mean  = mean(T_5Hz( (i-1)*N_block+1:i*N_block ) );
    w_mean  = mean(w_5Hz( (i-1)*N_block+1:i*N_block ) );
    H_10min(i) = c_p*p/(R_d*T_mean) * mean ( (T_5Hz((i-1)*N_block+1:i*N_block)-T_mean).* (w_5Hz((i-1)*N_block+1:i*N_block)-w_mean) );
end

%------------------------------------------------------------------------
% 1-min estimates:
N_avg = 1440;
N_block = 1 * 60* 5;

for i=1:N_avg
    t_1min(i) = mean(t_5Hz( (i-1)*N_block+1:i*N_block ) );
    T_mean  = mean(T_5Hz( (i-1)*N_block+1:i*N_block ) );
    w_mean  = mean(w_5Hz( (i-1)*N_block+1:i*N_block ) );
    H_1min(i) = c_p*p/(R_d*T_mean) * mean ( (T_5Hz((i-1)*N_block+1:i*N_block)-T_mean).* (w_5Hz((i-1)*N_block+1:i*N_block)-w_mean) );
end


figure(1)
plot(t_1min-4/24,H_1min,'b')
hold on
plot(t_10min-4/24,H_10min,'ro-','MarkerSize',3)
hold on
plot([t_5Hz(1)-4/24 t_5Hz(end)-4/24],[0 0],'k')
hold off
datetick('x',13)
ylim([-100 400])
xlabel('Eastern Standard Time [hhmmss]')
ylabel('H [W m^{-2}]')
title('Sensible heat flux, Amherst, 20090916')
legend('1 min','10 min')

figure(2)
plot(t_10min-4/24,H_10min,'r')
hold on
plot(t_30min-4/24,H_30min,'ko-','MarkerSize',3)
hold on
plot([t_5Hz(1)-4/24 t_5Hz(end)-4/24],[0 0],'k')
hold off
datetick('x',13)
ylim([-100 400])
xlabel('Eastern Standard Time [hhmmss]')
ylabel('H [W m^{-2}]')
title('Sensible heat flux, Amherst, 20090916')
legend('10 min','30 min')