Context Navigation

plt_map.pro @ 2

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

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

Download in other formats:

Original Format