Home

Context Navigation

cross_sections.ncl @ 1038

Last change on this file since 1038 was 983, checked in by hoffmann, 12 years ago
Bugfix in netcdf.f90 and improvements in palmplot.
File size: 78.1 KB

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

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

Download in other formats:

Original Format