source: trunk/src/SIMULS_IRCAAM/composite_hgt500_TROP_eof2_1030.m @ 23

Last change on this file since 23 was 23, checked in by pinsard, 15 years ago

replace fcommands save, load and print by function calls
File size: 10.5 KB
Line 
1%
2
3% initialisation
4clear;
5close all;
6
7status=load('hgt500_trop.mat');
8
9ncload('zg_d.TrNQIVIV.197106-09.nc');
10
11% selction du domaine geographique
12hgt500b=zeros(122,29,64,128);
13hgt500b(:,:,:,1:64)=hgt500(:,:,:,65:128);
14hgt500b(:,:,:,65:128)=hgt500(:,:,:,1:64);
15clear hgt500;
16hgt500=hgt500b;
17clear hgt500b;
18
19% on desaisonalise
20hgt1=permute(hgt500,[2 1 3 4]);
21hgt2(:,:,:)=nanmean(hgt1);
22clear hgt1;
23hgt3=reshape(hgt500, 3538,64,128);
24hgt4(:,:)=nanmean(hgt3);
25clear hgt3;
26
27for a=1:122;
28   for b=1:28;
29      for c=1:37;
30         for d=1:128;
31             hgt500(a,b,c,d)=hgt500(a,b,c,d)-(hgt2(a,c,d)-hgt4(c,d));
32         end;
33      end;
34   end;
35end;
36clear hgt2 hgt4;
37
38status=load('eof2_TR_1030.txt');
39
40% on met sous la forme annee mois
41hgt500_eof=reshape(eof2_TR_1030,122,29);
42clear eof2_TR_1030;
43
44% on calcule l ecart type sur la serie temporelle de l HGT500
45k=122*29;
46hgt500bis=reshape(hgt500_eof,k,1);
47hgt500bis=std(hgt500bis);
48
49% intialisation des tableaux
50hgt500_filtre_compomax=NaN*ones(50,31,64,128);
51hgt500_filtre_compomin=NaN*ones(50,31,64,128);
52
53compteurmax=0;
54compteurmin=0;
55
56% boucle sur les annees
57for a=1:29;
58   clear delta;
59   delta=zeros(122,1);
60   % boucle sur les jours : on prepare avant
61   for b=1:122;
62        if b>1;
63            delta(b)=hgt500_eof(b,a)-hgt500_eof(b-1,a);
64        end;
65   end;
66   % boucle sur les jours
67   for b=1:122;
68      % 1ere condition (laisser une marge de 10)
69      if (b>20 && b<112);
70% cas 1 : max deltas apres negetifs et deltas avant postifs
71         if (delta(b+1)<0 && delta(b)>0 && hgt500_eof(b,a)>(hgt500bis));
72            compteurmax=compteurmax+1;
73            hgt500_filtre_compomax(compteurmax,1,:,:)=hgt500(b-20,a,:,:);
74            hgt500_filtre_compomax(compteurmax,2,:,:)=hgt500(b-19,a,:,:);
75            hgt500_filtre_compomax(compteurmax,3,:,:)=hgt500(b-18,a,:,:);
76            hgt500_filtre_compomax(compteurmax,4,:,:)=hgt500(b-17,a,:,:);
77            hgt500_filtre_compomax(compteurmax,5,:,:)=hgt500(b-16,a,:,:);
78            hgt500_filtre_compomax(compteurmax,6,:,:)=hgt500(b-15,a,:,:);
79            hgt500_filtre_compomax(compteurmax,7,:,:)=hgt500(b-14,a,:,:);
80            hgt500_filtre_compomax(compteurmax,8,:,:)=hgt500(b-13,a,:,:);
81            hgt500_filtre_compomax(compteurmax,9,:,:)=hgt500(b-12,a,:,:);
82            hgt500_filtre_compomax(compteurmax,10,:,:)=hgt500(b-11,a,:,:);
83
84            hgt500_filtre_compomax(compteurmax,11,:,:)=hgt500(b-10,a,:,:);
85            hgt500_filtre_compomax(compteurmax,12,:,:)=hgt500(b-9,a,:,:);
86            hgt500_filtre_compomax(compteurmax,13,:,:)=hgt500(b-8,a,:,:);
87            hgt500_filtre_compomax(compteurmax,14,:,:)=hgt500(b-7,a,:,:);
88            hgt500_filtre_compomax(compteurmax,15,:,:)=hgt500(b-6,a,:,:);
89            hgt500_filtre_compomax(compteurmax,16,:,:)=hgt500(b-5,a,:,:);
90            hgt500_filtre_compomax(compteurmax,17,:,:)=hgt500(b-4,a,:,:);
91            hgt500_filtre_compomax(compteurmax,18,:,:)=hgt500(b-3,a,:,:);
92            hgt500_filtre_compomax(compteurmax,19,:,:)=hgt500(b-2,a,:,:);
93            hgt500_filtre_compomax(compteurmax,20,:,:)=hgt500(b-1,a,:,:);
94
95            hgt500_filtre_compomax(compteurmax,21,:,:)=hgt500(b,a,:,:);
96            hgt500_filtre_compomax(compteurmax,22,:,:)=hgt500(b+1,a,:,:);
97            hgt500_filtre_compomax(compteurmax,23,:,:)=hgt500(b+2,a,:,:);
98            hgt500_filtre_compomax(compteurmax,24,:,:)=hgt500(b+3,a,:,:);
99            hgt500_filtre_compomax(compteurmax,25,:,:)=hgt500(b+4,a,:,:);
100            hgt500_filtre_compomax(compteurmax,26,:,:)=hgt500(b+5,a,:,:);
101            hgt500_filtre_compomax(compteurmax,27,:,:)=hgt500(b+6,a,:,:);
102            hgt500_filtre_compomax(compteurmax,28,:,:)=hgt500(b+7,a,:,:);
103            hgt500_filtre_compomax(compteurmax,29,:,:)=hgt500(b+8,a,:,:);
104            hgt500_filtre_compomax(compteurmax,30,:,:)=hgt500(b+9,a,:,:);
105            hgt500_filtre_compomax(compteurmax,31,:,:)=hgt500(b+10,a,:,:);
106% cas 2 : min deltas apres positifs et deltas avant negatifs
107         elseif (delta(b+1)>0 && delta(b)<0 && hgt500_eof(b,a)<(hgt500bis));
108            compteurmin=compteurmin+1;
109
110            hgt500_filtre_compomin(compteurmin,1,:,:)=hgt500(b-20,a,:,:);
111            hgt500_filtre_compomin(compteurmin,2,:,:)=hgt500(b-19,a,:,:);
112            hgt500_filtre_compomin(compteurmin,3,:,:)=hgt500(b-18,a,:,:);
113            hgt500_filtre_compomin(compteurmin,4,:,:)=hgt500(b-17,a,:,:);
114            hgt500_filtre_compomin(compteurmin,5,:,:)=hgt500(b-16,a,:,:);
115            hgt500_filtre_compomin(compteurmin,6,:,:)=hgt500(b-15,a,:,:);
116            hgt500_filtre_compomin(compteurmin,7,:,:)=hgt500(b-14,a,:,:);
117            hgt500_filtre_compomin(compteurmin,8,:,:)=hgt500(b-13,a,:,:);
118            hgt500_filtre_compomin(compteurmin,9,:,:)=hgt500(b-12,a,:,:);
119            hgt500_filtre_compomin(compteurmin,10,:,:)=hgt500(b-11,a,:,:);
120
121            hgt500_filtre_compomin(compteurmin,11,:,:)=hgt500(b-10,a,:,:);
122            hgt500_filtre_compomin(compteurmin,12,:,:)=hgt500(b-9,a,:,:);
123            hgt500_filtre_compomin(compteurmin,13,:,:)=hgt500(b-8,a,:,:);
124            hgt500_filtre_compomin(compteurmin,14,:,:)=hgt500(b-7,a,:,:);
125            hgt500_filtre_compomin(compteurmin,15,:,:)=hgt500(b-6,a,:,:);
126            hgt500_filtre_compomin(compteurmin,16,:,:)=hgt500(b-5,a,:,:);
127            hgt500_filtre_compomin(compteurmin,17,:,:)=hgt500(b-4,a,:,:);
128            hgt500_filtre_compomin(compteurmin,18,:,:)=hgt500(b-3,a,:,:);
129            hgt500_filtre_compomin(compteurmin,19,:,:)=hgt500(b-2,a,:,:);
130            hgt500_filtre_compomin(compteurmin,20,:,:)=hgt500(b-1,a,:,:);
131
132            hgt500_filtre_compomin(compteurmin,21,:,:)=hgt500(b,a,:,:);
133            hgt500_filtre_compomin(compteurmin,22,:,:)=hgt500(b+1,a,:,:);
134            hgt500_filtre_compomin(compteurmin,23,:,:)=hgt500(b+2,a,:,:);
135            hgt500_filtre_compomin(compteurmin,24,:,:)=hgt500(b+3,a,:,:);
136            hgt500_filtre_compomin(compteurmin,25,:,:)=hgt500(b+4,a,:,:);
137            hgt500_filtre_compomin(compteurmin,26,:,:)=hgt500(b+5,a,:,:);
138            hgt500_filtre_compomin(compteurmin,27,:,:)=hgt500(b+6,a,:,:);
139            hgt500_filtre_compomin(compteurmin,28,:,:)=hgt500(b+7,a,:,:);
140            hgt500_filtre_compomin(compteurmin,29,:,:)=hgt500(b+8,a,:,:);
141            hgt500_filtre_compomin(compteurmin,30,:,:)=hgt500(b+9,a,:,:);
142            hgt500_filtre_compomin(compteurmin,31,:,:)=hgt500(b+10,a,:,:);
143         end;
144      end;
145   end;
146end;
147hgt500_filtre_compomin=hgt500_filtre_compomin(1:compteurmin,:,:,:);
148hgt500_filtre_compomax=hgt500_filtre_compomax(1:compteurmax,:,:,:);
149
150% on fait la somme des valeurs de hgt500_filtreentiel sur le nbre d occurences
151% on obtent un tableau de la forme (jours-avant-apres,lon,lat)
152hgt500_filtre_compomin1(:,:,:)=mean(hgt500_filtre_compomin);
153clear hgt500_filtre_compomin;
154hgt500_filtre_compomax1(:,:,:)=mean(hgt500_filtre_compomax);
155clear hgt500_filtre_compomax;
156
157x=[-10:1:10]';
158
159figure(1);
160orient('landscape');
161
162lat=lat';
163lon=lon-180;
164for e=1:7;
165       a=7-(e-1)+1;
166       f=2*(a-1)+1;
167       clear delta_hgt500_filtre_min;
168       delta_hgt500_filtre_min=zeros(64,128);
169       for b=1:64;
170          for c=1:128;
171             delta_hgt500_filtre_min(b,c)=hgt500_filtre_compomin1(f,b,c)-hgt500_filtre_compomax1(f,b,c);
172          end;
173       end;
174       for b=1:64;
175          for c=1:128;
176             if delta_hgt500_filtre_min(b,c)<-50;
177                  delta_hgt500_filtre_min(b,c)=-50;
178             elseif delta_hgt500_filtre_min(b,c)>50;
179                  delta_hgt500_filtre_min(b,c)=50;
180             end;
181          end;
182       end
183       subplot(7,1,8-e);
184       palettecomplet
185       delta_hgt500_filtre_min(1,1)=-50.0001;
186       delta_hgt500_filtre_min(1,2)=50.0001;
187       fin=-30+24*2.5;
188
189       clear fin
190       fin=-10+16*2.5;
191       aa=contourf(lon,lat,delta_hgt500_filtre_min,[-50:5:50],'LineStyle','none');
192       if e==1;
193       set (gca,'XTick',[-180:60:180],'XTickLabel',[' 180';'120W';' 60W';'   0';' 60E';'120E';' 180'],'fontname','Arial','fontsize',6);
194       else
195      set (gca,'XTick',[-180:60:180],'XTickLabel',[' ';' ';' ';' ';' ';' '],'fontname','Arial','fontsize',6);
196       end
197       set (gca,'YTick',[30:10:70],'YTickLabel',['30N';'40N';'50N';'60N';'70N'],'fontname','Arial','fontsize',6);
198       if e==4;
199          co=colorbar;
200          set (co,'xtick',[-20:10:20],'xticklabel',[-20:10:20],'fontname','Arial','fontsize',6);
201       end;
202       hold on;
203       cartemonde1;
204       set (gca,'tickDir','out');
205       axis equal
206       axis([-180 180 20 80]);
207       if e==7;
208          title('HGT500 composite EOF2 1030j simulation TROP','fontsize',8) ;
209       end
210       text(-270,10,['day = ',num2str(f-21)],'fontname','Arial','fontsize',10);
211       pos_vert=e/10;
212       set (gca,'position',[0.2 pos_vert 0.5 0.1])
213end;
214status=print('-depsc2','composite_hgt500_trop_eof2_1030a.eps');
215
216figure(2);
217orient('landscape');
218
219for e=1:7;
220       a=7-(e-1)+1+7;
221       f=2*(a-1)+1;
222       clear delta_hgt500_filtre_min;
223       delta_hgt500_filtre_min=zeros(64,128);
224       for b=1:64;
225          for c=1:128;
226             delta_hgt500_filtre_min(b,c)=hgt500_filtre_compomin1(f,b,c)-hgt500_filtre_compomax1(f,b,c);
227          end;
228       end;
229       for b=1:64;
230          for c=1:128;
231             if delta_hgt500_filtre_min(b,c)<-50;
232                  delta_hgt500_filtre_min(b,c)=-50;
233             elseif delta_hgt500_filtre_min(b,c)>50;
234                  delta_hgt500_filtre_min(b,c)=50;
235             end;
236          end;
237       end
238       subplot(7,1,8-e);
239       palettecomplet
240       delta_hgt500_filtre_min(1,1)=-50.0001;
241       delta_hgt500_filtre_min(1,2)=50.0001;
242       fin=-30+24*2.5;
243
244       clear fin
245       fin=-10+16*2.5;
246       aa=contourf(lon,lat,delta_hgt500_filtre_min,[-50:5:50],'LineStyle','none');
247       if e==1;
248       set (gca,'XTick',[-180:60:180],'XTickLabel',[' 180';'120W';' 60W';'   0';' 60E';'120E';' 180'],'fontname','Arial','fontsize',6);
249       else
250      set (gca,'XTick',[-180:60:180],'XTickLabel',[' ';' ';' ';' ';' ';' '],'fontname','Arial','fontsize',6);
251       end
252       set (gca,'YTick',[30:10:70],'YTickLabel',['30N';'40N';'50N';'60N';'70N'],'fontname','Arial','fontsize',6);
253       if e==4;
254          co=colorbar;
255          set (co,'xtick',[-20:10:20],'xticklabel',[-20:10:20],'fontname','Arial','fontsize',6);
256       end;
257       hold on;
258       cartemonde1;
259       set (gca,'tickDir','out');
260       axis equal
261       axis([-180 180 20 80]);
262       if e==7;
263          title('HGT500 composite EOF2 1030j simulation TROP','fontsize',8) ;
264       end
265       text(-270,10,['day = ',num2str(f-21)],'fontname','Arial','fontsize',10);
266       pos_vert=e/10;
267       set (gca,'position',[0.2 pos_vert 0.5 0.1])
268end;
269status=print('-depsc2','composite_hgt500_trop_eof2_1030b.eps');
Note: See TracBrowser for help on using the repository browser.