source: palm/trunk/SCRIPTS/NCL/cross_sections.ncl @ 424

Last change on this file since 424 was 418, checked in by heinze, 12 years ago

ncl_preferences is renamed .ncl.config.default

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