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

Last change on this file since 317 was 286, checked in by weinreis, 15 years ago

ncl-script cross_sections.ncl modified

File size: 59.5 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)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-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)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-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)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-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)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-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) = y_d(array_yl(ar))
1102         if (ar .GT. 0)
1103            do min_ar=0,3
1104               array_minor_yl(4*(ar-1)+min_ar)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-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) = x_d(array_xb(br))
1111         if (br .GT. 0)
1112            do min_br=0,3
1113               array_minor_xb(4*(br-1)+min_br)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-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      cs_res@tmYLValues = array_yl
1123      if (xyc .EQ. 1)then
1124         cs_res@tmYLLabels = array_yl_labels/norm_y
1125         cs_res@tmXBLabels = array_xb_labels/norm_x
1126         if (norm_x .NE. 1.)then
1127            cs_res@tiXAxisString = "x ["+unit_x+"]"
1128         else
1129            cs_res@tiXAxisString = "x [m]"
1130         end if
1131         if (norm_y .NE. 1.)then
1132            cs_res@tiYAxisString = "y ["+unit_y+"]"
1133         else
1134            cs_res@tiYAxisString = "y [m]"
1135         end if   
1136      end if
1137      if (xzc .EQ. 1)then
1138         cs_res@tmYLLabels = array_yl_labels/norm_z
1139         cs_res@tmXBLabels = array_xb_labels/norm_x
1140         if (norm_x .NE. 1.)then
1141            cs_res@tiXAxisString = "x ["+unit_x+"]"
1142         else
1143            cs_res@tiXAxisString = "x [m]"
1144         end if
1145         if (norm_z .NE. 1.)then
1146            cs_res@tiYAxisString = "z ["+unit_z+"]"
1147         else
1148            cs_res@tiYAxisString = "z [m]"
1149         end if
1150      end if
1151      if (yzc .EQ. 1)then
1152         cs_res@tmYLLabels = array_yl_labels/norm_z
1153         cs_res@tmXBLabels = array_xb_labels/norm_y
1154         if (norm_y .NE. 1.)then
1155            cs_res@tiXAxisString = "y ["+unit_y+"]"
1156         else
1157            cs_res@tiXAxisString = "y [m]"
1158         end if
1159         if (norm_z .NE. 1.)then
1160            cs_res@tiYAxisString = "z ["+unit_z+"]"
1161         else
1162            cs_res@tiYAxisString = "z [m]"
1163         end if
1164      end if
1165      cs_res@tmXBValues = array_xb     
1166      cs_res@tmYLMinorValues = array_minor_yl
1167      cs_res@tmXBMinorValues = array_minor_xb
1168   else
1169      if (axes_explicit .EQ. 0)then 
1170         cs_res@tmYLMinorPerMajor = 4
1171         cs_res@tmXBMinorPerMajor = 4
1172      else
1173         print(" ")
1174         print("'axes_explicit' is invalid and set to 0")
1175         print(" ")
1176         axes_explicit = 0
1177         cs_res@tmYLMinorPerMajor = 4
1178         cs_res@tmXBMinorPerMajor = 4
1179      end if
1180      if (xyc .EQ. 1)then
1181         cs_res@tiXAxisString = "x [m]"
1182         cs_res@tiYAxisString = "y [m]"
1183      end if
1184      if (xzc .EQ. 1)then
1185         cs_res@tiXAxisString = "x [m]"
1186         cs_res@tiYAxisString = "z [m]"
1187      end if
1188      if (yzc .EQ. 1)then
1189         cs_res@tiXAxisString = "y [m]"
1190         cs_res@tiYAxisString = "z [m]"
1191      end if
1192   end if
1193
1194   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
1195      if (xe .EQ. xs+1) then
1196         print(" ")
1197         print("range end for x-coordinate="+xe*delta_x+"m) must be at least two")
1198         print("more gridpoints(="+2*delta_x+"m) greater than start range="+xs*delta_x+"m)")
1199         print(" ")
1200         exit
1201      end if
1202   end if
1203   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
1204      if (ye .EQ. ys+1) then
1205         print(" ")
1206         print("range end for y-coordinate="+ye*delta_y+"m) must be at least two")
1207         print("more gridpoints(="+2*delta_y+"m) greater than start range="+ys*delta_y+"m)")
1208         print(" ")
1209         exit
1210      end if
1211   end if
1212   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
1213      if (ze .EQ. zs+1) then
1214         print(" ")
1215         print("range end for x-coordinate="+ze+") must be at least two")
1216         print("more gridpoints greater than start range="+zs+")")
1217         print(" ")
1218         exit
1219      end if
1220   end if
1221 
1222   if (xyc .EQ. 1) then
1223      no_layer = (ze-zs)+1
1224      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1225   end if
1226   if (xzc .EQ. 1) then
1227      no_layer = (ye-ys)+1
1228      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1229   end if
1230   if (yzc .EQ. 1) then
1231      no_layer = (xe-xs)+1
1232      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
1233   end if
1234
1235   MinVal = new(dim,float)
1236   MaxVal = new(dim,float)
1237   unit   = new(dim,string)
1238
1239   ; ****************************************************
1240   ; define inner and outer loops depending on "sort"
1241   ; ****************************************************       
1242   
1243   if ( xyc .eq. 1 ) then
1244      lays = zs
1245      laye = ze
1246   end if
1247   if ( xzc .eq. 1 ) then
1248      lays = ys
1249      laye = ye
1250   end if
1251   if ( yzc .eq. 1) then
1252      lays = xs
1253      laye = xe
1254   end if
1255
1256   if ( sort .eq. "time" ) then
1257      lis = start_time_step
1258      lie = end_time_step
1259      los = 0
1260      loe = laye-lays
1261   else
1262      lis = 0
1263      lie = laye-lays
1264      los = start_time_step
1265      loe = end_time_step
1266   end if
1267
1268   check = True
1269   v1=0
1270   v2=0
1271   no_var=0
1272   n=0
1273
1274   do varn=dim-1,0,1       
1275   
1276      if (vector .EQ. 1) then   
1277         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1278         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1279      end if
1280           
1281      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
1282         check = False
1283      else
1284         check = True
1285      end if 
1286
1287      if(check) then
1288     
1289         no_var=no_var+1
1290
1291         if (vector .EQ. 1) then
1292            if (","+vNam(varn)+"," .EQ. vec1) then
1293               v1=v1+1
1294            end if
1295            if (","+vNam(varn)+"," .EQ. vec2) then
1296               v2=v2+1
1297            end if
1298         end if
1299
1300         if (xyc .EQ. 1) then
1301            temp = f[:]->$vNam(varn)$
1302            data_att = f_att->$vNam(varn)$
1303            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)               
1304         end if
1305         if ( xzc .eq. 1 ) then
1306            temp = f[:]->$vNam(varn)$
1307            data_att = f_att->$vNam(varn)$
1308            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1309         end if
1310         if ( yzc .eq. 1) then
1311            temp = f[:]->$vNam(varn)$
1312            data_att = f_att->$vNam(varn)$
1313            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1314         end if
1315         if (check_vec1) then
1316            vect1=data(varn,:,:,:,:)
1317         end if
1318         if (check_vec2) then
1319            vect2=data(varn,:,:,:,:)
1320         end if
1321
1322         data!0 = "var"
1323         data!1 = "t"
1324         data!2 = "z"
1325         data!3 = "y"
1326         data!4 = "x" 
1327
1328         MinVal(varn) = min(data(varn,:,:,:,:))
1329         MaxVal(varn) = max(data(varn,:,:,:,:))
1330         
1331         unit(varn) = data_att@units
1332
1333      end if
1334
1335   end do
1336
1337   var_input=new(no_var,string)
1338   numb=0
1339   no_var=0
1340   
1341   do varn=dim-1,0,1   
1342   
1343      if (vector .EQ. 1) then   
1344         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1345         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1346      end if
1347           
1348      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
1349         check = False
1350      else
1351         check = True
1352      end if 
1353
1354      if(check) then     
1355         var_input(numb)=vNam(varn)
1356         numb=numb+1     
1357      end if
1358                 
1359      if (var .NE. "all") then
1360         check = isStrSubset( var,","+vNam(varn)+"," )
1361      end if
1362     
1363      if (check) then
1364         no_var = no_var+1
1365      end if
1366     
1367   end do
1368
1369   if (no_var .EQ. 0) then
1370      print(" ")
1371      print("The variables var='"+var+"' do not exist on your input file;")
1372      print("be sure to have one comma berfore and after each variable")
1373      print(" ")
1374      exit
1375   end if
1376
1377   if (vector .EQ. 1) then
1378      if (v1 .EQ. 0)then
1379         print(" ")
1380         print("Component 1 for the vector-plot ('vec1') must be one of the varibles on the input file:")
1381         print("- "+var_input)
1382         print("be sure to have one comma berfore and after the variable")
1383         print(" ")
1384         exit
1385      end if
1386
1387      if (v2 .EQ. 0)then
1388         print(" ")
1389         print("Component 2 for the vector-plot ('vec2') must be one of the varibles on the input file:")
1390         print("- "+var_input)
1391         print("be sure to have one comma berfore and after the variable")
1392         print(" ")
1393         exit
1394      end if
1395   end if
1396
1397   ; ***************************************************
1398   ; open workstation(s)
1399   ; ***************************************************
1400
1401   wks_ps  = gsn_open_wks(format_out,file_out)
1402   gsn_define_colormap(wks_ps,"rainbow+white")
1403
1404   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1405      plot=new((/no_time*no_layer*no_var+no_time*no_layer/),graphic)
1406   else
1407      plot=new((/no_time*no_layer*no_var/),graphic)
1408   end if
1409
1410   page = 0
1411   n=0
1412   print(" ")
1413   print("Plots sorted by '"+sort+"'")
1414   print(" ")
1415
1416   ; ***************************************************
1417   ; create plots
1418   ; ***************************************************
1419
1420   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1421      do lo = los, loe                                         
1422         do li = lis, lie
1423            if (xyc .EQ. 1)then
1424               if (sort .EQ. "time")then
1425                  level = "z=" + z_d(lo) + "m"
1426               else
1427                  level = "z=" + z_d(li) + "m"
1428               end if
1429            end if
1430            if (xzc .EQ. 1)then
1431               if (sort .EQ. "time")then
1432                  level = "y=" + y_d(lo) + "m"
1433               else
1434                  level = "y=" + y_d(li) + "m"
1435               end if
1436            end if
1437            if (yzc .EQ. 1)then
1438               if (sort .EQ. "time")then
1439                  level = "x=" + x_d(lo) + "m"
1440               else
1441                  level = "x=" + x_d(li) + "m"
1442               end if
1443            end if               
1444            vecres                  = True            ; vector only resources
1445            vecres@gsnDraw          = False           ; don't draw
1446            vecres@gsnFrame         = False           ; don't advance frame
1447            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1448            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1449            vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1450            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
1451            vecres@tmXBLabelFontHeightF   = font_size
1452            vecres@tmYLLabelFontHeightF   = font_size
1453            vecres@tiXAxisFontHeightF     = font_size
1454            vecres@tiYAxisFontHeightF     = font_size
1455            if (sort .EQ. "time")then
1456               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h   "+level
1457            else
1458               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level
1459            end if
1460            vecres@tiXAxisString    = " "
1461            if (xyc .EQ. 1)then 
1462               if (sort .EQ. "time")then                                         
1463                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1464               else
1465                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1466               end if
1467            end if
1468            if (xzc .EQ. 1) then
1469               if (sort .EQ. "time")then
1470                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
1471               else
1472                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
1473               end if
1474            end if
1475            if (yzc .EQ. 1) then
1476               if (sort .EQ. "time")then
1477                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
1478               else
1479                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
1480               end if
1481            end if
1482            n=n+1
1483         end do
1484      end do
1485   end if
1486 
1487   do varn=dim-1,0,1
1488
1489      if (vector .EQ. 1) then   
1490         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1491      end if
1492
1493      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
1494         check = False
1495      else
1496         check = True
1497      end if   
1498 
1499      if (var .NE. "all") then
1500         check = isStrSubset( var,","+vNam(varn)+"," )
1501      end if
1502   
1503      if(check) then
1504
1505         space=(MaxVal(varn)-MinVal(varn))/24
1506 
1507         cs_res@cnMinLevelValF = MinVal(varn)
1508         cs_res@cnMaxLevelValF = MaxVal(varn)
1509
1510         cs_res@cnLevelSpacingF = space
1511     
1512         ; ****************************************************
1513         ; loops over time and layer
1514         ; ****************************************************
1515
1516         
1517         do lo = los, loe                                       
1518            do li = lis, lie
1519
1520               ; ****************************************************
1521               ; xy cross section
1522               ; ****************************************************
1523
1524               if ( xyc .eq. 1 ) then
1525         
1526                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1527                  cs_res@gsnRightString = unit(varn)
1528                 
1529                  if ( sort .eq. "time" ) then
1530                     if ( z_d(zs+lo) .eq. -1)then
1531                        if (delta_z .EQ. -1) then
1532                           level = "-average"
1533                        else
1534                           level = "=" + z_d(zs+lo) + "m"
1535                        end if
1536                     else
1537                        level = "=" + z_d(zs+lo) + "m"
1538                     end if
1539                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level             
1540                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res)
1541                     if (vector .EQ. 1 .AND. check_vecp) then
1542                        vecres                  = True            ; vector only resources
1543                        vecres@gsnDraw          = False           ; don't draw
1544                        vecres@gsnFrame         = False           ; don't advance frame
1545                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1546                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1547                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref 
1548                        vecres@gsnRightString   = " "
1549                        vecres@gsnLeftString    = " "
1550                        vecres@tiXAxisString    = " "   
1551                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
1552                        overlay(plot(n), plot_vec)
1553                     end if                         
1554                  end if
1555                 
1556                  if ( sort .eq. "layer" ) then
1557                     if (z_d(zs+li) .eq. -1) then
1558                        if (delta_z .EQ. -1) then
1559                           level = "-average"
1560                        else
1561                           level = "=" + z_d(zs+li) + "m"
1562                        end if
1563                     else
1564                        level = "=" + z_d(zs+li) + "m"
1565                     end if
1566                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  z"+ level               
1567                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res)
1568                     if (vector .EQ. 1 .AND. check_vecp) then
1569                        vecres                  = True            ; vector only resources
1570                        vecres@gsnDraw          = False           ; don't draw
1571                        vecres@gsnFrame         = False           ; don't advance frame
1572                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1573                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1574                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1575                        vecres@gsnRightString   = " "             ; turn off right string
1576                        vecres@gsnLeftString    = " "             ; turn off left string
1577                        vecres@tiXAxisString    = " "   
1578                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
1579                        overlay(plot(n), plot_vec)
1580                     end if
1581                  end if
1582               end if
1583
1584               ; ****************************************************
1585               ; xz cross section
1586               ; ****************************************************
1587
1588               if ( xzc .eq. 1 ) then
1589                 
1590                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1591               
1592                  if ( sort .eq. "time" ) then
1593                     if ( y_d(ys+lo) .eq. -1 ) then
1594                        level = "-average"
1595                     else
1596                        level = "=" + y_d(ys+lo) + "m"
1597                     end if
1598                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
1599                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
1600                     if (vector .EQ. 1 .AND. check_vecp) then
1601                        vecres                  = True            ; vector only resources
1602                        vecres@gsnDraw          = False           ; don't draw
1603                        vecres@gsnFrame         = False           ; don't advance frame
1604                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1605                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1606                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1607                        vecres@gsnRightString   = " "             ; turn off right string
1608                        vecres@gsnLeftString    = " "             ; turn off left string
1609                        vecres@tiXAxisString    = " "   
1610                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
1611                        overlay(plot(n), plot_vec)
1612                     end if
1613                  end if
1614
1615                  if ( sort .eq. "layer" ) then
1616                     if ( y_d(ys+li) .eq. -1 ) then
1617                        level = "-average"
1618                     else
1619                        level = "=" + y_d(ys+li) + "m"
1620                     end if
1621                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
1622                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
1623                     if (vector .EQ. 1 .AND. check_vecp) then
1624                        vecres                  = True            ; vector only resources
1625                        vecres@gsnDraw          = False           ; don't draw
1626                        vecres@gsnFrame         = False           ; don't advance frame
1627                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1628                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1629                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1630                        vecres@gsnRightString   = " "             ; turn off right string
1631                        vecres@gsnLeftString    = " "             ; turn off left string
1632                        vecres@tiXAxisString    = " "   
1633                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
1634                        overlay(plot(n), plot_vec)
1635                     end if
1636                  end if                 
1637               end if
1638
1639               ; ****************************************************
1640               ; yz cross section
1641               ; ****************************************************
1642
1643               if ( yzc .eq. 1 ) then
1644                                 
1645                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1646
1647                  if ( sort .eq. "time" ) then
1648                     if ( x_d(xs+lo) .eq. -1 ) then
1649                        level = "-average"
1650                     else
1651                        level = "=" + x_d(xs+lo) + "m"
1652                     end if
1653                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
1654                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
1655                     if (vector .EQ. 1 .AND. check_vecp) then
1656                        vecres                  = True            ; vector only resources
1657                        vecres@gsnDraw          = False           ; don't draw
1658                        vecres@gsnFrame         = False           ; don't advance frame
1659                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1660                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1661                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1662                        vecres@gsnRightString   = " "             ; turn off right string
1663                        vecres@gsnLeftString    = " "             ; turn off left string
1664                        vecres@tiXAxisString    = " "   
1665                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
1666                        overlay(plot(n), plot_vec)
1667                     end if
1668                  end if
1669
1670                  if ( sort .eq. "layer" ) then
1671                     if ( x_d(xs+li) .eq. -1 ) then
1672                        level = "-average"
1673                     else
1674                        level = "=" + x_d(xs+li) + "m"
1675                     end if
1676                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
1677                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
1678                     if (vector .EQ. 1 .AND. check_vecp)then
1679                        vecres                  = True            ; vector only resources
1680                        vecres@gsnDraw          = False           ; don't draw
1681                        vecres@gsnFrame         = False           ; don't advance frame
1682                        vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1683                        vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1684                        vecres@vcRefLengthF     = 0.05            ; define length of vec ref
1685                        vecres@gsnRightString   = " "             ; turn off right string
1686                        vecres@gsnLeftString    = " "             ; turn off left string
1687                        vecres@tiXAxisString    = " "   
1688                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
1689                        overlay(plot(n), plot_vec)
1690                     end if
1691                  end if
1692               end if         
1693               n=n+1 
1694            end do     
1695         end do 
1696      end if
1697   end do
1698
1699   ; ***************************************************
1700   ; merge plots onto one page
1701   ; ***************************************************
1702   
1703   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec")then
1704      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1705         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),(/no_var+1,no_layer*no_time/),cs_resP)
1706         print(" ")
1707         print("Outputs to .eps or .epsi have only one frame")
1708         print(" ")
1709      else
1710         do np = 0,no_layer*no_time*(no_var+1)-1,no_rows*no_columns   
1711            if ( np + no_rows*no_columns .gt. (no_layer*no_time*(no_var+1))-1) then
1712               gsn_panel(wks_ps, plot(np:(no_layer*no_time*(no_var+1))-1),(/no_rows,no_columns/),cs_resP)
1713            else
1714               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
1715            end if
1716         end do
1717      end if
1718   else       
1719      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
1720         gsn_panel(wks_ps,plot(0:(no_time*no_layer*no_var)-1),(/no_var,no_layer*no_time/),cs_resP)
1721         print(" ")
1722         print("Outputs to .eps or .epsi have only one frame")
1723         print(" ")
1724      else
1725         do np = 0,no_layer*no_time*no_var-1,no_rows*no_columns 
1726            if ( np + no_rows*no_columns .gt. (no_layer*no_time*no_var)-1) then
1727               gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_var)-1),(/no_rows,no_columns/),cs_resP)
1728            else
1729               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
1730            end if
1731         end do
1732      end if
1733   end if
1734
1735   print(" ")
1736   print("Output to: " + file_out +"."+ format_out)
1737   print(" ")   
1738 
1739end
Note: See TracBrowser for help on using the repository browser.