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

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

Enable plot of profiles when data is of kind double in case of no_files > 1

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