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

Last change on this file since 2977 was 2837, checked in by gronemeier, 7 years ago

palmplot: config file of ncl scripts can be chosen in shell command

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