Home

Context Navigation

cross_sections.ncl @ 2334

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

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

Download in other formats:

Original Format