Home

Context Navigation

cross_sections.ncl

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

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

Download in other formats:

Original Format