Home

Context Navigation

cross_sections.ncl @ 1597

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

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

Download in other formats:

Original Format