% clamp_public_means.m
% means and stds of clamp AWS data 
% all variables are vectors 1 by x, with time sd(x) and location loc(x) 
% variables are: cl lA fname ...
%	sd ws wd rh0 rh1 t0 t1 p volt mslp loc ...
%	u v theta0 theta1 q0 q1 rhi0 rhi1 ...
%	cheight cslope clats clons sdmon
%       shfx lhfx

load C:/public/data/clamp/matlab_data/clamp_aws.mat


% set up mean variables as 5 by 4 arrays of NaN
length_data = NaN*ones(5,4);  
mean_ws = NaN*ones(5,4); vectormean_ws = NaN*ones(5,4); 
vectormean_wd = NaN*ones(5,4); mean_dir_constancy = NaN*ones(5,4); 
mean_t0 = NaN*ones(5,4); mean_t1 = NaN*ones(5,4);
mean_rhi0 = NaN*ones(5,4); mean_rhi1 = NaN*ones(5,4); 
mean_theta0 = NaN*ones(5,4); mean_theta1 = NaN*ones(5,4); 
mean_p = NaN*ones(5,4); mean_mslp = NaN*ones(5,4);
mean_shfx = NaN*ones(5,4); mean_lhfx = NaN*ones(5,4);

std_t1 = NaN*ones(5,4);  std_t1_t0 = NaN*ones(5,4);
std_ws = NaN*ones(5,4);  std_p = NaN*ones(5,4);
std_rhi1 = NaN*ones(5,4); 
std_shfx = NaN*ones(5,4); std_lhfx = NaN*ones(5,4);

% find certain periods and locations
% i=location (clamp 1 to 4) plus 5 = halley, 6 = ago77 
% j = season (1=djf, 2=mam, 3=jja, 4=son)
% where filternan used, checked that this does not reduce data set too much

stations = char('Halley', 'clamp 1','clamp 2','clamp 3','clamp 4','ago 77');
seasons = char('DJF','MAM','JJA','SON');
mons = [12 1 2; 3 4 5; 6 7 8; 9 10 11];

for i = 1:5
  for j = 1:4
   e=find(((sdmon==mons(j,1))|(sdmon==mons(j,2))|(sdmon==mons(j,3))) ...
	&(loc==i));
   if (length(e) ~= 0)
    length_data(i,j) = length(e);
    mean_ws(i,j) = mean(filternan(ws(e)));
    [vectormean_ws(i,j) vectormean_wd(i,j)] = ...
		 uv_to_vectormean(filternan(u(e)),filternan(v(e)));
    mean_dir_constancy(i,j) = dir_constancy(u(e),v(e));
    mean_t0(i,j) = mean(filternan(t0(e)));
    mean_t1(i,j) = mean(filternan(t1(e)));
    mean_rhi0(i,j) = mean(filternan(rhi0(e)));
    mean_rhi1(i,j) = mean(filternan(rhi1(e)));
    mean_theta0(i,j) = mean(filternan(theta0(e)));
    mean_theta1(i,j) = mean(filternan(theta1(e)));
    mean_q0(i,j) = mean(filternan(q0(e)));
    mean_q1(i,j) = mean(filternan(q1(e)));
    mean_p(i,j) =  mean(filternan(p(e)));
    mean_dthetadz(i,j) =  mean(filternan(dthetadz(e)));
    mean_shfx(i,j) =  mean(filternan(shfx(e)));
    mean_lhfx(i,j) =  mean(filternan(lhfx(e)));
    if (i==5) mean_cloud(i,j)= mean(filternan(cloud(e))); end;

    std_t1(i,j) = std(filternan(t1(e)));
    std_t1_t0(i,j) = std(filternan(t1(e)-t0(e)));
    std_ws(i,j) = std(filternan(ws(e)));
    std_q1(i,j) = std(filternan(q1(e)));
    std_p(i,j) =  std(filternan(p(e)));
    std_rhi1(i,j) =  std(filternan(rhi1(e)));
    std_dthetadz(i,j) =  std(filternan(dthetadz(e)));
    std_shfx(i,j) =  std(filternan(shfx(e)));
    if ((i~=3) | (j~=4)) std_lhfx(i,j) = std(filternan(lhfx(e))); end;   
    if (i==5) std_cloud(i,j)= std(filternan(cloud(e))); end;
   end;	%if loop

  end;
end;



