Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 580 was 566, checked in by heinze, 14 years ago

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