Home


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

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

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

allow plotting of data with very small time increments
File size: 77.4 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
424   if (nt .EQ. 1)then
425      delta_t=t_all(nt-1)/nt
426   else
427      delta_t_array = new(nt-1,double)
428
429      do i=0,nt-2
430         delta_t_array(i) = t_all(i+1)-t_all(i)
431      end do
432
433      delta_t = min(delta_t_array)
434      delete(delta_t_array)
435   end if
436
437
438   ; ****************************************************       
439   ; start of time step and different types of mistakes that could be done
440   ; ****************************************************
441
442   if (start_time_step .EQ. -1.) then           
443      start_time_step=t_all(0)/3600     
444   else
445      if (start_time_step .GT. t_all(nt-1)/3600)
446         print(" ")
447         print("'start_time_step' = "+ start_time_step +"h is greater than "+\
448               "last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
449         print(" ")
450         print("Select another 'start_time_step'")
451         print(" ")
452         exit
453      end if
454      if (start_time_step .LT. t_all(0)/3600)
455         print(" ")
456         print("'start_time_step' = "+ start_time_step +"h is lower than "+\
457               "first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
458         print(" ")
459         print("Select another 'start_time_step'")
460         print(" ")
461         exit
462      end if
463   end if
464   do i=0,nt-1   
465      if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
466          start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
467         st=i
468         break
469      else
470         st=0
471      end if
472   end do
473
474       
475   ; ****************************************************
476   ; end of time step and different types of mistakes that could be done
477   ; ****************************************************
478
479   if (end_time_step .EQ. -1.) then             
480      end_time_step = t_all(nt-1)/3600 
481   else 
482      if (end_time_step .LT. t_all(0)/3600)
483         print(" ")
484         print("'end_time_step' = "+end_time_step+ "h is lower than "+\
485               "first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
486         print(" ")
487         print("Select another 'end_time_step'")
488         print(" ")
489         exit
490      end if
491      if (end_time_step .GT. t_all(nt-1)/3600)
492         print(" ")
493         print("'end_time_step' = "+ end_time_step +"h is greater than "+\
494               "last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
495         print(" ")
496         print("Select another 'end_time_step'") 
497         print(" ")
498         exit
499      end if
500      if (end_time_step .LT. start_time_step/3600)
501         print(" ")
502         print("'end_time_step' = "+ end_time_step +"h is lower than "+\
503               "'start_time_step' = "+start_time_step+"h")
504         print(" ")
505         print("Select another 'start_time_step' or 'end_time_step'")
506         print(" ")
507         exit
508      end if
509   end if
510   do i=0,nt-1     
511      if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. \
512          end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
513         et=i
514         break
515      else
516         et=0
517      end if
518   end do
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   no_topo= 0
1424
1425   ;****************************************************
1426   ; Preparation of vector plots
1427   ; ****************************************************       
1428
1429   do varn=dim-1,0,1
1430
1431     if (vector .EQ. 1) then   
1432        check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
1433        if (check_vec1) then
1434            temp = f[:]->$vNam(varn)$
1435            data_att = f_att->$vNam(varn)$
1436            vect1=temp(:,zs:ze,ys:ye,xs:xe)
1437            delete(temp)
1438         end if
1439
1440         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
1441         if (check_vec2) then
1442            temp = f[:]->$vNam(varn)$
1443            data_att = f_att->$vNam(varn)$
1444            vect2=temp(:,zs:ze,ys:ye,xs:xe)
1445            delete(temp)
1446         end if
1447   
1448         if (","+vNam(varn)+"," .EQ. vec1) then
1449            v1=v1+1
1450         end if
1451         if (","+vNam(varn)+"," .EQ. vec2) then
1452             v2=v2+1
1453          end if
1454     end if
1455
1456   end do
1457
1458
1459
1460   do varn=dim-1,0,1     
1461           
1462      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1463           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .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            if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1488               dummy=0
1489            else
1490               data_att = f_att->$vNam(varn)$
1491            end if
1492            if(vNam(varn) .eq. "lwps_xy" .or. vNam(varn) .eq. "pras_xy"     \
1493              .or. vNam(varn) .eq. "prrs_xy" .or. vNam(varn) .eq. "qsws_xy" \
1494              .or. vNam(varn) .eq. "shfs_xy" .or. vNam(varn) .eq. "ts_xy"   \
1495              .or. vNam(varn) .eq. "us_xy" .or. vNam(varn) .eq. "z0s_xy"    \
1496              .or. vNam(varn) .eq. "lwp*_xy" .or. vNam(varn) .eq. "pra*_xy" \
1497              .or. vNam(varn) .eq. "prr*_xy" .or. vNam(varn) .eq. "qsws*_xy"\
1498              .or. vNam(varn) .eq. "shf*_xy" .or. vNam(varn) .eq. "t*_xy"   \
1499              .or. vNam(varn) .eq. "u*_xy" .or. vNam(varn) .eq. "z0*_xy") then
1500              ;these variables depend on zu1_xy and that's why they have
1501              ;only one z-layer
1502              data(varn,:,0,:,:)=temp(:,0,ys:ye,xs:xe)
1503              no_zu1=no_zu1+1
1504            else
1505              if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1506                  ;these variables depend on x and y
1507                  data(varn,0,0,:,:)=doubletofloat(temp(ys:ye,xs:xe))
1508                  no_topo=no_topo+1
1509               else
1510                  data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1511               end if
1512            end if
1513            delete(temp)               
1514         end if
1515         if ( xzc .eq. 1 ) then
1516            temp = f[:]->$vNam(varn)$
1517            data_att = f_att->$vNam(varn)$
1518            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1519            delete(temp)
1520         end if
1521         if ( yzc .eq. 1) then
1522            temp = f[:]->$vNam(varn)$
1523            data_att = f_att->$vNam(varn)$
1524            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe)
1525            delete(temp)
1526         end if
1527
1528         data!0 = "var"
1529         data!1 = "t"
1530         data!2 = "z"
1531         data!3 = "y"
1532         data!4 = "x" 
1533
1534         MinVal(varn) = min(data(varn,start_time_step:end_time_step,\
1535                                               0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1536         MaxVal(varn) = max(data(varn,start_time_step:end_time_step,\
1537                                               0:(ze-zs),0:(ye-ys),0:(xe-xs)))
1538         
1539         if (vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi")
1540            unit(varn) = "meters"
1541          else
1542            unit(varn) = data_att@units
1543            delete(data_att)
1544          end if
1545
1546      end if
1547
1548   end do
1549
1550   if (no_var .EQ. 0) then
1551      print(" ")
1552      print("The variables var='"+var+"' do not exist on your input file;")
1553      print("be sure to have one comma before and after each variable")
1554      print(" ")
1555      exit
1556   end if
1557
1558   var_input=new(no_var,string)
1559   numb=0
1560   no_var=0
1561   
1562   do varn=dim-1,0,1   
1563           
1564      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1565           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1566           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1567           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1568           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1569           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1570           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1571           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1572           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1573           vNam(varn) .eq. "ind_x_yz") then
1574         check = False
1575      else
1576         check = True
1577      end if 
1578
1579      if (var .NE. "all") then
1580         check = isStrSubset( var,","+vNam(varn)+"," )
1581      end if
1582
1583      if(check) then     
1584         var_input(numb)=vNam(varn)
1585         numb=numb+1     
1586      end if
1587     
1588      if (check) then
1589         no_var = no_var+1
1590      end if
1591     
1592   end do
1593
1594   if (no_var .EQ. 0) then
1595      print(" ")
1596      print("The variables var='"+var+"' do not exist on your input file;")
1597      print("be sure to have one comma before and after each variable")
1598      print(" ")
1599      exit
1600   end if
1601
1602   if (vector .EQ. 1) then
1603      if (v1 .EQ. 0)then
1604         print(" ")
1605         print("Component 1 for the vector-plot ('vec1') must be one "+\
1606               "of the variables on the input file:")
1607         print("- "+var_input)
1608         print("be sure to have one comma before and after the variable")
1609         print(" ")
1610         exit
1611      end if
1612
1613      if (v2 .EQ. 0)then
1614         print(" ")
1615         print("Component 2 for the vector-plot ('vec2') must be one "+\
1616               "of the variables on the input file:")
1617         print("- "+var_input)
1618         print("be sure to have one comma before and after the variable")
1619         print(" ")
1620         exit
1621      end if
1622   end if
1623
1624   ; ***************************************************
1625   ; open workstation(s)
1626   ; ***************************************************
1627
1628   if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
1629      format_out@wkPaperSize = "A4"
1630   end if
1631   if ( format_out .EQ. "png" ) then
1632      format_out@wkWidth  = 1000
1633      format_out@wkHeight = 1000
1634   end if
1635
1636   wks_ps  = gsn_open_wks(format_out,file_out)
1637   gsn_define_colormap(wks_ps,"rainbow+white")
1638
1639   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1640      plot=new((/no_time*no_layer/),graphic)
1641   else
1642      plot=new((/no_time*no_layer*(no_var-no_topo) + no_topo - \
1643                       no_time*(no_layer-1)*no_zu1/),graphic)
1644   end if
1645   dim_plot=dimsizes(plot)
1646
1647   page = 0
1648   n=0
1649   print(" ")
1650   print("Plots sorted by '"+sort+"'")
1651   print(" ")
1652
1653   ; ***************************************************
1654   ; create plots
1655   ; ***************************************************
1656
1657
1658   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1659      do lo = los, loe                                         
1660         do li = lis, lie
1661            if (xyc .EQ. 1)then
1662               if (sort .EQ. "time")then
1663                  if(z_d(lo) .eq. -1.d) then
1664                    level = "z-average"
1665                  else
1666                    level = "z=" + z_d(lo) + "m"
1667                  end if
1668               else
1669                  if(z_d(li) .eq. -1.d) then
1670                    level = "z-average"
1671                  else
1672                    level = "z=" + z_d(li) + "m"
1673                  end if
1674               end if
1675            end if
1676            if (xzc .EQ. 1)then
1677               if (sort .EQ. "time")then
1678                 if(y_d(lo) .eq. -1.d) then
1679                   level = "y-average"
1680                 else
1681                   level = "y=" + y_d(lo) + "m"
1682                 end if
1683               else
1684                  if(y_d(li) .eq. -1.d) then
1685                    level = "y-average"
1686                  else
1687                    level = "y=" + y_d(li) + "m"
1688                  end if
1689               end if
1690            end if
1691            if (yzc .EQ. 1)then
1692               if (sort .EQ. "time")then
1693                 if(x_d(lo) .eq. -1.d) then
1694                    level = "x-average"
1695                 else
1696                    level = "x=" + x_d(lo) + "m"
1697                 end if
1698               else
1699                 if(x_d(li) .eq. -1.d) then
1700                    level = "x-average"
1701                 else
1702                    level = "x=" + x_d(li) + "m"
1703                 end if
1704               end if
1705            end if               
1706            vecres                  = True            ; vector only resources
1707            vecres@gsnDraw          = False           ; don't draw
1708            vecres@gsnFrame         = False           ; don't advance frame
1709            vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
1710            vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
1711            vecres@vcRefLengthF     = 0.05            ; define length of
1712                                                      ; vec ref
1713            vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
1714            vecres@tmXBLabelFontHeightF   = font_size
1715            vecres@tmYLLabelFontHeightF   = font_size
1716            vecres@tiXAxisFontHeightF     = font_size
1717            vecres@tiYAxisFontHeightF     = font_size
1718            if (sort .EQ. "time")then
1719               vecres@gsnLeftString = "t=" + \
1720                             decimalPlaces(t_all(li)/3600,2,True) +"h   "+level
1721            else
1722               vecres@gsnLeftString = "t=" + \
1723                             decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level
1724            end if
1725            if (xyc .EQ. 1) then 
1726               vecres@tiXAxisString = "x (m)"
1727               vecres@tiYAxisString = "y (m)"   
1728               if (sort .EQ. "time")then
1729                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo-los,:,:),\
1730                                               vect2(li,lo-los,:,:),vecres)
1731               else
1732                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1733                                                vect2(lo,li-lis,:,:),vecres)
1734               end if
1735            end if
1736            if (xzc .EQ. 1) then
1737               vecres@tiXAxisString = "x (m)"
1738               vecres@tiYAxisString = "z (m)"
1739               if (sort .EQ. "time")then
1740                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
1741                                               vect2(li,:,lo-los,:),vecres)
1742               else
1743                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
1744                                                vect2(lo,:,li-lis,:),vecres)
1745               end if
1746            end if
1747            if (yzc .EQ. 1) then
1748               vecres@tiXAxisString = "y (m)"
1749               vecres@tiYAxisString = "z (m)"
1750               if (sort .EQ. "time")then
1751                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
1752                                                vect2(li,:,:,lo-los),vecres)
1753               else
1754                  plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
1755                                                vect2(lo,:,:,li-lis),vecres)
1756               end if
1757            end if
1758            n=n+1
1759         end do
1760      end do
1761   end if
1762 
1763   do varn=dim-1,0,1
1764
1765      if (vector .EQ. 1 ) then   
1766         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",")
1767      end if
1768     
1769      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or.          \
1770           vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or.          \
1771           vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or.     \
1772           vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or.         \
1773           vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or.        \
1774           vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or.   \
1775           vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or.  \
1776           vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. \
1777           vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or.     \
1778           vNam(varn) .eq. "ind_x_yz") then
1779         check = False
1780      else
1781         check = True
1782      end if   
1783 
1784      if (var .NE. "all") then
1785         check = isStrSubset( var,","+vNam(varn)+"," )
1786      end if
1787   
1788      if(check) then
1789
1790         space=(MaxVal(varn)-MinVal(varn))/24
1791 
1792         cs_res@cnMinLevelValF = MinVal(varn)
1793         cs_res@cnMaxLevelValF = MaxVal(varn)
1794
1795         cs_res@cnLevelSpacingF = space
1796     
1797         ; ****************************************************
1798         ; loops over time and layer
1799         ; ****************************************************
1800
1801         
1802         do lo = los, loe                                       
1803            do li = lis, lie
1804
1805               ; ****************************************************
1806               ; xy cross section
1807               ; ****************************************************
1808
1809               if ( xyc .eq. 1 ) then
1810         
1811                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1812                  cs_res@gsnRightString = unit(varn)
1813                 
1814                  if ( sort .eq. "time" ) then
1815                     if ( z_d(lo) .eq. -1)then
1816                        if (delta_z .EQ. -1) then
1817                           level = "-average"
1818                        else
1819                           level = "=" + z_d(lo) + "m"
1820                        end if
1821                     else
1822                        level = "=" + z_d(lo) + "m"
1823                     end if
1824
1825                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1826                        vNam(varn) .eq. "pras_xy"  .or. \
1827                        vNam(varn) .eq. "prrs_xy"  .or. \
1828                        vNam(varn) .eq. "qsws_xy"  .or. \
1829                        vNam(varn) .eq. "shfs_xy"  .or. \
1830                        vNam(varn) .eq. "ts_xy"    .or. \
1831                        vNam(varn) .eq. "us_xy"    .or. \
1832                        vNam(varn) .eq. "z0s_xy"   .or. \
1833                        vNam(varn) .eq. "lwp*_xy"  .or. \
1834                        vNam(varn) .eq. "pra*_xy"  .or. \
1835                        vNam(varn) .eq. "prr*_xy"  .or. \
1836                        vNam(varn) .eq. "qsws*_xy" .or. \
1837                        vNam(varn) .eq. "shf*_xy"  .or. \
1838                        vNam(varn) .eq. "t*_xy"    .or. \
1839                        vNam(varn) .eq. "u*_xy"    .or. \
1840                        vNam(varn) .eq. "z0*_xy") then
1841                         loe = 0
1842                         level = "=" + zu1(0) + "m"
1843                      else
1844                         loe = laye
1845                     end if
1846
1847                     if(vNam(varn) .eq. "zwwi"  .or. \
1848                        vNam(varn) .eq. "zusi") then
1849                         loe = 0
1850                         los = 0
1851                         lie  = 0
1852                         lis = 0
1853                         level = ""
1854                      else
1855                         lis = start_time_step
1856                         lie = end_time_step
1857                         los = lays
1858                         loe = laye
1859                     end if                 
1860
1861                     if(vNam(varn) .eq. "zwwi"  .or. \
1862                        vNam(varn) .eq. "zusi") then
1863                         cs_res@gsnCenterString = ""
1864                     else
1865                         cs_res@gsnCenterString = "t=" + \
1866                          decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level
1867                     end if
1868
1869                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1870                        ;nothing to be done
1871                     else
1872                         plot(n) = gsn_csm_contour(wks_ps,\
1873                                           data(varn,li,lo-los,:,:),cs_res)
1874                     end if
1875
1876
1877                     if (vector .EQ. 1 .AND. check_vecp) then
1878                        vecres                 = True           ; vector only
1879                                                                ; resources
1880                        vecres@gsnDraw         = False          ; don't draw
1881                        vecres@gsnFrame        = False          ; don't
1882                                                                ; advance frame
1883                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly
1884                                                                ; vectors
1885                        vecres@vcRefMagnitudeF = ref_mag        ; define
1886                                                                ; vector ref
1887                                                                ; mag
1888                        vecres@vcRefLengthF    = 0.05           ; define
1889                                                                ; length of
1890                                                                ; vec ref 
1891                        vecres@gsnRightString  = " "
1892                        vecres@gsnLeftString   = " "
1893                        vecres@tiXAxisString   = " "   
1894                        plot_vec=gsn_csm_vector(wks_ps,\
1895                              vect1(li,lo-los,:,:),vect2(li,lo-los,:,:),vecres)
1896                        overlay(plot(n), plot_vec)
1897                     end if                         
1898                  end if
1899                 
1900                  if ( sort .eq. "layer" ) then
1901                     if (z_d(li) .eq. -1) then
1902                        if (delta_z .EQ. -1) then
1903                           level = "-average"
1904                        else
1905                           level = "=" + z_d(li) + "m"
1906                        end if
1907                     else
1908                        level = "=" + z_d(li) + "m"
1909                     end if
1910       
1911                     if(vNam(varn) .eq. "lwps_xy"  .or. \
1912                        vNam(varn) .eq. "pras_xy"  .or. \
1913                        vNam(varn) .eq. "prrs_xy"  .or. \
1914                        vNam(varn) .eq. "qsws_xy"  .or. \
1915                        vNam(varn) .eq. "shfs_xy"  .or. \
1916                        vNam(varn) .eq. "ts_xy"    .or. \
1917                        vNam(varn) .eq. "us_xy"    .or. \
1918                        vNam(varn) .eq. "z0s_xy"   .or. \
1919                        vNam(varn) .eq. "lwp*_xy"  .or. \
1920                        vNam(varn) .eq. "pra*_xy"  .or. \
1921                        vNam(varn) .eq. "prr*_xy"  .or. \
1922                        vNam(varn) .eq. "qsws*_xy" .or. \
1923                        vNam(varn) .eq. "shf*_xy"  .or. \
1924                        vNam(varn) .eq. "t*_xy"    .or. \
1925                        vNam(varn) .eq. "u*_xy"    .or. \
1926                        vNam(varn) .eq. "z0*_xy") then
1927                         lie = 0
1928                         level = "=" + zu1(0) + "m"
1929                     else
1930                         lie = laye
1931                     end if
1932
1933                     if(vNam(varn) .eq. "zwwi"  .or. \
1934                        vNam(varn) .eq. "zusi") then
1935                         lie = 0
1936                         lis = 0
1937                         loe = 0
1938                         los = 0
1939                         level = ""
1940                     else
1941                         lis = lays
1942                         lie = laye
1943                         los = start_time_step
1944                         loe = end_time_step
1945                     end if
1946
1947                     if(vNam(varn) .eq. "zwwi"  .or. \
1948                        vNam(varn) .eq. "zusi") then
1949                         cs_res@gsnCenterString = ""
1950                     else
1951                         cs_res@gsnCenterString = "t=" + \
1952                          decimalPlaces(t_all(lo)/3600,2,True) +"h  z"+level
1953                     end if
1954
1955                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
1956                        ;nothing to be done
1957                     else
1958                        plot(n) = gsn_csm_contour(wks_ps,\
1959                                            data(varn,lo,li-lis,:,:),cs_res)
1960                     end if
1961
1962                     if (vector .EQ. 1 .AND. check_vecp) then
1963                        vecres                 = True           ; vector only
1964                                                                ; resources
1965                        vecres@gsnDraw         = False          ; don't draw
1966                        vecres@gsnFrame        = False          ; don't
1967                                                                ; advance frame
1968                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
1969                        vecres@vcRefMagnitudeF = ref_mag        ; define
1970                                                                ; vector ref
1971                                                                ; mag
1972                        vecres@vcRefLengthF    = 0.05           ; define
1973                                                                ; length of
1974                                                                ; vec ref
1975                        vecres@gsnRightString  = " "            ; turn off
1976                                                                ; right string
1977                        vecres@gsnLeftString   = " "            ; turn off
1978                                                                ; left string
1979                        vecres@tiXAxisString   = " "   
1980                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li-lis,:,:),\
1981                                                   vect2(lo,li-lis,:,:),vecres)
1982                        overlay(plot(n), plot_vec)
1983                     end if
1984                  end if
1985               end if
1986
1987               ; ****************************************************
1988               ; xz cross section
1989               ; ****************************************************
1990
1991               if ( xzc .eq. 1 ) then
1992                 
1993                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
1994               
1995                  if ( sort .eq. "time" ) then
1996                     if ( y_d(lo) .eq. -1 ) then
1997                        level = "-average"
1998                     else
1999                        level = "=" + y_d(lo) + "m"
2000                     end if
2001                     cs_res@gsnCenterString = "t=" + \
2002                           decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
2003
2004                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2005                        ;nothing to be done
2006                     else
2007                         plot(n) = gsn_csm_contour(wks_ps,\
2008                                              data(varn,li,:,lo-los,:),cs_res)
2009                     end if
2010
2011                     if (vector .EQ. 1 .AND. check_vecp) then
2012                        vecres                 = True           ; vector only
2013                                                                ; resources
2014                        vecres@gsnDraw         = False          ; don't draw
2015                        vecres@gsnFrame        = False          ; don't
2016                                                                ; advance frame
2017                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2018                        vecres@vcRefMagnitudeF = ref_mag        ; define
2019                                                                ; vector ref
2020                                                                ; mag
2021                        vecres@vcRefLengthF    = 0.05           ; define
2022                                                                ; length of
2023                                                                ; vec ref
2024                        vecres@gsnRightString  = " "            ; turn off
2025                                                                ; right string
2026                        vecres@gsnLeftString   = " "            ; turn off
2027                                                                ; left string
2028                        vecres@tiXAxisString   = " "   
2029                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo-los,:),\
2030                                                   vect2(li,:,lo-los,:),vecres)
2031                        overlay(plot(n), plot_vec)
2032                     end if
2033                  end if
2034
2035                  if ( sort .eq. "layer" ) then
2036                     if ( y_d(li) .eq. -1 ) then
2037                        level = "-average"
2038                     else
2039                        level = "=" + y_d(li) + "m"
2040                     end if
2041                     cs_res@gsnCenterString = "t=" + \
2042                           decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
2043
2044                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2045                        ;nothing to be done
2046                     else
2047                         plot(n) = gsn_csm_contour(wks_ps,\
2048                                              data(varn,lo,:,li-lis,:),cs_res) 
2049                     end if
2050
2051                     if (vector .EQ. 1 .AND. check_vecp) then
2052                        vecres                 = True           ; vector only
2053                                                                ; resources
2054                        vecres@gsnDraw         = False          ; don't draw
2055                        vecres@gsnFrame        = False          ; don't
2056                                                                ; advance frame
2057                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2058                        vecres@vcRefMagnitudeF = ref_mag        ; define
2059                                                                ; vector ref
2060                                                                ; mag
2061                        vecres@vcRefLengthF    = 0.05           ; define
2062                                                                ; length of
2063                                                                ; vec ref
2064                        vecres@gsnRightString  = " "            ; turn off
2065                                                                ; right string
2066                        vecres@gsnLeftString   = " "            ; turn off
2067                                                                ; left string
2068                        vecres@tiXAxisString   = " "   
2069                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li-lis,:),\
2070                                                   vect2(lo,:,li-lis,:),vecres)
2071                        overlay(plot(n), plot_vec)
2072                     end if
2073                  end if                 
2074               end if
2075
2076               ; ****************************************************
2077               ; yz cross section
2078               ; ****************************************************
2079
2080               if ( yzc .eq. 1 ) then
2081                                 
2082                  cs_res@gsnLeftString = "Plot of "+vNam(varn)
2083
2084                  if ( sort .eq. "time" ) then
2085                     if ( x_d(lo) .eq. -1 ) then
2086                        level = "-average"
2087                     else
2088                        level = "=" + x_d(lo) + "m"
2089                     end if
2090                     cs_res@gsnCenterString = "t=" + \
2091                          decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
2092
2093                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2094                        ;nothing to be done
2095                     else
2096                        plot(n) = gsn_csm_contour(wks_ps,\
2097                                          data(varn,li,:,:,lo-los),cs_res)
2098                     end if
2099
2100                     if (vector .EQ. 1 .AND. check_vecp) then
2101                        vecres                 = True           ; vector only
2102                                                                ; resources
2103                        vecres@gsnDraw         = False          ; don't draw
2104                        vecres@gsnFrame        = False          ; don't
2105                                                                ; advance frame
2106                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2107                        vecres@vcRefMagnitudeF = ref_mag        ; define
2108                                                                ; vector ref
2109                                                                ; mag
2110                        vecres@vcRefLengthF    = 0.05           ; define
2111                                                                ; length of
2112                                                                ; vec ref
2113                        vecres@gsnRightString  = " "            ; turn off
2114                                                                ; right string
2115                        vecres@gsnLeftString   = " "            ; turn off
2116                                                                ; left string
2117                        vecres@tiXAxisString   = " "   
2118                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo-los),\
2119                                                   vect2(li,:,:,lo-los),vecres)
2120                        overlay(plot(n), plot_vec)
2121                     end if
2122                  end if
2123
2124                  if ( sort .eq. "layer" ) then
2125                     if ( x_d(li) .eq. -1 ) then
2126                        level = "-average"
2127                     else
2128                        level = "=" + x_d(li) + "m"
2129                     end if
2130                     cs_res@gsnCenterString = "t=" + \
2131                           decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
2132
2133                     if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2134                        ;nothing to be done
2135                     else
2136                        plot(n) = gsn_csm_contour(wks_ps,\
2137                                          data(varn,lo,:,:,li-lis),cs_res)
2138                     end if
2139
2140                     if (vector .EQ. 1 .AND. check_vecp)then
2141                        vecres                 = True           ; vector only
2142                                                                ;resources
2143                        vecres@gsnDraw         = False          ; don't draw
2144                        vecres@gsnFrame        = False          ; don't
2145                                                                ; advance frame
2146                        vecres@vcGlyphStyle    = "CurlyVector"  ; curly vectors
2147                        vecres@vcRefMagnitudeF = ref_mag        ; define
2148                                                                ; vector ref
2149                                                                ; mag
2150                        vecres@vcRefLengthF    = 0.05           ; define
2151                                                                ; length of
2152                                                                ; vec ref
2153                        vecres@gsnRightString  = " "            ; turn off
2154                                                                ; right string
2155                        vecres@gsnLeftString   = " "            ; turn off
2156                                                                ; left string
2157                        vecres@tiXAxisString   = " "   
2158                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li-lis),\
2159                                                   vect2(lo,:,:,li-lis),vecres)
2160                        overlay(plot(n), plot_vec)
2161                     end if
2162                  end if
2163               end if         
2164               n=n+1 
2165            end do     
2166         end do 
2167      end if
2168   end do
2169
2170   ; ***************************************************
2171   ; merge plots onto one page
2172   ; ***************************************************   
2173
2174   no_frames = 0
2175
2176   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
2177      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
2178           no_time*no_layer*(no_var+1) .gt. no_rows*no_columns) then
2179         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),\
2180                                       (/no_var+1,no_layer*no_time/),cs_resP)
2181         print(" ")
2182         print("Outputs to .eps or .epsi have only one frame")
2183         print(" ")
2184      else
2185         do np = 0,dim_plot-1,no_rows*no_columns   
2186            if ( np + no_rows*no_columns .gt. dim_plot-1) then
2187               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2188                                                                      cs_resP)
2189               no_frames = no_frames + 1
2190            else
2191               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2192                                              (/no_rows,no_columns/),cs_resP)
2193               no_frames = no_frames + 1
2194            end if
2195         end do
2196      end if
2197   else       
2198      if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") \
2199                               .AND. dim_plot .gt. no_rows*no_columns) then
2200         gsn_panel(wks_ps,plot(0:dim_plot-1),(/dim_plot,1/),cs_resP)
2201         print(" ")
2202         print("Outputs to .eps or .epsi have only one frame")
2203         print(" ")
2204      else
2205         do np = 0,dim_plot-1,no_rows*no_columns 
2206            if ( np + no_rows*no_columns .gt. dim_plot-1) then
2207               gsn_panel(wks_ps, plot(np:dim_plot-1),(/no_rows,no_columns/),\
2208                                                                     cs_resP)
2209               no_frames = no_frames + 1
2210            else
2211               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),\
2212                                              (/no_rows,no_columns/),cs_resP)
2213               no_frames = no_frames + 1
2214            end if
2215         end do
2216      end if
2217   end if
2218
2219   if (format_out .EQ. "png" ) then
2220     png_output = new((/no_frames/), string)
2221     j = 0
2222     do i=0, no_frames-1
2223       j = i + 1
2224       if (j .LE. 9) then
2225         png_output(i) = file_out+".00000"+j+".png"
2226       end if
2227       if (j .GT. 9 .AND. j .LE. 99) then
2228         png_output(i) = file_out+".0000"+j+".png"
2229       end if
2230       if (j .GT. 99 .AND. j .LE. 999) then
2231         png_output(i) = file_out+".000"+j+".png"
2232       end if
2233       if (j .GT. 999) then
2234         png_output(i) = file_out+".00"+j+".png"
2235       end if
2236
2237       ;using imagemagick's convert for reducing the white
2238       ;space around the plot
2239       cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
2240              png_output(i) + " " + png_output(i)
2241       system(cmd)
2242     end do
2243
2244     print(" ")
2245     print("Output to: "+ png_output)
2246     print(" ")
2247   else
2248     print(" ")
2249     print("Output to: " + file_out +"."+ format_out)
2250     print(" ")
2251   end if
2252 
2253end
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |