/[MITgcm]/MITgcm_contrib/gmaze_pv/visu/eg_view_Timeserie_pl1.m
ViewVC logotype

Contents of /MITgcm_contrib/gmaze_pv/visu/eg_view_Timeserie_pl1.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Fri Oct 6 19:02:10 2006 UTC (19 years, 2 months ago) by gmaze
Branch: MAIN
CVS Tags: HEAD
update visu

1 %DEF 1 var per figure
2
3 % Map projection:
4 m_proj('mercator','long',subdomain.limlon,'lat',subdomain.limlat);
5 %m_proj('mercator','long',subdomain2.limlon,'lat',subdomain2.limlat);
6 %m_proj('mercator','long',subdomain.limlon,'lat',[25 40]);
7 %m_proj('mercator','long',[subdomain.limlon(1) 360-24],'lat',[25 50]);
8
9 % Which variables to plot:
10 wvar = [21 10 7 22];
11 wvar = [12];
12
13 for ip = 1 : length(wvar)
14
15 figur(10+ip);clf;drawnow;hold on;iw=1;jw=1;
16
17 % Variables loop:
18 % Default:
19 CBAR = 'v'; % Colorbar orientation
20 Tcontour = [17 19]; Tcolor = [0 0 0]; % Theta contours
21 Hcontour = -[0:200:600]; Hcolor = [0 0 0]; % MLD contours
22 unit = ''; % Default unit
23 load mapanom2 ; N = 256; c = [0 0]; cmap = jet; % Colormaping
24 uselogmap = 1; % Use the log colormap
25 showT = 1; % Iso-Theta contours
26 showW = 0; % Windstress arrows
27 showH = 0; colorW = 'w'; % Mixed Layer Depth
28 showE = 0; colorE = 'w'; % Ekman Layer Depth
29 showCLIM = 1; % Show CLIMODE region box
30 showQnet = 0 ; Qnetcontour = [-1000:100:1000]; % Show the Net heat flux
31 CONT = 0; % Contours instead of pcolor
32 CONTshlab = 0; % Show label for contours plot
33 CONTc = 0; % Highlighted contours instead of pcolor
34 CONTcshlab = 0; % Show label for contours plot
35 colorCOAST = [0 0 0]; % Land color
36 SHADE = 'flat'; % shading option
37
38 %if it==1, mini; end
39 switch wvar(ip)
40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1
41 case 1
42 C = Diso;
43 Clon = Qlon; Clat = Qlat;
44 tit = strcat('Depth of \sigma_\theta=',num2str(iso),'kg.m^{-3}');
45 showW = 0; % Windstress
46 cx = [-600 0];
47 unit = 'm';
48 titf = 'Diso';
49 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 11
50 case 11
51 C = Qiso;
52 %C(isnan(C)) = 10;
53 Clon = Qlon; Clat = Qlat;
54 tit = strcat('Full potential vorticity field on iso-\sigma=',num2str(iso));
55 colorCOAST = [1 1 1]*.5; % Land color
56 showW = 0; % Windstress
57 showH = 0;
58 showCLIM = 1;
59 showT = 1;
60 Tcontour = [22 22];
61 N = 256;
62 c = [1 5];
63 cx = [-(1+2*c(1)) 1+2*c(2)]*5e-10; cmap = mapanom; %cmap = jet;
64 unit = '1/m/s';
65 titf = strcat('PViso_',num2str(iso));
66 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 12
67 case 12
68 % First iso-ST have been computed in view_Timeserie
69 % Here is the 2nd one (supposed to be deeper than the 1st one)
70 iso2 = 25.35;
71 disp('Get 2nd iso-ST')
72 [Iiso mask] = subfct_getisoS(ST,iso2);
73 Diso1 = Diso;
74 Diso2 = ones(size(Iiso)).*NaN;
75 Qiso2 = ones(size(Iiso)).*NaN;
76 for ix = 1 : size(ST,3)
77 for iy = 1 : size(ST,2)
78 if ~isnan(Iiso(iy,ix)) & ~isnan( Q(Iiso(iy,ix),iy,ix) )
79 Diso2(iy,ix) = STdpt(Iiso(iy,ix));
80 Qiso2(iy,ix) = Q(Iiso(iy,ix),iy,ix);
81 end %if
82 end, end %for iy, ix
83 Diso1(isnan(squeeze(MASK(1,:,:)))) = NaN;
84 Diso2(isnan(squeeze(MASK(1,:,:)))) = NaN;
85 for ix = 1 : size(ST,3)
86 for iy = 1 : size(ST,2)
87 if isnan(Diso1(iy,ix)) & isnan(Diso2(iy,ix))
88 Hbiso(iy,ix) = -Inf;
89 elseif isnan(Diso1(iy,ix)) & ~isnan(Diso2(iy,ix))
90 Hbiso(iy,ix) = Inf;
91 elseif ~isnan(Diso1(iy,ix)) & ~isnan(Diso2(iy,ix))
92 Hbiso(iy,ix) = Diso1(iy,ix) - Diso2(iy,ix);
93 end
94 end, end %for iy, ix
95 Hbiso = Hbiso.*squeeze(MASK(1,:,:));
96 %figur(1);pcolor(Diso1);shading flat;colorbar
97 %figur(2);pcolor(Diso2);shading flat;colorbar
98 %figur(3);pcolor(Hbiso);shading flat;colorbar
99
100 C = Hbiso;
101 %C(isnan(C)) = NaN;
102 Clon = STlon; Clat = STlat;
103 tit = strvcat(strcat('Height between iso-\sigma=',num2str(iso),...
104 ' and: iso-\sigma=',num2str(iso2)),...
105 strcat('(White areas are outcrops for both iso-\sigma and red areas only for: ',...
106 num2str(iso),')'));
107 colorCOAST = [1 1 0]*.8; % Land color
108 showW = 1; colorW = 'k'; % Windstress
109 showH = 0; % Mixed layer depth
110 showCLIM = 1; % CLIMODE
111 showT = 0; Tcontour = [21 23]; % THETA
112 CONTc = 0; CONTcv = [0:0]; CONTcshlab = 1;
113 showQnet = 1; Qnetcontour=[[-1000:200:-400] [400:200:1000]];
114 cx = [0 300]; uselogmap = 0;
115 cmap = [[1 1 1]; jet ; [1 0 0]];
116 unit = 'm';
117 titf = strcat('Hbiso_',num2str(iso),'_',num2str(iso2));
118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2
119 case 2
120 C = QisoN; % C = Qiso;
121 Clon = Qlon; Clat = Qlat;
122 tit = strcat(snapshot,'/ Potential vorticity field: q = (-f/\rho . d\sigma_\theta/dz) / q_{ref}');
123 %tit = strcat(snapshot,'/ Potential vorticity field: q = - f . d\sigma_\theta/dz / \rho');
124 showW = 0; % Windstress
125 cx = [0 1]*10;
126 unit = char(1);
127 titf = 'PVisoN';
128
129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 21
130 case 21
131 C = squeeze(Q(1,:,:));
132 Clon = Qlon; Clat = Qlat;
133 tit = strcat('Surface potential vorticity field');
134 showW = 0; % Windstress
135 showH = 0;
136 showCLIM = 1;
137 N = 256;
138 c = [1 12]; cx = [-(1+2*c(1)) 1+2*c(2)]*1e-14; cmap = mapanom; % ecco2_bin1
139 c = [1 12]; cx = [-(1+2*c(1)) 1+2*c(2)]*1e-11; cmap = mapanom; % ecco2_bin2
140 unit = '1/m/s';
141 titf = 'PV.Lsurface';
142 %SHADE = 'interp';
143
144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 22
145 case 22
146 C = squeeze(Q(11,:,:));
147 Clon = Qlon; Clat = Qlat;
148 tit = strcat('Full potential vorticity field (-115m)');
149 showW = 0; % Windstress
150 N = 32;
151 c = [1 12];
152 c = [1 12]; cx = [-(1+2*c(1)) 1+2*c(2)]*1e-14; cmap = mapanom; % ecco2_bin1
153 c = [1 12]; cx = [-(1+2*c(1)) 1+2*c(2)]*1e-11; cmap = mapanom; % ecco2_bin2
154 unit = '1/m/s';
155 titf = 'PV.L115';
156 colorCOAST = [1 1 1]*.5;
157
158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3
159 case 3
160 C = JFz;
161 Clon = JFzlon; Clat = JFzlat;
162 %tit = strcat(snapshot,'/ Mechanical PV flux J^F_z and windstress');
163 tit = strvcat(['Mechanical PV flux J^F_z (positive upward), Ekman layer depth (green contours' ...
164 ' m)'],'and windstress (black arrows)');
165 %Tcolor = [0 0 0];
166 showW = 1; % Windstress
167 colorW = 'k';
168 showE = 1;
169 Econtour = [10 20:20:200];
170 N = 256;
171 c = [1 1];
172 cx = [-(1+2*c(1)) 1+2*c(2)]*1e-11; cmap = mapanom; %cmap = jet;
173 %cx = [-1 1]*10^(-11);
174 unit = 'kg/m^3/s^2';
175 titf = 'JFz';
176 showCLIM = 1;
177
178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4
179
180 case 4
181 C = JBz;
182 Clon = JBzlon; Clat = JBzlat;
183 tit = strcat(snapshot,'/ Diabatic PV flux J^B_z and windstress');
184 showW = 1; % Windstress
185 cx = [-1 1]*10^(-11);
186 unit = 'kg/m^3/s^2';
187 titf = 'JBz';
188
189 case 5
190 C = Qnet;
191 Clon = Qnetlon; Clat = Qnetlat;
192 tit = ['Net surface heat flux Q_{net} (positive downward), mixed layer depth (green contours,' ...
193 ' m) and windstress (black arrows)'];
194 tit = ['Net surface heat flux Q_{net} (positive downward), and windstress (black arrows)'];
195 showH = 0;
196 Hcontour = -[100 200:200:600];
197 showT = 0;
198 showW = 1; colorW = 'k';
199 N = 256;
200 c = [1 1]; cx = [-(1+2*c(1)) 1+2*c(2)]*200; cmap = mapanom;
201 %cx = [-1 1]*500;
202 cmap = mapanom;
203 unit = 'W/m^2';
204 titf = 'Qnet';
205 showCLIM = 1;
206 colorCOAST = [0 0 0];
207
208 case 6
209 C = JFz./JBz;
210 Clon = JFzlon; Clat = JFzlat;
211 tit = strcat(snapshot,'/ Ratio: J^F_z/J^B_z');
212 cx = [-1 1]*5;
213 unit = char(1);
214 titf = 'JFz_vs_JBz';
215
216 case 7
217 C = squeeze(ST(1,:,:));
218 C(isnan(C)) = 0;
219 Clon = STlon; Clat = STlat;
220 tit = strcat('Surface Potential density \sigma_\theta ');
221 showT = 0; %
222 cmap = flipud(hot);
223 CONT = 1; CONTv = [20:.2:30]; CONTshlab = 0;
224 CONTc = 1; CONTcv = [20:1:30]; CONTcshlab = 1;
225 cx = [23 28];
226 unit = 'kg/m^3';
227 titf = 'SIGMATHETA';
228 colorCOAST = [1 1 1]*.5;
229
230 case 8
231 C = squeeze(ZETA(1,:,:));
232 Clon = ZETAlon; Clat = ZETAlat;
233 tit = strcat('Surface relative vorticity');
234 showW = 0; % Windstress
235 showH = 0;
236 showCLIM = 1;
237 N = 256;
238 c = [0 0];
239 cx = [-(1+2*c(1)) 1+2*c(2)]*6e-5; cmap = mapanom;
240 unit = '1/s';
241 titf = 'ZETA';
242
243 case 9
244 C = abs(squeeze(ZETA(1,:,:))./squeeze(f(1,:,:)));
245 Clon = ZETAlon; Clat = ZETAlat;
246 tit = strcat('Absolute ratio between relative and planetary vorticity');
247 showW = 0; % Windstress
248 showH = 0;
249 showCLIM = 1;
250 N = 256;
251 c = [0 1];
252 cx = [0 1+3*c(2)];
253 cmap = flipud(hot); cmap = mapanom;
254 unit = '';
255 titf = 'ZETA_vs_f';
256
257 case 10
258 C = squeeze(T(1,:,:));
259 Clon = Tlon; Clat = Tlat;
260 tit = strcat('Surface Potential temperature \theta ');
261 showT = 0; %
262 N = 256; c = [0 0];
263 cmap = flipud(hot); cmap = jet;
264 CONT = 0; CONTv = [0:1:40]; CONTshlab = 0;
265 CONTc = 1; CONTcv = [0:1:40]; CONTcshlab = 1;
266 cx = [0 30];
267 unit = '^oK';
268 titf = 'THETA';
269 colorCOAST = [1 1 1]*.5;
270
271
272 end %switch what to plot
273
274 % Draw variable:
275 sp=subplot(iw,jw,1);hold on
276 if CONT ~= 1
277 m_pcolor(Clon,Clat,C);
278 shading(SHADE);
279 if uselogmap
280 colormap(logcolormap(N,c(1),c(2),cmap));
281 else
282 colormap(cmap);
283 end
284
285 if wvar(ip) == 10
286 if CONTc
287 clear cs h csh
288 EDW = [21:23];
289 [cs,h] = m_contour(Clon,Clat,C,CONTcv,'k');
290 csh = clabel(cs,h,'fontsize',8,'labelspacing',800);
291 set(csh,'visible','on');
292 for ih = 1 : length(h)
293 if find(EDW == get(h(ih),'Userdata'))
294 set(h(ih),'linewidth',1.5)
295 end %if
296 if find(EDW(2) == get(h(ih),'Userdata'))
297 set(h(ih),'linestyle','--')
298 end %if
299 end %for
300 for icsh = 1 : length(csh)
301 if find(EDW == get(csh(icsh),'userdata') )
302 set(csh(icsh),'visible','on');
303 end %if
304 end %for
305 end %if CONTc
306 end % if ST
307
308 else
309 clear cs h csh
310 [cs,h] = m_contourf(Clon,Clat,C,CONTv);
311 if uselogmap
312 colormap(mycolormap(logcolormap(N,c(1),c(2),cmap),length(CONTv)));
313 else
314 colormap(mycolormap(cmap,length(CONTv)));
315 end
316 csh = clabel(cs,h,'fontsize',8,'labelspacing',800);
317 set(csh,'visible','off');
318 if CONTshlab
319 set(csh,'visible','on');
320 end
321 if CONTc
322 for ih = 1 : length(h)
323 if find(CONTcv == get(h(ih),'CData'))
324 set(h(ih),'linewidth',1.5)
325 end
326 end
327 if CONTcshlab
328 for ih = 1 : length(csh)
329 if find(CONTcv == str2num( get(csh(ih),'string') ) )
330 set(csh(ih),'visible','on','color','k','fontsize',8);
331 set(csh(ih),'fontweight','bold','margin',1e-3);
332 end
333 end
334 end
335 end
336 end
337 caxis(cx);
338 ccol(ip) = colorbar(CBAR,'fontsize',10);
339 ctitle(ccol(ip),unit);
340 title(tit);
341 m_coast('patch',colorCOAST);
342 m_grid('xtick',360-[20:5:80],'ytick',[20:2:50]);
343 set(gcf,'name',titf);
344
345 if wvar == 5 % Qnet (Positions depend on map limits !)
346 yy=get(ccol,'ylabel');
347 set(yy,'string','Cooling Warming');
348 set(yy,'fontweight','bold');
349 end %if
350
351 if showT
352 clear cs h
353 [cs,h] = m_contour(Tlon,Tlat,squeeze(T(1,:,:)),Tcontour);
354 clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200);
355 for ih=1:length(h)
356 set(h(ih),'edgecolor',Tcolor,'linewidth',1.2);
357 end
358 end %if show THETA contours
359
360 if showQnet
361 clear cs h
362 CQnet = Qnet;
363 if wvar(ip) == 12
364 %CQnet(isnan(C)) = NaN;
365 end
366 [cs,h] = m_contour(Qnetlon,Qnetlat,CQnet,Qnetcontour);
367 Qnetmap = mycolormap(mapanom,length(Qnetcontour));
368 if ~isempty(cs)
369 clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',200);
370 for ih=1:length(h)
371 val = get(h(ih),'userdata');
372 set(h(ih),'edgecolor',Qnetmap( find(Qnetcontour == val) ,:),'linewidth',1);
373 end
374 end
375 end %if show Qnet contours
376
377 if showW
378 dx = 10*diff(Txlon(1:2)); dy = 8*diff(Txlat(1:2));
379 dx = 20*diff(Txlon(1:2)); dy = 10*diff(Txlat(1:2));
380 lo = [Txlon(1):dx:Txlon(length(Txlon))];
381 la = [Txlat(1):dy:Txlat(length(Txlat))];
382 [lo la] = meshgrid(lo,la);
383 Txn = interp2(Txlat,Txlon,Tx',la,lo);
384 Tyn = interp2(Txlat,Txlon,Ty',la,lo);
385 s = 2;
386 m_quiver(lo,la,Txn,Tyn,s,colorW,'linewidth',1.25);
387 % m_quiver(lo,la,-(1+sin(la*pi/180)).*Txn,(1+sin(la*pi/180)).*Tyn,s,'w');
388 m_quiver(360-82,37,1,0,s,'w','linewidth',1.25);
389 m_text(360-82,38,'1 N/m^2','color','w');
390 end %if show windstress
391
392 if showH
393 clear cs h
394 %[cs,h] = m_contour(MLDlon,MLDlat,MLD,Hcontour);
395 cm = flipud(mycolormap(jet,length(Hcontour)));
396 cm = mycolormap([linspace(0,0,20); linspace(1,.5,20) ;linspace(0,0,20)]',length(Hcontour));
397 cm = [0 1 0 ; 0 .6 0 ; 0 .2 0 ; 0 0 0];
398 for ii = 1 : length(Hcontour)
399 [cs,h] = m_contour(MLDlon,MLDlat,MLD,[1 1]*Hcontour(ii));
400 if ~isempty(cs)
401 % clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',300);
402 clabel(cs,h,'fontsize',8,'color',cm(ii,:),'labelspacing',600,'fontweight','bold');
403 for ih=1:length(h)
404 % set(h(ih),'edgecolor',Hcolor,'linewidth',1);
405 set(h(ih),'edgecolor',cm(ii,:),'linewidth',1.2);
406 end
407 end
408 end
409 end %if show Mixed Layer depth
410
411 if showE
412 clear cs h
413 %[cs,h] = m_contour(EKLlon,EKLlat,EKL,Econtour);
414 %cm = flipud(mycolormap(jet,length(Econtour)));
415 n = length(Econtour);
416 cm = flipud([linspace(0,0,n); linspace(1,0,n) ;linspace(0,0,n)]');
417 %cm = [0 1 0 ; 0 .6 0 ; 0 .2 0 ; 0 0 0];
418 for ii = 1 : length(Econtour)
419 [cs,h] = m_contour(EKLlon,EKLlat,EKL,[1 1]*Econtour(ii));
420 if ~isempty(cs)
421 % clabel(cs,h,'fontsize',8,'color',[0 0 0],'labelspacing',300);
422 cl=clabel(cs,h,'fontsize',8,'color',cm(ii,:),'labelspacing',600,'fontweight','bold');
423 for ih = 1 : length(h)
424 % set(h(ih),'edgecolor',Ecolor,'linewidth',1);
425 set(h(ih),'edgecolor',cm(ii,:),'linewidth',1.2);
426 end
427 end
428 end
429 end %if show Ekman Layer depth
430
431 if showCLIM
432 m_line(360-[71 62 62 71 71],[36 36 40.5 40.5 36],'color','r','linewidth',1.5)
433 end
434
435
436 if 1 % Show the date in big in the upper left corner
437 spp=subplot('position',[0 .95 .25 .05]);
438 p=patch([0 1 1 0],[0 0 1 1],'w');
439 set(spp,'ytick',[],'xtick',[]);
440 set(spp,'box','off');
441 dat = num2str(TIME(it,:));
442 dat = strcat(dat(1:4),'/',dat(5:6),'/',dat(7:8),':',dat(9:10),'H',dat(11:12));
443 text(0.1,.5,dat,'fontsize',16,...
444 'fontweight','bold','color','r','verticalalign','middle');
445 end
446
447 %%%%%%%%%%%%%%%%
448 drawnow
449 set(gcf,'position',[4 48 888 430]);
450 %videotimeline(num2str(zeros(size(TIME,1),1)),it,'b')
451 set(gcf,'color','white')
452 set(findobj('tag','m_grid_color'),'facecolor','none')
453 if prtimg
454 set(gcf,'paperposition',[0.6 6.5 25 14]);
455 %print(gcf,'-djpeg100',strcat(outimg,sla,titf,'.',snapshot,'.jpg'));
456 exportj(gcf,0,strcat(outimg,sla,titf,'.',snapshot));
457 end %if
458
459
460 end %for ip
461

  ViewVC Help
Powered by ViewVC 1.1.22