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

Last change on this file since 523 was 523, checked in by heinze, 14 years ago

Bugfix in profiles.ncl concerning no_files>1

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