Home

Context Navigation

cross_sections.ncl @ 981

Last change on this file since 981 was 957, checked in by heinze, 12 years ago
png output adjusted for NCL versions higher than 5.2.1
File size: 77.9 KB

Line
1	load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
2	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
3	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
4	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
5	load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
6
7	;***************************************************
8	; load .ncl.config or .ncl.config.default
9	;***************************************************
10
11	function which_script()
12	local script
13	begin
14	script="cross_section"
15	return(script)
16	end
17
18	if (isfilepresent("$PALM_BIN/../../.ncl.config")) then
19	loadscript("$PALM_BIN/../../.ncl.config")
20	else
21	if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then
22	loadscript( "$PALM_BIN/NCL/.ncl.config.default")
23	else
24	palm_bin_path = getenv("PALM_BIN")
25	print(" ")
26	print("Neither the personal configuration file '.ncl.config' exists in")
27	print("~/palm/current_version")
28	print("nor the default configuration file '.ncl.config.default' "+\
29	"exists in")
30	print(palm_bin_path + "/NCL")
31	print(" ")
32	exit
33	end if
34	end if
35
36	begin
37
38	;***************************************************
39	; Retrieving the NCL version used
40	;***************************************************
41
42	ncl_version_ch = systemfunc("ncl -V")
43	ncl_version = stringtofloat(ncl_version_ch)
44
45	; ***************************************************
46	; Retrieving the double quote character
47	; ***************************************************
48
49	dq=str_get_dq()
50
51	; ***************************************************
52	; set default parameter values and strings if not assigned
53	; in prompt or parameter list
54	; ***************************************************
55
56	; Selection of type of cross section is done in palmplot
57	if (.not. isvar("xyc"))then
58	xyc = 0
59	end if
60	if (.not. isvar("xzc"))then
61	xzc = 0
62	end if
63	if (.not. isvar("yzc"))then
64	yzc = 0
65	end if
66
67	if (file_1 .EQ. "File in") then
68	print(" ")
69	print("Declare input file 'file_1=' in '.ncl.config' or prompt")
70	print(" ")
71	exit
72	else
73	file_in = file_1
74	end if
75
76	if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND. \
77	format_out .NE. "eps" .AND. format_out .NE. "ps" .AND. \
78	format_out .NE. "epsi" .AND. format_out .NE. "ncgm" .AND. \
79	format_out .NE. "png") then
80	print(" ")
81	print("'format_out = "+format_out+"' is invalid and set to'x11'")
82	print(" ")
83	format_out="x11"
84	end if
85
86	if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then
87	print(" ")
88	print("Output of png files not available")
89	print("png output is avaiable with NCL version 5.2.0 and higher ")
90	print("NCL version used: " + ncl_version_ch)
91	print(" ")
92	exit
93	end if
94
95	if (sort .NE. "layer" .AND. sort .NE. "time") then
96	print(" ")
97	print("'sort'= "+sort+" is invalid and set to 'layer'")
98	print(" ")
99	sort = "layer"
100	end if
101
102	if (mode .NE. "Fill" .AND. mode .NE. "Line" .AND. mode .NE. "Both") then
103	print(" ")
104	print("'mode'= "+mode+" is invalid and set to 'Fill'")
105	print(" ")
106	mode = "Fill"
107	end if
108
109	if (fill_mode .NE. "AreaFill" .AND. fill_mode .NE. "RasterFill" .AND.\
110	fill_mode .NE. "CellFill") then
111	print(" ")
112	print("'fill_mode'= "+fill_mode+" is invalid and set to 'AreaFill'")
113	print(" ")
114	fill_mode = "AreaFill"
115	end if
116
117	if (shape .NE. 0 .AND. shape .NE. 1) then
118	print(" ")
119	print("'shape'= "+shape+" is invalid and set to 1")
120	print(" ")
121	shape = 1
122	end if
123
124	if (xyc .NE. 0 .AND. xyc .NE. 1)then
125	print(" ")
126	print("'xyc'= "+xyc+" is invalid and set to 0")
127	print(" ")
128	xyc = 0
129	end if
130
131	if (xzc .NE. 0 .AND. xzc .NE. 1)then
132	print(" ")
133	print("'xzc'= "+xzc+" is invalid and set to 0")
134	print(" ")
135	xyc = 0
136	end if
137
138	if (yzc .NE. 0 .AND. yzc .NE. 1)then
139	print(" ")
140	print("'yzc'= "+yzc+" is invalid and set to 0")
141	print(" ")
142	xyc = 0
143	end if
144
145	if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
146	print(" ")
147	print("Select one crossection (xyc=1 or xzc=1 or yzc=1)")
148	print(" ")
149	exit
150	end if
151	if (xyc .EQ. 1 ) then
152	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
153	print(" ")
154	print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
155	print(" ")
156	exit
157	end if
158	end if
159	if (xzc .EQ. 1) then
160	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
161	print(" ")
162	print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
163	print(" ")
164	exit
165	end if
166	end if
167	if (yzc .EQ. 1) then
168	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
169	print(" ")
170	print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
171	print(" ")
172	exit
173	end if
174	end if
175
176	if (vector .NE. 0 .AND. vector .NE. 1) then
177	print(" ")
178	print("Set 'vector' to 0 or 1")
179	print(" ")
180	exit
181	end if
182
183	if (norm_x .EQ. 0) then
184	print(" ")
185	print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
186	print(" ")
187	norm_x = 1.0
188	end if
189	if (norm_y .EQ. 0) then
190	print(" ")
191	print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
192	print(" ")
193	norm_y = 1.0
194	end if
195	if (norm_z .EQ. 0) then
196	print(" ")
197	print("Normalising with 0 is not allowed, 'norm_z' is set to 1.0")
198	print(" ")
199	norm_z = 1.0
200	end if
201
202	; ***************************************************
203	; open input file
204	; ***************************************************
205
206	file_in_1 = False
207	if (isStrSubset(file_in, ".nc"))then
208	start_f = -2
209	end_f = -2
210	file_in_1 = True
211	end if
212
213	if (start_f .EQ. -1)then
214	print(" ")
215	print("'start_f' must be one of the cyclic numbers (at least 0) "+\
216	"of your input file(s)")
217	print(" ")
218	exit
219	end if
220	if (end_f .EQ. -1)then
221	print(" ")
222	print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
223	"your input file(s)")
224	print(" ")
225	exit
226	end if
227
228	files=new(end_f-start_f+1,string)
229
230	if (file_in_1)then
231	if (isfilepresent(file_in))then
232	files(0)=file_in
233	else
234	print(" ")
235	print("1st input file: '"+file_in+"' does not exist")
236	print(" ")
237	exit
238	end if
239	else
240	if (start_f .EQ. 0)then
241	if (isfilepresent(file_in+".nc"))then
242	files(0)=file_in+".nc"
243	do i=1,end_f
244	if (isfilepresent(file_in+"."+i+".nc"))then
245	files(i)=file_in+"."+i+".nc"
246	else
247	print(" ")
248	print("Input file: '"+file_in+"."+i+".nc' does not exist")
249	print(" ")
250	exit
251	end if
252	end do
253	else
254	print(" ")
255	print("Input file: '"+file_in+".nc' does not exist")
256	print(" ")
257	exit
258	end if
259	else
260	do i=start_f,end_f
261	if (isfilepresent(file_in+"."+i+".nc"))then
262	files(i-start_f)=file_in+"."+i+".nc"
263	else
264	print(" ")
265	print("Input file: '"+file_in+"."+i+".nc' does not exist")
266	print(" ")
267	exit
268	end if
269	end do
270	end if
271	end if
272
273	f=addfiles(files,"r")
274	f_att=addfile(files(0),"r")
275	ListSetType(f,"cat")
276
277	vNam = getfilevarnames(f_att)
278	vType = getfilevartypes(f_att,vNam)
279
280	if ((all(vType .eq. "double"))) then ;distinction if data is double or float
281	check_vType = True
282	else
283	check_vType = False
284	end if
285
286	print(" ")
287	print("Variables in input file:")
288	print("- "+ vNam)
289	print(" ")
290	dim = dimsizes(vNam)
291
292	; ***************************************************
293	; check for kind of input file
294	; ***************************************************
295
296	check_3d = True
297	do varn=0,dim-1
298	checkxy = isStrSubset( vNam(varn),"xy")
299	checkxz = isStrSubset( vNam(varn),"xz")
300	checkyz = isStrSubset( vNam(varn),"yz")
301	if (checkxy .OR. checkxz .OR. checkyz) then
302	check_3d = False
303	break
304	end if
305	end do
306	if (.not. check_3d) then
307	if (xyc .EQ. 1)then
308	if (.not. checkxy) then
309	print(" ")
310	print("Your input file doesn't have values for xy cross-sections")
311	if (checkxz)then
312	print("Select another input file or xz cross-sections")
313	else
314	print("Select another input file or yz cross-sections")
315	end if
316	print(" ")
317	exit
318	else
319	print(" ")
320	print("Your input file contains xy data")
321	print(" ")
322	end if
323	end if
324	if (xzc .EQ. 1)then
325	if (.not. checkxz) then
326	print(" ")
327	print("Your input file doesn't have values for xz cross-sections")
328	if (checkxy)then
329	print("Select another input file or xy cross-sections")
330	else
331	print("Select another input file or yz cross-sections")
332	end if
333	print(" ")
334	exit
335	else
336	print(" ")
337	print("Your input file contains xz data")
338	print(" ")
339	end if
340	end if
341	if (yzc .EQ. 1)then
342	if (.not. checkyz) then
343	print(" ")
344	print("Your input file doesn't have values for yz cross-sections")
345	if (checkxy)then
346	print("Select another input file or xy cross-sections")
347	else
348	print("Select another input file or xz cross-sections")
349	end if
350	print(" ")
351	exit
352	else
353	print(" ")
354	print("Your input file contains yz data")
355	print(" ")
356	end if
357	end if
358	else
359	print(" ")
360	print("Your input file contains 3d or other data")
361	print(" ")
362	end if
363
364	if (dim .EQ. 0) then
365	print(" ")
366	print("There is no data on file")
367	print(" ")
368	exit
369	end if
370
371	; ***************************************************
372	; set up recourses
373	; ***************************************************
374
375	cs_res = True
376	vecres = True
377	cs_res@gsnYAxisIrregular2Linear = True
378	vecres@gsnYAxisIrregular2Linear = True
379
380	cs_res@gsnDraw = False
381	cs_res@gsnFrame = False
382	cs_res@gsnMaximize = True
383
384	cs_res@tmXBLabelFontHeightF = font_size
385	cs_res@tmYLLabelFontHeightF = font_size
386	cs_res@tiXAxisFontHeightF = font_size
387	cs_res@tiYAxisFontHeightF = font_size
388	cs_res@tmXBMode ="Automatic"
389	cs_res@tmYLMode ="Automatic"
390
391	cs_res@cnLevelSelectionMode = "ManualLevels"
392	cs_res@lbLabelFontHeightF = font_size_legend
393	cs_res@lbLabelStride = legend_label_stride
394
395
396	cs_resP = True
397	cs_resP@txString = f_att@title
398	cs_resP@txFuncCode = "~"
399	cs_resP@txFontHeightF = 0.015
400
401	if ( mode .eq. "Fill" ) then
402	cs_res@cnFillOn = True
403	cs_res@gsnSpreadColors = True
404	cs_res@cnFillMode = fill_mode
405	cs_res@cnLinesOn = False
406	cs_res@cnLineLabelsOn = False
407	end if
408
409	if ( mode .eq. "Both" ) then
410	cs_res@cnFillOn = True
411	cs_res@gsnSpreadColors = True
412	cs_res@cnFillMode = fill_mode
413	cs_res@cnLinesOn = True
414	cs_res@cnLineLabelsOn = True
415	end if
416
417	; *********************************************
418	; set up of start_time_step and end_time_step
419	; *********************************************
420
421	t_all = f[:]->time
422	nt = dimsizes(t_all)
423
424	if (nt .EQ. 1)then
425	delta_t=t_all(nt-1)/nt
426	else
427	delta_t_array = new(nt-1,double)
428
429	do i=0,nt-2
430	delta_t_array(i) = t_all(i+1)-t_all(i)
431	end do
432
433	delta_t = min(delta_t_array)
434	delete(delta_t_array)
435	end if
436
437
438	; ****************************************************
439	; start of time step and different types of mistakes that could be done
440	; ****************************************************
441
442	if (start_time_step .EQ. -1.) then
443	start_time_step=t_all(0)/3600
444	else
445	if (start_time_step .GT. t_all(nt-1)/3600)
446	print(" ")
447	print("'start_time_step' = "+ start_time_step +"h is greater than "+\
448	"last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
449	print(" ")
450	print("Select another 'start_time_step'")
451	print(" ")
452	exit
453	end if
454	if (start_time_step .LT. t_all(0)/3600)
455	print(" ")
456	print("'start_time_step' = "+ start_time_step +"h is lower than "+\
457	"first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
458	print(" ")
459	print("Select another 'start_time_step'")
460	print(" ")
461	exit
462	end if
463	end if
464	do i=0,nt-1
465	if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
466	start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
467	st=i
468	break
469	else
470	st=0
471	end if
472	end do
473
474
475	; ****************************************************
476	; end of time step and different types of mistakes that could be done
477	; ****************************************************
478
479	if (end_time_step .EQ. -1.) then
480	end_time_step = t_all(nt-1)/3600
481	else
482	if (end_time_step .LT. t_all(0)/3600)
483	print(" ")
484	print("'end_time_step' = "+end_time_step+ "h is lower than "+\
485	"first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
486	print(" ")
487	print("Select another 'end_time_step'")
488	print(" ")
489	exit
490	end if
491	if (end_time_step .GT. t_all(nt-1)/3600)
492	print(" ")
493	print("'end_time_step' = "+ end_time_step +"h is greater than "+\
494	"last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
495	print(" ")
496	print("Select another 'end_time_step'")
497	print(" ")
498	exit
499	end if
500	if (end_time_step .LT. start_time_step/3600)
501	print(" ")
502	print("'end_time_step' = "+ end_time_step +"h is lower than "+\
503	"'start_time_step' = "+start_time_step+"h")
504	print(" ")
505	print("Select another 'start_time_step' or 'end_time_step'")
506	print(" ")
507	exit
508	end if
509	end if
510	do i=0,nt-1
511	if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
512	end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
513	et=i
514	break
515	else
516	et=0
517	end if
518	end do
519
520	delete(start_time_step)
521	start_time_step=round(st,3)
522	delete(end_time_step)
523	end_time_step=round(et,3)
524
525	print(" ")
526	print("Output of time steps from "+t_all(start_time_step)/3600+" h = "+\
527	t_all(start_time_step)+" s => index = "+start_time_step)
528	print(" till "+t_all(end_time_step)/3600+" h = " +\
529	t_all(end_time_step)+" s => index = "+end_time_step)
530	print(" ")
531
532	no_time=(end_time_step-start_time_step)+1
533
534	; ****************************************************
535	; set up variables for vector plot if required
536	; ****************************************************
537
538	if (vector .EQ. 1) then
539	if (vec1 .EQ. "component1") then
540	print(" ")
541	print("Indicate Vector 1 ('vec1') for Vector-Plot or "+\
542	"set 'vector' to 0")
543	print(" ")
544	exit
545	end if
546	if (vec2 .EQ. "component2") then
547	print(" ")
548	print("Indicate Vector 2 ('vec2') for Vector-Plot or "+\
549	"set 'vector' to 0")
550	print(" ")
551	exit
552	end if
553	end if
554
555	check_vec1 = False
556	check_vec2 = False
557	check_vecp = False
558
559	; ****************************************************
560	; get data for plots
561	; ****************************************************
562
563	if (xyc .EQ. 1) then
564	do varn=0,dim-1
565	if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x")then
566	x_d = f_att->$vNam(varn)$
567	xdim = dimsizes(x_d)-1
568	delta_x = x_d(1)-x_d(0)
569	break
570	end if
571	end do
572	do varn=0,dim-1
573	if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then
574	y_d = f_att->$vNam(varn)$
575	ydim = dimsizes(y_d)-1
576	delta_y = y_d(1)-y_d(0)
577	break
578	end if
579	end do
580	do varn=0,dim-1
581	if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then
582	z_d = f_att->$vNam(varn)$
583	zdim = dimsizes(z_d)-1
584	delta_z = 0
585	break
586	else
587	if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then
588	z_d = f_att->$vNam(varn)$
589	zdim = dimsizes(z_d)-1
590	delta_z = -1.d
591	break
592	else
593	zdim = 0
594	delta_z = -1.d
595	end if
596	end if
597	if (vNam(varn) .eq. "zu1_xy") then
598	zu1 = f_att->$vNam(varn)$
599	end if
600	end do
601	end if
602	if (xzc .EQ. 1) then
603	do varn=0,dim-1
604	if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x") then
605	x_d = f_att->$vNam(varn)$
606	xdim = dimsizes(x_d)-1
607	delta_x = x_d(1)-x_d(0)
608	break
609	end if
610	end do
611	do varn=0,dim-1
612	if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. \
613	vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
614	z_d = f_att->$vNam(varn)$
615	zdim = dimsizes(z_d)-1
616	delta_z = z_d(1)-z_d(0)
617	break
618	end if
619	end do
620	do varn=0,dim-1
621	if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
622	y_d = f_att->$vNam(varn)$
623	ydim = dimsizes(y_d)-1
624	delta_y = y_d(1)-y_d(0)
625	break
626	else
627	if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
628	y_d = f_att->$vNam(varn)$
629	ydim = dimsizes(y_d)-1
630	delta_y = -1.d
631	break
632	else
633	ydim = 0
634	delta_y = -1.d
635	end if
636	end if
637	end do
638	end if
639	if (yzc .EQ. 1) then
640	do varn=0,dim-1
641	if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
642	y_d = f_att->$vNam(varn)$
643	ydim = dimsizes(y_d)-1
644	delta_y = y_d(1)-y_d(0)
645	break
646	end if
647	end do
648	do varn=0,dim-1
649	if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. \
650	vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
651	z_d = f_att->$vNam(varn)$
652	zdim = dimsizes(z_d)-1
653	delta_z = z_d(1)-z_d(0)
654	break
655	end if
656	end do
657	do varn=0,dim-1
658	if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
659	x_d = f_att->$vNam(varn)$
660	xdim = dimsizes(x_d)-1
661	delta_x = x_d(1)-x_d(0)
662	break
663	else
664	if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
665	x_d = f_att->$vNam(varn)$
666	xdim = dimsizes(x_d)-1
667	delta_x = -1.d
668	break
669	else
670	xdim = 0
671	delta_x = -1.d
672	end if
673	end if
674	end do
675	end if
676
677	; ****************************************************
678	; set up ranges of x-, y- and z-coordinates
679	; ****************************************************
680
681	if (xs .EQ. -1.d) then
682	xs = x_d(0)
683	if (delta_x .EQ. -1) then
684	print(" ")
685	print("You cannot choose a start value for x, "+\
686	"there are preset layers for x")
687	print(" ")
688	xstart=0
689	end if
690	else
691	if (delta_x .EQ. -1) then
692	print(" ")
693	print("You cannot choose a start value for x, "+\
694	"there are preset layers for x")
695	print(" ")
696	xstart=0
697	else
698	if (xs .LT. 0-delta_x/2) then
699	print(" ")
700	print("range start for x-coordinate = "+\
701	xs+"m is lower than first value = "+(0-delta_x/2)+"m")
702	print(" ")
703	exit
704	end if
705	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
706	if (xs .GE. xdim*delta_x) then
707	print(" ")
708	print("range start for x-coordinate = "+\
709	xs+"m is equal or greater than last value = "+\
710	xdim*delta_x+"m")
711	print(" ")
712	exit
713	end if
714	else
715	if (xs .GT. xdim*delta_x) then
716	print(" ")
717	print("range start for x-coordinate = "+\
718	xs+"m is greater than last value = "+\
719	xdim*delta_x+"m")
720	print(" ")
721	exit
722	end if
723	end if
724	end if
725	end if
726
727	if (ys .EQ. -1.d) then
728	ys = y_d(0)
729	if (delta_y .EQ. -1) then
730	print(" ")
731	print("You cannot choose a start value for y, "+\
732	"there are preset layers for y")
733	print(" ")
734	ystart=0
735	end if
736	else
737	if (delta_y .EQ. -1) then
738	print(" ")
739	print("You cannot choose a start value for y, "+\
740	"there are preset layers for y")
741	print(" ")
742	ystart=0
743	else
744	if (ys .LT. 0-delta_y/2) then
745	print(" ")
746	print("range start for y-coordinate = "+\
747	ys+"m is lower than first value = "+0-delta_y/2+"m")
748	print(" ")
749	exit
750	end if
751	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
752	if (ys .GE. ydim*delta_y) then
753	print(" ")
754	print("range start for y-coordinate = "+\
755	ys+"m is equal or greater than last value = "+\
756	ydim*delta_y+"m")
757	print(" ")
758	exit
759	end if
760	else
761	if (ys .GT. ydim*delta_y) then
762	print(" ")
763	print("range start for y-coordinate = "+\
764	ys+"m is greater than last value = "+ydim*delta_y+"m")
765	print(" ")
766	exit
767	end if
768	end if
769	end if
770	end if
771
772	if (zs .EQ. -1) then
773	zs = 0
774	if (delta_z .EQ. -1) then
775	print(" ")
776	print("You cannot choose a start value for z, "+\
777	"there are preset layers for z")
778	print(" ")
779	end if
780	else
781	if (delta_z .EQ. -1) then
782	print(" ")
783	print("You cannot choose a start value for z, "+\
784	"there are preset layers for z")
785	print(" ")
786	zs = 0
787	else
788	if (zs .LT. 0) then
789	print(" ")
790	print("range start for z-coordinate = "+\
791	zs+" is lower than first gridpoint = 0")
792	print(" ")
793	exit
794	end if
795	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
796	if (zs .GE. zdim) then
797	print(" ")
798	print("range start for z-coordinate = "+\
799	zs+" is equal or greater than last gridpoint = "+zdim)
800	print(" ")
801	exit
802	end if
803	else
804	if (zs .GT. zdim) then
805	print(" ")
806	print("range start for z-coordinate = "+\
807	zs+" is greater than last gridpoint = "+zdim)
808	print(" ")
809	exit
810	end if
811	end if
812	end if
813	end if
814
815	if (xe .EQ. -1) then
816	xe = x_d(xdim)
817	if (delta_x .EQ. -1) then
818	print(" ")
819	print("You cannot choose an end value for x, "+\
820	"there are preset layers for x")
821	print(" ")
822	xend=xdim
823	end if
824	else
825	if (delta_x .EQ. -1) then
826	print(" ")
827	print("You cannot choose an end value for x, "+\
828	"there are preset layers for x")
829	print(" ")
830	xend=xdim
831	else
832	if (xe .GT. xdim*delta_x) then
833	print(" ")
834	print("range end for x-coordinate = "+\
835	xe+"m is greater than last value = "+xdim*delta_x+"m")
836	print(" ")
837	exit
838	end if
839	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
840	if (xe .LE. 0-delta_x/2)
841	print(" ")
842	print("range end for x-coordinate = "+\
843	xe+"m is equal or lower than first value = "+\
844	(0-delta_x/2)+"m")
845	print(" ")
846	exit
847	end if
848	if (xe .LE. xs) then
849	print(" ")
850	print("range end for x-coordinate = "+\
851	xe+"m is equal or lower than start range = "+xs+"m")
852	print(" ")
853	exit
854	end if
855	if (xe .EQ. xs+1) then
856	print(" ")
857	print("range end for x-coordinate = "+\
858	xe+"m must be at least two more gridpoints "+\
859	"greater than start range = "+xs+"m")
860	print(" ")
861	exit
862	end if
863	else
864	if (xe .LT. 0-delta_x/2)
865	print(" ")
866	print("range end for x-coordinate = "+\
867	xe+"m is lower than first value = "+(0-delta_x/2)+"m")
868	print(" ")
869	exit
870	end if
871	if (xe .LT. xs) then
872	print(" ")
873	print("range end for x-coordinate = "+\
874	xe+"m is lower than start range = "+xs+"m")
875	print(" ")
876	exit
877	end if
878	end if
879	end if
880	end if
881
882	if (ye .EQ. -1) then
883	ye = y_d(ydim)
884	if (delta_y .EQ. -1) then
885	print(" ")
886	print("You cannot choose an end value for y, "+\
887	"there are preset layers for y")
888	print(" ")
889	yend=ydim
890	end if
891	else
892	if (delta_y .EQ. -1) then
893	print(" ")
894	print("You cannot choose an end value for y, "+\
895	"there are preset layers for y")
896	print(" ")
897	yend=ydim
898	else
899	if (ye .GT. ydim*delta_y) then
900	print(" ")
901	print("range end for y-coordinate = "+\
902	ye+"m is greater than last value = "+ydim*delta_y+"m")
903	print(" ")
904	exit
905	end if
906	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
907	if (ye .LE. 0-delta_y/2)
908	print(" ")
909	print("range end for y-coordinate = "+\
910	ye+"m is equal or lower than first value = "+\
911	(0-delta_y/2)+"m")
912	print(" ")
913	exit
914	end if
915	if (ye .LE. ys) then
916	print(" ")
917	print("range end for y-coordinate = "+\
918	ye+"m is equal or lower than start range = "+ys+"m")
919	print(" ")
920	exit
921	end if
922	if (ye .EQ. ys+1) then
923	print(" ")
924	print("range end for y-coordinate = "+\
925	ye+"m must be at least two more gridpoints greater "+\
926	"than start range = "+ys+"m")
927	print(" ")
928	exit
929	end if
930	else
931	if (ye .LT. 0-delta_y/2)
932	print(" ")
933	print("range end for y-coordinate = "+\
934	ye+"m is lower than first value = "+(0-delta_y/2)+"m")
935	print(" ")
936	exit
937	end if
938	if (ye .LT. ys) then
939	print(" ")
940	print("range end for y-coordinate = "+\
941	ye+"m is lower than start range = "+ys+"m")
942	print(" ")
943	exit
944	end if
945	end if
946	end if
947	end if
948
949	if (ze .EQ. -1) then
950	ze = zdim
951	if (delta_z .EQ. -1) then
952	print(" ")
953	print("You cannot choose an end value for z, "+\
954	"there are preset layers for z")
955	print(" ")
956	ze = zdim
957	end if
958	else
959	if (delta_z .EQ. -1) then
960	print(" ")
961	print("You cannot choose an end value for z, "+\
962	"there are preset layers for z")
963	print(" ")
964	ze = zdim
965	else
966	if (ze .GT. zdim) then
967	print(" ")
968	print("range end for z-coordinate = "+\
969	ze+" is greater than last gridpoint = "+zdim)
970	print(" ")
971	exit
972	end if
973	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
974	if (ze .LE. 0)
975	print(" ")
976	print("range end for z-coordinate = "+\
977	ze+" is equal or lower than first gridpoint = 0")
978	print(" ")
979	exit
980	end if
981	if (ze .LE. zs) then
982	print(" ")
983	print("range end for z-coordinate = "+\
984	ze+" is equal or lower than start range = "+zs)
985	print(" ")
986	exit
987	end if
988	if (ze .EQ. zs+1) then
989	print(" ")
990	print("range end for z-coordinate = "+\
991	ze+" must be at least two more gridpoints "+\
992	"greater than start range = "+zs)
993	print(" ")
994	exit
995	end if
996	else
997	if (ze .LT. 0)
998	print(" ")
999	print("range end for z-coordinate = "+\
1000	ze+" is lower than first gridpoint = 0")
1001	print(" ")
1002	exit
1003	end if
1004	if (ze .LT. zs) then
1005	print(" ")
1006	print("range end for z-coordinate = "+\
1007	ze+" is lower than start range = "+zs)
1008	print(" ")
1009	exit
1010	end if
1011	end if
1012	end if
1013	end if
1014
1015	if (delta_x .NE. -1) then
1016	do i=0,xdim
1017	if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
1018	xstart=i
1019	break
1020	end if
1021	end do
1022	do i=0,xdim
1023	if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
1024	xend=i
1025	break
1026	end if
1027	end do
1028	end if
1029	if (delta_y .NE. -1) then
1030	do i=0,ydim
1031	if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
1032	ystart=i
1033	break
1034	end if
1035	end do
1036	do i=0,ydim
1037	if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
1038	yend=i
1039	break
1040	end if
1041	end do
1042	end if
1043
1044	if( shape .eq. 1 ) then
1045	if (xyc .EQ. 1)then
1046	cs_res@vpWidthF = (xe-xs)/(ye-ys)
1047	cs_res@vpHeightF = 1
1048	vecres@vpWidthF = (xe-xs)/(ye-ys)
1049	vecres@vpHeightF = 1
1050	if (xe-xs .GT. ye-ys)then
1051	if (no_rows .LE. no_columns)then
1052	cs_res@gsnPaperOrientation = "landscape"
1053	vecres@gsnPaperOrientation = "landscape"
1054	cs_res@lbOrientation = "Horizontal"
1055	else
1056	cs_res@gsnPaperOrientation = "portrait"
1057	vecres@gsnPaperOrientation = "portrait"
1058	cs_res@lbOrientation = "Vertical"
1059	end if
1060	else
1061	if (no_rows .GE. no_columns)then
1062	cs_res@gsnPaperOrientation = "portrait"
1063	vecres@gsnPaperOrientation = "portrait"
1064	cs_res@lbOrientation = "Vertical"
1065	else
1066	cs_res@gsnPaperOrientation = "landscape"
1067	vecres@gsnPaperOrientation = "landscape"
1068	cs_res@lbOrientation = "Horizontal"
1069	end if
1070	end if
1071	end if
1072	if (xzc .EQ. 1)then
1073	cs_res@vpWidthF = (xe-xs)/(delta_x*(ze-zs))
1074	cs_res@vpHeightF = 1
1075	vecres@vpWidthF = (xe-xs)/(delta_x*(ze-zs))
1076	vecres@vpHeightF = 1
1077	if (xe-xs .GT. (delta_x*(ze-zs)))then
1078	if (no_rows .LE. no_columns)then
1079	cs_res@gsnPaperOrientation = "landscape"
1080	vecres@gsnPaperOrientation = "landscape"
1081	cs_res@lbOrientation = "Horizontal"
1082	else
1083	cs_res@gsnPaperOrientation = "portrait"
1084	vecres@gsnPaperOrientation = "portrait"
1085	cs_res@lbOrientation = "Vertical"
1086	end if
1087	else
1088	if (no_rows .GE. no_columns)then
1089	cs_res@gsnPaperOrientation = "portrait"
1090	vecres@gsnPaperOrientation = "portrait"
1091	cs_res@lbOrientation = "Vertical"
1092	else
1093	cs_res@gsnPaperOrientation = "landscape"
1094	vecres@gsnPaperOrientation = "landscape"
1095	cs_res@lbOrientation = "Horizontal"
1096	end if
1097	end if
1098	end if
1099	if (yzc .EQ. 1)then
1100	cs_res@vpWidthF = (ye-ys)/(delta_y*(ze-zs))
1101	cs_res@vpHeightF = 1
1102	vecres@vpWidthF = (ye-ys)/(delta_y*(ze-zs))
1103	vecres@vpHeightF = 1
1104	if (ye-ys .GT. (delta_y*(ze-zs)))then
1105	if (no_rows .LE. no_columns)then
1106	cs_res@gsnPaperOrientation = "landscape"
1107	vecres@gsnPaperOrientation = "landscape"
1108	cs_res@lbOrientation = "Horizontal"
1109	else
1110	cs_res@gsnPaperOrientation = "portrait"
1111	vecres@gsnPaperOrientation = "portrait"
1112	cs_res@lbOrientation = "Vertical"
1113	end if
1114	else
1115	if (no_rows .GE. no_columns)then
1116	cs_res@gsnPaperOrientation = "portrait"
1117	vecres@gsnPaperOrientation = "portrait"
1118	cs_res@lbOrientation = "Vertical"
1119	else
1120	cs_res@gsnPaperOrientation = "landscape"
1121	vecres@gsnPaperOrientation = "landscape"
1122	cs_res@lbOrientation = "Horizontal"
1123	end if
1124	end if
1125	end if
1126	end if
1127
1128	delete(xs)
1129	delete(xe)
1130	delete(ys)
1131	delete(ye)
1132
1133	xs=xstart
1134	xe=xend
1135	ys=ystart
1136	ye=yend
1137
1138	if (xyc .EQ. 1) then
1139	d= (ye-ys+1)/(major_ticks_y-1)
1140	e=(xe-xs+1)/(major_ticks_x-1)
1141	array_yl =new(major_ticks_y,integer)
1142	array_yl_labels=new(major_ticks_y,double)
1143	array_minor_yl =new((major_ticks_y-1)*4,double)
1144	array_xb =new(major_ticks_x,integer)
1145	array_xb_labels=new(major_ticks_x,double)
1146	array_minor_xb =new((major_ticks_x-1)*4,double)
1147	array_yl(0)=ys
1148	array_xb(0)=xs
1149	array_yl_labels(0)=y_d(ys)
1150	array_xb_labels(0)=x_d(xs)
1151	do ar=1,major_ticks_y-1
1152	array_yl(ar)=d*(ar-1)+d-1
1153	array_yl_labels(ar) = y_d(array_yl(ar))
1154	if (ar .GT. 0)
1155	do min_ar=0,3
1156	array_minor_yl(4*(ar-1)+min_ar)= y_d(array_yl(ar-1))+\
1157	(y_d(array_yl(ar))-y_d(array_yl(ar-1)))/5*(min_ar+1)
1158	end do
1159	end if
1160	end do
1161	do br=1,major_ticks_x-1
1162	array_xb(br)=e*(br-1)+e-1
1163	array_xb_labels(br) = x_d(array_xb(br))
1164	if (br .GT. 0)
1165	do min_br=0,3
1166	array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+\
1167	(x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
1168	end do
1169	end if
1170	end do
1171	end if
1172
1173	if (xzc .EQ. 1) then
1174	d=(ze-zs+1)/(major_ticks_y-1)
1175	e=(xe-xs+1)/(major_ticks_x-1)
1176	array_yl =new(major_ticks_y,integer)
1177	array_yl_labels=new(major_ticks_y,double)
1178	array_minor_yl =new((major_ticks_y-1)*4,double)
1179	array_xb =new(major_ticks_x,integer)
1180	array_xb_labels=new(major_ticks_x,double)
1181	array_minor_xb =new((major_ticks_x-1)*4,double)
1182	array_yl(0)=zs
1183	array_xb(0)=xs
1184	array_yl_labels(0)=z_d(zs)
1185	array_xb_labels(0)=x_d(xs)
1186	do ar=1,major_ticks_y-1
1187	array_yl(ar)=d*(ar-1)+d-1
1188	array_yl_labels(ar) = z_d(array_yl(ar))
1189	if (ar .GT. 0)
1190	do min_ar=0,3
1191	array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+\
1192	(z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
1193	end do
1194	end if
1195	end do
1196	do br=1,major_ticks_x-1
1197	array_xb(br)=e*(br-1)+e-1
1198	array_xb_labels(br) = x_d(array_xb(br))
1199	if (br .GT. 0)
1200	do min_br=0,3
1201	array_minor_xb(4*(br-1)+min_br)= x_d(array_xb(br-1))+\
1202	(x_d(array_xb(br))-x_d(array_xb(br-1)))/5*(min_br+1)
1203	end do
1204	end if
1205	end do
1206	end if
1207
1208	if (yzc .EQ. 1) then
1209	d=(ze-zs+1)/(major_ticks_y-1)
1210	e=(ye-ys+1)/(major_ticks_x-1)
1211	array_yl =new(major_ticks_y,integer)
1212	array_yl_labels=new(major_ticks_y,double)
1213	array_minor_yl =new((major_ticks_y-1)*4,double)
1214	array_xb =new(major_ticks_x,integer)
1215	array_xb_labels=new(major_ticks_x,double)
1216	array_minor_xb =new((major_ticks_x-1)*4,double)
1217	array_yl(0)=zs
1218	array_xb(0)=ys
1219	array_yl_labels(0)=z_d(zs)
1220	array_xb_labels(0)=y_d(ys)
1221	do ar=1,major_ticks_y-1
1222	array_yl(ar)=d*(ar-1)+d-1
1223	array_yl_labels(ar) = z_d(array_yl(ar))
1224	if (ar .GT. 0)
1225	do min_ar=0,3
1226	array_minor_yl(4*(ar-1)+min_ar)= z_d(array_yl(ar-1))+\
1227	(z_d(array_yl(ar))-z_d(array_yl(ar-1)))/5*(min_ar+1)
1228	end do
1229	end if
1230	end do
1231	do br=1,major_ticks_x-1
1232	array_xb(br)=e*(br-1)+e-1
1233	array_xb_labels(br) = y_d(array_xb(br))
1234	if (br .GT. 0)
1235	do min_br=0,3
1236	array_minor_xb(4*(br-1)+min_br)= y_d(array_xb(br-1))+\
1237	(y_d(array_xb(br))-y_d(array_xb(br-1)))/5*(min_br+1)
1238	end do
1239	end if
1240	end do
1241	end if
1242
1243	if (axes_explicit .EQ. 1)then
1244	cs_res@tmYLMode = "Explicit"
1245	cs_res@tmXBMode = "Explicit"
1246	if (xyc .EQ. 1)then
1247	cs_res@tmYLValues = y_d(array_yl)
1248	cs_res@tmXBValues = x_d(array_xb)
1249	cs_res@tmYLLabels = array_yl_labels/norm_y
1250	cs_res@tmXBLabels = array_xb_labels/norm_x
1251	if (norm_x .NE. 1.)then
1252	cs_res@tiXAxisString = "x ("+unit_x+")"
1253	else
1254	cs_res@tiXAxisString = "x (m)"
1255	end if
1256	if (norm_y .NE. 1.)then
1257	cs_res@tiYAxisString = "y ("+unit_y+")"
1258	else
1259	cs_res@tiYAxisString = "y (m)"
1260	end if
1261	end if
1262	if (xzc .EQ. 1)then
1263	cs_res@tmYLValues = z_d(array_yl)
1264	cs_res@tmXBValues = x_d(array_xb)
1265	cs_res@tmYLLabels = array_yl_labels/norm_z
1266	cs_res@tmXBLabels = array_xb_labels/norm_x
1267	if (norm_x .NE. 1.)then
1268	cs_res@tiXAxisString = "x ("+unit_x+")"
1269	else
1270	cs_res@tiXAxisString = "x (m)"
1271	end if
1272	if (norm_z .NE. 1.)then
1273	cs_res@tiYAxisString = "z ("+unit_z+")"
1274	else
1275	cs_res@tiYAxisString = "z (m)"
1276	end if
1277	end if
1278	if (yzc .EQ. 1)then
1279	cs_res@tmYLValues = z_d(array_yl)
1280	cs_res@tmXBValues = y_d(array_xb)
1281	cs_res@tmYLLabels = array_yl_labels/norm_z
1282	cs_res@tmXBLabels = array_xb_labels/norm_y
1283	if (norm_y .NE. 1.)then
1284	cs_res@tiXAxisString = "y ("+unit_y+")"
1285	else
1286	cs_res@tiXAxisString = "y (m)"
1287	end if
1288	if (norm_z .NE. 1.)then
1289	cs_res@tiYAxisString = "z ("+unit_z+")"
1290	else
1291	cs_res@tiYAxisString = "z (m)"
1292	end if
1293	end if
1294	cs_res@tmYLMinorValues = array_minor_yl
1295	cs_res@tmXBMinorValues = array_minor_xb
1296	else
1297	if (axes_explicit .EQ. 0)then
1298	cs_res@tmYLMinorPerMajor = 4
1299	cs_res@tmXBMinorPerMajor = 4
1300	else
1301	print(" ")
1302	print("'axes_explicit' is invalid and set to 0")
1303	print(" ")
1304	axes_explicit = 0
1305	cs_res@tmYLMinorPerMajor = 4
1306	cs_res@tmXBMinorPerMajor = 4
1307	end if
1308	if (xyc .EQ. 1)then
1309	cs_res@tiXAxisString = "x (m)"
1310	cs_res@tiYAxisString = "y (m)"
1311	end if
1312	if (xzc .EQ. 1)then
1313	cs_res@tiXAxisString = "x (m)"
1314	cs_res@tiYAxisString = "z (m)"
1315	end if
1316	if (yzc .EQ. 1)then
1317	cs_res@tiXAxisString = "y (m)"
1318	cs_res@tiYAxisString = "z (m)"
1319	end if
1320	end if
1321
1322	if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1323	if (xe .EQ. xs+1) then
1324	print(" ")
1325	print("range end for x-coordinate="+xe*delta_x+\
1326	"m must be at least two")
1327	print("more gridpoints(="+2*delta_x+"m) greater than start range="+\
1328	xs*delta_x+"m)")
1329	print(" ")
1330	exit
1331	end if
1332	end if
1333	if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1334	if (ye .EQ. ys+1) then
1335	print(" ")
1336	print("range end for y-coordinate="+ye*delta_y+\
1337	"m must be at least two")
1338	print("more gridpoints(="+2*delta_y+"m greater than start range="+\
1339	ys*delta_y+"m)")
1340	print(" ")
1341	exit
1342	end if
1343	end if
1344	if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1345	if (ze .EQ. zs+1) then
1346	print(" ")
1347	print("range end for z-coordinate="+ze+" must be at least two")
1348	print("more gridpoints greater than start range="+zs+")")
1349	print(" ")
1350	exit
1351	end if
1352	end if
1353
1354	if (xyc .EQ. 1) then
1355	no_layer = (ze-zs)+1
1356	if (check_vType) then
1357	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1358	else
1359	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1360	end if
1361	end if
1362	if (xzc .EQ. 1) then
1363	no_layer = (ye-ys)+1
1364	if (check_vType) then
1365	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1366	else
1367	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1368	end if
1369	end if
1370	if (yzc .EQ. 1) then
1371	no_layer = (xe-xs)+1
1372	if (check_vType) then
1373	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),double)
1374	else
1375	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1376	end if
1377	end if
1378
1379	if (check_vType) then
1380	MinVal = new(dim,double)
1381	MaxVal = new(dim,double)
1382	else
1383	MinVal = new(dim,float)
1384	MaxVal = new(dim,float)
1385	end if
1386	unit = new(dim,string)
1387
1388	; ****************************************************
1389	; define inner and outer loops depending on "sort"
1390	; ****************************************************
1391
1392	if ( xyc .eq. 1 ) then
1393	lays = zs
1394	laye = ze
1395	end if
1396	if ( xzc .eq. 1 ) then
1397	lays = ys
1398	laye = ye
1399	end if
1400	if ( yzc .eq. 1) then
1401	lays = xs
1402	laye = xe
1403	end if
1404
1405	if ( sort .eq. "time" ) then
1406	lis = start_time_step
1407	lie = end_time_step
1408	los = lays
1409	loe = laye
1410	else
1411	lis = lays
1412	lie = laye
1413	los = start_time_step
1414	loe = end_time_step
1415	end if
1416
1417	check = True
1418	v1 = 0
1419	v2 = 0
1420	no_var = 0
1421	n = 0
1422	no_zu1 = 0
1423	no_topo= 0
1424
1425	;****************************************************
1426	; Preparation of vector plots
1427	; ****************************************************
1428
1429	do varn=dim-1,0,1
1430
1431	if (vector .EQ. 1) then
1432	check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1433	if (check_vec1) then
1434	temp = f[:]->$vNam(varn)$
1435	data_att = f_att->$vNam(varn)$
1436	vect1=temp(:,zs:ze,ys:ye,xs:xe)
1437	delete(temp)
1438	end if
1439
1440	check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1441	if (check_vec2) then
1442	temp = f[:]->$vNam(varn)$
1443	data_att = f_att->$vNam(varn)$
1444	vect2=temp(:,zs:ze,ys:ye,xs:xe)
1445	delete(temp)
1446	end if
1447
1448	if (","+vNam(varn)+"," .EQ. vec1) then
1449	v1=v1+1
1450	end if
1451	if (","+vNam(varn)+"," .EQ. vec2) then
1452	v2=v2+1
1453	end if
1454	end if
1455
1456	end do
1457
1458
1459
1460	do varn=dim-1,0,1
1461
1462	if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. \
1463	vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. \
1464	vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. \
1465	vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. \
1466	vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. \
1467	vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. \
1468	vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. \
1469	vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1470	vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. \
1471	vNam(varn) .eq. "ind_x_yz") then
1472	check = False
1473	else
1474	check = True
1475	end if
1476
1477	if (var .NE. "all") then
1478	check = isStrSubset( var,","+vNam(varn)+"," )
1479	end if
1480
1481	if(check) then
1482
1483	no_var=no_var+1
1484
1485	if (xyc .EQ. 1) then
1486	temp = f[:]->$vNam(varn)$
1487	if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1488	dummy=0
1489	else
1490	data_att = f_att->$vNam(varn)$
1491	end if
1492	if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy" \
1493	.or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \
1494	.or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy" \
1495	.or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy" \
1496	.or. vNam(varn) .eq. "lwp_xy" .or. vNam(varn) .eq. "pra_xy" \
1497	.or. vNam(varn) .eq. "prr_xy" .or. vNam(varn) .eq. "qsws_xy"\
1498	.or. vNam(varn) .eq. "shf_xy" .or. vNam(varn) .eq. "t_xy" \
1499	.or. vNam(varn) .eq. "u_xy" .or. vNam(varn) .eq. "z0_xy") then
1500	;these variables depend on zu1_xy and that's why they have
1501	;only one z-layer
1502	data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe)
1503	no_zu1=no_zu1+1
1504	else
1505	if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1506	;these variables depend on x and y
1507	data(varn,0,0,:,:)=doubletofloat(temp(ys:ye,xs:xe))
1508	no_topo=no_topo+1
1509	else
1510	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1511	end if
1512	end if
1513	delete(temp)
1514	end if
1515	if ( xzc .eq. 1 ) then
1516	temp = f[:]->$vNam(varn)$
1517	data_att = f_att->$vNam(varn)$
1518	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1519	delete(temp)
1520	end if
1521	if ( yzc .eq. 1) then
1522	temp = f[:]->$vNam(varn)$
1523	data_att = f_att->$vNam(varn)$
1524	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1525	delete(temp)
1526	end if
1527
1528	data!0 = "var"
1529	data!1 = "t"
1530	data!2 = "z"
1531	data!3 = "y"
1532	data!4 = "x"
1533
1534	MinVal(varn) = min(data(varn,start_time_step:end_time_step,\
1535	0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1536	MaxVal(varn) = max(data(varn,start_time_step:end_time_step,\
1537	0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1538
1539	if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1540	unit(varn) = "meters"
1541	else
1542	unit(varn) = data_att@units
1543	delete(data_att)
1544	end if
1545
1546	end if
1547
1548	end do
1549
1550	if (no_var .EQ. 0) then
1551	print(" ")
1552	print("The variables var='"+var+"' do not exist on your input file;")
1553	print("be sure to have one comma before and after each variable")
1554	print(" ")
1555	exit
1556	end if
1557
1558	var_input=new(no_var,string)
1559	numb=0
1560	no_var=0
1561
1562	do varn=dim-1,0,1
1563
1564	if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. \
1565	vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. \
1566	vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. \
1567	vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. \
1568	vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. \
1569	vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. \
1570	vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. \
1571	vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1572	vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. \
1573	vNam(varn) .eq. "ind_x_yz") then
1574	check = False
1575	else
1576	check = True
1577	end if
1578
1579	if (var .NE. "all") then
1580	check = isStrSubset( var,","+vNam(varn)+"," )
1581	end if
1582
1583	if(check) then
1584	var_input(numb)=vNam(varn)
1585	numb=numb+1
1586	end if
1587
1588	if (check) then
1589	no_var = no_var+1
1590	end if
1591
1592	end do
1593
1594	if (no_var .EQ. 0) then
1595	print(" ")
1596	print("The variables var='"+var+"' do not exist on your input file;")
1597	print("be sure to have one comma before and after each variable")
1598	print(" ")
1599	exit
1600	end if
1601
1602	if (vector .EQ. 1) then
1603	if (v1 .EQ. 0)then
1604	print(" ")
1605	print("Component 1 for the vector-plot ('vec1') must be one "+\
1606	"of the variables on the input file:")
1607	print("- "+var_input)
1608	print("be sure to have one comma before and after the variable")
1609	print(" ")
1610	exit
1611	end if
1612
1613	if (v2 .EQ. 0)then
1614	print(" ")
1615	print("Component 2 for the vector-plot ('vec2') must be one "+\
1616	"of the variables on the input file:")
1617	print("- "+var_input)
1618	print("be sure to have one comma before and after the variable")
1619	print(" ")
1620	exit
1621	end if
1622	end if
1623
1624	; ***************************************************
1625	; open workstation(s)
1626	; ***************************************************
1627
1628	if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
1629	format_out@wkPaperSize = "A4"
1630	end if
1631	if ( format_out .EQ. "png" ) then
1632	format_out@wkWidth = 1000
1633	format_out@wkHeight = 1000
1634	end if
1635
1636	wks_ps = gsn_open_wks(format_out,file_out)
1637	gsn_define_colormap(wks_ps,"rainbow+white")
1638
1639	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1640	plot=new((/no_time*no_layer/),graphic)
1641	else
1642	plot=new((/no_timeno_layer(no_var-no_topo) + no_topo - \
1643	no_time(no_layer-1)no_zu1/),graphic)
1644	end if
1645	dim_plot=dimsizes(plot)
1646
1647	page = 0
1648	n=0
1649	print(" ")
1650	print("Plots sorted by '"+sort+"'")
1651	print(" ")
1652
1653	; ***************************************************
1654	; create plots
1655	; ***************************************************
1656
1657
1658	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1659	do lo = los, loe
1660	do li = lis, lie
1661	if (xyc .EQ. 1)then
1662	if (sort .EQ. "time")then
1663	if(z_d(lo) .eq. -1.d) then
1664	level = "z-average"
1665	else
1666	level = "z=" + z_d(lo) + "m"
1667	end if
1668	else
1669	if(z_d(li) .eq. -1.d) then
1670	level = "z-average"
1671	else
1672	level = "z=" + z_d(li) + "m"
1673	end if
1674	end if
1675	end if
1676	if (xzc .EQ. 1)then
1677	if (sort .EQ. "time")then
1678	if(y_d(lo) .eq. -1.d) then
1679	level = "y-average"
1680	else
1681	level = "y=" + y_d(lo) + "m"
1682	end if
1683	else
1684	if(y_d(li) .eq. -1.d) then
1685	level = "y-average"
1686	else
1687	level = "y=" + y_d(li) + "m"
1688	end if
1689	end if
1690	end if
1691	if (yzc .EQ. 1)then
1692	if (sort .EQ. "time")then
1693	if(x_d(lo) .eq. -1.d) then
1694	level = "x-average"
1695	else
1696	level = "x=" + x_d(lo) + "m"
1697	end if
1698	else
1699	if(x_d(li) .eq. -1.d) then
1700	level = "x-average"
1701	else
1702	level = "x=" + x_d(li) + "m"
1703	end if
1704	end if
1705	end if
1706	vecres = True ; vector only resources
1707	vecres@gsnDraw = False ; don't draw
1708	vecres@gsnFrame = False ; don't advance frame
1709	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1710	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1711	vecres@vcRefLengthF = 0.05 ; define length of
1712	; vec ref
1713	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
1714	vecres@tmXBLabelFontHeightF = font_size
1715	vecres@tmYLLabelFontHeightF = font_size
1716	vecres@tiXAxisFontHeightF = font_size
1717	vecres@tiYAxisFontHeightF = font_size
1718	if (sort .EQ. "time")then
1719	vecres@gsnLeftString = "t=" + \
1720	decimalPlaces(t_all(li)/3600,2,True) +"h "+level
1721	else
1722	vecres@gsnLeftString = "t=" + \
1723	decimalPlaces(t_all(lo)/3600,2,True) +"h "+level
1724	end if
1725	if (xyc .EQ. 1) then
1726	vecres@tiXAxisString = "x (m)"
1727	vecres@tiYAxisString = "y (m)"
1728	if (sort .EQ. "time")then
1729	plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),\
1730	vect2(li,lo-los,:,:),vecres)
1731	else
1732	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1733	vect2(lo,li-lis,:,:),vecres)
1734	end if
1735	end if
1736	if (xzc .EQ. 1) then
1737	vecres@tiXAxisString = "x (m)"
1738	vecres@tiYAxisString = "z (m)"
1739	if (sort .EQ. "time")then
1740	plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
1741	vect2(li,:,lo-los,:),vecres)
1742	else
1743	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
1744	vect2(lo,:,li-lis,:),vecres)
1745	end if
1746	end if
1747	if (yzc .EQ. 1) then
1748	vecres@tiXAxisString = "y (m)"
1749	vecres@tiYAxisString = "z (m)"
1750	if (sort .EQ. "time")then
1751	plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
1752	vect2(li,:,:,lo-los),vecres)
1753	else
1754	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
1755	vect2(lo,:,:,li-lis),vecres)
1756	end if
1757	end if
1758	n=n+1
1759	end do
1760	end do
1761	end if
1762
1763	do varn=dim-1,0,1
1764
1765	if (vector .EQ. 1 ) then
1766	check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1767	end if
1768
1769	if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. \
1770	vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. \
1771	vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. \
1772	vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. \
1773	vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. \
1774	vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. \
1775	vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. \
1776	vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1777	vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. \
1778	vNam(varn) .eq. "ind_x_yz") then
1779	check = False
1780	else
1781	check = True
1782	end if
1783
1784	if (var .NE. "all") then
1785	check = isStrSubset( var,","+vNam(varn)+"," )
1786	end if
1787
1788	if(check) then
1789
1790	space=(MaxVal(varn)-MinVal(varn))/24
1791
1792	cs_res@cnMinLevelValF = MinVal(varn)
1793	cs_res@cnMaxLevelValF = MaxVal(varn)
1794
1795	cs_res@cnLevelSpacingF = space
1796
1797	; ****************************************************
1798	; loops over time and layer
1799	; ****************************************************
1800
1801
1802	do lo = los, loe
1803	do li = lis, lie
1804
1805	; ****************************************************
1806	; xy cross section
1807	; ****************************************************
1808
1809	if ( xyc .eq. 1 ) then
1810
1811	cs_res@gsnLeftString = "Plot of "+vNam(varn)
1812	cs_res@gsnRightString = unit(varn)
1813
1814	if ( sort .eq. "time" ) then
1815	if ( z_d(lo) .eq. -1)then
1816	if (delta_z .EQ. -1) then
1817	level = "-average"
1818	else
1819	level = "=" + z_d(lo) + "m"
1820	end if
1821	else
1822	level = "=" + z_d(lo) + "m"
1823	end if
1824
1825	if(vNam(varn) .eq. "lwps_xy" .or. \
1826	vNam(varn) .eq. "pras_xy" .or. \
1827	vNam(varn) .eq. "prrs_xy" .or. \
1828	vNam(varn) .eq. "qsws_xy" .or. \
1829	vNam(varn) .eq. "shfs_xy" .or. \
1830	vNam(varn) .eq. "ts_xy" .or. \
1831	vNam(varn) .eq. "us_xy" .or. \
1832	vNam(varn) .eq. "z0s_xy" .or. \
1833	vNam(varn) .eq. "lwp*_xy" .or. \
1834	vNam(varn) .eq. "pra*_xy" .or. \
1835	vNam(varn) .eq. "prr*_xy" .or. \
1836	vNam(varn) .eq. "qsws*_xy" .or. \
1837	vNam(varn) .eq. "shf*_xy" .or. \
1838	vNam(varn) .eq. "t*_xy" .or. \
1839	vNam(varn) .eq. "u*_xy" .or. \
1840	vNam(varn) .eq. "z0*_xy") then
1841	loe = 0
1842	level = "=" + zu1(0) + "m"
1843	else
1844	loe = laye
1845	end if
1846
1847	if(vNam(varn) .eq. "zwwi" .or. \
1848	vNam(varn) .eq. "zusi") then
1849	loe = 0
1850	los = 0
1851	lie = 0
1852	lis = 0
1853	level = ""
1854	else
1855	lis = start_time_step
1856	lie = end_time_step
1857	los = lays
1858	loe = laye
1859	end if
1860
1861	if(vNam(varn) .eq. "zwwi" .or. \
1862	vNam(varn) .eq. "zusi") then
1863	cs_res@gsnCenterString = ""
1864	else
1865	cs_res@gsnCenterString = "t=" + \
1866	decimalPlaces(t_all(li)/3600,2,True) +"h z"+level
1867	end if
1868
1869	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1870	;nothing to be done
1871	else
1872	plot(n) = gsn_csm_contour(wks_ps,\
1873	data(varn,li,lo-los,:,:),cs_res)
1874	end if
1875
1876
1877	if (vector .EQ. 1 .AND. check_vecp) then
1878	vecres = True ; vector only
1879	; resources
1880	vecres@gsnDraw = False ; don't draw
1881	vecres@gsnFrame = False ; don't
1882	; advance frame
1883	vecres@vcGlyphStyle = "CurlyVector" ; curly
1884	; vectors
1885	vecres@vcRefMagnitudeF = ref_mag ; define
1886	; vector ref
1887	; mag
1888	vecres@vcRefLengthF = 0.05 ; define
1889	; length of
1890	; vec ref
1891	vecres@gsnRightString = " "
1892	vecres@gsnLeftString = " "
1893	vecres@tiXAxisString = " "
1894	plot_vec=gsn_csm_vector(wks_ps,\
1895	vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1896	overlay(plot(n), plot_vec)
1897	end if
1898	end if
1899
1900	if ( sort .eq. "layer" ) then
1901	if (z_d(li) .eq. -1) then
1902	if (delta_z .EQ. -1) then
1903	level = "-average"
1904	else
1905	level = "=" + z_d(li) + "m"
1906	end if
1907	else
1908	level = "=" + z_d(li) + "m"
1909	end if
1910
1911	if(vNam(varn) .eq. "lwps_xy" .or. \
1912	vNam(varn) .eq. "pras_xy" .or. \
1913	vNam(varn) .eq. "prrs_xy" .or. \
1914	vNam(varn) .eq. "qsws_xy" .or. \
1915	vNam(varn) .eq. "shfs_xy" .or. \
1916	vNam(varn) .eq. "ts_xy" .or. \
1917	vNam(varn) .eq. "us_xy" .or. \
1918	vNam(varn) .eq. "z0s_xy" .or. \
1919	vNam(varn) .eq. "lwp*_xy" .or. \
1920	vNam(varn) .eq. "pra*_xy" .or. \
1921	vNam(varn) .eq. "prr*_xy" .or. \
1922	vNam(varn) .eq. "qsws*_xy" .or. \
1923	vNam(varn) .eq. "shf*_xy" .or. \
1924	vNam(varn) .eq. "t*_xy" .or. \
1925	vNam(varn) .eq. "u*_xy" .or. \
1926	vNam(varn) .eq. "z0*_xy") then
1927	lie = 0
1928	level = "=" + zu1(0) + "m"
1929	else
1930	lie = laye
1931	end if
1932
1933	if(vNam(varn) .eq. "zwwi" .or. \
1934	vNam(varn) .eq. "zusi") then
1935	lie = 0
1936	lis = 0
1937	loe = 0
1938	los = 0
1939	level = ""
1940	else
1941	lis = lays
1942	lie = laye
1943	los = start_time_step
1944	loe = end_time_step
1945	end if
1946
1947	if(vNam(varn) .eq. "zwwi" .or. \
1948	vNam(varn) .eq. "zusi") then
1949	cs_res@gsnCenterString = ""
1950	else
1951	cs_res@gsnCenterString = "t=" + \
1952	decimalPlaces(t_all(lo)/3600,2,True) +"h z"+level
1953	end if
1954
1955	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1956	;nothing to be done
1957	else
1958	plot(n) = gsn_csm_contour(wks_ps,\
1959	data(varn,lo,li-lis,:,:),cs_res)
1960	end if
1961
1962	if (vector .EQ. 1 .AND. check_vecp) then
1963	vecres = True ; vector only
1964	; resources
1965	vecres@gsnDraw = False ; don't draw
1966	vecres@gsnFrame = False ; don't
1967	; advance frame
1968	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1969	vecres@vcRefMagnitudeF = ref_mag ; define
1970	; vector ref
1971	; mag
1972	vecres@vcRefLengthF = 0.05 ; define
1973	; length of
1974	; vec ref
1975	vecres@gsnRightString = " " ; turn off
1976	; right string
1977	vecres@gsnLeftString = " " ; turn off
1978	; left string
1979	vecres@tiXAxisString = " "
1980	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1981	vect2(lo,li-lis,:,:),vecres)
1982	overlay(plot(n), plot_vec)
1983	end if
1984	end if
1985	end if
1986
1987	; ****************************************************
1988	; xz cross section
1989	; ****************************************************
1990
1991	if ( xzc .eq. 1 ) then
1992
1993	cs_res@gsnLeftString = "Plot of "+vNam(varn)
1994
1995	if ( sort .eq. "time" ) then
1996	if ( y_d(lo) .eq. -1 ) then
1997	level = "-average"
1998	else
1999	level = "=" + y_d(lo) + "m"
2000	end if
2001	cs_res@gsnCenterString = "t=" + \
2002	decimalPlaces(t_all(li)/3600,2,True) + "h y"+ level
2003
2004	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2005	;nothing to be done
2006	else
2007	plot(n) = gsn_csm_contour(wks_ps,\
2008	data(varn,li,:,lo-los,:),cs_res)
2009	end if
2010
2011	if (vector .EQ. 1 .AND. check_vecp) then
2012	vecres = True ; vector only
2013	; resources
2014	vecres@gsnDraw = False ; don't draw
2015	vecres@gsnFrame = False ; don't
2016	; advance frame
2017	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
2018	vecres@vcRefMagnitudeF = ref_mag ; define
2019	; vector ref
2020	; mag
2021	vecres@vcRefLengthF = 0.05 ; define
2022	; length of
2023	; vec ref
2024	vecres@gsnRightString = " " ; turn off
2025	; right string
2026	vecres@gsnLeftString = " " ; turn off
2027	; left string
2028	vecres@tiXAxisString = " "
2029	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
2030	vect2(li,:,lo-los,:),vecres)
2031	overlay(plot(n), plot_vec)
2032	end if
2033	end if
2034
2035	if ( sort .eq. "layer" ) then
2036	if ( y_d(li) .eq. -1 ) then
2037	level = "-average"
2038	else
2039	level = "=" + y_d(li) + "m"
2040	end if
2041	cs_res@gsnCenterString = "t=" + \
2042	decimalPlaces(t_all(lo)/3600,2,True) + "h y"+ level
2043
2044	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2045	;nothing to be done
2046	else
2047	plot(n) = gsn_csm_contour(wks_ps,\
2048	data(varn,lo,:,li-lis,:),cs_res)
2049	end if
2050
2051	if (vector .EQ. 1 .AND. check_vecp) then
2052	vecres = True ; vector only
2053	; resources
2054	vecres@gsnDraw = False ; don't draw
2055	vecres@gsnFrame = False ; don't
2056	; advance frame
2057	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
2058	vecres@vcRefMagnitudeF = ref_mag ; define
2059	; vector ref
2060	; mag
2061	vecres@vcRefLengthF = 0.05 ; define
2062	; length of
2063	; vec ref
2064	vecres@gsnRightString = " " ; turn off
2065	; right string
2066	vecres@gsnLeftString = " " ; turn off
2067	; left string
2068	vecres@tiXAxisString = " "
2069	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
2070	vect2(lo,:,li-lis,:),vecres)
2071	overlay(plot(n), plot_vec)
2072	end if
2073	end if
2074	end if
2075
2076	; ****************************************************
2077	; yz cross section
2078	; ****************************************************
2079
2080	if ( yzc .eq. 1 ) then
2081
2082	cs_res@gsnLeftString = "Plot of "+vNam(varn)
2083
2084	if ( sort .eq. "time" ) then
2085	if ( x_d(lo) .eq. -1 ) then
2086	level = "-average"
2087	else
2088	level = "=" + x_d(lo) + "m"
2089	end if
2090	cs_res@gsnCenterString = "t=" + \
2091	decimalPlaces(t_all(li)/3600,2,True) + "h x"+ level
2092
2093	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2094	;nothing to be done
2095	else
2096	plot(n) = gsn_csm_contour(wks_ps,\
2097	data(varn,li,:,:,lo-los),cs_res)
2098	end if
2099
2100	if (vector .EQ. 1 .AND. check_vecp) then
2101	vecres = True ; vector only
2102	; resources
2103	vecres@gsnDraw = False ; don't draw
2104	vecres@gsnFrame = False ; don't
2105	; advance frame
2106	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
2107	vecres@vcRefMagnitudeF = ref_mag ; define
2108	; vector ref
2109	; mag
2110	vecres@vcRefLengthF = 0.05 ; define
2111	; length of
2112	; vec ref
2113	vecres@gsnRightString = " " ; turn off
2114	; right string
2115	vecres@gsnLeftString = " " ; turn off
2116	; left string
2117	vecres@tiXAxisString = " "
2118	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
2119	vect2(li,:,:,lo-los),vecres)
2120	overlay(plot(n), plot_vec)
2121	end if
2122	end if
2123
2124	if ( sort .eq. "layer" ) then
2125	if ( x_d(li) .eq. -1 ) then
2126	level = "-average"
2127	else
2128	level = "=" + x_d(li) + "m"
2129	end if
2130	cs_res@gsnCenterString = "t=" + \
2131	decimalPlaces(t_all(lo)/3600,2,True) + "h x"+ level
2132
2133	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2134	;nothing to be done
2135	else
2136	plot(n) = gsn_csm_contour(wks_ps,\
2137	data(varn,lo,:,:,li-lis),cs_res)
2138	end if
2139
2140	if (vector .EQ. 1 .AND. check_vecp)then
2141	vecres = True ; vector only
2142	;resources
2143	vecres@gsnDraw = False ; don't draw
2144	vecres@gsnFrame = False ; don't
2145	; advance frame
2146	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
2147	vecres@vcRefMagnitudeF = ref_mag ; define
2148	; vector ref
2149	; mag
2150	vecres@vcRefLengthF = 0.05 ; define
2151	; length of
2152	; vec ref
2153	vecres@gsnRightString = " " ; turn off
2154	; right string
2155	vecres@gsnLeftString = " " ; turn off
2156	; left string
2157	vecres@tiXAxisString = " "
2158	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
2159	vect2(lo,:,:,li-lis),vecres)
2160	overlay(plot(n), plot_vec)
2161	end if
2162	end if
2163	end if
2164	n=n+1
2165	end do
2166	end do
2167	end if
2168	end do
2169
2170	; ***************************************************
2171	; merge plots onto one page
2172	; ***************************************************
2173
2174	no_frames = 0
2175
2176	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2177	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
2178	no_timeno_layer(no_var+1) .gt. no_rows*no_columns) then
2179	gsn_panel(wks_ps,plot(0:(no_timeno_layer(no_var+1))-1),\
2180	(/no_var+1,no_layer*no_time/),cs_resP)
2181	print(" ")
2182	print("Outputs to .eps or .epsi have only one frame")
2183	print(" ")
2184	else
2185	do np = 0,dim_plot-1,no_rows*no_columns
2186	if ( np + no_rows*no_columns .gt. dim_plot-1) then
2187	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2188	cs_resP)
2189	no_frames = no_frames + 1
2190	else
2191	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2192	(/no_rows,no_columns/),cs_resP)
2193	no_frames = no_frames + 1
2194	end if
2195	end do
2196	end if
2197	else
2198	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
2199	.AND. dim_plot .gt. no_rows*no_columns) then
2200	gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
2201	print(" ")
2202	print("Outputs to .eps or .epsi have only one frame")
2203	print(" ")
2204	else
2205	do np = 0,dim_plot-1,no_rows*no_columns
2206	if ( np + no_rows*no_columns .gt. dim_plot-1) then
2207	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2208	cs_resP)
2209	no_frames = no_frames + 1
2210	else
2211	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2212	(/no_rows,no_columns/),cs_resP)
2213	no_frames = no_frames + 1
2214	end if
2215	end do
2216	end if
2217	end if
2218
2219	if (format_out .EQ. "png" ) then
2220	png_output = new((/no_frames/), string)
2221	j = 0
2222
2223	if (no_frames .eq. 1) then
2224	if (ncl_version .GE. 6 ) then
2225	png_output(0) = file_out+".png"
2226	else
2227	png_output(0) = file_out+".00000"+1+".png"
2228	end if
2229	;using imagemagick's convert for reducing the white
2230	;space around the plot
2231	cmd = "convert -geometry 1000x1000 -density 300 -trim " + \
2232	png_output(0) + " " + png_output(0)
2233	system(cmd)
2234	else
2235
2236	do i=0, no_frames-1
2237	j = i + 1
2238	if (j .LE. 9) then
2239	png_output(i) = file_out+".00000"+j+".png"
2240	end if
2241	if (j .GT. 9 .AND. j .LE. 99) then
2242	png_output(i) = file_out+".0000"+j+".png"
2243	end if
2244	if (j .GT. 99 .AND. j .LE. 999) then
2245	png_output(i) = file_out+".000"+j+".png"
2246	end if
2247	if (j .GT. 999) then
2248	png_output(i) = file_out+".00"+j+".png"
2249	end if
2250
2251	;using imagemagick's convert for reducing the white
2252	;space around the plot
2253	cmd = "convert -geometry 1000x1000 -density 300 -trim " + \
2254	png_output(i) + " " + png_output(i)
2255	system(cmd)
2256	end do
2257	end if
2258
2259	print(" ")
2260	print("Output to: "+ png_output)
2261	print(" ")
2262	else
2263	print(" ")
2264	print("Output to: " + file_out +"."+ format_out)
2265	print(" ")
2266	end if
2267
2268	end

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

Download in other formats:

Original Format