Home

Context Navigation

cross_sections.ncl @ 320

Last change on this file since 320 was 286, checked in by weinreis, 16 years ago
ncl-script cross_sections.ncl modified
File size: 59.5 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	end do
533	end if
534	if (xzc .EQ. 1) then
535	do varn=0,dim-1
536	if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x") then
537	x_d = f_att->$vNam(varn)$
538	xdim = dimsizes(x_d)-1
539	delta_x = x_d(1)-x_d(0)
540	break
541	end if
542	end do
543	do varn=0,dim-1
544	if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
545	z_d = f_att->$vNam(varn)$
546	zdim = dimsizes(z_d)-1
547	delta_z = z_d(1)-z_d(0)
548	break
549	end if
550	end do
551	do varn=0,dim-1
552	if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
553	y_d = f_att->$vNam(varn)$
554	ydim = dimsizes(y_d)-1
555	delta_y = y_d(1)-y_d(0)
556	break
557	else
558	if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
559	y_d = f_att->$vNam(varn)$
560	ydim = dimsizes(y_d)-1
561	delta_y = -1.d
562	break
563	else
564	ydim = 0
565	delta_y = -1.d
566	end if
567	end if
568	end do
569	end if
570	if (yzc .EQ. 1) then
571	do varn=0,dim-1
572	if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
573	y_d = f_att->$vNam(varn)$
574	ydim = dimsizes(y_d)-1
575	delta_y = y_d(1)-y_d(0)
576	break
577	end if
578	end do
579	do varn=0,dim-1
580	if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
581	z_d = f_att->$vNam(varn)$
582	zdim = dimsizes(z_d)-1
583	delta_z = z_d(1)-z_d(0)
584	break
585	end if
586	end do
587	do varn=0,dim-1
588	if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
589	x_d = f_att->$vNam(varn)$
590	xdim = dimsizes(x_d)-1
591	delta_x = x_d(1)-x_d(0)
592	break
593	else
594	if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
595	x_d = f_att->$vNam(varn)$
596	xdim = dimsizes(x_d)-1
597	delta_x = -1.d
598	break
599	else
600	xdim = 0
601	delta_x = -1.d
602	end if
603	end if
604	end do
605	end if
606
607	; ****************************************************
608	; set up ranges of x-, y- and z-coordinates
609	; ****************************************************
610
611	if (xs .EQ. -1.d) then
612	xs = x_d(0)
613	if (delta_x .EQ. -1) then
614	print(" ")
615	print("You cannot choose a start value for x, there are preseted layers for x")
616	print(" ")
617	xstart=0
618	end if
619	else
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	else
626	if (xs .LT. 0-delta_x/2) then
627	print(" ")
628	print("range start for x-coordinate = "+xs+"m is lower than first value = "+(0-delta_x/2)+"m")
629	print(" ")
630	exit
631	end if
632	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
633	if (xs .GE. xdim*delta_x) then
634	print(" ")
635	print("range start for x-coordinate = "+xs+"m is equal or greater than last value = "+xdim*delta_x+"m")
636	print(" ")
637	exit
638	end if
639	else
640	if (xs .GT. xdim*delta_x) then
641	print(" ")
642	print("range start for x-coordinate = "+xs+"m is greater than last value = "+xdim*delta_x+"m")
643	print(" ")
644	exit
645	end if
646	end if
647	end if
648	end if
649
650	if (ys .EQ. -1.d) then
651	ys = y_d(0)
652	if (delta_y .EQ. -1) then
653	print(" ")
654	print("You cannot choose a start value for y, there are preseted layers for y")
655	print(" ")
656	ystart=0
657	end if
658	else
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	else
665	if (ys .LT. 0-delta_y/2) then
666	print(" ")
667	print("range start for y-coordinate = "+ys+"m is lower than first value = "+0-delta_y/2+"m")
668	print(" ")
669	exit
670	end if
671	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
672	if (ys .GE. ydim*delta_y) then
673	print(" ")
674	print("range start for y-coordinate = "+ys+"m is equal or greater than last value = "+ydim*delta_y+"m")
675	print(" ")
676	exit
677	end if
678	else
679	if (ys .GT. ydim*delta_y) then
680	print(" ")
681	print("range start for y-coordinate = "+ys+"m is greater than last value = "+ydim*delta_y+"m")
682	print(" ")
683	exit
684	end if
685	end if
686	end if
687	end if
688
689	if (zs .EQ. -1) then
690	zs = 0
691	if (delta_z .EQ. -1) then
692	print(" ")
693	print("You cannot choose a start value for z, there are preseted layers for z")
694	print(" ")
695	end if
696	else
697	if (delta_z .EQ. -1) then
698	print(" ")
699	print("You cannot choose a start value for z, there are preseted layers for z")
700	print(" ")
701	else
702	if (zs .LT. 0) then
703	print(" ")
704	print("range start for z-coordinate = "+zs+" is lower than first gridpoint = 0")
705	print(" ")
706	exit
707	end if
708	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
709	if (zs .GE. zdim) then
710	print(" ")
711	print("range start for z-coordinate = "+zs+" is equal or greater than last gridpoint = "+zdim)
712	print(" ")
713	exit
714	end if
715	else
716	if (zs .GT. zdim) then
717	print(" ")
718	print("range start for z-coordinate = "+zs+" is greater than last gridpoint = "+zdim)
719	print(" ")
720	exit
721	end if
722	end if
723	end if
724	end if
725
726	if (xe .EQ. -1) then
727	xe = x_d(xdim)
728	if (delta_x .EQ. -1) then
729	print(" ")
730	print("You cannot choose an end value for x, there are preseted layers for x")
731	print(" ")
732	xend=xdim
733	end if
734	else
735	if (delta_x .EQ. -1) then
736	print(" ")
737	print("You cannot choose an end value for x, there are preseted layers for x")
738	print(" ")
739	xend=xdim
740	else
741	if (xe .GT. xdim*delta_x) then
742	print(" ")
743	print("range end for x-coordinate = "+xe+"m is greater than last value = "+xdim*delta_x+"m")
744	print(" ")
745	exit
746	end if
747	if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
748	if (xe .LE. 0-delta_x/2)
749	print(" ")
750	print("range end for x-coordinate = "+xe+"m is equal or lower than first value = "+(0-delta_x/2)+"m")
751	print(" ")
752	exit
753	end if
754	if (xe .LE. xs) then
755	print(" ")
756	print("range end for x-coordinate = "+xe+"m is equal or lower than start range = "+xs+"m")
757	print(" ")
758	exit
759	end if
760	if (xe .EQ. xs+1) then
761	print(" ")
762	print("range end for x-coordinate = "+xe+"m must be at least two more gridpoints greater than start range = "+xs+"m")
763	print(" ")
764	exit
765	end if
766	else
767	if (xe .LT. 0-delta_x/2)
768	print(" ")
769	print("range end for x-coordinate = "+xe+"m is lower than first value = "+(0-delta_x/2)+"m")
770	print(" ")
771	exit
772	end if
773	if (xe .LT. xs) then
774	print(" ")
775	print("range end for x-coordinate = "+xe+"m is lower than start range = "+xs+"m")
776	print(" ")
777	exit
778	end if
779	end if
780	end if
781	end if
782
783	if (ye .EQ. -1) then
784	ye = y_d(ydim)
785	if (delta_y .EQ. -1) then
786	print(" ")
787	print("You cannot choose an end value for y, there are preseted layers for y")
788	print(" ")
789	yend=ydim
790	end if
791	else
792	if (delta_y .EQ. -1) then
793	print(" ")
794	print("You cannot choose an end value for y, there are preseted layers for y")
795	print(" ")
796	yend=ydim
797	else
798	if (ye .GT. ydim*delta_y) then
799	print(" ")
800	print("range end for y-coordinate = "+ye+"m is greater than last value = "+ydim*delta_y+"m")
801	print(" ")
802	exit
803	end if
804	if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
805	if (ye .LE. 0-delta_y/2)
806	print(" ")
807	print("range end for y-coordinate = "+ye+"m is equal or lower than first value = "+(0-delta_y/2)+"m")
808	print(" ")
809	exit
810	end if
811	if (ye .LE. ys) then
812	print(" ")
813	print("range end for y-coordinate = "+ye+"m is equal or lower than start range = "+ys+"m")
814	print(" ")
815	exit
816	end if
817	if (ye .EQ. ys+1) then
818	print(" ")
819	print("range end for y-coordinate = "+ye+"m must be at least two more gridpoints greater than start range = "+ys+"m")
820	print(" ")
821	exit
822	end if
823	else
824	if (ye .LT. 0-delta_y/2)
825	print(" ")
826	print("range end for y-coordinate = "+ye+"m is lower than first value = "+(0-delta_y/2)+"m")
827	print(" ")
828	exit
829	end if
830	if (ye .LT. ys) then
831	print(" ")
832	print("range end for y-coordinate = "+ye+"m is lower than start range = "+ys+"m")
833	print(" ")
834	exit
835	end if
836	end if
837	end if
838	end if
839
840	if (ze .EQ. -1) then
841	ze = zdim
842	if (delta_z .EQ. -1) then
843	print(" ")
844	print("You cannot choose an end value for z, there are preseted layers for z")
845	print(" ")
846	ze = zdim
847	end if
848	else
849	if (delta_z .EQ. -1) then
850	print(" ")
851	print("You cannot choose an end value for z, there are preseted layers for z")
852	print(" ")
853	ze = zdim
854	else
855	if (ze .GT. zdim) then
856	print(" ")
857	print("range end for z-coordinate = "+ze+" is greater than last gridpoint = "+zdim)
858	print(" ")
859	exit
860	end if
861	if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
862	if (ze .LE. 0)
863	print(" ")
864	print("range end for z-coordinate = "+ze+" is equal or lower than first gridpoint = 0")
865	print(" ")
866	exit
867	end if
868	if (ze .LE. zs) then
869	print(" ")
870	print("range end for z-coordinate = "+ze+" is equal or lower than start range = "+zs)
871	print(" ")
872	exit
873	end if
874	if (ze .EQ. zs+1) then
875	print(" ")
876	print("range end for z-coordinate = "+ze+" must be at least two more gridpoints greater than start range = "+zs)
877	print(" ")
878	exit
879	end if
880	else
881	if (ze .LT. 0)
882	print(" ")
883	print("range end for z-coordinate = "+ze+" is lower than first gridpoint = 0")
884	print(" ")
885	exit
886	end if
887	if (ze .LT. zs) then
888	print(" ")
889	print("range end for z-coordinate = "+ze+" is lower than start range = "+zs)
890	print(" ")
891	exit
892	end if
893	end if
894	end if
895	end if
896
897	if (delta_x .NE. -1) then
898	do i=0,xdim
899	if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
900	xstart=i
901	break
902	end if
903	end do
904	do i=0,xdim
905	if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
906	xend=i
907	break
908	end if
909	end do
910	end if
911	if (delta_y .NE. -1) then
912	do i=0,ydim
913	if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
914	ystart=i
915	break
916	end if
917	end do
918	do i=0,ydim
919	if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
920	yend=i
921	break
922	end if
923	end do
924	end if
925
926	if( shape .eq. 1 ) then
927	if (xyc .EQ. 1)then
928	cs_res@vpWidthF = (xe-xs)/(ye-ys)
929	cs_res@vpHeightF = 1
930	vecres@vpWidthF = (xe-xs)/(ye-ys)
931	vecres@vpHeightF = 1
932	if (xe-xs .GT. ye-ys)then
933	if (no_rows .LE. no_columns)then
934	cs_res@gsnPaperOrientation = "landscape"
935	vecres@gsnPaperOrientation = "landscape"
936	cs_res@lbOrientation = "Horizontal"
937	else
938	cs_res@gsnPaperOrientation = "portrait"
939	vecres@gsnPaperOrientation = "portrait"
940	cs_res@lbOrientation = "Vertical"
941	end if
942	else
943	if (no_rows .GE. no_columns)then
944	cs_res@gsnPaperOrientation = "portrait"
945	vecres@gsnPaperOrientation = "portrait"
946	cs_res@lbOrientation = "Vertical"
947	else
948	cs_res@gsnPaperOrientation = "landscape"
949	vecres@gsnPaperOrientation = "landscape"
950	cs_res@lbOrientation = "Horizontal"
951	end if
952	end if
953	end if
954	if (xzc .EQ. 1)then
955	cs_res@vpWidthF = (xe-xs)/(delta_x*(ze-zs))
956	cs_res@vpHeightF = 1
957	vecres@vpWidthF = (xe-xs)/(delta_x*(ze-zs))
958	vecres@vpHeightF = 1
959	if (xe-xs .GT. (delta_x*(ze-zs)))then
960	if (no_rows .LE. no_columns)then
961	cs_res@gsnPaperOrientation = "landscape"
962	vecres@gsnPaperOrientation = "landscape"
963	cs_res@lbOrientation = "Horizontal"
964	else
965	cs_res@gsnPaperOrientation = "portrait"
966	vecres@gsnPaperOrientation = "portrait"
967	cs_res@lbOrientation = "Vertical"
968	end if
969	else
970	if (no_rows .GE. no_columns)then
971	cs_res@gsnPaperOrientation = "portrait"
972	vecres@gsnPaperOrientation = "portrait"
973	cs_res@lbOrientation = "Vertical"
974	else
975	cs_res@gsnPaperOrientation = "landscape"
976	vecres@gsnPaperOrientation = "landscape"
977	cs_res@lbOrientation = "Horizontal"
978	end if
979	end if
980	end if
981	if (yzc .EQ. 1)then
982	cs_res@vpWidthF = (ye-ys)/(delta_y*(ze-zs))
983	cs_res@vpHeightF = 1
984	vecres@vpWidthF = (ye-ys)/(delta_y*(ze-zs))
985	vecres@vpHeightF = 1
986	if (ye-ys .GT. (delta_y*(ze-zs)))then
987	if (no_rows .LE. no_columns)then
988	cs_res@gsnPaperOrientation = "landscape"
989	vecres@gsnPaperOrientation = "landscape"
990	cs_res@lbOrientation = "Horizontal"
991	else
992	cs_res@gsnPaperOrientation = "portrait"
993	vecres@gsnPaperOrientation = "portrait"
994	cs_res@lbOrientation = "Vertical"
995	end if
996	else
997	if (no_rows .GE. no_columns)then
998	cs_res@gsnPaperOrientation = "portrait"
999	vecres@gsnPaperOrientation = "portrait"
1000	cs_res@lbOrientation = "Vertical"
1001	else
1002	cs_res@gsnPaperOrientation = "landscape"
1003	vecres@gsnPaperOrientation = "landscape"
1004	cs_res@lbOrientation = "Horizontal"
1005	end if
1006	end if
1007	end if
1008	end if
1009
1010	delete(xs)
1011	delete(xe)
1012	delete(ys)
1013	delete(ye)
1014
1015	xs=xstart
1016	xe=xend
1017	ys=ystart
1018	ye=yend
1019
1020	if (xyc .EQ. 1) then
1021	d=(ye-ys+1)/(major_ticks_y-1)
1022	e=(xe-xs+1)/(major_ticks_x-1)
1023	array_yl =new(major_ticks_y,integer)
1024	array_yl_labels=new(major_ticks_y,double)
1025	array_minor_yl =new((major_ticks_y-1)*4,double)
1026	array_xb =new(major_ticks_x,integer)
1027	array_xb_labels=new(major_ticks_x,double)
1028	array_minor_xb =new((major_ticks_x-1)*4,double)
1029	array_yl(0)=ys
1030	array_xb(0)=xs
1031	array_yl_labels(0)=y_d(ys)
1032	array_xb_labels(0)=x_d(xs)
1033	do ar=1,major_ticks_y-1
1034	array_yl(ar)=d*(ar-1)+d-1
1035	array_yl_labels(ar) = y_d(array_yl(ar))
1036	if (ar .GT. 0)
1037	do min_ar=0,3
1038	array_minor_yl(4(ar-1)+min_ar)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-array_yl(ar-1))/5(min_ar+1)
1039	end do
1040	end if
1041	end do
1042	do br=1,major_ticks_x-1
1043	array_xb(br)=e*(br-1)+e-1
1044	array_xb_labels(br) = x_d(array_xb(br))
1045	if (br .GT. 0)
1046	do min_br=0,3
1047	array_minor_xb(4(br-1)+min_br)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-array_xb(br-1))/5(min_br+1)
1048	end do
1049	end if
1050	end do
1051	end if
1052
1053	if (xzc .EQ. 1) then
1054	d=(ze-zs+1)/(major_ticks_y-1)
1055	e=(xe-xs+1)/(major_ticks_x-1)
1056	array_yl =new(major_ticks_y,integer)
1057	array_yl_labels=new(major_ticks_y,double)
1058	array_minor_yl =new((major_ticks_y-1)*4,double)
1059	array_xb =new(major_ticks_x,integer)
1060	array_xb_labels=new(major_ticks_x,double)
1061	array_minor_xb =new((major_ticks_x-1)*4,double)
1062	array_yl(0)=zs
1063	array_xb(0)=xs
1064	array_yl_labels(0)=z_d(zs)
1065	array_xb_labels(0)=x_d(xs)
1066	do ar=1,major_ticks_y-1
1067	array_yl(ar)=d*(ar-1)+d-1
1068	array_yl_labels(ar) = z_d(array_yl(ar))
1069	if (ar .GT. 0)
1070	do min_ar=0,3
1071	array_minor_yl(4(ar-1)+min_ar)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-array_yl(ar-1))/5(min_ar+1)
1072	end do
1073	end if
1074	end do
1075	do br=1,major_ticks_x-1
1076	array_xb(br)=e*(br-1)+e-1
1077	array_xb_labels(br) = x_d(array_xb(br))
1078	if (br .GT. 0)
1079	do min_br=0,3
1080	array_minor_xb(4(br-1)+min_br)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-array_xb(br-1))/5(min_br+1)
1081	end do
1082	end if
1083	end do
1084	end if
1085
1086	if (yzc .EQ. 1) then
1087	d=(ze-zs+1)/(major_ticks_y-1)
1088	e=(ye-ys+1)/(major_ticks_x-1)
1089	array_yl =new(major_ticks_y,integer)
1090	array_yl_labels=new(major_ticks_y,double)
1091	array_minor_yl =new((major_ticks_y-1)*4,double)
1092	array_xb =new(major_ticks_x,integer)
1093	array_xb_labels=new(major_ticks_x,double)
1094	array_minor_xb =new((major_ticks_x-1)*4,double)
1095	array_yl(0)=zs
1096	array_xb(0)=ys
1097	array_yl_labels(0)=z_d(zs)
1098	array_xb_labels(0)=y_d(ys)
1099	do ar=1,major_ticks_y-1
1100	array_yl(ar)=d*(ar-1)+d-1
1101	array_yl_labels(ar) = y_d(array_yl(ar))
1102	if (ar .GT. 0)
1103	do min_ar=0,3
1104	array_minor_yl(4(ar-1)+min_ar)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-array_yl(ar-1))/5(min_ar+1)
1105	end do
1106	end if
1107	end do
1108	do br=1,major_ticks_x-1
1109	array_xb(br)=e*(br-1)+e-1
1110	array_xb_labels(br) = x_d(array_xb(br))
1111	if (br .GT. 0)
1112	do min_br=0,3
1113	array_minor_xb(4(br-1)+min_br)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-array_xb(br-1))/5(min_br+1)
1114	end do
1115	end if
1116	end do
1117	end if
1118
1119	if (axes_explicit .EQ. 1)then
1120	cs_res@tmYLMode = "Explicit"
1121	cs_res@tmXBMode = "Explicit"
1122	cs_res@tmYLValues = array_yl
1123	if (xyc .EQ. 1)then
1124	cs_res@tmYLLabels = array_yl_labels/norm_y
1125	cs_res@tmXBLabels = array_xb_labels/norm_x
1126	if (norm_x .NE. 1.)then
1127	cs_res@tiXAxisString = "x ["+unit_x+"]"
1128	else
1129	cs_res@tiXAxisString = "x [m]"
1130	end if
1131	if (norm_y .NE. 1.)then
1132	cs_res@tiYAxisString = "y ["+unit_y+"]"
1133	else
1134	cs_res@tiYAxisString = "y [m]"
1135	end if
1136	end if
1137	if (xzc .EQ. 1)then
1138	cs_res@tmYLLabels = array_yl_labels/norm_z
1139	cs_res@tmXBLabels = array_xb_labels/norm_x
1140	if (norm_x .NE. 1.)then
1141	cs_res@tiXAxisString = "x ["+unit_x+"]"
1142	else
1143	cs_res@tiXAxisString = "x [m]"
1144	end if
1145	if (norm_z .NE. 1.)then
1146	cs_res@tiYAxisString = "z ["+unit_z+"]"
1147	else
1148	cs_res@tiYAxisString = "z [m]"
1149	end if
1150	end if
1151	if (yzc .EQ. 1)then
1152	cs_res@tmYLLabels = array_yl_labels/norm_z
1153	cs_res@tmXBLabels = array_xb_labels/norm_y
1154	if (norm_y .NE. 1.)then
1155	cs_res@tiXAxisString = "y ["+unit_y+"]"
1156	else
1157	cs_res@tiXAxisString = "y [m]"
1158	end if
1159	if (norm_z .NE. 1.)then
1160	cs_res@tiYAxisString = "z ["+unit_z+"]"
1161	else
1162	cs_res@tiYAxisString = "z [m]"
1163	end if
1164	end if
1165	cs_res@tmXBValues = array_xb
1166	cs_res@tmYLMinorValues = array_minor_yl
1167	cs_res@tmXBMinorValues = array_minor_xb
1168	else
1169	if (axes_explicit .EQ. 0)then
1170	cs_res@tmYLMinorPerMajor = 4
1171	cs_res@tmXBMinorPerMajor = 4
1172	else
1173	print(" ")
1174	print("'axes_explicit' is invalid and set to 0")
1175	print(" ")
1176	axes_explicit = 0
1177	cs_res@tmYLMinorPerMajor = 4
1178	cs_res@tmXBMinorPerMajor = 4
1179	end if
1180	if (xyc .EQ. 1)then
1181	cs_res@tiXAxisString = "x [m]"
1182	cs_res@tiYAxisString = "y [m]"
1183	end if
1184	if (xzc .EQ. 1)then
1185	cs_res@tiXAxisString = "x [m]"
1186	cs_res@tiYAxisString = "z [m]"
1187	end if
1188	if (yzc .EQ. 1)then
1189	cs_res@tiXAxisString = "y [m]"
1190	cs_res@tiYAxisString = "z [m]"
1191	end if
1192	end if
1193
1194	if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1195	if (xe .EQ. xs+1) then
1196	print(" ")
1197	print("range end for x-coordinate="+xe*delta_x+"m) must be at least two")
1198	print("more gridpoints(="+2delta_x+"m) greater than start range="+xsdelta_x+"m)")
1199	print(" ")
1200	exit
1201	end if
1202	end if
1203	if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1204	if (ye .EQ. ys+1) then
1205	print(" ")
1206	print("range end for y-coordinate="+ye*delta_y+"m) must be at least two")
1207	print("more gridpoints(="+2delta_y+"m) greater than start range="+ysdelta_y+"m)")
1208	print(" ")
1209	exit
1210	end if
1211	end if
1212	if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1213	if (ze .EQ. zs+1) then
1214	print(" ")
1215	print("range end for x-coordinate="+ze+") must be at least two")
1216	print("more gridpoints greater than start range="+zs+")")
1217	print(" ")
1218	exit
1219	end if
1220	end if
1221
1222	if (xyc .EQ. 1) then
1223	no_layer = (ze-zs)+1
1224	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1225	end if
1226	if (xzc .EQ. 1) then
1227	no_layer = (ye-ys)+1
1228	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1229	end if
1230	if (yzc .EQ. 1) then
1231	no_layer = (xe-xs)+1
1232	data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1233	end if
1234
1235	MinVal = new(dim,float)
1236	MaxVal = new(dim,float)
1237	unit = new(dim,string)
1238
1239	; ****************************************************
1240	; define inner and outer loops depending on "sort"
1241	; ****************************************************
1242
1243	if ( xyc .eq. 1 ) then
1244	lays = zs
1245	laye = ze
1246	end if
1247	if ( xzc .eq. 1 ) then
1248	lays = ys
1249	laye = ye
1250	end if
1251	if ( yzc .eq. 1) then
1252	lays = xs
1253	laye = xe
1254	end if
1255
1256	if ( sort .eq. "time" ) then
1257	lis = start_time_step
1258	lie = end_time_step
1259	los = 0
1260	loe = laye-lays
1261	else
1262	lis = 0
1263	lie = laye-lays
1264	los = start_time_step
1265	loe = end_time_step
1266	end if
1267
1268	check = True
1269	v1=0
1270	v2=0
1271	no_var=0
1272	n=0
1273
1274	do varn=dim-1,0,1
1275
1276	if (vector .EQ. 1) then
1277	check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1278	check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1279	end if
1280
1281	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
1282	check = False
1283	else
1284	check = True
1285	end if
1286
1287	if(check) then
1288
1289	no_var=no_var+1
1290
1291	if (vector .EQ. 1) then
1292	if (","+vNam(varn)+"," .EQ. vec1) then
1293	v1=v1+1
1294	end if
1295	if (","+vNam(varn)+"," .EQ. vec2) then
1296	v2=v2+1
1297	end if
1298	end if
1299
1300	if (xyc .EQ. 1) then
1301	temp = f[:]->$vNam(varn)$
1302	data_att = f_att->$vNam(varn)$
1303	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1304	end if
1305	if ( xzc .eq. 1 ) then
1306	temp = f[:]->$vNam(varn)$
1307	data_att = f_att->$vNam(varn)$
1308	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1309	end if
1310	if ( yzc .eq. 1) then
1311	temp = f[:]->$vNam(varn)$
1312	data_att = f_att->$vNam(varn)$
1313	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1314	end if
1315	if (check_vec1) then
1316	vect1=data(varn,:,:,:,:)
1317	end if
1318	if (check_vec2) then
1319	vect2=data(varn,:,:,:,:)
1320	end if
1321
1322	data!0 = "var"
1323	data!1 = "t"
1324	data!2 = "z"
1325	data!3 = "y"
1326	data!4 = "x"
1327
1328	MinVal(varn) = min(data(varn,:,:,:,:))
1329	MaxVal(varn) = max(data(varn,:,:,:,:))
1330
1331	unit(varn) = data_att@units
1332
1333	end if
1334
1335	end do
1336
1337	var_input=new(no_var,string)
1338	numb=0
1339	no_var=0
1340
1341	do varn=dim-1,0,1
1342
1343	if (vector .EQ. 1) then
1344	check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1345	check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1346	end if
1347
1348	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
1349	check = False
1350	else
1351	check = True
1352	end if
1353
1354	if(check) then
1355	var_input(numb)=vNam(varn)
1356	numb=numb+1
1357	end if
1358
1359	if (var .NE. "all") then
1360	check = isStrSubset( var,","+vNam(varn)+"," )
1361	end if
1362
1363	if (check) then
1364	no_var = no_var+1
1365	end if
1366
1367	end do
1368
1369	if (no_var .EQ. 0) then
1370	print(" ")
1371	print("The variables var='"+var+"' do not exist on your input file;")
1372	print("be sure to have one comma berfore and after each variable")
1373	print(" ")
1374	exit
1375	end if
1376
1377	if (vector .EQ. 1) then
1378	if (v1 .EQ. 0)then
1379	print(" ")
1380	print("Component 1 for the vector-plot ('vec1') must be one of the varibles on the input file:")
1381	print("- "+var_input)
1382	print("be sure to have one comma berfore and after the variable")
1383	print(" ")
1384	exit
1385	end if
1386
1387	if (v2 .EQ. 0)then
1388	print(" ")
1389	print("Component 2 for the vector-plot ('vec2') must be one of the varibles on the input file:")
1390	print("- "+var_input)
1391	print("be sure to have one comma berfore and after the variable")
1392	print(" ")
1393	exit
1394	end if
1395	end if
1396
1397	; ***************************************************
1398	; open workstation(s)
1399	; ***************************************************
1400
1401	wks_ps = gsn_open_wks(format_out,file_out)
1402	gsn_define_colormap(wks_ps,"rainbow+white")
1403
1404	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1405	plot=new((/no_timeno_layerno_var+no_time*no_layer/),graphic)
1406	else
1407	plot=new((/no_timeno_layerno_var/),graphic)
1408	end if
1409
1410	page = 0
1411	n=0
1412	print(" ")
1413	print("Plots sorted by '"+sort+"'")
1414	print(" ")
1415
1416	; ***************************************************
1417	; create plots
1418	; ***************************************************
1419
1420	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1421	do lo = los, loe
1422	do li = lis, lie
1423	if (xyc .EQ. 1)then
1424	if (sort .EQ. "time")then
1425	level = "z=" + z_d(lo) + "m"
1426	else
1427	level = "z=" + z_d(li) + "m"
1428	end if
1429	end if
1430	if (xzc .EQ. 1)then
1431	if (sort .EQ. "time")then
1432	level = "y=" + y_d(lo) + "m"
1433	else
1434	level = "y=" + y_d(li) + "m"
1435	end if
1436	end if
1437	if (yzc .EQ. 1)then
1438	if (sort .EQ. "time")then
1439	level = "x=" + x_d(lo) + "m"
1440	else
1441	level = "x=" + x_d(li) + "m"
1442	end if
1443	end if
1444	vecres = True ; vector only resources
1445	vecres@gsnDraw = False ; don't draw
1446	vecres@gsnFrame = False ; don't advance frame
1447	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1448	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1449	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1450	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
1451	vecres@tmXBLabelFontHeightF = font_size
1452	vecres@tmYLLabelFontHeightF = font_size
1453	vecres@tiXAxisFontHeightF = font_size
1454	vecres@tiYAxisFontHeightF = font_size
1455	if (sort .EQ. "time")then
1456	vecres@gsnLeftString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h "+level
1457	else
1458	vecres@gsnLeftString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) +"h "+level
1459	end if
1460	vecres@tiXAxisString = " "
1461	if (xyc .EQ. 1)then
1462	if (sort .EQ. "time")then
1463	plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1464	else
1465	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1466	end if
1467	end if
1468	if (xzc .EQ. 1) then
1469	if (sort .EQ. "time")then
1470	plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
1471	else
1472	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
1473	end if
1474	end if
1475	if (yzc .EQ. 1) then
1476	if (sort .EQ. "time")then
1477	plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
1478	else
1479	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
1480	end if
1481	end if
1482	n=n+1
1483	end do
1484	end do
1485	end if
1486
1487	do varn=dim-1,0,1
1488
1489	if (vector .EQ. 1) then
1490	check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1491	end if
1492
1493	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
1494	check = False
1495	else
1496	check = True
1497	end if
1498
1499	if (var .NE. "all") then
1500	check = isStrSubset( var,","+vNam(varn)+"," )
1501	end if
1502
1503	if(check) then
1504
1505	space=(MaxVal(varn)-MinVal(varn))/24
1506
1507	cs_res@cnMinLevelValF = MinVal(varn)
1508	cs_res@cnMaxLevelValF = MaxVal(varn)
1509
1510	cs_res@cnLevelSpacingF = space
1511
1512	; ****************************************************
1513	; loops over time and layer
1514	; ****************************************************
1515
1516
1517	do lo = los, loe
1518	do li = lis, lie
1519
1520	; ****************************************************
1521	; xy cross section
1522	; ****************************************************
1523
1524	if ( xyc .eq. 1 ) then
1525
1526	cs_res@gsnLeftString = "Plot of "+vNam(varn)
1527	cs_res@gsnRightString = unit(varn)
1528
1529	if ( sort .eq. "time" ) then
1530	if ( z_d(zs+lo) .eq. -1)then
1531	if (delta_z .EQ. -1) then
1532	level = "-average"
1533	else
1534	level = "=" + z_d(zs+lo) + "m"
1535	end if
1536	else
1537	level = "=" + z_d(zs+lo) + "m"
1538	end if
1539	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h z"+level
1540	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res)
1541	if (vector .EQ. 1 .AND. check_vecp) then
1542	vecres = True ; vector only resources
1543	vecres@gsnDraw = False ; don't draw
1544	vecres@gsnFrame = False ; don't advance frame
1545	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1546	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1547	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1548	vecres@gsnRightString = " "
1549	vecres@gsnLeftString = " "
1550	vecres@tiXAxisString = " "
1551	plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1552	overlay(plot(n), plot_vec)
1553	end if
1554	end if
1555
1556	if ( sort .eq. "layer" ) then
1557	if (z_d(zs+li) .eq. -1) then
1558	if (delta_z .EQ. -1) then
1559	level = "-average"
1560	else
1561	level = "=" + z_d(zs+li) + "m"
1562	end if
1563	else
1564	level = "=" + z_d(zs+li) + "m"
1565	end if
1566	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h z"+ level
1567	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res)
1568	if (vector .EQ. 1 .AND. check_vecp) then
1569	vecres = True ; vector only resources
1570	vecres@gsnDraw = False ; don't draw
1571	vecres@gsnFrame = False ; don't advance frame
1572	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1573	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1574	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1575	vecres@gsnRightString = " " ; turn off right string
1576	vecres@gsnLeftString = " " ; turn off left string
1577	vecres@tiXAxisString = " "
1578	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1579	overlay(plot(n), plot_vec)
1580	end if
1581	end if
1582	end if
1583
1584	; ****************************************************
1585	; xz cross section
1586	; ****************************************************
1587
1588	if ( xzc .eq. 1 ) then
1589
1590	cs_res@gsnLeftString = "Plot of "+vNam(varn)
1591
1592	if ( sort .eq. "time" ) then
1593	if ( y_d(ys+lo) .eq. -1 ) then
1594	level = "-average"
1595	else
1596	level = "=" + y_d(ys+lo) + "m"
1597	end if
1598	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h y"+ level
1599	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
1600	if (vector .EQ. 1 .AND. check_vecp) then
1601	vecres = True ; vector only resources
1602	vecres@gsnDraw = False ; don't draw
1603	vecres@gsnFrame = False ; don't advance frame
1604	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1605	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1606	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1607	vecres@gsnRightString = " " ; turn off right string
1608	vecres@gsnLeftString = " " ; turn off left string
1609	vecres@tiXAxisString = " "
1610	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
1611	overlay(plot(n), plot_vec)
1612	end if
1613	end if
1614
1615	if ( sort .eq. "layer" ) then
1616	if ( y_d(ys+li) .eq. -1 ) then
1617	level = "-average"
1618	else
1619	level = "=" + y_d(ys+li) + "m"
1620	end if
1621	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h y"+ level
1622	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
1623	if (vector .EQ. 1 .AND. check_vecp) then
1624	vecres = True ; vector only resources
1625	vecres@gsnDraw = False ; don't draw
1626	vecres@gsnFrame = False ; don't advance frame
1627	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1628	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1629	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1630	vecres@gsnRightString = " " ; turn off right string
1631	vecres@gsnLeftString = " " ; turn off left string
1632	vecres@tiXAxisString = " "
1633	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
1634	overlay(plot(n), plot_vec)
1635	end if
1636	end if
1637	end if
1638
1639	; ****************************************************
1640	; yz cross section
1641	; ****************************************************
1642
1643	if ( yzc .eq. 1 ) then
1644
1645	cs_res@gsnLeftString = "Plot of "+vNam(varn)
1646
1647	if ( sort .eq. "time" ) then
1648	if ( x_d(xs+lo) .eq. -1 ) then
1649	level = "-average"
1650	else
1651	level = "=" + x_d(xs+lo) + "m"
1652	end if
1653	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h x"+ level
1654	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
1655	if (vector .EQ. 1 .AND. check_vecp) then
1656	vecres = True ; vector only resources
1657	vecres@gsnDraw = False ; don't draw
1658	vecres@gsnFrame = False ; don't advance frame
1659	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1660	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1661	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1662	vecres@gsnRightString = " " ; turn off right string
1663	vecres@gsnLeftString = " " ; turn off left string
1664	vecres@tiXAxisString = " "
1665	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
1666	overlay(plot(n), plot_vec)
1667	end if
1668	end if
1669
1670	if ( sort .eq. "layer" ) then
1671	if ( x_d(xs+li) .eq. -1 ) then
1672	level = "-average"
1673	else
1674	level = "=" + x_d(xs+li) + "m"
1675	end if
1676	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h x"+ level
1677	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
1678	if (vector .EQ. 1 .AND. check_vecp)then
1679	vecres = True ; vector only resources
1680	vecres@gsnDraw = False ; don't draw
1681	vecres@gsnFrame = False ; don't advance frame
1682	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1683	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1684	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1685	vecres@gsnRightString = " " ; turn off right string
1686	vecres@gsnLeftString = " " ; turn off left string
1687	vecres@tiXAxisString = " "
1688	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
1689	overlay(plot(n), plot_vec)
1690	end if
1691	end if
1692	end if
1693	n=n+1
1694	end do
1695	end do
1696	end if
1697	end do
1698
1699	; ***************************************************
1700	; merge plots onto one page
1701	; ***************************************************
1702
1703	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec")then
1704	if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1705	gsn_panel(wks_ps,plot(0:(no_timeno_layer(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP)
1706	print(" ")
1707	print("Outputs to .eps or .epsi have only one frame")
1708	print(" ")
1709	else
1710	do np = 0,no_layerno_time(no_var+1)-1,no_rows*no_columns
1711	if ( np + no_rowsno_columns .gt. (no_layerno_time*(no_var+1))-1) then
1712	gsn_panel(wks_ps, plot(np:(no_layerno_time(no_var+1))-1),(/no_rows,no_columns/),cs_resP)
1713	else
1714	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
1715	end if
1716	end do
1717	end if
1718	else
1719	if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1720	gsn_panel(wks_ps,plot(0:(no_timeno_layerno_var)-1),(/no_var,no_layer*no_time/),cs_resP)
1721	print(" ")
1722	print("Outputs to .eps or .epsi have only one frame")
1723	print(" ")
1724	else
1725	do np = 0,no_layerno_timeno_var-1,no_rows*no_columns
1726	if ( np + no_rowsno_columns .gt. (no_layerno_time*no_var)-1) then
1727	gsn_panel(wks_ps, plot(np:(no_layerno_timeno_var)-1),(/no_rows,no_columns/),cs_resP)
1728	else
1729	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
1730	end if
1731	end do
1732	end if
1733	end if
1734
1735	print(" ")
1736	print("Output to: " + file_out +"."+ format_out)
1737	print(" ")
1738
1739	end

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

Download in other formats:

Original Format