Home

Context Navigation

cross_sections.ncl @ 524

Last change on this file since 524 was 514, checked in by heinze, 15 years ago
Bugfix concerning the plot of xy cross sections which have only one layer (zu1_xy)
File size: 63.0 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. "lwps_xy" .or. vNam(varn) .eq. "pras_xy" \
1321	.or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \
1322	.or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy" \
1323	.or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy") then
1324	;these variables depend on zu1_xy and that's why they have only one z-layer
1325	data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe)
1326	no_zu1=no_zu1+1
1327	else
1328	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1329	end if
1330	delete(temp)
1331	end if
1332	if ( xzc .eq. 1 ) then
1333	temp = f[:]->$vNam(varn)$
1334	data_att = f_att->$vNam(varn)$
1335	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1336	delete(temp)
1337	end if
1338	if ( yzc .eq. 1) then
1339	temp = f[:]->$vNam(varn)$
1340	data_att = f_att->$vNam(varn)$
1341	data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1342	delete(temp)
1343	end if
1344	if (check_vec1) then
1345	vect1=data(varn,:,:,:,:)
1346	end if
1347	if (check_vec2) then
1348	vect2=data(varn,:,:,:,:)
1349	end if
1350
1351	data!0 = "var"
1352	data!1 = "t"
1353	data!2 = "z"
1354	data!3 = "y"
1355	data!4 = "x"
1356
1357	MinVal(varn) = min(data(varn,start_time_step:end_time_step,0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1358	MaxVal(varn) = max(data(varn,start_time_step:end_time_step,0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1359
1360	unit(varn) = data_att@units
1361	delete(data_att)
1362
1363	end if
1364
1365	end do
1366
1367	if (no_var .EQ. 0) then
1368	print(" ")
1369	print("The variables var='"+var+"' do not exist on your input file;")
1370	print("be sure to have one comma before and after each variable")
1371	print(" ")
1372	exit
1373	end if
1374
1375	var_input=new(no_var,string)
1376	numb=0
1377	no_var=0
1378
1379	do varn=dim-1,0,1
1380
1381	if (vector .EQ. 1) then
1382	check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1383	check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1384	end if
1385
1386	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
1387	check = False
1388	else
1389	check = True
1390	end if
1391
1392	if (var .NE. "all") then
1393	check = isStrSubset( var,","+vNam(varn)+"," )
1394	end if
1395
1396	if(check) then
1397	var_input(numb)=vNam(varn)
1398	numb=numb+1
1399	end if
1400
1401	if (check) then
1402	no_var = no_var+1
1403	end if
1404
1405	end do
1406
1407	if (no_var .EQ. 0) then
1408	print(" ")
1409	print("The variables var='"+var+"' do not exist on your input file;")
1410	print("be sure to have one comma before and after each variable")
1411	print(" ")
1412	exit
1413	end if
1414
1415	if (vector .EQ. 1) then
1416	if (v1 .EQ. 0)then
1417	print(" ")
1418	print("Component 1 for the vector-plot ('vec1') must be one of the varibles on the input file:")
1419	print("- "+var_input)
1420	print("be sure to have one comma before and after the variable")
1421	print(" ")
1422	exit
1423	end if
1424
1425	if (v2 .EQ. 0)then
1426	print(" ")
1427	print("Component 2 for the vector-plot ('vec2') must be one of the varibles on the input file:")
1428	print("- "+var_input)
1429	print("be sure to have one comma before and after the variable")
1430	print(" ")
1431	exit
1432	end if
1433	end if
1434
1435	; ***************************************************
1436	; open workstation(s)
1437	; ***************************************************
1438
1439	wks_ps = gsn_open_wks(format_out,file_out)
1440	gsn_define_colormap(wks_ps,"rainbow+white")
1441
1442	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1443	plot=new((/no_timeno_layerno_var - no_time(no_layer-1)no_zu1 \
1444	+ no_time*no_layer/),graphic)
1445	else
1446	plot=new((/no_timeno_layerno_var - no_time(no_layer-1)no_zu1/),graphic)
1447	end if
1448	dim_plot=dimsizes(plot)
1449
1450	page = 0
1451	n=0
1452	print(" ")
1453	print("Plots sorted by '"+sort+"'")
1454	print(" ")
1455
1456	; ***************************************************
1457	; create plots
1458	; ***************************************************
1459
1460
1461	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1462	do lo = los, loe
1463	do li = lis, lie
1464	if (xyc .EQ. 1)then
1465	if (sort .EQ. "time")then
1466	if(z_d(lo) .eq. -1.d) then
1467	level = "z-average"
1468	else
1469	level = "z=" + z_d(lo) + "m"
1470	end if
1471	else
1472	if(z_d(li) .eq. -1.d) then
1473	level = "z-average"
1474	else
1475	level = "z=" + z_d(li) + "m"
1476	end if
1477	end if
1478	end if
1479	if (xzc .EQ. 1)then
1480	if (sort .EQ. "time")then
1481	if(y_d(lo) .eq. -1.d) then
1482	level = "y-average"
1483	else
1484	level = "y=" + y_d(lo) + "m"
1485	end if
1486	else
1487	if(y_d(li) .eq. -1.d) then
1488	level = "y-average"
1489	else
1490	level = "y=" + y_d(li) + "m"
1491	end if
1492	end if
1493	end if
1494	if (yzc .EQ. 1)then
1495	if (sort .EQ. "time")then
1496	if(x_d(lo) .eq. -1.d) then
1497	level = "x-average"
1498	else
1499	level = "x=" + x_d(lo) + "m"
1500	end if
1501	else
1502	if(x_d(li) .eq. -1.d) then
1503	level = "x-average"
1504	else
1505	level = "x=" + x_d(li) + "m"
1506	end if
1507	end if
1508	end if
1509	vecres = True ; vector only resources
1510	vecres@gsnDraw = False ; don't draw
1511	vecres@gsnFrame = False ; don't advance frame
1512	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1513	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1514	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1515	vecres@gsnRightString = "Vector plot of "+vec1+" and "+vec2
1516	vecres@tmXBLabelFontHeightF = font_size
1517	vecres@tmYLLabelFontHeightF = font_size
1518	vecres@tiXAxisFontHeightF = font_size
1519	vecres@tiYAxisFontHeightF = font_size
1520	if (sort .EQ. "time")then
1521	vecres@gsnLeftString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h "+level
1522	else
1523	vecres@gsnLeftString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) +"h "+level
1524	end if
1525	vecres@tiXAxisString = " "
1526	if (xyc .EQ. 1)then
1527	if (sort .EQ. "time")then
1528	plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1529	else
1530	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),vect2(lo,li-lis,:,:),vecres)
1531	end if
1532	end if
1533	if (xzc .EQ. 1) then
1534	if (sort .EQ. "time")then
1535	plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),vect2(li,:,lo-los,:),vecres)
1536	else
1537	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),vect2(lo,:,li-lis,:),vecres)
1538	end if
1539	end if
1540	if (yzc .EQ. 1) then
1541	if (sort .EQ. "time")then
1542	plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),vect2(li,:,:,lo-los),vecres)
1543	else
1544	plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),vect2(lo,:,:,li-lis),vecres)
1545	end if
1546	end if
1547	n=n+1
1548	end do
1549	end do
1550	end if
1551
1552	do varn=dim-1,0,1
1553
1554	if (vector .EQ. 1) then
1555	check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1556	end if
1557
1558	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
1559	check = False
1560	else
1561	check = True
1562	end if
1563
1564	if (var .NE. "all") then
1565	check = isStrSubset( var,","+vNam(varn)+"," )
1566	end if
1567
1568	if(check) then
1569
1570	space=(MaxVal(varn)-MinVal(varn))/24
1571
1572	cs_res@cnMinLevelValF = MinVal(varn)
1573	cs_res@cnMaxLevelValF = MaxVal(varn)
1574
1575	cs_res@cnLevelSpacingF = space
1576
1577	; ****************************************************
1578	; loops over time and layer
1579	; ****************************************************
1580
1581
1582	do lo = los, loe
1583	do li = lis, lie
1584
1585	; ****************************************************
1586	; xy cross section
1587	; ****************************************************
1588
1589	if ( xyc .eq. 1 ) then
1590
1591	cs_res@gsnLeftString = "Plot of "+vNam(varn)
1592	cs_res@gsnRightString = unit(varn)
1593
1594	if ( sort .eq. "time" ) then
1595	if ( z_d(lo) .eq. -1)then
1596	if (delta_z .EQ. -1) then
1597	level = "-average"
1598	else
1599	level = "=" + z_d(lo) + "m"
1600	end if
1601	else
1602	level = "=" + z_d(lo) + "m"
1603	end if
1604
1605	if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy" \
1606	.or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \
1607	.or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy" \
1608	.or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy") then
1609	loe = 0
1610	level = "=" + zu1(0) + "m"
1611	end if
1612
1613	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h z"+level
1614	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo-los,:,:),cs_res)
1615	if (vector .EQ. 1 .AND. check_vecp) then
1616	vecres = True ; vector only resources
1617	vecres@gsnDraw = False ; don't draw
1618	vecres@gsnFrame = False ; don't advance frame
1619	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1620	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1621	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1622	vecres@gsnRightString = " "
1623	vecres@gsnLeftString = " "
1624	vecres@tiXAxisString = " "
1625	plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1626	overlay(plot(n), plot_vec)
1627	end if
1628	end if
1629
1630	if ( sort .eq. "layer" ) then
1631	if (z_d(li) .eq. -1) then
1632	if (delta_z .EQ. -1) then
1633	level = "-average"
1634	else
1635	level = "=" + z_d(li) + "m"
1636	end if
1637	else
1638	level = "=" + z_d(li) + "m"
1639	end if
1640
1641	if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy" \
1642	.or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \
1643	.or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy" \
1644	.or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy") then
1645	lie = 0
1646	level = "=" + zu1(0) + "m"
1647	end if
1648
1649	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h z"+ level
1650
1651	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li-lis,:,:),cs_res)
1652	if (vector .EQ. 1 .AND. check_vecp) then
1653	vecres = True ; vector only resources
1654	vecres@gsnDraw = False ; don't draw
1655	vecres@gsnFrame = False ; don't advance frame
1656	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1657	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1658	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1659	vecres@gsnRightString = " " ; turn off right string
1660	vecres@gsnLeftString = " " ; turn off left string
1661	vecres@tiXAxisString = " "
1662	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1663	overlay(plot(n), plot_vec)
1664	end if
1665	end if
1666	end if
1667
1668	; ****************************************************
1669	; xz cross section
1670	; ****************************************************
1671
1672	if ( xzc .eq. 1 ) then
1673
1674	cs_res@gsnLeftString = "Plot of "+vNam(varn)
1675
1676	if ( sort .eq. "time" ) then
1677	if ( y_d(lo) .eq. -1 ) then
1678	level = "-average"
1679	else
1680	level = "=" + y_d(lo) + "m"
1681	end if
1682	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h y"+ level
1683	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo-los,:),cs_res)
1684	if (vector .EQ. 1 .AND. check_vecp) then
1685	vecres = True ; vector only resources
1686	vecres@gsnDraw = False ; don't draw
1687	vecres@gsnFrame = False ; don't advance frame
1688	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1689	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1690	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1691	vecres@gsnRightString = " " ; turn off right string
1692	vecres@gsnLeftString = " " ; turn off left string
1693	vecres@tiXAxisString = " "
1694	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),vect2(li,:,lo-los,:),vecres)
1695	overlay(plot(n), plot_vec)
1696	end if
1697	end if
1698
1699	if ( sort .eq. "layer" ) then
1700	if ( y_d(li) .eq. -1 ) then
1701	level = "-average"
1702	else
1703	level = "=" + y_d(li) + "m"
1704	end if
1705	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h y"+ level
1706	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li-lis,:),cs_res)
1707	if (vector .EQ. 1 .AND. check_vecp) then
1708	vecres = True ; vector only resources
1709	vecres@gsnDraw = False ; don't draw
1710	vecres@gsnFrame = False ; don't advance frame
1711	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1712	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1713	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1714	vecres@gsnRightString = " " ; turn off right string
1715	vecres@gsnLeftString = " " ; turn off left string
1716	vecres@tiXAxisString = " "
1717	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),vect2(lo,:,li-lis,:),vecres)
1718	overlay(plot(n), plot_vec)
1719	end if
1720	end if
1721	end if
1722
1723	; ****************************************************
1724	; yz cross section
1725	; ****************************************************
1726
1727	if ( yzc .eq. 1 ) then
1728
1729	cs_res@gsnLeftString = "Plot of "+vNam(varn)
1730
1731	if ( sort .eq. "time" ) then
1732	if ( x_d(lo) .eq. -1 ) then
1733	level = "-average"
1734	else
1735	level = "=" + x_d(lo) + "m"
1736	end if
1737	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h x"+ level
1738	plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo-los),cs_res)
1739	if (vector .EQ. 1 .AND. check_vecp) then
1740	vecres = True ; vector only resources
1741	vecres@gsnDraw = False ; don't draw
1742	vecres@gsnFrame = False ; don't advance frame
1743	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1744	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1745	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1746	vecres@gsnRightString = " " ; turn off right string
1747	vecres@gsnLeftString = " " ; turn off left string
1748	vecres@tiXAxisString = " "
1749	plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),vect2(li,:,:,lo-los),vecres)
1750	overlay(plot(n), plot_vec)
1751	end if
1752	end if
1753
1754	if ( sort .eq. "layer" ) then
1755	if ( x_d(li) .eq. -1 ) then
1756	level = "-average"
1757	else
1758	level = "=" + x_d(li) + "m"
1759	end if
1760	cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h x"+ level
1761	plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li-lis),cs_res)
1762	if (vector .EQ. 1 .AND. check_vecp)then
1763	vecres = True ; vector only resources
1764	vecres@gsnDraw = False ; don't draw
1765	vecres@gsnFrame = False ; don't advance frame
1766	vecres@vcGlyphStyle = "CurlyVector" ; curly vectors
1767	vecres@vcRefMagnitudeF = ref_mag ; define vector ref mag
1768	vecres@vcRefLengthF = 0.05 ; define length of vec ref
1769	vecres@gsnRightString = " " ; turn off right string
1770	vecres@gsnLeftString = " " ; turn off left string
1771	vecres@tiXAxisString = " "
1772	plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),vect2(lo,:,:,li-lis),vecres)
1773	overlay(plot(n), plot_vec)
1774	end if
1775	end if
1776	end if
1777	n=n+1
1778	end do
1779	end do
1780	end if
1781	end do
1782
1783	; ***************************************************
1784	; merge plots onto one page
1785	; ***************************************************
1786	if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec")then
1787	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. no_timeno_layer(no_var+1) .gt. no_rows*no_columns) then
1788	gsn_panel(wks_ps,plot(0:(no_timeno_layer(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP)
1789	print(" ")
1790	print("Outputs to .eps or .epsi have only one frame")
1791	print(" ")
1792	else
1793	do np = 0,dim_plot-1,no_rows*no_columns
1794	if ( np + no_rows*no_columns .gt. dim_plot-1) then
1795	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),cs_resP)
1796	else
1797	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
1798	end if
1799	end do
1800	end if
1801	else
1802	if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. dim_plot .gt. no_rows*no_columns) then
1803	gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
1804	print(" ")
1805	print("Outputs to .eps or .epsi have only one frame")
1806	print(" ")
1807	else
1808	do np = 0,dim_plot-1,no_rows*no_columns
1809	if ( np + no_rows*no_columns .gt. dim_plot-1) then
1810	gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),cs_resP)
1811	else
1812	gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
1813	end if
1814	end do
1815	end if
1816	end if
1817
1818	print(" ")
1819	print("Output to: " + file_out +"."+ format_out)
1820	print(" ")
1821
1822	end

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

Download in other formats:

Original Format