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

Last change on this file since 992 was 963, checked in by hoffmann, 12 years ago

minor improvements in palmplot pr

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