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

Last change on this file since 758 was 758, checked in by heinze, 13 years ago

Bugfixes for NCL scripts timeseries.ncl and cross_sections.ncl concerning usage of norm_t and vector plot

File size: 75.6 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 ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1551           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1552           vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or.      \
1553           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1554           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1555           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1556           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1557           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1558           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1559           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1560           vNam(varn) .eq. "ind_x_yz") then
1561         check = False
1562      else
1563         check = True
1564      end if 
1565
1566      if (var .NE. "all") then
1567         check = isStrSubset( var,","+vNam(varn)+"," )
1568      end if
1569
1570      if(check) then     
1571         var_input(numb)=vNam(varn)
1572         numb=numb+1     
1573      end if
1574     
1575      if (check) then
1576         no_var = no_var+1
1577      end if
1578     
1579   end do
1580
1581   if (no_var .EQ. 0) then
1582      print(" ")
1583      print("The variables var='"+var+"' do not exist on your input file;")
1584      print("be sure to have one comma before and after each variable")
1585      print(" ")
1586      exit
1587   end if
1588
1589   if (vector .EQ. 1) then
1590      if (v1 .EQ. 0)then
1591         print(" ")
1592         print("Component 1 for the vector-plot ('vec1') must be one "+\
1593               "of the variables on the input file:")
1594         print("- "+var_input)
1595         print("be sure to have one comma before and after the variable")
1596         print(" ")
1597         exit
1598      end if
1599
1600      if (v2 .EQ. 0)then
1601         print(" ")
1602         print("Component 2 for the vector-plot ('vec2') must be one "+\
1603               "of the variables on the input file:")
1604         print("- "+var_input)
1605         print("be sure to have one comma before and after the variable")
1606         print(" ")
1607         exit
1608      end if
1609   end if
1610
1611   ; ***************************************************
1612   ; open workstation(s)
1613   ; ***************************************************
1614
1615   if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
1616      format_out@wkPaperSize = "A4"
1617   end if
1618   if ( format_out .EQ. "png" ) then
1619      format_out@wkWidth  = 1000
1620      format_out@wkHeight = 1000
1621   end if
1622
1623   wks_ps  = gsn_open_wks(format_out,file_out)
1624   gsn_define_colormap(wks_ps,"rainbow+white")
1625
1626   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1627      plot=new((/no_time*no_layer/),graphic)
1628   else
1629      plot=new((/no_time*no_layer*no_var - no_time*(no_layer-1)*no_zu1/)\
1630                                                                ,graphic)
1631   end if
1632   dim_plot=dimsizes(plot)
1633
1634   page = 0
1635   n=0
1636   print(" ")
1637   print("Plots sorted by '"+sort+"'")
1638   print(" ")
1639
1640   ; ***************************************************
1641   ; create plots
1642   ; ***************************************************
1643
1644
1645   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1646      do lo = los, loe                                         
1647         do li = lis, lie
1648            if (xyc .EQ. 1)then
1649               if (sort .EQ. "time")then
1650                  if(z_d(lo) .eq. -1.d) then
1651                    level = "z-average"
1652                  else
1653                    level = "z=" + z_d(lo) + "m"
1654                  end if
1655               else
1656                  if(z_d(li) .eq. -1.d) then
1657                    level = "z-average"
1658                  else
1659                    level = "z=" + z_d(li) + "m"
1660                  end if
1661               end if
1662            end if
1663            if (xzc .EQ. 1)then
1664               if (sort .EQ. "time")then
1665                 if(y_d(lo) .eq. -1.d) then
1666                   level = "y-average"
1667                 else
1668                   level = "y=" + y_d(lo) + "m"
1669                 end if
1670               else
1671                  if(y_d(li) .eq. -1.d) then
1672                    level = "y-average"
1673                  else
1674                    level = "y=" + y_d(li) + "m"
1675                  end if
1676               end if
1677            end if
1678            if (yzc .EQ. 1)then
1679               if (sort .EQ. "time")then
1680                 if(x_d(lo) .eq. -1.d) then
1681                    level = "x-average"
1682                 else
1683                    level = "x=" + x_d(lo) + "m"
1684                 end if
1685               else
1686                 if(x_d(li) .eq. -1.d) then
1687                    level = "x-average"
1688                 else
1689                    level = "x=" + x_d(li) + "m"
1690                 end if
1691               end if
1692            end if               
1693            vecres                  = True            ; vector only resources
1694            vecres@gsnDraw          = False           ; don't draw
1695            vecres@gsnFrame         = False           ; don't advance frame
1696            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1697            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1698            vecres@vcRefLengthF     = 0.05            ; define length of
1699                                                      ; vec ref
1700            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
1701            vecres@tmXBLabelFontHeightF   = font_size
1702            vecres@tmYLLabelFontHeightF   = font_size
1703            vecres@tiXAxisFontHeightF     = font_size
1704            vecres@tiYAxisFontHeightF     = font_size
1705            if (sort .EQ. "time")then
1706               vecres@gsnLeftString = "t=" + \
1707                             decimalPlaces(t_all(li)/3600,2,True) +"h   "+level
1708            else
1709               vecres@gsnLeftString = "t=" + \
1710                             decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level
1711            end if
1712            if (xyc .EQ. 1) then 
1713               vecres@tiXAxisString = "x (m)"
1714               vecres@tiYAxisString = "y (m)"   
1715               if (sort .EQ. "time")then
1716                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),\
1717                                               vect2(li,lo-los,:,:),vecres)
1718               else
1719                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1720                                                vect2(lo,li-lis,:,:),vecres)
1721               end if
1722            end if
1723            if (xzc .EQ. 1) then
1724               vecres@tiXAxisString = "x (m)"
1725               vecres@tiYAxisString = "z (m)"
1726               if (sort .EQ. "time")then
1727                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
1728                                               vect2(li,:,lo-los,:),vecres)
1729               else
1730                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
1731                                                vect2(lo,:,li-lis,:),vecres)
1732               end if
1733            end if
1734            if (yzc .EQ. 1) then
1735               vecres@tiXAxisString = "y (m)"
1736               vecres@tiYAxisString = "z (m)"
1737               if (sort .EQ. "time")then
1738                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
1739                                                vect2(li,:,:,lo-los),vecres)
1740               else
1741                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
1742                                                vect2(lo,:,:,li-lis),vecres)
1743               end if
1744            end if
1745            n=n+1
1746         end do
1747      end do
1748   end if
1749 
1750   do varn=dim-1,0,1
1751
1752      if (vector .EQ. 1 ) then   
1753         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1754      end if
1755     
1756      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1757           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1758           vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or.      \
1759           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1760           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1761           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1762           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1763           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1764           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1765           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1766           vNam(varn) .eq. "ind_x_yz") then
1767         check = False
1768      else
1769         check = True
1770      end if   
1771 
1772      if (var .NE. "all") then
1773         check = isStrSubset( var,","+vNam(varn)+"," )
1774      end if
1775   
1776      if(check) then
1777
1778         space=(MaxVal(varn)-MinVal(varn))/24
1779 
1780         cs_res@cnMinLevelValF = MinVal(varn)
1781         cs_res@cnMaxLevelValF = MaxVal(varn)
1782
1783         cs_res@cnLevelSpacingF = space
1784     
1785         ; ****************************************************
1786         ; loops over time and layer
1787         ; ****************************************************
1788
1789         
1790         do lo = los, loe                                       
1791            do li = lis, lie
1792
1793               ; ****************************************************
1794               ; xy cross section
1795               ; ****************************************************
1796
1797               if ( xyc .eq. 1 ) then
1798         
1799                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1800                  cs_res@gsnRightString = unit(varn)
1801                 
1802                  if ( sort .eq. "time" ) then
1803                     if ( z_d(lo) .eq. -1)then
1804                        if (delta_z .EQ. -1) then
1805                           level = "-average"
1806                        else
1807                           level = "=" + z_d(lo) + "m"
1808                        end if
1809                     else
1810                        level = "=" + z_d(lo) + "m"
1811                     end if
1812
1813                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1814                        vNam(varn) .eq. "pras_xy"  .or. \
1815                        vNam(varn) .eq. "prrs_xy"  .or. \
1816                        vNam(varn) .eq. "qsws_xy"  .or. \
1817                        vNam(varn) .eq. "shfs_xy"  .or. \
1818                        vNam(varn) .eq. "ts_xy"    .or. \
1819                        vNam(varn) .eq. "us_xy"    .or. \
1820                        vNam(varn) .eq. "z0s_xy"   .or. \
1821                        vNam(varn) .eq. "lwp*_xy"  .or. \
1822                        vNam(varn) .eq. "pra*_xy"  .or. \
1823                        vNam(varn) .eq. "prr*_xy"  .or. \
1824                        vNam(varn) .eq. "qsws*_xy" .or. \
1825                        vNam(varn) .eq. "shf*_xy"  .or. \
1826                        vNam(varn) .eq. "t*_xy"    .or. \
1827                        vNam(varn) .eq. "u*_xy"    .or. \
1828                        vNam(varn) .eq. "z0*_xy") then
1829                         loe = 0
1830                         level = "=" + zu1(0) + "m"
1831                     end if
1832                   
1833                     cs_res@gsnCenterString = "t=" + \
1834                       decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level
1835
1836                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1837                        ;nothing to be done
1838                     else
1839                         plot(n) = gsn_csm_contour(wks_ps,\
1840                                           data(varn,li,lo-los,:,:),cs_res)
1841                     end if
1842
1843
1844                     if (vector .EQ. 1 .AND. check_vecp) then
1845                        vecres                 = True           ; vector only
1846                                                                ; resources
1847                        vecres@gsnDraw         = False          ; don't draw
1848                        vecres@gsnFrame        = False          ; don't
1849                                                                ; advance frame
1850                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly
1851                                                                ; vectors
1852                        vecres@vcRefMagnitudeF = ref_mag        ; define
1853                                                                ; vector ref
1854                                                                ; mag
1855                        vecres@vcRefLengthF    = 0.05           ; define
1856                                                                ; length of
1857                                                                ; vec ref 
1858                        vecres@gsnRightString  = " "
1859                        vecres@gsnLeftString   = " "
1860                        vecres@tiXAxisString   = " "   
1861                        plot_vec=gsn_csm_vector(wks_ps,\
1862                              vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1863                        overlay(plot(n), plot_vec)
1864                     end if                         
1865                  end if
1866                 
1867                  if ( sort .eq. "layer" ) then
1868                     if (z_d(li) .eq. -1) then
1869                        if (delta_z .EQ. -1) then
1870                           level = "-average"
1871                        else
1872                           level = "=" + z_d(li) + "m"
1873                        end if
1874                     else
1875                        level = "=" + z_d(li) + "m"
1876                     end if
1877       
1878                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1879                        vNam(varn) .eq. "pras_xy"  .or. \
1880                        vNam(varn) .eq. "prrs_xy"  .or. \
1881                        vNam(varn) .eq. "qsws_xy"  .or. \
1882                        vNam(varn) .eq. "shfs_xy"  .or. \
1883                        vNam(varn) .eq. "ts_xy"    .or. \
1884                        vNam(varn) .eq. "us_xy"    .or. \
1885                        vNam(varn) .eq. "z0s_xy"   .or. \
1886                        vNam(varn) .eq. "lwp*_xy"  .or. \
1887                        vNam(varn) .eq. "pra*_xy"  .or. \
1888                        vNam(varn) .eq. "prr*_xy"  .or. \
1889                        vNam(varn) .eq. "qsws*_xy" .or. \
1890                        vNam(varn) .eq. "shf*_xy"  .or. \
1891                        vNam(varn) .eq. "t*_xy"    .or. \
1892                        vNam(varn) .eq. "u*_xy"    .or. \
1893                        vNam(varn) .eq. "z0*_xy") then
1894                         lie = 0
1895                         level = "=" + zu1(0) + "m"
1896                     end if
1897               
1898                     cs_res@gsnCenterString = "t=" + \
1899                           decimalPlaces(t_all(lo)/3600,2,True) + "h  z"+ level
1900
1901                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1902                        ;nothing to be done
1903                     else
1904                        plot(n) = gsn_csm_contour(wks_ps,\
1905                                            data(varn,lo,li-lis,:,:),cs_res)
1906                     end if
1907
1908                     if (vector .EQ. 1 .AND. check_vecp) then
1909                        vecres                 = True           ; vector only
1910                                                                ; resources
1911                        vecres@gsnDraw         = False          ; don't draw
1912                        vecres@gsnFrame        = False          ; don't
1913                                                                ; advance frame
1914                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
1915                        vecres@vcRefMagnitudeF = ref_mag        ; define
1916                                                                ; vector ref
1917                                                                ; mag
1918                        vecres@vcRefLengthF    = 0.05           ; define
1919                                                                ; length of
1920                                                                ; vec ref
1921                        vecres@gsnRightString  = " "            ; turn off
1922                                                                ; right string
1923                        vecres@gsnLeftString   = " "            ; turn off
1924                                                                ; left string
1925                        vecres@tiXAxisString   = " "   
1926                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1927                                                   vect2(lo,li-lis,:,:),vecres)
1928                        overlay(plot(n), plot_vec)
1929                     end if
1930                  end if
1931               end if
1932
1933               ; ****************************************************
1934               ; xz cross section
1935               ; ****************************************************
1936
1937               if ( xzc .eq. 1 ) then
1938                 
1939                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1940               
1941                  if ( sort .eq. "time" ) then
1942                     if ( y_d(lo) .eq. -1 ) then
1943                        level = "-average"
1944                     else
1945                        level = "=" + y_d(lo) + "m"
1946                     end if
1947                     cs_res@gsnCenterString = "t=" + \
1948                           decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
1949
1950                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1951                        ;nothing to be done
1952                     else
1953                         plot(n) = gsn_csm_contour(wks_ps,\
1954                                              data(varn,li,:,lo-los,:),cs_res)
1955                     end if
1956
1957                     if (vector .EQ. 1 .AND. check_vecp) then
1958                        vecres                 = True           ; vector only
1959                                                                ; resources
1960                        vecres@gsnDraw         = False          ; don't draw
1961                        vecres@gsnFrame        = False          ; don't
1962                                                                ; advance frame
1963                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
1964                        vecres@vcRefMagnitudeF = ref_mag        ; define
1965                                                                ; vector ref
1966                                                                ; mag
1967                        vecres@vcRefLengthF    = 0.05           ; define
1968                                                                ; length of
1969                                                                ; vec ref
1970                        vecres@gsnRightString  = " "            ; turn off
1971                                                                ; right string
1972                        vecres@gsnLeftString   = " "            ; turn off
1973                                                                ; left string
1974                        vecres@tiXAxisString   = " "   
1975                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
1976                                                   vect2(li,:,lo-los,:),vecres)
1977                        overlay(plot(n), plot_vec)
1978                     end if
1979                  end if
1980
1981                  if ( sort .eq. "layer" ) then
1982                     if ( y_d(li) .eq. -1 ) then
1983                        level = "-average"
1984                     else
1985                        level = "=" + y_d(li) + "m"
1986                     end if
1987                     cs_res@gsnCenterString = "t=" + \
1988                           decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
1989
1990                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1991                        ;nothing to be done
1992                     else
1993                         plot(n) = gsn_csm_contour(wks_ps,\
1994                                              data(varn,lo,:,li-lis,:),cs_res) 
1995                     end if
1996
1997                     if (vector .EQ. 1 .AND. check_vecp) then
1998                        vecres                 = True           ; vector only
1999                                                                ; resources
2000                        vecres@gsnDraw         = False          ; don't draw
2001                        vecres@gsnFrame        = False          ; don't
2002                                                                ; advance frame
2003                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2004                        vecres@vcRefMagnitudeF = ref_mag        ; define
2005                                                                ; vector ref
2006                                                                ; mag
2007                        vecres@vcRefLengthF    = 0.05           ; define
2008                                                                ; length of
2009                                                                ; vec ref
2010                        vecres@gsnRightString  = " "            ; turn off
2011                                                                ; right string
2012                        vecres@gsnLeftString   = " "            ; turn off
2013                                                                ; left string
2014                        vecres@tiXAxisString   = " "   
2015                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
2016                                                   vect2(lo,:,li-lis,:),vecres)
2017                        overlay(plot(n), plot_vec)
2018                     end if
2019                  end if                 
2020               end if
2021
2022               ; ****************************************************
2023               ; yz cross section
2024               ; ****************************************************
2025
2026               if ( yzc .eq. 1 ) then
2027                                 
2028                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
2029
2030                  if ( sort .eq. "time" ) then
2031                     if ( x_d(lo) .eq. -1 ) then
2032                        level = "-average"
2033                     else
2034                        level = "=" + x_d(lo) + "m"
2035                     end if
2036                     cs_res@gsnCenterString = "t=" + \
2037                          decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
2038
2039                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2040                        ;nothing to be done
2041                     else
2042                        plot(n) = gsn_csm_contour(wks_ps,\
2043                                          data(varn,li,:,:,lo-los),cs_res)
2044                     end if
2045
2046                     if (vector .EQ. 1 .AND. check_vecp) then
2047                        vecres                 = True           ; vector only
2048                                                                ; resources
2049                        vecres@gsnDraw         = False          ; don't draw
2050                        vecres@gsnFrame        = False          ; don't
2051                                                                ; advance frame
2052                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2053                        vecres@vcRefMagnitudeF = ref_mag        ; define
2054                                                                ; vector ref
2055                                                                ; mag
2056                        vecres@vcRefLengthF    = 0.05           ; define
2057                                                                ; length of
2058                                                                ; vec ref
2059                        vecres@gsnRightString  = " "            ; turn off
2060                                                                ; right string
2061                        vecres@gsnLeftString   = " "            ; turn off
2062                                                                ; left string
2063                        vecres@tiXAxisString   = " "   
2064                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
2065                                                   vect2(li,:,:,lo-los),vecres)
2066                        overlay(plot(n), plot_vec)
2067                     end if
2068                  end if
2069
2070                  if ( sort .eq. "layer" ) then
2071                     if ( x_d(li) .eq. -1 ) then
2072                        level = "-average"
2073                     else
2074                        level = "=" + x_d(li) + "m"
2075                     end if
2076                     cs_res@gsnCenterString = "t=" + \
2077                           decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
2078
2079                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2080                        ;nothing to be done
2081                     else
2082                        plot(n) = gsn_csm_contour(wks_ps,\
2083                                          data(varn,lo,:,:,li-lis),cs_res)
2084                     end if
2085
2086                     if (vector .EQ. 1 .AND. check_vecp)then
2087                        vecres                 = True           ; vector only
2088                                                                ;resources
2089                        vecres@gsnDraw         = False          ; don't draw
2090                        vecres@gsnFrame        = False          ; don't
2091                                                                ; advance frame
2092                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2093                        vecres@vcRefMagnitudeF = ref_mag        ; define
2094                                                                ; vector ref
2095                                                                ; mag
2096                        vecres@vcRefLengthF    = 0.05           ; define
2097                                                                ; length of
2098                                                                ; vec ref
2099                        vecres@gsnRightString  = " "            ; turn off
2100                                                                ; right string
2101                        vecres@gsnLeftString   = " "            ; turn off
2102                                                                ; left string
2103                        vecres@tiXAxisString   = " "   
2104                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
2105                                                   vect2(lo,:,:,li-lis),vecres)
2106                        overlay(plot(n), plot_vec)
2107                     end if
2108                  end if
2109               end if         
2110               n=n+1 
2111            end do     
2112         end do 
2113      end if
2114   end do
2115
2116   ; ***************************************************
2117   ; merge plots onto one page
2118   ; ***************************************************   
2119
2120   no_frames = 0
2121
2122   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2123      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
2124           no_time*no_layer*(no_var+1) .gt. no_rows*no_columns) then
2125         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),\
2126                                       (/no_var+1,no_layer*no_time/),cs_resP)
2127         print(" ")
2128         print("Outputs to .eps or .epsi have only one frame")
2129         print(" ")
2130      else
2131         do np = 0,dim_plot-1,no_rows*no_columns   
2132            if ( np + no_rows*no_columns .gt. dim_plot-1) then
2133               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2134                                                                      cs_resP)
2135               no_frames = no_frames + 1
2136            else
2137               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2138                                              (/no_rows,no_columns/),cs_resP)
2139               no_frames = no_frames + 1
2140            end if
2141         end do
2142      end if
2143   else       
2144      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
2145                               .AND. dim_plot .gt. no_rows*no_columns) then
2146         gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
2147         print(" ")
2148         print("Outputs to .eps or .epsi have only one frame")
2149         print(" ")
2150      else
2151         do np = 0,dim_plot-1,no_rows*no_columns 
2152            if ( np + no_rows*no_columns .gt. dim_plot-1) then
2153               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2154                                                                     cs_resP)
2155               no_frames = no_frames + 1
2156            else
2157               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2158                                              (/no_rows,no_columns/),cs_resP)
2159               no_frames = no_frames + 1
2160            end if
2161         end do
2162      end if
2163   end if
2164
2165   if (format_out .EQ. "png" ) then
2166     png_output = new((/no_frames/), string)
2167     j = 0
2168     do i=0, no_frames-1
2169       j = i + 1
2170       if (j .LE. 9) then
2171         png_output(i) = file_out+".00000"+j+".png"
2172       end if
2173       if (j .GT. 9 .AND. j .LE. 99) then
2174         png_output(i) = file_out+".0000"+j+".png"
2175       end if
2176       if (j .GT. 99 .AND. j .LE. 999) then
2177         png_output(i) = file_out+".000"+j+".png"
2178       end if
2179       if (j .GT. 999) then
2180         png_output(i) = file_out+".00"+j+".png"
2181       end if
2182
2183       ;using imagemagick's convert for reducing the white
2184       ;space around the plot
2185       cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
2186              png_output(i) + " " + png_output(i)
2187       system(cmd)
2188     end do
2189
2190     print(" ")
2191     print("Output to: "+ png_output)
2192     print(" ")
2193   else
2194     print(" ")
2195     print("Output to: " + file_out +"."+ format_out)
2196     print(" ")
2197   end if
2198 
2199end
Note: See TracBrowser for help on using the repository browser.