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

Last change on this file since 331 was 329, checked in by heinze, 16 years ago

Small change in formatting in message.f90. Bugfix in cross_sections.ncl in case of normalizing.

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