Changeset 951


Ignore:
Timestamp:
Jul 19, 2012 2:22:52 PM (10 years ago)
Author:
hoffmann
Message:

changes for using cross_profiles

Location:
palm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/NCL/.ncl.config.default

    r782 r951  
    828828     
    829829         
    830          var = "all"
     830         var = ",cross_profiles,"
    831831         
    832832         
     
    10491049      ; data type: integer
    10501050      ;
    1051       ; default:   1
     1051      ; default:   -1
    10521052      ;***************************************************
    10531053      if (.not. isvar("no_columns"))then         
    10541054     
    10551055     
    1056          no_columns = 1
     1056         no_columns = -1
    10571057         
    10581058         
     
    10651065      ; data type: integer
    10661066      ;
    1067       ; default:   2
     1067      ; default:   -1
    10681068      ;***************************************************
    10691069      if (.not. isvar("no_rows"))then   
    10701070     
    10711071       
    1072          no_rows = 2
     1072         no_rows = -1
    10731073         
    10741074         
  • palm/trunk/SCRIPTS/NCL/profiles.ncl

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