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

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

bugfix: palmplot pr can now handle *pr.nc files consisting exclusively of user profiles

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