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

Last change on this file since 1204 was 1127, checked in by hoffmann, 12 years ago

bugfix in palmplot pr

File size: 142.3 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      if ( no_files .EQ. 1 )  then
1542         resP@gsnPanelXWhiteSpacePercent = 6.0
1543         resP@gsnPanelYWhiteSpacePercent = 6.0
1544      else
1545         resP@gsnPanelXWhiteSpacePercent = 6.0
1546         resP@gsnPanelYWhiteSpacePercent = 0.0
1547      end if
1548      resP@txFont                     = "helvetica"
1549      resP@txString                   = f_att@title
1550      resP@txFuncCode                 = "~"
1551      resP@txFontHeightF              = 0.0105
1552
1553      ; ***************************************************
1554      ; set up graphics for plot
1555      ; ***************************************************
1556
1557      if (combine .EQ. 1) then
1558          plot_o = new(number_comb,graphic) 
1559          label=new(number_comb,string)
1560          color_o=new(number_comb,integer)
1561
1562          if (check_vType) then
1563            mini=new(number_comb,double)
1564            maxi=new(number_comb,double)
1565          else
1566            mini=new(number_comb,float)
1567            maxi=new(number_comb,float)
1568          end if
1569      end if
1570   
1571      if ( format_out .EQ. "pdf" .OR. format_out .EQ. "ps" ) then
1572          format_out@wkPaperSize = "A4"
1573      end if
1574      if ( format_out .EQ. "png" ) then
1575          format_out@wkWidth  = 1000
1576          format_out@wkHeight = 1000
1577      end if
1578
1579      if (wks_gsn_log) then
1580          wks=gsn_open_wks(format_out,file_out)
1581          gsn_define_colormap(wks,"rainbow+white")
1582          wks_gsn_log = False
1583      end if
1584
1585      ; ***************************************************
1586      ; read data and create plots
1587      ; ***************************************************
1588         
1589      do ti =0,np-1
1590          if( t_all(ti_in(ti)) .lt. 10^36) then
1591            start_time_step = ti
1592            break
1593          end if
1594      end do 
1595     
1596      if (log_z .EQ. 1) then     
1597          if (check_vType) then
1598            data   = new((/dim,np,dimz-1/),double)
1599            data_0 = new((/np,dimz-1/),double)
1600          else
1601            data   = new((/dim,np,dimz-1/),float)
1602            data_0 = new((/np,dimz-1/),float)
1603          end if
1604          data@_FillValue=9.96921e+36
1605          data_0 = 0.1
1606          t      = new((/np,dimz-1/),float)
1607          t      = 0.0
1608          unit   = new(dim,string)
1609          if (isvar("z_u"))then
1610            if (typeof(z_u) .EQ. "double")then
1611                z_v    = new((/dim,dimz/),double)
1612                z_     = new((/dim,dimz-1/),double)
1613            else
1614                if (typeof(z_u) .EQ. "float")then
1615                  z_v    = new((/dim,dimz/),float)
1616                  z_     = new((/dim,dimz-1/),float)
1617                end if
1618            end if
1619          else
1620            if (isvar("z_w"))then
1621                if (typeof(z_w) .EQ. "double")then
1622                  z_v    = new((/dim,dimz/),double)
1623                  z_     = new((/dim,dimz-1/),double)
1624                else
1625                  if (typeof(z_w) .EQ. "float")then
1626                      z_v    = new((/dim,dimz/),float)
1627                      z_     = new((/dim,dimz-1/),float)
1628                  end if
1629                end if
1630            end if
1631          end if
1632      else
1633          if (check_vType) then
1634            data   = new((/dim,np,dimz/),double)
1635            data_0 = new((/np,dimz/),double)
1636          else
1637            data   = new((/dim,np,dimz/),float)
1638            data_0 = new((/np,dimz/),float)
1639          end if
1640          data@_FillValue=9.96921e+36     
1641          data_0 = 0.0
1642          t      = new((/np,dimz/),float)
1643          t      = 0.0
1644          unit   = new(dim,string)
1645          if (isvar("z_u"))then
1646            if (typeof(z_u) .EQ. "double")then
1647                z_v    = new((/dim,dimz/),double)
1648                z_     = new((/dim,dimz/),double)
1649            else
1650                if (typeof(z_u) .EQ. "float")then
1651                  z_v    = new((/dim,dimz/),float)
1652                  z_     = new((/dim,dimz/),float)
1653                end if
1654            end if
1655          else
1656            if (isvar("z_w"))then
1657                if (typeof(z_w) .EQ. "double")then
1658                  z_v    = new((/dim,dimz/),double)
1659                  z_     = new((/dim,dimz/),double)
1660                else
1661                  if (typeof(z_w) .EQ. "float")then
1662                      z_v    = new((/dim,dimz/),float)
1663                      z_     = new((/dim,dimz/),float)
1664                  end if
1665                end if
1666            end if
1667          end if
1668      end if
1669
1670      end if
1671      ;-------above steps only for first file
1672
1673      ; ***************************************************
1674      ; indicate plot number
1675      ; ***************************************************
1676     
1677      if (combine .EQ. 1) then
1678          n = 1
1679      else
1680          n = 0
1681      end if
1682
1683      if (over .EQ. 1) then
1684          plot_u         = gsn_csm_xy(wks,t,data_0(:,:),res)
1685          plot_v         = gsn_csm_xy(wks,t,data_0(:,:),res)
1686          plot_w         = gsn_csm_xy(wks,t,data_0(:,:),res)
1687          plot_pt        = gsn_csm_xy(wks,t,data_0(:,:),res)     
1688          plot_vpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
1689          plot_lpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
1690          plot_q         = gsn_csm_xy(wks,t,data_0(:,:),res)
1691          plot_qv        = gsn_csm_xy(wks,t,data_0(:,:),res)
1692          plot_ql        = gsn_csm_xy(wks,t,data_0(:,:),res)
1693          plot_rho       = gsn_csm_xy(wks,t,data_0(:,:),res)
1694          plot_s         = gsn_csm_xy(wks,t,data_0(:,:),res)
1695          plot_sa        = gsn_csm_xy(wks,t,data_0(:,:),res)
1696          plot_e         = gsn_csm_xy(wks,t,data_0(:,:),res)
1697          plot_es        = gsn_csm_xy(wks,t,data_0(:,:),res)
1698          plot_km        = gsn_csm_xy(wks,t,data_0(:,:),res)
1699          plot_kh        = gsn_csm_xy(wks,t,data_0(:,:),res)
1700          plot_l         = gsn_csm_xy(wks,t,data_0(:,:),res)     
1701          plot_wpup      = gsn_csm_xy(wks,t,data_0(:,:),res)
1702          plot_wsus      = gsn_csm_xy(wks,t,data_0(:,:),res)
1703          plot_wu        = gsn_csm_xy(wks,t,data_0(:,:),res)
1704          plot_wpvp      = gsn_csm_xy(wks,t,data_0(:,:),res)
1705          plot_wsvs      = gsn_csm_xy(wks,t,data_0(:,:),res)
1706          plot_wv        = gsn_csm_xy(wks,t,data_0(:,:),res)
1707          plot_wpptp     = gsn_csm_xy(wks,t,data_0(:,:),res)
1708          plot_wspts     = gsn_csm_xy(wks,t,data_0(:,:),res)
1709          plot_wpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
1710          plot_wsptsBC   = gsn_csm_xy(wks,t,data_0(:,:),res)
1711          plot_wptBC     = gsn_csm_xy(wks,t,data_0(:,:),res)
1712          plot_wpvptp    = gsn_csm_xy(wks,t,data_0(:,:),res)
1713          plot_wsvpts    = gsn_csm_xy(wks,t,data_0(:,:),res)
1714          plot_wvpt      = gsn_csm_xy(wks,t,data_0(:,:),res)
1715          plot_wpqp      = gsn_csm_xy(wks,t,data_0(:,:),res)
1716          plot_wsqs      = gsn_csm_xy(wks,t,data_0(:,:),res)
1717          plot_wq        = gsn_csm_xy(wks,t,data_0(:,:),res)
1718          plot_wpqvp     = gsn_csm_xy(wks,t,data_0(:,:),res)
1719          plot_wsqvs     = gsn_csm_xy(wks,t,data_0(:,:),res)
1720          plot_wqv       = gsn_csm_xy(wks,t,data_0(:,:),res)
1721          plot_wpsp      = gsn_csm_xy(wks,t,data_0(:,:),res)
1722          plot_wsss      = gsn_csm_xy(wks,t,data_0(:,:),res)
1723          plot_ws        = gsn_csm_xy(wks,t,data_0(:,:),res)
1724          plot_wpsap     = gsn_csm_xy(wks,t,data_0(:,:),res)
1725          plot_wssas     = gsn_csm_xy(wks,t,data_0(:,:),res)
1726          plot_wsa       = gsn_csm_xy(wks,t,data_0(:,:),res)
1727          plot_wses      = gsn_csm_xy(wks,t,data_0(:,:),res)
1728          plot_us2       = gsn_csm_xy(wks,t,data_0(:,:),res)
1729          plot_vs2       = gsn_csm_xy(wks,t,data_0(:,:),res)
1730          plot_ws2       = gsn_csm_xy(wks,t,data_0(:,:),res)
1731          plot_pts2      = gsn_csm_xy(wks,t,data_0(:,:),res)
1732          plot_ws3       = gsn_csm_xy(wks,t,data_0(:,:),res)
1733          plot_Sw        = gsn_csm_xy(wks,t,data_0(:,:),res)
1734          plot_ws2pts    = gsn_csm_xy(wks,t,data_0(:,:),res)
1735          plot_wspts2    = gsn_csm_xy(wks,t,data_0(:,:),res)
1736          plot_wsususodz = gsn_csm_xy(wks,t,data_0(:,:),res)
1737          plot_wspsodz   = gsn_csm_xy(wks,t,data_0(:,:),res)
1738          plot_wpeodz    = gsn_csm_xy(wks,t,data_0(:,:),res)
1739     
1740          if (check_vType) then
1741            miniu         =  100000.d
1742            maxiu         = -100000.d
1743            miniv         =  100000.d
1744            maxiv         = -100000.d
1745            miniw         =  100000.d
1746            maxiw         = -100000.d
1747            minipt        =  100000.d
1748            maxipt        = -100000.d
1749            minivpt       =  100000.d
1750            maxivpt       = -100000.d
1751            minilpt       =  100000.d
1752            maxilpt       = -100000.d
1753            miniq         =  100000.d
1754            maxiq         = -100000.d
1755            miniqv        =  100000.d
1756            maxiqv        = -100000.d
1757            miniql        =  100000.d
1758            maxiql        = -100000.d
1759            minie         =  100000.d
1760            maxie         = -100000.d
1761            minies        =  100000.d
1762            maxies        = -100000.d
1763            minikm        =  100000.d
1764            maxikm        = -100000.d
1765            minikh        =  100000.d
1766            maxikh        = -100000.d
1767            miniwpup      =  100000.d
1768            maxiwpup      = -100000.d
1769            miniwsus      =  100000.d
1770            maxiwsus      = -100000.d
1771            miniwu        =  100000.d
1772            maxiwu        = -100000.d
1773            miniwpvp      =  100000.d
1774            maxiwpvp      = -100000.d
1775            miniwsvs      =  100000.d
1776            maxiwsvs      = -100000.d
1777            miniwv        =  100000.d
1778            maxiwv        = -100000.d
1779            miniwpptp     =  100000.d
1780            maxiwpptp     = -100000.d
1781            miniwspts     =  100000.d
1782            maxiwspts     = -100000.d
1783            miniwpt       =  100000.d
1784            maxiwpt       = -100000.d
1785            miniwsptsBC   =  100000.d
1786            maxiwsptsBC   = -100000.d
1787            miniwptBC     =  100000.d
1788            maxiwptBC     = -100000.d
1789            miniwpvptp    =  100000.d
1790            maxiwpvptp    = -100000.d
1791            miniwsvpts    =  100000.d
1792            maxiwsvpts    = -100000.d
1793            miniwvpt      =  100000.d
1794            maxiwvpt      = -100000.d
1795            miniwpqp      =  100000.d
1796            maxiwpqp      = -100000.d
1797            miniwsqs      =  100000.d
1798            maxiwsqs      = -100000.d
1799            miniwq        =  100000.d
1800            maxiwq        = -100000.d
1801            miniwpqvp     =  100000.d
1802            maxiwpqvp     = -100000.d
1803            miniwsqvs     =  100000.d
1804            maxiwsqvs     = -100000.d
1805            miniwqv       =  100000.d
1806            maxiwqv       = -100000.d
1807            miniwpsp      =  100000.d
1808            maxiwpsp      = -100000.d
1809            miniwsss      =  100000.d
1810            maxiwsss      = -100000.d
1811            miniws        =  100000.d
1812            maxiws        = -100000.d
1813            miniwpsap     =  100000.d
1814            maxiwpsap     = -100000.d
1815            miniwssas     =  100000.d
1816            maxiwssas     = -100000.d
1817            miniwsa       =  100000.d
1818            maxiwsa       = -100000.d
1819            minius2       =  100000.d
1820            maxius2       = -100000.d
1821            minivs2       =  100000.d
1822            maxivs2       = -100000.d
1823            miniws2       =  100000.d
1824            maxiws2       = -100000.d
1825            miniwsususodz =  100000.d
1826            maxiwsususodz = -100000.d
1827            miniwspsodz   =  100000.d
1828            maxiwspsodz   = -100000.d
1829            miniwpeodz    =  100000.d
1830            maxiwpeodz    = -100000.d
1831          else
1832            miniu         =  100000.
1833            maxiu         = -100000.
1834            miniv         =  100000.
1835            maxiv         = -100000.
1836            miniw         =  100000.
1837            maxiw         = -100000.
1838            minipt        =  100000.
1839            maxipt        = -100000.
1840            minivpt       =  100000.
1841            maxivpt       = -100000.
1842            minilpt       =  100000.
1843            maxilpt       = -100000.
1844            miniq         =  100000.
1845            maxiq         = -100000.
1846            miniqv        =  100000.
1847            maxiqv        = -100000.
1848            miniql        =  100000.
1849            maxiql        = -100000.
1850            minie         =  100000.
1851            maxie         = -100000.
1852            minies        =  100000.
1853            maxies        = -100000.
1854            minikm        =  100000.
1855            maxikm        = -100000.
1856            minikh        =  100000.
1857            maxikh        = -100000.
1858            miniwpup      =  100000.
1859            maxiwpup      = -100000.
1860            miniwsus      =  100000.
1861            maxiwsus      = -100000.
1862            miniwu        =  100000.
1863            maxiwu        = -100000.
1864            miniwpvp      =  100000.
1865            maxiwpvp      = -100000.
1866            miniwsvs      =  100000.
1867            maxiwsvs      = -100000.
1868            miniwv        =  100000.
1869            maxiwv        = -100000.
1870            miniwpptp     =  100000.
1871            maxiwpptp     = -100000.
1872            miniwspts     =  100000.
1873            maxiwspts     = -100000.
1874            miniwpt       =  100000.
1875            maxiwpt       = -100000.
1876            miniwsptsBC   =  100000.
1877            maxiwsptsBC   = -100000.
1878            miniwptBC     =  100000.
1879            maxiwptBC     = -100000.
1880            miniwpvptp    =  100000.
1881            maxiwpvptp    = -100000.
1882            miniwsvpts    =  100000.
1883            maxiwsvpts    = -100000.
1884            miniwvpt      =  100000.
1885            maxiwvpt      = -100000.
1886            miniwpqp      =  100000.
1887            maxiwpqp      = -100000.
1888            miniwsqs      =  100000.
1889            maxiwsqs      = -100000.
1890            miniwq        =  100000.
1891            maxiwq        = -100000.
1892            miniwpqvp     =  100000.
1893            maxiwpqvp     = -100000.
1894            miniwsqvs     =  100000.
1895            maxiwsqvs     = -100000.
1896            miniwqv       =  100000.
1897            maxiwqv       = -100000.
1898            miniwpsp      =  100000.
1899            maxiwpsp      = -100000.
1900            miniwsss      =  100000.
1901            maxiwsss      = -100000.
1902            miniws        =  100000.
1903            maxiws        = -100000.
1904            miniwpsap     =  100000.
1905            maxiwpsap     = -100000.
1906            miniwssas     =  100000.
1907            maxiwssas     = -100000.
1908            miniwsa       =  100000.
1909            maxiwsa       = -100000.
1910            minius2       =  100000.
1911            maxius2       = -100000.
1912            minivs2       =  100000.
1913            maxivs2       = -100000.
1914            miniws2       =  100000.
1915            maxiws2       = -100000.
1916            miniwsususodz =  100000.
1917            maxiwsususodz = -100000.
1918            miniwspsodz   =  100000.
1919            maxiwspsodz   = -100000.
1920            miniwpeodz    =  100000.
1921            maxiwpeodz    = -100000.
1922          end if
1923
1924      end if
1925
1926      if (prof3d .EQ. 1)then
1927
1928      if (end_x .EQ. -1) then
1929          end_x=dimx-2
1930      end if
1931      if (end_y .EQ. -1)then
1932          end_y=dimy-2
1933      end if
1934      if (start_x .LT. 0)then
1935          print(" ")
1936          print("'start_x' is lower than 0 and set to 0")
1937          print(" ")
1938          start_x=0
1939      end if
1940      if (start_x .GT. dimx-1)then
1941          print(" ")
1942          print("'start_x' is greater than available x-range and set to "+\
1943                "maximum of x-range (excluding ghostpoint)")
1944          print(" ")
1945          start_x=dimx-2
1946      end if
1947      if (end_x .EQ. dimx-1)then
1948          print(" ")
1949          print("'end_x' = "+end_x+" and includes the ghostpoint")
1950          print(" ")
1951      end if
1952      if (end_x .GT. dimx-1)then
1953          print(" ")
1954          print("'end_x' = "+end_x+" is greater than available x-range and set "+\
1955                "to maximum of x-range (excluding ghostpoint)")
1956          print(" ")
1957          end_x=dimx-2
1958      end if
1959      if (end_x .LT. start_x)then
1960          print(" ")
1961          print("'end_x' = "+end_x+" is lower than 'start_x' = "+start_x+\
1962                " and set to maximum of x-range (excluding ghostpoint)")
1963          print(" ")
1964          end_x=dimx-2
1965      end if
1966      if (start_y .LT. 0)then
1967          print(" ")
1968          print("'start_y' is lower than 0 and set to 0")
1969          print(" ")
1970          start_y=0
1971      end if
1972      if (start_y .GT. dimy-1)then
1973          print(" ")
1974          print("'start_y' is greater than available y-range and set to "+\
1975                "maximum of y-range (excluding ghostpoint)")
1976          print(" ")
1977          start_x=dimy-2
1978      end if
1979      if (end_y .EQ. dimy-1)then
1980          print(" ")
1981          print("'end_y' = "+end_y+" and includes the ghostpoint")
1982          print(" ")
1983      end if
1984      if (end_y .GT. dimy-1)then
1985          print(" ")
1986          print("'end_y' = "+end_y+" is greater than available y-range and "+\
1987                "set to maximum of y-range (excluding ghostpoint)")
1988          print(" ")
1989          end_x=dimy-2
1990      end if
1991      if (end_y .LT. start_y)then
1992          print(" ")
1993          print("'end_y' = "+end_y+" is lower than 'start_y' = "+start_y+\
1994                " and set to maximum of y-range (excluding ghostpoint)")
1995          print(" ")
1996          end_y=dimy-2
1997      end if
1998     
1999      end if
2000     
2001      n_o=0
2002      count_var=0
2003
2004      res@xyDashPattern = 1*nof
2005
2006      over_remind = False
2007      if ( over .eq. 1)then
2008        over_remind = True
2009      end if
2010
2011      do varn = 0,dim-1
2012
2013          check = True
2014         
2015          if (prof3d .EQ. 0) then
2016            if ( isStrSubset( vNam(varn), "time") .OR. \
2017                  isStrSubset( vNam(varn), "NORM")) then
2018                check = False
2019            end if
2020          else
2021            if ( isStrSubset( vNam(varn), "time") .OR.  \
2022                  isStrSubset( vNam(varn), "zusi") .OR.  \
2023                  isStrSubset( vNam(varn), "zwwi") .OR.  \
2024                  isStrSubset( vNam(varn), "x") .OR.     \
2025                  isStrSubset( vNam(varn), "xu") .OR.    \
2026                  isStrSubset( vNam(varn), "y") .OR.     \
2027                  isStrSubset( vNam(varn), "yv") .OR.    \
2028                  isStrSubset( vNam(varn), "zu_3d") .OR. \
2029                  isStrSubset( vNam(varn), "zw_3d")) then
2030                check = False
2031            end if
2032          end if
2033
2034          if (var .NE. "all") then
2035            check = isStrSubset( var,","+vNam(varn)+"," )
2036          end if
2037
2038          if (combine .EQ. 1) then         
2039            com=isStrSubset(c_var,","+vNam(varn)+"," )     
2040            if (com) then     
2041                if (prof3d .EQ. 0) then             
2042                  temp     = f[:]->$vNam(varn)$       
2043                  temp_att = f_att->$vNam(varn)$
2044                  if (log_z .EQ. 1) then
2045                      do j=0,np-1
2046                        data(varn,j,:) = temp(ti_in(j),1:dimz-1)
2047                      end do
2048                  else
2049                      do j=0,np-1
2050                        data(varn,j,:) = temp(ti_in(j),0:dimz-1)
2051                      end do
2052                  end if
2053                else
2054                  if (log_z .EQ. 1) then
2055                      do i=1,dimz-1
2056                        do j=0,np-1
2057                            temp = f[:]->$vNam(varn)$
2058                            temp_att = f_att->$vNam(varn)$
2059                            data_temp = temp(ti_in(j),i,\
2060                                            start_y:end_y,start_x:end_x)
2061                            data(varn,j,i-1) = dim_avg_Wrap(\
2062                                                      dim_avg_Wrap(data_temp))
2063                        end do
2064                      end do
2065                  else
2066                      do i=0,dimz-1
2067                        do j=0,np-1
2068                            temp = f[:]->$vNam(varn)$
2069                            temp_att = f_att->$vNam(varn)$
2070                            data_temp = temp(ti_in(j),i,\
2071                                                start_y:end_y,start_x:end_x)
2072                            data(varn,j,i) = dim_avg_Wrap(dim_avg_Wrap(data_temp))
2073                        end do
2074                      end do
2075                  end if
2076                  print(" ")
2077                  print("Variable for combine '"+vNam(varn)+"' is read")
2078                  print(" ")
2079                end if                 
2080                unit(varn) = temp_att@units
2081                if (n_o .GT. number_comb-1) then
2082                  print(" ")
2083                  print("Set 'number_comb' to the number of overlaying "+\
2084                        "variables ('c_var' = "+c_var+")")
2085                  print(" ")
2086                  exit
2087                end if
2088                if (abs(min(data(varn,:,:))) .GT. 10)then
2089                  min_value = abs(0.01*min(data(varn,:,:)))
2090                  max_value = abs(0.01*max(data(varn,:,:)))
2091                else
2092                  if (abs(min(data(varn,:,:))) .LT. 0.01 .AND. \
2093                      abs(max(data(varn,:,:))) .GT. 0.01)then
2094                      min_value = abs(0.1*max(data(varn,:,:)))
2095                      max_value = abs(0.1*max(data(varn,:,:)))
2096                  else
2097                      if (abs(max(data(varn,:,:))) .LT. 0.01 .AND. \
2098                          abs(min(data(varn,:,:))) .GT. 0.01)then
2099                        min_value = abs(0.1*min(data(varn,:,:)))
2100                        max_value = abs(0.1*min(data(varn,:,:)))
2101                      else
2102                        min_value = abs(0.1*min(data(varn,:,:)))
2103                        max_value = abs(0.1*max(data(varn,:,:)))
2104                      end if
2105                  end if
2106                end if
2107                if (min(data(varn,:,:)) .EQ. 0 .AND. \
2108                    max(data(varn,:,:)) .EQ. 0)then
2109                  min_value = 0.1
2110                  max_value = 0.1
2111                end if
2112                mini(n_o)=min(data(varn,:,:))
2113                maxi(n_o)=max(data(varn,:,:))
2114                n_o=n_o+1
2115            end if
2116          end if
2117
2118          if(check) then
2119           
2120            count_var=count_var+1
2121
2122            if (prof3d .EQ. 0) then       
2123                temp = f[:]->$vNam(varn)$
2124                temp_att = f_att->$vNam(varn)$
2125            else
2126                if (log_z .EQ. 1) then
2127                  do i=1,dimz-1
2128                      do j=0,np-1
2129                        temp= f[:]->$vNam(varn)$
2130                        temp_att = f_att->$vNam(varn)$
2131                        data_temp = temp(ti_in(j),i,start_y:end_y,start_x:end_x)
2132                        data(varn,j,i-1) = dim_avg_Wrap(dim_avg_Wrap(data_temp))
2133                        delete(data_temp)
2134                      end do
2135                  end do
2136                else
2137                  do i=0,dimz-1
2138                      do j=0,np-1
2139                        temp= f[:]->$vNam(varn)$
2140                        temp_att = f_att->$vNam(varn)$
2141                        data_temp = temp(ti_in(j),i,start_y:end_y,start_x:end_x)
2142                        data_temp!0 = "t"
2143                        data_temp!1 = "z"
2144                        data(varn,j,i) = dim_avg_Wrap(dim_avg_Wrap(data_temp))
2145                        delete(data_temp)
2146                      end do
2147                  end do
2148                end if
2149                print(" ")
2150                print("Variable '"+vNam(varn)+"' is read")
2151                print(" ")
2152                unit(varn) = temp_att@units
2153                a=getvardims(temp_att)
2154                b=dimsizes(a)
2155            end if
2156
2157            if (prof3d .EQ. 0) then
2158                if (log_z .EQ. 1) then
2159                  z = f_att->$vNam(varn+1)$(1:dimz-1)
2160                  unit(varn) = temp_att@units
2161                  do j=0,np-1
2162                      data(varn,j,:) = temp(ti_in(j),1:dimz-1)
2163                  end do
2164                else
2165                  z = f_att->$vNam(varn+1)$
2166                  unit(varn) = temp_att@units
2167                  do j=0,np-1
2168                      data(varn,j,:) = temp(ti_in(j),:)
2169                  end do
2170                end if
2171            else
2172                do i=0,b-1           
2173                  if (isStrSubset( a(i),"zu_3d" ))then
2174                      z_v(varn,:) = z_u
2175                      if (log_z .EQ. 1) then
2176                        z = z_v(varn,1:dimz-1)
2177                      else
2178                        z = z_v(varn,:)
2179                      end if
2180                  else
2181                      if (isStrSubset( a(i),"zw_3d" ))then
2182                        z_v(varn,:) = z_w
2183                        if (log_z .EQ. 1) then
2184                            z = z_v(varn,1:dimz-1)
2185                        else
2186                            z = z_v(varn,:)
2187                        end if
2188                      end if                   
2189                  end if
2190                end do           
2191            end if
2192
2193            if (isStrSubset(data@long_name," SR " ) ) then
2194              over = 0
2195            end if
2196           
2197            if (nof .EQ. 0) then
2198                z_(n,:)=z/norm_z
2199                z    = z_(n,:)
2200            else
2201                z=z/norm_z
2202            end if
2203
2204            delta_z = z(2) - z(1)
2205
2206            max_z_int=doubletoint(max_z/delta_z)
2207            min_z_int=doubletoint(min_z/delta_z)
2208
2209            if(max_z_int .eq. min_z_int)
2210                print(" ")
2211                print("Please increase 'max_z' or decrease 'min_z' so that "+\
2212                      "there are")
2213                print("at least two layers for the z-axis to plot")
2214                print(" ")
2215                exit
2216            end if
2217
2218            if(min_z_int .gt. max_z_int)
2219                print(" ")
2220                print("'min_z' is greater than 'max_z',")
2221                print("please change this")
2222                print(" ")
2223                exit
2224            end if
2225
2226            if (max_z_int .ge. dimz-1)
2227                max_z_int = dimz-1
2228                if (log_z .EQ. 1) then
2229                max_z_int = max_z_int - 1
2230                end if
2231            end if
2232
2233            if (min_z_int .lt. 0)
2234                min_z_int = 0
2235                if (log_z .EQ. 1) then
2236                  min_z_int = min_z_int + 1
2237                end if
2238            end if         
2239
2240            ;data can contain missing values in case of output of t=0h (inital profiles)
2241            ;where no output is possible
2242            n_not_ismissing = num(.not.ismissing(data(varn,:,:)))
2243
2244            if (n_not_ismissing .GT. 0 ) then   
2245   
2246              if (abs(min(data(varn,:,min_z_int:max_z_int))) .GT. 10)then
2247                  min_value = abs(0.001*min(data(varn,:,min_z_int:max_z_int)))
2248                  max_value = abs(0.001*max(data(varn,:,min_z_int:max_z_int)))
2249              else
2250                  if (abs(min(data(varn,:,min_z_int:max_z_int))) .LT. 0.01 .AND. \
2251                      abs(max(data(varn,:,min_z_int:max_z_int))) .GT. 0.01)then
2252                    min_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
2253                    max_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
2254                  else
2255                    if (abs(max(data(varn,:,min_z_int:max_z_int))) .LT. 0.01 .AND.\
2256                        abs(min(data(varn,:,min_z_int:max_z_int))) .GT. 0.01)then
2257                        min_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
2258                        max_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
2259                    else
2260                        min_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
2261                        max_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
2262                    end if
2263                  end if
2264              end if
2265           
2266              if (min(data(varn,:,min_z_int:max_z_int)) .EQ. 0 .AND. \
2267                  max(data(varn,:,min_z_int:max_z_int)) .EQ. 0)then
2268                  min_value = 0.1
2269                  max_value = 0.1
2270              end if
2271
2272            else     
2273              print(" ")
2274              print(vNam(varn) + " contains only missing values")
2275              print(" ")
2276            end if   
2277
2278            if (over .EQ. 0) then 
2279                res@gsnLeftString      = " "
2280                res@tiXAxisString      = "("+unit(varn)+")"
2281                res@gsnRightString     = " "
2282                res@trYMinF            = min_z
2283                res@trYMaxF            = max_z
2284                if (xs .EQ. -1) then
2285                  res@trXMinF            = min(data(varn,:,min_z_int:max_z_int))-\
2286                                                                          min_value
2287                else
2288                  res@trXMinF            = xs     
2289                end if
2290                if (xe .EQ. -1) then
2291                  res@trXMaxF            = max(data(varn,:,min_z_int:max_z_int))+\
2292                                                                          max_value
2293                else
2294                  res@trXMaxF            = xe 
2295                end if         
2296
2297                if (c_var_log .EQ. 1) then
2298                  if ((mod(com_i, no_rows * no_columns) .EQ. 0) .AND. (no_files .LT. 2))then
2299                    res@pmLegendDisplayMode  = "Always"
2300                  else
2301                    res@pmLegendDisplayMode  = "Never"
2302                  end if
2303                else
2304                  if ((mod(no_plots-1-n+combine, no_rows * no_columns) .EQ. 0) .AND. (no_files .LT. 2))then
2305                    res@pmLegendDisplayMode  = "Always"
2306                  else
2307                    res@pmLegendDisplayMode  = "Never"
2308                  end if
2309                end if
2310
2311                plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res) 
2312
2313            end if
2314     
2315       
2316            if (vNam(varn) .EQ. "u" ) then
2317                  miniu=min(data(varn,:,min_z_int:max_z_int))-min_value
2318                  maxiu=max(data(varn,:,min_z_int:max_z_int))+max_value
2319                  if (over .EQ. 1) then
2320                      res@xyDashPattern  = 0
2321                      plot_u = gsn_csm_xy(wks,data(varn,:,:),z,res)
2322                  else
2323                      res@gsnLeftString      = " "
2324                      res@tiXAxisString      = "("+unit(varn)+")"
2325                      res@gsnRightString     = " "                 
2326                      if (xs .EQ. -1) then 
2327                        res@trXMinF            = miniu
2328                      else
2329                        res@trXMinF            = xs
2330                      end if
2331                      if (xe .EQ. -1) then       
2332                        res@trXMaxF            = maxiu
2333                      else
2334                        res@trXMaxF            = xe 
2335                      end if               
2336                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res) 
2337                  end if
2338                end if
2339                if (vNam(varn) .EQ. "v") then
2340                  miniv=min(data(varn,:,min_z_int:max_z_int))-min_value
2341                  maxiv=max(data(varn,:,min_z_int:max_z_int))+max_value
2342                  if (over .EQ. 1) then
2343                      res@xyMonoDashPattern = True
2344                      res@xyDashPattern  = 1
2345                      plot_v = gsn_csm_xy(wks,data(varn,:,:),z,res)
2346                  else
2347                      res@gsnLeftString      = " "
2348                      res@tiXAxisString      = "("+unit(varn)+")"
2349                      res@gsnRightString     = " "                 
2350                      if (xs .EQ. -1) then
2351                        res@trXMinF            = miniv
2352                      else
2353                        res@trXMinF            = xs
2354                      end if
2355                      if (xe .EQ. -1) then
2356                        res@trXMaxF            = maxiv
2357                      else
2358                        res@trXMaxF            = xe 
2359                      end if               
2360                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2361                  end if 
2362                end if
2363                if (vNam(varn) .EQ. "w") then
2364                  miniw=min(data(varn,:,min_z_int:max_z_int))-min_value
2365                  maxiw=max(data(varn,:,min_z_int:max_z_int))+max_value
2366                  if (over .EQ. 1) then
2367                      res@xyDashPattern = 2
2368                      plot_w = gsn_csm_xy(wks,data(varn,:,:),z,res)
2369                  else
2370                      res@gsnLeftString      = " "
2371                      res@tiXAxisString      = "("+unit(varn)+")"
2372                      res@gsnRightString     = " "
2373                      if (xs .EQ. -1) then
2374                        res@trXMinF            = miniw
2375                      else
2376                        res@trXMinF            = xs
2377                      end if
2378                      if (xe .EQ. -1) then
2379                        res@trXMaxF            = maxiw
2380                      else
2381                        res@trXMaxF            = xe 
2382                      end if           
2383                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2384                  end if   
2385                end if
2386
2387                if (vNam(varn) .EQ. "pt") then
2388                  minipt=min(data(varn,:,min_z_int:max_z_int))-min_value
2389                  maxipt=max(data(varn,:,min_z_int:max_z_int))+max_value
2390                  if (over .EQ. 1) then
2391                      res@xyDashPattern  = 0
2392                      plot_pt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2393                  else
2394                      res@gsnLeftString      = " "
2395                      res@tiXAxisString      = "("+unit(varn)+")"
2396                      res@gsnRightString     = " "
2397                      if (xs .EQ. -1) then
2398                        res@trXMinF            = minipt
2399                      else
2400                        res@trXMinF            = xs
2401                      end if
2402                      if (xe .EQ. -1) then       
2403                        res@trXMaxF            = maxipt
2404                      else
2405                        res@trXMaxF            = xe 
2406                      end if
2407                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2408                  end if 
2409                end if
2410                if (vNam(varn) .EQ. "vpt") then
2411                  minivpt=min(data(varn,:,min_z_int:max_z_int))-min_value
2412                  maxivpt=max(data(varn,:,min_z_int:max_z_int))+max_value
2413                  if (over .EQ. 1) then
2414                      res@xyDashPattern  = 1
2415                      plot_vpt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2416                  else
2417                      res@gsnLeftString      = " "
2418                      res@tiXAxisString      = "("+unit(varn)+")"
2419                      res@gsnRightString     = " "
2420                      if (xs .EQ. -1) then
2421                        res@trXMinF            = minivpt
2422                      else
2423                        res@trXMinF            = xs
2424                      end if
2425                      if (xe .EQ. -1) then       
2426                        res@trXMaxF            = maxivpt
2427                      else
2428                        res@trXMaxF            = xe 
2429                      end if
2430                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2431                  end if
2432                end if
2433                if (vNam(varn) .EQ. "lpt") then
2434                  minilpt=min(data(varn,:,min_z_int:max_z_int))-min_value
2435                  maxilpt=max(data(varn,:,min_z_int:max_z_int))+max_value
2436                  if (over .EQ. 1) then
2437                      res@xyDashPattern  = 2
2438                      plot_lpt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2439                  else
2440                      res@gsnLeftString      = " "
2441                      res@tiXAxisString      = "("+unit(varn)+")"
2442                      res@gsnRightString     = " "
2443                      if (xs .EQ. -1) then
2444                        res@trXMinF            = minilpt
2445                      else
2446                        res@trXMinF            = xs
2447                      end if
2448                      if (xe .EQ. -1) then       
2449                        res@trXMaxF            = maxilpt
2450                      else
2451                        res@trXMaxF            = xe 
2452                      end if
2453                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2454                  end if
2455                end if
2456
2457                if (vNam(varn) .EQ. "q") then
2458                  miniq=min(data(varn,:,min_z_int:max_z_int))-min_value
2459                  maxiq=max(data(varn,:,min_z_int:max_z_int))+max_value
2460                  if (over .EQ. 1) then
2461                      res@xyDashPattern  = 0
2462                      plot_q = gsn_csm_xy(wks,data(varn,:,:),z,res)
2463                  else
2464                      res@gsnLeftString      = " "
2465                      res@tiXAxisString      = "("+unit(varn)+")"
2466                      res@gsnRightString     = " "
2467                      if (xs .EQ. -1) then
2468                        res@trXMinF            = miniq
2469                      else
2470                        res@trXMinF            = xs
2471                      end if
2472                      if (xe .EQ. -1) then       
2473                        res@trXMaxF            = maxiq
2474                      else
2475                        res@trXMaxF            = xe 
2476                      end if
2477                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2478                  end if
2479                end if
2480                if (vNam(varn) .EQ. "qv") then
2481                  miniqv=min(data(varn,:,min_z_int:max_z_int))-min_value
2482                  maxiqv=max(data(varn,:,min_z_int:max_z_int))+max_value
2483                  if (over .EQ. 1) then
2484                      res@xyDashPattern  = 1
2485                      plot_qv = gsn_csm_xy(wks,data(varn,:,:),z,res)
2486                  else
2487                      res@gsnLeftString      = " "
2488                      res@tiXAxisString      = "("+unit(varn)+")"
2489                      res@gsnRightString     = " "
2490                      if (xs .EQ. -1) then
2491                        res@trXMinF            = miniqv
2492                      else
2493                        res@trXMinF            = xs
2494                      end if
2495                      if (xe .EQ. -1) then         
2496                        res@trXMaxF            = maxiqv
2497                      else
2498                        res@trXMaxF            = xe 
2499                      end if
2500                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2501                  end if
2502                end if
2503                if (vNam(varn) .EQ. "ql") then
2504                  miniql=min(data(varn,:,min_z_int:max_z_int))-min_value
2505                  maxiql=max(data(varn,:,min_z_int:max_z_int))+max_value
2506                  if (over .EQ. 1) then
2507                      res@xyDashPattern  = 2
2508                      plot_ql = gsn_csm_xy(wks,data(varn,:,:),z,res)
2509                  else
2510                      res@gsnLeftString      = " "
2511                      res@tiXAxisString      = "("+unit(varn)+")"
2512                      res@gsnRightString     = " "
2513                      if (xs .EQ. -1) then
2514                        res@trXMinF            = miniql
2515                      else
2516                        res@trXMinF            = xs
2517                      end if
2518                      if (xe .EQ. -1) then       
2519                        res@trXMaxF            = maxiql
2520                      else
2521                        res@trXMaxF            = xe 
2522                      end if
2523                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2524                  end if
2525                end if
2526
2527                if (vNam(varn) .EQ. "e") then
2528                  minie=min(data(varn,:,min_z_int:max_z_int))-min_value
2529                  maxie=max(data(varn,:,min_z_int:max_z_int))+max_value
2530                  if (over .EQ. 1) then
2531                      res@xyDashPattern  = 0
2532                      plot_e = gsn_csm_xy(wks,data(varn,:,:),z,res)
2533                  else
2534                      res@gsnLeftString      = " "
2535                      res@tiXAxisString      = "("+unit(varn)+")"
2536                      res@gsnRightString     = " "
2537                      if (xs .EQ. -1) then
2538                        res@trXMinF            = minie
2539                      else
2540                        res@trXMinF            = xs
2541                      end if
2542                      if (xe .EQ. -1) then       
2543                        res@trXMaxF            = maxie
2544                      else
2545                        res@trXMaxF            = xe 
2546                      end if
2547                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2548                  end if
2549                end if
2550                if (vNam(varn) .EQ. "es" .OR. vNam(varn) .EQ. "e*") then
2551                  minies=min(data(varn,:,min_z_int:max_z_int))-min_value
2552                  maxies=max(data(varn,:,min_z_int:max_z_int))+max_value
2553                  if (over .EQ. 1) then
2554                      res@xyDashPattern  = 1
2555                      plot_es = gsn_csm_xy(wks,data(varn,:,:),z,res)
2556                  else
2557                      res@gsnLeftString      = " "
2558                      res@tiXAxisString      = "("+unit(varn)+")"
2559                      res@gsnRightString     = " "
2560                      if (xs .EQ. -1) then
2561                        res@trXMinF            = minies
2562                      else
2563                        res@trXMinF            = xs
2564                      end if
2565                      if (xe .EQ. -1) then       
2566                        res@trXMaxF            = maxies
2567                      else
2568                        res@trXMaxF            = xe 
2569                      end if
2570                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2571                  end if
2572                end if
2573
2574                if (vNam(varn) .EQ. "km") then
2575                  minikm=min(data(varn,:,min_z_int:max_z_int))-min_value
2576                  maxikm=max(data(varn,:,min_z_int:max_z_int))+max_value
2577                  if (over .EQ. 1) then
2578                      res@xyDashPattern  = 0
2579                      plot_km = gsn_csm_xy(wks,data(varn,:,:),z,res)
2580                  else
2581                      res@gsnLeftString      = " "
2582                      res@tiXAxisString      = "("+unit(varn)+")"
2583                      res@gsnRightString     = " "
2584                      if (xs .EQ. -1) then
2585                        res@trXMinF            = minikm
2586                      else
2587                        res@trXMinF            = xs
2588                      end if
2589                      if (xe .EQ. -1) then     
2590                        res@trXMaxF            = maxikm
2591                      else
2592                        res@trXMaxF            = xe 
2593                      end if
2594                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2595                  end if
2596                end if
2597                if (vNam(varn) .EQ. "kh") then
2598                  minikh=min(data(varn,:,min_z_int:max_z_int))-min_value
2599                  maxikh=max(data(varn,:,min_z_int:max_z_int))+max_value
2600                  if (over .EQ. 1) then
2601                      res@xyDashPattern  = 1
2602                      plot_kh = gsn_csm_xy(wks,data(varn,:,:),z,res)
2603                  else
2604                      res@gsnLeftString      = " "
2605                      res@tiXAxisString      = "("+unit(varn)+")"
2606                      res@gsnRightString     = " "
2607                      if (xs .EQ. -1) then
2608                        res@trXMinF            = minikh
2609                      else
2610                        res@trXMinF            = xs
2611                      end if
2612                      if (xe .EQ. -1) then     
2613                        res@trXMaxF            = maxikh
2614                      else
2615                        res@trXMaxF            = xe 
2616                      end if
2617                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2618                  end if
2619                end if
2620
2621                if (vNam(varn) .EQ. "wpup" .OR. vNam(varn) .EQ. "w"+dq+"u"+dq) then
2622                  miniwpup=min(data(varn,:,min_z_int:max_z_int))-min_value
2623                  maxiwpup=max(data(varn,:,min_z_int:max_z_int))+max_value
2624                  if (over .EQ. 1) then
2625                      res@xyDashPattern  = 0
2626                      plot_wpup = gsn_csm_xy(wks,data(varn,:,:),z,res)
2627                  else
2628                      res@gsnLeftString      = " "
2629                      res@tiXAxisString      = "("+unit(varn)+")"
2630                      res@gsnRightString     = " "
2631                      if (xs .EQ. -1) then
2632                        res@trXMinF            = miniwpup
2633                      else
2634                        res@trXMinF            = xs
2635                      end if
2636                      if (xe .EQ. -1) then
2637                        res@trXMaxF            = maxiwpup
2638                      else
2639                        res@trXMaxF            = xe 
2640                      end if
2641                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2642                  end if
2643                end if
2644                if (vNam(varn) .EQ. "wsus" .OR. vNam(varn) .EQ. "w*u*") then
2645                  miniwsus=min(data(varn,:,min_z_int:max_z_int))-min_value
2646                  maxiwsus=max(data(varn,:,min_z_int:max_z_int))+max_value
2647                  if (over .EQ. 1) then
2648                      res@xyDashPattern  = 1
2649                      plot_wsus = gsn_csm_xy(wks,data(varn,:,:),z,res)
2650                  else
2651                      res@gsnLeftString      = " "
2652                      res@tiXAxisString      = "("+unit(varn)+")"
2653                      res@gsnRightString     = " "
2654                      if (xs .EQ. -1) then
2655                        res@trXMinF            = miniwsus
2656                      else
2657                        res@trXMinF            = xs
2658                      end if
2659                      if (xe .EQ. -1) then     
2660                        res@trXMaxF            = maxiwsus
2661                      else
2662                        res@trXMaxF            = xe 
2663                      end if
2664                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2665                  end if
2666                end if
2667                if (vNam(varn) .EQ. "wu") then
2668                  miniwu=min(data(varn,:,min_z_int:max_z_int))-min_value
2669                  maxiwu=max(data(varn,:,min_z_int:max_z_int))+max_value
2670                  if (over .EQ. 1) then
2671                      res@xyDashPattern  = 2
2672                      plot_wu = gsn_csm_xy(wks,data(varn,:,:),z,res)
2673                  else
2674                      res@gsnLeftString      = " "
2675                      res@tiXAxisString      = "("+unit(varn)+")"
2676                      res@gsnRightString     = " "
2677                      if (xs .EQ. -1) then
2678                        res@trXMinF            = miniwu
2679                      else
2680                        res@trXMinF            = xs
2681                      end if
2682                      if (xe .EQ. -1) then
2683                        res@trXMaxF            = maxiwu
2684                      else
2685                        res@trXMaxF            = xe 
2686                      end if
2687                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2688                  end if
2689                end if
2690
2691                if (vNam(varn) .EQ. "wpvp" .OR. vNam(varn) .EQ. "w"+dq+"v"+dq) then
2692                  miniwpvp=min(data(varn,:,min_z_int:max_z_int))-min_value
2693                  maxiwpvp=max(data(varn,:,min_z_int:max_z_int))+max_value
2694                  if (over .EQ. 1) then
2695                      res@xyDashPattern  = 0
2696                      plot_wpvp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2697                  else
2698                      res@gsnLeftString      = " "
2699                      res@tiXAxisString      = "("+unit(varn)+")"
2700                      res@gsnRightString     = " "
2701                      if (xs .EQ. -1) then
2702                        res@trXMinF            = miniwpvp
2703                      else
2704                        res@trXMinF            = xs
2705                      end if
2706                      if (xe .EQ. -1) then     
2707                        res@trXMaxF            = maxiwpvp
2708                      else
2709                        res@trXMaxF            = xe 
2710                      end if
2711                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2712                  end if
2713                end if
2714                if (vNam(varn) .EQ. "wsvs" .OR. vNam(varn) .EQ. "w*v*") then
2715                  miniwsvs=min(data(varn,:,min_z_int:max_z_int))-min_value
2716                  maxiwsvs=max(data(varn,:,min_z_int:max_z_int))+max_value
2717                  if (over .EQ. 1) then
2718                      res@xyDashPattern  = 1
2719                      plot_wsvs = gsn_csm_xy(wks,data(varn,:,:),z,res)
2720                  else
2721                      res@gsnLeftString      = " "
2722                      res@tiXAxisString      = "("+unit(varn)+")"
2723                      res@gsnRightString     = " "
2724                      if (xs .EQ. -1) then
2725                        res@trXMinF            = miniwsvs
2726                      else
2727                        res@trXMinF            = xs
2728                      end if
2729                      if (xe .EQ. -1) then     
2730                        res@trXMaxF            = maxiwsvs
2731                      else
2732                        res@trXMaxF            = xe 
2733                      end if
2734                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2735                  end if
2736                end if
2737                if (vNam(varn) .EQ. "wv") then
2738                  miniwv=min(data(varn,:,min_z_int:max_z_int))-min_value
2739                  maxiwv=max(data(varn,:,min_z_int:max_z_int))+max_value
2740                  if (over .EQ. 1) then
2741                      res@xyDashPattern  = 2
2742                      plot_wv = gsn_csm_xy(wks,data(varn,:,:),z,res)
2743                  else
2744                      res@gsnLeftString      = " "
2745                      res@tiXAxisString      = "("+unit(varn)+")"
2746                      res@gsnRightString     = " "
2747                      if (xs .EQ. -1) then
2748                        res@trXMinF            = miniwv
2749                      else
2750                        res@trXMinF            = xs
2751                      end if
2752                      if (xe .EQ. -1) then       
2753                        res@trXMaxF            = maxiwv
2754                      else
2755                        res@trXMaxF            = xe 
2756                      end if
2757                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2758                  end if
2759                end if
2760
2761                if (vNam(varn) .EQ. "wpptp" .OR. vNam(varn) \
2762                  .EQ. "w"+dq+"pt"+dq) then
2763                  miniwpptp=min(data(varn,:,min_z_int:max_z_int))-min_value
2764                  maxiwpptp=max(data(varn,:,min_z_int:max_z_int))+max_value
2765                  if (over .EQ. 1) then
2766                      res@xyDashPattern  = 0
2767                      plot_wpptp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2768                  else
2769                      res@gsnLeftString      = " "
2770                      res@tiXAxisString      = "("+unit(varn)+")"
2771                      res@gsnRightString     = " "
2772                      if (xs .EQ. -1) then
2773                        res@trXMinF            = miniwpptp
2774                      else
2775                        res@trXMinF            = xs
2776                      end if
2777                      if (xe .EQ. -1) then       
2778                        res@trXMaxF            = maxiwpptp
2779                      else
2780                        res@trXMaxF            = xe 
2781                      end if
2782                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2783                  end if
2784                end if
2785                if (vNam(varn) .EQ. "wspts" .OR. vNam(varn) .EQ. "w*pt*") then
2786                  miniwspts=min(data(varn,:,min_z_int:max_z_int))-min_value
2787                  maxiwspts=max(data(varn,:,min_z_int:max_z_int))+max_value
2788                  if (over .EQ. 1) then
2789                      res@xyDashPattern  = 1
2790                      plot_wspts = gsn_csm_xy(wks,data(varn,:,:),z,res)
2791                  else
2792                      res@gsnLeftString      = " "
2793                      res@tiXAxisString      = "("+unit(varn)+")"
2794                      res@gsnRightString     = " "
2795                      if (xs .EQ. -1) then
2796                        res@trXMinF            = miniwspts
2797                      else
2798                        res@trXMinF            = xs
2799                      end if
2800                      if (xe .EQ. -1) then       
2801                        res@trXMaxF            = maxiwspts
2802                      else
2803                        res@trXMaxF            = xe 
2804                      end if
2805                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2806                  end if
2807                end if
2808                if (vNam(varn) .EQ. "wpt") then       
2809                  miniwpt=min(data(varn,:,min_z_int:max_z_int))-min_value
2810                  maxiwpt=max(data(varn,:,min_z_int:max_z_int))+max_value
2811                  if (over .EQ. 1) then
2812                      res@xyDashPattern  = 2
2813                      plot_wpt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2814                  else
2815                      res@gsnLeftString      = " "
2816                      res@tiXAxisString      = "("+unit(varn)+")"
2817                      res@gsnRightString     = " "
2818                      if (xs .EQ. -1) then
2819                        res@trXMinF            = miniwpt
2820                      else
2821                        res@trXMinF            = xs
2822                      end if
2823                      if (xe .EQ. -1) then       
2824                        res@trXMaxF            = maxiwpt
2825                      else
2826                        res@trXMaxF            = xe 
2827                      end if
2828                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2829                  end if
2830                end if
2831
2832                if (vNam(varn) .EQ. "wsptsBC".OR. vNam(varn) .EQ. "w*pt*BC" ) then
2833                  miniwsptsBC=min(data(varn,:,min_z_int:max_z_int))-min_value
2834                  maxiwsptsBC=max(data(varn,:,min_z_int:max_z_int))+max_value
2835                  if (over .EQ. 1) then
2836                      res@xyDashPattern  = 0
2837                      plot_wsptsBC = gsn_csm_xy(wks,data(varn,:,:),z,res)
2838                  else
2839                      res@gsnLeftString      = " "
2840                      res@tiXAxisString      = "("+unit(varn)+")"
2841                      res@gsnRightString     = " "
2842                      if (xs .EQ. -1) then
2843                        res@trXMinF            = miniwsptsBC
2844                      else
2845                        res@trXMinF            = xs
2846                      end if
2847                      if (xe .EQ. -1) then       
2848                        res@trXMaxF            = maxiwsptsBC
2849                      else
2850                        res@trXMaxF            = xe 
2851                      end if
2852                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2853                  end if
2854                end if             
2855                if (vNam(varn) .EQ. "wptBC") then
2856                  miniwptBC=min(data(varn,:,min_z_int:max_z_int))-min_value
2857                  maxiwptBC=max(data(varn,:,min_z_int:max_z_int))+max_value
2858                  if (over .EQ. 1) then
2859                      res@xyDashPattern  = 1
2860                      plot_wptBC = gsn_csm_xy(wks,data(varn,:,:),z,res)
2861                  else
2862                      res@gsnLeftString      = " "
2863                      res@tiXAxisString      = "("+unit(varn)+")"
2864                      res@gsnRightString     = " "
2865                      if (xs .EQ. -1) then
2866                        res@trXMinF            = miniwptBC
2867                      else
2868                        res@trXMinF            = xs
2869                      end if
2870                      if (xe .EQ. -1) then       
2871                        res@trXMaxF            = maxiwptBC
2872                      else
2873                        res@trXMaxF            = xe 
2874                      end if
2875                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2876                  end if
2877                end if
2878
2879                if (vNam(varn) .EQ. "wpvptp" .OR. vNam(varn) \
2880                    .EQ. "w"+dq+"vpt"+dq) then
2881                  miniwpvptp=min(data(varn,:,min_z_int:max_z_int))-min_value
2882                  maxiwpvptp=max(data(varn,:,min_z_int:max_z_int))+max_value
2883                  if (over .EQ. 1) then
2884                      res@xyDashPattern  = 0
2885                      plot_wpvptp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2886                  else
2887                      res@gsnLeftString      = " "
2888                      res@tiXAxisString      = "("+unit(varn)+")"
2889                      res@gsnRightString     = " "
2890                      if (xs .EQ. -1) then
2891                        res@trXMinF            = miniwpvptp
2892                      else
2893                        res@trXMinF            = xs
2894                      end if
2895                      if (xe .EQ. -1) then       
2896                        res@trXMaxF            = maxiwpvptp
2897                      else
2898                        res@trXMaxF            = xe 
2899                      end if
2900                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2901                  end if
2902                end if
2903                if (vNam(varn) .EQ. "wsvpts" .OR. vNam(varn) .EQ. "w*vpt*") then
2904                  miniwsvpts=min(data(varn,:,min_z_int:max_z_int))-min_value
2905                  maxiwsvpts=max(data(varn,:,min_z_int:max_z_int))+max_value
2906                  if (over .EQ. 1) then
2907                      res@xyDashPattern  = 1
2908                      plot_wsvpts = gsn_csm_xy(wks,data(varn,:,:),z,res)
2909                  else
2910                      res@gsnLeftString      = " "
2911                      res@tiXAxisString      = "("+unit(varn)+")"
2912                      res@gsnRightString     = " "
2913                      if (xs .EQ. -1) then
2914                        res@trXMinF            = miniwsvpts
2915                      else
2916                        res@trXMinF            = xs
2917                      end if
2918                      if (xe .EQ. -1) then       
2919                        res@trXMaxF            = maxiwsvpts
2920                      else
2921                        res@trXMaxF            = xe 
2922                      end if
2923                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2924                  end if
2925                end if
2926                if (vNam(varn) .EQ. "wvpt") then
2927                  miniwvpt=min(data(varn,:,min_z_int:max_z_int))-min_value
2928                  maxiwvpt=max(data(varn,:,min_z_int:max_z_int))+max_value
2929                  if (over .EQ. 1) then
2930                      res@xyDashPattern  = 2
2931                      plot_wvpt = gsn_csm_xy(wks,data(varn,:,:),z,res)
2932                  else
2933                      res@gsnLeftString      = " "
2934                      res@tiXAxisString      = "("+unit(varn)+")"
2935                      res@gsnRightString     = " "
2936                      if (xs .EQ. -1) then
2937                        res@trXMinF            = miniwvpt
2938                      else
2939                        res@trXMinF            = xs
2940                      end if
2941                      if (xe .EQ. -1) then       
2942                        res@trXMaxF            = maxiwvpt
2943                      else
2944                        res@trXMaxF            = xe 
2945                      end if
2946                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2947                  end if
2948                end if
2949
2950                if (vNam(varn) .EQ. "wpqp" .OR. vNam(varn) .EQ. "w"+dq+"q"+dq) then
2951                  miniwpqp=min(data(varn,:,min_z_int:max_z_int))-min_value
2952                  maxiwpqp=max(data(varn,:,min_z_int:max_z_int))+max_value
2953                  if (over .EQ. 1) then
2954                      res@xyDashPattern  = 0
2955                      plot_wpqp = gsn_csm_xy(wks,data(varn,:,:),z,res)
2956                  else
2957                      res@gsnLeftString      = " "
2958                      res@tiXAxisString      = "("+unit(varn)+")"
2959                      res@gsnRightString     = " "
2960                      if (xs .EQ. -1) then
2961                        res@trXMinF            = miniwpqp
2962                      else
2963                        res@trXMinF            = xs
2964                      end if
2965                      if (xe .EQ. -1) then       
2966                        res@trXMaxF            = maxiwpqp
2967                      else
2968                        res@trXMaxF            = xe 
2969                      end if
2970                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2971                  end if
2972                end if
2973                if (vNam(varn) .EQ. "wsqs".OR. vNam(varn) .EQ. "w*s*" ) then
2974                  miniwsqs=min(data(varn,:,min_z_int:max_z_int))-min_value
2975                  maxiwsqs=max(data(varn,:,min_z_int:max_z_int))+max_value
2976                  if (over .EQ. 1) then
2977                      res@xyDashPattern  = 1
2978                      plot_wsqs = gsn_csm_xy(wks,data(varn,:,:),z,res)
2979                  else
2980                      res@gsnLeftString      = " "
2981                      res@tiXAxisString      = "("+unit(varn)+")"
2982                      res@gsnRightString     = " "
2983                      if (xs .EQ. -1) then
2984                        res@trXMinF            = miniwsqs
2985                      else
2986                        res@trXMinF            = xs
2987                      end if
2988                      if (xe .EQ. -1) then       
2989                        res@trXMaxF            = maxiwsqs
2990                      else
2991                        res@trXMaxF            = xe 
2992                      end if
2993                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
2994                  end if
2995                end if
2996                if (vNam(varn) .EQ. "wq") then
2997                  miniwq=min(data(varn,:,min_z_int:max_z_int))-min_value
2998                  maxiwq=max(data(varn,:,min_z_int:max_z_int))+max_value
2999                  if (over .EQ. 1) then
3000                      res@xyDashPattern  = 2
3001                      plot_wq = gsn_csm_xy(wks,data(varn,:,:),z,res)
3002                  else
3003                      res@gsnLeftString      = " "
3004                      res@tiXAxisString      = "("+unit(varn)+")"
3005                      res@gsnRightString     = " "
3006                      if (xs .EQ. -1) then
3007                        res@trXMinF            = miniwq
3008                      else
3009                        res@trXMinF            = xs
3010                      end if
3011                      if (xe .EQ. -1) then       
3012                        res@trXMaxF            = maxiwq
3013                      else
3014                        res@trXMaxF            = xe 
3015                      end if
3016                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3017                  end if
3018                end if
3019
3020                if (vNam(varn) .EQ. "wpqvp" .OR. \
3021                  vNam(varn) .EQ. "w"+dq+"qv"+dq) then
3022                  miniwpqvp=min(data(varn,:,min_z_int:max_z_int))-min_value
3023                  maxiwpqvp=max(data(varn,:,min_z_int:max_z_int))+max_value
3024                  if (over .EQ. 1) then
3025                      res@xyDashPattern  = 0
3026                      plot_wpqvp = gsn_csm_xy(wks,data(varn,:,:),z,res)
3027                  else
3028                      res@gsnLeftString      = " "
3029                      res@tiXAxisString      = "("+unit(varn)+")"
3030                      res@gsnRightString     = " "
3031                      if (xs .EQ. -1) then
3032                        res@trXMinF            = miniwpqvp
3033                      else
3034                        res@trXMinF            = xs
3035                      end if
3036                      if (xe .EQ. -1) then       
3037                        res@trXMaxF            = maxiwpqvp
3038                      else
3039                        res@trXMaxF            = xe 
3040                      end if
3041                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3042                  end if
3043                end if
3044                if (vNam(varn) .EQ. "wsqvs" .OR. vNam(varn) .EQ. "w*qv*") then
3045                  miniwsqvs=min(data(varn,:,min_z_int:max_z_int))-min_value
3046                  maxiwsqvs=max(data(varn,:,min_z_int:max_z_int))+max_value
3047                  if (over .EQ. 1) then
3048                      res@xyDashPattern  = 1
3049                      plot_wsqvs = gsn_csm_xy(wks,data(varn,:,:),z,res)
3050                  else
3051                      res@gsnLeftString      = " "
3052                      res@tiXAxisString      = "("+unit(varn)+")"
3053                      res@gsnRightString     = " "
3054                      if (xs .EQ. -1) then
3055                        res@trXMinF            = miniwsqvs
3056                      else
3057                        res@trXMinF            = xs
3058                      end if
3059                      if (xe .EQ. -1) then       
3060                        res@trXMaxF            = maxiwsqvs
3061                      else
3062                        res@trXMaxF            = xe 
3063                      end if
3064                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3065                  end if
3066                end if
3067                if (vNam(varn) .EQ. "wqv") then
3068                  miniwqv=min(data(varn,:,min_z_int:max_z_int))-min_value
3069                  maxiwqv=max(data(varn,:,min_z_int:max_z_int))+max_value
3070                  if (over .EQ. 1) then
3071                      res@xyDashPattern  = 2
3072                      plot_wqv = gsn_csm_xy(wks,data(varn,:,:),z,res)
3073                  else
3074                      res@gsnLeftString      = " "
3075                      res@tiXAxisString      = "("+unit(varn)+")"
3076                      res@gsnRightString     = " "
3077                      if (xs .EQ. -1) then
3078                        res@trXMinF            = miniwqv
3079                      else
3080                        res@trXMinF            = xs
3081                      end if
3082                      if (xe .EQ. -1) then       
3083                        res@trXMaxF            = maxiwqv
3084                      else
3085                        res@trXMaxF            = xe 
3086                      end if
3087                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3088                  end if
3089                end if
3090
3091                if (vNam(varn) .EQ. "wpsp" .OR. vNam(varn) .EQ. "w"+dq+"s"+dq) then
3092                  miniwpsp=min(data(varn,:,min_z_int:max_z_int))-min_value
3093                  maxiwpsp=max(data(varn,:,min_z_int:max_z_int))+max_value
3094                  if (over .EQ. 1) then
3095                      res@xyDashPattern  = 0
3096                      plot_wpsp = gsn_csm_xy(wks,data(varn,:,:),z,res)
3097                  else
3098                      res@gsnLeftString      = " "
3099                      res@tiXAxisString      = "("+unit(varn)+")"
3100                      res@gsnRightString     = " "
3101                      if (xs .EQ. -1) then
3102                        res@trXMinF            = miniwpsp
3103                      else
3104                        res@trXMinF            = xs
3105                      end if
3106                      if (xe .EQ. -1) then       
3107                        res@trXMaxF            = maxiwpsp
3108                      else
3109                        res@trXMaxF            = xe 
3110                      end if
3111                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3112                  end if
3113                end if
3114                if (vNam(varn) .EQ. "wsss" .OR. vNam(varn) .EQ. "w*s*" ) then
3115                  miniwsss=min(data(varn,:,min_z_int:max_z_int))-min_value
3116                  maxiwsss=max(data(varn,:,min_z_int:max_z_int))+max_value
3117                  if (over .EQ. 1) then
3118                      res@xyDashPattern  = 1
3119                      plot_wsss = gsn_csm_xy(wks,data(varn,:,:),z,res)
3120                  else
3121                      res@gsnLeftString      = " "
3122                      res@tiXAxisString      = "("+unit(varn)+")"
3123                      res@gsnRightString     = " "
3124                      if (xs .EQ. -1) then
3125                        res@trXMinF            = miniwsss
3126                      else
3127                        res@trXMinF            = xs
3128                      end if
3129                      if (xe .EQ. -1) then       
3130                        res@trXMaxF            = maxiwsss
3131                      else
3132                        res@trXMaxF            = xe 
3133                      end if
3134                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3135                  end if
3136                end if
3137                if (vNam(varn) .EQ. "ws") then
3138                  miniws=min(data(varn,:,min_z_int:max_z_int))-min_value
3139                  maxiws=max(data(varn,:,min_z_int:max_z_int))+max_value
3140                  if (over .EQ. 1) then
3141                      res@xyDashPattern  = 2
3142                      plot_ws = gsn_csm_xy(wks,data(varn,:,:),z,res)
3143                  else
3144                      res@gsnLeftString      = " "
3145                      res@tiXAxisString      = "("+unit(varn)+")"
3146                      res@gsnRightString     = " "
3147                      if (xs .EQ. -1) then
3148                        res@trXMinF            = miniws
3149                      else
3150                        res@trXMinF            = xs
3151                      end if
3152                      if (xe .EQ. -1) then       
3153                        res@trXMaxF            = maxiws
3154                      else
3155                        res@trXMaxF            = xe 
3156                      end if
3157                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3158                  end if
3159                end if
3160
3161                if (vNam(varn) .EQ. "wpsap" .OR. \
3162                    vNam(varn) .EQ. "w"+dq+"sa"+dq) then
3163                  miniwpsap=min(data(varn,:,min_z_int:max_z_int))-min_value
3164                  maxiwpsap=max(data(varn,:,min_z_int:max_z_int))+max_value
3165                  if (over .EQ. 1) then
3166                      res@xyDashPattern  = 0
3167                      plot_wpsap = gsn_csm_xy(wks,data(varn,:,:),z,res)
3168                  else
3169                      res@gsnLeftString      = " "
3170                      res@tiXAxisString      = "("+unit(varn)+")"
3171                      res@gsnRightString     = " "
3172                      if (xs .EQ. -1) then
3173                        res@trXMinF            = miniwpsap
3174                      else
3175                        res@trXMinF            = xs
3176                      end if
3177                      if (xe .EQ. -1) then       
3178                        res@trXMaxF            = maxiwpsap
3179                      else
3180                        res@trXMaxF            = xe 
3181                      end if
3182                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3183                  end if
3184                end if
3185                if (vNam(varn) .EQ. "wssas" .OR. vNam(varn) .EQ. "w*sa*") then
3186                  miniwssas=min(data(varn,:,min_z_int:max_z_int))-min_value
3187                  maxiwssas=max(data(varn,:,min_z_int:max_z_int))+max_value
3188                  if (over .EQ. 1) then
3189                      res@xyDashPattern  = 1
3190                      plot_wssas = gsn_csm_xy(wks,data(varn,:,:),z,res)
3191                  else
3192                      res@gsnLeftString      = " "
3193                      res@tiXAxisString      = "("+unit(varn)+")"
3194                      res@gsnRightString     = " "
3195                      if (xs .EQ. -1) then
3196                        res@trXMinF            = miniwssas
3197                      else
3198                        res@trXMinF            = xs
3199                      end if
3200                      if (xe .EQ. -1) then       
3201                        res@trXMaxF            = maxiwssas
3202                      else
3203                        res@trXMaxF            = xe 
3204                      end if
3205                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3206                  end if
3207                end if
3208                if (vNam(varn) .EQ. "wsa") then
3209                  miniwsa=min(data(varn,:,min_z_int:max_z_int))-min_value
3210                  maxiwsa=max(data(varn,:,min_z_int:max_z_int))+max_value
3211                  if (over .EQ. 1) then
3212                      res@xyDashPattern  = 2
3213                      plot_wsa = gsn_csm_xy(wks,data(varn,:,:),z,res)
3214                  else
3215                      res@gsnLeftString      = " "
3216                      res@tiXAxisString      = "("+unit(varn)+")"
3217                      res@gsnRightString     = " "
3218                      if (xs .EQ. -1) then
3219                        res@trXMinF            = miniwsa
3220                      else
3221                        res@trXMinF            = xs
3222                      end if
3223                      if (xe .EQ. -1) then       
3224                        res@trXMaxF            = maxiwsa
3225                      else
3226                        res@trXMaxF            = xe 
3227                      end if
3228                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3229                  end if
3230                end if
3231
3232                if (vNam(varn) .EQ. "us2" .OR. vNam(varn) .EQ. "u*2") then
3233                  minius2=min(data(varn,:,min_z_int:max_z_int))-min_value
3234                  maxius2=max(data(varn,:,min_z_int:max_z_int))+max_value
3235                  if (over .EQ. 1) then
3236                      res@xyDashPattern  = 0
3237                      plot_us2 = gsn_csm_xy(wks,data(varn,:,:),z,res)
3238                  else
3239                      res@gsnLeftString      = " "
3240                      res@tiXAxisString      = "("+unit(varn)+")"
3241                      res@gsnRightString     = " "
3242                      if (xs .EQ. -1) then
3243                        res@trXMinF            = minius2
3244                      else
3245                        res@trXMinF            = xs
3246                      end if
3247                      if (xe .EQ. -1) then       
3248                        res@trXMaxF            = maxius2
3249                      else
3250                        res@trXMaxF            = xe 
3251                      end if
3252                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3253                  end if
3254                end if
3255                if (vNam(varn) .EQ. "vs2" .OR. vNam(varn) .EQ. "v*2") then
3256                  minivs2=min(data(varn,:,min_z_int:max_z_int))-min_value
3257                  maxivs2=max(data(varn,:,min_z_int:max_z_int))+max_value
3258                  if (over .EQ. 1) then
3259                      res@xyDashPattern  = 1
3260                      plot_vs2 = gsn_csm_xy(wks,data(varn,:,:),z,res)
3261                  else
3262                      res@gsnLeftString      = " "
3263                      res@tiXAxisString      = "("+unit(varn)+")"
3264                      res@gsnRightString     = " "
3265                      if (xs .EQ. -1) then
3266                        res@trXMinF            = minivs2
3267                      else
3268                        res@trXMinF            = xs
3269                      end if
3270                      if (xe .EQ. -1) then       
3271                        res@trXMaxF            = maxivs2
3272                      else
3273                        res@trXMaxF            = xe 
3274                      end if
3275                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3276                  end if
3277                end if
3278                if (vNam(varn) .EQ. "ws2" .OR. vNam(varn) .EQ. "w*2") then
3279                  miniws2=min(data(varn,:,min_z_int:max_z_int))-min_value
3280                  maxiws2=max(data(varn,:,min_z_int:max_z_int))+max_value
3281                  if (over .EQ. 1) then
3282                      res@xyDashPattern  = 2
3283                      plot_ws2 = gsn_csm_xy(wks,data(varn,:,:),z,res)
3284                  else
3285                      res@gsnLeftString      = " "
3286                      res@tiXAxisString      = "("+unit(varn)+")"
3287                      res@gsnRightString     = " "
3288                      if (xs .EQ. -1) then
3289                        res@trXMinF            = miniws2
3290                      else
3291                        res@trXMinF            = xs
3292                      end if
3293                      if (xe .EQ. -1) then       
3294                        res@trXMaxF            = maxiws2
3295                      else
3296                        res@trXMaxF            = xe 
3297                      end if
3298                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3299                  end if
3300                end if
3301
3302                if (vNam(varn) .EQ. "wsususodz" .OR. \
3303                    vNam(varn) .EQ. "w*u*u*:dz") then
3304                  miniwsususodz=min(data(varn,:,min_z_int:max_z_int))-min_value
3305                  maxiwsususodz=max(data(varn,:,min_z_int:max_z_int))+max_value
3306                  if (over .EQ. 1) then
3307                      res@xyDashPattern  = 0
3308                      plot_wsususodz = gsn_csm_xy(wks,data(varn,:,:),z,res)
3309                  else
3310                      res@gsnLeftString      = " "
3311                      res@tiXAxisString      = "("+unit(varn)+")"
3312                      res@gsnRightString     = " "
3313                      if (xs .EQ. -1) then
3314                        res@trXMinF            = miniwsususodz
3315                      else
3316                        res@trXMinF            = xs
3317                      end if
3318                      if (xe .EQ. -1) then       
3319                        res@trXMaxF            = maxiwsususodz
3320                      else
3321                        res@trXMaxF            = xe 
3322                      end if
3323                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3324                  end if 
3325                end if
3326                if (vNam(varn) .EQ. "wspsodz" .OR. vNam(varn) .EQ. "w*p*:dz") then
3327                  miniwspsodz=min(data(varn,:,min_z_int:max_z_int))-min_value
3328                  maxiwspsodz=max(data(varn,:,min_z_int:max_z_int))+max_value
3329                  if (over .EQ. 1) then
3330                      res@xyDashPattern  = 1
3331                      plot_wspsodz = gsn_csm_xy(wks,data(varn,:,:),z,res)
3332                  else
3333                      res@gsnLeftString      = " "
3334                      res@tiXAxisString      = "("+unit(varn)+")"
3335                      res@gsnRightString     = " "
3336                      if (xs .EQ. -1) then
3337                        res@trXMinF            = miniwspsodz
3338                      else
3339                        res@trXMinF            = xs
3340                      end if
3341                      if (xe .EQ. -1) then       
3342                        res@trXMaxF            = maxiwspsodz
3343                      else
3344                        res@trXMaxF            = xe 
3345                      end if
3346                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3347                  end if
3348                end if
3349                if (vNam(varn) .EQ. "wpeodz" .OR. \
3350                    vNam(varn) .EQ. "w"+dq+"p:dz") then
3351                  miniwpeodz=min(data(varn,:,min_z_int:max_z_int))-min_value
3352                  maxiwpeodz=max(data(varn,:,min_z_int:max_z_int))+max_value
3353                  if (over .EQ. 1) then
3354                      res@xyDashPattern  = 2
3355                      plot_wpeodz = gsn_csm_xy(wks,data(varn,:,:),z,res)
3356                  else
3357                      res@gsnLeftString      = " "
3358                      res@tiXAxisString      = "("+unit(varn)+")"
3359                      res@gsnRightString     = " "
3360                      if (xs .EQ. -1) then
3361                        res@trXMinF            = miniwpeodz
3362                      else
3363                        res@trXMinF            = xs
3364                      end if
3365                      if (xe .EQ. -1) then       
3366                        res@trXMaxF            = maxiwpeodz
3367                      else
3368                        res@trXMaxF            = xe 
3369                      end if
3370                      plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)
3371                  end if
3372                end if
3373
3374                ; ***************************************************
3375                ; legend for each plot
3376                ; ***************************************************
3377               
3378                if (over .EQ. 0) then
3379
3380                   label = " " + vNam(varn)
3381
3382                   lgres1                    = True
3383                   lgMonoDashIndex           = False
3384                   lgres1@lgDashIndexes      = (/0,1,2/)
3385                   lgres1@lgLabelFont        = "helvetica"   
3386                   lgres1@lgLabelFontHeightF = font_size_legend * 10.0 * 0.7
3387                   lgres1@vpWidthF           = 0.14           
3388                   lgres1@vpHeightF          = 0.1 / 3.0 + 0.01         
3389               
3390                   lbid = gsn_create_legend(wks,1,label,lgres1)       
3391
3392                   amres = True
3393                   amres@amParallelPosF   = -0.5                   
3394                   amres@amOrthogonalPosF = -0.5   
3395                   amres@amJust           = "TopLeft"
3396
3397                   annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3398     
3399                end if
3400
3401            if (no_files .GT. 1) then
3402;               print("nof="+nof+" und n="+n)
3403                multi_plot(nof,n)=plot(n)
3404                max_nof(nof,n)=max(data(varn,:,min_z_int:max_z_int))
3405                min_nof(nof,n)=min(data(varn,:,min_z_int:max_z_int))
3406                name(nof,n)   =vNam(varn)
3407                unit_(nof,n)  =unit(varn)
3408            end if
3409            if (over .EQ. 0) then
3410                n=n+1 
3411            end if 
3412            if (prof3d .EQ. 0)then   
3413                varn=varn+1
3414            end if
3415            delete(temp)
3416            delete(temp_att)
3417          end if         
3418      end do
3419      if (no_files .GT. 1) then
3420          delete(vNam)
3421          delete(files)
3422      end if
3423     
3424      end do
3425
3426      ;#########ENDE DO LOOP FOR no_files GT 1#############
3427
3428      if (isStrSubset(data@long_name," SR " ) .and. over_remind) then
3429          print(" ")
3430          print("If you have outputs of statistic regions you cannot overlay "+\
3431                "variables;")
3432          print("'over' is set to 0" )
3433          print(" ")
3434          over = 0
3435      end if
3436
3437      if (count_var .EQ. 0) then
3438          print(" ")
3439          print("The variables 'var="+var+"'" )
3440          print("do not exist on your input file;")
3441          print("be sure to have one comma before and after each variable")
3442          print(" ")
3443          exit       
3444      end if
3445     
3446      if (no_files .GT. 1) then
3447          multi_legend=new(6,string)
3448          string_len=new(6,integer)
3449          multi_dash=new(no_files,string)
3450          multi_legend(0)="  "+name_legend_1
3451          string_len(0)=strlen(multi_legend(0))
3452          multi_legend(1)="  "+name_legend_2
3453          string_len(1)=strlen(multi_legend(1))
3454          multi_legend(2)="  "+name_legend_3
3455          string_len(2)=strlen(multi_legend(2))
3456          multi_legend(3)="  "+name_legend_4
3457          string_len(3)=strlen(multi_legend(3))
3458          multi_legend(4)="  "+name_legend_5
3459          string_len(4)=strlen(multi_legend(4))
3460          multi_legend(5)="  "+name_legend_6
3461          string_len(5)=strlen(multi_legend(5))
3462          do ml=1,no_files
3463            multi_dash(ml-1)=ml-1
3464          end do
3465          delete(plot)
3466          plot = new(dim,graphic)
3467          do pl=0,n-1
3468            plot0 = new(1,graphic)
3469            res@trXMinF = min(min_nof(:,pl))
3470            res@trXMaxF = max(max_nof(:,pl))
3471            res@xyDashPattern=0
3472            res@gsnLeftString  = " "; name(0,pl)
3473            res@gsnRightString = " ";unit_(0,pl)
3474            res@tiXAxisString  = "(" + unit_(0,pl) + ")"
3475            data_0(:,:) = min(min_nof(:,pl))
3476
3477            if ( (mod(no_plots - pl, no_rows * no_columns) .EQ. 1) .OR. (no_rows * no_columns .EQ. 1) ) then
3478                res@pmLegendDisplayMode  = "Always"
3479            else
3480                res@pmLegendDisplayMode  = "Never"
3481            end if
3482
3483            plot0 = gsn_csm_xy(wks,data_0(:,:),z_(pl,:),res)
3484
3485            ; ***************************************************
3486            ; legend for combined plot (File 1, File 2, ...)
3487            ; ***************************************************
3488
3489            if ( (mod(no_plots - pl, no_rows * no_columns) .EQ. 1) .OR. (no_rows * no_columns .EQ. 1) )then
3490                lgres                    = True
3491                lgMonoDashIndex          = False
3492                lgres@lgLabelFont        = "helvetica"   
3493                lgres@lgLabelFontHeightF = font_size_legend * 10.0       
3494                lgres@vpWidthF           = max(string_len)*0.02;0.015           
3495                lgres@vpHeightF          = 0.035*no_files   
3496                lgres@lgDashIndexes      = multi_dash(no_files-1:0)
3497                lgres@lgLineDashSegLenF  = 0.075
3498                lbid = gsn_create_legend(\
3499                                    wks,no_files,multi_legend(no_files-1:0),lgres)
3500
3501                amres = True
3502                amres@amParallelPosF   = -0.5       
3503                amres@amOrthogonalPosF = -0.55
3504                amres@amJust           = "BottomLeft"
3505                annoid1 = gsn_add_annotation(plot0,lbid,amres)
3506            end if
3507
3508            do plo=0,no_files-1
3509                overlay(plot0,multi_plot(plo,pl))
3510                plot(pl)=plot0
3511            end do
3512            delete(plot0)
3513          end do
3514      end if
3515
3516      if (count_var .EQ. 0) then
3517          print(" ")
3518          print("Select a variable 'var=' or use the default value")
3519          print(" ")
3520          print("Your selection '"+var+"' does not exist on the input file")
3521          print(" ")
3522          exit
3523      end if
3524
3525      if (over .EQ. 1 ) then
3526
3527          overlay(plot_u,plot_v)
3528          overlay(plot_u,plot_w)
3529          u=0
3530          overlay(plot_pt,plot_vpt)
3531          overlay(plot_pt,plot_lpt)
3532          pt=0
3533          overlay(plot_q,plot_qv)
3534          overlay(plot_q,plot_ql)
3535          q=0
3536          overlay(plot_e,plot_es)
3537          e=0
3538          overlay(plot_km,plot_kh)
3539          km=0
3540          overlay(plot_wpup,plot_wsus)
3541          overlay(plot_wpup,plot_wu)
3542          wpup=0
3543          overlay(plot_wpvp,plot_wsvs)
3544          overlay(plot_wpvp,plot_wv)
3545          wpvp=0
3546          overlay(plot_wpptp,plot_wspts)
3547          overlay(plot_wpptp,plot_wpt)
3548          wpptp=0
3549          overlay(plot_wsptsBC,plot_wptBC)
3550          wsptsBC=0
3551          overlay(plot_wpvptp,plot_wsvpts)
3552          overlay(plot_wpvptp,plot_wvpt)
3553          wpvptp=0
3554          overlay(plot_wpqp,plot_wsqs)
3555          overlay(plot_wpqp,plot_wq)
3556          wpqp=0
3557          overlay(plot_wpqvp,plot_wsqvs)
3558          overlay(plot_wpqvp,plot_wqv)
3559          wpqvp=0
3560          overlay(plot_wpsp,plot_wsss)
3561          overlay(plot_wpsp,plot_ws)
3562          wpsp=0
3563          overlay(plot_wpsap,plot_wssas)
3564          overlay(plot_wpsap,plot_wsa)
3565          wpsap=0
3566          overlay(plot_us2,plot_vs2)
3567          overlay(plot_us2,plot_ws2)
3568          us2=0
3569          overlay(plot_wsususodz,plot_wspsodz)
3570          overlay(plot_wsususodz,plot_wpeodz)
3571          wsususodz=0
3572
3573      end if
3574
3575      if (over .EQ. 1) then
3576     
3577          do varn = 0,dim-1   
3578           
3579            check = True
3580         
3581            if (prof3d .EQ. 0) then
3582                if ( isStrSubset( vNam(varn), "time") .OR. \
3583                    isStrSubset( vNam(varn), "NORM")) then
3584                  check = False
3585                end if
3586            else
3587                if ( isStrSubset( vNam(varn), "time") .OR.  \
3588                    isStrSubset( vNam(varn), "zusi") .OR.  \
3589                    isStrSubset( vNam(varn), "zwwi") .OR.  \
3590                    isStrSubset( vNam(varn), "x") .OR.     \
3591                    isStrSubset( vNam(varn), "xu") .OR.    \
3592                    isStrSubset( vNam(varn), "y") .OR.     \
3593                    isStrSubset( vNam(varn), "yv") .OR.    \
3594                    isStrSubset( vNam(varn), "zu_3d") .OR. \
3595                    isStrSubset( vNam(varn), "zw_3d")) then
3596                  check = False
3597                end if
3598            end if
3599
3600            if (var .NE. "all") then     
3601                check = isStrSubset( var,","+vNam(varn)+"," )
3602            end if     
3603
3604            if (check)then
3605
3606                if (prof3d .EQ. 0) then
3607                  if (log_z .EQ. 1) then
3608                      z = f_att->$vNam(varn+1)$(1:dimz-1)
3609                  else
3610                      z = f_att->$vNam(varn+1)$               
3611                  end if
3612                else
3613                  do i=0,b-1           
3614                      if (isStrSubset( a(i),"zu_3d" ))then
3615                        z_v(varn,:) = z_u
3616                        if (log_z .EQ. 1) then
3617                            z = z_v(varn,1:dimz-1)
3618                        else
3619                            z = z_v(varn,:)
3620                        end if
3621                      else
3622                        if (isStrSubset( a(i),"zw_3d" ))then
3623                            z_v(varn,:) = z_w
3624                            if (log_z .EQ. 1) then
3625                              z = z_v(varn,1:dimz-1)
3626                            else
3627                              z = z_v(varn,:)
3628                            end if
3629                        end if                   
3630                      end if
3631                  end do           
3632                end if
3633
3634                z=z/norm_z
3635               
3636                if (abs(min(data(varn,:,min_z_int:max_z_int))) .GT. 10)then
3637                  min_value = abs(0.01*min(data(varn,:,min_z_int:max_z_int)))
3638                  max_value = abs(0.01*max(data(varn,:,min_z_int:max_z_int)))
3639                else
3640                  if (abs(min(data(varn,:,min_z_int:max_z_int))) .LT. 0.01 .AND. \
3641                      abs(max(data(varn,:,min_z_int:max_z_int))) .GT. 0.01)then
3642                      min_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
3643                      max_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
3644                  else
3645                      if (abs(max(data(varn,:,:))) .LT. 0.01 .AND. \
3646                          abs(min(data(varn,:,min_z_int:max_z_int))) .GT. 0.01)then
3647                        min_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
3648                        max_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
3649                      else
3650                        min_value = abs(0.1*min(data(varn,:,min_z_int:max_z_int)))
3651                        max_value = abs(0.1*max(data(varn,:,min_z_int:max_z_int)))
3652                      end if
3653                  end if
3654                end if
3655                if (min(data(varn,:,min_z_int:max_z_int)) .EQ. 0 .AND.\
3656                    max(data(varn,:,min_z_int:max_z_int)) .EQ. 0)then
3657                  min_value = 0.1
3658                  max_value = 0.1
3659                end if
3660
3661                res@gsnLeftString      = vNam(varn)
3662                res@tiXAxisString      = "("+unit(varn)+")"
3663                res@gsnRightString     = " "
3664                res@trYMinF            = min_z
3665                res@trYMaxF            = max_z 
3666                res@xyDashPattern = 0
3667
3668                if (xs .EQ. -1) then
3669                  res@trXMinF = min(data(varn,:,min_z_int:max_z_int))-min_value
3670                else
3671                  res@trXMinF = xs
3672                end if
3673                if (xe .EQ. -1) then
3674                  res@trXMaxF = max(data(varn,:,min_z_int:max_z_int))+max_value
3675                else
3676                  res@trXMaxF = xe 
3677                end if
3678                plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res) 
3679               
3680                if (vNam(varn) .EQ. "u" .OR. vNam(varn) .EQ. "v" .OR. \
3681                    vNam(varn) .EQ. "w") then
3682                  if (u .EQ. 0) then
3683                      res@gsnLeftString      = "u, v and w"
3684                      res@tiXAxisString      = "("+unit(varn)+")"
3685                      res@gsnRightString     = " "
3686                      if (xs .EQ. -1) then
3687                        res@trXMinF = min((/miniu,miniv,miniw/))
3688                      else
3689                        res@trXMinF = xs
3690                      end if
3691                      if (xe .EQ. -1) then
3692                        res@trXMaxF = max((/maxiu,maxiv,maxiw/))
3693                      else
3694                        res@trXMaxF = xe 
3695                      end if 
3696                      if (vNam(varn) .EQ. "v")then
3697                        res@xyDashPattern = 1
3698                      end if
3699                      if (vNam(varn) .EQ. "w")then
3700                        res@xyDashPattern = 2
3701                      end if
3702                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3703
3704                      ; ***************************************************
3705                      ; legend for overlaid plot
3706                      ; ***************************************************
3707         
3708                      lgres                    = True
3709                      lgMonoDashIndex          = False
3710                      lgres@lgLabelFont        = "helvetica"   
3711                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
3712                      lgres@vpWidthF           = 0.07           
3713                      lgres@vpHeightF          = 0.12         
3714                      lgres@lgDashIndexes      = (/0,1,2/)
3715                      lbid = gsn_create_legend(wks,3,(/"u","v","w"/),lgres)       
3716
3717                      amres = True
3718                      amres@amParallelPosF   = 0.88                 
3719                      amres@amOrthogonalPosF = 0.33           
3720                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3721                      overlay(plot(n),plot_u)
3722                      u=1
3723                  else
3724                      if (prof3d .EQ. 0)then
3725                        varn=varn+1
3726                      end if
3727                      continue
3728                  end if       
3729                end if 
3730       
3731                if (vNam(varn) .EQ. "pt" .OR. vNam(varn) .EQ. "vpt" .OR. \
3732                    vNam(varn) .EQ. "lpt") then
3733                  if (pt .EQ. 0) then
3734                      res@gsnLeftString      = "pt, vpt and lpt"
3735                      res@tiXAxisString      = "("+unit(varn)+")"
3736                      res@gsnRightString     = " "
3737                      if (xs .EQ. -1) then
3738                        res@trXMinF = min((/minipt,minivpt,minilpt/))
3739                      else
3740                        res@trXMinF = xs
3741                      end if
3742                      if (xe .EQ. -1) then
3743                        res@trXMaxF = max((/maxipt,maxivpt,maxilpt/))
3744                      else
3745                        res@trXMaxF = xe 
3746                      end if
3747                      if (vNam(varn) .EQ. "vpt")then
3748                        res@xyDashPattern = 1
3749                      end if
3750                      if (vNam(varn) .EQ. "lpt")then
3751                        res@xyDashPattern = 2
3752                      end if
3753                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3754
3755                      ; ***************************************************
3756                      ; legend for overlaid plot
3757                      ; ***************************************************
3758         
3759                      lgres                    = True
3760                      lgMonoDashIndex          = False
3761                      lgres@lgLabelFont        = "helvetica"   
3762                      lgres@lgLabelFontHeightF = font_size_legend*10.0         
3763                      lgres@vpWidthF           = 0.07           
3764                      lgres@vpHeightF          = 0.12         
3765                      lgres@lgDashIndexes      = (/0,1,2/)
3766                      lbid = gsn_create_legend(wks,3,(/"pt","vpt","lpt"/),lgres)
3767                      amres = True
3768                      amres@amParallelPosF   = 0.88         
3769                      amres@amOrthogonalPosF = 0.33           
3770                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3771                      overlay(plot(n),plot_pt)
3772                      pt=1
3773                  else
3774                      if (prof3d .EQ. 0)then
3775                        varn=varn+1
3776                      end if
3777                      continue       
3778                  end if
3779                end if           
3780                if (vNam(varn) .EQ. "q" .OR. vNam(varn) .EQ. "qv" .OR. \
3781                    vNam(varn) .EQ. "ql") then
3782                  if (q .EQ. 0) then
3783                      res@gsnLeftString      = "q, qv and ql"
3784                      res@tiXAxisString      = "("+unit(varn)+")"
3785                      res@gsnRightString     = " "
3786                      if (xs .EQ. -1) then
3787                        res@trXMinF = min((/miniq,miniqv,miniql/))
3788                      else
3789                        res@trXMinF = xs
3790                      end if
3791                      if (xe .EQ. -1) then
3792                        res@trXMaxF = max((/maxiq,maxiqv,maxiql/))
3793                      else
3794                        res@trXMaxF = xe 
3795                      end if
3796                      if (vNam(varn) .EQ. "qv")then
3797                        res@xyDashPattern = 1
3798                      end if
3799                      if (vNam(varn) .EQ. "ql")then
3800                        res@xyDashPattern = 2
3801                      end if
3802                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3803
3804                      ; ***************************************************
3805                      ; legend for overlaid plot
3806                      ; ***************************************************
3807         
3808                      lgres                    = True
3809                      lgMonoDashIndex          = False
3810                      lgres@lgLabelFont        = "helvetica"   
3811                      lgres@lgLabelFontHeightF = font_size_legend*10.0         
3812                      lgres@vpWidthF           = 0.07         
3813                      lgres@vpHeightF          = 0.12         
3814                      lgres@lgDashIndexes      = (/0,1,2/)
3815                      lbid = gsn_create_legend(wks,3,(/"q","qv","ql"/),lgres)
3816
3817                      amres = True
3818                      amres@amParallelPosF   = 0.88                 
3819                      amres@amOrthogonalPosF = 0.33           
3820                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3821                      overlay(plot(n),plot_q)
3822                      q=1
3823                  else
3824                      if (prof3d .EQ. 0)then
3825                        varn=varn+1
3826                      end if
3827                      continue   
3828                  end if
3829                end if   
3830             
3831                if (vNam(varn) .EQ. "e" .OR. vNam(varn) .EQ. "es" .OR. \
3832                    vNam(varn) .EQ. "e*" ) then
3833                  if (e .EQ. 0) then
3834                      res@gsnLeftString      = "e and e*"
3835                      res@tiXAxisString      = "("+unit(varn)+")"
3836                      res@gsnRightString     = " "
3837                      if (xs .EQ. -1) then
3838                        res@trXMinF = min((/minie,minies/))
3839                      else
3840                        res@trXMinF = xs
3841                      end if
3842                      if (xe .EQ. -1) then
3843                        res@trXMaxF = max((/maxie,maxies/))
3844                      else
3845                        res@trXMaxF = xe 
3846                      end if
3847                      if (vNam(varn) .EQ. "es")then
3848                        res@xyDashPattern = 1
3849                      end if
3850                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3851
3852                      ; ***************************************************
3853                      ; legend for overlaid plot
3854                      ; ***************************************************
3855         
3856                      lgres                    = True
3857                      lgMonoDashIndex          = False
3858                      lgres@lgLabelFont        = "helvetica"   
3859                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
3860                      lgres@vpWidthF           = 0.07           
3861                      lgres@vpHeightF          = 0.08         
3862                      lgres@lgDashIndexes      = (/0,1,2/)
3863                      lbid = gsn_create_legend(wks,2,(/"e","e*"/),lgres)       
3864
3865                      amres = True
3866                      amres@amParallelPosF   = 0.88                 
3867                      amres@amOrthogonalPosF = 0.365           
3868                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3869                      overlay(plot(n),plot_e)
3870                      e=1
3871                  else
3872                      if (prof3d .EQ. 0)then
3873                        varn=varn+1
3874                      end if
3875                      continue   
3876                  end if
3877                end if           
3878                if (vNam(varn) .EQ. "km" .OR. vNam(varn) .EQ. "kh") then
3879                  if (km .EQ. 0) then
3880                      res@gsnLeftString      = "km and kh"
3881                      res@tiXAxisString      = "("+unit(varn)+")"
3882                      res@gsnRightString     = " "
3883                      if (xs .EQ. -1) then
3884                        res@trXMinF = min((/minikm,minikh/))
3885                      else
3886                        res@trXMinF = xs
3887                      end if
3888                      if (xe .EQ. -1) then
3889                        res@trXMaxF = max((/maxikm,maxikh/))
3890                      else
3891                        res@trXMaxF = xe 
3892                      end if
3893                      if (vNam(varn) .EQ. "kh")then
3894                        res@xyDashPattern = 1
3895                      end if
3896                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3897
3898                      ; ***************************************************
3899                      ; legend for overlaid plot
3900                      ; ***************************************************
3901         
3902                      lgres                    = True
3903                      lgMonoDashIndex          = False
3904                      lgres@lgLabelFont        = "helvetica"   
3905                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
3906                      lgres@vpWidthF           = 0.07           
3907                      lgres@vpHeightF          = 0.08         
3908                      lgres@lgDashIndexes      = (/0,1,2/)
3909                      lbid = gsn_create_legend(wks,2,(/"km","kh"/),lgres)       
3910
3911                      amres = True
3912                      amres@amParallelPosF   = 0.88                 
3913                      amres@amOrthogonalPosF = 0.365           
3914                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3915                      overlay(plot(n),plot_km)
3916                      km=1
3917                  else
3918                      if (prof3d .EQ. 0)then
3919                        varn=varn+1
3920                      end if
3921                      continue   
3922                  end if
3923                end if           
3924             
3925                if (vNam(varn) .EQ. "wpup" .OR. vNam(varn) .EQ. "wsus" .OR.      \
3926                    vNam(varn) .EQ. "wu" .OR. vNam(varn) .EQ. "w"+dq+"u"+dq .OR. \
3927                    vNam(varn) .EQ. "w*u*") then
3928                  if (wpup .EQ. 0) then
3929                      res@gsnLeftString      = "w"+dq+"u"+dq+", w*u* and wu"
3930                      res@tiXAxisString      = "("+unit(varn)+")"
3931                      res@gsnRightString     = " "
3932                      if (xs .EQ. -1) then
3933                        res@trXMinF = min((/miniwpup,miniwsus,miniwu/))
3934                      else
3935                        res@trXMinF = xs
3936                      end if
3937                      if (xe .EQ. -1) then
3938                        res@trXMaxF = max((/maxiwpup,maxiwsus,maxiwu/))
3939                      else
3940                        res@trXMaxF = xe 
3941                      end if 
3942                      if (vNam(varn) .EQ. "wsus")then
3943                        res@xyDashPattern = 1
3944                      end if
3945                      if (vNam(varn) .EQ. "wu")then
3946                        res@xyDashPattern = 2
3947                      end if
3948                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
3949
3950                      ; ***************************************************
3951                      ; legend for overlaid plot
3952                      ; ***************************************************
3953         
3954                      lgres                    = True
3955                      lgMonoDashIndex          = False
3956                      lgres@lgLabelFont        = "helvetica"   
3957                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
3958                      lgres@vpWidthF           = 0.08           
3959                      lgres@vpHeightF          = 0.12         
3960                      lgres@lgDashIndexes      = (/0,1,2/)
3961                      lbid = gsn_create_legend(\
3962                                    wks,3,(/"w"+dq+"u"+dq,"w*u*","wu"/),lgres)
3963
3964                      amres = True
3965                      amres@amParallelPosF   = 0.88                 
3966                      amres@amOrthogonalPosF = 0.33           
3967                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
3968                      overlay(plot(n),plot_wpup)
3969                      wpup=1
3970                  else
3971                      if (prof3d .EQ. 0)then
3972                        varn=varn+1
3973                      end if
3974                      continue   
3975                  end if
3976                end if
3977                if (vNam(varn) .EQ. "wpvp" .OR. vNam(varn) .EQ. "wsvs" .OR.     \
3978                  vNam(varn) .EQ. "wv" .OR. vNam(varn) .EQ. "w"+dq+"v"+dq .OR. \
3979                  vNam(varn) .EQ. "w*v*") then
3980                  if (wpvp .EQ. 0) then
3981                      res@gsnLeftString      = "w"+dq+"v"+dq+", w*v* and wv"
3982                      res@tiXAxisString      = "("+unit(varn)+")"
3983                      res@gsnRightString     = " "
3984                      if (xs .EQ. -1) then
3985                        res@trXMinF = min((/miniwpvp,miniwsvs,miniwv/))
3986                      else
3987                        res@trXMinF = xs
3988                      end if
3989                      if (xe .EQ. -1) then
3990                        res@trXMaxF = max((/maxiwpvp,maxiwsvs,maxiwv/))
3991                      else
3992                        res@trXMaxF = xe 
3993                      end if
3994                      if (vNam(varn) .EQ. "wsvs")then
3995                        res@xyDashPattern = 1
3996                      end if
3997                      if (vNam(varn) .EQ. "wv")then
3998                        res@xyDashPattern = 2
3999                      end if
4000                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4001
4002                      ; ***************************************************
4003                      ; legend for overlaid plot
4004                      ; ***************************************************
4005         
4006                      lgres                    = True
4007                      lgMonoDashIndex          = False
4008                      lgres@lgLabelFont        = "helvetica"   
4009                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4010                      lgres@vpWidthF           = 0.08         
4011                      lgres@vpHeightF          = 0.12         
4012                      lgres@lgDashIndexes      = (/0,1,2/)
4013                      lbid = gsn_create_legend(\
4014                                wks,3,(/"w"+dq+"v"+dq,"w*v*","wv"/),lgres)       
4015
4016                      amres = True
4017                      amres@amParallelPosF   = 0.88                 
4018                      amres@amOrthogonalPosF = 0.33           
4019                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4020                      overlay(plot(n),plot_wpvp)
4021                      wpvp=1
4022                  else
4023                      if (prof3d .EQ. 0)then
4024                        varn=varn+1
4025                      end if
4026                      continue   
4027                  end if
4028                end if
4029                if (vNam(varn) .EQ. "wpptp" .OR. vNam(varn) .EQ. "wspts" .OR.     \
4030                    vNam(varn) .EQ. "wpt" .OR. vNam(varn) .EQ. "w"+dq+"pt"+dq .OR.\
4031                    vNam(varn) .EQ. "w*pt*") then
4032                  if (wpptp .EQ. 0) then                 
4033                      res@gsnLeftString      = "w"+dq+"pt"+dq+", w*pt* and wpt"
4034                      res@tiXAxisString      = "("+unit(varn)+")"
4035                      res@gsnRightString     = " "
4036                      if (xs .EQ. -1) then
4037                        res@trXMinF = min((/miniwpptp,miniwspts,miniwpt/))
4038                      else
4039                        res@trXMinF = xs
4040                      end if
4041                      if (xe .EQ. -1) then
4042                        res@trXMaxF = max((/maxiwpptp,maxiwspts,maxiwpt/))
4043                      else
4044                        res@trXMaxF = xe 
4045                      end if
4046                      if (vNam(varn) .EQ. "wspts" .OR. vNam(varn) .EQ. "w*pt*")then
4047                        res@xyDashPattern = 1
4048                      end if
4049                      if (vNam(varn) .EQ. "wpt")then
4050                        res@xyDashPattern = 2
4051                      end if
4052                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4053
4054                      ; ***************************************************
4055                      ; legend for overlaid plot
4056                      ; ***************************************************
4057         
4058                      lgres                    = True
4059                      lgMonoDashIndex          = False
4060                      lgres@lgLabelFont        = "helvetica"   
4061                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4062                      lgres@vpWidthF           = 0.09           
4063                      lgres@vpHeightF          = 0.12         
4064                      lgres@lgDashIndexes      = (/0,1,2/)             
4065                      lbid = gsn_create_legend(\
4066                                  wks,3,(/"w"+dq+"pt"+dq,"w*pt*","wpt"/),lgres)
4067
4068                      amres = True
4069                      amres@amParallelPosF   = 0.88                 
4070                      amres@amOrthogonalPosF = 0.33           
4071                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4072                      overlay(plot(n),plot_wpptp)
4073                      wpptp=1
4074                  else
4075                      if (prof3d .EQ. 0)then
4076                        varn=varn+1
4077                      end if
4078                      continue   
4079                  end if
4080                end if
4081                if (vNam(varn) .EQ. "wsptsBC" .OR. vNam(varn) .EQ. "wptBC" .OR.\
4082                    vNam(varn) .EQ. "w*pt*BC") then
4083                  if (wsptsBC .EQ. 0) then
4084                      res@gsnLeftString      = "w*pt*BC and wptBC"
4085                      res@tiXAxisString      = "("+unit(varn)+")"
4086                      res@gsnRightString     = " "
4087                      if (xs .EQ. -1) then
4088                        res@trXMinF = min((/miniwsptsBC,miniwptBC/))
4089                      else
4090                        res@trXMinF = xs
4091                      end if
4092                      if (xe .EQ. -1) then
4093                        res@trXMaxF = max((/maxiwsptsBC,maxiwptBC/))
4094                      else
4095                        res@trXMaxF = xe 
4096                      end if
4097                      if (vNam(varn) .EQ. "wptBC")then
4098                        res@xyDashPattern = 1
4099                      end if
4100                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4101
4102                      ; ***************************************************
4103                      ; legend for overlaid plot
4104                      ; ***************************************************
4105         
4106                      lgres                    = True
4107                      lgMonoDashIndex          = False
4108                      lgres@lgLabelFont        = "helvetica"   
4109                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4110                      lgres@vpWidthF           = 0.1           
4111                      lgres@vpHeightF          = 0.12         
4112                      lgres@lgDashIndexes      = (/0,1,2/)
4113                      lbid = gsn_create_legend(\
4114                                            wks,3,(/"w*pt*BC","wptBC"/),lgres)
4115
4116                      amres = True
4117                      amres@amParallelPosF   = 0.88                 
4118                      amres@amOrthogonalPosF = 0.33           
4119                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4120                      overlay(plot(n),plot_wsptsBC)
4121                      wsptsBC=1
4122                  else
4123                      if (prof3d .EQ. 0)then
4124                        varn=varn+1
4125                      end if
4126                      continue   
4127                  end if 
4128                end if             
4129                if (vNam(varn) .EQ. "wpvptp" .OR. vNam(varn) .EQ. "wsvpts" .OR. \
4130                    vNam(varn) .EQ. "wvpt" .OR. vNam(varn) .EQ. \
4131                    "w"+dq+"vpt"+dq .OR. vNam(varn) .EQ. "w*vpt*") then
4132                  if (wpvptp .EQ. 0) then
4133                      res@gsnLeftString      = "w"+dq+"vpt"+dq+", w*vpt* and wvpt"
4134                      res@tiXAxisString      = "("+unit(varn)+")"
4135                      res@gsnRightString     = " "
4136                      if (xs .EQ. -1) then
4137                        res@trXMinF = min((/miniwpvptp,miniwsvpts,miniwvpt/))
4138                      else
4139                        res@trXMinF = xs
4140                      end if
4141                      if (xe .EQ. -1) then
4142                        res@trXMaxF = max((/maxiwpvptp,maxiwsvpts,maxiwvpt/))
4143                      else
4144                        res@trXMaxF = xe 
4145                      end if
4146                      if (vNam(varn) .EQ. "wsvpts")then
4147                        res@xyDashPattern = 1
4148                      end if
4149                      if (vNam(varn) .EQ. "wvpt")then
4150                        res@xyDashPattern = 2
4151                      end if
4152                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4153
4154                      ; ***************************************************
4155                      ; legend for overlaid plot
4156                      ; ***************************************************
4157         
4158                      lgres                    = True
4159                      lgMonoDashIndex          = False
4160                      lgres@lgLabelFont        = "helvetica"   
4161                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4162                      lgres@vpWidthF           = 0.1           
4163                      lgres@vpHeightF          = 0.12         
4164                      lgres@lgDashIndexes      = (/0,1,2/)
4165                      lbid = gsn_create_legend(\
4166                                wks,3,(/"w"+dq+"vpt"+dq,"w*vpt*","wvpt"/),lgres)
4167                      amres = True
4168                      amres@amParallelPosF   = 0.88                 
4169                      amres@amOrthogonalPosF = 0.33           
4170                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4171                      overlay(plot(n),plot_wpvptp)
4172                      wpvptp=1
4173                  else
4174                      if (prof3d .EQ. 0)then
4175                        varn=varn+1
4176                      end if
4177                      continue   
4178                  end if
4179                end if
4180                if (vNam(varn) .EQ. "wpqp" .OR. vNam(varn) .EQ. "wsqs" .OR.      \
4181                    vNam(varn) .EQ. "wq" .OR. vNam(varn) .EQ. "w"+dq+"q"+dq .OR. \
4182                    vNam(varn) .EQ. "w*q*") then
4183                  if (wpqp .EQ. 0) then
4184                      res@gsnLeftString      = "w"+dq+"q"+dq+", w*q* and wq"
4185                      res@tiXAxisString      = "("+unit(varn)+")"
4186                      res@gsnRightString     = " "
4187                      if (xs .EQ. -1) then
4188                        res@trXMinF = min((/miniwpqp,miniwsqs,miniwq/))
4189                      else
4190                        res@trXMinF = xs
4191                      end if
4192                      if (xe .EQ. -1) then
4193                        res@trXMaxF = max((/maxiwpqp,maxiwsqs,maxiwq/))
4194                      else
4195                        res@trXMaxF = xe 
4196                      end if 
4197                      if (vNam(varn) .EQ. "wsqs")then
4198                        res@xyDashPattern = 1
4199                      end if
4200                      if (vNam(varn) .EQ. "wq")then
4201                        res@xyDashPattern = 2
4202                      end if
4203                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4204
4205                      ; ***************************************************
4206                      ; legend for overlaid plot
4207                      ; ***************************************************
4208         
4209                      lgres                    = True
4210                      lgMonoDashIndex          = False
4211                      lgres@lgLabelFont        = "helvetica"   
4212                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4213                      lgres@vpWidthF           = 0.08           
4214                      lgres@vpHeightF          = 0.12         
4215                      lgres@lgDashIndexes      = (/0,1,2/)
4216                      lbid = gsn_create_legend(\
4217                                wks,3,(/"w"+dq+"q"+dq,"w*q*","wq"/),lgres)       
4218
4219                      amres = True
4220                      amres@amParallelPosF   = 0.88                 
4221                      amres@amOrthogonalPosF = 0.33           
4222                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4223                      overlay(plot(n),plot_wpqp)
4224                      wpqp=1
4225                  else
4226                      if (prof3d .EQ. 0)then
4227                        varn=varn+1
4228                      end if
4229                      continue   
4230                  end if
4231                end if
4232                if (vNam(varn) .EQ. "wpqvp" .OR. vNam(varn) .EQ. "wsqvs" .OR.     \
4233                    vNam(varn) .EQ. "wqv" .OR. vNam(varn) .EQ. "w"+dq+"qv"+dq .OR.\
4234                    vNam(varn) .EQ. "w*qv*") then
4235                  if (wpqvp .EQ. 0) then
4236                      res@gsnLeftString      ="w"+dq+"qv"+dq+" , w*qv* and wqv"
4237                      res@tiXAxisString      = "("+unit(varn)+")"
4238                      res@gsnRightString     = " "
4239                      if (xs .EQ. -1) then
4240                        res@trXMinF = min((/miniwpqp,miniwsqvs,miniwqv/))
4241                      else
4242                        res@trXMinF = xs
4243                      end if
4244                      if (xe .EQ. -1) then
4245                        res@trXMaxF = max((/maxiwpqp,maxiwsqvs,maxiwqv/))
4246                      else
4247                        res@trXMaxF = xe 
4248                      end if
4249                      if (vNam(varn) .EQ. "wsqvs")then
4250                        res@xyDashPattern = 1
4251                      end if
4252                      if (vNam(varn) .EQ. "wqv")then
4253                        res@xyDashPattern = 2
4254                      end if
4255                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4256
4257                      ; ***************************************************
4258                      ; legend for overlaid plot
4259                      ; ***************************************************
4260         
4261                      lgres                    = True
4262                      lgMonoDashIndex          = False
4263                      lgres@lgLabelFont        = "helvetica"   
4264                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4265                      lgres@vpWidthF           = 0.09           
4266                      lgres@vpHeightF          = 0.12         
4267                      lgres@lgDashIndexes      = (/0,1,2/)
4268                      lbid = gsn_create_legend(\
4269                                    wks,3,(/"w"+dq+"qv"+dq,"w*qv*","wqv"/),lgres)
4270
4271                      amres = True
4272                      amres@amParallelPosF   = 0.88                 
4273                      amres@amOrthogonalPosF = 0.33           
4274                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4275                      overlay(plot(n),plot_wpqvp)
4276                      wpqvp=1
4277                  else
4278                      if (prof3d .EQ. 0)then
4279                        varn=varn+1
4280                      end if
4281                      continue   
4282                  end if
4283                end if
4284                if (vNam(varn) .EQ. "wpsp" .OR. vNam(varn) .EQ. "wsss" .OR.     \
4285                    vNam(varn) .EQ. "ws" .OR. vNam(varn) .EQ. "w"+dq+"s"+dq .OR.\
4286                    vNam(varn) .EQ. "w*s*") then
4287                  if (wpsp .EQ. 0) then
4288                      res@gsnLeftString      = "w"+dq+"s"+dq+", w*s* and ws"
4289                      res@tiXAxisString      = "("+unit(varn)+")"
4290                      res@gsnRightString     = " "
4291                      if (xs .EQ. -1) then
4292                        res@trXMinF = min((/miniwpsp,miniwsss,miniws/))
4293                      else
4294                        res@trXMinF = xs
4295                      end if
4296                      if (xe .EQ. -1) then
4297                        res@trXMaxF = max((/maxiwpsp,maxiwsss,maxiws/))
4298                      else
4299                        res@trXMaxF = xe 
4300                      end if
4301                      if (vNam(varn) .EQ. "wsss")then
4302                        res@xyDashPattern = 1
4303                      end if
4304                      if (vNam(varn) .EQ. "ws")then
4305                        res@xyDashPattern = 2
4306                      end if
4307                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4308
4309                      ; ***************************************************
4310                      ; legend for overlaid plot
4311                      ; ***************************************************
4312         
4313                      lgres                    = True
4314                      lgMonoDashIndex          = False
4315                      lgres@lgLabelFont        = "helvetica"   
4316                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4317                      lgres@vpWidthF           = 0.08           
4318                      lgres@vpHeightF          = 0.12         
4319                      lgres@lgDashIndexes      = (/0,1,2/)
4320                      lbid = gsn_create_legend(\
4321                              wks,3,(/"w"+dq+"s"+dq,"w*s*","ws"/),lgres)       
4322
4323                      amres = True
4324                      amres@amParallelPosF   = 0.88                 
4325                      amres@amOrthogonalPosF = 0.33           
4326                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4327                      overlay(plot(n),plot_wpsp)
4328                      wpsp=1
4329                  else
4330                      if (prof3d .EQ. 0)then
4331                        varn=varn+1
4332                      end if
4333                      continue   
4334                  end if
4335                end if
4336                if (vNam(varn) .EQ. "wpsap" .OR.vNam(varn) .EQ. "wssas" .OR.      \
4337                    vNam(varn) .EQ. "wsa" .OR. vNam(varn) .EQ. "w"+dq+"sa"+dq .OR.\
4338                    vNam(varn) .EQ. "w*sa*") then
4339                  if (wpsap .EQ. 0) then
4340                      res@gsnLeftString      = "w"+dq+"sa"+dq+", w*sa* and wsa"
4341                      res@tiXAxisString      = "("+unit(varn)+")"
4342                      res@gsnRightString     = " "
4343                      if (xs .EQ. -1) then
4344                        res@trXMinF = min((/miniwpsap,miniwssas,miniwsa/))
4345                      else
4346                        res@trXMinF = xs
4347                      end if
4348                      if (xe .EQ. -1) then
4349                        res@trXMaxF = max((/maxiwpsap,maxiwssas,maxiwsa/))
4350                      else
4351                        res@trXMaxF = xe 
4352                      end if
4353                      if (vNam(varn) .EQ. "wssas")then
4354                        res@xyDashPattern = 1
4355                      end if
4356                      if (vNam(varn) .EQ. "wsa")then
4357                        res@xyDashPattern = 2
4358                      end if
4359                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4360
4361                      ; ***************************************************
4362                      ; legend for overlaid plot
4363                      ; ***************************************************
4364         
4365                      lgres                    = True
4366                      lgMonoDashIndex          = False
4367                      lgres@lgLabelFont        = "helvetica"   
4368                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4369                      lgres@vpWidthF           = 0.09           
4370                      lgres@vpHeightF          = 0.12         
4371                      lgres@lgDashIndexes      = (/0,1,2/)
4372                      lbid = gsn_create_legend(\
4373                                wks,3,(/"w"+dq+"sa"+dq,"w*sa*","wsa"/),lgres)
4374                      amres = True
4375                      amres@amParallelPosF   = 0.88                 
4376                      amres@amOrthogonalPosF = 0.33           
4377                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4378                      overlay(plot(n),plot_wpsap)
4379                      wpsap=1
4380                  else
4381                      if (prof3d .EQ. 0)then
4382                        varn=varn+1
4383                      end if
4384                      continue   
4385                  end if
4386                end if
4387             
4388                if (vNam(varn) .EQ. "us2" .OR. vNam(varn) .EQ. "vs2" .OR. \
4389                    vNam(varn) .EQ. "ws2" .OR. vNam(varn) .EQ. "u*2" .OR. \
4390                    vNam(varn) .EQ. "v*2" .OR. vNam(varn) .EQ. "w*2" ) then
4391                  if (us2 .EQ. 0) then
4392                      res@gsnLeftString      = "u*2, v*2 and w*2"
4393                      res@tiXAxisString      = "("+unit(varn)+")"
4394                      res@gsnRightString     = " "
4395                      if (xs .EQ. -1) then
4396                        res@trXMinF = min((/minius2,minivs2,miniws2/))
4397                      else
4398                        res@trXMinF = xs
4399                      end if
4400                      if (xe .EQ. -1) then
4401                        res@trXMaxF = max((/maxius2,maxivs2,maxiws2/))
4402                      else
4403                        res@trXMaxF = xe 
4404                      end if
4405                      if (vNam(varn) .EQ. "vs2")then
4406                        res@xyDashPattern = 1
4407                      end if
4408                      if (vNam(varn) .EQ. "ws2")then
4409                        res@xyDashPattern = 2
4410                      end if
4411                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4412
4413                      ; ***************************************************
4414                      ; legend for overlaid plot
4415                      ; ***************************************************
4416         
4417                      lgres                    = True
4418                      lgMonoDashIndex          = False
4419                      lgres@lgLabelFont        = "helvetica"   
4420                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4421                      lgres@vpWidthF           = 0.07           
4422                      lgres@vpHeightF          = 0.12         
4423                      lgres@lgDashIndexes      = (/0,1,2/)
4424                      lbid = gsn_create_legend(wks,3,(/"u*2","v*2","w*2"/),lgres)
4425                      amres = True
4426                      amres@amParallelPosF   = 0.88                 
4427                      amres@amOrthogonalPosF = 0.33           
4428                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4429                      overlay(plot(n),plot_us2)
4430                      us2=1
4431                  else
4432                      if (prof3d .EQ. 0)then
4433                        varn=varn+1
4434                      end if
4435                      continue   
4436                  end if
4437                end if
4438             
4439                if (vNam(varn) .EQ. "wsususodz" .OR. \
4440                    vNam(varn) .EQ. "wspsodz" .OR.   \
4441                    vNam(varn) .EQ. "wpeodz" .OR.    \
4442                    vNam(varn) .EQ. "w*u*u*:dz" .OR. \
4443                    vNam(varn) .EQ. "w*p*:dz" .OR.   \
4444                    vNam(varn) .EQ. "w"+dq+"e:dz") then
4445                  if (wsususodz .EQ. 0) then
4446                      res@gsnLeftString      = "w*u*u*:dz, w*p*:dz and w"+dq+"e:dz"
4447                      res@tiXAxisString      = "("+unit(varn)+")"
4448                      res@gsnRightString     = " "
4449                      if (xs .EQ. -1) then
4450                        res@trXMinF = min((/miniwsususodz,\
4451                                                      miniwspsodz,miniwpeodz/))
4452                      else
4453                        res@trXMinF = xs
4454                      end if
4455                      if (xe .EQ. -1) then
4456                        res@trXMaxF = max((/maxiwsususodz,maxiwspsodz,\
4457                                                                  maxiwpeodz/))
4458                      else
4459                        res@trXMaxF = xe 
4460                      end if
4461                      if (vNam(varn) .EQ. "wspsodz")then
4462                        res@xyDashPattern = 1
4463                      end if
4464                      if (vNam(varn) .EQ. "wpeodz")then
4465                        res@xyDashPattern = 2
4466                      end if
4467                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z,res)
4468
4469                      ; ***************************************************
4470                      ; legend for overlaid plot
4471                      ; ***************************************************
4472         
4473                      lgres                    = True
4474                      lgMonoDashIndex          = False
4475                      lgres@lgLabelFont        = "helvetica"   
4476                      lgres@lgLabelFontHeightF = font_size_legend*10.0           
4477                      lgres@vpWidthF           = 0.12           
4478                      lgres@vpHeightF          = 0.12         
4479                      lgres@lgDashIndexes      = (/0,1,2/)
4480                      lbid = gsn_create_legend(\
4481                              wks,3,(/"w*u*u*:dz","w*p*:dz","w"+dq+"e:dz"/),lgres)
4482                      amres = True
4483                      amres@amParallelPosF   = 0.88                 
4484                      amres@amOrthogonalPosF = 0.33           
4485                      annoid1 = gsn_add_annotation(plot(n),lbid,amres)
4486                      overlay(plot(n),plot_wsususodz)
4487                      wsususodz=1
4488                  else
4489                      if (prof3d .EQ. 0)then
4490                        varn=varn+1
4491                      end if
4492                      continue   
4493                  end if
4494                end if     
4495                n=n+1
4496                if (prof3d .EQ. 0)then
4497                  varn=varn+1
4498                end if   
4499            end if
4500          end do
4501      end if
4502     
4503      if (com_i .NE. 0) then
4504          delete(com_var_avail)
4505      end if
4506      com_var_avail=new(count_var,string)   
4507
4508      if (combine .EQ. 1) then
4509          co=0     
4510          n_o=0
4511          do varn = 0,dim-1
4512   
4513            check = True
4514         
4515            if (prof3d .EQ. 0) then
4516                if ( isStrSubset( vNam(varn), "time") .OR. \
4517                    isStrSubset( vNam(varn), "NORM")) then
4518                  check = False
4519                end if
4520            else
4521                if ( isStrSubset( vNam(varn), "time") .OR.  \
4522                    isStrSubset( vNam(varn), "zusi") .OR.  \
4523                    isStrSubset( vNam(varn), "zwwi") .OR.  \
4524                    isStrSubset( vNam(varn), "x") .OR.     \
4525                    isStrSubset( vNam(varn), "xu") .OR.    \
4526                    isStrSubset( vNam(varn), "y") .OR.     \
4527                    isStrSubset( vNam(varn), "yv") .OR.    \
4528                    isStrSubset( vNam(varn), "zu_3d") .OR. \
4529                    isStrSubset( vNam(varn), "zw_3d")) then
4530                  check = False
4531                end if
4532            end if
4533
4534            if (var .NE. "all") then         
4535                check = isStrSubset( var,","+vNam(varn)+"," )
4536            end if     
4537
4538            if (check)then
4539               
4540                if (prof3d .EQ. 0) then
4541                  if (log_z .EQ. 1) then
4542                      z = f_att->$vNam(varn+1)$(1:dimz-1)
4543                  else
4544                      z = f_att->$vNam(varn+1)$               
4545                  end if
4546                else
4547                  do i=0,b-1           
4548                      if (isStrSubset( a(i),"zu_3d" ))then
4549                        z_v(varn,:) = z_u
4550                        if (log_z .EQ. 1) then
4551                            z = z_v(varn,1:dimz-1)
4552                        else
4553                            z = z_v(varn,:)
4554                        end if
4555                      else
4556                        if (isStrSubset( a(i),"zw_3d" ))then
4557                            z_v(varn,:) = z_w
4558                            if (log_z .EQ. 1) then
4559                              z = z_v(varn,1:dimz-1)
4560                            else
4561                              z = z_v(varn,:)
4562                            end if
4563                        end if                   
4564                      end if
4565                  end do           
4566                end if
4567
4568                z=z/norm_z
4569               
4570                com_var_avail(n_o)=vNam(varn)
4571               
4572                com=isStrSubset( c_var,","+vNam(varn)+"," )
4573                if (com)then
4574                  co = co+1           
4575                  if (n_o .EQ. 1) then
4576                      res@xyDashPattern  = 1                                   
4577                  else           
4578                      if (n_o .EQ. 2) then
4579                        res@xyDashPattern  = 2
4580                      else
4581                        res@xyDashPattern  = 0
4582                        res@gsnLeftString  = " " ; "Combined Plot of "+c_var_name
4583                        res@tiXAxisString  = "("+unit(varn)+")"
4584                        res@gsnRightString = " "
4585                        if (xs .EQ. -1) then
4586                            res@trXMinF = min(mini)
4587                        else
4588                            res@trXMinF = xs
4589                        end if
4590                        if (xe .EQ. -1) then
4591                            res@trXMaxF = max(maxi)
4592                        else
4593                            res@trXMaxF = xe 
4594                        end if
4595                      end if
4596                  end if
4597                  label(n_o)=" " + vNam(varn)
4598                  color_o(n_o)=237
4599                  plot_o(n_o)=gsn_csm_xy(wks,data(varn,:,:),z,res)
4600                  n_o=n_o+1
4601                end if
4602                if (prof3d .EQ. 0)then
4603                  varn=varn+1
4604                end if           
4605            end if
4606          end do
4607     
4608          if(number_comb .EQ. 2)then
4609            if (co .EQ. 2)then
4610                overlay(plot_o(0),plot_o(1))
4611            else
4612                print(" ")
4613                print("combining is not possible,")
4614                print("'c_var'(= "+c_var+") must include two variables of "+\
4615                      "the general plots = ")
4616                print("- "+com_var_avail)
4617                print("be sure to have one comma before and after the variable")
4618                print(" ")
4619                exit 
4620            end if
4621          end if
4622          if(number_comb .EQ. 3)then
4623            if (co .EQ. 3)then
4624                overlay(plot_o(0),plot_o(1))
4625                overlay(plot_o(0),plot_o(2))
4626            else
4627                print(" ")
4628                print("combining is not possible,")
4629                print("'c_var'(= "+c_var+") must include three variables of "+\
4630                      "the general plots = ")
4631                print("- "+com_var_avail)
4632                print("be sure to have one comma before and after the variable")
4633                print(" ")
4634                exit
4635            end if
4636          end if
4637
4638          ; ***************************************************
4639          ; legend for combined plot
4640          ; ***************************************************
4641         
4642          lgres                    = True
4643          lgMonoDashIndex          = False
4644          lgres@lgLineDashSegLenF  = 0.075
4645          lgres@lgDashIndexes      = (/0,1,2/)
4646          lgres@lgLabelFont        = "helvetica"   
4647          lgres@vpWidthF           = 0.14           
4648          lgres@vpHeightF          = 0.1 * number_comb / 3.0 + (number_comb + 1) * 0.01       
4649          lgres@lgLabelFontHeightF = 10.0 * font_size_legend * ((number_comb \
4650                                     - 2.0) * 0.35 + 0.65)
4651
4652          lbid = gsn_create_legend(wks,number_comb,label,lgres)       
4653
4654          amres = True
4655          amres@amParallelPosF   = -0.5             
4656          amres@amOrthogonalPosF = -0.5   
4657          amres@amJust           = "TopLeft"
4658
4659          annoid1 = gsn_add_annotation(plot_o(0),lbid,amres)
4660       
4661          plot(0) = plot_o(0)
4662
4663      end if
4664
4665      if (c_var_log .EQ. 1) then
4666          c_var_plot(com_i) = plot(0)
4667      end if
4668
4669   end do
4670
4671   ; ***************************************************
4672   ; merge plots onto one page
4673   ; ***************************************************
4674
4675   if (c_var_log .EQ. 1) then
4676      delete(plot_)
4677      plot_=new(max_com_i+1,graphic)
4678      plot_=c_var_plot
4679      n = max_com_i + 1
4680   else
4681      do m=0,n-1
4682         plot_(m)=plot(n-1-m)
4683      end do
4684   end if
4685
4686   no_frames = 0
4687
4688   if ((format_out .EQ. "eps" .OR. format_out .EQ. "epsi") .AND. \
4689       n .gt. no_rows*no_columns) then
4690      gsn_panel(wks,plot_,(/n,1/),resP)
4691      print(" ")
4692      print("Outputs to .eps or .epsi have only one frame")
4693      print(" ")
4694   else   
4695      do i = 0,n-1, no_rows*no_columns
4696         if( (i+no_rows*no_columns) .gt. (n-1)) then
4697            gsn_panel(wks,plot_(i:n-1),(/no_rows,no_columns/),resP)
4698            no_frames = no_frames + 1 
4699         else
4700            gsn_panel(wks,plot_(i:i+no_rows*no_columns-1),\
4701                                                (/no_rows,no_columns/),resP)
4702            no_frames = no_frames + 1 
4703         end if
4704      end do
4705   end if
4706
4707
4708   if (format_out .EQ. "png" ) then
4709     png_output = new((/no_frames/), string)
4710     j = 0
4711
4712     if (no_frames .eq. 1) then
4713        if (ncl_version .GE. 6 ) then
4714           png_output(0) = file_out+".png"
4715        else
4716           png_output(0) = file_out+".00000"+1+".png"
4717        end if
4718        ;using imagemagick's convert for reducing the white
4719        ;space around the plot
4720        cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
4721              png_output(0) + " " + png_output(0)
4722       system(cmd)
4723     else
4724
4725       do i=0, no_frames-1
4726         j = i + 1
4727         if (j .LE. 9) then
4728           png_output(i) = file_out+".00000"+j+".png"
4729         end if
4730         if (j .GT. 9 .AND. j .LE. 99) then
4731           png_output(i) = file_out+".0000"+j+".png"
4732         end if
4733         if (j .GT. 99 .AND. j .LE. 999) then
4734           png_output(i) = file_out+".000"+j+".png"
4735         end if
4736         if (j .GT. 999) then
4737           png_output(i) = file_out+".00"+j+".png"
4738         end if
4739
4740         ;using imagemagick's convert for reducing the white
4741         ;space around the plot
4742         cmd = "convert -geometry 1000x1000 -density 300 -trim " +  \
4743                png_output(i) + " " + png_output(i)
4744         system(cmd)
4745       end do
4746     end if
4747
4748     print(" ")
4749     print("Output to: "+ png_output)
4750     print(" ")
4751   else
4752     print(" ")
4753     print("Output to: " + file_out +"."+ format_out)
4754     print(" ")
4755   end if
4756
4757end
Note: See TracBrowser for help on using the repository browser.