Home

Context Navigation

cross_sections.ncl @ 513

Last change on this file since 513 was 513, checked in by heinze, 14 years ago
Using the NCL scripts by means of the shell script palmplot
File size: 62.3 KB

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

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

Download in other formats:

Original Format