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

Last change on this file since 742 was 585, checked in by heinze, 14 years ago

Bugfix: enable plot of data if it is of kind double instead of kind float

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