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

Last change on this file since 983 was 983, checked in by hoffmann, 10 years ago

Bugfix in netcdf.f90 and improvements in palmplot.

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