%

%%%%%%%%%%%%%%% FICHIER pluie IRD
% ouverture / lecture du fichier
clear;

% 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;

status=load('eof12_rossby_30100.txt');
eof12_rossby_30100=status;
clear status;

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

% 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;

% boucle sur les annees
for a=1:12;
   clear delta;
   delta=zeros(122,1);
   % boucle sur les jours : on prepare avant
   for b=1:122;
        if b>1;
            delta(b)=olr_eof(b,a)-olr_eof(b-1,a);
        end;
   end;
   % boucle sur les jours
   for b=1:122;
      % 1ere condition (laisser une marge de 10)
      if (b>20 && b<112);
% 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_ird_30100',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;