Home

Context Navigation

cross_sections.ncl @ 1380

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

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

Download in other formats:

Original Format