source: trunk/src/mode_sahelien/composite_pluie_eof23_1030.m @ 15

%%%%%%%%%%%%%%% FICHIER pluie IRD
% ouverture / lecture du fichier
clear;
cd('Pluie_ird_25');

% on selectionne les mois d ete (attention aux annees bixestiles!)

siz1=31+28+31+30+31+1;
siz1bis=31+28+31+30+31+30+31+31+30;
pluie=zeros(41,17,122,12);

clear v fid; fid=fopen('pluvio.79.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,1)=v(:,:,siz1:siz1bis); % 1979

clear v fid; fid=fopen('pluvio.80.haut','r','b'); siz=41*17*366;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,366);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:366;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,2)=v(:,:,siz1+1:siz1bis+1); % 1980

clear v fid; fid=fopen('pluvio.81.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,3)=v(:,:,siz1:siz1bis); % 1981

clear v fid; fid=fopen('pluvio.82.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,4)=v(:,:,siz1:siz1bis); % 1982

clear v fid; fid=fopen('pluvio.83.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;


pluie(:,:,:,5)=v(:,:,siz1:siz1bis); % 1983

clear v fid; fid=fopen('pluvio.84.haut','r','b'); siz=41*17*366;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,366);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:366;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,6)=v(:,:,siz1+1:siz1bis+1); % 1984

clear v fid; fid=fopen('pluvio.85.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,7)=v(:,:,siz1:siz1bis); % 1985

clear v fid; fid=fopen('pluvio.86.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,8)=v(:,:,siz1:siz1bis); % 1986

clear v fid; fid=fopen('pluvio.87.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,9)=v(:,:,siz1:siz1bis); % 1987

clear v fid; fid=fopen('pluvio.88.haut','r','b'); siz=41*17*366;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,366);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:366;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,10)=v(:,:,siz1:siz1bis); % 1988

clear v fid; fid=fopen('pluvio.89.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,11)=v(:,:,siz1:siz1bis); % 1989

clear v fid; fid=fopen('pluvio.90.haut','r','b'); siz=41*17*365;
v=fread(fid,siz,'float'); fclose(fid);
v=reshape(v,41,17,365);
clear a b c;
for a=1:41;
   for b=1:17;
       for c=1:365;
          if v(a,b,c)>990;
              v(a,b,c)=NaN;        
          elseif v(a,b,c)<0;
              v(a,b,c)=NaN;
          end;
       end;
   end;
end;
pluie(:,:,:,12)=v(:,:,siz1:siz1bis); % 1990

clear siz sizbis;
clear v;
cd('..');
save pluie.mat pluie;
clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load pluie.mat;
% pluie dimension : lon lat jours annee (41 17 122 12) => jour annee lat lon  (17 41 122 12)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load eof23_1030.txt


% on met sous la forme annee mois
eof23_1030=reshape(eof23_1030,122,28);
olr_eof=eof23_1030(:,1:12); % selection des annees avant 1990 
clear eof23_1030

% on calcule l ecart type sur la serie temporelle de l'olr
k=122*12; 
olrbis=reshape(olr_eof,k,1);
olrbis=std(olrbis); 

% intialisation des tableaux
pluie_filtre_compomax=NaN*ones(30,17,41); 
pluie_filtre_compomin=NaN*ones(30,17,41);
pluie=permute(pluie,[3 4 2 1]);
compteurmax=0; 
compteurmin=0;

for a=1:12; % boucle sur les annees
   clear delta; 
   delta=zeros(122,1);
   for b=1:122; % boucle sur les jours : on prepare avant   
        if b>1;
            delta(b)=olr_eof(b,a)-olr_eof(b-1,a);
        end;
   end;
   for b=1:122; % boucle sur les jours
      if (b>20 && b<112); % 1ere condition (laisser une marge de 10) 
% cas 1 : max deltas apres negetifs et deltas avant postifs   
         if (delta(b+1)<0 && delta(b)<0 && delta(b-1)>0 && delta(b-2)>0 && olr_eof(b,a)>olrbis); 
            compteurmax=compteurmax+1;
            pluie_filtre_compomax(compteurmax,:,:)=pluie(b,a,:,:);
% cas 2 : min deltas apres positifs et deltas avant negatifs   
         elseif (delta(b+1)>0 && delta(b)>0 && delta(b-1)<0 && delta(b-2)<0 && olr_eof(b,a)<(-olrbis));
            compteurmin=compteurmin+1;
            pluie_filtre_compomin(compteurmin,:,:)=pluie(b,a,:,:);
         end;
      end;
   end;         
end;


% on fait la somme des valeurs de pluie_filtreentiel sur le nbre d occurences
% on obtent un tableau de la forme (jours-avant-apres,lon,lat)
pluie_filtre_compomin1(:,:)=nanmean(pluie_filtre_compomin); clear pluie_filtre_compomin;
pluie_filtre_compomax1(:,:)=nanmean(pluie_filtre_compomax); clear pluie_filtre_compomax;

clear pluie_min pluie_max;

pluie_min=pluie_filtre_compomin1;
pluie_max=pluie_filtre_compomax1;

diff=pluie_min-pluie_max;

sauvegrads('composite_pluie_eof23',diff,[-70 2.5 -10 2.5]);

for b=1:17;
    for c=1:41;
         if diff(b,c)<-5;
             diff(b,c)=5;
          elseif diff(b,c)>5;
             diff(b,c)=5;
         end;
    end;
end
diff(1,1)=-5.01;
diff(1,2)=5.01;

palette
lat=[-10:2.5:30]; lon=[-70:2.5:30]';

contourf(lon,lat,diff,[-5:0.5:5],'LineStyle','none');
caxis([-5 5])
co=colorbar('horiz');            
set(co,'position',[0.15 0.26 0.3 0.02]);
set(co,'xtick',[-4:2:4],'xticklabel',[-4:2:4],'fontname','Arial','fontsize',4);

hold on;