Home

Context Navigation

cross_sections.ncl @ 1673

Last change on this file since 1673 was 1626, checked in by heinze, 9 years ago
Bugfix to plot xy-cross sections in case variables contain only one level (like lwp_xy, shf_xy, ...)
File size: 78.7 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. "zwwi" .or. \
1857	vNam(varn) .eq. "zusi") then
1858	loe = 0
1859	los = 0
1860	lie = 0
1861	lis = 0
1862	level = ""
1863	end if
1864
1865	if(vNam(varn) .eq. "lwps_xy" .or. \
1866	vNam(varn) .eq. "pras_xy" .or. \
1867	vNam(varn) .eq. "prrs_xy" .or. \
1868	vNam(varn) .eq. "qsws_xy" .or. \
1869	vNam(varn) .eq. "shfs_xy" .or. \
1870	vNam(varn) .eq. "ts_xy" .or. \
1871	vNam(varn) .eq. "us_xy" .or. \
1872	vNam(varn) .eq. "z0s_xy" .or. \
1873	vNam(varn) .eq. "lwp*_xy" .or. \
1874	vNam(varn) .eq. "pra*_xy" .or. \
1875	vNam(varn) .eq. "prr*_xy" .or. \
1876	vNam(varn) .eq. "qsws*_xy" .or. \
1877	vNam(varn) .eq. "shf*_xy" .or. \
1878	vNam(varn) .eq. "t*_xy" .or. \
1879	vNam(varn) .eq. "u*_xy" .or. \
1880	vNam(varn) .eq. "z0*_xy") then
1881	loe = 0
1882	level = "=" + zu1(0) + "m"
1883	else
1884	lis = start_time_step
1885	lie = end_time_step
1886	los = lays
1887	loe = laye
1888	end if
1889
1890	if(vNam(varn) .eq. "zwwi" .or. \
1891	vNam(varn) .eq. "zusi") then
1892	cs_res@gsnCenterString = ""
1893	else
1894	cs_res@gsnCenterString = "t=" + \
1895	decimalPlaces(t_all(li)/3600,2,True) +"h z"+level
1896	end if
1897
1898	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1899	;nothing to be done
1900	else
1901	plot(n) = gsn_csm_contour(wks_ps,\
1902	data(varn,li,lo-los,:,:),cs_res)
1903	end if
1904
1905
1906	if (vector .EQ. 1 .AND. check_vecp) then
1907	vecres = True ; vector only
1908	; resources
1909	vecres@gsnDraw = False ; don't draw
1910	vecres@gsnFrame = False ; don't
1911	; advance frame
1912	vecres@vcGlyphStyle = "CurlyVector" ; curly
1913	; vectors
1914	vecres@vcRefMagnitudeF = ref_mag ; define
1915	; vector ref
1916	; mag
1917	vecres@vcRefLengthF = 0.05 ; define
1918	; length of
1919	; vec ref
1920	vecres@gsnRightString = " "
1921	vecres@gsnLeftString = " "
1922	vecres@tiXAxisString = " "
1923	plot_vec=gsn_csm_vector(wks_ps,\
1924	vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1925	overlay(plot(n), plot_vec)
1926	end if
1927	end if
1928
1929	if ( sort .eq. "layer" ) then
1930	if (z_d(li) .eq. -1) then
1931	if (delta_z .EQ. -1) then
1932	level = "-average"
1933	else
1934	level = "=" + z_d(li) + "m"
1935	end if
1936	else
1937	level = "=" + z_d(li) + "m"
1938	end if
1939
1940	if(vNam(varn) .eq. "zwwi" .or. \
1941	vNam(varn) .eq. "zusi") then
1942	lie = 0
1943	lis = 0
1944	loe = 0
1945	los = 0
1946	level = ""
1947	end if
1948
1949	if(vNam(varn) .eq. "lwps_xy" .or. \
1950	vNam(varn) .eq. "pras_xy" .or. \
1951	vNam(varn) .eq. "prrs_xy" .or. \
1952	vNam(varn) .eq. "qsws_xy" .or. \
1953	vNam(varn) .eq. "shfs_xy" .or. \
1954	vNam(varn) .eq. "ts_xy" .or. \
1955	vNam(varn) .eq. "us_xy" .or. \
1956	vNam(varn) .eq. "z0s_xy" .or. \
1957	vNam(varn) .eq. "lwp*_xy" .or. \
1958	vNam(varn) .eq. "pra*_xy" .or. \
1959	vNam(varn) .eq. "prr*_xy" .or. \
1960	vNam(varn) .eq. "qsws*_xy" .or. \
1961	vNam(varn) .eq. "shf*_xy" .or. \
1962	vNam(varn) .eq. "t*_xy" .or. \
1963	vNam(varn) .eq. "u*_xy" .or. \
1964	vNam(varn) .eq. "z0*_xy") then
1965	lie = 0
1966	level = "=" + zu1(0) + "m"
1967	else
1968	lis = lays
1969	lie = laye
1970	los = start_time_step
1971	loe = end_time_step
1972	end if
1973
1974	if(vNam(varn) .eq. "zwwi" .or. \
1975	vNam(varn) .eq. "zusi") then
1976	cs_res@gsnCenterString = ""
1977	else
1978	cs_res@gsnCenterString = "t=" + \
1979	decimalPlaces(t_all(lo)/3600,2,True) +"h z"+level
1980	end if
1981
1982	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1983	;nothing to be done
1984	else
1985	plot(n) = gsn_csm_contour(wks_ps,\
1986	data(varn,lo,li-lis,:,:),cs_res)
1987	end if
1988
1989	if (vector .EQ. 1 .AND. check_vecp) then
1990	vecres = True ; vector only
1991	; resources
1992	vecres@gsnDraw = False ; don't draw
1993	vecres@gsnFrame = False ; don't
1994	; advance frame
1995	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1996	vecres@vcRefMagnitudeF = ref_mag ; define
1997	; vector ref
1998	; mag
1999	vecres@vcRefLengthF = 0.05 ; define
2000	; length of
2001	; vec ref
2002	vecres@gsnRightString = " " ; turn off
2003	; right string
2004	vecres@gsnLeftString = " " ; turn off
2005	; left string
2006	vecres@tiXAxisString = " "
2007	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
2008	vect2(lo,li-lis,:,:),vecres)
2009	overlay(plot(n), plot_vec)
2010	end if
2011	end if
2012	end if
2013
2014	; ****************************************************
2015	; xz cross section
2016	; ****************************************************
2017
2018	if ( xzc .eq. 1 ) then
2019
2020	cs_res@gsnLeftString = "Plot of "+vNam(varn)
2021
2022	if ( sort .eq. "time" ) then
2023	if ( y_d(lo) .eq. -1 ) then
2024	level = "-average"
2025	else
2026	level = "=" + y_d(lo) + "m"
2027	end if
2028	cs_res@gsnCenterString = "t=" + \
2029	decimalPlaces(t_all(li)/3600,2,True) + "h y"+ level
2030
2031	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2032	;nothing to be done
2033	else
2034	plot(n) = gsn_csm_contour(wks_ps,\
2035	data(varn,li,:,lo-los,:),cs_res)
2036	end if
2037
2038	if (vector .EQ. 1 .AND. check_vecp) then
2039	vecres = True ; vector only
2040	; resources
2041	vecres@gsnDraw = False ; don't draw
2042	vecres@gsnFrame = False ; don't
2043	; advance frame
2044	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
2045	vecres@vcRefMagnitudeF = ref_mag ; define
2046	; vector ref
2047	; mag
2048	vecres@vcRefLengthF = 0.05 ; define
2049	; length of
2050	; vec ref
2051	vecres@gsnRightString = " " ; turn off
2052	; right string
2053	vecres@gsnLeftString = " " ; turn off
2054	; left string
2055	vecres@tiXAxisString = " "
2056	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
2057	vect2(li,:,lo-los,:),vecres)
2058	overlay(plot(n), plot_vec)
2059	end if
2060	end if
2061
2062	if ( sort .eq. "layer" ) then
2063	if ( y_d(li) .eq. -1 ) then
2064	level = "-average"
2065	else
2066	level = "=" + y_d(li) + "m"
2067	end if
2068	cs_res@gsnCenterString = "t=" + \
2069	decimalPlaces(t_all(lo)/3600,2,True) + "h y"+ level
2070
2071	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2072	;nothing to be done
2073	else
2074	plot(n) = gsn_csm_contour(wks_ps,\
2075	data(varn,lo,:,li-lis,:),cs_res)
2076	end if
2077
2078	if (vector .EQ. 1 .AND. check_vecp) then
2079	vecres = True ; vector only
2080	; resources
2081	vecres@gsnDraw = False ; don't draw
2082	vecres@gsnFrame = False ; don't
2083	; advance frame
2084	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
2085	vecres@vcRefMagnitudeF = ref_mag ; define
2086	; vector ref
2087	; mag
2088	vecres@vcRefLengthF = 0.05 ; define
2089	; length of
2090	; vec ref
2091	vecres@gsnRightString = " " ; turn off
2092	; right string
2093	vecres@gsnLeftString = " " ; turn off
2094	; left string
2095	vecres@tiXAxisString = " "
2096	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
2097	vect2(lo,:,li-lis,:),vecres)
2098	overlay(plot(n), plot_vec)
2099	end if
2100	end if
2101	end if
2102
2103	; ****************************************************
2104	; yz cross section
2105	; ****************************************************
2106
2107	if ( yzc .eq. 1 ) then
2108
2109	cs_res@gsnLeftString = "Plot of "+vNam(varn)
2110
2111	if ( sort .eq. "time" ) then
2112	if ( x_d(lo) .eq. -1 ) then
2113	level = "-average"
2114	else
2115	level = "=" + x_d(lo) + "m"
2116	end if
2117	cs_res@gsnCenterString = "t=" + \
2118	decimalPlaces(t_all(li)/3600,2,True) + "h x"+ level
2119
2120	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2121	;nothing to be done
2122	else
2123	plot(n) = gsn_csm_contour(wks_ps,\
2124	data(varn,li,:,:,lo-los),cs_res)
2125	end if
2126
2127	if (vector .EQ. 1 .AND. check_vecp) then
2128	vecres = True ; vector only
2129	; resources
2130	vecres@gsnDraw = False ; don't draw
2131	vecres@gsnFrame = False ; don't
2132	; advance frame
2133	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
2134	vecres@vcRefMagnitudeF = ref_mag ; define
2135	; vector ref
2136	; mag
2137	vecres@vcRefLengthF = 0.05 ; define
2138	; length of
2139	; vec ref
2140	vecres@gsnRightString = " " ; turn off
2141	; right string
2142	vecres@gsnLeftString = " " ; turn off
2143	; left string
2144	vecres@tiXAxisString = " "
2145	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
2146	vect2(li,:,:,lo-los),vecres)
2147	overlay(plot(n), plot_vec)
2148	end if
2149	end if
2150
2151	if ( sort .eq. "layer" ) then
2152	if ( x_d(li) .eq. -1 ) then
2153	level = "-average"
2154	else
2155	level = "=" + x_d(li) + "m"
2156	end if
2157	cs_res@gsnCenterString = "t=" + \
2158	decimalPlaces(t_all(lo)/3600,2,True) + "h x"+ level
2159
2160	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2161	;nothing to be done
2162	else
2163	plot(n) = gsn_csm_contour(wks_ps,\
2164	data(varn,lo,:,:,li-lis),cs_res)
2165	end if
2166
2167	if (vector .EQ. 1 .AND. check_vecp)then
2168	vecres = True ; vector only
2169	;resources
2170	vecres@gsnDraw = False ; don't draw
2171	vecres@gsnFrame = False ; don't
2172	; advance frame
2173	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
2174	vecres@vcRefMagnitudeF = ref_mag ; define
2175	; vector ref
2176	; mag
2177	vecres@vcRefLengthF = 0.05 ; define
2178	; length of
2179	; vec ref
2180	vecres@gsnRightString = " " ; turn off
2181	; right string
2182	vecres@gsnLeftString = " " ; turn off
2183	; left string
2184	vecres@tiXAxisString = " "
2185	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
2186	vect2(lo,:,:,li-lis),vecres)
2187	overlay(plot(n), plot_vec)
2188	end if
2189	end if
2190	end if
2191	n=n+1
2192	end do
2193	end do
2194	end if
2195	end do
2196
2197	; ***************************************************
2198	; merge plots onto one page
2199	; ***************************************************
2200
2201	no_frames = 0
2202
2203	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2204	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
2205	no_timeno_layer(no_var+1) .gt. no_rows*no_columns) then
2206	gsn_panel(wks_ps,plot(0:(no_timeno_layer(no_var+1))-1),\
2207	(/no_var+1,no_layer*no_time/),cs_resP)
2208	print(" ")
2209	print("Outputs to .eps or .epsi have only one frame")
2210	print(" ")
2211	else
2212	do np = 0,dim_plot-1,no_rows*no_columns
2213	if ( np + no_rows*no_columns .gt. dim_plot-1) then
2214	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2215	cs_resP)
2216	no_frames = no_frames + 1
2217	else
2218	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2219	(/no_rows,no_columns/),cs_resP)
2220	no_frames = no_frames + 1
2221	end if
2222	end do
2223	end if
2224	else
2225	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
2226	.AND. dim_plot .gt. no_rows*no_columns) then
2227	gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
2228	print(" ")
2229	print("Outputs to .eps or .epsi have only one frame")
2230	print(" ")
2231	else
2232	do np = 0,dim_plot-1,no_rows*no_columns
2233	if ( np + no_rows*no_columns .gt. dim_plot-1) then
2234	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2235	cs_resP)
2236	no_frames = no_frames + 1
2237	else
2238	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2239	(/no_rows,no_columns/),cs_resP)
2240	no_frames = no_frames + 1
2241	end if
2242	end do
2243	end if
2244	end if
2245
2246	if (format_out .EQ. "png" ) then
2247	png_output = new((/no_frames/), string)
2248	j = 0
2249
2250	if (no_frames .eq. 1) then
2251	if (ncl_version .GE. 6 ) then
2252	png_output(0) = file_out+".png"
2253	else
2254	png_output(0) = file_out+".00000"+1+".png"
2255	end if
2256	;using imagemagick's convert for reducing the white
2257	;space around the plot
2258	cmd = "convert -geometry 1000x1000 -density 300 -trim " + \
2259	png_output(0) + " " + png_output(0)
2260	system(cmd)
2261	else
2262
2263	do i=0, no_frames-1
2264	j = i + 1
2265	if (j .LE. 9) then
2266	png_output(i) = file_out+".00000"+j+".png"
2267	end if
2268	if (j .GT. 9 .AND. j .LE. 99) then
2269	png_output(i) = file_out+".0000"+j+".png"
2270	end if
2271	if (j .GT. 99 .AND. j .LE. 999) then
2272	png_output(i) = file_out+".000"+j+".png"
2273	end if
2274	if (j .GT. 999) then
2275	png_output(i) = file_out+".00"+j+".png"
2276	end if
2277
2278	;using imagemagick's convert for reducing the white
2279	;space around the plot
2280	cmd = "convert -geometry 1000x1000 -density 300 -trim " + \
2281	png_output(i) + " " + png_output(i)
2282	system(cmd)
2283	end do
2284	end if
2285
2286	print(" ")
2287	print("Output to: "+ png_output)
2288	print(" ")
2289	else
2290	print(" ")
2291	print("Output to: " + file_out +"."+ format_out)
2292	print(" ")
2293	end if
2294
2295	end

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

Download in other formats:

Original Format