Context Navigation

plt_map.pro @ 26

Last change on this file since 26 was 26, checked in by kolasinski, 16 years ago
Add varexp in data_read, plt_map, mesh_from_file, mesh_orca and mesh_pcmdi
File size: 51.6 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	sw_diffg = 0
159	def_grid, cmd
160	IF debug_w THEN print, ' after def_grid in plt_map: cmd.grid, varexp, vargrid ',cmd.grid, varexp, vargrid
161
162	; vertical average if 3D field
163	vert_switch = 0
164	IF plttyp EQ 'plt' AND vert_type ne '0' THEN vert_switch = 1
165	IF plttyp EQ 'pltt' AND vert_type ne '0' THEN vert_switch = 1
166	IF plttyp EQ 'plt1d' AND vert_type ne '0' THEN vert_switch = 1
167	IF plt1dtyp EQ 'z' THEN vert_switch = 0
168
169	; read all data if netcdf output
170
171	IF debug_w THEN print, ' force_all_data_read in pltmap: ', force_all_data_read
172
173	IF cmd.out NE 'cdf' THEN BEGIN
174	all_data_str = ''
175	ENDIF ELSE BEGIN
176	all_data_str = ',/ALL_DATA'
177	END
178
179	; ===========
180	; Read data
181	; ===========
182
183	really_1m_st = 1
184	IF cmd.timave EQ '1y' AND strmid(cmd.plt, 0, 2) EQ 'st' THEN BEGIN
185	@densit_pltmap_read
186	ENDIF ELSE BEGIN
187	pltcmd = 'field = data_read(cmd,'''+hotyp+''','''+plttyp+''','+string(dimplot)+','+string(iover)+all_data_str+', ZMTYP = '''+cmd.plt+''')'
188	printf, nulhis, strcompress(pltcmd)
189	printf, nulhis, ' '
190	IF execute(pltcmd) EQ 0 THEN stop
191	ENDELSE
192
193	; data specific treatments
194	IF cmd.var EQ 'sla' THEN field.origin = 'diff'
195	IF cmd.var EQ 'sosstsst' AND cmd.exp EQ 'HadISST1' THEN BEGIN
196	; limit data to -1.8 C
197	idx = where(field.data EQ valmask)
198	field.data = field.data > (-1.8)
199	IF idx[0] NE -1 THEN field.data(idx) = valmask
200	ENDIF
201
202	IF n_elements(field.data) EQ 1 THEN return
203
204	; change cmd.var if macro
205	CASE STRMID(cmd.var, 0, 2) OF
206	'@@': BEGIN
207	cmd.var = field.name
208	cmd.var = STRMID(cmd.var, 2, 15)
209	IF debug_w THEN print, 'cmd.var : ', cmd.var
210	END
211	ELSE:
212	ENDCASE
213	IF vecplot EQ 1 THEN BEGIN ; vectors case
214	CASE STRMID(cmd.var2, 0, 2) OF
215	'@@': begin
216	cmd.var2 = STRMID(cmd.var2, 2, 15)
217	END
218	ELSE:
219	ENDCASE
220	ENDIF
221
222	; spectum plot
223	IF spectplt EQ 1 THEN plttyp = 'spec'
224	; wavelet plot
225	IF waveplt EQ 1 THEN plttyp = 'wavelet'
226
227	xlogax = ''
228
229	; density pojection plot (splot=1)
230	IF splot EQ 1 THEN BEGIN
231	; keep grid attributes
232	izminmesh_b = izminmesh
233	izmaxmesh_b = izmaxmesh
234	jpk_b = jpk
235	jpkglo_b = jpkglo
236	gdept_b = gdept
237	gdepw_b = gdepw
238	e3t_b = e3t
239	e3w_b = e3w
240	tmask_b = tmask
241	; IF cmd.var EQ 'vosigvol' THEN xlogax = ',xlog = 1'
242	; modify vertical grid -> sigma
243	n_sig = (sig_max - sig_min)/sig_del + 1
244	z_sig = sig_min+findgen(n_sig)*sig_del
245	izminmesh = 0
246	izmaxmesh = n_sig-1
247	jpk = long(izmaxmesh-izminmesh+1)
248	jpkglo = jpk
249	gdept = z_sig
250	gdepw = gdept
251	e3t = shift(gdept, 1)-gdept
252	e3t[0] = e3t[1]
253	e3w = e3t
254	prof1 = sig_min
255	prof2 = sig_max
256	; grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz ; added PDW 12/5/04
257	; tmask = mask ; added PDW 12/5/04
258	tmask = reform(reform(tmask[, , 0], jpi*jpj)#replicate(1, jpk), jpi, jpj, jpk)
259	domdef
260	END
261
262	IF debug_w THEN print, ' cmd after data_read in plt_map = ', cmd
263	;
264	; ======================
265	; 2. window attributes
266	; ======================
267	;
268
269	; find field attributes
270	; ---------------------
271
272	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}
273
274	; extremum
275	fldext = fld_pltext(cmd.var, plttyp, dimplot, hotyp)
276
277	fldatt.name = fldext.name
278	fldatt.min = fldext.min
279	fldatt.max = fldext.max
280	fldatt.homin = fldext.homin
281	fldatt.homax = fldext.homax
282	fldatt.min1d = fldext.min1d
283	fldatt.max1d = fldext.max1d
284	fldatt.dmax = fldext.dmax
285	fldatt.assos = fldext.assos
286
287	; interval + change of unit
288	fldint = fld_pltint(cmd.var, plttyp, dimplot, hotyp)
289
290	; help, fldint, /struct
291
292	fldatt.int = fldint.int
293	fldatt.mult = fldint.mult
294	fldatt.add = fldint.add
295	fldatt.unit = fldint.unit
296	fldatt.mid = fldint.mid
297
298	; print, 'plt_map 1 :', fldatt.mult, fldatt.add, fldatt.min, fldatt.max
299
300	IF fldatt.mult NE 1. OR fldatt.add NE 0. THEN field.units = fldatt.unit
301	fld = field.data*fldatt.mult + fldatt.add
302	IF (where(field.data EQ valmask))[0] NE -1 THEN $
303	fld[where(field.data EQ valmask)] = valmask
304
305	; field 2 if needed
306
307	IF plttyp EQ 'yfx' OR vecplot EQ 1 THEN BEGIN
308	fldatt2 = fldatt
309
310	fldext2 = fld_pltext(cmd.var2, plttyp, dimplot, hotyp)
311
312	fldatt2.name = fldext2.name
313	fldatt2.min = fldext2.min
314	fldatt2.max = fldext2.max
315	fldatt2.homin = fldext2.homin
316	fldatt2.homax = fldext2.homax
317	fldatt2.min1d = fldext2.min1d
318	fldatt2.max1d = fldext2.max1d
319	fldatt2.dmax = fldext2.dmax
320
321	; interval + change of unit
322	fldint2 = fld_pltint(cmd.var2, plttyp, dimplot, hotyp)
323
324	fldatt2.int = fldint2.int
325	fldatt2.mult = fldint2.mult
326	fldatt2.add = fldint2.add
327	fldatt2.unit = fldint2.unit
328	fldatt2.mid = fldint2.mid
329
330	IF fldatt2.mult NE 1. OR fldatt2.add NE 0. THEN field.units2 = fldatt2.unit
331	fld2 = field.data2*fldatt2.mult + fldatt2.add
332	IF (where(field.data2 EQ valmask))[0] NE -1 THEN $
333	fld2[where(field.data2 EQ valmask)] = valmask
334
335	ENDIF
336
337	; If running std dev change min/max
338	IF run_stddev GT 0 THEN BEGIN
339	fldatt.homin = 0.
340	fldatt.homax = fldext.dmax
341	ENDIF
342
343	; legend text
344	; -----------
345
346	; IF plttyp EQ "pltt" THEN BEGIN
347	CASE strmid(cmd.trend, 0, 1) OF
348	'1': trd_txt = ' trend'
349	'2': trd_txt = ' drift'
350	'3': trd_txt = ' inverse trend'
351	'4': trd_txt = ' anomaly'
352	'5': trd_txt = ' filter'
353	'7': trd_txt = ' month sel'
354	ELSE: trd_txt = ''
355	ENDCASE
356	; ENDIF ELSE trd_txt = ''
357	IF field_int EQ 1 AND dimplot EQ 1 THEN BEGIN
358	trd_txt = trd_txt+' integral'
359	ENDIF
360
361	varlegend = field.legend+trd_txt+' ('+field.units+')'
362
363	CASE plttyp OF
364	'yfx': varlegend2 = field.legend2+trd_txt+' ('+field.units2+')'
365	ELSE: BEGIN
366	IF run_stddev GT 0 THEN trd_txt = trd_txt+' running Std Dev ['+string(run_stddev, format = '(I3)')+']'
367	END
368	ENDCASE
369
370	; date text
371
372	CASE strmid(cmd.timave, 0, 2) OF
373	'1m': BEGIN
374	mn = def_month(cmd.timave, cmd.date1)
375	IF strmid(cmd.timave, 0, 3) EQ '1mm' THEN $
376	date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $
377	ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2)
378	END
379	'3m': BEGIN
380	mn = def_month(cmd.timave, cmd.date1)
381	IF strmid(cmd.timave, 0, 3) EQ '3mm' THEN $
382	date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $
383	ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2)
384	END
385	ELSE:date_txt = cmd.timave+' '+cmd.date1
386	ENDCASE
387
388	CASE plttyp OF
389	'pltt':BEGIN
390	date_txt = cmd.timave
391	IF strmid(cmd.timave, 1, 2) EQ 'mm' THEN $
392	date_txt = cmd.timave+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')'
393	END
394	'yfx': BEGIN
395	IF hotyp NE '-' THEN date_txt = cmd.timave+' '+cmd.date1+' '+cmd.spec
396	END
397	'wavelet': BEGIN
398	date_txt = cmd.timave+' Wavelet'
399	IF strmid(cmd.timave, 1, 2) EQ 'mm' THEN $
400	date_txt = cmd.timave+' Wavelet ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')'
401	END
402	ELSE:
403	ENDCASE
404
405	; min/max/int
406
407	; defined if fraction of x/y.range is added to plot domain
408
409	; IF free_1d_minmax EQ 'no' THEN fraction = 0. ELSE fraction = 1.0
410	fraction = 0.
411
412	CASE dimplot OF
413	1: BEGIN
414	CASE plttyp OF
415	'pltt': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END
416	'wavelet': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END
417	'plt1d': BEGIN & minc = fldatt.min1d & maxc = fldatt.max1d & END
418	'yfx': BEGIN ; y=f(x) scatter plot
419	IF long(strmid(cmd.trend, 0, 1)) GT 0 THEN BEGIN
420	IF run_stddev EQ 0 THEN BEGIN
421	minc = -fldatt.dmax & maxc = fldatt.dmax
422	minc2 = -fldatt2.dmax & maxc2 = fldatt2.dmax
423	ENDIF ELSE BEGIN
424	minc = 0. & maxc = fldatt.dmax
425	minc2 = 0. & maxc2 = fldatt2.dmax
426	ENDELSE
427	ENDIF ELSE BEGIN
428	minc = fldatt.min1d & maxc = fldatt.max1d
429	minc2 = fldatt2.min1d & maxc2 = fldatt2.max1d
430	ENDELSE
431	END
432	ELSE: BEGIN & minc = 0. & maxc = 0. & END
433	ENDCASE
434	fldatt.int = !VALUES.F_NAN
435	END
436	2: BEGIN
437	CASE plttyp OF
438	'pltt': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END
439	ELSE: BEGIN & minc = fldatt.min & maxc = fldatt.max & END
440	ENDCASE
441	END
442	ENDCASE
443
444	IF cmd.var NE fld_prev OR cmd.trend NE '0' THEN BEGIN
445	; field.origin EQ 'diff' OR
446	print, ''
447	CASE dimplot OF
448	2: BEGIN
449	print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max/int=', $
450	minc, maxc, fldatt.int
451	IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $
452	print, ' --> Modify field using : ', 'mult/add=' , $
453	fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit
454	END
455	1: BEGIN
456	print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max=',minc, maxc
457	IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $
458	print, ' --> Modify field using : ', 'mult/add=' , $
459	fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit
460	END
461	0: BEGIN
462	print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max=',minc, maxc
463	IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $
464	print, ' --> Modify field using : ', 'mult/add=' , $
465	fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit
466	print, ' '+cmd.var2, ' plot attributes ('+plttyp+') : min/max=',minc2, maxc2
467	IF ( fldatt2.mult NE 1.0 OR fldatt2.add NE 0.0 ) THEN $
468	print, ' --> Modify field using : ', 'mult/add=' , $
469	fldatt2.mult, fldatt2.add, ' to obtain : ', fldatt2.unit
470	END
471	ELSE:
472	ENDCASE
473	print, ''
474
475	ENDIF
476	IF ( fldatt.mult EQ -1.0) THEN BEGIN
477	temp_m = minc
478	minc = -maxc
479	maxc = -temp_m
480	ENDIF
481	IF finite(fldatt.int) EQ 0 THEN BEGIN
482	intcmd = ''
483	colbarfor = ''
484	fmt = '(f6.1)'
485	ENDIF ELSE BEGIN
486	intcmd = ',int = '+string(fldatt.int)
487	IF fldatt.int LT 0.1 THEN BEGIN
488	fmt = '(f6.2)'
489	ENDIF ELSE IF fldatt.int LT 1 THEN BEGIN
490	fmt = '(f5.1)'
491	ENDIF ELSE BEGIN
492	fmt = '(f5.0)'
493	ENDELSE
494	IF long(fldatt.int) NE 0 THEN BEGIN
495	IF fldatt.int/long(fldatt.int) NE 1 THEN fmt = '(f5.1)'
496	ENDIF
497	colbarfor = ', format = '''+fmt+''''
498	ENDELSE
499
500	; if window > 1 or overlay > 1 noerase
501
502	IF win[2] NE 1 OR iover GT 1 THEN BEGIN
503	noerase = 1
504	ENDIF ELSE BEGIN
505	noerase = 0
506	ENDELSE
507
508	; choose window number
509
510	CASE multi_win OF
511	1: window_number = ', window='+string(iplot)
512	ELSE: window_number = ''
513	ENDCASE
514
515	; color label control
516
517	labstr = ''
518	IF iover EQ 1 THEN BEGIN
519	nofill = 1-shading
520	ENDIF ELSE BEGIN
521	nofill = 1
522	ENDELSE
523
524	nocoltxt = ''
525	IF nofill EQ 1 THEN nocoltxt = ',NOCOULEUR=1'
526
527	IF dimplot NE 2 THEN nofill = 1
528
529	nocolorbar = 0
530	IF col_palette EQ 'no' THEN nocolorbar = 1
531	IF pal_type EQ '2dom' THEN nocolorbar = 1
532	IF iover GT 1 THEN nocolorbar = 1
533	IF dimplot NE 2 THEN nocolorbar = 1
534
535	IF cmd.out NE 'cdf' THEN BEGIN
536	readpal = 0
537	IF nofill EQ 0 AND dimplot GT 1 THEN readpal = 1
538	IF waveplt EQ 1 THEN readpal = 2
539	IF field.origin EQ 'div' AND plttyp NE 'pltt' THEN readpal = 3
540	;; Necessary for correlation plots (overlay of 2D fields)
541	IF dimplot GT 1 AND nover GT 1 AND iover GT 1 THEN readpal = 1
542	IF readpal GE 1 THEN BEGIN
543	lec_pal_gmt, cmd.var, c_anot_str, fmt, found, readpal
544	colbarfor = ', format = '''+fmt+''''
545	IF found EQ 1 THEN BEGIN
546	labstr = ',label=3'
547	; if mincmd/naxcmd not defined then use one from palette
548	; if define take min of 2 mins and max of 2 maxs
549	min_palette = min(levels_gmt)
550	max_palette = max(levels_gmt)
551	IF finite(minc) NE 0 THEN minc = min([minc, min_palette]) ELSE minc = min_palette
552	IF finite(maxc) NE 0 THEN maxc = max([maxc, max_palette]) ELSE maxc = max_palette
553	print, ' '+cmd.var, ' plot attributes ('+plttyp+') modified via color palette: min/max=', $
554	minc, maxc
555	ENDIF
556	ENDIF
557	ENDIF
558	colbar = colbarfor
559
560	mincmd = ''
561	maxcmd = ''
562
563	IF finite(minc) NE 0 THEN BEGIN
564	mincmd = ','+string(minc)
565	maxcmd = ','+string(maxc)
566	ENDIF
567
568	; contour annotation + other color bar controls
569
570	IF n_elements(c_anot_str) EQ 0 THEN BEGIN
571	c_anot = ''
572	ENDIF ELSE BEGIN
573	c_anot = ',c_annotation=c_anot_str'
574	; sampling of colorbar annotation
575	IF col_palette EQ 'yes' THEN BEGIN
576	sampling = win[0]*win[1]
577	; idx_cb = lindgen((ncont_gmt-2)/sampling)*sampling+1
578	idx_cb0 = lindgen(ncont_gmt-1)+1
579	idx_cb = lindgen(ncont_gmt-2)+2
580	idx_cb1 = lindgen(ncont_gmt-2)+1
581	colbar = colbar+',discret=coul_gmt[idx_cb1],div=n_elements(idx_cb0)-1,cb_label=levels_gmt[idx_cb0]'
582	; c_anot = ',c_annotation=levels_gmt[idx_cb0]'
583	ENDIF
584	ENDELSE
585	;
586	; vertical grid type
587	;
588	IF debug_w THEN print, 'splot=', splot
589	CASE mesh_type OF
590	'oce': BEGIN
591	IF splot EQ 1 THEN BEGIN
592	type_yz = ',type_yz =''sigma'''
593	zoom_txt = ', zoom='+string(sig_max)
594	ENDIF ELSE BEGIN
595	type_yz = ',type_yz=''m'''
596	zoom_txt = ', zoom='+string(zoom_z)
597	ENDELSE
598	END
599	'atm': BEGIN & type_yz = ',type_yz=''hPa''' & zoom_txt = ', zoom='+string(hpa_max) & END
600	ELSE:type_yz = ''
601	ENDCASE
602
603	; continents color, thickness
604	; c_cont = 100
605	; continent fill
606	; real continent drawn
607	strcont =''
608	IF mesh_type NE 'oce' THEN BEGIN
609	IF cont_fill EQ 0 THEN strcont = ',/CONT_NOFILL, nite=0'
610	IF cont_real GE 1 THEN strcont = strcont+', REALCONT = '+string(strtrim(cont_real, 2))+', COAST_THICK = 2'
611	ENDIF
612
613	; calendar type
614
615	cal_typ = ''
616	IF calendar_type GT 1 THEN cal_typ = ',ndayspm=calendar_type'
617
618	; titles options
619	; By default, SAXO puts a title and a subtitle for the first plot
620	; Then, if you overlay another field, title and subtitle are set to ''
621	CASE title_type OF
622	"T": tit_str = ',subtitle='''','+marge_option
623	"S": tit_str = ',title='''','+marge_option
624	"TS": tit_str = ','+marge_option
625	"off":tit_str = ',title='''',subtitle='''','+marge_option
626	ENDCASE
627
628	; fill_space option
629
630	filltxt = ''
631	IF fill_space EQ 1 THEN filltxt = ',/rempli'
632
633	; common command string to plots
634
635	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+', /MEMEINDICES'
636
637	; add contour_options not overlay
638
639	IF iover EQ 1 THEN com_strplt = com_strplt+contour_options
640
641	; run specific fixes
642
643	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
644	fld = -fld
645	print, ' *** WARNING: Multiply tauu by -1 ', cmd.exp, cmd.var, cmd.timave
646	ENDIF
647
648	IF debug_w THEN print, ' cmd before making output in plt_map = ', cmd
649
650
651	; ============================
652	; 3. make output (data/plot)
653	; ============================
654	;
655
656	CASE cmd.out OF
657	'data': write_data = long((iplot-1)*10+win[2]) ; write 1D data ascii in trends.pro
658	'tcdf': write_data = -long((iplot-1)*10+win[2]) ; write 1D data netcdf in trends.pro
659	'cdf': write_data = 0
660	ELSE: write_data = 0
661	ENDCASE
662
663	; make output
664
665	CASE cmd.out OF
666	'tv': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ tvnplot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
667	erase
668
669	CASE (size(fld))[0] OF
670	2: fldtv = fld
671	3: BEGIN
672	niveau = xquestion('3d data : which level ?', '1', /chkwidget)
673	fldtv = fld[, , niveau]
674	END
675	4: BEGIN
676	niveau = xquestion('3d data : which level ?', '1', /chkwidget)
677	timeplot = xquestion('4d data : which date ?', '1', /chkwidget)
678	fldtv = fld[, , niveau, timeplot]
679	END
680	ELSE:
681	ENDCASE
682
683	CASE cmd.grid OF
684	; 'T': tvnplot, fldtv*tmask
685	'T': tvplus, fldtv, min = fldatt.min, max = fldatt.max
686	'U': tvplus, fldtv, min = fldatt.min, max = fldatt.max
687	'V': tvplus, fldtv, min = fldatt.min, max = fldatt.max
688	'W': tvplus, fldtv, min = fldatt.min, max = fldatt.max
689	ELSE: tvplus, fldtv, min = fldatt.min, max = fldatt.max
690	ENDCASE
691	END
692
693	'cdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
694
695	; build structures containing grids
696	cdf_description = nc_build(cmd, field, field.direc, vargrid)
697	; build netCDF file (here if whole data)
698	IF hotyp EQ '-' OR hotyp EQ 'xyt' THEN BEGIN
699	fldcdf = {data:field.data, units:field.units, short_name:cmd.var, long_name:field.legend, missing_value:valmask, direc:field.direc}
700	file_out = cmd.var+'_'+strtrim(string(FORMAT = '(I2.2)', iplot), 2)+'.nc'
701	pltcmd ='nc_put, fldcdf, file_out, ncdf_db = hom_idl'+cdf_description
702	printf, nulhis, strcompress(pltcmd)
703	res = execute(pltcmd)
704	ENDIF
705	END
706	'tcdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ 1D netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
707	mask_z, fld, cmd, boite_plt1d, dimplot
708	fld = checkfield(fld, 'pltt', type = datyp.hotyp, boite = boite_plt1d)
709	IF debug_w THEN print, 'size(fld)', size(fld)
710	IF debug_w THEN print, 'boite_plt1d', boite_plt1d
711	IF cmd.trend GT 0 THEN BEGIN
712	fld = trends(fld, cmd.trend, datyp.hotyp)
713	add_txt = trd_txt
714	ENDIF ELSE print, 'You need to have a trend to make 1D netcdf output'
715	END
716	ELSE: BEGIN ; $$$$$$$$$$$$$$$$$$$$$$$$$$ make plot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
717
718	CASE plttyp OF
719	'plt': BEGIN
720	;
721	; Map plot
722	; --------
723	;
724	print, ' '+cmd.var+' data to plot min and max : ', $
725	min(fld(where (fld NE -valmask)), /nan), max(fld(where (fld NE valmask)), /nan)
726	print, ''
727
728	map_cmd = ''
729
730	niveau = 1
731
732	IF vecplot EQ 0 THEN BEGIN
733	; eddy field ? (if so, remove zonal mean)
734	IF strpos(cmd.plt, '@z') GT 1 THEN BEGIN
735	fldzm = moyenne(fld, 'x', boite=[20,380,-90,90], /NAN)
736	fldzm = replicate(1, (size(fld))(1))#fldzm
737	fld = fld-fldzm
738	filleg = ' Eddy '
739	ENDIF ELSE BEGIN
740	filleg = ''
741	ENDELSE
742
743	IF hotyp EQ 'xyt' THEN date_txt = 'monthly'
744
745	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
746	;; For 2D plots with overlays (correlation and proba for instance)
747	;; Scratch the titles and define contours style for the second plot
748	IF nover GT 1 AND iover GT 1 THEN BEGIN
749	contour_style = ',c_thick=2,c_linestyle=2,c_charthick=2'
750	CASE title_type OF
751	'T': pltcmd = pltcmd+',title='''''+contour_style
752	'S': pltcmd = pltcmd+',subtitle='''''+contour_style
753	'TS': pltcmd = pltcmd+',title='''',subtitle='''''+contour_style
754	'off': pltcmd = pltcmd+contour_style
755	END
756	ENDIF
757	ENDIF ELSE BEGIN
758
759	fld = {a:fld, g:strmid(cmd.grid, 0, 1)}
760	fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)}
761	fldv = {data1: fld, data2: fld2}
762	vargrid = strmid(cmd.grid, 2, 1)
763
764	pltcmd = 'ajoutvect,fldv'+vect_sample+com_strplt+strcont+map_cmd
765
766	ENDELSE
767
768	IF debug_w THEN print, strcompress(pltcmd)
769	printf, nulhis, strcompress(pltcmd)
770	res = execute(pltcmd(0))
771
772	END
773	'pltz': BEGIN
774	;
775	; Vertical section/mean plot
776	; --------------------------
777	;
778	g = gphit
779	t = tmask
780
781	; if already 1D, reform fld
782	sz = size(fld)
783	IF sz[0] EQ 2 THEN BEGIN
784	IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2])
785	ENDIF
786	IF sz[0] EQ 3 THEN BEGIN
787	IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3])
788	IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3])
789	IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2])
790	ENDIF
791
792	; mask fld
793	mask_z, fld, cmd, boite_pltz, dimplot, legz
794	printf, nulhis, ' boite_pltz = ', boite_pltz
795
796	print, ' '+cmd.var+' data to plot min and max : ', $
797	min(fld(where (fld NE valmask)), /NAN), max(fld(where (fld NE valmask)), /NAN)
798	print, ''
799	IF vecplot EQ 0 THEN BEGIN
800
801	IF cmd.var EQ 'vozonbsf' THEN xindstr = ', /xindex' ELSE xindstr = ''
802
803	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
804
805
806	ENDIF ELSE BEGIN
807
808	fld = {a:fld, g:strmid(cmd.grid, 0, 1)}
809	fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)}
810	fldv = {data1: fld, data2: fld2}
811	vargrid = strmid(cmd.grid, 2, 1)
812
813	pltcmd = 'ajoutvectz,fldv'+com_strplt+strcont+',type='''+strmid(cmd.plt, 0, 2)+''', boite=boite_pltz'
814
815	ENDELSE
816
817	IF debug_w THEN print, pltcmd
818	printf, nulhis, strcompress(pltcmd)
819	res = execute(pltcmd(0))
820
821
822	; overlay bowl if sig_bowl=1
823	IF sig_bowl EQ 1 THEN begin
824	; define line color, thickness and type
825	iover = 2
826	overc = overlay_type(iover, dimplot)
827	; type of latitude axis
828	IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $
829	sin = ',/sin' ELSE sin = ''
830	plt1dtyp = strmid(pltztyp, 0, 1)
831	noerase = 1
832	com_strplt = ', NOERASE = '+string(noerase)
833	pltcmd = 'plt1d,field.bowl,/'+plt1dtyp+', boite=boite_pltz'+overc+sin+com_strplt
834	IF debug_w THEN print, pltcmd
835
836	res = execute(pltcmd(0))
837	ENDIF
838
839	gphit = g
840	tmask = t
841
842	END
843	'plt1d': BEGIN
844	;
845	; 1D, x, y, z plot
846	; ----------------
847	;
848	g = gphit
849	t = tmask
850	; if already 1D, reform fld
851	sz = size(fld)
852	IF sz[0] EQ 2 THEN BEGIN
853	IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2])
854	ENDIF
855	IF sz[0] EQ 3 THEN BEGIN
856	IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3])
857	IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3])
858	IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2])
859	ENDIF
860	; mask fld
861	mask_z, fld, cmd, boite_plt1d, dimplot, legz
862
863	IF debug_w THEN print, boite_plt1d
864	niveau = 1
865	; define line color, thickness and type
866	overc = overlay_type(iover, dimplot)
867	; type of latitude axis
868	IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $
869	sin = ',/sin' ELSE sin = ''
870
871	print, ' '+cmd.var+' data to plot min and max : ', $
872	min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask)))
873	print, ''
874	take_log = ',take_log=0'
875	IF cmd.var EQ 'vosigvol' AND n_elements(str_sep(cmd.exp,'-')) EQ 1 THEN BEGIN ; don't take log if difference plot
876	take_log = ',take_log=1'
877	ENDIF
878	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
879	printf, nulhis, 'boite_plt1d=',boite_plt1d
880	printf, nulhis, strcompress(pltcmd)
881	IF debug_w THEN print, pltcmd
882	res = execute(pltcmd(0))
883
884	@legend_overlay
885
886	gphit = g
887	tmask = t
888
889	END
890	'pltt': BEGIN
891	;
892	; hovmoeller
893	; -----------
894	;
895	g = gphit
896	t = tmask
897
898	; mask fld
899	mask_z, fld, cmd, boite_pltt, dimplot, legz
900	; define line color, thickness and type
901	overc = overlay_type(iover, dimplot)
902	; define additional text for pltt
903	@add_txt_pltt
904	; repeat for seasonal mean
905	CASE cmd.timave OF
906	'1mm': rep_txt = ',repeat_c=nb_cycles'
907	ELSE: rep_txt = ''
908	ENDCASE
909
910	print, ' '+cmd.var+' data to plot min and max : ', $
911	min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask)))
912	print, ''
913	vardate = 'toto'
914	grille, mask, glam, gphi, gdep, nx, ny,nz
915	IF st_rms EQ 1 THEN BEGIN ; time series of rms deviation on a sigma surface
916	@densit_pltmap_plot
917	ENDIF ELSE BEGIN
918	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'
919	ENDELSE
920
921	printf, nulhis, 'boite_pltt=',boite_pltt
922	printf, nulhis, strcompress(pltcmd)
923	IF debug_w THEN print, pltcmd
924	res = execute(pltcmd(0))
925	IF iover EQ 4 THEN print, 'There might be a pb with the legends !'
926	; legend if plot=t
927	IF hotyp EQ 't' THEN BEGIN
928	; positions of the legend depend on nb_cycles
929	; it is different from 1 only with 1mm case
930	IF cmd.timave NE '1mm' THEN nb_cycles = 1
931	@legend_overlay
932	ENDIF
933	gphit = g
934	tmask = t
935	END
936
937	'yfx': BEGIN
938	;
939	; Scatter plot y=f(x)
940	; -------------------
941	;
942	t = tmask
943	boxx = ''
944	boxy = ''
945	add_txt = ''
946	; arrange/average data if needed
947	; space and time plots fld
948	IF strpos(cmd.plt, 'xy_') NE -1 THEN BEGIN
949	mask_z, fld, cmd, boite_plt, dimplot, boxx
950	boxx = ' ['+boxx+']'
951	print, ' Averaging (space) '+cmd.var+' in'+boxx
952	fld = checkfield(fld, 'plt', type = 'xy', boite = boite_plt)
953	ENDIF
954	IF datyp.hotyp NE '-' THEN BEGIN
955	mask_z, fld, cmd, boite_plt1d, dimplot, boxx
956	boxx = ' ['+boxx+']'
957	vargrid = vargrid1
958	IF debug_w THEN print, ' domdef and vargrid : ', boite_plt1d[0:3], vargrid
959	domdef, boite_plt1d[0:3]
960	print, ' Averaging (time serie plot) cmd.var '+cmd.var+' in'+boxx
961	IF debug_w THEN print, 'fld size', size(fld)
962	fld = checkfield(fld, 'plt1d', type = datyp.hotyp, boite = boite_plt1d)
963	IF cmd.trend GT 0 THEN BEGIN
964	fld_flag = 1
965	fld = trends(fld, cmd.trend, datyp.hotyp)
966	add_txt = trd_txt
967	ENDIF
968	ENDIF
969	; space and time plots fld2
970	IF sw_diffg EQ 1 THEN BEGIN
971	jptb = jpt
972	IF debug_w THEN print, ' calling def_grid, cmd2 dans plt_map '
973	def_grid, cmd2
974	jpt = jptb
975	ENDIF
976	IF strpos(cmd2.plt, 'xy_') NE -1 THEN BEGIN
977	mask_z, fld2, cmd2, boite_plt2, dimplot, boxy
978	boxy = ' ['+boxy+']'
979	print, ' Averaging (space) '+cmd2.var+' in'+boxy
980	fld2 = checkfield(fld2, 'plt', type = 'xy', boite = boite_plt2)
981	ENDIF
982	IF datyp2.hotyp NE '-' THEN BEGIN
983	run_stddev = 0
984	IF strpos(cmd2.plt, '@r') GE 1 THEN BEGIN
985	run_stddev = strmid(cmd2.plt, strpos(cmd2.plt, '@r')+2, strlen(cmd2.plt))
986	cmd2.plt = strmid(cmd2.plt, 0, strpos(cmd2.plt, '@r'))
987	ENDIF
988	mask_z, fld2, cmd2, boite_plt1d2, dimplot, boxy
989	boxy = ' ['+boxy+']'
990	vargrid = vargrid2
991	IF debug_w THEN print, ' domdef and vargrid : ', boite_plt1d2[0:3], vargrid
992	domdef, boite_plt1d2[0:3]
993	print, ' Averaging (time serie plot) cmd2.var '+cmd2.var+' in'+boxy
994	IF debug_w THEN print, ' fld2 size', size(fld2)
995	fld2 = checkfield(fld2, 'plt1d', type = datyp2.hotyp, direc = 'xy', boite = boite_plt1d2, /timearray)
996	IF cmd2.trend GT 0 THEN BEGIN
997	fld_flag = 2
998	fld2 = trends(fld2, cmd2.trend, datyp2.hotyp)
999	add_txt = trd_txt
1000	ENDIF
1001	ENDIF
1002	IF sw_diffg EQ 1 THEN BEGIN
1003	sw_diffg = 0
1004	jptb = jpt
1005	IF debug_w THEN print, ' calling def_grid, cmd dans plt_map '
1006	def_grid, cmd
1007	jpt = jptb
1008	ENDIF
1009	coeff = linfit(fld2, fld, CHISQ = linerr)
1010	print, ' Linfit coef + error (full serie)=', coeff(0), coeff(1), linerr
1011
1012	; plot domain
1013	boxyfx = def_box(cmd.plt, dimplot, legz, time_stride)
1014	; define line color and type
1015	overc = overlay_type(iover, dimplot)
1016	vardate = date_txt
1017	varname = varlegend
1018	varunit = ' '+cmd.var+boxx+' =f('+cmd.var2+boxy+') '+add_txt
1019	pltcmdstd = 'pltsc,fld,fld2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+overc+com_strplt
1020
1021	IF hotyp EQ 't' THEN BEGIN ; time scatter plot
1022	IF mean_sc_only EQ 0 OR mean_sc_only EQ 4 THEN BEGIN
1023	; number of time colors
1024	IF strpos(symbol_families, 'x') NE -1 THEN BEGIN
1025	coding = str_sep(symbol_families, 'x')
1026	ncolors = long(coding(0))
1027	ntimes = long(coding(1))
1028	ENDIF ELSE BEGIN
1029	ncolors = long(symbol_families)
1030	ntimes = 1
1031	ENDELSE
1032	IF ncolors GE 2 THEN BEGIN ; time color coding
1033	ic = 0
1034	jpto = jpt
1035	jpt = jpt/ncolors
1036	lincoef = fltarr(ncolors)
1037	linerro = fltarr(ncolors)
1038	linprob = fltarr(ncolors)
1039	WHILE ic LE ncolors-1 DO BEGIN
1040	idx0 = (floor(findgen(jpto/(ncolorsntimes)))ncolorsntimes)+ic(ntimes)
1041	idx = idx0
1042	jl = 1
1043	WHILE jl LE ntimes-1 DO BEGIN
1044	idx = [idx, idx0+jl]
1045	jl = jl+1
1046	ENDWHILE
1047	fldp = fld(idx)
1048	fldp2 = fld2(idx)
1049	coeff = linfit(fldp2(where (fldp2() GT 0.)), fldp(where (fldp2() GT 0.)), CHISQ = errlin, PROB = prblin)
1050	linerro(ic) = errlin
1051	linprob(ic) = prblin
1052	lincoef(ic) = coeff(1)
1053	; print, ' Period '+strtrim(string(ic), 2)+' Linear fit slope, chisq, prob =', coeff(1)1000., errlin1000., prblin
1054	IF mean_sc_only EQ 0 THEN BEGIN
1055	pltcmd = 'pltsc,fldp,fldp2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d='+string(ic)+',COLOR='+string(symbol_color(ic))+', STY1D='+string(symbol_style(ic)-1)+', SYMSIZE='+string(symbol_size)
1056	printf, nulhis, strcompress(pltcmd)
1057	IF debug_w THEN print, ' ',pltcmd
1058	res = execute(pltcmd)
1059	ENDIF
1060	IF mean_sc_only EQ 4 AND ic EQ ncolors-1 THEN BEGIN
1061	print, ' Linfit slope + error (June-Nov)=', mean(lincoef(5:10)), mean(linerro(5:10))
1062	print, ' Linfit slope + error (Dec-May)=', mean([lincoef(0:4), lincoef(11)]), mean([linerro(0:4), linerro(11)])
1063	vardate = 'toto'
1064	varname = 'Coupling strength '+cmd.date1+' '+cmd.spec
1065	jpt = ncolors
1066	time = lindgen(12)301+ $
1067	julday(1,1,01, ndayspm = calendar_type)
1068	hotyp = 't'
1069	minmax = ',0,25'
1070	pltcmd = 'pltt,lincoef*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=0,COLOR=1, thick=4, STY1D=0'
1071
1072	printf, nulhis, strcompress(pltcmd)
1073	IF debug_w THEN print, ' ', pltcmd
1074	res = execute(pltcmd)
1075	pltcmd = 'pltt,(lincoef-linerro)*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=0'
1076	printf, nulhis, strcompress(pltcmd)
1077	IF debug_w THEN print, ' ',pltcmd
1078	res = execute(pltcmd)
1079	pltcmd = 'pltt,(lincoef+linerro)*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=0'
1080	printf, nulhis, strcompress(pltcmd)
1081	IF debug_w THEN print, ' ',pltcmd
1082	res = execute(pltcmd)
1083	ENDIF
1084	ic = ic + 1
1085	ENDWHILE
1086	ENDIF ELSE BEGIN ; one color
1087	IF debug_w THEN print, ' ',pltcmdstd
1088	res = execute(pltcmdstd(0))
1089	ENDELSE
1090	ENDIF
1091	ENDIF ELSE BEGIN ; standard scatter plot
1092	IF debug_w THEN print, ' ',pltcmdstd
1093	res = execute(pltcmdstd(0))
1094	ENDELSE
1095
1096	; add mean SC in 4x3 plots
1097
1098	IF symbol_families EQ '4x3' THEN BEGIN
1099
1100	IF cmd.trend EQ '412' THEN BEGIN
1101	fldrem_t1 = fldrem_t1-mean(fldrem_t1)
1102	fldrem_t2 = fldrem_t2-mean(fldrem_t2)
1103	ENDIF ELSE BEGIN
1104	running = 12L
1105	lenght = (size(fld))[1]
1106	fldrem_t1 = fltarr(running)
1107	fldrem_t2 = fltarr(running)
1108	FOR t1 = 0, running-1 DO BEGIN
1109	fldrem_t1(t1) = mean(fld(long(findgen(lenght/running)*running+t1)))
1110	fldrem_t2(t1) = mean(fld2(long(findgen(lenght/running)*running+t1)))
1111	ENDFOR
1112	ENDELSE
1113
1114	print, ' Seasonal cycle var/stddev fld1 ', (moment(fldrem_t1))[1], sqrt((moment(fldrem_t1))[1])
1115	print, ' Seasonal cycle var/stddev fld2 ', (moment(fldrem_t2))[1], sqrt((moment(fldrem_t2))[1])
1116
1117	fldrem_t1 = [fldrem_t1, fldrem_t1(0)]
1118	fldrem_t2 = [fldrem_t2, fldrem_t2(0)]
1119	sw_ov1d = mean_sc_only
1120	IF mean_sc_only EQ 2 AND cmd.trend EQ '412' THEN BEGIN
1121	; SC of std dev
1122	stat = spectra(fld, jpt/12, 20, 15, 0)
1123	stat2 = spectra(fld2, jpt/12, 20, 15, 0)
1124	print, ' Stddev of monthly stddevs fld', sqrt((moment(stat.sc_std-(moment(stat.sc_std))[0]))[1])
1125	print, ' Stddev of monthly stddevs fld2', sqrt((moment(stat2.sc_std-(moment(stat2.sc_std))[0]))[1])
1126	stdsc = [stat.sc_std, stat.sc_std(0)]
1127	stdsc2 = [stat2.sc_std, stat2.sc_std(0)]
1128	pltcmd = 'pltsc,stdsc,stdsc2,0. ,maxc,0.,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d=1-min(1,sw_ov1d), STY1D=-1, THICKN=3, SYMSIZE='+string(symbol_size)+', FRACTION=fraction'
1129	printf, nulhis, strcompress(pltcmd)
1130	IF debug_w THEN print, pltcmd
1131	res = execute(pltcmd(0))
1132	pltcmd = 'xyouts,stdsc2-(maxc2)/(.3win(0)win(1)),stdsc,string(long(findgen(12))+1),charsize=1.3,alignment=0.5'
1133	printf, nulhis, strcompress(pltcmd)
1134	IF debug_w THEN print, pltcmd
1135	res = execute(pltcmd(0))
1136	ENDIF ELSE BEGIN
1137	ov1d_val = 0
1138	IF mean_sc_only EQ 0 THEN ov1d_val = 1
1139	pltcmd = 'pltsc,fldrem_t1,fldrem_t2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d=ov1d_val, STY1D=-1, THICKN=5, SYMSIZE='+string(symbol_size)+', FRACTION=fraction'
1140	printf, nulhis, strcompress(pltcmd)
1141	IF debug_w THEN print, pltcmd
1142	res = execute(pltcmd(0))
1143	IF mean_sc_only EQ 1 THEN BEGIN
1144	; add month info win(0)*win(1)
1145	pltcmd = 'xyouts,fldrem_t2-(maxc2-minc2)/(0.5win(0)win(1)),fldrem_t1,string(long(findgen(12))+1),charsize=1.3,alignment=0.5'
1146	printf, nulhis, strcompress(pltcmd)
1147	IF debug_w THEN print, pltcmd
1148	res = execute(pltcmd(0))
1149	ENDIF
1150	ENDELSE
1151	ENDIF
1152
1153	; overplot sigma contours in T-S plane
1154	IF (cmd.var EQ 'sosstsst' and cmd.var2 EQ 'sosaline') OR (cmd.var EQ 'votemper' and cmd.var2 EQ 'vosaline') THEN BEGIN
1155	IF noerase EQ 0 THEN BEGIN
1156	min_t = minc-2.0
1157	max_t = maxc+2.0
1158	min_s = minc2-2.0
1159	max_s = maxc2+2.0
1160	delta_t = 0.5
1161	delta_s = 0.05
1162	t_vec = findgen((max_t-min_t)/delta_t+1)*delta_t+min_t
1163	s_vec = findgen((max_s-min_s)/delta_s+1)*delta_s+min_s
1164	t_ = t ; since t is already defined
1165	t = t_vec
1166	FOR i = 1, (max_s-min_s)/delta_s DO t = [[t], [t_vec]]
1167	s = s_vec
1168	FOR i = 1, (max_t-min_t)/delta_t DO s = [[s], [s_vec]]
1169	s = transpose(s)
1170	sr=sqrt(abs(s))
1171	r1=((((6.536332E-9t-1.120083E-6)t+1.001685E-4)*t $
1172	-9.095290E-3)t+6.793952E-2)t+999.842594
1173	r2=(((5.3875E-9t-8.2467E-7)t+7.6438E-5)t-4.0899E-3)t+8.24493E-1
1174	r3=(-1.6546E-6t+1.0227E-4)t-5.72466E-3
1175	rhopn = ( ( 4.8314E-4s + r3sr +r2)*s +r1) -1000.0
1176	min_sig = round(min(rhopn)+0.5)-1.
1177	max_sig = round(max(rhopn)-0.5)+1.
1178	delta_sig = 1. ; define contour interval here!
1179	levels = findgen((max_sig-min_sig)/delta_sig+1)*delta_sig+min_sig
1180	labels = levels
1181	labels(*) = 1
1182	contour, rhopn, s, t, /overplot, levels = levels, c_labels = labels, c_linestyle = 0
1183	t = t_
1184	ENDIF
1185	ENDIF
1186
1187	printf, nulhis, 'boite_yfx=',boxyfx
1188	printf, nulhis, strcompress(pltcmd)
1189
1190	tmask = t
1191	; force read grid for next window
1192	cmd_prev.grid = 'toto'
1193
1194	END
1195	'spec': BEGIN
1196	;
1197	; Spectrum plot y=fft(x)
1198	; ---------------------
1199	;
1200	t = tmask
1201	; plot domain
1202	mask_z, fld, cmd, boite_pltspec, dimplot, legspec
1203	; define line color and type
1204	overc = overlay_type(iover, dimplot)
1205	vardate = date_txt
1206	varname = varlegend
1207	varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec
1208
1209	; make mean in box
1210	fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex)
1211	; apply trends
1212	IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't')
1213
1214	; make spectrum
1215	time_interval = time[1]-time[0]
1216
1217	print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360
1218	print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360
1219	print, ' spectrum window (years)/chunks = ', spec_win/360, long(time(jpt-1)-time(0)+time_interval)/spec_win
1220	print, ' '
1221
1222	; number of time windowsš
1223	nt_win = long((time(jpt-1)-time(0)+time_interval)/spec_win)
1224	idx_win_size = long(spec_win/time_interval)
1225	IF idx_win_size*nt_win NE jpt THEN BEGIN
1226	print, ' *** Warning : spec_win must divide lenght of time serie ', idx_win_size, jpt
1227	return
1228	ENDIF
1229	; mean of idx_win_size chunks
1230	spect = spectrum(fld[0:idx_win_size-1], time_interval)
1231	FOR it = 2, nt_win DO BEGIN
1232	idx1 = idx_win_size*(it-1)
1233	idx2 = idx1 + idx_win_size - 1
1234	specc = spectrum(fld[idx1:idx2], time_interval)
1235	spect = spect + specc
1236	ENDFOR
1237	spect = spect/nt_win
1238	; build new time array
1239	idx = where (spect[0, *] LE max_spec)
1240	jpt = n_elements(idx)
1241	fld = reverse(reform(spect[1, idx], jpt))
1242	time = findgen(jpt)
1243	time = reverse(reform(spect[0, idx], jpt))
1244	IF max(time) GT 20*360 THEN BEGIN
1245	time = time/360
1246	xtitle = 'Years'
1247	ENDIF ELSE IF max(time) GT 5*360 THEN BEGIN
1248	time = time/30
1249	xtitle = 'Months'
1250	ENDIF ELSE xtitle = 'Days'
1251	; time/space filter ?
1252	IF strpos(cmd.plt, '@f') GT 1 THEN BEGIN
1253	filter = long(strmid(cmd.plt, strpos(cmd.plt, '@f')+3, strlen(cmd.plt)-strpos(cmd.plt, '@f')-3))
1254	fildirec = strmid(cmd.plt, strpos(cmd.plt, '@f')+2, 1)
1255	IF fildirec EQ 't' THEN BEGIN
1256	xtitle = xtitle+' (filter '+strtrim(string(filter), 2)+' warning : discrete filter)'
1257	fld = smooth(fld, filter)
1258	fld[0:filter/2-1] = 0.
1259	fld[(size(fld))[1]-filter/2-1:(size(fld))[1]-1] = 0.
1260	ENDIF
1261	ENDIF
1262	; plot min/max
1263	min1 = min(fld)
1264	max1 = max(fld)
1265	min2 = 0
1266	max2 = max(time)
1267	!x.range = [min2-abs(max2-min2)/5.,max2+abs(max2-min2)/5.]
1268	!y.range = [min1-abs(max1-min1)/5.,max1+abs(max1-min1)/5.]
1269
1270	; plot
1271	ytitle = 'Power spectrum (window='+strtrim(string(spec_win/360), 2)+'y)'
1272	pltcmd = 'splot,time,fld,xstyle=1,ystyle=1,title=varunit,xtitle=xtitle,ytitle=ytitle'+overc+com_strplt
1273	printf, nulhis, 'boite_pltspec=',boite_pltspec
1274	printf, nulhis, strcompress(pltcmd)
1275	IF debug_w THEN print, pltcmd
1276	res = execute(pltcmd(0))
1277
1278	tmask = t
1279
1280	END
1281	'wavelet': BEGIN
1282	;
1283	; Wavelet plot
1284	; ------------
1285	;
1286	t = tmask
1287	; plot domain
1288	mask_z, fld, cmd, boite_pltspec, dimplot, legspec
1289	; define line color and type
1290	vardate = date_txt
1291	varname = varlegend + ' '+date_txt
1292	varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec
1293	; make mean in box
1294	fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex)
1295	; apply trends
1296	IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't')
1297
1298	; make spectrum
1299	time_interval = time[1]-time[0]
1300
1301	print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360
1302	print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360
1303	print, ' '
1304	; make wavelet
1305	wave = wavelet(fld,time_interval,period=period,coi=coi,/pad,signif=signif)
1306	power = abs(wave)^2
1307	nscale = n_elements(period)
1308	period = period/365 ; to have period in years
1309	; compute significance
1310	signif = rebin(transpose(signif),jpt,nscale)
1311	; plot wavelet
1312	mincmd = ','+string(min(power))
1313	maxcmd = ','+string(max(power))
1314	boite_pltspec = [0, 0, min(period), max(period)]
1315	; domdef, boite_pltspec
1316	lat1r = lat1
1317	lat2r = lat2
1318	lat1 = 0
1319	lat2 = max_spec/360
1320	key_onearth = 0
1321	pltcmd = 'plttg,power,period'+mincmd+maxcmd+intcmd+',boite=boite_pltspec,reverse_y=1,nocontour=1'+com_strplt
1322	printf, nulhis, 'boite_pltspec=',boite_pltspec
1323	printf, nulhis, strcompress(pltcmd)
1324	IF debug_w THEN print, pltcmd
1325	res = execute(pltcmd(0))
1326	contour,abs(wave)^2/signif,time,period, /overplot,level=1.0,c_annot='95%'
1327	plot,time,coi/365, noclip = 0, /noerase
1328	tmask = t
1329	lat1 = lat1r
1330	lat2 = lat2r
1331	key_onearth = 1
1332	END
1333	ELSE: BEGIN
1334	print, ' unknown projection plot ', cmd.plt, ' / plot type = ', plttyp
1335
1336	END
1337	ENDCASE
1338
1339	END
1340
1341	ENDCASE
1342
1343	fld_prev = cmd.var
1344	;
1345	; reset incompatible options of plt_def
1346
1347	cmd.trend = trend_typp
1348
1349	;
1350	; reset vertical grid after density bining
1351
1352	IF splot EQ 1 THEN BEGIN
1353
1354	izminmesh = izminmesh_b
1355	izmaxmesh = izmaxmesh_b
1356	jpk = jpk_b
1357	jpkglo = jpkglo_b
1358	gdept = gdept_b
1359	gdepw = gdepw_b
1360	e3t = e3t_b
1361	e3w = e3w_b
1362	tmask = tmask_b
1363
1364	ENDIF
1365	IF debug_w THEN print, ' ...EXIT plt_map'
1366
1367	END

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

Download in other formats:

Original Format