%Makes climatology of data according to the month
%Input has to have time as first dimension
%The output has same dimension of data_in
%data_mon is format 20101010
%method is 'clim' or something else
%method2 is 'daily' or 'season'
%FORMAT: [data_new,monthnum,std_new] = make_clim(data_in,data_mon,method,method2)

function [data_new,monthnum,std_new] = make_clim(data_in,data_mon,method,method2)
addpath /data/mgoes/AX8/Time_transform

if nargin < 3, method = [];method2=[];end
if nargin < 4, method2 = [];end

if size(data_in,1)~=size(data_mon,1); error('First dimension is time');end

%ndays = [31 29 31 30 31 30 31 31 30 31 30 31];
%NOW ADD DAILY
[nmonths,incday] = size(data_mon);
allmonths = nan*ones(1,nmonths);
alldays   = nan*ones(1,nmonths);                         %DAILY
for i=1:nmonths
    aux_mon = str2num(data_mon(i,5:6));
    if isempty(aux_mon),continue;end
    allmonths(i) = aux_mon;%str2num(data_mon(i,5:6));
    alldays(i) = str2num(data_mon(i,5:8));                %DAILY
end

if strcmp(method,'clim')
 for i =1:12
 acha = find(allmonths == i);
 if strcmp(method2,'season')
    ii = [i-1 i i+1]; ii(ii<=0) = 12+ii(ii<=0); ii(ii>12) = ii(ii>12) -12;
  acha = find(allmonths == ii(1) | allmonths == ii(2)...
      | allmonths == ii(3));
 end 
 monthnum(i) = i;
 data_aux = nanmean(data_in(acha,:,:),1);
 data_new(i,:,:) = data_aux(1,:,:);
 std_aux = nanstd(data_in(acha,:,:),[],1);
 std_new(i,:,:) = std_aux(1,:,:);             %standard deviation
 %std_new(i,:,:) = std_aux(1,:,:)/sqrt(length(acha));%standard error
 end
return 
end

% if strcmp(method,'season')
%  for i =1:12
%      ii = [i-1 i i+1]; ii(ii<=0) = 12+ii(ii<=0); ii(ii>12) = ii(ii>12) -12;
%  acha = find(allmonths == ii(1) | allmonths == ii(2) | allmonths == ii(3));
%  monthnum(i) = i;
%  data_aux = nanmean(data_in(acha,:,:),1);
%  data_new(i,:,:) = data_aux(1,:,:);
%  std_aux = nanstd(data_in(acha,:,:),[],1);
%  std_new(i,:,:) = std_aux(1,:,:);             %standard deviation
%  %std_new(i,:,:) = std_aux(1,:,:)/sqrt(length(acha));%standard error
%  end
% return 
% end



data_new = zeros(size(data_in));
monthnum = zeros(1,incday);
std_new  = nan*ones(size(data_in));

%DAILY%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(method2,'daily')
    
    day366 = unique(alldays);
for i=1:length(day366)
    acha = find(alldays==day366(i));
    data_aux = nanmean(data_in(acha,:,:),1);
 data_new(acha,:,:) = data_aux(ones(1,length(acha)),:,:);
 std_aux = nanstd(data_in(acha,:,:),[],1);
 std_new(acha,:,:) = std_aux(ones(1,length(acha)),:,:);
end

    
else
%MONTHLY
for i =1:12
 acha = find(allmonths == i);
   monthnum(acha) = i;
if strcmp(method2,'season')
    ii = [i-1 i i+1]; ii(ii<=0) = 12+ii(ii<=0); ii(ii>12) = ii(ii>12) -12;
  acha = find(allmonths == ii(1) | allmonths == ii(2)...
      | allmonths == ii(3));
 end 
 data_aux = nanmean(data_in(acha,:,:),1);
 data_new(acha,:,:) = data_aux(ones(1,length(acha)),:,:);
 std_aux = nanstd(data_in(acha,:,:),[],1);
 std_new(acha,:,:) = std_aux(ones(1,length(acha)),:,:);
end

end %daily



%year =  str2double(data_mon(i,1:4));
%month = str2double(data_mon(i,5:6));
% if incday>6
%    day   = str2double(data_mon(i,7:8));
% end
%
%frac_month = month/13;
%frac_month_p = (month+1)/13;
%frac_day = (frac_month_p-frac_month)*day/(ndays(month)+1);
%new_data(i)=year + frac_month + frac_day;
%end
