Context Navigation

plt_map.pro @ 34

Last change on this file since 34 was 34, checked in by kolasinski, 16 years ago
Add x/ychartxt read from plt_def
File size: 51.8 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	; axis charsize option (x/ychartxt read from plt_def)
634
635	; common command string to plots
636
637	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)
638
639	; add contour_options not overlay
640
641	IF iover EQ 1 THEN com_strplt = com_strplt+contour_options
642
643	; run specific fixes
644
645	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
646	fld = -fld
647	print, ' *** WARNING: Multiply tauu by -1 ', cmd.exp, cmd.var, cmd.timave
648	ENDIF
649
650	IF debug_w THEN print, ' cmd before making output in plt_map = ', cmd
651
652
653	; ============================
654	; 3. make output (data/plot)
655	; ============================
656	;
657
658	CASE cmd.out OF
659	'data': write_data = long((iplot-1)*10+win[2]) ; write 1D data ascii in trends.pro
660	'tcdf': write_data = -long((iplot-1)*10+win[2]) ; write 1D data netcdf in trends.pro
661	'cdf': write_data = 0
662	ELSE: write_data = 0
663	ENDCASE
664
665	; make output
666
667	CASE cmd.out OF
668	'tv': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ tvnplot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
669	erase
670
671	CASE (size(fld))[0] OF
672	2: fldtv = fld
673	3: BEGIN
674	niveau = xquestion('3d data : which level ?', '1', /chkwidget)
675	fldtv = fld[, , niveau]
676	END
677	4: BEGIN
678	niveau = xquestion('3d data : which level ?', '1', /chkwidget)
679	timeplot = xquestion('4d data : which date ?', '1', /chkwidget)
680	fldtv = fld[, , niveau, timeplot]
681	END
682	ELSE:
683	ENDCASE
684
685	CASE cmd.grid OF
686	; 'T': tvnplot, fldtv*tmask
687	'T': tvplus, fldtv, min = fldatt.min, max = fldatt.max
688	'U': tvplus, fldtv, min = fldatt.min, max = fldatt.max
689	'V': tvplus, fldtv, min = fldatt.min, max = fldatt.max
690	'W': tvplus, fldtv, min = fldatt.min, max = fldatt.max
691	ELSE: tvplus, fldtv, min = fldatt.min, max = fldatt.max
692	ENDCASE
693	END
694
695	'cdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
696
697	; build structures containing grids
698	cdf_description = nc_build(cmd, field, field.direc, vargrid)
699	; build netCDF file (here if whole data)
700	IF hotyp EQ '-' OR hotyp EQ 'xyt' THEN BEGIN
701	fldcdf = {data:field.data, units:field.units, short_name:cmd.var, long_name:field.legend, missing_value:valmask, direc:field.direc}
702	file_out = cmd.var+'_'+strtrim(string(FORMAT = '(I2.2)', iplot), 2)+'.nc'
703	pltcmd ='nc_put, fldcdf, file_out, ncdf_db = hom_idl'+cdf_description
704	printf, nulhis, strcompress(pltcmd)
705	res = execute(pltcmd)
706	ENDIF
707	END
708	'tcdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ 1D netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
709	mask_z, fld, cmd, boite_plt1d, dimplot
710	fld = checkfield(fld, 'pltt', type = datyp.hotyp, boite = boite_plt1d)
711	IF debug_w THEN print, 'size(fld)', size(fld)
712	IF debug_w THEN print, 'boite_plt1d', boite_plt1d
713	IF cmd.trend GT 0 THEN BEGIN
714	fld = trends(fld, cmd.trend, datyp.hotyp)
715	add_txt = trd_txt
716	ENDIF ELSE print, 'You need to have a trend to make 1D netcdf output'
717	END
718	ELSE: BEGIN ; $$$$$$$$$$$$$$$$$$$$$$$$$$ make plot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
719
720	CASE plttyp OF
721	'plt': BEGIN
722	;
723	; Map plot
724	; --------
725	;
726	print, ' '+cmd.var+' data to plot min and max : ', $
727	min(fld(where (fld NE -valmask)), /nan), max(fld(where (fld NE valmask)), /nan)
728	print, ''
729
730	map_cmd = ''
731
732	niveau = 1
733
734	IF vecplot EQ 0 THEN BEGIN
735	; eddy field ? (if so, remove zonal mean)
736	IF strpos(cmd.plt, '@z') GT 1 THEN BEGIN
737	fldzm = moyenne(fld, 'x', boite=[20,380,-90,90], /NAN)
738	fldzm = replicate(1, (size(fld))(1))#fldzm
739	fld = fld-fldzm
740	filleg = ' Eddy '
741	ENDIF ELSE BEGIN
742	filleg = ''
743	ENDELSE
744
745	IF hotyp EQ 'xyt' THEN date_txt = 'monthly'
746
747	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
748	;; For 2D plots with overlays (correlation and proba for instance)
749	;; Scratch the titles and define contours style for the second plot
750	IF nover GT 1 AND iover GT 1 THEN BEGIN
751	contour_style = ',c_thick=2,c_linestyle=2,c_charthick=2'
752	CASE title_type OF
753	'T': pltcmd = pltcmd+',title='''''+contour_style
754	'S': pltcmd = pltcmd+',subtitle='''''+contour_style
755	'TS': pltcmd = pltcmd+',title='''',subtitle='''''+contour_style
756	'off': pltcmd = pltcmd+contour_style
757	END
758	ENDIF
759	ENDIF ELSE BEGIN
760
761	fld = {a:fld, g:strmid(cmd.grid, 0, 1)}
762	fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)}
763	fldv = {data1: fld, data2: fld2}
764	vargrid = strmid(cmd.grid, 2, 1)
765
766	pltcmd = 'ajoutvect,fldv'+vect_sample+com_strplt+strcont+map_cmd
767
768	ENDELSE
769
770	IF debug_w THEN print, strcompress(pltcmd)
771	printf, nulhis, strcompress(pltcmd)
772	res = execute(pltcmd(0))
773
774	END
775	'pltz': BEGIN
776	;
777	; Vertical section/mean plot
778	; --------------------------
779	;
780	g = gphit
781	t = tmask
782
783	; if already 1D, reform fld
784	sz = size(fld)
785	IF sz[0] EQ 2 THEN BEGIN
786	IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2])
787	ENDIF
788	IF sz[0] EQ 3 THEN BEGIN
789	IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3])
790	IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3])
791	IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2])
792	ENDIF
793
794	; mask fld
795	mask_z, fld, cmd, boite_pltz, dimplot, legz
796	printf, nulhis, ' boite_pltz = ', boite_pltz
797
798	print, ' '+cmd.var+' data to plot min and max : ', $
799	min(fld(where (fld NE valmask)), /NAN), max(fld(where (fld NE valmask)), /NAN)
800	print, ''
801	IF vecplot EQ 0 THEN BEGIN
802
803	IF cmd.var EQ 'vozonbsf' THEN xindstr = ', /xindex' ELSE xindstr = ''
804
805	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
806
807
808	ENDIF ELSE BEGIN
809
810	fld = {a:fld, g:strmid(cmd.grid, 0, 1)}
811	fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)}
812	fldv = {data1: fld, data2: fld2}
813	vargrid = strmid(cmd.grid, 2, 1)
814
815	pltcmd = 'ajoutvectz,fldv'+com_strplt+strcont+',type='''+strmid(cmd.plt, 0, 2)+''', boite=boite_pltz'
816
817	ENDELSE
818
819	IF debug_w THEN print, pltcmd
820	printf, nulhis, strcompress(pltcmd)
821	res = execute(pltcmd(0))
822
823
824	; overlay bowl if sig_bowl=1
825	IF sig_bowl EQ 1 THEN begin
826	; define line color, thickness and type
827	iover = 2
828	overc = overlay_type(iover, dimplot)
829	; type of latitude axis
830	IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $
831	sin = ',/sin' ELSE sin = ''
832	plt1dtyp = strmid(pltztyp, 0, 1)
833	noerase = 1
834	com_strplt = ', NOERASE = '+string(noerase)
835	pltcmd = 'plt1d,field.bowl,/'+plt1dtyp+', boite=boite_pltz'+overc+sin+com_strplt
836	IF debug_w THEN print, pltcmd
837
838	res = execute(pltcmd(0))
839	ENDIF
840
841	gphit = g
842	tmask = t
843
844	END
845	'plt1d': BEGIN
846	;
847	; 1D, x, y, z plot
848	; ----------------
849	;
850	g = gphit
851	t = tmask
852	; if already 1D, reform fld
853	sz = size(fld)
854	IF sz[0] EQ 2 THEN BEGIN
855	IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2])
856	ENDIF
857	IF sz[0] EQ 3 THEN BEGIN
858	IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3])
859	IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3])
860	IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2])
861	ENDIF
862	; mask fld
863	mask_z, fld, cmd, boite_plt1d, dimplot, legz
864
865	IF debug_w THEN print, boite_plt1d
866	niveau = 1
867	; define line color, thickness and type
868	overc = overlay_type(iover, dimplot)
869	; type of latitude axis
870	IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $
871	sin = ',/sin' ELSE sin = ''
872
873	print, ' '+cmd.var+' data to plot min and max : ', $
874	min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask)))
875	print, ''
876	take_log = ',take_log=0'
877	IF cmd.var EQ 'vosigvol' AND n_elements(str_sep(cmd.exp,'-')) EQ 1 THEN BEGIN ; don't take log if difference plot
878	take_log = ',take_log=1'
879	ENDIF
880	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
881	printf, nulhis, 'boite_plt1d=',boite_plt1d
882	printf, nulhis, strcompress(pltcmd)
883	IF debug_w THEN print, pltcmd
884	res = execute(pltcmd(0))
885
886	@legend_overlay
887
888	gphit = g
889	tmask = t
890
891	END
892	'pltt': BEGIN
893	;
894	; hovmoeller
895	; -----------
896	;
897	g = gphit
898	t = tmask
899
900	; mask fld
901	mask_z, fld, cmd, boite_pltt, dimplot, legz
902	; define line color, thickness and type
903	overc = overlay_type(iover, dimplot)
904	; define additional text for pltt
905	@add_txt_pltt
906	; repeat for seasonal mean
907	CASE cmd.timave OF
908	'1mm': rep_txt = ',repeat_c=nb_cycles'
909	ELSE: rep_txt = ''
910	ENDCASE
911
912	print, ' '+cmd.var+' data to plot min and max : ', $
913	min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask)))
914	print, ''
915	vardate = 'toto'
916	grille, mask, glam, gphi, gdep, nx, ny,nz
917	IF st_rms EQ 1 THEN BEGIN ; time series of rms deviation on a sigma surface
918	@densit_pltmap_plot
919	ENDIF ELSE BEGIN
920	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'
921	ENDELSE
922
923	printf, nulhis, 'boite_pltt=',boite_pltt
924	printf, nulhis, strcompress(pltcmd)
925	IF debug_w THEN print, pltcmd
926	res = execute(pltcmd(0))
927	IF iover EQ 4 THEN print, 'There might be a pb with the legends !'
928	; legend if plot=t
929	IF hotyp EQ 't' THEN BEGIN
930	; positions of the legend depend on nb_cycles
931	; it is different from 1 only with 1mm case
932	IF cmd.timave NE '1mm' THEN nb_cycles = 1
933	@legend_overlay
934	ENDIF
935	gphit = g
936	tmask = t
937	END
938
939	'yfx': BEGIN
940	;
941	; Scatter plot y=f(x)
942	; -------------------
943	;
944	t = tmask
945	boxx = ''
946	boxy = ''
947	add_txt = ''
948	; arrange/average data if needed
949	; space and time plots fld
950	IF strpos(cmd.plt, 'xy_') NE -1 THEN BEGIN
951	mask_z, fld, cmd, boite_plt, dimplot, boxx
952	boxx = ' ['+boxx+']'
953	print, ' Averaging (space) '+cmd.var+' in'+boxx
954	fld = checkfield(fld, 'plt', type = 'xy', boite = boite_plt)
955	ENDIF
956	IF datyp.hotyp NE '-' THEN BEGIN
957	mask_z, fld, cmd, boite_plt1d, dimplot, boxx
958	boxx = ' ['+boxx+']'
959	vargrid = vargrid1
960	IF debug_w THEN print, ' domdef and vargrid : ', boite_plt1d[0:3], vargrid
961	domdef, boite_plt1d[0:3]
962	print, ' Averaging (time serie plot) cmd.var '+cmd.var+' in'+boxx
963	IF debug_w THEN print, 'fld size', size(fld)
964	fld = checkfield(fld, 'plt1d', type = datyp.hotyp, boite = boite_plt1d)
965	IF cmd.trend GT 0 THEN BEGIN
966	fld_flag = 1
967	fld = trends(fld, cmd.trend, datyp.hotyp)
968	add_txt = trd_txt
969	ENDIF
970	ENDIF
971	; space and time plots fld2
972	IF sw_diffg EQ 1 THEN BEGIN
973	jptb = jpt
974	IF debug_w THEN print, ' calling def_grid, cmd2 dans plt_map '
975	def_grid, cmd2
976	jpt = jptb
977	ENDIF
978	IF strpos(cmd2.plt, 'xy_') NE -1 THEN BEGIN
979	mask_z, fld2, cmd2, boite_plt2, dimplot, boxy
980	boxy = ' ['+boxy+']'
981	print, ' Averaging (space) '+cmd2.var+' in'+boxy
982	fld2 = checkfield(fld2, 'plt', type = 'xy', boite = boite_plt2)
983	ENDIF
984	IF datyp2.hotyp NE '-' THEN BEGIN
985	run_stddev = 0
986	IF strpos(cmd2.plt, '@r') GE 1 THEN BEGIN
987	run_stddev = strmid(cmd2.plt, strpos(cmd2.plt, '@r')+2, strlen(cmd2.plt))
988	cmd2.plt = strmid(cmd2.plt, 0, strpos(cmd2.plt, '@r'))
989	ENDIF
990	mask_z, fld2, cmd2, boite_plt1d2, dimplot, boxy
991	boxy = ' ['+boxy+']'
992	vargrid = vargrid2
993	IF debug_w THEN print, ' domdef and vargrid : ', boite_plt1d2[0:3], vargrid
994	domdef, boite_plt1d2[0:3]
995	print, ' Averaging (time serie plot) cmd2.var '+cmd2.var+' in'+boxy
996	IF debug_w THEN print, ' fld2 size', size(fld2)
997	fld2 = checkfield(fld2, 'plt1d', type = datyp2.hotyp, direc = 'xy', boite = boite_plt1d2, /timearray)
998	IF cmd2.trend GT 0 THEN BEGIN
999	fld_flag = 2
1000	fld2 = trends(fld2, cmd2.trend, datyp2.hotyp)
1001	add_txt = trd_txt
1002	ENDIF
1003	ENDIF
1004	IF sw_diffg EQ 1 THEN BEGIN
1005	sw_diffg = 0
1006	jptb = jpt
1007	IF debug_w THEN print, ' calling def_grid, cmd dans plt_map '
1008	def_grid, cmd
1009	jpt = jptb
1010	ENDIF
1011	coeff = linfit(fld2, fld, CHISQ = linerr, PROB = proberr, SIGMA = sigmaerr)
1012	print, ' Linfit coef (a+bx) = ', coeff(0), coeff(1), ' errbar = ',sigmaerr, ' (proba = ',proberr, ')'
1013
1014	; plot domain
1015	boxyfx = def_box(cmd.plt, dimplot, legz, time_stride)
1016	; define line color and type
1017	overc = overlay_type(iover, dimplot)
1018	vardate = date_txt
1019	varname = varlegend
1020	varunit = ' '+cmd.var+boxx+' =f('+cmd.var2+boxy+') '+add_txt
1021	pltcmdstd = 'pltsc,fld,fld2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+overc+com_strplt
1022
1023	IF hotyp EQ 't' THEN BEGIN ; time scatter plot
1024	IF mean_sc_only EQ 0 OR mean_sc_only EQ 4 THEN BEGIN
1025	; number of time colors
1026	IF strpos(symbol_families, 'x') NE -1 THEN BEGIN
1027	coding = str_sep(symbol_families, 'x')
1028	ncolors = long(coding(0))
1029	ntimes = long(coding(1))
1030	ENDIF ELSE BEGIN
1031	ncolors = long(symbol_families)
1032	ntimes = 1
1033	ENDELSE
1034	IF ncolors GE 2 THEN BEGIN ; time color coding
1035	ic = 0
1036	jpto = jpt
1037	jpt = jpt/ncolors
1038	lincoef = fltarr(ncolors)
1039	linerro = fltarr(ncolors)
1040	linprob = fltarr(ncolors)
1041	WHILE ic LE ncolors-1 DO BEGIN
1042	idx0 = (floor(findgen(jpto/(ncolorsntimes)))ncolorsntimes)+ic(ntimes)
1043	idx = idx0
1044	jl = 1
1045	WHILE jl LE ntimes-1 DO BEGIN
1046	idx = [idx, idx0+jl]
1047	jl = jl+1
1048	ENDWHILE
1049	fldp = fld(idx)
1050	fldp2 = fld2(idx)
1051	coeff = linfit(fldp2(where (fldp2() GT 0.)), fldp(where (fldp2() GT 0.)), CHISQ = errlin, PROB = prblin)
1052	linerro(ic) = errlin
1053	linprob(ic) = prblin
1054	lincoef(ic) = coeff(1)
1055	; print, ' Period '+strtrim(string(ic), 2)+' Linear fit slope, chisq, prob =', coeff(1)1000., errlin1000., prblin
1056	IF mean_sc_only EQ 0 THEN BEGIN
1057	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)
1058	printf, nulhis, strcompress(pltcmd)
1059	IF debug_w THEN print, ' ',pltcmd
1060	res = execute(pltcmd)
1061	ENDIF
1062	IF mean_sc_only EQ 4 AND ic EQ ncolors-1 THEN BEGIN
1063	print, ' Linfit slope + error (June-Nov)=', mean(lincoef(5:10)), mean(linerro(5:10))
1064	print, ' Linfit slope + error (Dec-May)=', mean([lincoef(0:4), lincoef(11)]), mean([linerro(0:4), linerro(11)])
1065	vardate = 'toto'
1066	varname = 'Coupling strength '+cmd.date1+' '+cmd.spec
1067	jpt = ncolors
1068	time = lindgen(12)301+ $
1069	julday(1,1,01, ndayspm = calendar_type)
1070	hotyp = 't'
1071	minmax = ',0,25'
1072	pltcmd = 'pltt,lincoef*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=0,COLOR=1, thick=4, STY1D=0'
1073
1074	printf, nulhis, strcompress(pltcmd)
1075	IF debug_w THEN print, ' ', pltcmd
1076	res = execute(pltcmd)
1077	pltcmd = 'pltt,(lincoef-linerro)*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=0'
1078	printf, nulhis, strcompress(pltcmd)
1079	IF debug_w THEN print, ' ',pltcmd
1080	res = execute(pltcmd)
1081	pltcmd = 'pltt,(lincoef+linerro)*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=0'
1082	printf, nulhis, strcompress(pltcmd)
1083	IF debug_w THEN print, ' ',pltcmd
1084	res = execute(pltcmd)
1085	ENDIF
1086	ic = ic + 1
1087	ENDWHILE
1088	ENDIF ELSE BEGIN ; one color
1089	IF debug_w THEN print, ' ',pltcmdstd
1090	res = execute(pltcmdstd(0))
1091	ENDELSE
1092	ENDIF
1093	ENDIF ELSE BEGIN ; standard scatter plot
1094	IF debug_w THEN print, ' ',pltcmdstd
1095	res = execute(pltcmdstd(0))
1096	ENDELSE
1097
1098	; add mean SC in 4x3 plots
1099
1100	IF symbol_families EQ '4x3' THEN BEGIN
1101
1102	IF cmd.trend EQ '412' THEN BEGIN
1103	fldrem_t1 = fldrem_t1-mean(fldrem_t1)
1104	fldrem_t2 = fldrem_t2-mean(fldrem_t2)
1105	ENDIF ELSE BEGIN
1106	running = 12L
1107	lenght = (size(fld))[1]
1108	fldrem_t1 = fltarr(running)
1109	fldrem_t2 = fltarr(running)
1110	FOR t1 = 0, running-1 DO BEGIN
1111	fldrem_t1(t1) = mean(fld(long(findgen(lenght/running)*running+t1)))
1112	fldrem_t2(t1) = mean(fld2(long(findgen(lenght/running)*running+t1)))
1113	ENDFOR
1114	ENDELSE
1115
1116	print, ' Seasonal cycle var/stddev fld1 ', (moment(fldrem_t1))[1], sqrt((moment(fldrem_t1))[1])
1117	print, ' Seasonal cycle var/stddev fld2 ', (moment(fldrem_t2))[1], sqrt((moment(fldrem_t2))[1])
1118
1119	fldrem_t1 = [fldrem_t1, fldrem_t1(0)]
1120	fldrem_t2 = [fldrem_t2, fldrem_t2(0)]
1121	sw_ov1d = mean_sc_only
1122	IF mean_sc_only EQ 2 AND cmd.trend EQ '412' THEN BEGIN
1123	; SC of std dev
1124	stat = spectra(fld, jpt/12, 20, 15, 0)
1125	stat2 = spectra(fld2, jpt/12, 20, 15, 0)
1126	print, ' Stddev of monthly stddevs fld', sqrt((moment(stat.sc_std-(moment(stat.sc_std))[0]))[1])
1127	print, ' Stddev of monthly stddevs fld2', sqrt((moment(stat2.sc_std-(moment(stat2.sc_std))[0]))[1])
1128	stdsc = [stat.sc_std, stat.sc_std(0)]
1129	stdsc2 = [stat2.sc_std, stat2.sc_std(0)]
1130	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'
1131	printf, nulhis, strcompress(pltcmd)
1132	IF debug_w THEN print, pltcmd
1133	res = execute(pltcmd(0))
1134	pltcmd = 'xyouts,stdsc2-(maxc2)/(.3win(0)win(1)),stdsc,string(long(findgen(12))+1),charsize=1.3,alignment=0.5'
1135	printf, nulhis, strcompress(pltcmd)
1136	IF debug_w THEN print, pltcmd
1137	res = execute(pltcmd(0))
1138	ENDIF ELSE BEGIN
1139	ov1d_val = 0
1140	IF mean_sc_only EQ 0 THEN ov1d_val = 1
1141	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'
1142	printf, nulhis, strcompress(pltcmd)
1143	IF debug_w THEN print, pltcmd
1144	res = execute(pltcmd(0))
1145	IF mean_sc_only EQ 1 THEN BEGIN
1146	; add month info win(0)*win(1)
1147	pltcmd = 'xyouts,fldrem_t2-(maxc2-minc2)/(0.5win(0)win(1)),fldrem_t1,string(long(findgen(12))+1),charsize=1.3,alignment=0.5'
1148	printf, nulhis, strcompress(pltcmd)
1149	IF debug_w THEN print, pltcmd
1150	res = execute(pltcmd(0))
1151	ENDIF
1152	ENDELSE
1153	ENDIF
1154
1155	; overplot sigma contours in T-S plane
1156	IF (cmd.var EQ 'sosstsst' and cmd.var2 EQ 'sosaline') OR (cmd.var EQ 'votemper' and cmd.var2 EQ 'vosaline') THEN BEGIN
1157	IF noerase EQ 0 THEN BEGIN
1158	min_t = minc-2.0
1159	max_t = maxc+2.0
1160	min_s = minc2-2.0
1161	max_s = maxc2+2.0
1162	delta_t = 0.5
1163	delta_s = 0.05
1164	t_vec = findgen((max_t-min_t)/delta_t+1)*delta_t+min_t
1165	s_vec = findgen((max_s-min_s)/delta_s+1)*delta_s+min_s
1166	t_ = t ; since t is already defined
1167	t = t_vec
1168	FOR i = 1, (max_s-min_s)/delta_s DO t = [[t], [t_vec]]
1169	s = s_vec
1170	FOR i = 1, (max_t-min_t)/delta_t DO s = [[s], [s_vec]]
1171	s = transpose(s)
1172	sr=sqrt(abs(s))
1173	r1=((((6.536332E-9t-1.120083E-6)t+1.001685E-4)*t $
1174	-9.095290E-3)t+6.793952E-2)t+999.842594
1175	r2=(((5.3875E-9t-8.2467E-7)t+7.6438E-5)t-4.0899E-3)t+8.24493E-1
1176	r3=(-1.6546E-6t+1.0227E-4)t-5.72466E-3
1177	rhopn = ( ( 4.8314E-4s + r3sr +r2)*s +r1) -1000.0
1178	min_sig = round(min(rhopn)+0.5)-1.
1179	max_sig = round(max(rhopn)-0.5)+1.
1180	delta_sig = 1. ; define contour interval here!
1181	levels = findgen((max_sig-min_sig)/delta_sig+1)*delta_sig+min_sig
1182	labels = levels
1183	labels(*) = 1
1184	contour, rhopn, s, t, /overplot, levels = levels, c_labels = labels, c_linestyle = 0
1185	t = t_
1186	ENDIF
1187	ENDIF
1188
1189	printf, nulhis, 'boite_yfx=',boxyfx
1190	printf, nulhis, strcompress(pltcmd)
1191
1192	tmask = t
1193	; force read grid for next window
1194	cmd_prev.grid = 'toto'
1195
1196	END
1197	'spec': BEGIN
1198	;
1199	; Spectrum plot y=fft(x)
1200	; ---------------------
1201	;
1202	t = tmask
1203	; plot domain
1204	mask_z, fld, cmd, boite_pltspec, dimplot, legspec
1205	; define line color and type
1206	overc = overlay_type(iover, dimplot)
1207	vardate = date_txt
1208	varname = varlegend
1209	varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec
1210
1211	; make mean in box
1212	fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex)
1213	; apply trends
1214	IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't')
1215
1216	; make spectrum
1217	time_interval = time[1]-time[0]
1218
1219	print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360
1220	print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360
1221	print, ' spectrum window (years)/chunks = ', spec_win/360, long(time(jpt-1)-time(0)+time_interval)/spec_win
1222	print, ' '
1223
1224	; number of time windowsš
1225	nt_win = long((time(jpt-1)-time(0)+time_interval)/spec_win)
1226	idx_win_size = long(spec_win/time_interval)
1227	IF idx_win_size*nt_win NE jpt THEN BEGIN
1228	print, ' *** Warning : spec_win must divide lenght of time serie ', idx_win_size, jpt
1229	return
1230	ENDIF
1231	; mean of idx_win_size chunks
1232	spect = spectrum(fld[0:idx_win_size-1], time_interval)
1233	FOR it = 2, nt_win DO BEGIN
1234	idx1 = idx_win_size*(it-1)
1235	idx2 = idx1 + idx_win_size - 1
1236	specc = spectrum(fld[idx1:idx2], time_interval)
1237	spect = spect + specc
1238	ENDFOR
1239	spect = spect/nt_win
1240	; build new time array
1241	idx = where (spect[0, *] LE max_spec)
1242	jpt = n_elements(idx)
1243	fld = reverse(reform(spect[1, idx], jpt))
1244	time = findgen(jpt)
1245	time = reverse(reform(spect[0, idx], jpt))
1246	IF max(time) GT 20*360 THEN BEGIN
1247	time = time/360
1248	xtitle = 'Years'
1249	ENDIF ELSE IF max(time) GT 5*360 THEN BEGIN
1250	time = time/30
1251	xtitle = 'Months'
1252	ENDIF ELSE xtitle = 'Days'
1253	; time/space filter ?
1254	IF strpos(cmd.plt, '@f') GT 1 THEN BEGIN
1255	filter = long(strmid(cmd.plt, strpos(cmd.plt, '@f')+3, strlen(cmd.plt)-strpos(cmd.plt, '@f')-3))
1256	fildirec = strmid(cmd.plt, strpos(cmd.plt, '@f')+2, 1)
1257	IF fildirec EQ 't' THEN BEGIN
1258	xtitle = xtitle+' (filter '+strtrim(string(filter), 2)+' warning : discrete filter)'
1259	fld = smooth(fld, filter)
1260	fld[0:filter/2-1] = 0.
1261	fld[(size(fld))[1]-filter/2-1:(size(fld))[1]-1] = 0.
1262	ENDIF
1263	ENDIF
1264	; plot min/max
1265	min1 = min(fld)
1266	max1 = max(fld)
1267	min2 = 0
1268	max2 = max(time)
1269	!x.range = [min2-abs(max2-min2)/5.,max2+abs(max2-min2)/5.]
1270	!y.range = [min1-abs(max1-min1)/5.,max1+abs(max1-min1)/5.]
1271
1272	; plot
1273	ytitle = 'Power spectrum (window='+strtrim(string(spec_win/360), 2)+'y)'
1274	pltcmd = 'splot,time,fld,xstyle=1,ystyle=1,title=varunit,xtitle=xtitle,ytitle=ytitle'+overc+com_strplt
1275	printf, nulhis, 'boite_pltspec=',boite_pltspec
1276	printf, nulhis, strcompress(pltcmd)
1277	IF debug_w THEN print, pltcmd
1278	res = execute(pltcmd(0))
1279
1280	tmask = t
1281
1282	END
1283	'wavelet': BEGIN
1284	;
1285	; Wavelet plot
1286	; ------------
1287	;
1288	t = tmask
1289	; plot domain
1290	mask_z, fld, cmd, boite_pltspec, dimplot, legspec
1291	; define line color and type
1292	vardate = date_txt
1293	varname = varlegend + ' '+date_txt
1294	varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec
1295	; make mean in box
1296	fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex)
1297	; apply trends
1298	IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't')
1299
1300	; make spectrum
1301	time_interval = time[1]-time[0]
1302
1303	print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360
1304	print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360
1305	print, ' '
1306	; make wavelet
1307	wave = wavelet(fld,time_interval,period=period,coi=coi,/pad,signif=signif)
1308	power = abs(wave)^2
1309	nscale = n_elements(period)
1310	period = period/365 ; to have period in years
1311	; compute significance
1312	signif = rebin(transpose(signif),jpt,nscale)
1313	; plot wavelet
1314	mincmd = ','+string(min(power))
1315	maxcmd = ','+string(max(power))
1316	boite_pltspec = [0, 0, min(period), max(period)]
1317	; domdef, boite_pltspec
1318	lat1r = lat1
1319	lat2r = lat2
1320	lat1 = 0
1321	lat2 = max_spec/360
1322	key_onearth = 0
1323	pltcmd = 'plttg,power,period'+mincmd+maxcmd+intcmd+',boite=boite_pltspec,reverse_y=1,nocontour=1'+com_strplt
1324	printf, nulhis, 'boite_pltspec=',boite_pltspec
1325	printf, nulhis, strcompress(pltcmd)
1326	IF debug_w THEN print, pltcmd
1327	res = execute(pltcmd(0))
1328	contour,abs(wave)^2/signif,time,period, /overplot,level=1.0,c_annot='95%'
1329	plot,time,coi/365, noclip = 0, /noerase
1330	tmask = t
1331	lat1 = lat1r
1332	lat2 = lat2r
1333	key_onearth = 1
1334	END
1335	ELSE: BEGIN
1336	print, ' unknown projection plot ', cmd.plt, ' / plot type = ', plttyp
1337
1338	END
1339	ENDCASE
1340
1341	END
1342
1343	ENDCASE
1344
1345	fld_prev = cmd.var
1346	;
1347	; reset incompatible options of plt_def
1348
1349	cmd.trend = trend_typp
1350
1351	;
1352	; reset vertical grid after density bining
1353
1354	IF splot EQ 1 THEN BEGIN
1355
1356	izminmesh = izminmesh_b
1357	izmaxmesh = izmaxmesh_b
1358	jpk = jpk_b
1359	jpkglo = jpkglo_b
1360	gdept = gdept_b
1361	gdepw = gdepw_b
1362	e3t = e3t_b
1363	e3w = e3w_b
1364	tmask = tmask_b
1365
1366	ENDIF
1367	IF debug_w THEN print, ' ...EXIT plt_map'
1368
1369	END

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

Download in other formats:

Original Format