Context Navigation

plt_map.pro @ 141

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

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

Download in other formats:

Original Format