Home

Context Navigation

cross_sections.ncl @ 526

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

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

Download in other formats:

Original Format