source: palm/trunk/SCRIPTS/NCL/profiles.ncl @ 874

Last change on this file since 874 was 874, checked in by heinze, 12 years ago

Bugfix: correction of x-axis string in case of no_files>1

File size: 171.1 KB
Line 
1load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
2load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
3load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
4
5;***************************************************
6; load .ncl.config or .ncl.config.default
7;***************************************************
8
9function which_script()
10local script
11begin
12   script="profiles"
13   return(script)
14end
15   
16if (isfilepresent("$PALM_BIN/../../.ncl.config")) then
17   loadscript("$PALM_BIN/../../.ncl.config")
18else
19  if (isfilepresent("$PALM_BIN/NCL/.ncl.config.default")) then
20     loadscript( "$PALM_BIN/NCL/.ncl.config.default")
21  else
22      palm_bin_path = getenv("PALM_BIN")
23      print(" ")
24      print("Neither the personal configuration file '.ncl.config' exists in")
25      print("~/palm/current_version")
26      print("nor the default configuration file '.ncl.config.default' "+\
27            "exists in")
28      print(palm_bin_path + "/NCL")
29      print(" ")
30      exit
31   end if
32end if   
33
34   
35begin
36
37   ;***************************************************
38   ; Retrieving the NCL version used
39   ;***************************************************
40   
41   ncl_version_ch = systemfunc("ncl -V")
42   ncl_version    = stringtofloat(ncl_version_ch)
43
44   ;***************************************************
45   ; Retrieving the double quote character
46   ;***************************************************
47   
48   dq=str_get_dq()
49
50   ;***************************************************
51   ; set up default parameter values and strings
52   ;***************************************************
53
54   if (no_files .LT. 1 .OR. no_files .GT. 6) then
55      print(" ")
56      print("Assign 'no_files' between 1 and 6") 
57      print(" ")
58      exit
59   end if
60   
61   file_in   = new(no_files,string)
62   file_in_1 = new(no_files,logical)
63   start_f_t = new(no_files,integer)
64   end_f_t   = new(no_files,integer)
65   
66   if (file_1 .EQ. "File in") then
67      print(" ")
68      print("Declare 1st input file 'file_1=' in '.ncl.config' or prompt")
69      print(" ")
70      exit
71   else
72      file_in(0) = file_1   
73   end if     
74   file_in_1(0) = False
75   if (isStrSubset(file_in(0), ".nc"))then
76      start_f = -2
77      end_f = -2
78      file_in_1(0) = True     
79   end if 
80   if (start_f .EQ. -1)then
81      print(" ")
82      print("'start_f' must be one of the cyclic numbers (at least 0) "+\
83            "of your input file(s)")
84      print(" ") 
85      exit
86   end if
87   if (end_f .EQ. -1)then
88      print(" ")
89      print("'end_f' must be one of the cyclic numbers (at least 0) of "+\
90            "your input file(s)")
91      print(" ") 
92      exit
93   end if           
94   start_f_t(0) = start_f
95   end_f_t(0) = end_f   
96 
97   if (no_files .GT. 1) then
98      if (file_2 .EQ. "File in") then
99         print(" ")
100         print("Declare 2nd input file 'file_2=' in '.ncl.config' or prompt")
101         print(" ")
102         exit
103      else
104         file_in(1) = file_2   
105      end if     
106      file_in_1(1) = False
107      if (isStrSubset(file_in(1), ".nc"))then
108         start_f_2 = -2
109         end_f_2 = -2
110         file_in_1(1) = True     
111      end if
112      if (start_f_2 .EQ. -1)then
113         print(" ")
114         print("'start_f_2' must be one of the cyclic numbers (at least 0) "+\
115              "of your input file(s)")
116         print(" ") 
117         exit
118      end if
119      if (end_f_2 .EQ. -1)then
120         print(" ")
121         print("'end_f_2' must be one of the cyclic numbers (at least 0) "+\
122               "of your input file(s)")
123         print(" ") 
124         exit
125      end if     
126      start_f_t(1) = start_f_2
127      end_f_t(1) = end_f_2   
128   end if   
129   
130   if (no_files .GT. 2) then
131      if (file_3 .EQ. "File in") then
132         print(" ")
133         print("Declare 3rd input file 'file_3=' in '.ncl.config' or prompt")
134         print(" ")
135         exit
136      else
137         file_in(2) = file_3   
138      end if     
139      file_in_1(2) = False
140      if (isStrSubset(file_in(2), ".nc"))then
141         start_f_3 = -2
142         end_f_3 = -2
143         file_in_1(2) = True     
144      end if
145      if (start_f_3 .EQ. -1)then
146         print(" ")
147         print("'start_f_3' must be one of the cyclic numbers (at least 0) "+\
148               "of your input file(s)")
149         print(" ") 
150         exit
151      end if
152      if (end_f_3 .EQ. -1)then
153         print(" ")
154         print("'end_f_3' must be one of the cyclic numbers (at least 0) "+\
155               "of your input file(s)")
156         print(" ") 
157         exit
158      end if     
159      start_f_t(2) = start_f_3
160      end_f_t(2) = end_f_3 
161   end if
162   
163   if (no_files .GT. 3) then
164      if (file_4 .EQ. "File in") then
165         print(" ")
166         print("Declare 4th input file 'file_4=' in '.ncl.config' or prompt")
167         print(" ")
168         exit
169      else
170         file_in(3) = file_4   
171      end if     
172      file_in_1(3) = False
173      if (isStrSubset(file_in(3), ".nc"))then
174         start_f_4 = -2
175         end_f_4 = -2
176         file_in_1(3) = True     
177      end if
178      if (start_f_4 .EQ. -1)then
179         print(" ")
180         print("'start_f_4' must be one of the cyclic numbers (at least 0) "+\
181               "of your input file(s)")
182         print(" ") 
183         exit
184      end if
185      if (end_f_4 .EQ. -1)then
186         print(" ")
187         print("'end_f_4' must be one of the cyclic numbers (at least 0) "+\
188               "of your input file(s)")
189         print(" ") 
190         exit
191      end if     
192      start_f_t(3) = start_f_4
193      end_f_t(3) = end_f_4
194   end if
195   
196   if (no_files .GT. 4) then
197      if (file_5 .EQ. "File in") then
198         print(" ")
199         print("Declare 5th input file 'file_5=' in '.ncl.config' or prompt")
200         print(" ")
201         exit
202      else
203         file_in(4) = file_5   
204      end if     
205      file_in_1(4) = False
206      if (isStrSubset(file_in(4), ".nc"))then
207         start_f_5 = -2
208         end_f_5 = -2
209         file_in_1(4) = True     
210      end if
211      if (start_f_5 .EQ. -1)then
212         print(" ")
213         print("'start_f_5' must be one of the cyclic numbers (at least 0) "+\
214               "of your input file(s)")
215         print(" ") 
216         exit
217      end if
218      if (end_f_5 .EQ. -1)then
219         print(" ")
220         print("'end_f_5' must be one of the cyclic numbers (at least 0) "+\
221               "of your input file(s)")
222         print(" ") 
223         exit
224      end if     
225      start_f_t(4) = start_f_5
226      end_f_t(4) = end_f_5   
227   end if
228   
229   if (no_files .GT. 5) then
230      if (file_6 .EQ. "File in") then
231         print(" ")
232         print("Declare 6th input file 'file_6=' in '.ncl.config' or prompt")
233         print(" ")
234         exit
235      else
236         file_in(5) = file_6   
237      end if     
238      file_in_1(5) = False
239      if (isStrSubset(file_in(5), ".nc"))then
240         start_f_6 = -2
241         end_f_6 = -2
242         file_in_1(5) = True     
243      end if
244      if (start_f_6 .EQ. -1)then
245         print(" ")
246         print("'start_f_6' must be one of the cyclic numbers (at least 0) "+\
247               "of your input file(s)")
248         print(" ") 
249         exit
250      end if
251      if (end_f_6 .EQ. -1)then
252         print(" ")
253         print("'end_f_6' must be one of the cyclic numbers (at least 0) "+\
254               "of your input file(s)")
255         print(" ") 
256         exit
257      end if     
258      start_f_t(5) = start_f_6
259      end_f_t(5) = end_f_6   
260   end if
261
262   
263   if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND.  \
264      format_out .NE. "eps" .AND. format_out .NE. "ps" .AND.    \
265      format_out .NE. "epsi" .AND. format_out .NE. "ncgm" .AND. \
266      format_out .NE. "png")then
267      print(" ")
268      print("'format_out = "+format_out+"' is invalid and set to'x11'")
269      print(" ")
270      format_out="x11"
271   end if 
272
273   if (ncl_version .LE. 5.1 .AND. format_out .EQ. "png") then
274      print(" ")
275      print("Output of png files not available")
276      print("png output is avaiable with NCL version 5.2.0 and higher ")
277      print("NCL version used: " + ncl_version_ch)
278      print(" ")
279      exit
280   end if
281
282   if (over .NE. 0 .AND. over .NE. 1) then
283      print(" ")
284      print("'over'= "+over+" is invalid and set to 0")
285      print(" ")
286      over = 0
287   end if
288   if (no_files .GT. 1) then
289      over = 0
290      print(" ")
291      print("If you have more than one input file - you cannot overlay "+\
292            "variables: over is set to 0")
293      print(" ")
294   end if
295   
296   if (combine .NE. 0 .AND. combine .NE. 1)then
297      print(" ")
298      print("Your 'combine'= "+combine+" is invalid and set to 0")
299      print(" ")
300      combine = 0
301   end if
302   if (no_files .GT. 1) then
303      combine = 0
304      print(" ")
305      print("If you have more than one input file you cannot combine "+\
306            "variables: combine is set to 0")
307      print(" ")
308   end if   
309   if (combine .EQ. 1 .AND. number_comb .EQ. -1) then
310      print(" ")
311      print("Set 'number_comb' to 2 or 3 or combine to 0")
312      print(" ")
313      exit
314   end if
315   if (combine .EQ. 1)then
316      if (number_comb .EQ. -1) then
317         print(" ")
318         print("Set 'number_comb' to 2 or 3 or combine to 0")
319         print(" ")
320         exit
321      end if
322      if (number_comb .NE. 2 .AND. number_comb .NE. 3) then
323         print(" ")
324         print("Set 'number_comb' to 2 or 3 or combine to 0")
325         print(" ")
326         exit
327      end if
328   end if 
329   if (combine .EQ. 1 .AND. c_var .EQ. "c_variables") then
330      print(" ")
331      print("Select variables for overlaying ('c_var') or set combine to 0")
332      print(" ")
333      exit
334   end if
335
336   if (log_z .NE. 0 .AND. log_z .NE. 1)then
337      print(" ")
338      print("'log_z'= "+log_z+" is invalid and set to 0")
339      print(" ")
340      log_z = 0
341   end if   
342 
343   if (norm_z .EQ. 0) then
344      print(" ")
345      print("Normalising with 0 is not allowed, 'norm_z' is set to 1.0")
346      print(" ")
347      norm_z = 1.0
348   end if
349
350   ;***************************************************
351   ; open input file
352   ;***************************************************
353
354   do nof=0,no_files-1
355
356   files=new(end_f_t(nof)-start_f_t(nof)+1,string)
357   if (file_in_1(nof))then
358      if (isfilepresent(file_in(nof)))then
359         files(0)=file_in(nof)
360      else
361         print(" ")
362         print("1st input file: '"+file_in(nof)+"' does not exist")
363         print(" ")
364         exit
365      end if
366   else
367      if (start_f_t(nof) .EQ. 0)then
368         if (isfilepresent(file_in(nof)+".nc"))then
369            files(0)=file_in(nof)+".nc"
370            do i=1,end_f_t(nof)
371               if (isfilepresent(file_in(nof)+"."+i+".nc"))then   
372                  files(i)=file_in(nof)+"."+i+".nc"
373               else
374                  print(" ")
375                  print("Input file: '"+file_in(nof)+"."+i+".nc' does not "+\
376                        "exist")
377                  print(" ")
378                  exit 
379               end if       
380            end do         
381         else
382            print(" ")
383            print("Input file: '"+file_in(nof)+".nc' does not exist")
384            print(" ")
385            exit
386         end if
387      else
388         do i=start_f_t(nof),end_f_t(nof)
389            if (isfilepresent(file_in(nof)+"."+i+".nc"))then   
390               files(i-start_f_t(nof))=file_in(nof)+"."+i+".nc"
391            else
392               print(" ")
393               print("Input file: '"+file_in(nof)+"."+i+".nc' does not exist")
394               print(" ")
395               exit 
396            end if
397         end do
398      end if
399   end if
400
401   f=addfiles(files,"r")
402   f_att=addfile(files(0),"r")
403   ListSetType(f,"cat")
404   
405   vNam  = getfilevarnames(f_att)
406   vType = getfilevartypes(f_att,vNam)
407
408   if ((all(vType .eq. "double"))) then ;distinction if data is double or float
409      check_vType = True
410   else
411      check_vType = False
412   end if
413   
414   if (nof .EQ. 0)then
415      vNam0=vNam
416   end if
417   if (nof .NE. 0)then
418      if (dim0 .NE. dim)then
419         print(" ")
420         print("There are 'no_files'="+no_files+" input files but they do "+\
421               "not contain the same variables")
422         print(" ")
423         exit
424      else
425         do i=0,dim0-1
426            if (vNam0(i) .NE. vNam(i))then
427               print(" ")
428               print("There are 'no_files'="+no_files+" input files but "+\
429                     "they do not contain the same variables")
430               print(" ")
431               exit
432            end if
433         end do
434      end if
435   end if
436   nof=nof+1
437   print(" ")
438   print("Variables in input file "+nof+":")
439   print("- "+ vNam)
440   print(" ")
441   nof=nof-1
442   dim = dimsizes(vNam)
443   dim0=dim
444 
445   if (dim .EQ. 0) then
446      print(" ")
447      print("There is no data on file")
448      print(" ")
449   end if
450   
451   prof3d = 0
452   do varn = dim-1,0,1
453      if ( isStrSubset( vNam(varn), "zu_3d") .OR. \
454           isStrSubset( vNam(varn), "zw_3d")) then
455         prof3d = 1
456         break
457      end if
458   end do
459 
460   if (var .NE. "all") then
461      ;rearrange the order of the variables in vNam so that the variables
462      ;specified by "var" are in the top of vNam
463
464      vNam_static = new((/dim/),string)
465      vNam_temp   = new((/dim/),string)
466
467      var_char    = stringtocharacter(var)
468      no_char     = dimsizes(var_char)
469      comma       = 0
470       
471      do j=0,no_char-1
472         if(var_char(j) .eq. ",")
473            comma = comma + 1
474         end if
475      end do
476
477      if(comma .le. 1)
478          print(" ")
479          print("The variables 'var="+var+"'" )
480          print("do not exist on your input file;")
481          print("be sure to have one comma before and after each variable")
482          print(" ")
483          exit
484      end if
485
486      if(comma .gt. dim)
487          print(" ")
488          print("The variables 'var="+var+"'" )
489          print("do not exist on your input file;")
490          print("be sure to have one comma before and after each variable")
491          print(" ")
492          exit
493      end if
494
495      indices = new((/comma/),integer)
496      comma   = 0
497       
498      do j=0,no_char-1
499         if(var_char(j) .eq. ",")
500           indices(comma) = j
501           comma          = comma + 1
502         end if
503      end do
504
505      do j=0,comma-2
506         vNam_temp(j) = charactertostring(\
507                                   var_char(indices(j)+1:indices(j+1)-1))
508      end do
509
510      do j=0,comma-2
511        count_check = 0
512        do varn=0,dim-1
513           if(vNam_temp(j) .ne. vNam(varn))
514             count_check=count_check+1
515           end if
516        end do   
517
518        if(count_check .eq. dim)
519            print(" ")
520            print("The variables 'var="+var+"'" )
521            print("do not exist on your input file;")
522            print("be sure to have one comma before and after each variable")
523            print(" ")
524            exit
525        end if
526
527        vNam_static(j) = vNam_temp(comma - 2 - j)
528      end do
529
530      vNam_temp   = vNam_static
531      vNam_static = ""
532
533      if (prof3d .EQ. 0) then
534         i = 0
535         do j=0,dim-1,2
536            vNam_static(j)   = vNam_temp(i)
537            vNam_static(j+1) = "z" + vNam_static(j)
538            i = i+1
539            if(i .eq. comma-1)
540                break
541            end if
542         end do
543
544         delete(vNam_temp)
545         count_variable = 2*(comma-1)
546      else
547         i = 0
548          do j=0,dim-1,1
549             vNam_static(j) = vNam_temp(i)
550             i = i+1
551             if(i .eq. comma-1)
552                 break
553             end if
554          end do
555
556          delete(vNam_temp)
557          count_variable = comma-1
558      end if
559
560      counter  = 0
561      counter2 = 0
562      do j=0,dim-1
563         counter = 0
564         do i = 0, count_variable - 1       
565            if(vNam_static(i) .ne. vNam(j))
566               counter = counter + 1
567            end if
568         end do
569     
570         if (counter .eq. count_variable)
571            vNam_static(count_variable + counter2) = vNam(j)
572            counter2 = counter2 + 1
573         end if
574      end do
575     
576      vNam = vNam_static
577      delete(vNam_static)
578    end if
579
580   ;---------below steps only for first file -> nof=0
581   if (nof .EQ. 0) then
582
583      plot = new(dim,graphic)
584      plot_ = new(dim,graphic)
585      if (no_files .GT. 1) then
586         multi_plot = new((/no_files,dim/),graphic)
587        if (check_vType) then
588          max_nof = new((/no_files,dim/),double)
589          min_nof = new((/no_files,dim/),double)
590        else
591          max_nof = new((/no_files,dim/),float)
592          min_nof = new((/no_files,dim/),float)
593        end if
594         name    = new((/no_files,dim/),string)
595         unit_   = new((/no_files,dim/),string)
596      end if
597
598   if (prof3d .EQ. 0) then
599 
600   do varn=0,dim-1
601      if ( isStrSubset( vNam(varn), "time") .OR. \
602           isStrSubset( vNam(varn), "NORM")) then
603         continue
604      end if
605      if (vNam(varn) .EQ. "u" .OR. isStrSubset(vNam(varn), "u_"))then
606         z_u = f_att->$vNam(varn+1)$
607         break
608      else
609         if (vNam(varn) .EQ. "v" .OR. isStrSubset(vNam(varn), "v_"))then 
610            z_u = f_att->$vNam(varn+1)$
611            break
612         else
613            if(vNam(varn) .EQ. "pt" .OR. isStrSubset(vNam(varn), "pt_"))then
614               z_u = f_att->$vNam(varn+1)$
615               break
616            else
617               if(vNam(varn) .EQ. "vpt" .OR. isStrSubset(vNam(varn), "vpt_"))then
618                  z_u = f_att->$vNam(varn+1)$
619                  break
620               else
621                  if(vNam(varn) .EQ. "lpt" .OR. isStrSubset(vNam(varn), "lpt_"))then
622                     z_u = f_att->$vNam(varn+1)$
623                     break
624                  else   
625                     if(vNam(varn) .EQ. "q" .OR. isStrSubset(vNam(varn), "q_") )then
626                        z_u = f_att->$vNam(varn+1)$
627                        break
628                     else
629                        if(vNam(varn) .EQ. "qv" .OR. isStrSubset(vNam(varn), "qv_"))then
630                           z_u = f_att->$vNam(varn+1)$
631                           break
632                        else
633                           if(vNam(varn) .EQ. "ql" .OR. isStrSubset(vNam(varn), "ql_"))then
634                              z_u = f_att->$vNam(varn+1)$
635                              break
636                           else
637                              if(vNam(varn) .EQ. "rho" .OR. isStrSubset(vNam(varn), "rho_"))then
638                                 z_u = f_att->$vNam(varn+1)$
639                                 break
640                              else
641                                 if(vNam(varn) .EQ. "s" .OR. isStrSubset(vNam(varn), "s_") )then
642                                    z_u = f_att->$vNam(varn+1)$
643                                    break
644                                 else
645                                    if(vNam(varn) .EQ. "sa" .OR. isStrSubset(vNam(varn), "sa_"))then
646                                       z_u = f_att->$vNam(varn+1)$
647                                       break
648                                    else
649                                       if(vNam(varn) .EQ. "e" .OR. isStrSubset(vNam(varn), "e_"))then
650                                         
651                                          z_u = f_att->$vNam(varn+1)$
652                                          break
653                                       else
654                                          if(vNam(varn) .EQ. "es" .OR. \
655                                           isStrSubset(vNam(varn), "es_") .OR. vNam(varn) .EQ. "e*" .OR. isStrSubset(vNam(varn), "e*_"))then
656                                             z_u = f_att->$vNam(varn+1)$
657                                             break
658                                          else
659                                             if(vNam(varn) .EQ. "km" .OR. isStrSubset(vNam(varn), "km_"))then
660                                                z_u = f_att->$vNam(varn+1)$
661                                                break
662                                             else
663                                                if(vNam(varn) .EQ. "kh" .OR. isStrSubset(vNam(varn), "kh_"))then
664                                                   z_u = f_att->$vNam(varn+1)$
665                                                   break
666                                                else
667                                                   if(vNam(varn) .EQ. "l" .OR. isStrSubset(vNam(varn), "l_"))then
668                                                      z_u = f_att->$vNam(varn+1)$
669                                                      break
670                                                   else
671                                                      if(vNam(varn) .EQ. "us2" .OR. isStrSubset(vNam(varn), "us2_")\
672                                                         .OR. vNam(varn) .EQ. "u*2" .OR. isStrSubset(vNam(varn), "u*2_"))then
673                                                         z_u = f_att->$vNam(varn+1)$
674                                                         break
675                                                      else
676                                                         if(vNam(varn) .EQ. "vs2" .OR. isStrSubset(vNam(varn), "vs2_")\
677                                                            .OR. vNam(varn) .EQ. "v*2" .OR. isStrSubset(vNam(varn), "v*2_"))then
678                                                            z_u = f_att->$vNam(varn+1)$
679                                                            break
680                                                         else
681                                                            if(vNam(varn) .EQ. "pts2" .OR. isStrSubset(vNam(varn), "pts2_")\
682                                                               .OR. vNam(varn) .EQ. "pt*2" .OR. isStrSubset(vNam(varn), "pt*2_"))then
683                                                               z_u = f_att->$vNam(varn+1)$
684                                                               break
685                                                            else
686                                                               if(vNam(varn) .EQ. "wsususodz" .OR. isStrSubset(vNam(varn), "wsususodz_")\
687                                                                  .OR. vNam(varn) .EQ. "w*u*u*:dz" .OR. isStrSubset(vNam(varn), "w*u*u*:dz_"))then
688                                                                  z_u = f_att->$vNam(varn+1)$
689                                                                  break
690                                                               else
691                                                                  if(vNam(varn) .EQ. "wspsodz" .OR. isStrSubset(vNam(varn), "wspsodz_")\
692                                                                     .OR. vNam(varn) .EQ. "w*p*:dz" .OR. isStrSubset(vNam(varn), "w*p*:dz_"))then
693                                                                     z_u = f_att->$vNam(varn+1)$
694                                                                     break
695                                                                  else
696                                                                     if(vNam(varn) .EQ. "wpeodz" .OR. isStrSubset(vNam(varn), "wpeodz_")\
697                                                                        .OR. vNam(varn) .EQ. "w"+dq+"e:dz" .OR. isStrSubset(vNam(varn), "w"+dq+"e:dz_"))then
698                                                                        z_u = f_att->$vNam(varn+1)$
699                                                                        break   
700                                                                     else
701                                                                        if(vNam(varn) .EQ. "qs2" .OR. isStrSubset(vNam(varn), "qs2_")\
702                                                                           .OR. vNam(varn) .EQ. "q*2" .OR. isStrSubset(vNam(varn), "q*2_"))then
703                                                                           z_u = f_att->$vNam(varn+1)$
704                                                                           break 
705                                                                       end if               
706                                                                     end if
707                                                                  end if
708                                                               end if
709                                                            end if
710                                                         end if
711                                                      end if
712                                                   end if
713                                                end if
714                                             end if
715                                          end if
716                                       end if
717                                    end if
718                                 end if
719                              end if
720                           end if
721                        end if
722                     end if
723                  end if
724               end if
725            end if
726         end if
727      end if
728   end do
729
730   do varn=0,dim-1
731      if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "NORM")) then
732         continue
733      end if
734      if (vNam(varn) .EQ. "w" .OR. isStrSubset(vNam(varn), "w_"))then
735         z_w = f_att->$vNam(varn+1)$
736         break
737      else
738         if (vNam(varn) .EQ. "wpup" .OR. isStrSubset(vNam(varn), "wpup_")\
739            .OR. vNam(varn) .EQ. "w"+dq+"u"+dq .OR. isStrSubset(vNam(varn), "w"+dq+"u"+dq+"_"))then
740            z_w = f_att->$vNam(varn+1)$
741            break
742         else
743            if(vNam(varn) .EQ. "wsus" .OR. isStrSubset(vNam(varn), "wsus_")\
744               .OR. vNam(varn) .EQ. "w*u*" .OR. isStrSubset(vNam(varn), "w*u*_"))then
745               z_w = f_att->$vNam(varn+1)$
746               break
747            else
748               if(vNam(varn) .EQ. "wu" .OR. isStrSubset(vNam(varn), "wu_"))then
749                  z_w = f_att->$vNam(varn+1)$
750                  break
751               else
752                  if(vNam(varn) .EQ. "wpvp" .OR. isStrSubset(vNam(varn), "wpvp_")\
753                     .OR. vNam(varn) .EQ. "w"+dq+"v"+dq .OR. isStrSubset(vNam(varn), "w"+dq+"v"+dq+"_"))then
754                     z_w = f_att->$vNam(varn+1)$
755                     break
756                  else   
757                     if(vNam(varn) .EQ. "wsvs" .OR. isStrSubset(vNam(varn), "wsvs_")\
758                        .OR. vNam(varn) .EQ. "w*v*" .OR. isStrSubset(vNam(varn), "w*v*_"))then
759                        z_w = f_att->$vNam(varn+1)$
760                        break
761                     else
762                        if(vNam(varn) .EQ. "wv" .OR. isStrSubset(vNam(varn), "wv_"))then
763                           z_w = f_att->$vNam(varn+1)$
764                           break
765                        else
766                           if(vNam(varn) .EQ. "wptpp" .OR. isStrSubset(vNam(varn), "wptpp_")\
767                              .OR. vNam(varn) .EQ. "w"+dq+"pt"+dq .OR. isStrSubset(vNam(varn), "w"+dq+"pt"+dq+"_"))then
768                              z_w = f_att->$vNam(varn+1)$
769                              break
770                           else
771                              if(vNam(varn) .EQ. "wspts" .OR. isStrSubset(vNam(varn), "wspts_")\
772                                 .OR. vNam(varn) .EQ. "w*pt*" .OR. isStrSubset(vNam(varn), "w*pt*_"))then
773                                 z_w = f_att->$vNam(varn+1)$
774                                 break
775                              else
776                                 if(vNam(varn) .EQ. "wpt" .OR. isStrSubset(vNam(varn), "wpt_"))then
777                                    z_w = f_att->$vNam(varn+1)$
778                                    break
779                                 else
780                                    if(vNam(varn) .EQ. "wsptsBC" .OR. isStrSubset(vNam(varn), "wsptsBC_")\
781                                       .OR. vNam(varn) .EQ. "w*pt*BC" .OR. isStrSubset(vNam(varn), "w*pt*BC_"))then
782                                       z_w = f_att->$vNam(varn+1)$
783                                       break
784                                    else
785                                       if(vNam(varn) .EQ. "wptBC" .OR. isStrSubset(vNam(varn), "wptBC_"))then
786                                          z_w = f_att->$vNam(varn+1)$
787                                          break
788                                       else
789                                          if(vNam(varn) .EQ. "wpvptp" .OR. isStrSubset(vNam(varn), "wpvptp_")\
790                                              .OR. vNam(varn) .EQ. "w"+dq+"vpt"+dq .OR. isStrSubset(vNam(varn), "w"+dq+"vpt"+dq+"_"))then
791                                             z_w = f_att->$vNam(varn+1)$
792                                             break
793                                          else
794                                             if(vNam(varn) .EQ. "wsvpts" .OR. isStrSubset(vNam(varn), "wsvpts_")\
795                                               .OR. vNam(varn) .EQ. "w*vpt*" .OR. isStrSubset(vNam(varn), "w*vpt*_"))then
796                                                z_w = f_att->$vNam(varn+1)$
797                                                break
798                                             else
799                                                if(vNam(varn) .EQ. "wvpt" .OR. isStrSubset(vNam(varn), "wvpt_"))then
800                                                   z_w = f_att->$vNam(varn+1)$
801                                                   break
802                                                else
803                                                   if(vNam(varn) .EQ. "wpqp" .OR. isStrSubset(vNam(varn), "wpqp_")\
804                                                      .OR. vNam(varn) .EQ. "w"+dq+"q"+dq .OR. isStrSubset(vNam(varn), "w"+dq+"q"+dq+"_"))then
805                                                      z_w = f_att->$vNam(varn+1)$
806                                                      break
807                                                   else
808                                                      if(vNam(varn) .EQ. "wsqs" .OR. isStrSubset(vNam(varn), "wsqs_")\
809                                                         .OR. vNam(varn) .EQ. "w*q*" .OR. isStrSubset(vNam(varn), "w*q*_"))then
810                                                         z_w = f_att->$vNam(varn+1)$
811                                                         break
812                                                      else
813                                                         if(vNam(varn) .EQ. "wq" .OR. isStrSubset(vNam(varn), "wq_"))then
814                                                            z_w = f_att->$vNam(varn+1)$
815                                                            break
816                                                         else
817                                                            if(vNam(varn) .EQ. "wpqvp" .OR. isStrSubset(vNam(varn), "wpqvp_")\
818                                                               .OR. vNam(varn) .EQ. "w"+dq+"qv"+dq .OR. isStrSubset(vNam(varn), "w"+dq+"qv"+dq+"_"))then
819                                                               z_w = f_att->$vNam(varn+1)$
820                                                               break
821                                                            else
822                                                               if(vNam(varn) .EQ. "wsqvs" .OR. isStrSubset(vNam(varn), "wsqvs_")\
823                                                                  .OR. vNam(varn) .EQ. "w*qv*" .OR. isStrSubset(vNam(varn), "w*qv*_"))then
824                                                                  z_w = f_att->$vNam(varn+1)$
825                                                                  break
826                                                               else
827                                                                  if(vNam(varn) .EQ. "wqv" .OR. isStrSubset(vNam(varn), "wqv_"))then
828                                                                     z_w = f_att->$vNam(varn+1)$
829                                                                     break
830                                                                  else
831                                                                     if(vNam(varn) .EQ. "wpsp" .OR. isStrSubset(vNam(varn), "wpsp_")\
832                                                                        .OR. vNam(varn) .EQ. "w"+dq+"s"+dq .OR. isStrSubset(vNam(varn), "w"+dq+"s"+dq+"_"))then
833                                                                        z_w = f_att->$vNam(varn+1)$
834                                                                        break
835                                                                     else
836                                                                        if(vNam(varn) .EQ. "wsss" .OR. isStrSubset(vNam(varn), "wsss_")\
837                                                                            .OR. vNam(varn) .EQ. "w*s*" .OR. isStrSubset(vNam(varn), "w*s*_"))then
838                                                                           z_w = f_att->$vNam(varn+1)$
839                                                                           break
840                                                                        else
841                                                                           if(vNam(varn) .EQ. "ws" .OR. isStrSubset(vNam(varn), "ws_"))then
842                                                                              z_w = f_att->$vNam(varn+1)$
843                                                                              break
844                                                                           else
845                                                                              if(vNam(varn) .EQ. "wpsap" .OR. isStrSubset(vNam(varn), "wpsap_")\
846                                                                                  .OR. vNam(varn) .EQ. "w"+dq+"sa"+dq .OR. isStrSubset(vNam(varn), "w"+dq+"sa"+dq+"_"))then
847                                                                                 z_w = f_att->$vNam(varn+1)$
848                                                                                 break
849                                                                              else
850                                                                                 if(vNam(varn) .EQ. "wssas" .OR. isStrSubset(vNam(varn), "wssas_")\
851                                                                                     .OR. vNam(varn) .EQ. "w*sa*" .OR. isStrSubset(vNam(varn), "w*sa*_"))then
852                                                                                    z_w = f_att->$vNam(varn+1)$
853                                                                                    break
854                                                                                 else
855                                                                                    if(vNam(varn) .EQ. "wsa" .OR. isStrSubset(vNam(varn), "wsa_"))then
856                                                                                       z_w = f_att->$vNam(varn+1)$
857                                                                                       break
858                                                                                    else
859                                                                                       if(vNam(varn) .EQ. "wses" .OR. isStrSubset(vNam(varn), "wses_")\
860                                                                                          .OR. vNam(varn) .EQ. "w*e*" .OR. isStrSubset(vNam(varn), "w*e*_"))then
861                                                                                          z_w = f_att->$vNam(varn+1)$
862                                                                                          break
863                                                                                       else
864                                                                                          if(vNam(varn) .EQ. "ws2" .OR. isStrSubset(vNam(varn), "ws2_")\
865                                                                                             .OR. vNam(varn) .EQ. "w*2" .OR. isStrSubset(vNam(varn), "w*2_"))then
866                                                                                             z_w = f_att->$vNam(varn+1)$
867                                                                                             break
868                                                                                          else
869                                                                                             if(vNam(varn) .EQ. "ws3" .OR. isStrSubset(vNam(varn), "ws3_")\
870                                                                                                .OR. vNam(varn) .EQ. "w*3" .OR. isStrSubset(vNam(varn), "w*3_"))then
871                                                                                                z_w = f_att->$vNam(varn+1)$
872                                                                                                break
873                                                                                             else
874                                                                                                if(vNam(varn) .EQ. "Sw" .OR. isStrSubset(vNam(varn), "Sw_"))then
875                                                                                                   z_w = f_att->$vNam(varn+1)$
876                                                                                                   break
877                                                                                                else
878                                                                                                   if(vNam(varn) .EQ. "ws2pts".OR. isStrSubset(vNam(varn), "ws2pts_")\
879                                                                                                      .OR. vNam(varn) .EQ. "w*2pt*" .OR. isStrSubset(vNam(varn), "w*2pt*_"))then
880                                                                                                      z_w = f_att->$vNam(varn+1)$
881                                                                                                      break
882                                                                                                   else
883                                                                                                      if(vNam(varn) .EQ. "wspts2" .OR. isStrSubset(vNam(varn), "wspts2_")\
884                                                                                                         .OR. vNam(varn) .EQ. "w*pt*2" .OR. isStrSubset(vNam(varn), "w*pt*2_"))then
885                                                                                                         z_w = f_att->$vNam(varn+1)$
886                                                                                                         break                                           
887                                                                                                      end if
888                                                                                                   end if
889                                                                                                end if
890                                                                                             end if
891                                                                                          end if
892                                                                                       end if
893                                                                                    end if
894                                                                                 end if
895                                                                              end if
896                                                                           end if
897                                                                        end if   
898                                                                     end if
899                                                                  end if
900                                                               end if
901                                                            end if
902                                                         end if
903                                                      end if
904                                                   end if
905                                                end if
906                                             end if
907                                          end if
908                                       end if
909                                    end if
910                                 end if
911                              end if
912                           end if
913                        end if
914                     end if
915                  end if
916               end if
917            end if
918         end if
919      end if
920   end do
921   
922   if (.not. isvar("z_u") .AND. .not. isvar ("z_w")) then
923      co=0
924      do varn=0,dim-1     
925         if ( isStrSubset( vNam(varn), "time") .OR. \
926              isStrSubset( vNam(varn), "NORM")) then
927            check = False
928         else
929            check = True
930            if (var .NE. "all") then
931               check = isStrSubset( var,","+vNam(varn)+"," )
932            end if         
933         end if
934         if (check)then
935            co=co+1
936            z = f_att->$vNam(varn+1)$
937            if (getvardims(z) .EQ. "zu")then
938               z_u = z
939            else
940               if (getvardims(z) .EQ. "zw")then
941                  z_w = z
942               end if
943            end if
944            dimz = dimsizes(z)
945            break
946         end if   
947      end do
948      if (co .EQ. 0) then
949         print(" ")
950         print("The variables 'var="+var+"'" )
951         print("do not exist on your input file;")
952         print("be sure to have one comma before and after each variable")
953         print(" ")
954         exit           
955      end if
956   end if
957
958   if (isvar("z_u") ) then
959      dimz  = dimsizes(z_u)
960   else
961      if (isvar("z_w"))then
962         dimz  = dimsizes(z_w)
963      end if
964   end if
965
966   else
967
968      do varn = dim-1,0,1
969         if (vNam(varn) .EQ. "zu_3d")then
970            z_u = f_att->zu_3d 
971            dimz  = dimsizes(z_u)         
972         else
973            if (vNam(varn) .EQ. "zw_3d")then
974               z_w = f_att->zw_3d
975               dimz  = dimsizes(z_w)
976            end if
977         end if
978         if (vNam(varn) .EQ. "x")then
979            x = f_att->x
980            dimx=dimsizes(x)
981         else
982            if (vNam(varn) .EQ. "xu")then
983               x = f_att->xu
984               dimx=dimsizes(x)
985            end if
986         end if
987         if (vNam(varn) .EQ. "y")then
988            y = f_att->y
989            dimy=dimsizes(y)
990         else
991            if (vNam(varn) .EQ. "yv")then
992               y = f_att->yv
993               dimy=dimsizes(y)
994            end if
995         end if
996      end do
997
998   end if
999   
1000   if(.not. isvar("z_u") .AND. .not. isvar("z_w"))then
1001      print(" ")
1002      print("Program aborts - there are no z-variables available")
1003      print("Be sure if 'plot_3d' is set correctly")
1004      print(" ")
1005      exit
1006   end if
1007 
1008
1009   t_all = f[:]->time
1010   nt    = dimsizes(t_all)
1011
1012   if (nt .EQ. 1)then
1013      delta_t=t_all(nt-1)/nt
1014   else
1015      delta_t_array = new(nt-1,double)
1016
1017      do i=0,nt-2
1018         delta_t_array(i) = t_all(i+1)-t_all(i)
1019      end do
1020
1021      delta_t = min(delta_t_array)
1022      delete(delta_t_array)
1023   end if
1024
1025
1026   ; ****************************************************       
1027   ; start of time step and different types of mistakes that could be done
1028   ; ****************************************************
1029   
1030   if (start_time_step .EQ. -1.) then           
1031      start_time_step=t_all(0)/3600     
1032   else
1033      if (start_time_step .GT. t_all(nt-1)/3600)then
1034         print(" ")
1035         print("'start_time_step' = "+ start_time_step +"h is greater "+\
1036               "than last time step = " \
1037               + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
1038         print(" ")
1039         print("Select another 'start_time_step'")
1040         print(" ")
1041         exit
1042      end if
1043      if (start_time_step .LT. t_all(0)/3600)then
1044         print(" ")
1045         print("'start_time_step' = "+ start_time_step +"h is lower "+\
1046               "than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
1047         print(" ")
1048         exit
1049      end if
1050   end if
1051
1052   do i=0,nt-1   
1053      if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND.\
1054          start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
1055         st=i
1056         break
1057      else
1058         st=0
1059      end if
1060   end do
1061   
1062   ; ****************************************************
1063   ; end of time step and different types of mistakes that could be done
1064   ; ****************************************************
1065
1066   if (end_time_step .EQ. -1.) then             
1067      end_time_step = t_all(nt-1)/3600 
1068   else
1069      if (end_time_step .GT. t_all(nt-1)/3600)then
1070         print(" ")
1071         print("'end_time_step' = "+ end_time_step +"h is greater "+\
1072               "than last time step = " +\
1073                t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
1074         print(" ")
1075         print("Select another 'end_time_step'") 
1076         print(" ")
1077         exit
1078      end if
1079      if (end_time_step .LT. start_time_step/3600)then
1080         print(" ")
1081         print("'end_time_step' = "+ end_time_step +"h is lower "+\
1082               "than 'start_time_step' = "+start_time_step+"h")
1083         print(" ")
1084         print("Select another 'start_time_step' or 'end_time_step'")
1085         print(" ")
1086         exit
1087      end if
1088   end if
1089
1090   do i=0,nt-1     
1091      if (end_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND.\
1092          end_time_step .LT. (t_all(i)+delta_t/2)/3600)then
1093         et=i
1094         break
1095       else
1096         et=0
1097      end if
1098   end do
1099 
1100   delete(start_time_step)
1101   start_time_step=round(st,3)
1102   delete(end_time_step)
1103   end_time_step=round(et,3)
1104
1105   no_time = end_time_step-start_time_step+1
1106   
1107   ; ****************************************************
1108   ; time_stride and different types of mistakes that could be done
1109   ; ****************************************************
1110 
1111   if (time_stride .LT. 1) then
1112      print(" ")
1113      print("'time_stride' has to be positive and is set to 1")
1114      print(" ")
1115      time_stride = 1
1116   end if
1117
1118   if (time_stride .GT. no_time) then
1119      print(" ")
1120      print("'time_stride' is greater than number of available "+\
1121           "time steps,")
1122      print("=> 'time_stride' is set to 1")
1123      time_stride = 1
1124   end if
1125
1126   ti_in = ispan(start_time_step,end_time_step,time_stride) ;ti_in contents
1127                                                            ;the time indices
1128                                                            ;to plot
1129   np    = dimsizes(ti_in) 
1130
1131   print(" ")
1132   print("Output of time steps from "+t_all(start_time_step)/3600+\
1133         " h = "+t_all(start_time_step)+" s => index = "+start_time_step)
1134   print("                     till "+t_all(ti_in(np-1))/3600+" h = "\
1135        +t_all(ti_in(np-1))+" s => index = "+end_time_step)
1136   print("                     with temporal stride = "+time_stride)
1137   print(" ")
1138
1139
1140   ; ****************************************************
1141   ; set up legend and colors
1142   ; ****************************************************
1143   
1144   legend_label=new(np,string)
1145   do p=0, np-1
1146      legend_label(p)=sprintf("%6.2f", t_all(ti_in(p))/3600)
1147   end do
1148
1149   ; ***************************************************
1150   ; set up recourses
1151   ; ***************************************************
1152
1153   res                         = True
1154   res@gsnDraw                 = False
1155   res@gsnFrame                = False
1156   res@txFont                  = "helvetica"
1157   res@tiMainFont              = "helvetica"
1158   res@tiXAxisFont             = "helvetica"
1159   res@tiYAxisFont             = "helvetica"
1160   res@tmXBLabelFont           = "helvetica"
1161   res@tmYLLabelFont           = "helvetica"
1162   res@lgLabelFont             = "helvetica"
1163   res@tmLabelAutoStride       = True
1164   if (legend .EQ. 1)then
1165      res@pmLegendDisplayMode  = "Always"
1166   end if
1167
1168   res@pmLegendSide            = "Top"
1169   res@xyExplicitLegendLabels  = legend_label
1170   res@pmLegendParallelPosF    = 1.15
1171   res@pmLegendOrthogonalPosF  = -1.0
1172   res@pmLegendWidthF          = 0.12
1173   res@pmLegendHeightF         = 0.05*np
1174   res@lgLabelFontHeightF      = font_size_legend
1175   res@lgTitleString           = "Time (h)"
1176   res@lgTitleFontHeightF      = font_size   
1177   res@txFontHeightF           = font_size
1178   res@tiXAxisFontHeightF      = font_size
1179   res@tiYAxisFontHeightF      = font_size
1180   res@tmXBLabelFontHeightF    = font_size
1181   res@tmYLLabelFontHeightF    = font_size
1182   res@tiXAxisString            = " "
1183   if ( black .eq. 0 ) then 
1184      res@xyLineColors = -(ispan(-237,-2,235/np))
1185   end if
1186   if (norm_z .EQ. 1)then
1187      res@tiYAxisString = "Height (m)"
1188   else   
1189      res@tiYAxisString = "Height / "+norm_z+" (m)"
1190   end if
1191   
1192   if (log_z .EQ. 1) then
1193      res@trYLog = True
1194   end if
1195
1196   if (dash .EQ. 0 ) then
1197      res@xyMonoDashPattern = True
1198   else
1199      res@xyMonoDashPattern = False
1200      if (no_files .GT. 1)
1201         res@xyMonoDashPattern = True 
1202         print(" ")
1203         print("If you use more than one file, patterns for different "+\
1204               "timesteps cannot be used")
1205         print(" ")
1206      end if       
1207   end if
1208
1209   res@tmXBMinorPerMajor = 4
1210   res@tmYLMinorPerMajor = 4
1211
1212   resP                        = True
1213   resP@txFont                 = "helvetica"
1214   resP@txString               = f_att@title
1215   resP@txFuncCode             = "~"
1216   resP@txFontHeightF          = 0.015
1217
1218   ; ***************************************************
1219   ; set up graphics for plot
1220   ; ***************************************************
1221
1222   if (combine .EQ. 1) then
1223      plot_o = new(number_comb,graphic)   
1224      label=new(number_comb,string)
1225      color_o=new(number_comb,integer)
1226
1227      if (check_vType) then
1228        mini=new(number_comb,double)
1229        maxi=new(number_comb,double)
1230      else
1231        mini=new(number_comb,float)
1232        maxi=new(number_comb,float)
1233      end if
1234   end if
1235 
1236   if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
1237      format_out@wkPaperSize = "A4"
1238   end if
1239   if ( format_out .EQ. "png" ) then
1240      format_out@wkWidth  = 1000
1241      format_out@wkHeight = 1000
1242   end if
1243
1244   wks=gsn_open_wks(format_out,file_out)
1245   gsn_define_colormap(wks,"rainbow+white")
1246
1247   ; ***************************************************
1248   ; set up minimum and maximum height
1249   ; ***************************************************
1250
1251   if (log_z .EQ. 1)then
1252      if (min_z .EQ. -1)then
1253         if (isvar("z_u"))then
1254            min_z=z_u(1)
1255         else
1256            min_z=z_w(1)
1257         end if       
1258      else
1259         if (isvar("z_w"))then
1260            if (min_z .GE. max(z_w) ) then
1261               print(" ")
1262               print("Minimum of height ('min_z'="+min_z+") is greater "+\
1263                     "than available heights (="+max(z_w)+")")
1264               print(" ")
1265               exit
1266            end if         
1267         else
1268            if (min_z .GE. max(z_u) ) then
1269               print(" ")
1270               print("Minimum of height ('min_z'="+min_z+") is greater "+\
1271                     "than available heights (="+max(z_u)+")")
1272               print(" ")
1273               exit
1274            end if
1275         end if
1276         if (isvar("z_u"))then   
1277            if (min_z .LT. z_u(1) ) then
1278               print(" ")
1279               print("Begin height 'min_z' at least at level k=1 (="+\
1280                     z_u(1)+"m) due to the logarithmic scale of the y-axis")
1281               print(" ")
1282               exit
1283            end if
1284         else
1285            if (min_z .LT. z_w(1) ) then
1286               print(" ")
1287               print("Begin height 'min_z' at least at level k=1 (="+\
1288                     z_w(1)+"m) due to the logarithmic scale of the y-axis")
1289               print(" ")
1290               exit
1291            end if   
1292         end if
1293      end if
1294   else
1295      if (isvar("z_u"))then
1296         if (min_z .EQ. -1)then
1297            min_z=z_u(0)
1298         end if
1299      else
1300         if (min_z .EQ. -1)then
1301            min_z=z_w(0)
1302         end if   
1303      end if
1304   
1305      if (isvar("z_w"))then
1306         if (min_z .GE. max(z_w) ) then
1307            print(" ")
1308            print("Minimum of height ('min_z'="+min_z+") is greater than "+\
1309                  "available heights (="+max(z_w)+")")
1310            print(" ")
1311            exit
1312         end if
1313      else
1314         if (min_z .EQ. -1)then
1315            min_z=z_u(0)
1316         end if
1317         if (min_z .GE. max(z_u) ) then
1318            print(" ")
1319            print("Minimum of height ('min_z'="+min_z+") is greater than "+\
1320                  "available heights (="+max(z_u)+")")
1321            print(" ")
1322            exit
1323         end if
1324      end if
1325   end if
1326   
1327   if (isvar("z_w"))then
1328      if (max_z .EQ. -1)then
1329         max_z=max(z_w)
1330      end if
1331   else
1332      if (max_z .EQ. -1)then
1333         max_z=max(z_u)
1334      end if   
1335   end if
1336   
1337   if (isvar("z_w"))then
1338      if (max_z .GT. max(z_w) ) then
1339         print(" ")
1340         print("Maximum of height ('max_z'="+max_z+") is greater than "+\
1341               "available heights (="+max(z_w)+")")
1342         print(" ")
1343         exit
1344      end if
1345   end if
1346
1347   min_z=min_z/norm_z
1348   max_z=max_z/norm_z
1349
1350   ; ***************************************************
1351   ; read data and create plots
1352   ; ***************************************************
1353     
1354   do ti =0,np-1
1355      if( t_all(ti_in(ti)) .lt. 10^36) then
1356         start_time_step = ti
1357         break
1358      end if
1359   end do 
1360   
1361   if (log_z .EQ. 1) then     
1362      if (check_vType) then
1363         data   = new((/dim,np,dimz-1/),double)
1364         data_0 = new((/np,dimz-1/),double)
1365      else
1366         data   = new((/dim,np,dimz-1/),float)
1367         data_0 = new((/np,dimz-1/),float)
1368      end if
1369      data@_FillValue=9.96921e+36
1370      data_0 = 0.1
1371      t      = new((/np,dimz-1/),float)
1372      t      = 0.0
1373      unit   = new(dim,string)
1374      if (isvar("z_u"))then
1375         if (typeof(z_u) .EQ. "double")then
1376            z_v    = new((/dim,dimz/),double)
1377            z_     = new((/dim,dimz-1/),double)
1378         else
1379            if (typeof(z_u) .EQ. "float")then
1380               z_v    = new((/dim,dimz/),float)
1381               z_     = new((/dim,dimz-1/),float)
1382            end if
1383         end if
1384      else
1385         if (isvar("z_w"))then
1386            if (typeof(z_w) .EQ. "double")then
1387               z_v    = new((/dim,dimz/),double)
1388               z_     = new((/dim,dimz-1/),double)
1389            else
1390               if (typeof(z_w) .EQ. "float")then
1391                  z_v    = new((/dim,dimz/),float)
1392                  z_     = new((/dim,dimz-1/),float)
1393               end if
1394            end if
1395         end if
1396      end if
1397   else
1398      if (check_vType) then
1399         data   = new((/dim,np,dimz/),double)
1400         data_0 = new((/np,dimz/),double)
1401      else
1402         data   = new((/dim,np,dimz/),float)
1403         data_0 = new((/np,dimz/),float)
1404      end if
1405      data@_FillValue=9.96921e+36     
1406      data_0 = 0.0
1407      t      = new((/np,dimz/),float)
1408      t      = 0.0
1409      unit   = new(dim,string)
1410      if (isvar("z_u"))then
1411         if (typeof(z_u) .EQ. "double")then
1412            z_v    = new((/dim,dimz/),double)
1413            z_     = new((/dim,dimz/),double)
1414         else
1415            if (typeof(z_u) .EQ. "float")then
1416               z_v    = new((/dim,dimz/),float)
1417               z_     = new((/dim,dimz/),float)
1418            end if
1419         end if
1420      else
1421         if (isvar("z_w"))then
1422            if (typeof(z_w) .EQ. "double")then
1423               z_v    = new((/dim,dimz/),double)
1424               z_     = new((/dim,dimz/),double)
1425            else
1426               if (typeof(z_w) .EQ. "float")then
1427                  z_v    = new((/dim,dimz/),float)
1428                  z_     = new((/dim,dimz/),float)
1429               end if
1430            end if
1431         end if
1432      end if
1433   end if
1434
1435   end if
1436   ;-------above steps only for first file
1437
1438   ; ***************************************************
1439   ; indicate plot number
1440   ; ***************************************************
1441   
1442   if (combine .EQ. 1) then
1443      n = 1
1444   else
1445      n = 0
1446   end if
1447
1448   if (over .EQ. 1) then
1449      plot_u         = gsn_csm_xy(wks,t,data_0(:,:),res)
1450      plot_v         = gsn_csm_xy(wks,t,data_0(:,:),res)
1451      plot_w         = gsn_csm_xy(wks,t,data_0(:,:),res)
1452      plot_pt        = gsn_csm_xy(wks,t,data_0(:,:),res)     
1453      plot_vpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
1454      plot_lpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
1455      plot_q         = gsn_csm_xy(wks,t,data_0(:,:),res)
1456      plot_qv        = gsn_csm_xy(wks,t,data_0(:,:),res)
1457      plot_ql        = gsn_csm_xy(wks,t,data_0(:,:),res)
1458      plot_rho       = gsn_csm_xy(wks,t,data_0(:,:),res)
1459      plot_s         = gsn_csm_xy(wks,t,data_0(:,:),res)
1460      plot_sa        = gsn_csm_xy(wks,t,data_0(:,:),res)
1461      plot_e         = gsn_csm_xy(wks,t,data_0(:,:),res)
1462      plot_es        = gsn_csm_xy(wks,t,data_0(:,:),res)
1463      plot_km        = gsn_csm_xy(wks,t,data_0(:,:),res)
1464      plot_kh        = gsn_csm_xy(wks,t,data_0(:,:),res)
1465      plot_l         = gsn_csm_xy(wks,t,data_0(:,:),res)     
1466      plot_wpup      = gsn_csm_xy(wks,t,data_0(:,:),res)
1467      plot_wsus      = gsn_csm_xy(wks,t,data_0(:,:),res)
1468      plot_wu        = gsn_csm_xy(wks,t,data_0(:,:),res)
1469      plot_wpvp      = gsn_csm_xy(wks,t,data_0(:,:),res)
1470      plot_wsvs      = gsn_csm_xy(wks,t,data_0(:,:),res)
1471      plot_wv        = gsn_csm_xy(wks,t,data_0(:,:),res)
1472      plot_wpptp     = gsn_csm_xy(wks,t,data_0(:,:),res)
1473      plot_wspts     = gsn_csm_xy(wks,t,data_0(:,:),res)
1474      plot_wpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
1475      plot_wsptsBC   = gsn_csm_xy(wks,t,data_0(:,:),res)
1476      plot_wptBC     = gsn_csm_xy(wks,t,data_0(:,:),res)
1477      plot_wpvptp    = gsn_csm_xy(wks,t,data_0(:,:),res)
1478      plot_wsvpts    = gsn_csm_xy(wks,t,data_0(:,:),res)
1479      plot_wvpt      = gsn_csm_xy(wks,t,data_0(:,:),res)
1480      plot_wpqp      = gsn_csm_xy(wks,t,data_0(:,:),res)
1481      plot_wsqs      = gsn_csm_xy(wks,t,data_0(:,:),res)
1482      plot_wq        = gsn_csm_xy(wks,t,data_0(:,:),res)
1483      plot_wpqvp     = gsn_csm_xy(wks,t,data_0(:,:),res)
1484      plot_wsqvs     = gsn_csm_xy(wks,t,data_0(:,:),res)
1485      plot_wqv       = gsn_csm_xy(wks,t,data_0(:,:),res)
1486      plot_wpsp      = gsn_csm_xy(wks,t,data_0(:,:),res)
1487      plot_wsss      = gsn_csm_xy(wks,t,data_0(:,:),res)
1488      plot_ws        = gsn_csm_xy(wks,t,data_0(:,:),res)
1489      plot_wpsap     = gsn_csm_xy(wks,t,data_0(:,:),res)
1490      plot_wssas     = gsn_csm_xy(wks,t,data_0(:,:),res)
1491      plot_wsa       = gsn_csm_xy(wks,t,data_0(:,:),res)
1492      plot_wses      = gsn_csm_xy(wks,t,data_0(:,:),res)
1493      plot_us2       = gsn_csm_xy(wks,t,data_0(:,:),res)
1494      plot_vs2       = gsn_csm_xy(wks,t,data_0(:,:),res)
1495      plot_ws2       = gsn_csm_xy(wks,t,data_0(:,:),res)
1496      plot_pts2      = gsn_csm_xy(wks,t,data_0(:,:),res)
1497      plot_ws3       = gsn_csm_xy(wks,t,data_0(:,:),res)
1498      plot_Sw        = gsn_csm_xy(wks,t,data_0(:,:),res)
1499      plot_ws2pts    = gsn_csm_xy(wks,t,data_0(:,:),res)
1500      plot_wspts2    = gsn_csm_xy(wks,t,data_0(:,:),res)
1501      plot_wsususodz = gsn_csm_xy(wks,t,data_0(:,:),res)
1502      plot_wspsodz   = gsn_csm_xy(wks,t,data_0(:,:),res)
1503      plot_wpeodz    = gsn_csm_xy(wks,t,data_0(:,:),res)
1504 
1505      if (check_vType) then
1506        miniu         =  100000.d
1507        maxiu         = -100000.d
1508        miniv         =  100000.d
1509        maxiv         = -100000.d
1510        miniw         =  100000.d
1511        maxiw         = -100000.d
1512        minipt        =  100000.d
1513        maxipt        = -100000.d
1514        minivpt       =  100000.d
1515        maxivpt       = -100000.d
1516        minilpt       =  100000.d
1517        maxilpt       = -100000.d
1518        miniq         =  100000.d
1519        maxiq         = -100000.d
1520        miniqv        =  100000.d
1521        maxiqv        = -100000.d
1522        miniql        =  100000.d
1523        maxiql        = -100000.d
1524        minie         =  100000.d
1525        maxie         = -100000.d
1526        minies        =  100000.d
1527        maxies        = -100000.d
1528        minikm        =  100000.d
1529        maxikm        = -100000.d
1530        minikh        =  100000.d
1531        maxikh        = -100000.d
1532        miniwpup      =  100000.d
1533        maxiwpup      = -100000.d
1534        miniwsus      =  100000.d
1535        maxiwsus      = -100000.d
1536        miniwu        =  100000.d
1537        maxiwu        = -100000.d
1538        miniwpvp      =  100000.d
1539        maxiwpvp      = -100000.d
1540        miniwsvs      =  100000.d
1541        maxiwsvs      = -100000.d
1542        miniwv        =  100000.d
1543        maxiwv        = -100000.d
1544        miniwpptp     =  100000.d
1545        maxiwpptp     = -100000.d
1546        miniwspts     =  100000.d
1547        maxiwspts     = -100000.d
1548        miniwpt       =  100000.d
1549        maxiwpt       = -100000.d
1550        miniwsptsBC   =  100000.d
1551        maxiwsptsBC   = -100000.d
1552        miniwptBC     =  100000.d
1553        maxiwptBC     = -100000.d
1554        miniwpvptp    =  100000.d
1555        maxiwpvptp    = -100000.d
1556        miniwsvpts    =  100000.d
1557        maxiwsvpts    = -100000.d
1558        miniwvpt      =  100000.d
1559        maxiwvpt      = -100000.d
1560        miniwpqp      =  100000.d
1561        maxiwpqp      = -100000.d
1562        miniwsqs      =  100000.d
1563        maxiwsqs      = -100000.d
1564        miniwq        =  100000.d
1565        maxiwq        = -100000.d
1566        miniwpqvp     =  100000.d
1567        maxiwpqvp     = -100000.d
1568        miniwsqvs     =  100000.d
1569        maxiwsqvs     = -100000.d
1570        miniwqv       =  100000.d
1571        maxiwqv       = -100000.d
1572        miniwpsp      =  100000.d
1573        maxiwpsp      = -100000.d
1574        miniwsss      =  100000.d
1575        maxiwsss      = -100000.d
1576        miniws        =  100000.d
1577        maxiws        = -100000.d
1578        miniwpsap     =  100000.d
1579        maxiwpsap     = -100000.d
1580        miniwssas     =  100000.d
1581        maxiwssas     = -100000.d
1582        miniwsa       =  100000.d
1583        maxiwsa       = -100000.d
1584        minius2       =  100000.d
1585        maxius2       = -100000.d
1586        minivs2       =  100000.d
1587        maxivs2       = -100000.d
1588        miniws2       =  100000.d
1589        maxiws2       = -100000.d
1590        miniwsususodz =  100000.d
1591        maxiwsususodz = -100000.d
1592        miniwspsodz   =  100000.d
1593        maxiwspsodz   = -100000.d
1594        miniwpeodz    =  100000.d
1595        maxiwpeodz    = -100000.d
1596      else
1597        miniu         =  100000.
1598        maxiu         = -100000.
1599        miniv         =  100000.
1600        maxiv         = -100000.
1601        miniw         =  100000.
1602        maxiw         = -100000.
1603        minipt        =  100000.
1604        maxipt        = -100000.
1605        minivpt       =  100000.
1606        maxivpt       = -100000.
1607        minilpt       =  100000.
1608        maxilpt       = -100000.
1609        miniq         =  100000.
1610        maxiq         = -100000.
1611        miniqv        =  100000.
1612        maxiqv        = -100000.
1613        miniql        =  100000.
1614        maxiql        = -100000.
1615        minie         =  100000.
1616        maxie         = -100000.
1617        minies        =  100000.
1618        maxies        = -100000.
1619        minikm        =  100000.
1620        maxikm        = -100000.
1621        minikh        =  100000.
1622        maxikh        = -100000.
1623        miniwpup      =  100000.
1624        maxiwpup      = -100000.
1625        miniwsus      =  100000.
1626        maxiwsus      = -100000.
1627        miniwu        =  100000.
1628        maxiwu        = -100000.
1629        miniwpvp      =  100000.
1630        maxiwpvp      = -100000.
1631        miniwsvs      =  100000.
1632        maxiwsvs      = -100000.
1633        miniwv        =  100000.
1634        maxiwv        = -100000.
1635        miniwpptp     =  100000.
1636        maxiwpptp     = -100000.
1637        miniwspts     =  100000.
1638        maxiwspts     = -100000.
1639        miniwpt       =  100000.
1640        maxiwpt       = -100000.
1641        miniwsptsBC   =  100000.
1642        maxiwsptsBC   = -100000.
1643        miniwptBC     =  100000.
1644        maxiwptBC     = -100000.
1645        miniwpvptp    =  100000.
1646        maxiwpvptp    = -100000.
1647        miniwsvpts    =  100000.
1648        maxiwsvpts    = -100000.
1649        miniwvpt      =  100000.
1650        maxiwvpt      = -100000.
1651        miniwpqp      =  100000.
1652        maxiwpqp      = -100000.
1653        miniwsqs      =  100000.
1654        maxiwsqs      = -100000.
1655        miniwq        =  100000.
1656        maxiwq        = -100000.
1657        miniwpqvp     =  100000.
1658        maxiwpqvp     = -100000.
1659        miniwsqvs     =  100000.
1660        maxiwsqvs     = -100000.
1661        miniwqv       =  100000.
1662        maxiwqv       = -100000.
1663        miniwpsp      =  100000.
1664        maxiwpsp      = -100000.
1665        miniwsss      =  100000.
1666        maxiwsss      = -100000.
1667        miniws        =  100000.
1668        maxiws        = -100000.
1669        miniwpsap     =  100000.
1670        maxiwpsap     = -100000.
1671        miniwssas     =  100000.
1672        maxiwssas     = -100000.
1673        miniwsa       =  100000.
1674        maxiwsa       = -100000.
1675        minius2       =  100000.
1676        maxius2       = -100000.
1677        minivs2       =  100000.
1678        maxivs2       = -100000.
1679        miniws2       =  100000.
1680        maxiws2       = -100000.
1681        miniwsususodz =  100000.
1682        maxiwsususodz = -100000.
1683        miniwspsodz   =  100000.
1684        maxiwspsodz   = -100000.
1685        miniwpeodz    =  100000.
1686        maxiwpeodz    = -100000.
1687      end if
1688
1689   end if
1690
1691   if (prof3d .EQ. 1)then
1692
1693   if (end_x .EQ. -1) then
1694      end_x=dimx-2
1695   end if
1696   if (end_y .EQ. -1)then
1697      end_y=dimy-2
1698   end if
1699   if (start_x .LT. 0)then
1700      print(" ")
1701      print("'start_x' is lower than 0 and set to 0")
1702      print(" ")
1703      start_x=0
1704   end if
1705   if (start_x .GT. dimx-1)then
1706      print(" ")
1707      print("'start_x' is greater than available x-range and set to "+\
1708            "maximum of x-range (excluding ghostpoint)")
1709      print(" ")
1710      start_x=dimx-2
1711   end if
1712   if (end_x .EQ. dimx-1)then
1713      print(" ")
1714      print("'end_x' = "+end_x+" and includes the ghostpoint")
1715      print(" ")
1716   end if
1717   if (end_x .GT. dimx-1)then
1718      print(" ")
1719      print("'end_x' = "+end_x+" is greater than available x-range and set "+\
1720            "to maximum of x-range (excluding ghostpoint)")
1721      print(" ")
1722      end_x=dimx-2
1723   end if
1724   if (end_x .LT. start_x)then
1725      print(" ")
1726      print("'end_x' = "+end_x+" is lower than 'start_x' = "+start_x+\
1727            " and set to maximum of x-range (excluding ghostpoint)")
1728      print(" ")
1729      end_x=dimx-2
1730   end if
1731   if (start_y .LT. 0)then
1732      print(" ")
1733      print("'start_y' is lower than 0 and set to 0")
1734      print(" ")
1735      start_y=0
1736   end if
1737   if (start_y .GT. dimy-1)then
1738      print(" ")
1739      print("'start_y' is greater than available y-range and set to "+\
1740            "maximum of y-range (excluding ghostpoint)")
1741      print(" ")
1742      start_x=dimy-2
1743   end if
1744   if (end_y .EQ. dimy-1)then
1745      print(" ")
1746      print("'end_y' = "+end_y+" and includes the ghostpoint")
1747      print(" ")
1748   end if
1749   if (end_y .GT. dimy-1)then
1750      print(" ")
1751      print("'end_y' = "+end_y+" is greater than available y-range and "+\
1752            "set to maximum of y-range (excluding ghostpoint)")
1753      print(" ")
1754      end_x=dimy-2
1755   end if
1756   if (end_y .LT. start_y)then
1757      print(" ")
1758      print("'end_y' = "+end_y+" is lower than 'start_y' = "+start_y+\
1759            " and set to maximum of y-range (excluding ghostpoint)")
1760      print(" ")
1761      end_y=dimy-2
1762   end if
1763   
1764   end if
1765 
1766   n_o=0
1767   count_var=0
1768
1769   res@xyDashPattern = 1*nof
1770
1771   over_remind = False
1772   if ( over .eq. 1)then
1773     over_remind = True
1774   end if
1775   
1776   do varn = 0,dim-1
1777
1778      check = True
1779     
1780      if (prof3d .EQ. 0) then
1781         if ( isStrSubset( vNam(varn), "time") .OR. \
1782              isStrSubset( vNam(varn), "NORM")) then
1783            check = False
1784         end if
1785      else
1786         if ( isStrSubset( vNam(varn), "time") .OR.  \
1787              isStrSubset( vNam(varn), "zusi") .OR.  \
1788              isStrSubset( vNam(varn), "zwwi") .OR.  \
1789              isStrSubset( vNam(varn), "x") .OR.     \
1790              isStrSubset( vNam(varn), "xu") .OR.    \
1791              isStrSubset( vNam(varn), "y") .OR.     \
1792              isStrSubset( vNam(varn), "yv") .OR.    \
1793              isStrSubset( vNam(varn), "zu_3d") .OR. \
1794              isStrSubset( vNam(varn), "zw_3d")) then
1795            check = False
1796         end if
1797      end if
1798
1799      if (var .NE. "all") then
1800         check = isStrSubset( var,","+vNam(varn)+"," )
1801      end if
1802
1803      if (combine .EQ. 1) then         
1804         com=isStrSubset(c_var,","+vNam(varn)+"," )     
1805         if (com) then     
1806            if (prof3d .EQ. 0) then                 
1807               temp     = f[:]->$vNam(varn)$         
1808               temp_att = f_att->$vNam(varn)$
1809               if (log_z .EQ. 1) then
1810                  do j=0,np-1
1811                     data(varn,j,:) = temp(ti_in(j),1:dimz-1)
1812                  end do
1813               else
1814                  do j=0,np-1
1815                     data(varn,j,:) = temp(ti_in(j),0:dimz-1)
1816                  end do
1817               end if
1818            else
1819               if (log_z .EQ. 1) then
1820                  do i=1,dimz-1
1821                     do j=0,np-1
1822                        temp = f[:]->$vNam(varn)$
1823                        temp_att = f_att->$vNam(varn)$
1824                        data_temp = temp(ti_in(j),i,\
1825                                         start_y:end_y,start_x:end_x)
1826                        data(varn,j,i-1) = dim_avg_Wrap(\
1827                                                   dim_avg_Wrap(data_temp))
1828                     end do
1829                  end do
1830               else
1831                  do i=0,dimz-1
1832                     do j=0,np-1
1833                        temp = f[:]->$vNam(varn)$
1834                        temp_att = f_att->$vNam(varn)$
1835                        data_temp = temp(ti_in(j),i,\
1836                                            start_y:end_y,start_x:end_x)
1837                        data(varn,j,i) = dim_avg_Wrap(dim_avg_Wrap(data_temp))
1838                     end do
1839                  end do
1840               end if
1841               print(" ")
1842               print("Variable for combine '"+vNam(varn)+"' is read")
1843               print(" ")
1844            end if                 
1845            unit(varn) = temp_att@units
1846            if (n_o .GT. number_comb-1) then
1847               print(" ")
1848               print("Set 'number_comb' to the number of overlaying "+\
1849                     "variables ('c_var' = "+c_var+")")
1850               print(" ")
1851               exit
1852            end if
1853            if (abs(min(data(varn,:,:))) .GT. 10)then
1854               min_value = abs(0.01*min(data(varn,:,:)))
1855               max_value = abs(0.01*max(data(varn,:,:)))
1856            else
1857               if (abs(min(data(varn,:,:))) .LT. 0.01 .AND. \
1858                   abs(max(data(varn,:,:))) .GT. 0.01)then
1859                  min_value = abs(0.1*max(data(varn,:,:)))
1860                  max_value = abs(0.1*max(data(varn,:,:)))
1861               else
1862                  if (abs(max(data(varn,:,:))) .LT. 0.01 .AND. \
1863                      abs(min(data(varn,:,:))) .GT. 0.01)then
1864                     min_value = abs(0.1*min(data(varn,:,:)))
1865                     max_value = abs(0.1*min(data(varn,:,:)))
1866                  else
1867                     min_value = abs(0.1*min(data(varn,:,:)))
1868                     max_value = abs(0.1*max(data(varn,:,:)))
1869                  end if
1870               end if
1871            end if
1872            if (min(data(varn,:,:)) .EQ. 0 .AND. \
1873                max(data(varn,:,:)) .EQ. 0)then
1874               min_value = 0.1
1875               max_value = 0.1
1876            end if
1877            mini(n_o)=min(data(varn,:,:))
1878            maxi(n_o)=max(data(varn,:,:))
1879            n_o=n_o+1
1880         end if
1881      end if
1882
1883      if(check) then
1884         
1885         count_var=count_var+1
1886
1887         if (prof3d .EQ. 0) then         
1888            temp = f[:]->$vNam(varn)$
1889            temp_att = f_att->$vNam(varn)$
1890         else
1891             if (log_z .EQ. 1) then
1892               do i=1,dimz-1
1893                  do j=0,np-1
1894                     temp= f[:]->$vNam(varn)$
1895                     temp_att = f_att->$vNam(varn)$
1896                     data_temp = temp(ti_in(j),i,start_y:end_y,start_x:end_x)
1897                     data(varn,j,i-1) = dim_avg_Wrap(dim_avg_Wrap(data_temp))
1898                     delete(data_temp)
1899                  end do
1900               end do
1901            else
1902               do i=0,dimz-1
1903                  do j=0,np-1
1904                     temp= f[:]->$vNam(varn)$
1905                     temp_att = f_att->$vNam(varn)$
1906                     data_temp = temp(ti_in(j),i,start_y:end_y,start_x:end_x)
1907                     data_temp!0 = "t"
1908                     data_temp!1 = "z"
1909                     data(varn,j,i) = dim_avg_Wrap(dim_avg_Wrap(data_temp))
1910                     delete(data_temp)
1911                  end do
1912               end do
1913            end if
1914            print(" ")
1915            print("Variable '"+vNam(varn)+"' is read")
1916            print(" ")
1917            unit(varn) = temp_att@units
1918            a=getvardims(temp_att)
1919            b=dimsizes(a)
1920         end if
1921
1922         if (prof3d .EQ. 0) then
1923            if (log_z .EQ. 1) then
1924               z = f_att->$vNam(varn+1)$(1:dimz-1)
1925               unit(varn) = temp_att@units
1926               do j=0,np-1
1927                  data(varn,j,:) = temp(ti_in(j),1:dimz-1)
1928               end do
1929            else
1930               z = f_att->$vNam(varn+1)$
1931               unit(varn) = temp_att@units
1932               do j=0,np-1
1933                  data(varn,j,:) = temp(ti_in(j),:)
1934               end do
1935            end if
1936         else
1937            do i=0,b-1           
1938               if (isStrSubset( a(i),"zu_3d" ))then
1939                  z_v(varn,:) = z_u
1940                  if (log_z .EQ. 1) then
1941                     z = z_v(varn,1:dimz-1)
1942                  else
1943                     z = z_v(varn,:)
1944                  end if
1945               else
1946                  if (isStrSubset( a(i),"zw_3d" ))then
1947                     z_v(varn,:) = z_w
1948                     if (log_z .EQ. 1) then
1949                        z = z_v(varn,1:dimz-1)
1950                     else
1951                        z = z_v(varn,:)
1952                     end if
1953                  end if                   
1954               end if
1955            end do           
1956         end if
1957
1958         if (isStrSubset(data@long_name," SR " ) ) then
1959           over = 0
1960         end if
1961         
1962         if (nof .EQ. 0) then
1963            z_(n,:)=z/norm_z
1964            z    = z_(n,:)
1965         else
1966            z=z/norm_z
1967         end if
1968
1969         delta_z = z(2) - z(1)
1970
1971         max_z_int=doubletoint(max_z/delta_z)
1972         min_z_int=doubletoint(min_z/delta_z)
1973
1974         if(max_z_int .eq. min_z_int)
1975             print(" ")
1976             print("Please increase 'max_z' or decrease 'min_z' so that "+\
1977                   "there are")
1978             print("at least two layers for the z-axis to plot")
1979             print(" ")
1980             exit
1981         end if
1982
1983         if(min_z_int .gt. max_z_int)
1984            print(" ")
1985            print("'min_z' is greater than 'max_z',")
1986            print("please change this")
1987            print(" ")
1988            exit
1989         end if
1990
1991         if (max_z_int .ge. dimz-1)
1992            max_z_int = dimz-1
1993            if (log_z .EQ. 1) then
1994             max_z_int = max_z_int - 1
1995            end if
1996         end if
1997
1998         if (min_z_int .lt. 0)
1999            min_z_int = 0
2000            if (log_z .EQ. 1) then
2001              min_z_int = min_z_int + 1
2002            end if
2003         end if         
2004
2005
2006         ;data can contain missing values in case of output of t=0h (inital profiles)
2007         ;where no output is possible
2008         n_not_ismissing = num(.not.ismissing(data(varn,:,:)))
2009
2010         if (n_not_ismissing .GT. 0 ) then   
2011 
2012           if (abs(min(data(varn,:,min_z_int:max_z_int))) .GT. 10)then
2013              min_value = abs(0.001*min(data(varn,:,min_z_int:max_z_int)))
2014              max_value = abs(0.001*max(data(varn,:,min_z_int:max_z_int)))
2015           else
2016              if (abs(min(data(varn,:,min_z_int:max_z_int))) .LT. 0.01 .AND. \
2017                  abs(max(data(varn,:,min_z_int:max_z_int))) .GT. 0.01)then
2018                 min_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
2019                 max_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
2020              else
2021                 if (abs(max(data(varn,:,min_z_int:max_z_int))) .LT. 0.01 .AND.\
2022                     abs(min(data(varn,:,min_z_int:max_z_int))) .GT. 0.01)then
2023                    min_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
2024                    max_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
2025                 else
2026                    min_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
2027                    max_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
2028                 end if
2029              end if
2030           end if
2031       
2032           if (min(data(varn,:,min_z_int:max_z_int)) .EQ. 0 .AND. \
2033               max(data(varn,:,min_z_int:max_z_int)) .EQ. 0)then
2034              min_value = 0.1
2035              max_value = 0.1
2036           end if
2037
2038         else     
2039           print(" ")
2040           print(vNam(varn) + " contains only missing values")
2041           print(" ")
2042         end if   
2043
2044         if (over .EQ. 0) then 
2045            res@gsnLeftString      = vNam(varn)
2046            res@tiXAxisString      = "("+unit(varn)+")"
2047            res@gsnRightString     = " "
2048            res@trYMinF            = min_z
2049            res@trYMaxF            = max_z
2050            if (xs .EQ. -1) then
2051               res@trXMinF            = min(data(varn,:,min_z_int:max_z_int))-\
2052                                                                      min_value
2053            else
2054               res@trXMinF            = xs     
2055            end if
2056            if (xe .EQ. -1) then
2057               res@trXMaxF            = max(data(varn,:,min_z_int:max_z_int))+\
2058                                                                      max_value
2059            else
2060               res@trXMaxF            = xe 
2061            end if         
2062            plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res) 
2063         end if
2064 
2065     
2066         if (vNam(varn) .EQ. "u" ) then
2067               miniu=min(data(varn,:,min_z_int:max_z_int))-min_value
2068               maxiu=max(data(varn,:,min_z_int:max_z_int))+max_value
2069               if (over .EQ. 1) then
2070                  res@xyDashPattern  = 0
2071                  plot_u = gsn_csm_xy(wks,data(varn,:,:),z,res)
2072               else
2073                  res@gsnLeftString      = vNam(varn)
2074                  res@tiXAxisString      = "("+unit(varn)+")"
2075                  res@gsnRightString     = " "                 
2076                  if (xs .EQ. -1) then 
2077                     res@trXMinF            = miniu
2078                  else
2079                     res@trXMinF            = xs
2080                  end if
2081                  if (xe .EQ. -1) then       
2082                     res@trXMaxF            = maxiu
2083                  else
2084                     res@trXMaxF            = xe 
2085                  end if               
2086                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res) 
2087               end if
2088            end if
2089            if (vNam(varn) .EQ. "v") then
2090               miniv=min(data(varn,:,min_z_int:max_z_int))-min_value
2091               maxiv=max(data(varn,:,min_z_int:max_z_int))+max_value
2092               if (over .EQ. 1) then
2093                  res@xyMonoDashPattern = True
2094                  res@xyDashPattern  = 1
2095                  plot_v = gsn_csm_xy(wks,data(varn,:,:),z,res)
2096               else
2097                  res@gsnLeftString      = vNam(varn)
2098                  res@tiXAxisString      = "("+unit(varn)+")"
2099                  res@gsnRightString     = " "                 
2100                  if (xs .EQ. -1) then
2101                     res@trXMinF            = miniv
2102                  else
2103                     res@trXMinF            = xs
2104                  end if
2105                  if (xe .EQ. -1) then
2106                     res@trXMaxF            = maxiv
2107                  else
2108                     res@trXMaxF            = xe 
2109                  end if               
2110                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2111               end if 
2112            end if
2113            if (vNam(varn) .EQ. "w") then
2114               miniw=min(data(varn,:,min_z_int:max_z_int))-min_value
2115               maxiw=max(data(varn,:,min_z_int:max_z_int))+max_value
2116               if (over .EQ. 1) then
2117                  res@xyDashPattern = 2
2118                  plot_w = gsn_csm_xy(wks,data(varn,:,:),z,res)
2119               else
2120                  res@gsnLeftString      = vNam(varn)
2121                  res@tiXAxisString      = "("+unit(varn)+")"
2122                  res@gsnRightString     = " "
2123                  if (xs .EQ. -1) then
2124                     res@trXMinF            = miniw
2125                  else
2126                     res@trXMinF            = xs
2127                  end if
2128                  if (xe .EQ. -1) then
2129                     res@trXMaxF            = maxiw
2130                  else
2131                     res@trXMaxF            = xe 
2132                  end if           
2133                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2134               end if   
2135            end if
2136
2137            if (vNam(varn) .EQ. "pt") then
2138               minipt=min(data(varn,:,min_z_int:max_z_int))-min_value
2139               maxipt=max(data(varn,:,min_z_int:max_z_int))+max_value
2140               if (over .EQ. 1) then
2141                  res@xyDashPattern  = 0
2142                  plot_pt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2143               else
2144                  res@gsnLeftString      = vNam(varn)
2145                  res@tiXAxisString      = "("+unit(varn)+")"
2146                  res@gsnRightString     = " "
2147                  if (xs .EQ. -1) then
2148                     res@trXMinF            = minipt
2149                  else
2150                     res@trXMinF            = xs
2151                  end if
2152                  if (xe .EQ. -1) then       
2153                     res@trXMaxF            = maxipt
2154                  else
2155                     res@trXMaxF            = xe 
2156                  end if
2157                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2158               end if 
2159            end if
2160            if (vNam(varn) .EQ. "vpt") then
2161               minivpt=min(data(varn,:,min_z_int:max_z_int))-min_value
2162               maxivpt=max(data(varn,:,min_z_int:max_z_int))+max_value
2163               if (over .EQ. 1) then
2164                  res@xyDashPattern  = 1
2165                  plot_vpt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2166               else
2167                  res@gsnLeftString      = vNam(varn)
2168                  res@tiXAxisString      = "("+unit(varn)+")"
2169                  res@gsnRightString     = " "
2170                  if (xs .EQ. -1) then
2171                     res@trXMinF            = minivpt
2172                  else
2173                     res@trXMinF            = xs
2174                  end if
2175                  if (xe .EQ. -1) then       
2176                     res@trXMaxF            = maxivpt
2177                  else
2178                     res@trXMaxF            = xe 
2179                  end if
2180                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2181               end if
2182            end if
2183            if (vNam(varn) .EQ. "lpt") then
2184               minilpt=min(data(varn,:,min_z_int:max_z_int))-min_value
2185               maxilpt=max(data(varn,:,min_z_int:max_z_int))+max_value
2186               if (over .EQ. 1) then
2187                  res@xyDashPattern  = 2
2188                  plot_lpt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2189               else
2190                  res@gsnLeftString      = vNam(varn)
2191                  res@tiXAxisString      = "("+unit(varn)+")"
2192                  res@gsnRightString     = " "
2193                  if (xs .EQ. -1) then
2194                     res@trXMinF            = minilpt
2195                  else
2196                     res@trXMinF            = xs
2197                  end if
2198                  if (xe .EQ. -1) then       
2199                     res@trXMaxF            = maxilpt
2200                  else
2201                     res@trXMaxF            = xe 
2202                  end if
2203                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2204               end if
2205            end if
2206
2207            if (vNam(varn) .EQ. "q") then
2208               miniq=min(data(varn,:,min_z_int:max_z_int))-min_value
2209               maxiq=max(data(varn,:,min_z_int:max_z_int))+max_value
2210               if (over .EQ. 1) then
2211                  res@xyDashPattern  = 0
2212                  plot_q = gsn_csm_xy(wks,data(varn,:,:),z,res)
2213               else
2214                  res@gsnLeftString      = vNam(varn)
2215                  res@tiXAxisString      = "("+unit(varn)+")"
2216                  res@gsnRightString     = " "
2217                  if (xs .EQ. -1) then
2218                     res@trXMinF            = miniq
2219                  else
2220                     res@trXMinF            = xs
2221                  end if
2222                  if (xe .EQ. -1) then       
2223                     res@trXMaxF            = maxiq
2224                  else
2225                     res@trXMaxF            = xe 
2226                  end if
2227                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2228               end if
2229            end if
2230            if (vNam(varn) .EQ. "qv") then
2231               miniqv=min(data(varn,:,min_z_int:max_z_int))-min_value
2232               maxiqv=max(data(varn,:,min_z_int:max_z_int))+max_value
2233               if (over .EQ. 1) then
2234                  res@xyDashPattern  = 1
2235                  plot_qv = gsn_csm_xy(wks,data(varn,:,:),z,res)
2236               else
2237                  res@gsnLeftString      = vNam(varn)
2238                  res@tiXAxisString      = "("+unit(varn)+")"
2239                  res@gsnRightString     = " "
2240                  if (xs .EQ. -1) then
2241                     res@trXMinF            = miniqv
2242                  else
2243                     res@trXMinF            = xs
2244                  end if
2245                  if (xe .EQ. -1) then         
2246                     res@trXMaxF            = maxiqv
2247                  else
2248                     res@trXMaxF            = xe 
2249                  end if
2250                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2251               end if
2252            end if
2253            if (vNam(varn) .EQ. "ql") then
2254               miniql=min(data(varn,:,min_z_int:max_z_int))-min_value
2255               maxiql=max(data(varn,:,min_z_int:max_z_int))+max_value
2256               if (over .EQ. 1) then
2257                  res@xyDashPattern  = 2
2258                  plot_ql = gsn_csm_xy(wks,data(varn,:,:),z,res)
2259               else
2260                  res@gsnLeftString      = vNam(varn)
2261                  res@tiXAxisString      = "("+unit(varn)+")"
2262                  res@gsnRightString     = " "
2263                  if (xs .EQ. -1) then
2264                     res@trXMinF            = miniql
2265                  else
2266                     res@trXMinF            = xs
2267                  end if
2268                  if (xe .EQ. -1) then       
2269                     res@trXMaxF            = maxiql
2270                  else
2271                     res@trXMaxF            = xe 
2272                  end if
2273                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2274               end if
2275            end if
2276
2277            if (vNam(varn) .EQ. "e") then
2278               minie=min(data(varn,:,min_z_int:max_z_int))-min_value
2279               maxie=max(data(varn,:,min_z_int:max_z_int))+max_value
2280               if (over .EQ. 1) then
2281                  res@xyDashPattern  = 0
2282                  plot_e = gsn_csm_xy(wks,data(varn,:,:),z,res)
2283               else
2284                  res@gsnLeftString      = vNam(varn)
2285                  res@tiXAxisString      = "("+unit(varn)+")"
2286                  res@gsnRightString     = " "
2287                  if (xs .EQ. -1) then
2288                     res@trXMinF            = minie
2289                  else
2290                     res@trXMinF            = xs
2291                  end if
2292                  if (xe .EQ. -1) then       
2293                     res@trXMaxF            = maxie
2294                  else
2295                     res@trXMaxF            = xe 
2296                  end if
2297                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2298               end if
2299            end if
2300            if (vNam(varn) .EQ. "es" .OR. vNam(varn) .EQ. "e*") then
2301               minies=min(data(varn,:,min_z_int:max_z_int))-min_value
2302               maxies=max(data(varn,:,min_z_int:max_z_int))+max_value
2303               if (over .EQ. 1) then
2304                  res@xyDashPattern  = 1
2305                  plot_es = gsn_csm_xy(wks,data(varn,:,:),z,res)
2306               else
2307                  res@gsnLeftString      = vNam(varn)
2308                  res@tiXAxisString      = "("+unit(varn)+")"
2309                  res@gsnRightString     = " "
2310                  if (xs .EQ. -1) then
2311                     res@trXMinF            = minies
2312                  else
2313                     res@trXMinF            = xs
2314                  end if
2315                  if (xe .EQ. -1) then       
2316                     res@trXMaxF            = maxies
2317                  else
2318                     res@trXMaxF            = xe 
2319                  end if
2320                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2321               end if
2322            end if
2323
2324            if (vNam(varn) .EQ. "km") then
2325               minikm=min(data(varn,:,min_z_int:max_z_int))-min_value
2326               maxikm=max(data(varn,:,min_z_int:max_z_int))+max_value
2327               if (over .EQ. 1) then
2328                  res@xyDashPattern  = 0
2329                  plot_km = gsn_csm_xy(wks,data(varn,:,:),z,res)
2330               else
2331                  res@gsnLeftString      = vNam(varn)
2332                  res@tiXAxisString      = "("+unit(varn)+")"
2333                  res@gsnRightString     = " "
2334                  if (xs .EQ. -1) then
2335                     res@trXMinF            = minikm
2336                  else
2337                     res@trXMinF            = xs
2338                  end if
2339                  if (xe .EQ. -1) then     
2340                     res@trXMaxF            = maxikm
2341                  else
2342                     res@trXMaxF            = xe 
2343                  end if
2344                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2345               end if
2346            end if
2347            if (vNam(varn) .EQ. "kh") then
2348               minikh=min(data(varn,:,min_z_int:max_z_int))-min_value
2349               maxikh=max(data(varn,:,min_z_int:max_z_int))+max_value
2350               if (over .EQ. 1) then
2351                  res@xyDashPattern  = 1
2352                  plot_kh = gsn_csm_xy(wks,data(varn,:,:),z,res)
2353               else
2354                  res@gsnLeftString      = vNam(varn)
2355                  res@tiXAxisString      = "("+unit(varn)+")"
2356                  res@gsnRightString     = " "
2357                  if (xs .EQ. -1) then
2358                     res@trXMinF            = minikh
2359                  else
2360                     res@trXMinF            = xs
2361                  end if
2362                  if (xe .EQ. -1) then     
2363                     res@trXMaxF            = maxikh
2364                  else
2365                     res@trXMaxF            = xe 
2366                  end if
2367                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2368               end if
2369            end if
2370
2371            if (vNam(varn) .EQ. "wpup" .OR. vNam(varn) .EQ. "w"+dq+"u"+dq) then
2372               miniwpup=min(data(varn,:,min_z_int:max_z_int))-min_value
2373               maxiwpup=max(data(varn,:,min_z_int:max_z_int))+max_value
2374               if (over .EQ. 1) then
2375                  res@xyDashPattern  = 0
2376                  plot_wpup = gsn_csm_xy(wks,data(varn,:,:),z,res)
2377               else
2378                  res@gsnLeftString      = vNam(varn)
2379                  res@tiXAxisString      = "("+unit(varn)+")"
2380                  res@gsnRightString     = " "
2381                  if (xs .EQ. -1) then
2382                     res@trXMinF            = miniwpup
2383                  else
2384                     res@trXMinF            = xs
2385                  end if
2386                  if (xe .EQ. -1) then
2387                     res@trXMaxF            = maxiwpup
2388                  else
2389                     res@trXMaxF            = xe 
2390                  end if
2391                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2392               end if
2393            end if
2394            if (vNam(varn) .EQ. "wsus" .OR. vNam(varn) .EQ. "w*u*") then
2395               miniwsus=min(data(varn,:,min_z_int:max_z_int))-min_value
2396               maxiwsus=max(data(varn,:,min_z_int:max_z_int))+max_value
2397               if (over .EQ. 1) then
2398                  res@xyDashPattern  = 1
2399                  plot_wsus = gsn_csm_xy(wks,data(varn,:,:),z,res)
2400               else
2401                  res@gsnLeftString      = vNam(varn)
2402                  res@tiXAxisString      = "("+unit(varn)+")"
2403                  res@gsnRightString     = " "
2404                  if (xs .EQ. -1) then
2405                     res@trXMinF            = miniwsus
2406                  else
2407                     res@trXMinF            = xs
2408                  end if
2409                  if (xe .EQ. -1) then     
2410                     res@trXMaxF            = maxiwsus
2411                  else
2412                     res@trXMaxF            = xe 
2413                  end if
2414                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2415               end if
2416            end if
2417            if (vNam(varn) .EQ. "wu") then
2418               miniwu=min(data(varn,:,min_z_int:max_z_int))-min_value
2419               maxiwu=max(data(varn,:,min_z_int:max_z_int))+max_value
2420               if (over .EQ. 1) then
2421                  res@xyDashPattern  = 2
2422                  plot_wu = gsn_csm_xy(wks,data(varn,:,:),z,res)
2423               else
2424                  res@gsnLeftString      = vNam(varn)
2425                  res@tiXAxisString      = "("+unit(varn)+")"
2426                  res@gsnRightString     = " "
2427                  if (xs .EQ. -1) then
2428                     res@trXMinF            = miniwu
2429                  else
2430                     res@trXMinF            = xs
2431                  end if
2432                  if (xe .EQ. -1) then
2433                     res@trXMaxF            = maxiwu
2434                  else
2435                     res@trXMaxF            = xe 
2436                  end if
2437                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2438               end if
2439            end if
2440
2441            if (vNam(varn) .EQ. "wpvp" .OR. vNam(varn) .EQ. "w"+dq+"v"+dq) then
2442               miniwpvp=min(data(varn,:,min_z_int:max_z_int))-min_value
2443               maxiwpvp=max(data(varn,:,min_z_int:max_z_int))+max_value
2444               if (over .EQ. 1) then
2445                  res@xyDashPattern  = 0
2446                  plot_wpvp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2447               else
2448                  res@gsnLeftString      = vNam(varn)
2449                  res@tiXAxisString      = "("+unit(varn)+")"
2450                  res@gsnRightString     = " "
2451                  if (xs .EQ. -1) then
2452                     res@trXMinF            = miniwpvp
2453                  else
2454                     res@trXMinF            = xs
2455                  end if
2456                  if (xe .EQ. -1) then     
2457                     res@trXMaxF            = maxiwpvp
2458                  else
2459                     res@trXMaxF            = xe 
2460                  end if
2461                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2462               end if
2463            end if
2464            if (vNam(varn) .EQ. "wsvs" .OR. vNam(varn) .EQ. "w*v*") then
2465               miniwsvs=min(data(varn,:,min_z_int:max_z_int))-min_value
2466               maxiwsvs=max(data(varn,:,min_z_int:max_z_int))+max_value
2467               if (over .EQ. 1) then
2468                  res@xyDashPattern  = 1
2469                  plot_wsvs = gsn_csm_xy(wks,data(varn,:,:),z,res)
2470               else
2471                  res@gsnLeftString      = vNam(varn)
2472                  res@tiXAxisString      = "("+unit(varn)+")"
2473                  res@gsnRightString     = " "
2474                  if (xs .EQ. -1) then
2475                     res@trXMinF            = miniwsvs
2476                  else
2477                     res@trXMinF            = xs
2478                  end if
2479                  if (xe .EQ. -1) then     
2480                     res@trXMaxF            = maxiwsvs
2481                  else
2482                     res@trXMaxF            = xe 
2483                  end if
2484                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2485               end if
2486            end if
2487            if (vNam(varn) .EQ. "wv") then
2488               miniwv=min(data(varn,:,min_z_int:max_z_int))-min_value
2489               maxiwv=max(data(varn,:,min_z_int:max_z_int))+max_value
2490               if (over .EQ. 1) then
2491                  res@xyDashPattern  = 2
2492                  plot_wv = gsn_csm_xy(wks,data(varn,:,:),z,res)
2493               else
2494                  res@gsnLeftString      = vNam(varn)
2495                  res@tiXAxisString      = "("+unit(varn)+")"
2496                  res@gsnRightString     = " "
2497                  if (xs .EQ. -1) then
2498                     res@trXMinF            = miniwv
2499                  else
2500                     res@trXMinF            = xs
2501                  end if
2502                  if (xe .EQ. -1) then       
2503                     res@trXMaxF            = maxiwv
2504                  else
2505                     res@trXMaxF            = xe 
2506                  end if
2507                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2508               end if
2509            end if
2510
2511            if (vNam(varn) .EQ. "wpptp" .OR. vNam(varn) \
2512               .EQ. "w"+dq+"pt"+dq) then
2513               miniwpptp=min(data(varn,:,min_z_int:max_z_int))-min_value
2514               maxiwpptp=max(data(varn,:,min_z_int:max_z_int))+max_value
2515               if (over .EQ. 1) then
2516                  res@xyDashPattern  = 0
2517                  plot_wpptp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2518               else
2519                  res@gsnLeftString      = vNam(varn)
2520                  res@tiXAxisString      = "("+unit(varn)+")"
2521                  res@gsnRightString     = " "
2522                  if (xs .EQ. -1) then
2523                     res@trXMinF            = miniwpptp
2524                  else
2525                     res@trXMinF            = xs
2526                  end if
2527                  if (xe .EQ. -1) then       
2528                     res@trXMaxF            = maxiwpptp
2529                  else
2530                     res@trXMaxF            = xe 
2531                  end if
2532                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2533               end if
2534            end if
2535            if (vNam(varn) .EQ. "wspts" .OR. vNam(varn) .EQ. "w*pt*") then
2536               miniwspts=min(data(varn,:,min_z_int:max_z_int))-min_value
2537               maxiwspts=max(data(varn,:,min_z_int:max_z_int))+max_value
2538               if (over .EQ. 1) then
2539                  res@xyDashPattern  = 1
2540                  plot_wspts = gsn_csm_xy(wks,data(varn,:,:),z,res)
2541               else
2542                  res@gsnLeftString      = vNam(varn)
2543                  res@tiXAxisString      = "("+unit(varn)+")"
2544                  res@gsnRightString     = " "
2545                  if (xs .EQ. -1) then
2546                     res@trXMinF            = miniwspts
2547                  else
2548                     res@trXMinF            = xs
2549                  end if
2550                  if (xe .EQ. -1) then       
2551                     res@trXMaxF            = maxiwspts
2552                  else
2553                     res@trXMaxF            = xe 
2554                  end if
2555                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2556               end if
2557            end if
2558            if (vNam(varn) .EQ. "wpt") then       
2559               miniwpt=min(data(varn,:,min_z_int:max_z_int))-min_value
2560               maxiwpt=max(data(varn,:,min_z_int:max_z_int))+max_value
2561               if (over .EQ. 1) then
2562                  res@xyDashPattern  = 2
2563                  plot_wpt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2564               else
2565                  res@gsnLeftString      = vNam(varn)
2566                  res@tiXAxisString      = "("+unit(varn)+")"
2567                  res@gsnRightString     = " "
2568                  if (xs .EQ. -1) then
2569                     res@trXMinF            = miniwpt
2570                  else
2571                     res@trXMinF            = xs
2572                  end if
2573                  if (xe .EQ. -1) then       
2574                     res@trXMaxF            = maxiwpt
2575                  else
2576                     res@trXMaxF            = xe 
2577                  end if
2578                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2579               end if
2580            end if
2581
2582            if (vNam(varn) .EQ. "wsptsBC".OR. vNam(varn) .EQ. "w*pt*BC" ) then
2583               miniwsptsBC=min(data(varn,:,min_z_int:max_z_int))-min_value
2584               maxiwsptsBC=max(data(varn,:,min_z_int:max_z_int))+max_value
2585               if (over .EQ. 1) then
2586                  res@xyDashPattern  = 0
2587                  plot_wsptsBC = gsn_csm_xy(wks,data(varn,:,:),z,res)
2588               else
2589                  res@gsnLeftString      = vNam(varn)
2590                  res@tiXAxisString      = "("+unit(varn)+")"
2591                  res@gsnRightString     = " "
2592                  if (xs .EQ. -1) then
2593                     res@trXMinF            = miniwsptsBC
2594                  else
2595                     res@trXMinF            = xs
2596                  end if
2597                  if (xe .EQ. -1) then       
2598                     res@trXMaxF            = maxiwsptsBC
2599                  else
2600                     res@trXMaxF            = xe 
2601                  end if
2602                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2603               end if
2604            end if             
2605            if (vNam(varn) .EQ. "wptBC") then
2606               miniwptBC=min(data(varn,:,min_z_int:max_z_int))-min_value
2607               maxiwptBC=max(data(varn,:,min_z_int:max_z_int))+max_value
2608               if (over .EQ. 1) then
2609                  res@xyDashPattern  = 1
2610                  plot_wptBC = gsn_csm_xy(wks,data(varn,:,:),z,res)
2611               else
2612                  res@gsnLeftString      = vNam(varn)
2613                  res@tiXAxisString      = "("+unit(varn)+")"
2614                  res@gsnRightString     = " "
2615                  if (xs .EQ. -1) then
2616                     res@trXMinF            = miniwptBC
2617                  else
2618                     res@trXMinF            = xs
2619                  end if
2620                  if (xe .EQ. -1) then       
2621                     res@trXMaxF            = maxiwptBC
2622                  else
2623                     res@trXMaxF            = xe 
2624                  end if
2625                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2626               end if
2627            end if
2628
2629            if (vNam(varn) .EQ. "wpvptp" .OR. vNam(varn) \
2630                .EQ. "w"+dq+"vpt"+dq) then
2631               miniwpvptp=min(data(varn,:,min_z_int:max_z_int))-min_value
2632               maxiwpvptp=max(data(varn,:,min_z_int:max_z_int))+max_value
2633               if (over .EQ. 1) then
2634                  res@xyDashPattern  = 0
2635                  plot_wpvptp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2636               else
2637                  res@gsnLeftString      = vNam(varn)
2638                  res@tiXAxisString      = "("+unit(varn)+")"
2639                  res@gsnRightString     = " "
2640                  if (xs .EQ. -1) then
2641                     res@trXMinF            = miniwpvptp
2642                  else
2643                     res@trXMinF            = xs
2644                  end if
2645                  if (xe .EQ. -1) then       
2646                     res@trXMaxF            = maxiwpvptp
2647                  else
2648                     res@trXMaxF            = xe 
2649                  end if
2650                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2651               end if
2652            end if
2653            if (vNam(varn) .EQ. "wsvpts" .OR. vNam(varn) .EQ. "w*vpt*") then
2654               miniwsvpts=min(data(varn,:,min_z_int:max_z_int))-min_value
2655               maxiwsvpts=max(data(varn,:,min_z_int:max_z_int))+max_value
2656               if (over .EQ. 1) then
2657                  res@xyDashPattern  = 1
2658                  plot_wsvpts = gsn_csm_xy(wks,data(varn,:,:),z,res)
2659               else
2660                  res@gsnLeftString      = vNam(varn)
2661                  res@tiXAxisString      = "("+unit(varn)+")"
2662                  res@gsnRightString     = " "
2663                  if (xs .EQ. -1) then
2664                     res@trXMinF            = miniwsvpts
2665                  else
2666                     res@trXMinF            = xs
2667                  end if
2668                  if (xe .EQ. -1) then       
2669                     res@trXMaxF            = maxiwsvpts
2670                  else
2671                     res@trXMaxF            = xe 
2672                  end if
2673                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2674               end if
2675            end if
2676            if (vNam(varn) .EQ. "wvpt") then
2677               miniwvpt=min(data(varn,:,min_z_int:max_z_int))-min_value
2678               maxiwvpt=max(data(varn,:,min_z_int:max_z_int))+max_value
2679               if (over .EQ. 1) then
2680                  res@xyDashPattern  = 2
2681                  plot_wvpt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2682               else
2683                  res@gsnLeftString      = vNam(varn)
2684                  res@tiXAxisString      = "("+unit(varn)+")"
2685                  res@gsnRightString     = " "
2686                  if (xs .EQ. -1) then
2687                     res@trXMinF            = miniwvpt
2688                  else
2689                     res@trXMinF            = xs
2690                  end if
2691                  if (xe .EQ. -1) then       
2692                     res@trXMaxF            = maxiwvpt
2693                  else
2694                     res@trXMaxF            = xe 
2695                  end if
2696                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2697               end if
2698            end if
2699
2700            if (vNam(varn) .EQ. "wpqp" .OR. vNam(varn) .EQ. "w"+dq+"q"+dq) then
2701               miniwpqp=min(data(varn,:,min_z_int:max_z_int))-min_value
2702               maxiwpqp=max(data(varn,:,min_z_int:max_z_int))+max_value
2703               if (over .EQ. 1) then
2704                  res@xyDashPattern  = 0
2705                  plot_wpqp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2706               else
2707                  res@gsnLeftString      = vNam(varn)
2708                  res@tiXAxisString      = "("+unit(varn)+")"
2709                  res@gsnRightString     = " "
2710                  if (xs .EQ. -1) then
2711                     res@trXMinF            = miniwpqp
2712                  else
2713                     res@trXMinF            = xs
2714                  end if
2715                  if (xe .EQ. -1) then       
2716                     res@trXMaxF            = maxiwpqp
2717                  else
2718                     res@trXMaxF            = xe 
2719                  end if
2720                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2721               end if
2722            end if
2723            if (vNam(varn) .EQ. "wsqs".OR. vNam(varn) .EQ. "w*s*" ) then
2724               miniwsqs=min(data(varn,:,min_z_int:max_z_int))-min_value
2725               maxiwsqs=max(data(varn,:,min_z_int:max_z_int))+max_value
2726               if (over .EQ. 1) then
2727                  res@xyDashPattern  = 1
2728                  plot_wsqs = gsn_csm_xy(wks,data(varn,:,:),z,res)
2729               else
2730                  res@gsnLeftString      = vNam(varn)
2731                  res@tiXAxisString      = "("+unit(varn)+")"
2732                  res@gsnRightString     = " "
2733                  if (xs .EQ. -1) then
2734                     res@trXMinF            = miniwsqs
2735                  else
2736                     res@trXMinF            = xs
2737                  end if
2738                  if (xe .EQ. -1) then       
2739                     res@trXMaxF            = maxiwsqs
2740                  else
2741                     res@trXMaxF            = xe 
2742                  end if
2743                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2744               end if
2745            end if
2746            if (vNam(varn) .EQ. "wq") then
2747               miniwq=min(data(varn,:,min_z_int:max_z_int))-min_value
2748               maxiwq=max(data(varn,:,min_z_int:max_z_int))+max_value
2749               if (over .EQ. 1) then
2750                  res@xyDashPattern  = 2
2751                  plot_wq = gsn_csm_xy(wks,data(varn,:,:),z,res)
2752               else
2753                  res@gsnLeftString      = vNam(varn)
2754                  res@tiXAxisString      = "("+unit(varn)+")"
2755                  res@gsnRightString     = " "
2756                  if (xs .EQ. -1) then
2757                     res@trXMinF            = miniwq
2758                  else
2759                     res@trXMinF            = xs
2760                  end if
2761                  if (xe .EQ. -1) then       
2762                     res@trXMaxF            = maxiwq
2763                  else
2764                     res@trXMaxF            = xe 
2765                  end if
2766                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2767               end if
2768            end if
2769
2770            if (vNam(varn) .EQ. "wpqvp" .OR. \
2771               vNam(varn) .EQ. "w"+dq+"qv"+dq) then
2772               miniwpqvp=min(data(varn,:,min_z_int:max_z_int))-min_value
2773               maxiwpqvp=max(data(varn,:,min_z_int:max_z_int))+max_value
2774               if (over .EQ. 1) then
2775                  res@xyDashPattern  = 0
2776                  plot_wpqvp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2777               else
2778                  res@gsnLeftString      = vNam(varn)
2779                  res@tiXAxisString      = "("+unit(varn)+")"
2780                  res@gsnRightString     = " "
2781                  if (xs .EQ. -1) then
2782                     res@trXMinF            = miniwpqvp
2783                  else
2784                     res@trXMinF            = xs
2785                  end if
2786                  if (xe .EQ. -1) then       
2787                     res@trXMaxF            = maxiwpqvp
2788                  else
2789                     res@trXMaxF            = xe 
2790                  end if
2791                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2792               end if
2793            end if
2794            if (vNam(varn) .EQ. "wsqvs" .OR. vNam(varn) .EQ. "w*qv*") then
2795               miniwsqvs=min(data(varn,:,min_z_int:max_z_int))-min_value
2796               maxiwsqvs=max(data(varn,:,min_z_int:max_z_int))+max_value
2797               if (over .EQ. 1) then
2798                  res@xyDashPattern  = 1
2799                  plot_wsqvs = gsn_csm_xy(wks,data(varn,:,:),z,res)
2800               else
2801                  res@gsnLeftString      = vNam(varn)
2802                  res@tiXAxisString      = "("+unit(varn)+")"
2803                  res@gsnRightString     = " "
2804                  if (xs .EQ. -1) then
2805                     res@trXMinF            = miniwsqvs
2806                  else
2807                     res@trXMinF            = xs
2808                  end if
2809                  if (xe .EQ. -1) then       
2810                     res@trXMaxF            = maxiwsqvs
2811                  else
2812                     res@trXMaxF            = xe 
2813                  end if
2814                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2815               end if
2816            end if
2817            if (vNam(varn) .EQ. "wqv") then
2818               miniwqv=min(data(varn,:,min_z_int:max_z_int))-min_value
2819               maxiwqv=max(data(varn,:,min_z_int:max_z_int))+max_value
2820               if (over .EQ. 1) then
2821                  res@xyDashPattern  = 2
2822                  plot_wqv = gsn_csm_xy(wks,data(varn,:,:),z,res)
2823               else
2824                  res@gsnLeftString      = vNam(varn)
2825                  res@tiXAxisString      = "("+unit(varn)+")"
2826                  res@gsnRightString     = " "
2827                  if (xs .EQ. -1) then
2828                     res@trXMinF            = miniwqv
2829                  else
2830                     res@trXMinF            = xs
2831                  end if
2832                  if (xe .EQ. -1) then       
2833                     res@trXMaxF            = maxiwqv
2834                  else
2835                     res@trXMaxF            = xe 
2836                  end if
2837                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2838               end if
2839            end if
2840
2841            if (vNam(varn) .EQ. "wpsp" .OR. vNam(varn) .EQ. "w"+dq+"s"+dq) then
2842               miniwpsp=min(data(varn,:,min_z_int:max_z_int))-min_value
2843               maxiwpsp=max(data(varn,:,min_z_int:max_z_int))+max_value
2844               if (over .EQ. 1) then
2845                  res@xyDashPattern  = 0
2846                  plot_wpsp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2847               else
2848                  res@gsnLeftString      = vNam(varn)
2849                  res@tiXAxisString      = "("+unit(varn)+")"
2850                  res@gsnRightString     = " "
2851                  if (xs .EQ. -1) then
2852                     res@trXMinF            = miniwpsp
2853                  else
2854                     res@trXMinF            = xs
2855                  end if
2856                  if (xe .EQ. -1) then       
2857                     res@trXMaxF            = maxiwpsp
2858                  else
2859                     res@trXMaxF            = xe 
2860                  end if
2861                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2862               end if
2863            end if
2864            if (vNam(varn) .EQ. "wsss" .OR. vNam(varn) .EQ. "w*s*" ) then
2865               miniwsss=min(data(varn,:,min_z_int:max_z_int))-min_value
2866               maxiwsss=max(data(varn,:,min_z_int:max_z_int))+max_value
2867               if (over .EQ. 1) then
2868                  res@xyDashPattern  = 1
2869                  plot_wsss = gsn_csm_xy(wks,data(varn,:,:),z,res)
2870               else
2871                  res@gsnLeftString      = vNam(varn)
2872                  res@tiXAxisString      = "("+unit(varn)+")"
2873                  res@gsnRightString     = " "
2874                  if (xs .EQ. -1) then
2875                     res@trXMinF            = miniwsss
2876                  else
2877                     res@trXMinF            = xs
2878                  end if
2879                  if (xe .EQ. -1) then       
2880                     res@trXMaxF            = maxiwsss
2881                  else
2882                     res@trXMaxF            = xe 
2883                  end if
2884                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2885               end if
2886            end if
2887            if (vNam(varn) .EQ. "ws") then
2888               miniws=min(data(varn,:,min_z_int:max_z_int))-min_value
2889               maxiws=max(data(varn,:,min_z_int:max_z_int))+max_value
2890               if (over .EQ. 1) then
2891                  res@xyDashPattern  = 2
2892                  plot_ws = gsn_csm_xy(wks,data(varn,:,:),z,res)
2893               else
2894                  res@gsnLeftString      = vNam(varn)
2895                  res@tiXAxisString      = "("+unit(varn)+")"
2896                  res@gsnRightString     = " "
2897                  if (xs .EQ. -1) then
2898                     res@trXMinF            = miniws
2899                  else
2900                     res@trXMinF            = xs
2901                  end if
2902                  if (xe .EQ. -1) then       
2903                     res@trXMaxF            = maxiws
2904                  else
2905                     res@trXMaxF            = xe 
2906                  end if
2907                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2908               end if
2909            end if
2910
2911            if (vNam(varn) .EQ. "wpsap" .OR. \
2912                vNam(varn) .EQ. "w"+dq+"sa"+dq) then
2913               miniwpsap=min(data(varn,:,min_z_int:max_z_int))-min_value
2914               maxiwpsap=max(data(varn,:,min_z_int:max_z_int))+max_value
2915               if (over .EQ. 1) then
2916                  res@xyDashPattern  = 0
2917                  plot_wpsap = gsn_csm_xy(wks,data(varn,:,:),z,res)
2918               else
2919                  res@gsnLeftString      = vNam(varn)
2920                  res@tiXAxisString      = "("+unit(varn)+")"
2921                  res@gsnRightString     = " "
2922                  if (xs .EQ. -1) then
2923                     res@trXMinF            = miniwpsap
2924                  else
2925                     res@trXMinF            = xs
2926                  end if
2927                  if (xe .EQ. -1) then       
2928                     res@trXMaxF            = maxiwpsap
2929                  else
2930                     res@trXMaxF            = xe 
2931                  end if
2932                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2933               end if
2934            end if
2935            if (vNam(varn) .EQ. "wssas" .OR. vNam(varn) .EQ. "w*sa*") then
2936               miniwssas=min(data(varn,:,min_z_int:max_z_int))-min_value
2937               maxiwssas=max(data(varn,:,min_z_int:max_z_int))+max_value
2938               if (over .EQ. 1) then
2939                  res@xyDashPattern  = 1
2940                  plot_wssas = gsn_csm_xy(wks,data(varn,:,:),z,res)
2941               else
2942                  res@gsnLeftString      = vNam(varn)
2943                  res@tiXAxisString      = "("+unit(varn)+")"
2944                  res@gsnRightString     = " "
2945                  if (xs .EQ. -1) then
2946                     res@trXMinF            = miniwssas
2947                  else
2948                     res@trXMinF            = xs
2949                  end if
2950                  if (xe .EQ. -1) then       
2951                     res@trXMaxF            = maxiwssas
2952                  else
2953                     res@trXMaxF            = xe 
2954                  end if
2955                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2956               end if
2957            end if
2958            if (vNam(varn) .EQ. "wsa") then
2959               miniwsa=min(data(varn,:,min_z_int:max_z_int))-min_value
2960               maxiwsa=max(data(varn,:,min_z_int:max_z_int))+max_value
2961               if (over .EQ. 1) then
2962                  res@xyDashPattern  = 2
2963                  plot_wsa = gsn_csm_xy(wks,data(varn,:,:),z,res)
2964               else
2965                  res@gsnLeftString      = vNam(varn)
2966                  res@tiXAxisString      = "("+unit(varn)+")"
2967                  res@gsnRightString     = " "
2968                  if (xs .EQ. -1) then
2969                     res@trXMinF            = miniwsa
2970                  else
2971                     res@trXMinF            = xs
2972                  end if
2973                  if (xe .EQ. -1) then       
2974                     res@trXMaxF            = maxiwsa
2975                  else
2976                     res@trXMaxF            = xe 
2977                  end if
2978                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2979               end if
2980            end if
2981
2982            if (vNam(varn) .EQ. "us2" .OR. vNam(varn) .EQ. "u*2") then
2983               minius2=min(data(varn,:,min_z_int:max_z_int))-min_value
2984               maxius2=max(data(varn,:,min_z_int:max_z_int))+max_value
2985               if (over .EQ. 1) then
2986                  res@xyDashPattern  = 0
2987                  plot_us2 = gsn_csm_xy(wks,data(varn,:,:),z,res)
2988               else
2989                  res@gsnLeftString      = vNam(varn)
2990                  res@tiXAxisString      = "("+unit(varn)+")"
2991                  res@gsnRightString     = " "
2992                  if (xs .EQ. -1) then
2993                     res@trXMinF            = minius2
2994                  else
2995                     res@trXMinF            = xs
2996                  end if
2997                  if (xe .EQ. -1) then       
2998                     res@trXMaxF            = maxius2
2999                  else
3000                     res@trXMaxF            = xe 
3001                  end if
3002                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3003               end if
3004            end if
3005            if (vNam(varn) .EQ. "vs2" .OR. vNam(varn) .EQ. "v*2") then
3006               minivs2=min(data(varn,:,min_z_int:max_z_int))-min_value
3007               maxivs2=max(data(varn,:,min_z_int:max_z_int))+max_value
3008               if (over .EQ. 1) then
3009                  res@xyDashPattern  = 1
3010                  plot_vs2 = gsn_csm_xy(wks,data(varn,:,:),z,res)
3011               else
3012                  res@gsnLeftString      = vNam(varn)
3013                  res@tiXAxisString      = "("+unit(varn)+")"
3014                  res@gsnRightString     = " "
3015                  if (xs .EQ. -1) then
3016                     res@trXMinF            = minivs2
3017                  else
3018                     res@trXMinF            = xs
3019                  end if
3020                  if (xe .EQ. -1) then       
3021                     res@trXMaxF            = maxivs2
3022                  else
3023                     res@trXMaxF            = xe 
3024                  end if
3025                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3026               end if
3027            end if
3028            if (vNam(varn) .EQ. "ws2" .OR. vNam(varn) .EQ. "w*2") then
3029               miniws2=min(data(varn,:,min_z_int:max_z_int))-min_value
3030               maxiws2=max(data(varn,:,min_z_int:max_z_int))+max_value
3031               if (over .EQ. 1) then
3032                  res@xyDashPattern  = 2
3033                  plot_ws2 = gsn_csm_xy(wks,data(varn,:,:),z,res)
3034               else
3035                  res@gsnLeftString      = vNam(varn)
3036                  res@tiXAxisString      = "("+unit(varn)+")"
3037                  res@gsnRightString     = " "
3038                  if (xs .EQ. -1) then
3039                     res@trXMinF            = miniws2
3040                  else
3041                     res@trXMinF            = xs
3042                  end if
3043                  if (xe .EQ. -1) then       
3044                     res@trXMaxF            = maxiws2
3045                  else
3046                     res@trXMaxF            = xe 
3047                  end if
3048                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3049               end if
3050            end if
3051
3052            if (vNam(varn) .EQ. "wsususodz" .OR. \
3053                vNam(varn) .EQ. "w*u*u*:dz") then
3054               miniwsususodz=min(data(varn,:,min_z_int:max_z_int))-min_value
3055               maxiwsususodz=max(data(varn,:,min_z_int:max_z_int))+max_value
3056               if (over .EQ. 1) then
3057                  res@xyDashPattern  = 0
3058                  plot_wsususodz = gsn_csm_xy(wks,data(varn,:,:),z,res)
3059               else
3060                  res@gsnLeftString      = vNam(varn)
3061                  res@tiXAxisString      = "("+unit(varn)+")"
3062                  res@gsnRightString     = " "
3063                  if (xs .EQ. -1) then
3064                     res@trXMinF            = miniwsususodz
3065                  else
3066                     res@trXMinF            = xs
3067                  end if
3068                  if (xe .EQ. -1) then       
3069                     res@trXMaxF            = maxiwsususodz
3070                  else
3071                     res@trXMaxF            = xe 
3072                  end if
3073                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3074               end if 
3075            end if
3076            if (vNam(varn) .EQ. "wspsodz" .OR. vNam(varn) .EQ. "w*p*:dz") then
3077               miniwspsodz=min(data(varn,:,min_z_int:max_z_int))-min_value
3078               maxiwspsodz=max(data(varn,:,min_z_int:max_z_int))+max_value
3079               if (over .EQ. 1) then
3080                  res@xyDashPattern  = 1
3081                  plot_wspsodz = gsn_csm_xy(wks,data(varn,:,:),z,res)
3082               else
3083                  res@gsnLeftString      = vNam(varn)
3084                  res@tiXAxisString      = "("+unit(varn)+")"
3085                  res@gsnRightString     = " "
3086                  if (xs .EQ. -1) then
3087                     res@trXMinF            = miniwspsodz
3088                  else
3089                     res@trXMinF            = xs
3090                  end if
3091                  if (xe .EQ. -1) then       
3092                     res@trXMaxF            = maxiwspsodz
3093                  else
3094                     res@trXMaxF            = xe 
3095                  end if
3096                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3097               end if
3098            end if
3099            if (vNam(varn) .EQ. "wpeodz" .OR. \
3100                vNam(varn) .EQ. "w"+dq+"p:dz") then
3101               miniwpeodz=min(data(varn,:,min_z_int:max_z_int))-min_value
3102               maxiwpeodz=max(data(varn,:,min_z_int:max_z_int))+max_value
3103               if (over .EQ. 1) then
3104                  res@xyDashPattern  = 2
3105                  plot_wpeodz = gsn_csm_xy(wks,data(varn,:,:),z,res)
3106               else
3107                  res@gsnLeftString      = vNam(varn)
3108                  res@tiXAxisString      = "("+unit(varn)+")"
3109                  res@gsnRightString     = " "
3110                  if (xs .EQ. -1) then
3111                     res@trXMinF            = miniwpeodz
3112                  else
3113                     res@trXMinF            = xs
3114                  end if
3115                  if (xe .EQ. -1) then       
3116                     res@trXMaxF            = maxiwpeodz
3117                  else
3118                     res@trXMaxF            = xe 
3119                  end if
3120                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3121               end if
3122            end if
3123         if (no_files .GT. 1) then
3124            print("nof="+nof+" und n="+n)
3125            multi_plot(nof,n)=plot(n)
3126            max_nof(nof,n)=max(data(varn,:,min_z_int:max_z_int))
3127            min_nof(nof,n)=min(data(varn,:,min_z_int:max_z_int))
3128            name(nof,n)   =vNam(varn)
3129            unit_(nof,n)  =unit(varn)
3130         end if
3131         if (over .EQ. 0) then
3132            n=n+1 
3133         end if 
3134         if (prof3d .EQ. 0)then   
3135            varn=varn+1
3136         end if
3137         delete(temp)
3138         delete(temp_att)
3139      end if         
3140   end do
3141   if (no_files .GT. 1) then
3142      delete(vNam)
3143      delete(files)
3144   end if
3145   
3146   end do
3147   ;#########ENDE DO LOOP FOR no_files GT 1#############
3148
3149   if (isStrSubset(data@long_name," SR " ) .and. over_remind) then
3150      print(" ")
3151      print("If you have outputs of statistic regions you cannot overlay "+\
3152            "variables;")
3153      print("'over' is set to 0" )
3154      print(" ")
3155      over = 0
3156   end if
3157
3158   if (count_var .EQ. 0) then
3159      print(" ")
3160      print("The variables 'var="+var+"'" )
3161      print("do not exist on your input file;")
3162      print("be sure to have one comma before and after each variable")
3163      print(" ")
3164      exit       
3165   end if
3166   
3167   if (no_files .GT. 1) then
3168      multi_legend=new(6,string)
3169      string_len=new(6,integer)
3170      multi_dash=new(no_files,string)
3171      multi_legend(0)="  "+name_legend_1
3172      string_len(0)=strlen(multi_legend(0))
3173      multi_legend(1)="  "+name_legend_2
3174      string_len(1)=strlen(multi_legend(1))
3175      multi_legend(2)="  "+name_legend_3
3176      string_len(2)=strlen(multi_legend(2))
3177      multi_legend(3)="  "+name_legend_4
3178      string_len(3)=strlen(multi_legend(3))
3179      multi_legend(4)="  "+name_legend_5
3180      string_len(4)=strlen(multi_legend(4))
3181      multi_legend(5)="  "+name_legend_6
3182      string_len(5)=strlen(multi_legend(5))
3183      do ml=1,no_files
3184         multi_dash(ml-1)=ml-1
3185      end do
3186      delete(plot)
3187      plot = new(dim,graphic)
3188      do pl=0,n-1
3189         plot0 = new(1,graphic)
3190         res@trXMinF = min(min_nof(:,pl))
3191         res@trXMaxF = max(max_nof(:,pl))
3192         res@gsnLeftString  = name(0,pl)
3193         res@gsnRightString = unit_(0,pl)
3194         res@tiXAxisString  = "("+unit_(0,pl)+")"
3195         
3196         data_0(:,:) = min(min_nof(:,pl))
3197         plot0 = gsn_csm_xy(wks,data_0(:,:),z_(pl,:),res)
3198
3199         ; ***************************************************
3200         ; legend for combined plot
3201         ; ***************************************************
3202
3203         lgres                    = True
3204         lgMonoDashIndex          = False
3205         lgres@lgLabelFont        = "helvetica"   
3206         lgres@lgLabelFontHeightF = font_size_legend*10.0       
3207         lgres@vpWidthF           = max(string_len)*0.015           
3208         lgres@vpHeightF          = 0.03*no_files         
3209         lgres@lgDashIndexes      = multi_dash(no_files-1:0)
3210         lbid = gsn_create_legend(\
3211                            wks,no_files,multi_legend(no_files-1:0),lgres)
3212
3213         amres = True
3214         amres@amParallelPosF   = max(string_len)*0.012+0.78               
3215         amres@amOrthogonalPosF = -0.0315*no_files+0.431         
3216         annoid1 = gsn_add_annotation(plot0,lbid,amres)
3217
3218         do plo=0,no_files-1
3219            overlay(plot0,multi_plot(plo,pl))
3220            plot(pl)=plot0
3221         end do
3222         delete(plot0)
3223      end do
3224   end if
3225
3226   if (count_var .EQ. 0) then
3227      print(" ")
3228      print("Select a variable 'var=' or use the default value")
3229      print(" ")
3230      print("Your selection '"+var+"' does not exist on the input file")
3231      print(" ")
3232      exit
3233   end if
3234
3235   if (over .EQ. 1 ) then
3236
3237      overlay(plot_u,plot_v)
3238      overlay(plot_u,plot_w)
3239      u=0
3240      overlay(plot_pt,plot_vpt)
3241      overlay(plot_pt,plot_lpt)
3242      pt=0
3243      overlay(plot_q,plot_qv)
3244      overlay(plot_q,plot_ql)
3245      q=0
3246      overlay(plot_e,plot_es)
3247      e=0
3248      overlay(plot_km,plot_kh)
3249      km=0
3250      overlay(plot_wpup,plot_wsus)
3251      overlay(plot_wpup,plot_wu)
3252      wpup=0
3253      overlay(plot_wpvp,plot_wsvs)
3254      overlay(plot_wpvp,plot_wv)
3255      wpvp=0
3256      overlay(plot_wpptp,plot_wspts)
3257      overlay(plot_wpptp,plot_wpt)
3258      wpptp=0
3259      overlay(plot_wsptsBC,plot_wptBC)
3260      wsptsBC=0
3261      overlay(plot_wpvptp,plot_wsvpts)
3262      overlay(plot_wpvptp,plot_wvpt)
3263      wpvptp=0
3264      overlay(plot_wpqp,plot_wsqs)
3265      overlay(plot_wpqp,plot_wq)
3266      wpqp=0
3267      overlay(plot_wpqvp,plot_wsqvs)
3268      overlay(plot_wpqvp,plot_wqv)
3269      wpqvp=0
3270      overlay(plot_wpsp,plot_wsss)
3271      overlay(plot_wpsp,plot_ws)
3272      wpsp=0
3273      overlay(plot_wpsap,plot_wssas)
3274      overlay(plot_wpsap,plot_wsa)
3275      wpsap=0
3276      overlay(plot_us2,plot_vs2)
3277      overlay(plot_us2,plot_ws2)
3278      us2=0
3279      overlay(plot_wsususodz,plot_wspsodz)
3280      overlay(plot_wsususodz,plot_wpeodz)
3281      wsususodz=0
3282
3283   end if
3284
3285   if (over .EQ. 1) then
3286   
3287      do varn = 0,dim-1   
3288         
3289         check = True
3290     
3291         if (prof3d .EQ. 0) then
3292            if ( isStrSubset( vNam(varn), "time") .OR. \
3293                 isStrSubset( vNam(varn), "NORM")) then
3294               check = False
3295            end if
3296         else
3297            if ( isStrSubset( vNam(varn), "time") .OR.  \
3298                 isStrSubset( vNam(varn), "zusi") .OR.  \
3299                 isStrSubset( vNam(varn), "zwwi") .OR.  \
3300                 isStrSubset( vNam(varn), "x") .OR.     \
3301                 isStrSubset( vNam(varn), "xu") .OR.    \
3302                 isStrSubset( vNam(varn), "y") .OR.     \
3303                 isStrSubset( vNam(varn), "yv") .OR.    \
3304                 isStrSubset( vNam(varn), "zu_3d") .OR. \
3305                 isStrSubset( vNam(varn), "zw_3d")) then
3306               check = False
3307            end if
3308         end if
3309
3310         if (var .NE. "all") then     
3311            check = isStrSubset( var,","+vNam(varn)+"," )
3312         end if     
3313
3314         if (check)then
3315
3316            if (prof3d .EQ. 0) then
3317               if (log_z .EQ. 1) then
3318                  z = f_att->$vNam(varn+1)$(1:dimz-1)
3319               else
3320                  z = f_att->$vNam(varn+1)$               
3321               end if
3322            else
3323               do i=0,b-1           
3324                  if (isStrSubset( a(i),"zu_3d" ))then
3325                     z_v(varn,:) = z_u
3326                     if (log_z .EQ. 1) then
3327                        z = z_v(varn,1:dimz-1)
3328                     else
3329                        z = z_v(varn,:)
3330                     end if
3331                  else
3332                     if (isStrSubset( a(i),"zw_3d" ))then
3333                        z_v(varn,:) = z_w
3334                        if (log_z .EQ. 1) then
3335                           z = z_v(varn,1:dimz-1)
3336                        else
3337                           z = z_v(varn,:)
3338                        end if
3339                     end if                   
3340                  end if
3341               end do           
3342            end if
3343
3344            z=z/norm_z
3345           
3346            if (abs(min(data(varn,:,min_z_int:max_z_int))) .GT. 10)then
3347               min_value = abs(0.01*min(data(varn,:,min_z_int:max_z_int)))
3348               max_value = abs(0.01*max(data(varn,:,min_z_int:max_z_int)))
3349            else
3350               if (abs(min(data(varn,:,min_z_int:max_z_int))) .LT. 0.01 .AND. \
3351                   abs(max(data(varn,:,min_z_int:max_z_int))) .GT. 0.01)then
3352                  min_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
3353                  max_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
3354               else
3355                  if (abs(max(data(varn,:,:))) .LT. 0.01 .AND. \
3356                      abs(min(data(varn,:,min_z_int:max_z_int))) .GT. 0.01)then
3357                     min_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
3358                     max_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
3359                  else
3360                     min_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
3361                     max_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
3362                  end if
3363               end if
3364            end if
3365            if (min(data(varn,:,min_z_int:max_z_int)) .EQ. 0 .AND.\
3366                max(data(varn,:,min_z_int:max_z_int)) .EQ. 0)then
3367               min_value = 0.1
3368               max_value = 0.1
3369            end if
3370
3371            res@gsnLeftString      = vNam(varn)
3372            res@tiXAxisString      = "("+unit(varn)+")"
3373            res@gsnRightString     = " "
3374            res@trYMinF            = min_z
3375            res@trYMaxF            = max_z 
3376            res@xyDashPattern = 0
3377
3378            if (xs .EQ. -1) then
3379               res@trXMinF = min(data(varn,:,min_z_int:max_z_int))-min_value
3380            else
3381               res@trXMinF = xs
3382            end if
3383            if (xe .EQ. -1) then
3384               res@trXMaxF = max(data(varn,:,min_z_int:max_z_int))+max_value
3385            else
3386               res@trXMaxF = xe 
3387            end if
3388            plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res) 
3389           
3390            if (vNam(varn) .EQ. "u" .OR. vNam(varn) .EQ. "v" .OR. \
3391                vNam(varn) .EQ. "w") then
3392               if (u .EQ. 0) then
3393                  res@gsnLeftString      = "u, v and w"
3394                  res@tiXAxisString      = "("+unit(varn)+")"
3395                  res@gsnRightString     = " "
3396                  if (xs .EQ. -1) then
3397                     res@trXMinF = min((/miniu,miniv,miniw/))
3398                  else
3399                     res@trXMinF = xs
3400                  end if
3401                  if (xe .EQ. -1) then
3402                     res@trXMaxF = max((/maxiu,maxiv,maxiw/))
3403                  else
3404                     res@trXMaxF = xe 
3405                  end if 
3406                  if (vNam(varn) .EQ. "v")then
3407                     res@xyDashPattern = 1
3408                  end if
3409                  if (vNam(varn) .EQ. "w")then
3410                     res@xyDashPattern = 2
3411                  end if
3412                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3413
3414                  ; ***************************************************
3415                  ; legend for overlaid plot
3416                  ; ***************************************************
3417     
3418                  lgres                    = True
3419                  lgMonoDashIndex          = False
3420                  lgres@lgLabelFont        = "helvetica"   
3421                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3422                  lgres@vpWidthF           = 0.07           
3423                  lgres@vpHeightF          = 0.12         
3424                  lgres@lgDashIndexes      = (/0,1,2/)
3425                  lbid = gsn_create_legend(wks,3,(/"u","v","w"/),lgres)       
3426
3427                  amres = True
3428                  amres@amParallelPosF   = 0.88             
3429                  amres@amOrthogonalPosF = 0.33           
3430                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3431                  overlay(plot(n),plot_u)
3432                  u=1
3433               else
3434                  if (prof3d .EQ. 0)then
3435                     varn=varn+1
3436                  end if
3437                  continue
3438               end if       
3439            end if 
3440     
3441            if (vNam(varn) .EQ. "pt" .OR. vNam(varn) .EQ. "vpt" .OR. \
3442                vNam(varn) .EQ. "lpt") then
3443               if (pt .EQ. 0) then
3444                  res@gsnLeftString      = "pt, vpt and lpt"
3445                  res@tiXAxisString      = "("+unit(varn)+")"
3446                  res@gsnRightString     = " "
3447                  if (xs .EQ. -1) then
3448                     res@trXMinF = min((/minipt,minivpt,minilpt/))
3449                  else
3450                     res@trXMinF = xs
3451                  end if
3452                  if (xe .EQ. -1) then
3453                     res@trXMaxF = max((/maxipt,maxivpt,maxilpt/))
3454                  else
3455                     res@trXMaxF = xe 
3456                  end if
3457                  if (vNam(varn) .EQ. "vpt")then
3458                     res@xyDashPattern = 1
3459                  end if
3460                  if (vNam(varn) .EQ. "lpt")then
3461                     res@xyDashPattern = 2
3462                  end if
3463                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3464
3465                  ; ***************************************************
3466                  ; legend for overlaid plot
3467                  ; ***************************************************
3468     
3469                  lgres                    = True
3470                  lgMonoDashIndex          = False
3471                  lgres@lgLabelFont        = "helvetica"   
3472                  lgres@lgLabelFontHeightF = font_size_legend*10.0         
3473                  lgres@vpWidthF           = 0.07           
3474                  lgres@vpHeightF          = 0.12         
3475                  lgres@lgDashIndexes      = (/0,1,2/)
3476                  lbid = gsn_create_legend(wks,3,(/"pt","vpt","lpt"/),lgres)
3477                  amres = True
3478                  amres@amParallelPosF   = 0.88     
3479                  amres@amOrthogonalPosF = 0.33           
3480                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3481                  overlay(plot(n),plot_pt)
3482                  pt=1
3483               else
3484                  if (prof3d .EQ. 0)then
3485                     varn=varn+1
3486                  end if
3487                  continue       
3488               end if
3489            end if           
3490            if (vNam(varn) .EQ. "q" .OR. vNam(varn) .EQ. "qv" .OR. \
3491                vNam(varn) .EQ. "ql") then
3492               if (q .EQ. 0) then
3493                  res@gsnLeftString      = "q, qv and ql"
3494                  res@tiXAxisString      = "("+unit(varn)+")"
3495                  res@gsnRightString     = " "
3496                  if (xs .EQ. -1) then
3497                     res@trXMinF = min((/miniq,miniqv,miniql/))
3498                  else
3499                     res@trXMinF = xs
3500                  end if
3501                  if (xe .EQ. -1) then
3502                     res@trXMaxF = max((/maxiq,maxiqv,maxiql/))
3503                  else
3504                     res@trXMaxF = xe 
3505                  end if
3506                  if (vNam(varn) .EQ. "qv")then
3507                     res@xyDashPattern = 1
3508                  end if
3509                  if (vNam(varn) .EQ. "ql")then
3510                     res@xyDashPattern = 2
3511                  end if
3512                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3513
3514                  ; ***************************************************
3515                  ; legend for overlaid plot
3516                  ; ***************************************************
3517     
3518                  lgres                    = True
3519                  lgMonoDashIndex          = False
3520                  lgres@lgLabelFont        = "helvetica"   
3521                  lgres@lgLabelFontHeightF = font_size_legend*10.0         
3522                  lgres@vpWidthF           = 0.07         
3523                  lgres@vpHeightF          = 0.12         
3524                  lgres@lgDashIndexes      = (/0,1,2/)
3525                  lbid = gsn_create_legend(wks,3,(/"q","qv","ql"/),lgres)
3526
3527                  amres = True
3528                  amres@amParallelPosF   = 0.88             
3529                  amres@amOrthogonalPosF = 0.33           
3530                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3531                  overlay(plot(n),plot_q)
3532                  q=1
3533               else
3534                  if (prof3d .EQ. 0)then
3535                     varn=varn+1
3536                  end if
3537                  continue   
3538               end if
3539            end if   
3540           
3541            if (vNam(varn) .EQ. "e" .OR. vNam(varn) .EQ. "es" .OR. \
3542                vNam(varn) .EQ. "e*" ) then
3543               if (e .EQ. 0) then
3544                  res@gsnLeftString      = "e and e*"
3545                  res@tiXAxisString      = "("+unit(varn)+")"
3546                  res@gsnRightString     = " "
3547                  if (xs .EQ. -1) then
3548                     res@trXMinF = min((/minie,minies/))
3549                  else
3550                     res@trXMinF = xs
3551                  end if
3552                  if (xe .EQ. -1) then
3553                     res@trXMaxF = max((/maxie,maxies/))
3554                  else
3555                     res@trXMaxF = xe 
3556                  end if
3557                  if (vNam(varn) .EQ. "es")then
3558                     res@xyDashPattern = 1
3559                  end if
3560                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3561
3562                  ; ***************************************************
3563                  ; legend for overlaid plot
3564                  ; ***************************************************
3565     
3566                  lgres                    = True
3567                  lgMonoDashIndex          = False
3568                  lgres@lgLabelFont        = "helvetica"   
3569                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3570                  lgres@vpWidthF           = 0.07           
3571                  lgres@vpHeightF          = 0.08         
3572                  lgres@lgDashIndexes      = (/0,1,2/)
3573                  lbid = gsn_create_legend(wks,2,(/"e","e*"/),lgres)       
3574
3575                  amres = True
3576                  amres@amParallelPosF   = 0.88             
3577                  amres@amOrthogonalPosF = 0.365           
3578                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3579                  overlay(plot(n),plot_e)
3580                  e=1
3581               else
3582                  if (prof3d .EQ. 0)then
3583                     varn=varn+1
3584                  end if
3585                  continue   
3586               end if
3587            end if           
3588            if (vNam(varn) .EQ. "km" .OR. vNam(varn) .EQ. "kh") then
3589               if (km .EQ. 0) then
3590                  res@gsnLeftString      = "km and kh"
3591                  res@tiXAxisString      = "("+unit(varn)+")"
3592                  res@gsnRightString     = " "
3593                  if (xs .EQ. -1) then
3594                     res@trXMinF = min((/minikm,minikh/))
3595                  else
3596                     res@trXMinF = xs
3597                  end if
3598                  if (xe .EQ. -1) then
3599                     res@trXMaxF = max((/maxikm,maxikh/))
3600                  else
3601                     res@trXMaxF = xe 
3602                  end if
3603                  if (vNam(varn) .EQ. "kh")then
3604                     res@xyDashPattern = 1
3605                  end if
3606                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3607
3608                  ; ***************************************************
3609                  ; legend for overlaid plot
3610                  ; ***************************************************
3611     
3612                  lgres                    = True
3613                  lgMonoDashIndex          = False
3614                  lgres@lgLabelFont        = "helvetica"   
3615                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3616                  lgres@vpWidthF           = 0.07           
3617                  lgres@vpHeightF          = 0.08         
3618                  lgres@lgDashIndexes      = (/0,1,2/)
3619                  lbid = gsn_create_legend(wks,2,(/"km","kh"/),lgres)       
3620
3621                  amres = True
3622                  amres@amParallelPosF   = 0.88             
3623                  amres@amOrthogonalPosF = 0.365           
3624                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3625                  overlay(plot(n),plot_km)
3626                  km=1
3627               else
3628                  if (prof3d .EQ. 0)then
3629                     varn=varn+1
3630                  end if
3631                  continue   
3632               end if
3633            end if           
3634           
3635            if (vNam(varn) .EQ. "wpup" .OR. vNam(varn) .EQ. "wsus" .OR.      \
3636                vNam(varn) .EQ. "wu" .OR. vNam(varn) .EQ. "w"+dq+"u"+dq .OR. \
3637                vNam(varn) .EQ. "w*u*") then
3638               if (wpup .EQ. 0) then
3639                  res@gsnLeftString      = "w"+dq+"u"+dq+", w*u* and wu"
3640                  res@tiXAxisString      = "("+unit(varn)+")"
3641                  res@gsnRightString     = " "
3642                  if (xs .EQ. -1) then
3643                     res@trXMinF = min((/miniwpup,miniwsus,miniwu/))
3644                  else
3645                     res@trXMinF = xs
3646                  end if
3647                  if (xe .EQ. -1) then
3648                     res@trXMaxF = max((/maxiwpup,maxiwsus,maxiwu/))
3649                  else
3650                     res@trXMaxF = xe 
3651                  end if 
3652                  if (vNam(varn) .EQ. "wsus")then
3653                     res@xyDashPattern = 1
3654                  end if
3655                  if (vNam(varn) .EQ. "wu")then
3656                     res@xyDashPattern = 2
3657                  end if
3658                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3659
3660                  ; ***************************************************
3661                  ; legend for overlaid plot
3662                  ; ***************************************************
3663     
3664                  lgres                    = True
3665                  lgMonoDashIndex          = False
3666                  lgres@lgLabelFont        = "helvetica"   
3667                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3668                  lgres@vpWidthF           = 0.08           
3669                  lgres@vpHeightF          = 0.12         
3670                  lgres@lgDashIndexes      = (/0,1,2/)
3671                  lbid = gsn_create_legend(\
3672                                wks,3,(/"w"+dq+"u"+dq,"w*u*","wu"/),lgres)
3673
3674                  amres = True
3675                  amres@amParallelPosF   = 0.88             
3676                  amres@amOrthogonalPosF = 0.33           
3677                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3678                  overlay(plot(n),plot_wpup)
3679                  wpup=1
3680               else
3681                  if (prof3d .EQ. 0)then
3682                     varn=varn+1
3683                  end if
3684                  continue   
3685               end if
3686            end if
3687            if (vNam(varn) .EQ. "wpvp" .OR. vNam(varn) .EQ. "wsvs" .OR.     \
3688               vNam(varn) .EQ. "wv" .OR. vNam(varn) .EQ. "w"+dq+"v"+dq .OR. \
3689               vNam(varn) .EQ. "w*v*") then
3690               if (wpvp .EQ. 0) then
3691                  res@gsnLeftString      = "w"+dq+"v"+dq+", w*v* and wv"
3692                  res@tiXAxisString      = "("+unit(varn)+")"
3693                  res@gsnRightString     = " "
3694                  if (xs .EQ. -1) then
3695                     res@trXMinF = min((/miniwpvp,miniwsvs,miniwv/))
3696                  else
3697                     res@trXMinF = xs
3698                  end if
3699                  if (xe .EQ. -1) then
3700                     res@trXMaxF = max((/maxiwpvp,maxiwsvs,maxiwv/))
3701                  else
3702                     res@trXMaxF = xe 
3703                  end if
3704                  if (vNam(varn) .EQ. "wsvs")then
3705                     res@xyDashPattern = 1
3706                  end if
3707                  if (vNam(varn) .EQ. "wv")then
3708                     res@xyDashPattern = 2
3709                  end if
3710                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3711
3712                  ; ***************************************************
3713                  ; legend for overlaid plot
3714                  ; ***************************************************
3715     
3716                  lgres                    = True
3717                  lgMonoDashIndex          = False
3718                  lgres@lgLabelFont        = "helvetica"   
3719                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3720                  lgres@vpWidthF           = 0.08         
3721                  lgres@vpHeightF          = 0.12         
3722                  lgres@lgDashIndexes      = (/0,1,2/)
3723                  lbid = gsn_create_legend(\
3724                             wks,3,(/"w"+dq+"v"+dq,"w*v*","wv"/),lgres)       
3725
3726                  amres = True
3727                  amres@amParallelPosF   = 0.88             
3728                  amres@amOrthogonalPosF = 0.33           
3729                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3730                  overlay(plot(n),plot_wpvp)
3731                  wpvp=1
3732               else
3733                  if (prof3d .EQ. 0)then
3734                     varn=varn+1
3735                  end if
3736                  continue   
3737               end if
3738            end if
3739            if (vNam(varn) .EQ. "wpptp" .OR. vNam(varn) .EQ. "wspts" .OR.     \
3740                vNam(varn) .EQ. "wpt" .OR. vNam(varn) .EQ. "w"+dq+"pt"+dq .OR.\
3741                vNam(varn) .EQ. "w*pt*") then
3742               if (wpptp .EQ. 0) then                 
3743                  res@gsnLeftString      = "w"+dq+"pt"+dq+", w*pt* and wpt"
3744                  res@tiXAxisString      = "("+unit(varn)+")"
3745                  res@gsnRightString     = " "
3746                  if (xs .EQ. -1) then
3747                     res@trXMinF = min((/miniwpptp,miniwspts,miniwpt/))
3748                  else
3749                     res@trXMinF = xs
3750                  end if
3751                  if (xe .EQ. -1) then
3752                     res@trXMaxF = max((/maxiwpptp,maxiwspts,maxiwpt/))
3753                  else
3754                     res@trXMaxF = xe 
3755                  end if
3756                  if (vNam(varn) .EQ. "wspts" .OR. vNam(varn) .EQ. "w*pt*")then
3757                     res@xyDashPattern = 1
3758                  end if
3759                  if (vNam(varn) .EQ. "wpt")then
3760                     res@xyDashPattern = 2
3761                  end if
3762                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3763
3764                  ; ***************************************************
3765                  ; legend for overlaid plot
3766                  ; ***************************************************
3767     
3768                  lgres                    = True
3769                  lgMonoDashIndex          = False
3770                  lgres@lgLabelFont        = "helvetica"   
3771                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3772                  lgres@vpWidthF           = 0.09           
3773                  lgres@vpHeightF          = 0.12         
3774                  lgres@lgDashIndexes      = (/0,1,2/)             
3775                  lbid = gsn_create_legend(\
3776                               wks,3,(/"w"+dq+"pt"+dq,"w*pt*","wpt"/),lgres)
3777
3778                  amres = True
3779                  amres@amParallelPosF   = 0.88             
3780                  amres@amOrthogonalPosF = 0.33           
3781                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3782                  overlay(plot(n),plot_wpptp)
3783                  wpptp=1
3784               else
3785                  if (prof3d .EQ. 0)then
3786                     varn=varn+1
3787                  end if
3788                  continue   
3789               end if
3790            end if
3791            if (vNam(varn) .EQ. "wsptsBC" .OR. vNam(varn) .EQ. "wptBC" .OR.\
3792                vNam(varn) .EQ. "w*pt*BC") then
3793               if (wsptsBC .EQ. 0) then
3794                  res@gsnLeftString      = "w*pt*BC and wptBC"
3795                  res@tiXAxisString      = "("+unit(varn)+")"
3796                  res@gsnRightString     = " "
3797                  if (xs .EQ. -1) then
3798                     res@trXMinF = min((/miniwsptsBC,miniwptBC/))
3799                  else
3800                     res@trXMinF = xs
3801                  end if
3802                  if (xe .EQ. -1) then
3803                     res@trXMaxF = max((/maxiwsptsBC,maxiwptBC/))
3804                  else
3805                     res@trXMaxF = xe 
3806                  end if
3807                  if (vNam(varn) .EQ. "wptBC")then
3808                     res@xyDashPattern = 1
3809                  end if
3810                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3811
3812                  ; ***************************************************
3813                  ; legend for overlaid plot
3814                  ; ***************************************************
3815     
3816                  lgres                    = True
3817                  lgMonoDashIndex          = False
3818                  lgres@lgLabelFont        = "helvetica"   
3819                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3820                  lgres@vpWidthF           = 0.1           
3821                  lgres@vpHeightF          = 0.12         
3822                  lgres@lgDashIndexes      = (/0,1,2/)
3823                  lbid = gsn_create_legend(\
3824                                        wks,3,(/"w*pt*BC","wptBC"/),lgres)
3825
3826                  amres = True
3827                  amres@amParallelPosF   = 0.88             
3828                  amres@amOrthogonalPosF = 0.33           
3829                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3830                  overlay(plot(n),plot_wsptsBC)
3831                  wsptsBC=1
3832               else
3833                  if (prof3d .EQ. 0)then
3834                     varn=varn+1
3835                  end if
3836                  continue   
3837               end if 
3838            end if             
3839            if (vNam(varn) .EQ. "wpvptp" .OR. vNam(varn) .EQ. "wsvpts" .OR. \
3840                vNam(varn) .EQ. "wvpt" .OR. vNam(varn) .EQ. \
3841                "w"+dq+"vpt"+dq .OR. vNam(varn) .EQ. "w*vpt*") then
3842               if (wpvptp .EQ. 0) then
3843                  res@gsnLeftString      = "w"+dq+"vpt"+dq+", w*vpt* and wvpt"
3844                  res@tiXAxisString      = "("+unit(varn)+")"
3845                  res@gsnRightString     = " "
3846                  if (xs .EQ. -1) then
3847                     res@trXMinF = min((/miniwpvptp,miniwsvpts,miniwvpt/))
3848                  else
3849                     res@trXMinF = xs
3850                  end if
3851                  if (xe .EQ. -1) then
3852                     res@trXMaxF = max((/maxiwpvptp,maxiwsvpts,maxiwvpt/))
3853                  else
3854                     res@trXMaxF = xe 
3855                  end if
3856                  if (vNam(varn) .EQ. "wsvpts")then
3857                     res@xyDashPattern = 1
3858                  end if
3859                  if (vNam(varn) .EQ. "wvpt")then
3860                     res@xyDashPattern = 2
3861                  end if
3862                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3863
3864                  ; ***************************************************
3865                  ; legend for overlaid plot
3866                  ; ***************************************************
3867     
3868                  lgres                    = True
3869                  lgMonoDashIndex          = False
3870                  lgres@lgLabelFont        = "helvetica"   
3871                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3872                  lgres@vpWidthF           = 0.1           
3873                  lgres@vpHeightF          = 0.12         
3874                  lgres@lgDashIndexes      = (/0,1,2/)
3875                  lbid = gsn_create_legend(\
3876                             wks,3,(/"w"+dq+"vpt"+dq,"w*vpt*","wvpt"/),lgres)
3877                  amres = True
3878                  amres@amParallelPosF   = 0.88             
3879                  amres@amOrthogonalPosF = 0.33           
3880                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3881                  overlay(plot(n),plot_wpvptp)
3882                  wpvptp=1
3883               else
3884                  if (prof3d .EQ. 0)then
3885                     varn=varn+1
3886                  end if
3887                  continue   
3888               end if
3889            end if
3890            if (vNam(varn) .EQ. "wpqp" .OR. vNam(varn) .EQ. "wsqs" .OR.      \
3891                vNam(varn) .EQ. "wq" .OR. vNam(varn) .EQ. "w"+dq+"q"+dq .OR. \
3892                vNam(varn) .EQ. "w*q*") then
3893               if (wpqp .EQ. 0) then
3894                  res@gsnLeftString      = "w"+dq+"q"+dq+", w*q* and wq"
3895                  res@tiXAxisString      = "("+unit(varn)+")"
3896                  res@gsnRightString     = " "
3897                  if (xs .EQ. -1) then
3898                     res@trXMinF = min((/miniwpqp,miniwsqs,miniwq/))
3899                  else
3900                     res@trXMinF = xs
3901                  end if
3902                  if (xe .EQ. -1) then
3903                     res@trXMaxF = max((/maxiwpqp,maxiwsqs,maxiwq/))
3904                  else
3905                     res@trXMaxF = xe 
3906                  end if 
3907                  if (vNam(varn) .EQ. "wsqs")then
3908                     res@xyDashPattern = 1
3909                  end if
3910                  if (vNam(varn) .EQ. "wq")then
3911                     res@xyDashPattern = 2
3912                  end if
3913                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3914
3915                  ; ***************************************************
3916                  ; legend for overlaid plot
3917                  ; ***************************************************
3918     
3919                  lgres                    = True
3920                  lgMonoDashIndex          = False
3921                  lgres@lgLabelFont        = "helvetica"   
3922                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3923                  lgres@vpWidthF           = 0.08           
3924                  lgres@vpHeightF          = 0.12         
3925                  lgres@lgDashIndexes      = (/0,1,2/)
3926                  lbid = gsn_create_legend(\
3927                             wks,3,(/"w"+dq+"q"+dq,"w*q*","wq"/),lgres)       
3928
3929                  amres = True
3930                  amres@amParallelPosF   = 0.88             
3931                  amres@amOrthogonalPosF = 0.33           
3932                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3933                  overlay(plot(n),plot_wpqp)
3934                  wpqp=1
3935               else
3936                  if (prof3d .EQ. 0)then
3937                     varn=varn+1
3938                  end if
3939                  continue   
3940               end if
3941            end if
3942            if (vNam(varn) .EQ. "wpqvp" .OR. vNam(varn) .EQ. "wsqvs" .OR.     \
3943                vNam(varn) .EQ. "wqv" .OR. vNam(varn) .EQ. "w"+dq+"qv"+dq .OR.\
3944                vNam(varn) .EQ. "w*qv*") then
3945               if (wpqvp .EQ. 0) then
3946                  res@gsnLeftString      ="w"+dq+"qv"+dq+" , w*qv* and wqv"
3947                  res@tiXAxisString      = "("+unit(varn)+")"
3948                  res@gsnRightString     = " "
3949                  if (xs .EQ. -1) then
3950                     res@trXMinF = min((/miniwpqp,miniwsqvs,miniwqv/))
3951                  else
3952                     res@trXMinF = xs
3953                  end if
3954                  if (xe .EQ. -1) then
3955                     res@trXMaxF = max((/maxiwpqp,maxiwsqvs,maxiwqv/))
3956                  else
3957                     res@trXMaxF = xe 
3958                  end if
3959                  if (vNam(varn) .EQ. "wsqvs")then
3960                     res@xyDashPattern = 1
3961                  end if
3962                  if (vNam(varn) .EQ. "wqv")then
3963                     res@xyDashPattern = 2
3964                  end if
3965                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3966
3967                  ; ***************************************************
3968                  ; legend for overlaid plot
3969                  ; ***************************************************
3970     
3971                  lgres                    = True
3972                  lgMonoDashIndex          = False
3973                  lgres@lgLabelFont        = "helvetica"   
3974                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
3975                  lgres@vpWidthF           = 0.09           
3976                  lgres@vpHeightF          = 0.12         
3977                  lgres@lgDashIndexes      = (/0,1,2/)
3978                  lbid = gsn_create_legend(\
3979                                wks,3,(/"w"+dq+"qv"+dq,"w*qv*","wqv"/),lgres)
3980
3981                  amres = True
3982                  amres@amParallelPosF   = 0.88             
3983                  amres@amOrthogonalPosF = 0.33           
3984                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3985                  overlay(plot(n),plot_wpqvp)
3986                  wpqvp=1
3987               else
3988                  if (prof3d .EQ. 0)then
3989                     varn=varn+1
3990                  end if
3991                  continue   
3992               end if
3993            end if
3994            if (vNam(varn) .EQ. "wpsp" .OR. vNam(varn) .EQ. "wsss" .OR.     \
3995                vNam(varn) .EQ. "ws" .OR. vNam(varn) .EQ. "w"+dq+"s"+dq .OR.\
3996                vNam(varn) .EQ. "w*s*") then
3997               if (wpsp .EQ. 0) then
3998                  res@gsnLeftString      = "w"+dq+"s"+dq+", w*s* and ws"
3999                  res@tiXAxisString      = "("+unit(varn)+")"
4000                  res@gsnRightString     = " "
4001                  if (xs .EQ. -1) then
4002                     res@trXMinF = min((/miniwpsp,miniwsss,miniws/))
4003                  else
4004                     res@trXMinF = xs
4005                  end if
4006                  if (xe .EQ. -1) then
4007                     res@trXMaxF = max((/maxiwpsp,maxiwsss,maxiws/))
4008                  else
4009                     res@trXMaxF = xe 
4010                  end if
4011                  if (vNam(varn) .EQ. "wsss")then
4012                     res@xyDashPattern = 1
4013                  end if
4014                  if (vNam(varn) .EQ. "ws")then
4015                     res@xyDashPattern = 2
4016                  end if
4017                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4018
4019                  ; ***************************************************
4020                  ; legend for overlaid plot
4021                  ; ***************************************************
4022     
4023                  lgres                    = True
4024                  lgMonoDashIndex          = False
4025                  lgres@lgLabelFont        = "helvetica"   
4026                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
4027                  lgres@vpWidthF           = 0.08           
4028                  lgres@vpHeightF          = 0.12         
4029                  lgres@lgDashIndexes      = (/0,1,2/)
4030                  lbid = gsn_create_legend(\
4031                           wks,3,(/"w"+dq+"s"+dq,"w*s*","ws"/),lgres)       
4032
4033                  amres = True
4034                  amres@amParallelPosF   = 0.88             
4035                  amres@amOrthogonalPosF = 0.33           
4036                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4037                  overlay(plot(n),plot_wpsp)
4038                  wpsp=1
4039               else
4040                  if (prof3d .EQ. 0)then
4041                     varn=varn+1
4042                  end if
4043                  continue   
4044               end if
4045            end if
4046            if (vNam(varn) .EQ. "wpsap" .OR.vNam(varn) .EQ. "wssas" .OR.      \
4047                vNam(varn) .EQ. "wsa" .OR. vNam(varn) .EQ. "w"+dq+"sa"+dq .OR.\
4048                vNam(varn) .EQ. "w*sa*") then
4049               if (wpsap .EQ. 0) then
4050                  res@gsnLeftString      = "w"+dq+"sa"+dq+", w*sa* and wsa"
4051                  res@tiXAxisString      = "("+unit(varn)+")"
4052                  res@gsnRightString     = " "
4053                  if (xs .EQ. -1) then
4054                     res@trXMinF = min((/miniwpsap,miniwssas,miniwsa/))
4055                  else
4056                     res@trXMinF = xs
4057                  end if
4058                  if (xe .EQ. -1) then
4059                     res@trXMaxF = max((/maxiwpsap,maxiwssas,maxiwsa/))
4060                  else
4061                     res@trXMaxF = xe 
4062                  end if
4063                  if (vNam(varn) .EQ. "wssas")then
4064                     res@xyDashPattern = 1
4065                  end if
4066                  if (vNam(varn) .EQ. "wsa")then
4067                     res@xyDashPattern = 2
4068                  end if
4069                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4070
4071                  ; ***************************************************
4072                  ; legend for overlaid plot
4073                  ; ***************************************************
4074     
4075                  lgres                    = True
4076                  lgMonoDashIndex          = False
4077                  lgres@lgLabelFont        = "helvetica"   
4078                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
4079                  lgres@vpWidthF           = 0.09           
4080                  lgres@vpHeightF          = 0.12         
4081                  lgres@lgDashIndexes      = (/0,1,2/)
4082                  lbid = gsn_create_legend(\
4083                             wks,3,(/"w"+dq+"sa"+dq,"w*sa*","wsa"/),lgres)
4084                  amres = True
4085                  amres@amParallelPosF   = 0.88             
4086                  amres@amOrthogonalPosF = 0.33           
4087                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4088                  overlay(plot(n),plot_wpsap)
4089                  wpsap=1
4090               else
4091                  if (prof3d .EQ. 0)then
4092                     varn=varn+1
4093                  end if
4094                  continue   
4095               end if
4096            end if
4097         
4098            if (vNam(varn) .EQ. "us2" .OR. vNam(varn) .EQ. "vs2" .OR. \
4099                vNam(varn) .EQ. "ws2" .OR. vNam(varn) .EQ. "u*2" .OR. \
4100                vNam(varn) .EQ. "v*2" .OR. vNam(varn) .EQ. "w*2" ) then
4101               if (us2 .EQ. 0) then
4102                  res@gsnLeftString      = "u*2, v*2 and w*2"
4103                  res@tiXAxisString      = "("+unit(varn)+")"
4104                  res@gsnRightString     = " "
4105                  if (xs .EQ. -1) then
4106                     res@trXMinF = min((/minius2,minivs2,miniws2/))
4107                  else
4108                     res@trXMinF = xs
4109                  end if
4110                  if (xe .EQ. -1) then
4111                     res@trXMaxF = max((/maxius2,maxivs2,maxiws2/))
4112                  else
4113                     res@trXMaxF = xe 
4114                  end if
4115                  if (vNam(varn) .EQ. "vs2")then
4116                     res@xyDashPattern = 1
4117                  end if
4118                  if (vNam(varn) .EQ. "ws2")then
4119                     res@xyDashPattern = 2
4120                  end if
4121                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4122
4123                  ; ***************************************************
4124                  ; legend for overlaid plot
4125                  ; ***************************************************
4126     
4127                  lgres                    = True
4128                  lgMonoDashIndex          = False
4129                  lgres@lgLabelFont        = "helvetica"   
4130                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
4131                  lgres@vpWidthF           = 0.07           
4132                  lgres@vpHeightF          = 0.12         
4133                  lgres@lgDashIndexes      = (/0,1,2/)
4134                  lbid = gsn_create_legend(wks,3,(/"u*2","v*2","w*2"/),lgres)
4135                  amres = True
4136                  amres@amParallelPosF   = 0.88             
4137                  amres@amOrthogonalPosF = 0.33           
4138                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4139                  overlay(plot(n),plot_us2)
4140                  us2=1
4141               else
4142                  if (prof3d .EQ. 0)then
4143                     varn=varn+1
4144                  end if
4145                  continue   
4146               end if
4147            end if
4148           
4149            if (vNam(varn) .EQ. "wsususodz" .OR. \
4150                vNam(varn) .EQ. "wspsodz" .OR.   \
4151                vNam(varn) .EQ. "wpeodz" .OR.    \
4152                vNam(varn) .EQ. "w*u*u*:dz" .OR. \
4153                vNam(varn) .EQ. "w*p*:dz" .OR.   \
4154                vNam(varn) .EQ. "w"+dq+"e:dz") then
4155               if (wsususodz .EQ. 0) then
4156                  res@gsnLeftString      = "w*u*u*:dz, w*p*:dz and w"+dq+"e:dz"
4157                  res@tiXAxisString      = "("+unit(varn)+")"
4158                  res@gsnRightString     = " "
4159                  if (xs .EQ. -1) then
4160                     res@trXMinF = min((/miniwsususodz,\
4161                                                   miniwspsodz,miniwpeodz/))
4162                  else
4163                     res@trXMinF = xs
4164                  end if
4165                  if (xe .EQ. -1) then
4166                     res@trXMaxF = max((/maxiwsususodz,maxiwspsodz,\
4167                                                               maxiwpeodz/))
4168                  else
4169                     res@trXMaxF = xe 
4170                  end if
4171                  if (vNam(varn) .EQ. "wspsodz")then
4172                     res@xyDashPattern = 1
4173                  end if
4174                  if (vNam(varn) .EQ. "wpeodz")then
4175                     res@xyDashPattern = 2
4176                  end if
4177                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4178
4179                  ; ***************************************************
4180                  ; legend for overlaid plot
4181                  ; ***************************************************
4182     
4183                  lgres                    = True
4184                  lgMonoDashIndex          = False
4185                  lgres@lgLabelFont        = "helvetica"   
4186                  lgres@lgLabelFontHeightF = font_size_legend*10.0           
4187                  lgres@vpWidthF           = 0.12           
4188                  lgres@vpHeightF          = 0.12         
4189                  lgres@lgDashIndexes      = (/0,1,2/)
4190                  lbid = gsn_create_legend(\
4191                           wks,3,(/"w*u*u*:dz","w*p*:dz","w"+dq+"e:dz"/),lgres)
4192                  amres = True
4193                  amres@amParallelPosF   = 0.88             
4194                  amres@amOrthogonalPosF = 0.33           
4195                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4196                  overlay(plot(n),plot_wsususodz)
4197                  wsususodz=1
4198               else
4199                  if (prof3d .EQ. 0)then
4200                     varn=varn+1
4201                  end if
4202                  continue   
4203               end if
4204            end if     
4205            n=n+1
4206            if (prof3d .EQ. 0)then
4207               varn=varn+1
4208            end if   
4209         end if
4210      end do
4211   end if
4212   
4213   com_var_avail=new(count_var,string)   
4214
4215   if (combine .EQ. 1) then
4216      co=0     
4217      n_o=0
4218      do varn = 0,dim-1
4219 
4220         check = True
4221     
4222         if (prof3d .EQ. 0) then
4223            if ( isStrSubset( vNam(varn), "time") .OR. \
4224                 isStrSubset( vNam(varn), "NORM")) then
4225               check = False
4226            end if
4227         else
4228            if ( isStrSubset( vNam(varn), "time") .OR.  \
4229                 isStrSubset( vNam(varn), "zusi") .OR.  \
4230                 isStrSubset( vNam(varn), "zwwi") .OR.  \
4231                 isStrSubset( vNam(varn), "x") .OR.     \
4232                 isStrSubset( vNam(varn), "xu") .OR.    \
4233                 isStrSubset( vNam(varn), "y") .OR.     \
4234                 isStrSubset( vNam(varn), "yv") .OR.    \
4235                 isStrSubset( vNam(varn), "zu_3d") .OR. \
4236                 isStrSubset( vNam(varn), "zw_3d")) then
4237               check = False
4238            end if
4239         end if
4240
4241         if (var .NE. "all") then         
4242            check = isStrSubset( var,","+vNam(varn)+"," )
4243         end if     
4244
4245         if (check)then
4246           
4247            if (prof3d .EQ. 0) then
4248               if (log_z .EQ. 1) then
4249                  z = f_att->$vNam(varn+1)$(1:dimz-1)
4250               else
4251                  z = f_att->$vNam(varn+1)$               
4252               end if
4253            else
4254               do i=0,b-1           
4255                  if (isStrSubset( a(i),"zu_3d" ))then
4256                     z_v(varn,:) = z_u
4257                     if (log_z .EQ. 1) then
4258                        z = z_v(varn,1:dimz-1)
4259                     else
4260                        z = z_v(varn,:)
4261                     end if
4262                  else
4263                     if (isStrSubset( a(i),"zw_3d" ))then
4264                        z_v(varn,:) = z_w
4265                        if (log_z .EQ. 1) then
4266                           z = z_v(varn,1:dimz-1)
4267                        else
4268                           z = z_v(varn,:)
4269                        end if
4270                     end if                   
4271                  end if
4272               end do           
4273            end if
4274
4275            z=z/norm_z
4276           
4277            com_var_avail(n_o)=vNam(varn)
4278           
4279            com=isStrSubset( c_var,","+vNam(varn)+"," )
4280            if (com)then
4281               co = co+1           
4282               if (n_o .EQ. 1) then
4283                  res@xyDashPattern  = 1                                   
4284               else           
4285                  if (n_o .EQ. 2) then
4286                     res@xyDashPattern  = 2
4287                  else
4288                     res@xyDashPattern  = 0
4289                     res@gsnLeftString  = "Combined Plot of "+c_var
4290                     res@tiXAxisString      = "("+unit(varn)+")"
4291                     res@gsnRightString     = " "
4292                     if (xs .EQ. -1) then
4293                        res@trXMinF = min(mini)
4294                     else
4295                        res@trXMinF = xs
4296                     end if
4297                     if (xe .EQ. -1) then
4298                        res@trXMaxF = max(maxi)
4299                     else
4300                        res@trXMaxF = xe 
4301                     end if
4302                  end if
4303               end if
4304               label(n_o)=vNam(varn)
4305               color_o(n_o)=237
4306               plot_o(n_o)=gsn_csm_xy(wks,data(varn,:,:),z,res)
4307               n_o=n_o+1
4308            end if
4309            if (prof3d .EQ. 0)then
4310               varn=varn+1
4311            end if           
4312         end if
4313      end do
4314   
4315      if(number_comb .EQ. 2)then
4316         if (co .EQ. 2)then
4317            overlay(plot_o(0),plot_o(1))
4318         else
4319            print(" ")
4320            print("combining is not possible,")
4321            print("'c_var'(= "+c_var+") must include two variables of "+\
4322                  "the general plots = ")
4323            print("- "+com_var_avail)
4324            print("be sure to have one comma before and after the variable")
4325            print(" ")
4326            exit 
4327         end if
4328      end if
4329      if(number_comb .EQ. 3)then
4330         if (co .EQ. 3)then
4331            overlay(plot_o(0),plot_o(1))
4332            overlay(plot_o(0),plot_o(2))
4333         else
4334            print(" ")
4335            print("combining is not possible,")
4336            print("'c_var'(= "+c_var+") must include three variables of "+\
4337                  "the general plots = ")
4338            print("- "+com_var_avail)
4339            print("be sure to have one comma before and after the variable")
4340            print(" ")
4341            exit
4342         end if
4343      end if
4344
4345      ; ***************************************************
4346      ; legend for combined plot
4347      ; ***************************************************
4348     
4349      lgres                    = True
4350      lgMonoDashIndex          = False
4351      lgres@lgDashIndexes      = (/0,1,2/)
4352      lgres@lgLabelFont        = "helvetica"   
4353      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4354      lgres@vpWidthF           = 0.12           
4355      lgres@vpHeightF          = 0.1           
4356 
4357      lbid = gsn_create_legend(wks,number_comb,label,lgres)       
4358
4359      amres = True
4360      amres@amParallelPosF   = 0.65                 
4361      amres@amOrthogonalPosF = -0.2           
4362      annoid1 = gsn_add_annotation(plot_o(0),lbid,amres)
4363   
4364      plot(0) = plot_o(0)
4365
4366   end if
4367
4368   ; ***************************************************
4369   ; merge plots onto one page
4370   ; ***************************************************
4371
4372   do m=0,n-1
4373      plot_(m)=plot(n-1-m)
4374   end do
4375
4376   no_frames = 0
4377
4378   if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
4379       n .gt. no_rows*no_columns) then
4380      gsn_panel(wks,plot_,(/n,1/),resP)
4381      print(" ")
4382      print("Outputs to .eps or .epsi have only one frame")
4383      print(" ")
4384   else   
4385      do i = 0,n-1, no_rows*no_columns
4386         if( (i+no_rows*no_columns) .gt. (n-1)) then
4387            gsn_panel(wks,plot_(i:n-1),(/no_rows,no_columns/),resP)
4388            no_frames = no_frames + 1 
4389         else
4390            gsn_panel(wks,plot_(i:i+no_rows*no_columns-1),\
4391                                                (/no_rows,no_columns/),resP)
4392            no_frames = no_frames + 1 
4393         end if
4394      end do
4395   end if
4396
4397    if (format_out .EQ. "png" ) then
4398     png_output = new((/no_frames/), string)
4399     j = 0
4400     do i=0, no_frames-1
4401       j = i + 1
4402       if (j .LE. 9) then
4403         png_output(i) = file_out+".00000"+j+".png"
4404       end if
4405       if (j .GT. 9 .AND. j .LE. 99) then
4406         png_output(i) = file_out+".0000"+j+".png"
4407       end if
4408       if (j .GT. 99 .AND. j .LE. 999) then
4409         png_output(i) = file_out+".000"+j+".png"
4410       end if
4411       if (j .GT. 999) then
4412         png_output(i) = file_out+".00"+j+".png"
4413       end if
4414
4415       ;using imagemagick's convert for reducing the white
4416       ;space around the plot
4417       cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
4418              png_output(i) + " " + png_output(i)
4419       system(cmd)
4420     end do
4421
4422     print(" ")
4423     print("Output to: "+ png_output)
4424     print(" ")
4425   else
4426     print(" ")
4427     print("Output to: " + file_out +"."+ format_out)
4428     print(" ")
4429   end if
4430
4431end
Note: See TracBrowser for help on using the repository browser.