Home

Context Navigation

cross_sections.ncl @ 532

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

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

Download in other formats:

Original Format