function [out,stdev,samples] = climatology(in,time,type)
% computes mean climatology 'out', standard dev 'stdev', and number of samples 'samples'
% of vector 'in' with matlab time vector 'time'
% type specifies over which timeinterval to average: 1: yearly, 2: monthly, 3: daily, 4: hourly, 5: hourly by month
% all NaN values are omitted from dataset

time_vec = datevec(time);
if type == 1        % yearly
    min = time_vec(1,1); max=time_vec(length(time),1);
elseif type == 2    % monthly
    min = 1; max = 12;
elseif type == 3    % daily
    min = 1; max = 365;
elseif type == 4    % hourly
    min = 0; max = 23;
end

if type < 4
for i = min:max
out(i) = mean( in( find(time_vec(:,type)==i & ~isnan(in)) ) );
stdev(i) = std( in( find(time_vec(:,type)==i & ~isnan(in)) ) );
samples(i) = length( in( find(time_vec(:,type)==i & ~isnan(in)) ) );
end

elseif type == 4
for i = min:max
out(i+1) = mean( in( find(time_vec(:,type)==i & ~isnan(in)) ) );
stdev(i+1) = std( in( find(time_vec(:,type)==i & ~isnan(in)) ) );
samples(i+1) = length( in( find(time_vec(:,type)==i & ~isnan(in)) ) );
end

elseif type == 5
    for i = 1:12
        for j=0:23
            out(i,j+1) = mean(  in( find( time_vec(:,2)==i & time_vec(:,4)==j & ~isnan(in)) ) );
          stdev(i,j+1) =  std(  in( find( time_vec(:,2)==i & time_vec(:,4)==j & ~isnan(in)) ) );
        samples(i,j+1) =length( in( find( time_vec(:,2)==i & time_vec(:,4)==j & ~isnan(in)) ) );
        end
    end
end