Home

Context Navigation

cross_sections.ncl @ 357

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

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

Download in other formats:

Original Format