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

Last change on this file since 1248 was 1248, checked in by heinze, 9 years ago

NCL function getfilevarnames changed with NCL version 6.1.1 and higher -> adaptation required

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