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

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

bugfixes in palmplot pr

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