Context Navigation

plt_map.pro @ 68

Last change on this file since 68 was 41, checked in by ericg, 16 years ago
Misc modifs, incl. read_from_grid in yfx, debug_w prints, removal of mesh_lmdz and small bug corrections
File size: 38.0 KB

Line
1	PRO plt_map, cmd, iplot, win, iover, landscape
2
3	@common
4	@com_eg
5	;;
6	;;-------------------------------------
7	;;
8	;; make window plot
9	;;
10	;; EG 19-2-99
11	;;-------------------------------------
12	;
13	;
14	; cmd = command line of window
15	; win = position of window (petit[,,,])
16	; iover = overlay index
17	IF debug_w THEN print, ' '
18	IF debug_w THEN print, ' ENTER plt_map...'
19	;
20	; empty window case for on=2
21	;
22	ioverchk = iover
23
24	cmd_wrk = cmd
25
26	IF cmd.on EQ 2 THEN return
27	;
28	; incompatible options of plt_def and post_it
29	;
30	trend_typp = cmd.trend
31	; CASE cmd.timave OF
32	; '1mm': cmd.trend = '0'
33	; '3mm': cmd.trend = '0'
34	; ELSE:
35	; ENDCASE
36	;
37	; =============
38	; 1.read data
39	; =============
40	;
41	; type of plot
42
43	; time series of rms deviation on a sigma surface?
44	st_rms = 0
45	IF strpos(cmd.plt, 'st') NE -1 AND strpos(cmd.var, '@s') GT 4 THEN BEGIN
46	st_rms = 1
47	sig_surface = strmid(cmd.var, strpos(cmd.var, '@s')+2, strlen(cmd.var))
48	cmd.var = strmid(cmd.var, 0, strpos(cmd.var, '@s'))
49	legendsurf = ' on sigma='+sig_surface
50	ENDIF
51
52	CASE cmd.plt OF
53	'x': cmd.plt = cmd.plt+'_global'
54	'y': cmd.plt = cmd.plt+'_global'
55	'z': cmd.plt = cmd.plt+'_global'
56	's': cmd.plt = cmd.plt+'_global'
57	'yz': cmd.plt = cmd.plt+'_global'
58	ELSE:
59	ENDCASE
60
61	; Definition of plttyp, hotyp, dimplot
62	datyp = def_dptyp(cmd)
63
64	plttyp = datyp.plttyp
65	hotyp = datyp.hotyp
66	dimplot = datyp.dimplot
67	pltztyp = datyp.pltztyp
68	plt1dtyp = datyp.plt1dtyp
69	splot = datyp.splot
70
71	hotypchk = hotyp
72
73	; stat computations:
74
75	CASE cmd.plt OF
76	'xyt': BEGIN & plttyp = 'plt' & hotyp = 'xyt' & dimplot = 2 & END
77	ELSE:
78	ENDCASE
79
80	; xy plot on density surface (cmd.var = <var>@s<sig_surf>
81	legsurf = ''
82	sig_surf = 0.
83	IF strpos(cmd.var, '@s') GT 4 THEN BEGIN
84	splot = 1
85	sig_surf = strmid(cmd.var, strpos(cmd.var, '@s')+2, strlen(cmd.var))
86	cmd.var = strmid(cmd.var, 0, strpos(cmd.var, '@s'))
87	legsurf = ' on sigma='+sig_surf
88	ENDIF
89
90	; contour of coasts
91
92	CASE splot OF
93	1: c_cont = 0
94	0: c_cont = (!d.n_colors-1) < 255
95	ENDCASE
96
97	; scatter plot y=f(x)
98
99	IF strpos(cmd.var, '=f(') GE 1 THEN BEGIN
100	CASE cmd.plt OF
101	'xy': cmd.plt = cmd.plt+'_global'
102	ELSE:
103	ENDCASE
104	plttyp = 'yfx'
105	dimplot = 1
106	ENDIF
107
108	; vector plot
109
110	IF strpos(cmd.var, '[') EQ 0 THEN BEGIN
111	vecplot = 1
112	; sampling of vector
113	IF vector_sample GE 2 THEN BEGIN
114	vect_sample = ',UNVECTSUR=strtrim(string(vector_sample),2)'
115	ENDIF ELSE vect_sample = ''
116	ENDIF ELSE BEGIN
117	vecplot = 0
118	ENDELSE
119
120	; spectrum
121
122	spectplt = 0
123	IF strpos(cmd.plt, '@s') GE 1 THEN BEGIN
124	IF hotyp EQ 't' THEN spectplt = 1 ELSE BEGIN
125	print, ' *** spectrum [@s] applies to 1D time serie only (t_<box>)'
126	stop
127	ENDELSE
128	cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@s'))
129	ENDIF
130
131	; wavelet
132
133	waveplt = 0
134	IF strpos(cmd.plt, '@w') GE 1 THEN BEGIN
135	IF hotyp EQ 't' THEN waveplt = 1 ELSE BEGIN
136	print, ' *** wavelet [@w] applies to 1D time serie only (t_<box>)'
137	stop
138	ENDELSE
139	cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@w'))
140	ENDIF
141
142	; running std dev
143
144	run_stddev = 0
145	IF strpos(cmd.plt, '@r') GE 1 THEN BEGIN
146	IF hotyp NE 't' THEN BEGIN
147	print, ' *** running [@r<win>] std dev applies to 1D time serie only (t_<box> & 1m@t412)'
148	stop
149	ENDIF
150	run_stddev = strmid(cmd.plt, strpos(cmd.plt, '@r')+2, strlen(cmd.plt))
151	cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@r'))
152	ENDIF
153
154
155	; define grids
156	; read if new type + triangule
157	;
158	cmd1_back = cmd
159	sw_diffg = 0
160	def_grid, cmd
161	IF debug_w THEN print, ' after def_grid in plt_map: cmd.grid/varexp/vargrid ',cmd.grid,' ', varexp, ' ', vargrid
162
163	; vertical average if 3D field
164	vert_switch = 0
165	IF plttyp EQ 'plt' AND vert_type ne '0' THEN vert_switch = 1
166	IF plttyp EQ 'pltt' AND vert_type ne '0' THEN vert_switch = 1
167	IF plttyp EQ 'plt1d' AND vert_type ne '0' THEN vert_switch = 1
168	IF plt1dtyp EQ 'z' THEN vert_switch = 0
169
170	; read all data if netcdf output
171
172	IF debug_w THEN print, ' force_all_data_read in pltmap: ', force_all_data_read
173
174	IF cmd.out NE 'cdf' THEN BEGIN
175	all_data_str = ''
176	ENDIF ELSE BEGIN
177	all_data_str = ',/ALL_DATA'
178	END
179
180	; ===========
181	; Read data
182	; ===========
183
184	really_1m_st = 1
185	IF cmd.timave EQ '1y' AND strmid(cmd.plt, 0, 2) EQ 'st' THEN BEGIN
186	@densit_pltmap_read
187	ENDIF ELSE BEGIN
188	pltcmd = 'field = data_read(cmd,'''+hotyp+''','''+plttyp+''','+string(dimplot)+','+string(iover)+all_data_str+', ZMTYP = '''+cmd.plt+''')'
189	printf, nulhis, strcompress(pltcmd)
190	printf, nulhis, ' '
191	IF execute(pltcmd) EQ 0 THEN stop
192	ENDELSE
193
194	; data specific treatments
195	IF cmd.var EQ 'sla' THEN field.origin = 'diff'
196	IF cmd.var EQ 'sosstsst' AND cmd.exp EQ 'HadISST1' THEN BEGIN
197	; limit data to -1.8 C
198	idx = where(field.data EQ valmask)
199	field.data = field.data > (-1.8)
200	IF idx[0] NE -1 THEN field.data(idx) = valmask
201	ENDIF
202
203	IF n_elements(field.data) EQ 1 THEN return
204
205	; change cmd.var if macro
206	CASE STRMID(cmd.var, 0, 2) OF
207	'@@': BEGIN
208	cmd.var = field.name
209	cmd.var = STRMID(cmd.var, 2, 15)
210	IF debug_w THEN print, 'cmd.var : ', cmd.var
211	END
212	ELSE:
213	ENDCASE
214	IF vecplot EQ 1 THEN BEGIN ; vectors case
215	CASE STRMID(cmd.var2, 0, 2) OF
216	'@@': begin
217	cmd.var2 = STRMID(cmd.var2, 2, 15)
218	END
219	ELSE:
220	ENDCASE
221	ENDIF
222
223	; spectum plot
224	IF spectplt EQ 1 THEN plttyp = 'spec'
225	; wavelet plot
226	IF waveplt EQ 1 THEN plttyp = 'wavelet'
227
228	xlogax = ''
229
230	; density pojection plot (splot=1)
231	IF splot EQ 1 THEN BEGIN
232	; keep grid attributes
233	izminmesh_b = izminmesh
234	izmaxmesh_b = izmaxmesh
235	jpk_b = jpk
236	jpkglo_b = jpkglo
237	gdept_b = gdept
238	gdepw_b = gdepw
239	e3t_b = e3t
240	e3w_b = e3w
241	tmask_b = tmask
242	; IF cmd.var EQ 'vosigvol' THEN xlogax = ',xlog = 1'
243	; modify vertical grid -> sigma
244	n_sig = (sig_max - sig_min)/sig_del + 1
245	z_sig = sig_min+findgen(n_sig)*sig_del
246	izminmesh = 0
247	izmaxmesh = n_sig-1
248	jpk = long(izmaxmesh-izminmesh+1)
249	jpkglo = jpk
250	gdept = z_sig
251	gdepw = gdept
252	e3t = shift(gdept, 1)-gdept
253	e3t[0] = e3t[1]
254	e3w = e3t
255	prof1 = sig_min
256	prof2 = sig_max
257	; grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz ; added PDW 12/5/04
258	; tmask = mask ; added PDW 12/5/04
259	tmask = reform(reform(tmask[, , 0], jpi*jpj)#replicate(1, jpk), jpi, jpj, jpk)
260	domdef
261	END
262
263	IF debug_w THEN print, ' cmd after data_read in plt_map = ', cmd
264	;
265	; ======================
266	; 2. window attributes
267	; ======================
268	;
269
270	; find field attributes
271	; ---------------------
272
273	fldatt = {name:'', assos:'', min: 0.0, max: 0.0, int: 0.0, mult: 0.0, add: 0.0, unit: '', mid: 0.0, homin:0.0, homax:0.0, min1d:0.0, max1d:0.0, dmax:0.0}
274
275	; extremum
276	fldext = fld_pltext(cmd.var, plttyp, dimplot, hotyp)
277
278	fldatt.name = fldext.name
279	fldatt.min = fldext.min
280	fldatt.max = fldext.max
281	fldatt.homin = fldext.homin
282	fldatt.homax = fldext.homax
283	fldatt.min1d = fldext.min1d
284	fldatt.max1d = fldext.max1d
285	fldatt.dmax = fldext.dmax
286	fldatt.assos = fldext.assos
287
288	; interval + change of unit
289	fldint = fld_pltint(cmd.var, plttyp, dimplot, hotyp)
290
291	; help, fldint, /struct
292
293	fldatt.int = fldint.int
294	fldatt.mult = fldint.mult
295	fldatt.add = fldint.add
296	fldatt.unit = fldint.unit
297	fldatt.mid = fldint.mid
298
299	; print, 'plt_map 1 :', fldatt.mult, fldatt.add, fldatt.min, fldatt.max
300
301	IF fldatt.mult NE 1. OR fldatt.add NE 0. THEN field.units = fldatt.unit
302	fld = field.data*fldatt.mult + fldatt.add
303	IF (where(field.data EQ valmask))[0] NE -1 THEN $
304	fld[where(field.data EQ valmask)] = valmask
305
306	; field 2 if needed
307
308	IF plttyp EQ 'yfx' OR vecplot EQ 1 THEN BEGIN
309	fldatt2 = fldatt
310
311	fldext2 = fld_pltext(cmd.var2, plttyp, dimplot, hotyp)
312
313	fldatt2.name = fldext2.name
314	fldatt2.min = fldext2.min
315	fldatt2.max = fldext2.max
316	fldatt2.homin = fldext2.homin
317	fldatt2.homax = fldext2.homax
318	fldatt2.min1d = fldext2.min1d
319	fldatt2.max1d = fldext2.max1d
320	fldatt2.dmax = fldext2.dmax
321
322	; interval + change of unit
323	fldint2 = fld_pltint(cmd.var2, plttyp, dimplot, hotyp)
324
325	fldatt2.int = fldint2.int
326	fldatt2.mult = fldint2.mult
327	fldatt2.add = fldint2.add
328	fldatt2.unit = fldint2.unit
329	fldatt2.mid = fldint2.mid
330
331	IF fldatt2.mult NE 1. OR fldatt2.add NE 0. THEN field.units2 = fldatt2.unit
332	fld2 = field.data2*fldatt2.mult + fldatt2.add
333	IF (where(field.data2 EQ valmask))[0] NE -1 THEN $
334	fld2[where(field.data2 EQ valmask)] = valmask
335
336	ENDIF
337
338	; If running std dev change min/max
339	IF run_stddev GT 0 THEN BEGIN
340	fldatt.homin = 0.
341	fldatt.homax = fldext.dmax
342	ENDIF
343
344	; legend text
345	; -----------
346
347	; IF plttyp EQ "pltt" THEN BEGIN
348	CASE strmid(cmd.trend, 0, 1) OF
349	'1': trd_txt = ' trend'
350	'2': trd_txt = ' drift'
351	'3': trd_txt = ' inverse trend'
352	'4': trd_txt = ' anomaly'
353	'5': trd_txt = ' filter'
354	'7': trd_txt = ' month sel'
355	ELSE: trd_txt = ''
356	ENDCASE
357	; ENDIF ELSE trd_txt = ''
358	IF field_int EQ 1 AND dimplot EQ 1 THEN BEGIN
359	trd_txt = trd_txt+' integral'
360	ENDIF
361
362	varlegend = field.legend+trd_txt+' ('+field.units+')'
363
364	CASE plttyp OF
365	'yfx': varlegend2 = field.legend2+trd_txt+' ('+field.units2+')'
366	ELSE: BEGIN
367	IF run_stddev GT 0 THEN trd_txt = trd_txt+' running Std Dev ['+string(run_stddev, format = '(I3)')+']'
368	END
369	ENDCASE
370
371	; date text
372
373	CASE strmid(cmd.timave, 0, 2) OF
374	'1m': BEGIN
375	mn = def_month(cmd.timave, cmd.date1)
376	IF strmid(cmd.timave, 0, 3) EQ '1mm' THEN $
377	date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $
378	ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2)
379	END
380	'3m': BEGIN
381	mn = def_month(cmd.timave, cmd.date1)
382	IF strmid(cmd.timave, 0, 3) EQ '3mm' THEN $
383	date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $
384	ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2)
385	END
386	ELSE:date_txt = cmd.timave+' '+cmd.date1
387	ENDCASE
388
389	CASE plttyp OF
390	'pltt':BEGIN
391	date_txt = cmd.timave
392	IF strmid(cmd.timave, 1, 2) EQ 'mm' THEN $
393	date_txt = cmd.timave+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')'
394	END
395	'yfx': BEGIN
396	IF hotyp NE '-' THEN date_txt = cmd.timave+' '+cmd.date1+' '+cmd.spec
397	END
398	'wavelet': BEGIN
399	date_txt = cmd.timave+' Wavelet'
400	IF strmid(cmd.timave, 1, 2) EQ 'mm' THEN $
401	date_txt = cmd.timave+' Wavelet ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')'
402	END
403	ELSE:
404	ENDCASE
405
406	; min/max/int
407
408	; defined if fraction of x/y.range is added to plot domain
409
410	; IF free_1d_minmax EQ 'no' THEN fraction = 0. ELSE fraction = 1.0
411	fraction = 0.
412
413	CASE dimplot OF
414	1: BEGIN
415	CASE plttyp OF
416	'pltt': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END
417	'wavelet': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END
418	'plt1d': BEGIN & minc = fldatt.min1d & maxc = fldatt.max1d & END
419	'yfx': BEGIN ; y=f(x) scatter plot
420	IF long(strmid(cmd.trend, 0, 1)) GT 0 THEN BEGIN
421	IF run_stddev EQ 0 THEN BEGIN
422	minc = -fldatt.dmax & maxc = fldatt.dmax
423	minc2 = -fldatt2.dmax & maxc2 = fldatt2.dmax
424	ENDIF ELSE BEGIN
425	minc = 0. & maxc = fldatt.dmax
426	minc2 = 0. & maxc2 = fldatt2.dmax
427	ENDELSE
428	ENDIF ELSE BEGIN
429	minc = fldatt.min1d & maxc = fldatt.max1d
430	minc2 = fldatt2.min1d & maxc2 = fldatt2.max1d
431	ENDELSE
432	END
433	ELSE: BEGIN & minc = 0. & maxc = 0. & END
434	ENDCASE
435	fldatt.int = !VALUES.F_NAN
436	END
437	2: BEGIN
438	CASE plttyp OF
439	'pltt': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END
440	ELSE: BEGIN & minc = fldatt.min & maxc = fldatt.max & END
441	ENDCASE
442	END
443	ENDCASE
444
445	IF cmd.var NE fld_prev OR cmd.trend NE '0' THEN BEGIN
446	; field.origin EQ 'diff' OR
447	print, ''
448	CASE dimplot OF
449	2: BEGIN
450	print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max/int=', $
451	minc, maxc, fldatt.int
452	IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $
453	print, ' --> Modify field using : ', 'mult/add=' , $
454	fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit
455	END
456	1: BEGIN
457	print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max=',minc, maxc
458	IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $
459	print, ' --> Modify field using : ', 'mult/add=' , $
460	fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit
461	END
462	0: BEGIN
463	print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max=',minc, maxc
464	IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $
465	print, ' --> Modify field using : ', 'mult/add=' , $
466	fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit
467	print, ' '+cmd.var2, ' plot attributes ('+plttyp+') : min/max=',minc2, maxc2
468	IF ( fldatt2.mult NE 1.0 OR fldatt2.add NE 0.0 ) THEN $
469	print, ' --> Modify field using : ', 'mult/add=' , $
470	fldatt2.mult, fldatt2.add, ' to obtain : ', fldatt2.unit
471	END
472	ELSE:
473	ENDCASE
474	print, ''
475
476	ENDIF
477	IF ( fldatt.mult EQ -1.0) THEN BEGIN
478	temp_m = minc
479	minc = -maxc
480	maxc = -temp_m
481	ENDIF
482	IF finite(fldatt.int) EQ 0 THEN BEGIN
483	intcmd = ''
484	colbarfor = ''
485	fmt = '(f6.1)'
486	ENDIF ELSE BEGIN
487	intcmd = ',int = '+string(fldatt.int)
488	IF fldatt.int LT 0.1 THEN BEGIN
489	fmt = '(f6.2)'
490	ENDIF ELSE IF fldatt.int LT 1 THEN BEGIN
491	fmt = '(f5.1)'
492	ENDIF ELSE BEGIN
493	fmt = '(f5.0)'
494	ENDELSE
495	IF long(fldatt.int) NE 0 THEN BEGIN
496	IF fldatt.int/long(fldatt.int) NE 1 THEN fmt = '(f5.1)'
497	ENDIF
498	colbarfor = ', format = '''+fmt+''''
499	ENDELSE
500
501	; if window > 1 or overlay > 1 noerase
502
503	IF win[2] NE 1 OR iover GT 1 THEN BEGIN
504	noerase = 1
505	ENDIF ELSE BEGIN
506	noerase = 0
507	ENDELSE
508
509	; choose window number
510
511	CASE multi_win OF
512	1: window_number = ', window='+string(iplot)
513	ELSE: window_number = ''
514	ENDCASE
515
516	; color label control
517
518	labstr = ''
519	IF iover EQ 1 THEN BEGIN
520	nofill = 1-shading
521	ENDIF ELSE BEGIN
522	nofill = 1
523	ENDELSE
524
525	nocoltxt = ''
526	IF nofill EQ 1 THEN nocoltxt = ',NOCOULEUR=1'
527
528	IF dimplot NE 2 THEN nofill = 1
529
530	nocolorbar = 0
531	IF col_palette EQ 'no' THEN nocolorbar = 1
532	IF pal_type EQ '2dom' THEN nocolorbar = 1
533	IF iover GT 1 THEN nocolorbar = 1
534	IF dimplot NE 2 THEN nocolorbar = 1
535
536	IF cmd.out NE 'cdf' THEN BEGIN
537	readpal = 0
538	IF nofill EQ 0 AND dimplot GT 1 THEN readpal = 1
539	IF waveplt EQ 1 THEN readpal = 2
540	IF field.origin EQ 'div' AND plttyp NE 'pltt' THEN readpal = 3
541	;; Necessary for correlation plots (overlay of 2D fields)
542	IF dimplot GT 1 AND nover GT 1 AND iover GT 1 THEN readpal = 1
543	IF readpal GE 1 THEN BEGIN
544	lec_pal_gmt, cmd.var, c_anot_str, fmt, found, readpal
545	colbarfor = ', format = '''+fmt+''''
546	IF found EQ 1 THEN BEGIN
547	labstr = ',label=3'
548	; if mincmd/naxcmd not defined then use one from palette
549	; if define take min of 2 mins and max of 2 maxs
550	min_palette = min(levels_gmt)
551	max_palette = max(levels_gmt)
552	IF finite(minc) NE 0 THEN minc = min([minc, min_palette]) ELSE minc = min_palette
553	IF finite(maxc) NE 0 THEN maxc = max([maxc, max_palette]) ELSE maxc = max_palette
554	print, ' '+cmd.var, ' plot attributes ('+plttyp+') modified via color palette: min/max=', $
555	minc, maxc
556	ENDIF
557	ENDIF
558	ENDIF
559	colbar = colbarfor
560
561	mincmd = ''
562	maxcmd = ''
563
564	IF finite(minc) NE 0 THEN BEGIN
565	mincmd = ','+string(minc)
566	maxcmd = ','+string(maxc)
567	ENDIF
568
569	; contour annotation + other color bar controls
570
571	IF n_elements(c_anot_str) EQ 0 THEN BEGIN
572	c_anot = ''
573	ENDIF ELSE BEGIN
574	c_anot = ',c_annotation=c_anot_str'
575	; sampling of colorbar annotation
576	IF col_palette EQ 'yes' THEN BEGIN
577	sampling = win[0]*win[1]
578	; idx_cb = lindgen((ncont_gmt-2)/sampling)*sampling+1
579	idx_cb0 = lindgen(ncont_gmt-1)+1
580	idx_cb = lindgen(ncont_gmt-2)+2
581	idx_cb1 = lindgen(ncont_gmt-2)+1
582	colbar = colbar+',discret=coul_gmt[idx_cb1],div=n_elements(idx_cb0)-1,cb_label=levels_gmt[idx_cb0]'
583	; c_anot = ',c_annotation=levels_gmt[idx_cb0]'
584	ENDIF
585	ENDELSE
586	;
587	; vertical grid type
588	;
589	IF debug_w THEN print, ' splot=', splot
590	CASE mesh_type OF
591	'oce': BEGIN
592	IF splot EQ 1 THEN BEGIN
593	type_yz = ',type_yz =''sigma'''
594	zoom_txt = ', zoom='+string(sig_max)
595	ENDIF ELSE BEGIN
596	type_yz = ',type_yz=''m'''
597	zoom_txt = ', zoom='+string(zoom_z)
598	ENDELSE
599	END
600	'atm': BEGIN & type_yz = ',type_yz=''hPa''' & zoom_txt = ', zoom='+string(hpa_max) & END
601	ELSE:type_yz = ''
602	ENDCASE
603
604	; continents color, thickness
605	; c_cont = 100
606	; continent fill
607	; real continent drawn
608	strcont =''
609	IF mesh_type NE 'oce' THEN BEGIN
610	IF cont_fill EQ 0 THEN strcont = ',/CONT_NOFILL, nite=0'
611	IF cont_real GE 1 THEN strcont = strcont+', REALCONT = '+string(strtrim(cont_real, 2))+', COAST_THICK = 2'
612	ENDIF
613
614	; calendar type
615
616	cal_typ = ''
617	IF calendar_type GT 1 THEN cal_typ = ',ndayspm=calendar_type'
618
619	; titles options
620	; By default, SAXO puts a title and a subtitle for the first plot
621	; Then, if you overlay another field, title and subtitle are set to ''
622	CASE title_type OF
623	"T": tit_str = ',subtitle='''','+marge_option
624	"S": tit_str = ',title='''','+marge_option
625	"TS": tit_str = ','+marge_option
626	"off":tit_str = ',title='''',subtitle='''','+marge_option
627	ENDCASE
628
629	; fill_space option
630
631	filltxt = ''
632	IF fill_space EQ 1 THEN filltxt = ',/rempli'
633
634	; axis charsize option (x/ychartxt read from plt_def)
635
636	; common command string to plots
637
638	com_strplt = ',petit = ['+string(win[0])+','+string(win[1])+','+string(win[2])+']'+filltxt+nocoltxt+', LANDSCAPE = '+string(landscape)+', NOCOLORBAR = '+string(nocolorbar)+', NOERASE = '+string(noerase)+look+labstr+c_anot+colbar+cal_typ+',cont_thick=2'+type_yz+window_number+tit_str+xlogax+', xcharsize='+string(xchartxt)+', ycharsize='+string(ychartxt)
639
640	; add contour_options not overlay
641
642	IF iover EQ 1 THEN com_strplt = com_strplt+contour_options
643
644	; run specific fixes
645
646	IF (strpos(cmd.exp, 'MIHR') NE -1) AND cmd.var EQ 'tauu' AND strpos(cmd.timave, '1m') NE -1 AND strpos(cmd.timave, '1mm') EQ -1 THEN BEGIN
647	fld = -fld
648	print, ' *** WARNING: Multiply tauu by -1 ', cmd.exp, cmd.var, cmd.timave
649	ENDIF
650
651	IF debug_w THEN print, ' cmd before making output in plt_map = ', cmd
652
653
654	; ============================
655	; 3. make output (data/plot)
656	; ============================
657	;
658
659	CASE cmd.out OF
660	'data': write_data = long((iplot-1)*10+win[2]) ; write 1D data ascii in trends.pro
661	'tcdf': write_data = -long((iplot-1)*10+win[2]) ; write 1D data netcdf in trends.pro
662	'cdf': write_data = 0
663	ELSE: write_data = 0
664	ENDCASE
665
666	; make output
667
668	CASE cmd.out OF
669	'tv': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ tvnplot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
670	erase
671
672	CASE (size(fld))[0] OF
673	2: fldtv = fld
674	3: BEGIN
675	niveau = xquestion('3d data : which level ?', '1', /chkwidget)
676	fldtv = fld[, , niveau]
677	END
678	4: BEGIN
679	niveau = xquestion('3d data : which level ?', '1', /chkwidget)
680	timeplot = xquestion('4d data : which date ?', '1', /chkwidget)
681	fldtv = fld[, , niveau, timeplot]
682	END
683	ELSE:
684	ENDCASE
685
686	CASE cmd.grid OF
687	; 'T': tvnplot, fldtv*tmask
688	'T': tvplus, fldtv, min = fldatt.min, max = fldatt.max
689	'U': tvplus, fldtv, min = fldatt.min, max = fldatt.max
690	'V': tvplus, fldtv, min = fldatt.min, max = fldatt.max
691	'W': tvplus, fldtv, min = fldatt.min, max = fldatt.max
692	ELSE: tvplus, fldtv, min = fldatt.min, max = fldatt.max
693	ENDCASE
694	END
695
696	'cdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
697
698	; build structures containing grids
699	cdf_description = nc_build(cmd, field, field.direc, vargrid)
700	; build netCDF file (here if whole data)
701	IF hotyp EQ '-' OR hotyp EQ 'xyt' THEN BEGIN
702	fldcdf = {data:field.data, units:field.units, short_name:cmd.var, long_name:field.legend, missing_value:valmask, direc:field.direc}
703	file_out = cmd.var+'_'+strtrim(string(FORMAT = '(I2.2)', iplot), 2)+'.nc'
704	pltcmd ='nc_put, fldcdf, file_out, ncdf_db = hom_idl'+cdf_description
705	printf, nulhis, strcompress(pltcmd)
706	res = execute(pltcmd)
707	ENDIF
708	END
709	'tcdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ 1D netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
710	mask_z, fld, cmd, boite_plt1d, dimplot
711	fld = checkfield(fld, 'pltt', type = datyp.hotyp, boite = boite_plt1d)
712	IF debug_w THEN print, 'size(fld)', size(fld)
713	IF debug_w THEN print, 'boite_plt1d', boite_plt1d
714	IF cmd.trend GT 0 THEN BEGIN
715	fld = trends(fld, cmd.trend, datyp.hotyp)
716	add_txt = trd_txt
717	ENDIF ELSE print, 'You need to have a trend to make 1D netcdf output'
718	END
719	ELSE: BEGIN ; $$$$$$$$$$$$$$$$$$$$$$$$$$ make plot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
720
721	CASE plttyp OF
722	'plt': BEGIN
723	;
724	; Map plot
725	; --------
726	;
727	print, ' '+cmd.var+' data to plot min and max : ', $
728	min(fld(where (fld NE -valmask)), /nan), max(fld(where (fld NE valmask)), /nan)
729	print, ''
730
731	map_cmd = ''
732
733	niveau = 1
734
735	IF vecplot EQ 0 THEN BEGIN
736	; eddy field ? (if so, remove zonal mean)
737	IF strpos(cmd.plt, '@z') GT 1 THEN BEGIN
738	fldzm = moyenne(fld, 'x', boite=[20,380,-90,90], /NAN)
739	fldzm = replicate(1, (size(fld))(1))#fldzm
740	fld = fld-fldzm
741	filleg = ' Eddy '
742	ENDIF ELSE BEGIN
743	filleg = ''
744	ENDELSE
745
746	IF hotyp EQ 'xyt' THEN date_txt = 'monthly'
747
748	pltcmd = 'plt,{a:fld, d:'''+date_txt+''', n:'''+filleg+varlegend+' '+legbox+legsurf+''', u:'''+field.units+''', g:'''+vargrid+'''}'+mincmd+maxcmd+intcmd+com_strplt+strcont+map_cmd
749	;; For 2D plots with overlays (correlation and proba for instance)
750	;; Scratch the titles and define contours style for the second plot
751	IF nover GT 1 AND iover GT 1 THEN BEGIN
752	contour_style = ',c_thick=2,c_linestyle=2,c_charthick=2'
753	CASE title_type OF
754	'T': pltcmd = pltcmd+',title='''''+contour_style
755	'S': pltcmd = pltcmd+',subtitle='''''+contour_style
756	'TS': pltcmd = pltcmd+',title='''',subtitle='''''+contour_style
757	'off': pltcmd = pltcmd+contour_style
758	END
759	ENDIF
760	ENDIF ELSE BEGIN
761
762	fld = {a:fld, g:strmid(cmd.grid, 0, 1)}
763	fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)}
764	fldv = {data1: fld, data2: fld2}
765	vargrid = strmid(cmd.grid, 2, 1)
766
767	pltcmd = 'ajoutvect,fldv'+vect_sample+com_strplt+strcont+map_cmd
768
769	ENDELSE
770
771	IF debug_w THEN print, strcompress(pltcmd)
772	printf, nulhis, strcompress(pltcmd)
773	res = execute(pltcmd(0))
774
775	END
776	'pltz': BEGIN
777	;
778	; Vertical section/mean plot
779	; --------------------------
780	;
781	g = gphit
782	t = tmask
783
784	; if already 1D, reform fld
785	sz = size(fld)
786	IF sz[0] EQ 2 THEN BEGIN
787	IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2])
788	ENDIF
789	IF sz[0] EQ 3 THEN BEGIN
790	IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3])
791	IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3])
792	IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2])
793	ENDIF
794
795	; mask fld
796	mask_z, fld, cmd, boite_pltz, dimplot, legz
797	printf, nulhis, ' boite_pltz = ', boite_pltz
798
799	print, ' '+cmd.var+' data to plot min and max : ', $
800	min(fld(where (fld NE valmask)), /NAN), max(fld(where (fld NE valmask)), /NAN)
801	print, ''
802	IF vecplot EQ 0 THEN BEGIN
803
804	IF cmd.var EQ 'vozonbsf' THEN xindstr = ', /xindex' ELSE xindstr = ''
805
806	pltcmd = 'pltz,{a:fld, d:'''+date_txt+''', n:'''+varlegend+' '+legz+''', u:'''+field.units+''', g:'''+vargrid+'''}'+mincmd+maxcmd+intcmd+',/'+pltztyp+', boite=boite_pltz'+zoom_txt+com_strplt+xindstr
807
808
809	ENDIF ELSE BEGIN
810
811	fld = {a:fld, g:strmid(cmd.grid, 0, 1)}
812	fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)}
813	fldv = {data1: fld, data2: fld2}
814	vargrid = strmid(cmd.grid, 2, 1)
815
816	pltcmd = 'ajoutvectz,fldv'+com_strplt+strcont+',type='''+strmid(cmd.plt, 0, 2)+''', boite=boite_pltz'
817
818	ENDELSE
819
820	IF debug_w THEN print, pltcmd
821	printf, nulhis, strcompress(pltcmd)
822	res = execute(pltcmd(0))
823
824
825	; overlay bowl if sig_bowl=1
826	IF sig_bowl EQ 1 THEN begin
827	; define line color, thickness and type
828	iover = 2
829	overc = overlay_type(iover, dimplot)
830	; type of latitude axis
831	IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $
832	sin = ',/sin' ELSE sin = ''
833	plt1dtyp = strmid(pltztyp, 0, 1)
834	noerase = 1
835	com_strplt = ', NOERASE = '+string(noerase)
836	pltcmd = 'plt1d,field.bowl,/'+plt1dtyp+', boite=boite_pltz'+overc+sin+com_strplt
837	IF debug_w THEN print, pltcmd
838
839	res = execute(pltcmd(0))
840	ENDIF
841
842	gphit = g
843	tmask = t
844
845	END
846	'plt1d': BEGIN
847	;
848	; 1D, x, y, z plot
849	; ----------------
850	;
851	g = gphit
852	t = tmask
853	; if already 1D, reform fld
854	sz = size(fld)
855	IF sz[0] EQ 2 THEN BEGIN
856	IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2])
857	ENDIF
858	IF sz[0] EQ 3 THEN BEGIN
859	IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3])
860	IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3])
861	IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2])
862	ENDIF
863	; mask fld
864	mask_z, fld, cmd, boite_plt1d, dimplot, legz
865
866	IF debug_w THEN print, boite_plt1d
867	niveau = 1
868	; define line color, thickness and type
869	overc = overlay_type(iover, dimplot)
870	; type of latitude axis
871	IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $
872	sin = ',/sin' ELSE sin = ''
873
874	print, ' '+cmd.var+' data to plot min and max : ', $
875	min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask)))
876	print, ''
877	take_log = ',take_log=0'
878	IF cmd.var EQ 'vosigvol' AND n_elements(str_sep(cmd.exp,'-')) EQ 1 THEN BEGIN ; don't take log if difference plot
879	take_log = ',take_log=1'
880	ENDIF
881	pltcmd = 'plt1d,{a:fld, d:'''+date_txt+''', n:'''+varlegend+' '+legz+''', u:'''+field.units+''', g:'''+vargrid+'''},'''+plt1dtyp+''''+mincmd+maxcmd+intcmd+', boite=boite_plt1d'+overc+sin+com_strplt+take_log
882	printf, nulhis, 'boite_plt1d=',boite_plt1d
883	printf, nulhis, strcompress(pltcmd)
884	IF debug_w THEN print, pltcmd
885	res = execute(pltcmd(0))
886
887	@legend_overlay
888
889	gphit = g
890	tmask = t
891
892	END
893	'pltt': BEGIN
894	;
895	; hovmoeller
896	; -----------
897	;
898	g = gphit
899	t = tmask
900
901	; mask fld
902	mask_z, fld, cmd, boite_pltt, dimplot, legz
903	; define line color, thickness and type
904	overc = overlay_type(iover, dimplot)
905	; define additional text for pltt
906	@add_txt_pltt
907	; repeat for seasonal mean
908	CASE cmd.timave OF
909	'1mm': rep_txt = ',repeat_c=nb_cycles'
910	ELSE: rep_txt = ''
911	ENDCASE
912
913	print, ' '+cmd.var+' data to plot min and max : ', $
914	min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask)))
915	print, ''
916	vardate = 'toto'
917	grille, mask, glam, gphi, gdep, nx, ny,nz
918	IF st_rms EQ 1 THEN BEGIN ; time series of rms deviation on a sigma surface
919	@densit_pltmap_plot
920	ENDIF ELSE BEGIN
921	pltcmd = 'pltt,{a:fld, n:'''+date_txt+' '+varlegend+' '+legbox+''', u:'''+field.units+filleg+''', g:'''+vargrid+'''}, '''+hotyp+''''+mincmd+maxcmd+intcmd+',boite=boite_pltt'+com_strplt+overc+filtxt+',TREND_TYPE='+cmd.trend+rep_txt ; +',/integration'
922	ENDELSE
923
924	printf, nulhis, 'boite_pltt=',boite_pltt
925	printf, nulhis, strcompress(pltcmd)
926	IF debug_w THEN print, pltcmd
927	res = execute(pltcmd(0))
928	IF iover EQ 4 THEN print, 'There might be a pb with the legends !'
929	; legend if plot=t
930	IF hotyp EQ 't' THEN BEGIN
931	; positions of the legend depend on nb_cycles
932	; it is different from 1 only with 1mm case
933	IF cmd.timave NE '1mm' THEN nb_cycles = 1
934	@legend_overlay
935	ENDIF
936	gphit = g
937	tmask = t
938	END
939
940	'yfx': BEGIN
941	@yfx
942	END
943	'spec': BEGIN
944	;
945	; Spectrum plot y=fft(x)
946	; ---------------------
947	;
948	t = tmask
949	; plot domain
950	mask_z, fld, cmd, boite_pltspec, dimplot, legspec
951	; define line color and type
952	overc = overlay_type(iover, dimplot)
953	vardate = date_txt
954	varname = varlegend
955	varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec
956
957	; make mean in box
958	fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex)
959	; apply trends
960	IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't')
961
962	; make spectrum
963	time_interval = time[1]-time[0]
964
965	print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360
966	print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360
967	print, ' spectrum window (years)/chunks = ', spec_win/360, long(time(jpt-1)-time(0)+time_interval)/spec_win
968	print, ' '
969
970	; number of time windowsš
971	nt_win = long((time(jpt-1)-time(0)+time_interval)/spec_win)
972	idx_win_size = long(spec_win/time_interval)
973	IF idx_win_size*nt_win NE jpt THEN BEGIN
974	print, ' *** Warning : spec_win must divide lenght of time serie ', idx_win_size, jpt
975	return
976	ENDIF
977	; mean of idx_win_size chunks
978	spect = spectrum(fld[0:idx_win_size-1], time_interval)
979	FOR it = 2, nt_win DO BEGIN
980	idx1 = idx_win_size*(it-1)
981	idx2 = idx1 + idx_win_size - 1
982	specc = spectrum(fld[idx1:idx2], time_interval)
983	spect = spect + specc
984	ENDFOR
985	spect = spect/nt_win
986	; build new time array
987	idx = where (spect[0, *] LE max_spec)
988	jpt = n_elements(idx)
989	fld = reverse(reform(spect[1, idx], jpt))
990	time = findgen(jpt)
991	time = reverse(reform(spect[0, idx], jpt))
992	IF max(time) GT 20*360 THEN BEGIN
993	time = time/360
994	xtitle = 'Years'
995	ENDIF ELSE IF max(time) GT 5*360 THEN BEGIN
996	time = time/30
997	xtitle = 'Months'
998	ENDIF ELSE xtitle = 'Days'
999	; time/space filter ?
1000	IF strpos(cmd.plt, '@f') GT 1 THEN BEGIN
1001	filter = long(strmid(cmd.plt, strpos(cmd.plt, '@f')+3, strlen(cmd.plt)-strpos(cmd.plt, '@f')-3))
1002	fildirec = strmid(cmd.plt, strpos(cmd.plt, '@f')+2, 1)
1003	IF fildirec EQ 't' THEN BEGIN
1004	xtitle = xtitle+' (filter '+strtrim(string(filter), 2)+' warning : discrete filter)'
1005	fld = smooth(fld, filter)
1006	fld[0:filter/2-1] = 0.
1007	fld[(size(fld))[1]-filter/2-1:(size(fld))[1]-1] = 0.
1008	ENDIF
1009	ENDIF
1010	; plot min/max
1011	min1 = min(fld)
1012	max1 = max(fld)
1013	min2 = 0
1014	max2 = max(time)
1015	!x.range = [min2-abs(max2-min2)/5.,max2+abs(max2-min2)/5.]
1016	!y.range = [min1-abs(max1-min1)/5.,max1+abs(max1-min1)/5.]
1017
1018	; plot
1019	ytitle = 'Power spectrum (window='+strtrim(string(spec_win/360), 2)+'y)'
1020	pltcmd = 'splot,time,fld,xstyle=1,ystyle=1,title=varunit,xtitle=xtitle,ytitle=ytitle'+overc+com_strplt
1021	printf, nulhis, 'boite_pltspec=',boite_pltspec
1022	printf, nulhis, strcompress(pltcmd)
1023	IF debug_w THEN print, pltcmd
1024	res = execute(pltcmd(0))
1025
1026	tmask = t
1027
1028	END
1029	'wavelet': BEGIN
1030	;
1031	; Wavelet plot
1032	; ------------
1033	;
1034	t = tmask
1035	; plot domain
1036	mask_z, fld, cmd, boite_pltspec, dimplot, legspec
1037	; define line color and type
1038	vardate = date_txt
1039	varname = varlegend + ' '+date_txt
1040	varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec
1041	; make mean in box
1042	fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex)
1043	; apply trends
1044	IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't')
1045
1046	; make spectrum
1047	time_interval = time[1]-time[0]
1048
1049	print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360
1050	print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360
1051	print, ' '
1052	; make wavelet
1053	wave = wavelet(fld,time_interval,period=period,coi=coi,/pad,signif=signif)
1054	power = abs(wave)^2
1055	nscale = n_elements(period)
1056	period = period/365 ; to have period in years
1057	; compute significance
1058	signif = rebin(transpose(signif),jpt,nscale)
1059	; plot wavelet
1060	mincmd = ','+string(min(power))
1061	maxcmd = ','+string(max(power))
1062	boite_pltspec = [0, 0, min(period), max(period)]
1063	; domdef, boite_pltspec
1064	lat1r = lat1
1065	lat2r = lat2
1066	lat1 = 0
1067	lat2 = max_spec/360
1068	key_onearth = 0
1069	pltcmd = 'plttg,power,period'+mincmd+maxcmd+intcmd+',boite=boite_pltspec,reverse_y=1,nocontour=1'+com_strplt
1070	printf, nulhis, 'boite_pltspec=',boite_pltspec
1071	printf, nulhis, strcompress(pltcmd)
1072	IF debug_w THEN print, pltcmd
1073	res = execute(pltcmd(0))
1074	contour,abs(wave)^2/signif,time,period, /overplot,level=1.0,c_annot='95%'
1075	plot,time,coi/365, noclip = 0, /noerase
1076	tmask = t
1077	lat1 = lat1r
1078	lat2 = lat2r
1079	key_onearth = 1
1080	END
1081	ELSE: BEGIN
1082	print, ' unknown projection plot ', cmd.plt, ' / plot type = ', plttyp
1083
1084	END
1085	ENDCASE
1086
1087	END
1088
1089	ENDCASE
1090
1091	fld_prev = cmd.var
1092	;
1093	; reset incompatible options of plt_def
1094
1095	cmd.trend = trend_typp
1096
1097	;
1098	; reset vertical grid after density bining
1099
1100	IF splot EQ 1 THEN BEGIN
1101
1102	izminmesh = izminmesh_b
1103	izmaxmesh = izmaxmesh_b
1104	jpk = jpk_b
1105	jpkglo = jpkglo_b
1106	gdept = gdept_b
1107	gdepw = gdepw_b
1108	e3t = e3t_b
1109	e3w = e3w_b
1110	tmask = tmask_b
1111
1112	ENDIF
1113	IF debug_w THEN print, ' ...EXIT plt_map'
1114
1115	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format