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

Last change on this file since 1264 was 1248, checked in by heinze, 11 years ago

NCL function getfilevarnames changed with NCL version 6.1.1 and higher -> adaptation required

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