Index: palm/trunk/SCRIPTS/NCL/.ncl_preferences
===================================================================
--- palm/trunk/SCRIPTS/NCL/.ncl_preferences	(revision 157)
+++ palm/trunk/SCRIPTS/NCL/.ncl_preferences	(revision 161)
@@ -19,11 +19,11 @@
 ___no_lines__________[1;2;...]__________________________________default = 2_______________[ALL]
 2
-___'var'_____________[E,dt,umax,;...]___________________________default = variables_______[ALL]
+___'var'_____________[,E,dt,umax,;...]__________________________default = variables_______[ALL]
 variables
 ___combine___________[0] or [1]_________________________________default = 0_______________[PR]
 0
-___number_comb_______[0] or [1]_________________________________default = 0_______________[PR]
+___number_comb_______[2] or [3]_________________________________default = 0_______________[PR]
 0
-___'c_var'___________[umax,vmax,wmax,;...]______________________default = c_variables_____[PR]
+___'c_var'___________[,us2,vs2,ws2,;...]________________________default = c_variables_____[PR]
 c_variables
 ___dash______________[0] or [1]_________________________________default = 0_______________[PR]
@@ -33,6 +33,6 @@
 ___min_z_____________[0;10;...]_________________________________default = 0_______________[PR]
 0
-___max_z_____________[100;1000;...]_____________________________default = max(z_pr)_______[PR] 
-max(z_pr)
+___max_z_____________[100;1000;...]_____________________________default = max(z_u)________[PR] 
+max(z_u)
 ___'sort'____________[time] or [layer]__________________________default = time____________[CS]
 time
@@ -63,11 +63,13 @@
 ___vector____________[0] or [1]_________________________________default = 0_______________[CS]
 0
-___'vec1'____________[u,;v,;...]________________________________default = vec1____________[CS]
+___'vec1'____________[,u,;,v,;...]______________________________default = vec1____________[CS]
 vec1
-___'vec2'____________[u,;v,;...]________________________________default = vec2____________[CS]
+___'vec2'____________[,u,;,v,;...]______________________________default = vec2____________[CS]
 vec2
-__'plotvec'__________[u,;v,;...]________________________________default = plotvec_________[CS]
+__'plotvec'__________[,u,;,v,;...]______________________________default = plotvec_________[CS]
 plotvec
 ___ref_mag___________[0.1;0.02;...]_____________________________default = 0.05____________[CS]
 0.05
+___over______________[0] or [1]_________________________________default = 0_______________[TS+PR]
+0
 ___end of file ".ncl_preferences"_____________________________________________________________
Index: palm/trunk/SCRIPTS/NCL/cross_sections.ncl
===================================================================
--- palm/trunk/SCRIPTS/NCL/cross_sections.ncl	(revision 157)
+++ palm/trunk/SCRIPTS/NCL/cross_sections.ncl	(revision 161)
@@ -12,5 +12,5 @@
 
    if (isfilepresent("~/.ncl_preferences")) then
-      parameter = asciiread("~/.ncl_preferences",73,"string")
+      parameter = asciiread("~/.ncl_preferences",75,"string")
       delete(parameter@_FillValue)
    else
@@ -208,5 +208,10 @@
    print(" ")
    dim   = dimsizes(vNam)
-  
+   if (dim .EQ. 0) then
+      print(" ")
+      print("There are no data on file")
+      print(" ")
+   end if
+
    ; ***************************************************
    ; set up recourses
@@ -259,26 +264,9 @@
    end if
 
-   ; ***************************************************
-   ; open workstation(s) 
-   ; ***************************************************
-
-   wks_ps  = gsn_open_wks(format_out,file_out)
-   gsn_define_colormap(wks_ps,"rainbow+white")
-
-   plot=new((/no_columns*no_lines/),graphic)
-
-   page = 0
-
-   nc1 = no_columns
-   nl1 = no_lines
-
-   print("plots sorted by '"+sort+"'")
-   print(" ")
-
    ; *********************************************
    ; set up of start_time_step and end_time_step
    ; *********************************************
 
-   t = f->time
+   t    = f->time
    nt = dimsizes(t)
      
@@ -288,7 +276,7 @@
 
    if ( .not. isvar("start_time_step") ) then   	
-      start_time_step = 1
+      start_time_step = 0
       if (parameter(13) .NE. "1") then
-         if (parameter(13) .LE. "0")
+         if (parameter(13) .LT. "1")
             print(" ")
             print("Begin with time step 1")
@@ -302,5 +290,5 @@
             exit
          end if
-         start_time_step = stringtointeger(parameter(13))  
+         start_time_step = stringtointeger(parameter(13))-1  
       end if
    else 
@@ -317,4 +305,5 @@
          exit
       end if
+      start_time_step = start_time_step - 1
    end if
 
@@ -324,5 +313,5 @@
 
    if ( .not. isvar("end_time_step") ) then     	
-      end_time_step = nt
+      end_time_step = nt-1
       if (parameter(15) .NE. "nt") then
          if (parameter(15) .LE. "0")
@@ -344,5 +333,5 @@
             exit
          end if
-         end_time_step = stringtointeger(parameter(15))  
+         end_time_step = stringtointeger(parameter(15))-1  
       end if    
    else
@@ -365,5 +354,8 @@
          exit
       end if
-   end if
+      end_time_step = end_time_step-1 
+   end if
+
+   no_time=(end_time_step-start_time_step)+1
 
    ; ****************************************************
@@ -374,10 +366,4 @@
       if ( .not. isvar("plotvec") ) then
          plotvec = parameter(69)
-         if (parameter(69) .EQ. "plotvec") then
-            print(" ")
-            print("Please indicate the variable where the vector plot shall overlay ('plotvec')")
-            print(" ")
-            exit
-         end if
       end if
       if ( .not. isvar("vec1") ) then
@@ -409,22 +395,21 @@
    ; ****************************************************
 
-   do varn=0,dim-1  
-
-      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" ) then
+   do varn=dim-1,0,1  
+
+      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x" .or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then
          check = False
-      end if
-      if ( vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi") then
-         check = False
-      end if      
+      else
+         check = True
+      end if    
       if (  isvar("var") ) then  
-         check = isStrSubset( var,vNam(varn)+",")
+         check = isStrSubset( var,","+vNam(varn)+",")
       end if
       if (parameter(21) .NE. "variables") then
          var = parameter(21) 
-         check = isStrSubset( var,vNam(varn)+"," ) 
+         check = isStrSubset( var,","+vNam(varn)+"," ) 
       end if   
 
       if(check) then
-       
+         print(vNam(varn))
          data_all = f->$vNam(varn)$
 
@@ -536,4 +521,14 @@
    end do
 
+   if (xyc .EQ. 1) then
+      no_layer = (ze-zs)+1
+   end if
+   if (xzc .EQ. 1) then
+      no_layer = (ye-ys)+1
+   end if
+   if (yzc .EQ. 1) then
+      no_layer = (xe-xs)+1
+   end if
+
    data  = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,(xe-xs)+1/),float)
 
@@ -570,36 +565,38 @@
    v1=0
    v2=0
-
-   do varn=0,dim-1  
+   no_var=0
+   n=0
+
+   do varn=dim-1,0,1  
       
       data_all = f->$vNam(varn)$ 
 
       if (vector .EQ. 1) then    
-         check_vec1 = isStrSubset( vec1,vNam(varn)+",")
-         check_vec2 = isStrSubset( vec2,vNam(varn)+",")
+         check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
+         check_vec2 = isStrSubset( vec2,","+vNam(varn)+",")
       end if
            
-      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" ) then
+      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then
          check = False
-      end if
-      if ( vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi") then
-         check = False
-      end if      
+      else
+         check = True
+      end if  
       if (  isvar("var") ) then  
-         check = isStrSubset( var,vNam(varn)+",")
-      end if
-  
+         check = isStrSubset( var,","+vNam(varn)+",")
+      end if  
       if (parameter(21) .NE. "variables") then
          var = parameter(21) 
-         check = isStrSubset( var,vNam(varn)+"," ) 
+         check = isStrSubset( var,","+vNam(varn)+"," ) 
       end if  
 
       if(check) then 
 
+         no_var=no_var+1
+
          if (vector .EQ. 1) then
-            if (vNam(varn)+"," .EQ. vec1) then
+            if (","+vNam(varn)+"," .EQ. vec1) then
                v1=v1+1
             end if
-            if (vNam(varn)+"," .EQ. vec2) then
+            if (","+vNam(varn)+"," .EQ. vec2) then
                v2=v2+1
             end if
@@ -608,9 +605,7 @@
          data(varn,:,:,:,:)=data_all(0:nt-1,zs:ze,ys:ye,xs:xe) 
          if (check_vec1) then
-            print("in vec1")
             vect1=data_all
          end if
          if (check_vec2) then
-             print("in vec2")
             vect2=data_all
          end if
@@ -626,5 +621,7 @@
       delete(data_all)
 
-   end do  
+   end do
+  
+   plot_panel=new((/no_time*no_layer*no_var/),graphic)
 
    if (vector .EQ. 1) then
@@ -645,4 +642,19 @@
 
    ; ***************************************************
+   ; open workstation(s) 
+   ; ***************************************************
+
+   wks_ps  = gsn_open_wks(format_out,file_out)
+   gsn_define_colormap(wks_ps,"rainbow+white")
+
+   plot=new((/no_time*no_layer*no_var/),graphic)
+
+   page = 0
+   n=0
+
+   print("plots sorted by '"+sort+"'")
+   print(" ")
+
+   ; ***************************************************
    ; create plots 
    ; ***************************************************
@@ -650,91 +662,32 @@
    check = True
 
-   do varn=0,dim-1 
+   do varn=dim-1,0,1 
 
       if (vector .EQ. 1) then    
-         check_vecp = isStrSubset( plotvec,vNam(varn)+",") 
+         check_vecp = isStrSubset( plotvec,","+vNam(varn)+",") 
       end if 
 
-      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" ) then
-         check = False
-      end if
-      if ( vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi") then
-         check = False
-      end if      
+      if ( vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "zu" .or. vNam(varn) .eq. "zwwi" .or. vNam(varn) .eq. "zusi" .or. vNam(varn) .eq. "time" .or. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d" .or. vNam(varn) .eq. "x".or. vNam(varn) .eq. "y" .or. vNam(varn) .eq. "zu_xy" .or. vNam(varn) .eq. "zw_xy" .or. vNam(varn) .eq. "zu1_xy" .or. vNam(varn) .eq. "ind_z_xy" .or. vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz" .or. vNam(varn) .eq. "ind_y_xz" .or. vNam(varn) .eq. "x_yz" .or. vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "ind_x_yz") then
+         check = False 
+      else
+         check = True
+      end if   
       if (  isvar("var") ) then  
-         check = isStrSubset( var,vNam(varn)+",")
+         check = isStrSubset( var,","+vNam(varn)+",")
       end if
       if (parameter(21) .NE. "variables") then
          var = parameter(21) 
-         check = isStrSubset( var,vNam(varn)+"," ) 
+         check = isStrSubset( var,","+vNam(varn)+"," ) 
       end if 
     
       if(check) then
-      
+    
          ; ****************************************************
          ; loops over time and layer
          ; **************************************************** 
 
-         do lo = los, loe           				; lo: loop outer
-            do li = lis, lie, no_lines*no_columns		; li: loop inner 
-               
-               if ( sort .eq. "time" ) then
-                  ta = li
-                  tb = min( (/li + no_lines*no_columns - 1, end_time_step/) )
-                  if ( xyc .eq. 1 ) then
-                     xa = xs
-                     xb = xe
-                     ya = ys
-                     yb = ye
-                     za = lo
-                     zb = lo
-                  end if
-                  if ( xzc .eq. 1 ) then
-                     xa = xs
-                     xb = xe
-                     ya = lo
-                     yb = lo
-                     za = zs
-                     zb = ze
-                  end if
-                  if ( yzc .eq. 1 ) then
-                     xa = lo
-                     xb = lo
-                     ya = ys
-                     yb = ye
-                     za = zs
-                     zb = ze
-                  end if
-               end if
-
-               if ( sort .eq. "layer" ) then
-                  ta = lo
-                  tb = lo
-                  if ( xyc .eq. 1 ) then
-                     xa = xs 
-                     xb = xe
-                     ya = ys
-                     yb = ye
-                     za = li
-                     zb = min( (/li + no_lines*no_columns - 1, ze/) )  
-                  end if
-                  if ( xzc .eq. 1 ) then
-                     xa = xs 
-                     xb = xe
-                     ya = li
-                     yb = min( (/li + no_lines*no_columns - 1, ye/) ) 
-                     za = zs
-                     zb = ze
-                  end if
-                  if ( yzc .eq. 1 ) then
-                     xa = li 
-                     xb = min( (/li + no_lines*no_columns - 1, xe/) ) 
-                     ya = ys
-                     yb = ye
-                     za = zs
-                     zb = ze
-                  end if
-               end if   
-
+         do lo = los, loe           				
+            do li = lis, lie	                        
+                  
                ; ****************************************************
                ; xy cross section
@@ -743,8 +696,8 @@
                if ( xyc .eq. 1 ) then 
                 
-                  cs_res@tiXAxisString   = "x [m]"
-                  cs_res@tiYAxisString   = "y [m]"
+                  cs_res@tiXAxisString = "x [m]"
+                  cs_res@tiYAxisString = "y [m]"
                   cs_res@gsnLeftString = "Plot of "+vNam(varn)
-                  
+      
                   if ( sort .eq. "time" ) then
                      if ( data&z(lo) .eq. -1 ) then
@@ -753,50 +706,71 @@
                         level = "=" + data&z(lo) + "m"
                      end if
-                     do n = 0,tb-ta
-                        cs_res@gsnCenterString = "time=" + data&t(li+n-1) +"s  z"+level
-                        plot1 = gsn_csm_contour(wks_ps,data(varn,li+n-1,lo,:,:),cs_res)  
+                     cs_res@gsnCenterString = "time=" + data&t(li) +"s  z"+level
+                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,lo,:,:),cs_res) 
+                     if (vector .EQ. 1) then
                         if (check_vecp)
                            vecres                  = True            ; vector only resources
   			   vecres@gsnDraw          = False           ; don't draw
-  			   vecres@gsnFrame         = False           ; don't advance frame
-  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
-  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
-  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  	                   vecres@gsnFrame         = False           ; don't advance frame
+  	                   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  		           vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  		           vecres@vcRefLengthF     = 0.05            ; define length of vec ref
   			   vecres@gsnRightString   = " "             ; turn off right string
   			   vecres@gsnLeftString    = " "             ; turn off left string
   			   vecres@tiXAxisString    = " "   
-                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
-                           overlay(plot1, plot_vec)
-                        end if                      
-                        plot(n) = plot1                         
-                     end do
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                           overlay(plot(n), plot_vec)
+                        else  
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                           plot(n) = plot_vec
+                        end if
+                     end if                         
                   end if
 
                   if ( sort .eq. "layer" ) then
-                     do n = 0, zb-za
-                        if ( data&z(li+n) .eq. -1 ) then
-                           level = "-average"
-                        else 
-                           level = "=" + data&z(li+n) + "m"
+                     if ( data&z(li) .eq. -1 ) then
+                        level = "-average"
+                     else 
+                        level = "=" + data&z(li) + "m"
+                     end if
+                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  z"+ level
+                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,li,:,:),cs_res)
+                     if (vector .EQ. 1) then
+                        if (check_vecp)
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = " "             ; turn off right string
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+                           overlay(plot(n), plot_vec)
+                        else  
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+                           plot(n) = plot_vec
                         end if
-                        cs_res@gsnCenterString = "t=" + data&t(lo-1) + "s  z"+ level
-                        plot1 = gsn_csm_contour(wks_ps,data(varn,lo-1,li+n,:,:),cs_res)
-                        if (check_vecp)
-                           vecres                  = True            ; vector only resources
-  			   vecres@gsnDraw          = False           ; don't draw
-  			   vecres@gsnFrame         = False           ; don't advance frame
-  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
-  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
-  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
-  			   vecres@gsnRightString   = " "             ; turn off right string
-  			   vecres@gsnLeftString    = " "             ; turn off left string
-  			   vecres@tiXAxisString    = " "   
-                           plot_vec=gsn_csm_vector(wks_ps,vect1(lo-1,li+n,:,:),vect2(lo-1,li+n,:,:),vecres)
-                           overlay(plot1, plot_vec)
-                        end if 
-                        plot(n) = plot1
-                     end do
+                     end if
                   end if
-
                end if
 
@@ -817,7 +791,7 @@
                         level = "=" + data&y(lo) + "m"
                      end if
-                     do n = 0, tb-ta
-                        cs_res@gsnCenterString = "t=" + data&t(li+n-1) + "s  y"+ level
-                        plot1 = gsn_csm_contour(wks_ps,data(varn,li+n-1,:,lo,:),cs_res)
+                     cs_res@gsnCenterString = "t=" + data&t(li) + "s  y"+ level
+                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
+                     if (vector .EQ. 1) then
                         if (check_vecp)
                            vecres                  = True            ; vector only resources
@@ -830,37 +804,58 @@
   			   vecres@gsnLeftString    = " "             ; turn off left string
   			   vecres@tiXAxisString    = " "   
-                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
-                           overlay(plot1, plot_vec)
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                           overlay(plot(n), plot_vec)
+                        else  
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                           plot(n) = plot_vec
                         end if
-                        plot(n) = plot1  
-                     end do
+                     end if
                   end if
 
                   if ( sort .eq. "layer" ) then
-                     do n = 0, yb-ya
-                        if ( data&z(li+n) .eq. -1 ) then
-                           level = "-average"
-                        else 
-                           level = "=" + data&y(li+n) + "m"
+                     if ( data&z(li) .eq. -1 ) then
+                        level = "-average"
+                     else 
+                        level = "=" + data&y(li) + "m"
+                     end if
+                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  y"+ level
+                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
+                     if (vector .EQ. 1) then
+                        if (check_vecp)
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = " "             ; turn off right string
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+                           overlay(plot(n), plot_vec)
+                        else  
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+                           plot(n) = plot_vec
                         end if
-                        cs_res@gsnCenterString = "t=" + data&t(lo-1) + "s  y"+ level
-                        plot1 = gsn_csm_contour(wks_ps,data(varn,lo-1,:,li+n,:),cs_res)
-                        if (check_vecp)
-                           vecres                  = True            ; vector only resources
-  			   vecres@gsnDraw          = False           ; don't draw
-  			   vecres@gsnFrame         = False           ; don't advance frame
-  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
-  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
-  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
-  			   vecres@gsnRightString   = " "             ; turn off right string
-  			   vecres@gsnLeftString    = " "             ; turn off left string
-  			   vecres@tiXAxisString    = " "   
-                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
-                           overlay(plot1, plot_vec)
-                        end if
-                        plot(n) = plot1
-                     end do
+                     end if
                   end if                  
-                  
                end if
 
@@ -881,7 +876,7 @@
                         level = "=" + data&z(lo) + "m"
                      end if
-                     do n = 0, tb-ta
-                        cs_res@gsnCenterString = "t=" + data&t(li+n-1) + "s  x"+ level
-                        plot1 = gsn_csm_contour(wks_ps,data(varn,li+n-1,:,:,lo),cs_res)
+                     cs_res@gsnCenterString = "t=" + data&t(li) + "s  x"+ level
+                     plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
+    	             if (vector .EQ. 1) then
                         if (check_vecp)
                            vecres                  = True            ; vector only resources
@@ -894,54 +889,83 @@
   			   vecres@gsnLeftString    = " "             ; turn off left string
   			   vecres@tiXAxisString    = " "   
-                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
-                           overlay(plot1, plot_vec)
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                           overlay(plot(n), plot_vec)
+                        else  
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                           plot(n) = plot_vec
                         end if
-                        plot(n) = plot1  
-                     end do
+                     end if
                   end if
 
                   if ( sort .eq. "layer" ) then
-                     do n = 0, xb-xa
-                        if ( data&z(li+n) .eq. -1 ) then
-                           level = "-average"
-                        else 
-                           level = "=" + data&z(li+n) + "m"
+                     if ( data&z(li) .eq. -1 ) then
+                        level = "-average"
+                     else 
+                        level = "=" + data&z(li) + "m"
+                     end if
+                     cs_res@gsnCenterString = "t=" + data&t(lo) + "s  x"+ level
+                     plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
+                     if (vector .EQ. 1) then
+                        if (check_vecp)
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = " "             ; turn off right string
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+                           overlay(plot(n), plot_vec)
+                        else  
+                           vecres                  = True            ; vector only resources
+  			   vecres@gsnDraw          = False           ; don't draw
+  			   vecres@gsnFrame         = False           ; don't advance frame
+  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
+  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
+  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
+  			   vecres@gsnRightString   = "Vector plot of "+vec1+" and "+vec2
+  			   vecres@gsnLeftString    = " "             ; turn off left string
+  			   vecres@tiXAxisString    = " "   
+                           plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+                           plot(n) = plot_vec
                         end if
-                        cs_res@gsnCenterString = "t=" + data&t(lo-1) + "s  x"+ level
-                        plot1 = gsn_csm_contour(wks_ps,data(varn,lo-1,:,:,li+n),cs_res)
-                        if (check_vecp)
-                           vecres                  = True            ; vector only resources
-  			   vecres@gsnDraw          = False           ; don't draw
-  			   vecres@gsnFrame         = False           ; don't advance frame
-  			   vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
-  			   vecres@vcRefMagnitudeF  = ref_mag         ; define vector ref mag
-  			   vecres@vcRefLengthF     = 0.05            ; define length of vec ref
-  			   vecres@gsnRightString   = " "             ; turn off right string
-  			   vecres@gsnLeftString    = " "             ; turn off left string
-  			   vecres@tiXAxisString    = " "   
-                           plot_vec=gsn_csm_vector(wks_ps,vect1(li+n-1,lo,:,:),vect2(li+n-1,lo,:,:),vecres)
-                           overlay(plot1, plot_vec)
-                        end if
-                        plot(n) = plot1
-                     end do
+                     end if 
                   end if
-
-               end if
-
-	       ; ***************************************************
-               ; merge plots onto one page 
-               ; ***************************************************
-         
-               gsn_panel(wks_ps, plot(0:n-1),(/no_lines,no_columns/),cs_resP)
-
-               ; To hold no_lines and no_columns constant, it must be defined again:
-
-               no_lines   = nl1
-               no_columns = nc1
-
+               end if	      
+               n=n+1   
             end do 	
          end do  
       end if 
-   end do   
+   end do
+
+   ; ***************************************************
+   ; merge plots onto one page 
+   ; ***************************************************
+             
+   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
+      gsn_panel(wks_ps,plot(0:(no_time*no_layer*no_var)-1),(/no_var,no_layer*no_time/),cs_resP)
+      print(" ")
+      print("Outputs to .eps or .epsi have only one frame")
+      print(" ")
+   else
+      do np = 0,no_layer*no_time*no_var-1,no_lines*no_columns   
+         if ( np + no_lines*no_columns .gt. (no_layer*no_time*no_var)-1) then 
+            gsn_panel(wks_ps, plot(np:(no_layer*no_time*no_var)-1),(/no_lines,no_columns/),cs_resP)
+         else
+            gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
+         end if
+      end do
+   end if
 
    print(" ")
Index: palm/trunk/SCRIPTS/NCL/profiles.ncl
===================================================================
--- palm/trunk/SCRIPTS/NCL/profiles.ncl	(revision 157)
+++ palm/trunk/SCRIPTS/NCL/profiles.ncl	(revision 161)
@@ -10,5 +10,5 @@
 
    if (isfilepresent("~/.ncl_preferences")) then
-      parameter = asciiread("~/.ncl_preferences",73,"string")
+      parameter = asciiread("~/.ncl_preferences",75,"string")
       delete(parameter@_FillValue)
    else
@@ -71,6 +71,10 @@
    if (combine .EQ. 1) then
       if( .not. isvar("c_var") ) then
-         c_var=" "
-         if (parameter(27) .NE. "c_variables") then
+         if (parameter(27) .EQ. "c_variables") then
+            print(" ") 
+            print("Please select variables for overlaying ('c_var')")
+            print(" ")
+            exit
+         else
             c_var=parameter(27)
          end if
@@ -101,5 +105,17 @@
       end if
    end if
-
+   if ( .not. isvar("over") ) then     		        ; switches overlaying plots on
+      over = 0
+      if (stringtointeger(parameter(73)) .NE. 0) then
+         over = stringtointeger(parameter(73))
+         if (stringtointeger(parameter(73)) .NE. 1) then
+            print(" ")
+            print("Please set 'over' to 0 or 1")
+            print(" ")
+            exit
+         end if   
+      end if
+   end if
+  
    ; ***************************************************
    ; open input file
@@ -108,14 +124,342 @@
    f=addfile( file_in,"r")
    
-   vNam  = getfilevarnames(f) 
+   vNam = getfilevarnames(f) 
    print(" ")
    print("Variable on netCDF file: " + vNam)
    print(" ")
-   dim   = dimsizes(vNam)
-   z_pr  = f->zpt
-   dimz  = dimsizes(z_pr)
+   dim = dimsizes(vNam)
+   if (dim .EQ. 0) then
+      print(" ")
+      print("There are no data on file")
+      print(" ")
+   end if
+
+   do varn = dim-1,0,1
+      if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "NORM")) then
+         varn=varn-1
+         continue
+      end if
+      if (vNam(varn) .EQ. "u")
+         z_u = f->zu
+         break
+      else
+         if (vNam(varn) .EQ. "v")then
+            z_u = f->zv
+            break
+         else
+            if(vNam(varn) .EQ. "pt")then
+               z_u = f->zpt
+               break
+            else
+               if(vNam(varn) .EQ. "vpt")then
+                  z_u = f->zvpt 
+                  break
+               else
+                  if(vNam(varn) .EQ. "lpt")then
+                     z_u = f->zlpt 
+                     break
+                  else   
+                     if(vNam(varn) .EQ. "q")then
+                        z_u = f->zq 
+                        break
+                     else
+                        if(vNam(varn) .EQ. "qv")then
+                           z_u = f->zqv 
+                           break
+                        else
+                           if(vNam(varn) .EQ. "ql")then
+                              z_u = f->zql 
+                              break
+                           else
+                              if(vNam(varn) .EQ. "rho")then
+                                 z_u = f->zrho 
+                                 break
+                              else
+                                 if(vNam(varn) .EQ. "s")then
+                                    z_u = f->zs 
+                                    break
+                                 else
+                                    if(vNam(varn) .EQ. "sa")then
+                                       z_u = f->zsa 
+                                       break
+                                    else
+                                       if(vNam(varn) .EQ. "e")then
+                                          z_u = f->ze 
+                                          break
+                                       else
+                                          if(vNam(varn) .EQ. "es")then
+                                             z_u = f->zes 
+                                             break
+                                          else
+                                             if(vNam(varn) .EQ. "km")then
+                                                z_u = f->zkm 
+                                                break
+                                             else
+                                                if(vNam(varn) .EQ. "kh")then
+                                                   z_u = f->zkh 
+                                                   break
+                                                else
+                                                   if(vNam(varn) .EQ. "l")then
+                                                      z_u = f->zl 
+                                                      break
+                                                   else
+                                                      if(vNam(varn) .EQ. "us2")then
+                                                         z_u = f->zus2 
+                                                         break
+                                                      else
+                                                         if(vNam(varn) .EQ. "vs2")then
+                                                            z_u = f->zvs2 
+                                                            break
+                                                         else
+                                                            if(vNam(varn) .EQ. "pts2")then
+                                                               z_u = f->zpts2 
+                                                               break
+                                                            else
+                                                               if(vNam(varn) .EQ. "wsususodz")then
+                                                                  z_u = f->zwsususodz 
+                                                                  break
+                                                               else
+                                                                  if(vNam(varn) .EQ. "wspsodz")then
+                                                                     z_u = f->zwspsodz 
+                                                                     break
+                                                                  else
+                                                                     if(vNam(varn) .EQ. "wpeodz")then
+                                                                        z_u = f->zwpeodz
+                                                                        break                                                                        
+                                                                     end if
+                                                                  end if
+                                                               end if
+                                                            end if
+                                                         end if
+                                                      end if
+                                                   end if
+                                                end if
+                                             end if
+                                          end if
+                                       end if
+                                    end if
+                                 end if
+                              end if
+                           end if
+                        end if
+                     end if
+                  end if
+               end if
+            end if 
+         end if
+      end if
+      varn=varn-1
+   end do
+
+   do varn=dim-1,0,1
+      if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "NORM")) then
+         varn=varn-1
+         continue
+      end if 
+      if (vNam(varn) .EQ. "w")
+         z_w = f->zw
+         break
+      else
+         if (vNam(varn) .EQ. "wpup")then
+            z_w = f->zwpup
+            break
+         else
+            if(vNam(varn) .EQ. "wsus")then
+               z_w = f->zwsus
+               break
+            else
+               if(vNam(varn) .EQ. "wu")then
+                  z_w = f->zwu 
+                  break
+               else
+                  if(vNam(varn) .EQ. "wpvp")then
+                     z_w = f->zwpvp 
+                     break
+                  else   
+                     if(vNam(varn) .EQ. "wsvs")then
+                        z_w = f->zwsvs 
+                        break
+                     else
+                        if(vNam(varn) .EQ. "wv")then
+                           z_w = f->zwv 
+                           break
+                        else
+                           if(vNam(varn) .EQ. "wptpp")then
+                              z_w = f->zwptpp 
+                              break
+                           else
+                              if(vNam(varn) .EQ. "wspts")then
+                                 z_w = f->zwspts 
+                                 break
+                              else
+                                 if(vNam(varn) .EQ. "wpt")then
+                                    z_w = f->zwpt
+                                    break
+                                 else
+                                    if(vNam(varn) .EQ. "wsptsBC")then
+                                       z_w = f->zwsptsBC 
+                                       break
+                                    else
+                                       if(vNam(varn) .EQ. "wptBC")then
+                                          z_w = f->zwptBC 
+                                          break
+                                       else
+                                          if(vNam(varn) .EQ. "wpvptp")then
+                                             z_w = f->zwpvptp 
+                                             break
+                                          else
+                                             if(vNam(varn) .EQ. "wsvpts")then
+                                                z_w = f->zwsvpts
+                                                break
+                                             else
+                                                if(vNam(varn) .EQ. "wvpt")then
+                                                   z_w = f->zwvpt 
+                                                   break
+                                                else
+                                                   if(vNam(varn) .EQ. "wpqp")then
+                                                      z_w = f->zwpqp 
+                                                      break
+                                                   else
+                                                      if(vNam(varn) .EQ. "wsqs")then
+                                                         z_w = f->zwsqs 
+                                                         break
+                                                      else
+                                                         if(vNam(varn) .EQ. "wq")then
+                                                            z_w = f->zwq 
+                                                            break
+                                                         else
+                                                            if(vNam(varn) .EQ. "wpqvp")then
+                                                               z_w = f->zwpqvp 
+                                                               break
+                                                            else
+                                                               if(vNam(varn) .EQ. "wsqvs")then
+                                                                  z_w = f->zwsqvs 
+                                                                  break
+                                                               else
+                                                                  if(vNam(varn) .EQ. "wqv")then
+                                                                     z_w = f->zwqv 
+                                                                     break
+                                                                  else
+                                                                     if(vNam(varn) .EQ. "wpsp")then
+                                                                        z_w = f->zwpsp
+                                                                        break 
+                                                                     else
+                                                                        if(vNam(varn) .EQ. "wsss")then
+                                                                           z_w = f->zwsss
+                                                                           break 
+                                                                        else
+                                                                           if(vNam(varn) .EQ. "ws")then
+                                                                              z_w = f->zws
+                                                                              break 
+                                                                           else
+                                                                              if(vNam(varn) .EQ. "wpsap")then
+                                                                                 z_w = f->zwpsap
+                                                                                 break 
+                                                                              else
+                                                                                 if(vNam(varn) .EQ. "wssas")then
+                                                                                    z_w = f->zwssas
+                                                                                    break 
+                                                                                 else
+                                                                                    if(vNam(varn) .EQ. "wsa")then
+                                                                                       z_w = f->zwsa
+                                                                                       break 
+                                                                                    else
+                                                                                       if(vNam(varn) .EQ. "wses")then
+                                                                                          z_w = f->zwses
+                                                                                          break 
+                                                                                       else
+                                                                                          if(vNam(varn) .EQ. "ws2")then
+                                                                                             z_w = f->zws2
+                                                                                             break 
+                                                                                          else
+                                                                                             if(vNam(varn) .EQ. "ws3")then
+                                                                                                z_w = f->zws3
+                                                                                                break 
+                                                                                             else
+                                                                                                if(vNam(varn) .EQ. "Sw")then
+                                                                                                   z_w = f->zSw
+                                                                                                   break 
+                                                                                                else
+                                                                                                   if(vNam(varn) .EQ. "ws2pts")then
+                                                                                                      z_w = f->zws2pts
+                                                                                                      break 
+                                                                                                   else
+                                                                                                      if(vNam(varn) .EQ. "wspts2")then
+                                                                                                         z_w = f->zwspts2
+                                                                                                         break                                            
+                                                                                                      end if
+                                                                                                   end if 
+                                                                                                end if
+                                                                                             end if
+                                                                                          end if
+                                                                                       end if
+                                                                                    end if
+                                                                                 end if
+                                                                              end if
+                                                                           end if
+                                                                        end if    
+                                                                     end if
+                                                                  end if
+                                                               end if
+                                                            end if
+                                                         end if
+                                                      end if
+                                                   end if
+                                                end if
+                                             end if
+                                          end if
+                                       end if
+                                    end if
+                                 end if
+                              end if
+                           end if
+                        end if
+                     end if
+                  end if
+               end if
+            end if 
+         end if
+      end if
+      varn=varn-1
+   end do
+   if ( .not. isvar("z_u") ) then
+      zu = 0
+   else
+      zu = 1
+   end if
+   if ( .not. isvar("z_w") ) then
+      zw = 0
+   else 
+      zw = 1
+   end if
+
+   if (zu .EQ. 0 .AND. zw .EQ. 0) then
+      do varn=0,dim-1      
+         if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "NORM")) then
+            check = False
+         else
+            if (.not. isvar("var")) then
+               check = True
+               if (parameter(21) .NE. "variables") then
+                  var=parameter(21)
+                  check = isStrSubset( var,","+vNam(varn)+"," )
+               end if
+            else          
+               check = isStrSubset( var,","+vNam(varn)+"," )
+            end if
+         end if
+         if (check)
+            z = f->$vNam(varn+1)$
+            dimz = dimsizes(z)
+            break
+         end if   
+      end do
+   end if
+ 
+   dimz  = dimsizes(z_w)
+   dimz  = dimsizes(z_u)
    t_all = f->time
    nt    = dimsizes(t_all)
-
+ 
    ; ****************************************************	
    ; start of time step and different types of mistakes that could be done
@@ -125,7 +469,7 @@
       start_time_step = 1
       if (parameter(13) .NE. "1") then
-         if (parameter(13) .LE. "0")
-            print(" ")
-            print("Begin with time step 1")
+         if (parameter(13) .LE. "1")
+            print(" ")
+            print("Begin at least with time step 2")
             print(" ")
             exit
@@ -137,10 +481,10 @@
             exit
          end if
-         start_time_step = stringtointeger(parameter(13))  
+         start_time_step = stringtointeger(parameter(13))-1  
       end if
    else 
       if (start_time_step .LE. 0)
          print(" ")
-         print("Begin with time step 1")
+         print("Begin at least with time step 1")
          print(" ")
          exit
@@ -152,4 +496,5 @@
          exit
       end if
+      start_time_step = start_time_step - 1 
    end if
 
@@ -200,4 +545,5 @@
             exit
       end if
+      end_time_step = end_time_step - 1
    end if
 
@@ -206,5 +552,5 @@
    ; ****************************************************
    
-   legend_label=new(nt-1,double)
+   legend_label=new(nt,double)
    do p=start_time_step,end_time_step
       legend_label(p-start_time_step)=t_all(p)
@@ -224,5 +570,4 @@
    res@gsnFrame                = False
    res@gsnPaperOrientation     = "portrait"
-   ;  res@gsnmaximize             = True
    res@gsnPaperWidth           = 8.27
    res@gsnPaperHeight          = 11.69
@@ -239,11 +584,15 @@
    res@pmLegendSide            = "Top"
    res@xyExplicitLegendLabels  = legend_label
-   res@pmLegendParallelPosF    = 1.4
-   res@pmLegendOrthogonalPosF  = -1.085
+   res@pmLegendParallelPosF    = 1.15
+   res@pmLegendOrthogonalPosF  = -1.0
    res@pmLegendWidthF          = 0.12
    res@pmLegendHeightF         = 0.3
    res@lgLabelFontHeightF     = .02 
-   ;  res@XBLabelConstantSpacingF = 1.0
-
+   res@txFontHeightF      = 0.02
+   res@tiXAxisFontHeightF = 0.02
+   res@tiYAxisFontHeightF = 0.02
+   res@tiXAxisString      = " "    
+   res@tiYAxisString      = "Height [z]"
+   
    if ( dash .eq. 0 ) then
       res@xyMonoDashPattern       = True  
@@ -257,66 +606,36 @@
 
    ; ***************************************************
-   ; set up colors and recourses for combined plot
-   ; ***************************************************
-   
-   if (combine .EQ. 1) then
-
-      if (.not. isvar("number_comb")) then
-         m=0
-         m=stringtointeger(parameter(25))
-         if(m .EQ. 0) then
-            print(" ")
-            print("Please indicate the number of variables you would like to combine ('number_comb')")
-            print(" ")
-            exit
-         end if
-      end if
-    
-      ores=True
-      
-      colors=new(m*(nt-1),integer)
-      do j=0,m-1
-         colors(j+j*(nt-2):(j+1)*(nt-1)-1) = ispan(2,237,235/np)
-      end do   
-      ores@xyLineColors = colors
-
-      dash_oplot=new(m*(nt-1),integer)
-      do j=0,m-1
-         dash_oplot(j+j*(nt-2):(j+1)*(nt-1)-1)=j
-      end do
-      ores@xyDashPatterns = dash_oplot
-
-      ores                         = True
-      ores@gsnDraw                 = False
-      ores@gsnFrame                = False
-      ores@gsnPaperOrientation     = "portrait"
-      ;  ores@gsnmaximize             = True
-      ores@gsnPaperWidth           = 8.27
-      ores@gsnPaperHeight          = 11.69
-      ores@gsnPaperMargin          = 0.79
-      ores@txFont                  = "helvetica"
-      ores@tiMainFont              = "helvetica"
-      ores@tiXAxisFont             = "helvetica"
-      ores@tiYAxisFont             = "helvetica"
-      ores@tmXBLabelFont           = "helvetica"
-      ores@tmYLLabelFont           = "helvetica"
-      ores@lgLabelFont             = "helvetica"
-      ores@tmLabelAutoStride       = True
-      ores@pmLegendDisplayMode     = "Always"
-      ores@pmLegendSide            = "Top"
-      ores@pmLegendParallelPosF    = 1.4
-      ores@pmLegendOrthogonalPosF  = -1.085
-      ores@pmLegendWidthF          = 0.15
-      ores@pmLegendHeightF         = 0.60
-      ores@lgLabelFontHeightF     = .02 
-      ;  ores@XBLabelConstantSpacingF = 1.0
- 
-   end if
-
-   ; ***************************************************
    ; set up graphics for plot
    ; ***************************************************
 
-   plot  = new(dim,graphic)
+   plot = new(dim,graphic)
+   plot_ = new(dim,graphic)
+
+   if (combine .EQ. 1) then
+      if ( .not. isvar("number_comb") ) then        	
+         if (parameter(25) .EQ. "0") then
+            print(" ")
+            print("Please set 'number_comb' to 2 or 3 if you would like to overlay 2 or 3 variables in one plot")
+            print(" ")
+            exit
+         else
+            number_comb=stringtointeger(parameter(25))
+            plot_o = new(number_comb,graphic)   
+         end if
+      else
+         if(number_comb .EQ. 2 .OR. number_comb .EQ. 3) then
+            plot_o = new(number_comb,graphic)     
+         else
+            print(" ")
+            print("Please set 'number_comb' to 2 or 3 if you would like to overlay 2 or 3 variables in one plot")
+            print(" ")
+            exit 
+         end if    
+      end if
+      label=new(number_comb,string)
+      color_o=new(number_comb,integer)
+      mini=new(number_comb,float)
+      maxi=new(number_comb,float)
+   end if
 
    wks=gsn_open_wks(format_out,file_out)
@@ -324,13 +643,13 @@
  
    ; ***************************************************
-   ; indicate plot number
+   ; indicate plot number 
    ; ***************************************************
-
+   
    if (combine .EQ. 1) then
-      n = 1 
+      n = 1
    else
       n = 0
-   end if 
-  
+   end if
+
    ; ***************************************************
    ; set up minimum and maximum height
@@ -340,7 +659,7 @@
       min_z=0
       if (stringtointeger(parameter(33)) .NE. 0) then
-         if (stringtointeger(parameter(33)) .GE. max(z_pr) ) then
-            print(" ")
-            print("Minimum of height ('min_z'="+stringtointeger(parameter(33))+") is greater than available heights (="+max(z_pr)+")")
+         if (stringtointeger(parameter(33)) .GE. max(z_u) ) then
+            print(" ")
+            print("Minimum of height ('min_z'="+stringtointeger(parameter(33))+") is greater than available heights (="+max(z_u)+")")
             print(" ")
             exit
@@ -355,7 +674,7 @@
       end if
    else
-      if (min_z .GE. max(z_pr) ) then
+      if (min_z .GE. max(z_u) ) then
          print(" ")
-         print("Minimum of height ('min_z'="+min_z+") is greater than available heights (="+max(z_pr)+")")
+         print("Minimum of height ('min_z'="+min_z+") is greater than available heights (="+max(z_u)+")")
          print(" ")
          exit
@@ -370,9 +689,9 @@
 
    if (.not. isvar("max_z"))
-      max_z=max(z_pr)
-      if ((parameter(35)) .NE. "max(z_pr)") then
-         if (stringtofloat(parameter(35)) .GE. max(z_pr) ) then
-            print(" ")
-            print("Maximum of height ('max_z'="+parameter(35)+") is greater than available heights (="+max(z_pr)+")")
+      max_z=max(z_u)
+      if ((parameter(35)) .NE. "max(z_u)") then
+         if (stringtofloat(parameter(35)) .GE. max(z_u) ) then
+            print(" ")
+            print("Maximum of height ('max_z'="+parameter(35)+") is greater than available heights (="+max(z_u)+")")
             print(" ")
             exit
@@ -389,5 +708,5 @@
       if (max_z .GE. max(z_pr) ) then
          print(" ")
-         print("Maximum of height ('max_z'="+max_z+") is greater than available heights (="+max(z_pr)+")")
+         print("Maximum of height ('max_z'="+max_z+") is greater than available heights (="+max(z_u)+")")
          print(" ")
          exit
@@ -411,17 +730,166 @@
       end if
    end do  
-   
-   if (combine .EQ. 1) then
-      data_o=new((/m,end_time_step,dimz/),float)
-      data_o!0 = "e"
-      legend_label_oplot=new(m*(nt-1),string)
-      mini=new(m,float)
-      maxi=new(m,float)
-      e=-1
-   end if  
-
-   do varn = 0, dim-1   
-     
-      if ( isStrSubset( vNam(varn), "NORM") .or. isStrSubset( vNam(varn), "time") )
+ 
+   data   = new((/dim,(end_time_step-start_time_step)+1,dimz/),float) 
+   data_0 = new((/(end_time_step-start_time_step)+1,dimz/),float)
+   data_0 = 0.0
+   t      = new((/(end_time_step-start_time_step)+1,dimz/),float) 
+   t      = 0.0
+   unit   = new(dim,string)
+
+   if (over .EQ. 1) then 
+      plot_u         = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniu = 1.E27
+      maxiu =-1.E27
+      plot_v         = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniv = 1.E27
+      maxiv =-1.E27
+      plot_w         = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniw = 1.E27
+      maxiw =-1.E27
+      plot_pt        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minipt = 1.E27
+      maxipt =-1.E27
+      plot_vpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minivpt = 1.E27
+      maxivpt =-1.E27
+      plot_lpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minilpt = 1.E27
+      maxilpt =-1.E27
+      plot_q         = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniq = 1.E27
+      maxiq =-1.E27
+      plot_qv        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniqv = 1.E27
+      maxiqv =-1.E27
+      plot_ql        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniql = 1.E27
+      maxiql =-1.E27
+      plot_rho       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_s         = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_sa        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_e         = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minie = 1.E27
+      maxie =-1.E27
+      plot_es        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minie = 1.E27
+      maxie =-1.E27
+      plot_km        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minie = 1.E27
+      maxie =-1.E27
+      plot_kh        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minie = 1.E27
+      maxie =-1.E27
+      plot_l         = gsn_csm_xy(wks,t,data_0(:,:),res)     
+      plot_wpup      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpup = 1.E27
+      maxiwpup =-1.E27
+      plot_wsus      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsus = 1.E27
+      maxiwsus =-1.E27
+      plot_wu        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwu = 1.E27
+      maxiwu =-1.E27
+      plot_wpvp      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpvp = 1.E27
+      maxiwpvp =-1.E27
+      plot_wsvs      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsvs = 1.E27
+      maxiwsvs =-1.E27
+      plot_wv        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwv = 1.E27
+      maxiwv =-1.E27
+      plot_wpptp     = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpptp = 1.E27
+      maxiwpptp =-1.E27
+      plot_wspts     = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwspts = 1.E27
+      maxiwspts =-1.E27
+      plot_wpt       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpt = 1.E27
+      maxiwpt =-1.E27
+      plot_wsptsBC   = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsptsBC = 1.E27
+      maxiwsptsBC =-1.E27 
+      plot_wptBC     = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwptBC = 1.E27
+      maxiwptBC =-1.E27 
+      plot_wpvptp    = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpvptp = 1.E27
+      maxiwpvptp =-1.E27
+      plot_wsvpts    = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsvpts = 1.E27
+      maxiewsvpts=-1.E27
+      plot_wvpt      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwvpt = 1.E27
+      maxiwvpt =-1.E27 
+      plot_wpqp      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpqp = 1.E27
+      maxiwpqp =-1.E27
+      plot_wsqs      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsqs = 1.E27
+      maxiwsqs =-1.E27
+      plot_wq        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwq = 1.E27
+      maxiwq =-1.E27 
+      plot_wpqvp     = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpqvp = 1.E27
+      maxiwpqvp =-1.E27
+      plot_wsqvs     = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsqvs = 1.E27
+      maxiwsqvs =-1.E27
+      plot_wqv       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwqv = 1.E27
+      maxiwqv =-1.E27
+      plot_wpsp      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpsp = 1.E27
+      maxiwpsp =-1.E27
+      plot_wsss      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsss = 1.E27
+      maxiwsss =-1.E27
+      plot_ws        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniws = 1.E27
+      maxiws =-1.E27
+      plot_wpsap     = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpsap = 1.E27
+      maxiwpsap =-1.E27
+      plot_wssas     = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwssas = 1.E27
+      maxiwssas =-1.E27 
+      plot_wsa       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsa = 1.E27
+      maxiwsa =-1.E27
+      plot_wses      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_us2       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minius2 = 1.E27
+      maxius2 =-1.E27
+      plot_vs2       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      minivs2 = 1.E27
+      maxivs2 =-1.E27
+      plot_ws2       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniws2 = 1.E27
+      maxiws2 =-1.E27
+      plot_pts2      = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_ws3       = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_Sw        = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_ws2pts    = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_wspts2    = gsn_csm_xy(wks,t,data_0(:,:),res)
+      plot_wsususodz = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwsususodz = 1.E27
+      maxiwsususodz =-1.E27
+      plot_wspsodz   = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwspsodz = 1.E27
+      maxiwspsodz =-1.E27
+      plot_wpeodz    = gsn_csm_xy(wks,t,data_0(:,:),res)
+      miniwpeodz = 1.E27
+      maxiwpeodz =-1.E27
+   end if
+ 
+   n_o=0
+   do varn = 0,dim-1 
+        
+      temp = f->$vNam(varn)$
+      
+      if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "NORM")) then
          check = False
       else
@@ -430,96 +898,1597 @@
             if (parameter(21) .NE. "variables") then
                var=parameter(21)
-               check = isStrSubset( var,vNam(varn)+"," )
+               check = isStrSubset( var,","+vNam(varn)+"," )
             end if
          else          
-            check = isStrSubset( var,vNam(varn)+"," )
+            check = isStrSubset( var,","+vNam(varn)+"," )
          end if
       end if
+
+      if (combine .EQ. 1) then         
+         com=isStrSubset(c_var,","+vNam(varn)+"," )      
+         if (com)                       
+            check = False
+            data(varn,:,:) = temp(start_time_step:end_time_step,0:dimz-1)
+            unit(varn) = temp@units
+            if (n_o .GT. number_comb-1) then
+               print(" ")
+               print("Please set 'number_comb' to the right number of overlaying variables ('c_var')")
+               print(" ")
+               exit
+            end if 
+            mini(n_o)=min(data(varn,:,:))
+            maxi(n_o)=max(data(varn,:,:)) 
+            varn=varn+1
+            n_o=n_o+1
+         end if
+      end if
+
+      if(check) then
+               
+         z = f->$vNam(varn+1)$
+  
+         unit(varn) = temp@units
+         data(varn,:,:) = temp(start_time_step:end_time_step,0:dimz-1) 
+ 
+         if (over .EQ. 0) then  
+            res@gsnLeftString      = vNam(varn)
+            res@gsnRightString     = unit(varn)
+            res@trYMinF            = min_z
+            res@trYMaxF            = max_z  
+            res@trXMinF            = min(data(varn,:,:))
+            res@trXMaxF            = max(data(varn,:,:))        
+            plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)  
+         end if       
+         
+            if (vNam(varn) .EQ. "u") then 
+               miniu=min(data(varn,:,:))
+               maxiu=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_u = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniu
+                  res@trXMaxF            = maxiu                  
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)  
+               end if
+            end if
+            if (vNam(varn) .EQ. "v") then
+               miniv=min(data(varn,:,:))
+               maxiv=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_v = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniv
+                  res@trXMaxF            = maxiv                   
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               end if  
+            end if
+            if (vNam(varn) .EQ. "w") then
+               miniw=min(data(varn,:,:))
+               maxiw=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2
+                  plot_w = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniw
+                  res@trXMaxF            = maxiw                   
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if   
+            end if
+
+            if (vNam(varn) .EQ. "pt") then
+               minipt=min(data(varn,:,:))
+               maxipt=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_pt = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minipt
+                  res@trXMaxF            = maxipt
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if  
+            end if
+            if (vNam(varn) .EQ. "vpt") then
+               minivpt=min(data(varn,:,:))
+               maxivpt=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_vpt = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minivpt
+                  res@trXMaxF            = maxivpt 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               end if 
+            end if
+            if (vNam(varn) .EQ. "lpt") then
+               minilpt=min(data(varn,:,:))
+               maxilpt=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2
+                  plot_lpt = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minilpt
+                  res@trXMaxF            = maxilpt 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if 
+            end if
+
+            if (vNam(varn) .EQ. "q") then
+               miniq=min(data(varn,:,:))
+               maxiq=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_q = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniq
+                  res@trXMaxF            = maxiq 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if
+            end if
+            if (vNam(varn) .EQ. "qv") then
+               miniqv=min(data(varn,:,:))
+               maxiqv=max(data(varn,:,:))
+               if (over .EQ. 1) then 
+                  res@xyDashPattern  = 1
+                  plot_qv = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniqv
+                  res@trXMaxF            = maxiqv
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "ql") then
+               miniql=min(data(varn,:,:))
+               maxiql=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2
+                  plot_ql = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniql
+                  res@trXMaxF            = maxiql 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "rho") then                
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:))
+                  plot_q = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:)) 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "s") then                
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:))
+                  plot_s = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:)) 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "sa") then                 
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:)) 
+                  plot_sa = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:))  
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "e") then
+               minie=min(data(varn,:,:))
+               maxie=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_e = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minie
+                  res@trXMaxF            = maxie 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if 
+            end if
+            if (vNam(varn) .EQ. "es") then
+               minies=min(data(varn,:,:))
+               maxies=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_es = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minies
+                  res@trXMaxF            = maxies
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "km") then
+               minikm=min(data(varn,:,:))
+               maxikm=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_km = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minikm
+                  res@trXMaxF            = maxikm
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if 
+            end if
+            if (vNam(varn) .EQ. "kh") then
+               minikh=min(data(varn,:,:))
+               maxikh=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_kh = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minikh
+                  res@trXMaxF            = maxikh 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "l") then                 
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:))
+                  plot_l = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:)) 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wpup") then
+               miniwpup=min(data(varn,:,:))
+               maxiwpup=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wpup = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpup
+                  res@trXMaxF            = maxiwpup
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+            if (vNam(varn) .EQ. "wsus") then
+               miniwsus=min(data(varn,:,:))
+               maxiwsus=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wsus = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsus
+                  res@trXMaxF            = maxiwsus 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "wu") then
+               miniwu=min(data(varn,:,:))
+               maxiwu=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2 
+                  plot_wu = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwu
+                  res@trXMaxF            = maxiwu
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wpvp") then
+               miniwpvp=min(data(varn,:,:))
+               maxiwpvp=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wpvp = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpvp
+                  res@trXMaxF            = maxiwpvp
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+            if (vNam(varn) .EQ. "wsvs") then
+               miniwsvs=min(data(varn,:,:))
+               maxiwsvs=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wsvs = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsvs
+                  res@trXMaxF            = maxiwsvs 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "wv") then
+               miniwv=min(data(varn,:,:))
+               maxiwv=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2 
+                  plot_wv = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwv
+                  res@trXMaxF            = maxiwv 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wpptp") then
+               miniwpptp=min(data(varn,:,:))
+               maxiwpptp=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wpptp = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpptp
+                  res@trXMaxF            = maxiwpptp 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+            if (vNam(varn) .EQ. "wspts") then
+               miniwspts=min(data(varn,:,:))
+               maxiwspts=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wspts = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwspts
+                  res@trXMaxF            = maxiwspts 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "wpt") then
+               miniwpt=min(data(varn,:,:))
+               maxiwpt=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2 
+                  plot_wpt = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpt
+                  res@trXMaxF            = maxiwpt 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wsptsBC") then
+               miniwsptsBC=min(data(varn,:,:))
+               maxiwsptsBC=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wsptsBC = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsptsBC
+                  res@trXMaxF            = maxiwsptsBC 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if 
+            end if             
+            if (vNam(varn) .EQ. "wptBC") then
+               miniwptBC=min(data(varn,:,:))
+               maxiwptBC=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wptBC = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwptBC
+                  res@trXMaxF            = maxiwptBC 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wpvptp") then
+               miniwpvptp=min(data(varn,:,:))
+               maxiwpvptp=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wpvptp = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpvptp
+                  res@trXMaxF            = maxiwpvptp 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+            if (vNam(varn) .EQ. "wsvpts") then
+               miniwsvpts=min(data(varn,:,:))
+               maxiwsvpts=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wsvpts = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsvpts
+                  res@trXMaxF            = maxiwsvpts 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "wvpt") then
+               miniwvpt=min(data(varn,:,:))
+               maxiwvpt=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2 
+                  plot_wvpt = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwvpt
+                  res@trXMaxF            = maxiwvpt 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wpqp") then
+               miniwpqp=min(data(varn,:,:))
+               maxiwpqp=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wpqp = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpqp
+                  res@trXMaxF            = maxiwpqp 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+            if (vNam(varn) .EQ. "wsqs") then
+               miniwsqs=min(data(varn,:,:))
+               maxiwsqs=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wsqs = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsqs
+                  res@trXMaxF            = maxiwsqs 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "wq") then
+               miniwq=min(data(varn,:,:))
+               maxiwq=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2
+                  plot_wq = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwq
+                  res@trXMaxF            = maxiwq 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wpqvp") then
+               miniwpqvp=min(data(varn,:,:))
+               maxiwpqvp=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0 
+                  plot_wpqvp = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpqvp
+                  res@trXMaxF            = maxiwpqvp 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if 
+            end if
+            if (vNam(varn) .EQ. "wsqvs") then
+               miniwsqvs=min(data(varn,:,:))
+               maxiwsqvs=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wsqvs = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsqvs
+                  res@trXMaxF            = maxiwsqvs 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "wqv") then
+               miniwqv=min(data(varn,:,:))
+               maxiwqv=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2
+                  plot_wqv = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwqv
+                  res@trXMaxF            = maxiwqv 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wpsp") then
+               miniwpsp=min(data(varn,:,:))
+               maxiwpsp=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wpsp = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpsp
+                  res@trXMaxF            = maxiwpsp 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if 
+            end if 
+            if (vNam(varn) .EQ. "wsss") then
+               miniwsss=min(data(varn,:,:))
+               maxiwsss=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1 
+                  plot_wsss = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsss
+                  res@trXMaxF            = maxiwsss 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "ws") then
+               miniws=min(data(varn,:,:))
+               maxiws=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2
+                  plot_ws = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniws
+                  res@trXMaxF            = maxiws 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if 
+            end if
+
+            if (vNam(varn) .EQ. "wpsap") then
+               miniwpsap=min(data(varn,:,:))
+               maxiwpsap=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wpsap = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpsap
+                  res@trXMaxF            = maxiwpsap 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if 
+            end if
+            if (vNam(varn) .EQ. "wssas") then
+               miniwssas=min(data(varn,:,:))
+               maxiwssas=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wssas = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwssas
+                  res@trXMaxF            = maxiwssas 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "wsa") then
+               miniwsa=min(data(varn,:,:))
+               maxiwsa=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2
+                  plot_wsa = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsa
+                  res@trXMaxF            = maxiwsa 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wses") then                
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:))
+                  plot_wses = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min(data(varn,:,:))
+                  res@trXMaxF            = max(data(varn,:,:)) 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "us2") then
+               minius2=min(data(varn,:,:))
+               maxius2=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_us2 = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minius2
+                  res@trXMaxF            = maxius2 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if 
+            end if
+            if (vNam(varn) .EQ. "vs2") then
+               minivs2=min(data(varn,:,:))
+               maxivs2=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_vs2 = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minivs2
+                  res@trXMaxF            = maxivs2 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               end if
+            end if
+            if (vNam(varn) .EQ. "ws2") then
+               miniws2=min(data(varn,:,:))
+               maxiws2=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2 
+                  plot_ws2 = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniws2
+                  res@trXMaxF            = maxiws2 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "pts2") then
+               minipts2=min(data(varn,:,:))
+               maxipts2=max(data(varn,:,:))                 
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minipts2
+                  res@trXMaxF            = maxipts2
+                  plot_pts2 = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = minipts2
+                  res@trXMaxF            = maxipts2 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "ws3") then
+               miniws3=min(data(varn,:,:))
+               maxiws3=max(data(varn,:,:))                 
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniws3
+                  res@trXMaxF            = maxiws3
+                  plot_ws3 = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniws3
+                  res@trXMaxF            = maxiws3 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "Sw") then
+               miniSw=min(data(varn,:,:))
+               maxiSw=max(data(varn,:,:))                 
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniSw
+                  res@trXMaxF            = maxiSw
+                  plot_Sw = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniSw
+                  res@trXMaxF            = maxiSw 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "ws2pts") then
+               miniws2pts=min(data(varn,:,:))
+               maxiws2pts=max(data(varn,:,:))                 
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniws2pts
+                  res@trXMaxF            = maxiws2pts
+                  plot_ws2pts = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniws2pts
+                  res@trXMaxF            = maxiws2pts 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wspts2") then
+               miniwspts2=min(data(varn,:,:))
+               maxiwspts2=max(data(varn,:,:))                 
+               if (over .EQ. 1) then
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwspts2
+                  res@trXMaxF            = maxiwspts2
+                  plot_wspts2 = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwspts2
+                  res@trXMaxF            = maxiwspts2 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "wsususodz") then
+               miniwsususodz=min(data(varn,:,:))
+               maxiwsususodz=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 0
+                  plot_wsususodz = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwsususodz
+                  res@trXMaxF            = maxiwsususodz 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+               end if  
+            end if
+            if (vNam(varn) .EQ. "wspsodz") then
+               miniwspsodz=min(data(varn,:,:))
+               maxiwspsodz=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 1
+                  plot_wspsodz = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwspsodz
+                  res@trXMaxF            = maxiwspsodz 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
+               end if 
+            end if
+            if (vNam(varn) .EQ. "wpeodz") then
+               miniwpeodz=min(data(varn,:,:))
+               maxiwpeodz=max(data(varn,:,:))
+               if (over .EQ. 1) then
+                  res@xyDashPattern  = 2 
+                  plot_wpeodz = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
+               else
+                  res@gsnLeftString      = vNam(varn)
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = miniwpeodz
+                  res@trXMaxF            = maxiwpeodz 
+                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+               end if
+            end if 
+         if (over .EQ. 0) then
+            n=n+1  
+         end if      
+         varn=varn+1
+      end if      
+      delete(temp) 
+   end do
+   
+   if (over .EQ. 1 ) then
+
+      overlay(plot_u,plot_v)
+      overlay(plot_u,plot_w)
+      u=0
+      overlay(plot_pt,plot_vpt)
+      overlay(plot_pt,plot_lpt)
+      pt=0
+      overlay(plot_q,plot_qv)
+      overlay(plot_q,plot_ql)
+      q=0
+      overlay(plot_e,plot_es)
+      e=0
+      overlay(plot_km,plot_kh)
+      km=0
+      overlay(plot_wpup,plot_wsus)
+      overlay(plot_wpup,plot_wu)
+      wpup=0
+      overlay(plot_wpvp,plot_wsvs)
+      overlay(plot_wpvp,plot_wv)
+      wpvp=0
+      overlay(plot_wpptp,plot_wspts)
+      overlay(plot_wpptp,plot_wpt)
+      wpptp=0
+      overlay(plot_wsptsBC,plot_wptBC)
+      wsptsBC=0
+      overlay(plot_wpvptp,plot_wsvpts)
+      overlay(plot_wpvptp,plot_wvpt)
+      wpvptp=0
+      overlay(plot_wpqp,plot_wsqs)
+      overlay(plot_wpqp,plot_wq)
+      wpqp=0
+      overlay(plot_wpqvp,plot_wsqvs)
+      overlay(plot_wpqvp,plot_wqv)
+      wpqvp=0
+      overlay(plot_wpsp,plot_wsss)
+      overlay(plot_wpsp,plot_ws)
+      wpsp=0
+      overlay(plot_wpsap,plot_wssas)
+      overlay(plot_wpsap,plot_wsa)
+      wpsap=0
+      overlay(plot_us2,plot_vs2)
+      overlay(plot_us2,plot_ws2)
+      us2=0
+      overlay(plot_wsususodz,plot_wspsodz)
+      overlay(plot_wsususodz,plot_wpeodz)
+      wsususodz=0
+
+   end if
+
+   if (over .EQ. 1) then 
+   
+      do varn = 0,dim-1    
+      
+         if ( isStrSubset( vNam(varn), "time") .OR. isStrSubset( vNam(varn), "NORM")) then
+            check = False
+         else
+            if (.not. isvar("var")) then
+               check = True
+               if (parameter(21) .NE. "variables") then
+                  var=parameter(21)
+                  check = isStrSubset( var,","+vNam(varn)+"," )
+               end if
+            else          
+               check = isStrSubset( var,","+vNam(varn)+"," )
+            end if
+         end if  
+
+         if (check)
+           
+            z = f->$vNam(varn+1)$
+           
+            ;res@gsnLeftString      = vNam(varn)
+            res@gsnRightString     = unit(varn)
+            res@trYMinF            = min_z
+            res@trYMaxF            = max_z  
+            res@trXMinF            = min(data(varn,:,:))
+            res@trXMaxF            = max(data(varn,:,:))
+            plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)  
+            
+            if (vNam(varn) .EQ. "u" .OR. vNam(varn) .EQ. "v" .OR. vNam(varn) .EQ. "w") then
+               if (u .EQ. 0) then
+                  res@gsnLeftString      = "u, v and w"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniu,miniv,miniw/))
+                  res@trXMaxF            = max((/maxiu,maxiv,maxiw/))
+                  if (vNam(varn) .EQ. "w") then 
+                     plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+                  else
+                     plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+                  end if
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"u","v","w"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_u)
+                  u=1
+               else
+                  varn=varn+1
+                  continue
+               end if        
+            end if  
      
-      if (combine .EQ. 1) then
-         if (c_var .EQ. " ") then
-            print(" ")
-            print("Please indicate the variables you would like to combine ('c_var')")
-            print(" ")
-            exit
-         end if               
-         com=isStrSubset(c_var,vNam(varn)+",")      
-         if (com) then
-            e=e+1
+            if (vNam(varn) .EQ. "pt" .OR. vNam(varn) .EQ. "vpt" .OR. vNam(varn) .EQ. "lpt") then
+               if (pt .EQ. 0) then
+                  res@gsnLeftString      = "pt, vpt and lpt"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/minipt,minivpt,minilpt/))
+                  res@trXMaxF            = max((/maxipt,maxivpt,maxilpt/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"pt","vpt","lpt"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_pt)
+                  pt=1
+               else
+                  varn=varn+1
+                  continue       
+               end if
+            end if            
+            if (vNam(varn) .EQ. "q" .OR. vNam(varn) .EQ. "qv" .OR. vNam(varn) .EQ. "ql") then
+               if (q .EQ. 0) then
+                  res@gsnLeftString      = "q, qv and ql"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniq,miniqv,miniql/))
+                  res@trXMaxF            = max((/maxiq,maxiqv,maxiql/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"q","qv","ql"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_q)
+                  q=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if           
+            if (vNam(varn) .EQ. "rho") then
+               plot(n) = plot_rho
+            end if
+            if (vNam(varn) .EQ. "s") then
+               plot(n) = plot_s
+            end if
+            if (vNam(varn) .EQ. "sa") then
+               plot(n) = plot_sa
+            end if
+            if (vNam(varn) .EQ. "e" .OR. vNam(varn) .EQ. "es") then
+               if (e .EQ. 0) then
+                  res@gsnLeftString      = "e and es"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/minie,minies/))
+                  res@trXMaxF            = max((/maxie,maxies/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"e","es"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_e)
+                  e=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if            
+            if (vNam(varn) .EQ. "km" .OR. vNam(varn) .EQ. "kh") then
+               if (km .EQ. 0) then
+                  res@gsnLeftString      = "km and kh"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/minikm,minikh/))
+                  res@trXMaxF            = max((/maxikm,maxikh/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"km","kh"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_km) 
+                  km=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if           
+            if (vNam(varn) .EQ. "l") then
+               plot(n) = plot_l
+            end if
+            if (vNam(varn) .EQ. "wpup" .OR. vNam(varn) .EQ. "wsus" .OR. vNam(varn) .EQ. "wu") then
+               if (wpup .EQ. 0) then
+                  res@gsnLeftString      = "wpup, wsus and wu"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwpup,miniwsus,miniwu/))
+                  res@trXMaxF            = max((/maxiwpup,maxiwsus,maxiwu/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wpup","wsus","wu"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wpup)
+                  wpup=1
+               else
+                  varn=varn+1
+                  continue    
+               end if
+            end if
+            if (vNam(varn) .EQ. "wpvp" .OR. vNam(varn) .EQ. "wsvs" .OR.vNam(varn) .EQ. "wv") then
+               if (wpvp .EQ. 0) then
+                  res@gsnLeftString      = "wpvp, wsus and wv"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwpvp,miniwsvs,miniwv/))
+                  res@trXMaxF            = max((/maxiwpvp,maxiwsvs,maxiwv/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wpvp","wsvs","wv"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wpvp)
+                  wpup=0
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if
+            if (vNam(varn) .EQ. "wpptp" .OR. vNam(varn) .EQ. "wspts" .OR. vNam(varn) .EQ. "wpt") then
+               if (wpptp .EQ. 0) then
+                  res@gsnLeftString      = "wpptp, wspts and wv"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwpptp,miniwspts,miniwpt/))
+                  res@trXMaxF            = max((/maxiwpptp,maxiwspts,maxiwpt/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wpptp","wspts","wpt"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wpptp)
+                  wpptp=1
+               else
+                  varn=varn+1
+                  continue    
+               end if
+            end if
+            if (vNam(varn) .EQ. "wsptsBC" .OR. vNam(varn) .EQ. "wptBC") then
+               if (wsptsBC .EQ. 0) then
+                  res@gsnLeftString      = "wsptsBC and wptBC"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwsptsBC,miniwptBC/))
+                  res@trXMaxF            = max((/maxiwsptsBC,maxiwptBC/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wsptsBC","wptBC"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wsptsBC)
+                  wsptsBC=1
+               else
+                  varn=varn+1
+                  continue   
+               end if  
+            end if             
+            if (vNam(varn) .EQ. "wpvptp" .OR. vNam(varn) .EQ. "wsvpts" .OR. vNam(varn) .EQ. "wvpt") then
+               if (wpvptp .EQ. 0) then
+                  res@gsnLeftString      = "wpvptp, wsvpts and wv"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwpvptp,miniwsvpts,miniwvpt/))
+                  res@trXMaxF            = max((/maxiwpvptp,maxiwsvpts,maxiwvpt/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wpvptp","wsvpts","wvpt"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wpvptp)
+                  wpvptp=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if
+            if (vNam(varn) .EQ. "wpqp" .OR. vNam(varn) .EQ. "wsqs" .OR. vNam(varn) .EQ. "wq") then
+               if (wpqp .EQ. 0) then 
+                  res@gsnLeftString      = "wpqp, wsqs and wq"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwpqp,miniwsqs,miniwq/))
+                  res@trXMaxF            = max((/maxiwpqp,maxiwsqs,maxiwq/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wpqp","wsqs","wq"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wpqp)
+                  wpqp=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if
+            if (vNam(varn) .EQ. "wpqvp" .OR. vNam(varn) .EQ. "wsqvs" .OR. vNam(varn) .EQ. "wqv") then
+               if (wpqvp .EQ. 0) then
+                  res@gsnLeftString      = "wpqvp, wsqvs and wqv"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwpqp,miniwsqvs,miniwqv/))
+                  res@trXMaxF            = max((/maxiwpqp,maxiwsqvs,maxiwqv/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wpqvp","wsqvs","wqv"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wpqvp) 
+                  wpqvp=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if
+            if (vNam(varn) .EQ. "wpsp" .OR. vNam(varn) .EQ. "wsss" .OR. vNam(varn) .EQ. "ws") then
+               if (wpsp .EQ. 0) then
+                  res@gsnLeftString      = "wpsp, wsss and ws"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwpsp,miniwsss,miniws/))
+                  res@trXMaxF            = max((/maxiwpsp,maxiwsss,maxiws/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wpsp","wsss","ws"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wpsp)
+                  wpsp=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if 
+            if (vNam(varn) .EQ. "wpsap" .OR.vNam(varn) .EQ. "wssas" .OR. vNam(varn) .EQ. "wsa") then
+               if (wpsap .EQ. 0) then
+                  res@gsnLeftString      = "wpsap, wssas and wsa"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwpsap,miniwssas,miniwsa/))
+                  res@trXMaxF            = max((/maxiwpsap,maxiwssas,maxiwsa/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wpsap","wssas","wsa"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wpsap)
+                  wpsap=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if
+            if (vNam(varn) .EQ. "wses") then
+               plot(n) = plot_wses
+            end if
+            if (vNam(varn) .EQ. "us2" .OR. vNam(varn) .EQ. "vs2" .OR. vNam(varn) .EQ. "ws2") then
+               if (us2 .EQ. 0) then
+                  res@gsnLeftString      = "us2, vs2 and ws2"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/minius2,minivs2,miniws2/))
+                  res@trXMaxF            = max((/maxius2,maxivs2,maxiws2/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"us2","vs2","ws2"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_us2)
+                  us2=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if
+            if (vNam(varn) .EQ. "pts2") then
+               plot(n) = plot_pts2 
+            end if
+            if (vNam(varn) .EQ. "ws3") then
+               plot(n) = plot_ws3
+            end if
+            if (vNam(varn) .EQ. "Sw") then
+               plot(n) = plot_Sw
+            end if
+            if (vNam(varn) .EQ. "ws2pts") then
+               plot(n) = plot_ws2pts
+            end if
+            if (vNam(varn) .EQ. "wspts2") then
+               plot(n) = plot_wspts2 
+            end if
+            if (vNam(varn) .EQ. "wsususodz" .OR. vNam(varn) .EQ. "wspsodz" .OR. vNam(varn) .EQ. "wpeodz" ) then
+               if (wsususodz .EQ. 0) then 
+                  res@gsnLeftString      = "wsususodz, wspsodz and ws2"
+                  res@gsnRightString     = unit(varn)
+                  res@trXMinF            = min((/miniwsususodz,miniwspsodz,miniwpeodz/))
+                  res@trXMaxF            = max((/maxiwsususodz,maxiwspsodz,maxiwpeodz/))
+
+                  plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.12           
+                  lgres@vpHeightF          = 0.1          
+                  lgres@lgDashIndexes      = (/0,1,2/)
+                  lbid = gsn_create_legend(wks,3,(/"wsususodz","wspsodz","wpeodz"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.65		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot(n),lbid,amres)                  
+                  overlay(plot(n),plot_wsususodz)
+                  wsususodz=1
+               else
+                  varn=varn+1
+                  continue   
+               end if
+            end if     
+            n=n+1    
+            varn=varn+1
+         end if
+      end do 
+   end if     
+
+   if (combine .EQ. 1) then
+     
+      check=True      
+      n_o=0
+      do varn = 0,dim-1 
+  
+         if ( isStrSubset( vNam(varn), "NORM") .or. isStrSubset( vNam(varn), "time") )
             check = False
-            temp = f->$vNam(varn)$            
-            data_o!1 = temp!0
-            data_o!2 = temp!1           
-            data_o(e,:,:) = temp(start_time_step:end_time_step,:)
-            mini(e) = min(data_o(e,:,:)) 
-            maxi(e) = max(data_o(e,:,:))
-            
-            do j=e,e
-               do p=start_time_step,end_time_step
-                  legend_label_oplot(e*end_time_step+p-start_time_step)=" "+flt2string(doubletofloat(t_all(p)))+" "+vNam(varn)
-               end do
-            end do
-           
-            delete(temp)
-            varn = varn + 1
+         else
+            if (.not. isvar("var")) then
+               check = True
+               if (parameter(21) .NE. "variables") then
+                  var=parameter(21)
+                  check = isStrSubset( var,","+vNam(varn)+"," )
+               end if
+            else          
+               check = isStrSubset( var,","+vNam(varn)+"," )
+            end if
+         end if
+     
+         if(check) then
+
+            z = f->$vNam(varn+1)$ 
+            com=isStrSubset( c_var,","+vNam(varn)+"," )
+     
+            if (com) 
+                             
+               if (n_o .EQ. 1) then
+                  res@xyDashPattern  = 1 
+                                  
+               else            
+                  if (n_o .EQ. 2) then
+                     res@xyDashPattern  = 2
+                  else
+                     res@xyDashPattern  = 0
+                     res@gsnLeftString  = c_var
+                     res@gsnRightString = unit(varn)
+                     res@trXMinF        = min(mini)
+                     res@trXMaxF        = max(maxi)
+                  end if
+               end if
+               label(n_o)=vNam(varn)
+               color_o(n_o)=237
+               plot_o(n_o)=gsn_csm_xy(wks,data(varn,:,:),z,res) 
+               n_o=n_o+1
+            end if          
+            varn=varn+1   
+         end if
+      end do
+   
+      if(number_comb .EQ. 2)then
+         overlay(plot_o(0),plot_o(1))
+      end if
+      if(number_comb .EQ. 3)then
+         overlay(plot_o(0),plot_o(1))
+         overlay(plot_o(0),plot_o(2))
+      end if
+
+      ; ***************************************************
+      ; legend for combined plot
+      ; ***************************************************
+      
+      lgres                    = True
+      lgMonoDashIndex          = False
+      lgres@lgDashIndexes      = (/0,1,2/)
+      lgres@lgLabelFont        = "helvetica"    
+      lgres@lgLabelFontHeightF = .1            
+      lgres@vpWidthF           = 0.12           
+      lgres@vpHeightF          = 0.1           
+  
+      lbid = gsn_create_legend(wks,number_comb,label,lgres)       
+
+      amres = True
+      amres@amParallelPosF   = 0.65		    
+      amres@amOrthogonalPosF = -0.2            
+      annoid1 = gsn_add_annotation(plot_o(0),lbid,amres) 
     
-         end if        
-      end if  
- 
-      if(check) then
-
-         temp = f->$vNam(varn)$      
-         data = temp(start_time_step:end_time_step,:)   
-         testx= data*100  
-         print("         plot of " + vNam(varn))
-         data!1 = "z"
-  
-         res@gsnLeftString      = "PROFILE plot of"
-         res@gsnRightString     = vNam(varn)
-         res@txFontHeightF      = 0.02
-         res@tiXAxisFontHeightF = 0.02
-         res@tiYAxisFontHeightF = 0.02
-         res@tiXAxisString      = " "    
-         res@tiYAxisString      = "Height [z]"
-         res@trXMinF            = min(data(:,:))
-         res@trXMaxF            = max(data(:,:))
-         res@trYMinF            = min_z
-         res@trYMaxF            = max_z    
-        
-         plot(n) = gsn_csm_xy(wks,data,z_pr,res)  
-   
-         delete(temp)
-         delete(data)
-         n = n + 1 
-         varn = varn + 1
-       
-      end if      
-
-   end do
-
-   if (combine .EQ. 1) then
-      data_all=new((/m*end_time_step,dimz/),float)
-      do j=0,e
-         data_all(start_time_step-1+j*end_time_step:(j+1)*end_time_step-1,:) = data_o(j,:,:)
-      end do  
-
-      ores@gsnLeftString      = "combined PROFILE plot of" 
-      ores@gsnRightString     = c_var
-      ores@txFontHeightF      = 0.02
-      ores@tiXAxisFontHeightF = 0.02
-      ores@tiYAxisFontHeightF = 0.02
-      ores@tiXAxisString      = " "    
-      ores@tiYAxisString      = "Height [z]"
-      ores@trXMinF            = min(mini(:))
-      ores@trXMaxF            = max(maxi(:))
-      ores@trYMinF            = min_z
-      ores@trYMaxF            = max_z   
-      ores@xyExplicitLegendLabels  = legend_label_oplot
-
-      print(" ")
-      print("combined plot of " + c_var) 
- 
-      plot(0) = gsn_csm_xy(wks,data_all,z_pr,ores) 
-      
+      plot(0) = plot_o(0)
+
    end if
 
@@ -527,12 +2496,20 @@
    ; merge plots onto one page 
    ; ***************************************************
-  
-   do i = 0, n-1, no_lines*no_columns 
-      if( (i+no_lines*no_columns) .gt. (n-1)) then
-        gsn_panel(wks,plot(i:n-1),(/no_lines,no_columns/),resP)
-      else
-        gsn_panel(wks,plot(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
-      end if
+ print(n)
+   do m=0,n-1
+      plot_(m)=plot(n-1-m)
    end do
+
+   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
+      gsn_panel(wks,plot_(0:n-1),(/n,1/),resP)
+   else    
+      do i = 0,n-1, no_lines*no_columns
+         if( (i+no_lines*no_columns) .gt. (n-1)) then
+            gsn_panel(wks,plot_(i:n-1),(/no_lines,no_columns/),resP)
+         else
+           gsn_panel(wks,plot_(i:i+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
+         end if
+      end do
+   end if
 
    print(" ")
Index: palm/trunk/SCRIPTS/NCL/timeseries.ncl
===================================================================
--- palm/trunk/SCRIPTS/NCL/timeseries.ncl	(revision 157)
+++ palm/trunk/SCRIPTS/NCL/timeseries.ncl	(revision 161)
@@ -11,5 +11,5 @@
    
    if (isfilepresent("~/.ncl_preferences")) then
-      parameter = asciiread("~/.ncl_preferences",73,"string")
+      parameter = asciiread("~/.ncl_preferences",75,"string")
       delete(parameter@_FillValue)
    else
@@ -23,5 +23,5 @@
    ; set up default parameter values and strings if not assigned in prompt or parameter list
    ; ***************************************************
-   
+  
    if ( .not. isvar("file_in") ) then     		; path+name of input file      
       if (parameter(7) .EQ. "input file") then
@@ -60,4 +60,16 @@
    if ( .not. isvar("var") ) then   			; variable name
       check = True
+   end if
+   if ( .not. isvar("over") ) then     		        ; switches overlaying plots on
+      over = 0
+      if (stringtointeger(parameter(73)) .NE. 0) then
+         over = stringtointeger(parameter(73))
+         if (stringtointeger(parameter(73)) .NE. 1) then
+            print(" ")
+            print("Please set 'over' to 0 or 1")
+            print(" ")
+            exit
+         end if   
+      end if
    end if 
   
@@ -68,11 +80,17 @@
    f  = addfile(file_in , "r" )
 
-   vNam  = getfilevarnames(f) 
+   vNam  = getfilevarnames(f)
    print(" ")
    print("Variable on netCDF file: " + vNam)
    print(" ")
-   dim   = dimsizes(vNam)
+   dim = dimsizes(vNam)
+   if (dim .EQ. 0) then
+      print(" ")
+      print("There are no data on file")
+      print(" ")
+   end if
+
    t_all = f->time
-   nt    = dimsizes(t_all)
+   nt  = dimsizes(t_all)
 
    ; ****************************************************	
@@ -81,7 +99,7 @@
 
    if ( .not. isvar("start_time_step") ) then   	
-      start_time_step = 1
-      if (parameter(13) .LE. "1") then
-         if (parameter(13) .EQ. "0")
+      start_time_step = 0
+      if (parameter(13) .NE. "1") then
+         if (parameter(13) .LE. "0")
             print(" ")
             print("Begin with time step 1")
@@ -95,5 +113,5 @@
             exit
          end if
-         start_time_step = stringtointeger(parameter(13))  
+         start_time_step = stringtointeger(parameter(13))-1  
       end if
    else 
@@ -110,4 +128,5 @@
          exit
       end if
+      start_time_step = start_time_step - 1
    end if
 
@@ -137,5 +156,5 @@
             exit
          end if
-         end_time_step = stringtointeger(parameter(15))  
+         end_time_step = stringtointeger(parameter(15))-1 
       end if    
    else
@@ -153,9 +172,10 @@
       end if
       if (end_time_step .LT. start_time_step)
-            print(" ")
-            print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
-            print(" ")
-            exit
-      end if
+         print(" ")
+         print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
+         print(" ")
+         exit
+      end if
+      end_time_step = end_time_step - 1
    end if  
 
@@ -181,6 +201,9 @@
    res@tiYAxisFont             = "helvetica"
    res@tmXBLabelFont           = "helvetica"
-   res@tmYLLabelFont           = "helvetica"  
-   res@xyLineColors            = (/237/)
+   res@tmYLLabelFont           = "helvetica" 
+   res@xyLineColors            = (/237/) 
+   
+   res@lgLabelFontHeightF     = .02 
+
    resP                        = True
    resP@txFont                 = "helvetica"
@@ -191,4 +214,5 @@
    res@vpWidthF=4
 
+   txres = True
 
    ; ***************************************************
@@ -198,51 +222,634 @@
    wks_ps = gsn_open_wks(format_out,file_out)		
    gsn_define_colormap(wks_ps,"rainbow+white") 
-   plot_ps=new(dim,graphic)                            	 
+   plot_ps=new(dim,graphic) 
+                           	 
    n=0 
- 
-   do varn = 0, dim-1  
+   minE=1.E27
+   maxE=-1.E27
+   minus=1.E27
+   maxus=-1.E27
+   minu=1.E27
+   maxu=-1.E27
+   minz=1.E27
+   maxz=-1.E27
+   minw=1.E27
+   maxw=-1.E27
+   minp=1.E27
+   maxp=-1.E27
+   mins=1.E27
+   maxs=-1.E27
+
+   data   = new((/dim,(end_time_step-start_time_step)+1/),float)
+   unit   = new(dim,string)
+   data_0 = new((end_time_step-start_time_step)+1,float)
+   data_0 = 0.0
+   mini   = new(dim,float)
+   maxi   = new(dim,float)
+   
+   if (over .EQ. 1) then 
+      plot_E       = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_Es      = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_us      = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_ws      = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_umax    = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_vmax    = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_wmax    = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_z_i_wpt = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_z_i_pt  = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_wpptp0  = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_wpptp   = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_wpt     = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_pt_0_   = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_pt_zp_  = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_splptx  = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_splpty  = gsn_csm_xy(wks_ps,t,data_0(:),res)
+      plot_splptz  = gsn_csm_xy(wks_ps,t,data_0(:),res)
+   end if
+  
+   do varn = dim-1,0,1 
+
+      if( isStrSubset (vNam(varn), "time") )
+         check = False
+      else
+         check = True
+      end if 
       if( isvar("var") ) then     		
-         check = isStrSubset( var,vNam(varn)+"," )
+         check = isStrSubset( var,","+vNam(varn)+"," )
       end if
       if (parameter(21) .NE. "variables") then
          var = parameter(21) 
-         check = isStrSubset( var,vNam(varn)+"," ) 
+         check = isStrSubset( var,","+vNam(varn)+"," )
       end if 
     
+     
+      if(check) then
+       
+         data_all = f ->$vNam(varn)$ 
+         unit(varn) = data_all@units
+                
+         data(varn,:)=data_all(start_time_step:end_time_step)
+         
+         if (over .EQ. 1) then
+
+            mini(varn) = min(data(varn,:))
+            maxi(varn) = max(data(varn,:))
+           
+            if (vNam(varn) .EQ. "E" .OR. vNam(varn) .EQ. "Es") then
+               if (mini(varn) .EQ. maxi(varn)) then
+                  if (min(data(varn,:)) .EQ. 0)then
+                     mini(varn)= mini(varn)-1.
+                     maxi(varn)= maxi(varn)+1.
+                  end if
+                  if (min(data(varn,:)) .LT. 0)then
+                     mini(varn)= mini(varn)-1.+(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
+                  end if
+                  if (min(data(varn,:)) .GT. 0)then
+                     mini(varn)= mini(varn)-1.-(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
+                  end if
+               end if
+               minE=min((/minE,mini(varn)/))
+               maxE=max((/maxE,maxi(varn)/))
+            end if
+
+            if (vNam(varn) .EQ. "us" .OR. vNam(varn) .EQ. "ws") then
+               if (mini(varn) .EQ. maxi(varn)) then
+                  if (min(data(varn,:)) .EQ. 0)then
+                     mini(varn)= mini(varn)-1.
+                     maxi(varn)= maxi(varn)+1.
+                  end if
+                  if (min(data(varn,:)) .LT. 0)then
+                     mini(varn)= mini(varn)-1.+(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
+                  end if
+                  if (min(data(varn,:)) .GT. 0)then
+                     mini(varn)= mini(varn)-1.-(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
+                  end if
+               end if
+               minus=min((/minus,mini(varn)/))
+               maxus=max((/maxus,maxi(varn)/))
+            end if
+
+            if (vNam(varn) .EQ. "umax" .OR. vNam(varn) .EQ. "vmax" .OR. vNam(varn) .EQ. "wmax") then
+               if (mini(varn) .EQ. maxi(varn)) then
+                  if (mini(varn) .EQ. 0)then
+                     mini(varn)= mini(varn)-1.
+                     maxi(varn)= maxi(varn)+1.
+                  end if
+                  if (mini(varn) .LT. 0)then
+                     mini(varn)= mini(varn)-1.+(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
+                  end if
+                  if (mini(varn) .GT. 0)then
+                     mini(varn)= mini(varn)-1.-(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
+                  end if
+               end if
+               minu=min((/minu,mini(varn)/))
+               maxu=max((/maxu,maxi(varn)/))
+            end if
+
+            if (vNam(varn) .EQ. "z_i_wpt" .OR. vNam(varn) .EQ. "z_i_pt") then
+               if (mini(varn) .EQ. maxi(varn)) then
+                  if (min(data(varn,:)) .EQ. 0)then
+                     mini(varn)= mini(varn)-1.
+                     maxi(varn)= maxi(varn)+1.
+                  end if
+                  if (min(data(varn,:)) .LT. 0)then
+                     mini(varn)= mini(varn)-1.+(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
+                  end if
+                  if (min(data(varn,:)) .GT. 0)then
+                     mini(varn)= mini(varn)-1.-(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
+                  end if
+               end if
+               minz=min((/minz,mini(varn)/))
+               maxz=max((/maxz,maxi(varn)/))
+            end if
+
+            if (vNam(varn) .EQ. "wpptp0" .OR. vNam(varn) .EQ. "wpptp" .OR. vNam(varn) .EQ. "wpt") then
+               if (mini(varn) .EQ. maxi(varn)) then
+                  if (min(data(varn,:)) .EQ. 0)then
+                     mini(varn)= mini(varn)-1.
+                     maxi(varn)= maxi(varn)+1.
+                  end if
+                  if (min(data(varn,:)) .LT. 0)then
+                     mini(varn)= mini(varn)-1.+(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
+                  end if
+                  if (min(data(varn,:)) .GT. 0)then
+                     mini(varn)= mini(varn)-1.-(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
+                  end if
+               end if
+               minw=min((/minw,mini(varn)/))
+               maxw=max((/maxw,maxi(varn)/))
+            end if
+
+            if (vNam(varn) .EQ. "pt_0_" .OR. vNam(varn) .EQ. "pt_zp_") then
+               if (mini(varn) .EQ. maxi(varn)) then
+                  if (min(data(varn,:)) .EQ. 0)then
+                     mini(varn)= mini(varn)-1.
+                     maxi(varn)= maxi(varn)+1.
+                  end if
+                  if (min(data(varn,:)) .LT. 0)then
+                     mini(varn)= mini(varn)-1.+(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
+                  end if
+                  if (min(data(varn,:)) .GT. 0)then
+                     mini(varn)= mini(varn)-1.-(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
+                  end if
+               end if
+               minp=min((/minp,mini(varn)/))
+               maxp=max((/maxp,maxi(varn)/))
+            end if
+
+            if (vNam(varn) .EQ. "splptx" .OR. vNam(varn) .EQ. "splpty" .OR. vNam(varn) .EQ. "splptz") then
+               if (mini(varn) .EQ. maxi(varn)) then
+                  if (min(data(varn,:)) .EQ. 0)then
+                     mini(varn)= mini(varn)-1.
+                     maxi(varn)= maxi(varn)+1.
+                  end if
+                  if (min(data(varn,:)) .LT. 0)then
+                     mini(varn)= mini(varn)-1.+(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.-(maxi(varn))/2
+                  end if
+                  if (min(data(varn,:)) .GT. 0)then
+                     mini(varn)= mini(varn)-1.-(mini(varn))/2
+                     maxi(varn)= maxi(varn)+1.+(maxi(varn))/2
+                  end if
+               end if
+               mins=min((/mins,mini(varn)/))
+               maxs=max((/maxs,maxi(varn)/))
+            end if
+            
+         end if
+      end if
+   end do
+
+   do varn = dim-1,0,1 
+      
       if( isStrSubset (vNam(varn), "time") )
          check = False
-      end if
+      else
+         check = True
+      end if     
+      if( isvar("var") ) then     		
+         check = isStrSubset( var,","+vNam(varn)+"," )
+      end if
+      if (parameter(21) .NE. "variables") then
+         var = parameter(21) 
+         check = isStrSubset( var,","+vNam(varn)+"," ) 
+      end if 
   
       if(check) then
+
+        if (over .EQ. 1) then
  
-         n=n+1
-         data = f ->$vNam(varn)$          
-         print("plot of " + vNam(varn))
-
-         res@gsnRightString       = vNam(varn)
-         res@tiXAxisString        = " time [s] "
-         res@tiXAxisFontHeightF   = 0.07
-         res@txFontHeightF        = 0.07
-         res@tiYAxisFontHeightF   = 0.07
-      
-         plot_ps(n) = gsn_csm_xy(wks_ps,t,data(start_time_step:end_time_step),res)
+            res@gsnLeftString       = "overlayed plot"
+            res@gsnRightString      = unit(varn)
+            res@tiXAxisString        = " time [s] "
+            res@tiYAxisString        = " " 
+            res@tiXAxisFontHeightF   = 0.07
+            res@txFontHeightF        = 0.07
+            res@tiYAxisFontHeightF   = 0.07
+            res@xyLabelMode          = "Custom"
+            
+
+            if (vNam(varn) .EQ. "E")
+               E=0
+               res@xyLineColors     = (/237/)
+               res@xyExplicitLabels =  "E"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 237
+               res@trYMaxF           = minE
+               res@trYMinF           = maxE
+               plot_E = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if
+            if (vNam(varn) .EQ. "Es")
+               Es=0
+               res@xyLineColors            = (/144/)
+               res@xyExplicitLabels =  "Es"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 144
+               plot_Es = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
+            end if 
+             
+            if (vNam(varn) .EQ. "us")
+               us=0
+               res@xyLineColors            = (/237/)
+               res@xyExplicitLabels =  "us"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 237
+               res@trYMaxF           = minus
+               res@trYMinF           = maxus
+               plot_us = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if             
+            if (vNam(varn) .EQ. "ws")
+               ws=0
+               res@xyLineColors            = (/144/)
+               res@xyExplicitLabels =  "ws"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 144
+               plot_ws = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if
+             
+            if (vNam(varn) .EQ. "umax")
+               u=0
+               res@xyLineColors            = (/237/)
+               res@xyExplicitLabels =  "umax"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 237
+               res@trYMaxF           = minu
+               res@trYMinF           = maxu
+               plot_umax = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if
+            if (vNam(varn) .EQ. "vmax")
+               v=0
+               res@xyLineColors            = (/144/)
+               res@xyExplicitLabels =  "vmax"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 144
+               plot_vmax = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
+            end if
+            if (vNam(varn) .EQ. "wmax")
+               w=0
+               res@xyLineColors            = (/80/)
+               res@xyExplicitLabels =  "wmax"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 80
+               plot_wmax = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if 
    
-      end if
+            if (vNam(varn) .EQ. "z_i_wpt")
+               zw=0
+               res@xyLineColors            = (/237/) 
+               res@xyExplicitLabels =  "z_i_wpt"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 237
+               res@trYMaxF           = minz
+               res@trYMinF           = maxz 
+               plot_z_i_wpt = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if
+            if (vNam(varn) .EQ. "z_i_pt")
+               z=0
+               res@xyLineColors            = (/144/)  
+               res@xyExplicitLabels =  "z_i_pt"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 144 
+               plot_z_i_pt = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
+            end if
+
+            if (vNam(varn) .EQ. "wpptp0")
+               w0=0
+               res@xyLineColors            = (/237/)
+               res@xyExplicitLabels =  "wpptp0"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 237
+               res@trYMaxF           = minw
+               res@trYMinF           = maxw
+               plot_wpptp0 = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if
+            if (vNam(varn) .EQ. "wpptp")
+               wp=0
+               res@xyLineColors            = (/144/)
+               res@xyExplicitLabels =  "wpptp"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 144
+               plot_wpptp = gsn_csm_xy(wks_ps,t,data(varn,:),res)  
+            end if
+            if (vNam(varn) .EQ. "wpt")
+               wt=0
+               res@xyLineColors            = (/80/)
+               res@xyExplicitLabels =  "wpt"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 80
+               plot_wpt = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if
+
+            if (vNam(varn) .EQ. "pt_0_")
+               p=0
+               res@xyLineColors            = (/237/)
+               res@xyExplicitLabels =  "pt_0_"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 237
+               res@trYMaxF           = minp
+               res@trYMinF           = maxp
+               plot_pt_0_ = gsn_csm_xy(wks_ps,t,data(varn,:),res)  
+            end if
+            if (vNam(varn) .EQ. "pt_zp_")
+               pz=0
+               res@xyLineColors            = (/144/)
+               res@xyExplicitLabels =  "pt_zp_"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 144
+               plot_pt_zp_ = gsn_csm_xy(wks_ps,t,data(varn,:),res) 
+            end if
+
+            if (vNam(varn) .EQ. "splptx")
+               x=0
+               res@xyLineColors            = (/237/)
+               res@xyExplicitLabels =  "splptx"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 237
+               res@trYMaxF           = mins
+               res@trYMinF           = maxs
+               plot_splptx = gsn_csm_xy(wks_ps,t,data(varn,:),res)    
+            end if
+            if (vNam(varn) .EQ. "splpty")
+               y=0
+               res@xyLineColors            = (/144/)
+               res@xyExplicitLabels =  "splpty"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 144
+               plot_splpty = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+            end if
+            if (vNam(varn) .EQ. "splptz")
+               z=0
+               res@xyLineColors            = (/80/)
+               res@xyExplicitLabels =  "splptz"
+               res@xyLineLabelFontHeightF = 0.05 
+               res@xyLineLabelFontColor   = 80
+               plot_splptz = gsn_csm_xy(wks_ps,t,data(varn,:),res)         
+            end if
+
+         end if
+      end if
+   end do
+   
+   do varn = dim-1,0,1
  
+      if( isStrSubset (vNam(varn), "time") )
+         check = False
+      else 
+         check = True
+      end if    
+      if( isvar("var") ) then     		
+         check = isStrSubset( var,","+vNam(varn)+"," )
+      end if
+      if (parameter(21) .NE. "variables") then
+         var = parameter(21) 
+         check = isStrSubset( var,","+vNam(varn)+"," ) 
+      end if 
+    
+      if(check) then
+        
+         if (over .EQ. 1) then       
+   	
+            if (vNam(varn) .EQ. "E" .AND. Es .NE. 1) then 
+               E=1   
+               overlay(plot_E,plot_Es)
+               n=n+1
+               plot_ps(n) = plot_E               
+            end if
+            if (vNam(varn) .EQ. "Es" .AND. E .NE. 1) then
+               Es=1
+               overlay(plot_E,plot_Es)
+               n=n+1
+               plot_ps(n) = plot_E           
+            end if
+
+            if (vNam(varn) .EQ. "us" .AND. ws .NE. 1) then
+               us=1
+               overlay(plot_us,plot_ws)
+               n=n+1
+               plot_ps(n) = plot_us
+            end if
+            if (vNam(varn) .EQ. "ws" .AND. us .NE. 1) then
+               ws=1
+               overlay(plot_us,plot_ws)
+               n=n+1
+               plot_ps(n) = plot_us
+            end if
+         
+            if (vNam(varn) .EQ. "umax" .AND. v .NE. 1)
+               if (w .NE. 1) then 
+                  u=1          
+                  overlay(plot_umax,plot_vmax)
+                  overlay(plot_umax,plot_wmax)
+                  n=n+1
+                  plot_ps(n) = plot_umax
+               end if
+            end if
+            if (vNam(varn) .EQ. "vmax" .AND. u .NE. 1)
+               if (w .NE. 1) then
+                  v=1  
+                  overlay(plot_umax,plot_vmax)
+                  overlay(plot_umax,plot_wmax)
+                  n=n+1
+                  plot_ps(n) = plot_umax
+               end if
+            end if
+            if (vNam(varn) .EQ. "wmax" .AND. v .NE. 1)
+               if(u .NE. 1) then  
+                  w=1        
+                  overlay(plot_umax,plot_vmax)
+                  overlay(plot_umax,plot_wmax)
+                  n=n+1
+                  plot_ps(n) = plot_umax
+               end if
+            end if
+            
+            if (vNam(varn) .EQ. "z_i_wpt" .AND. z .NE. 1) then
+               zw=1       
+               overlay(plot_z_i_wpt,plot_z_i_pt)
+               n=n+1
+               plot_ps(n) = plot_z_i_wpt           
+            end if 
+            if (vNam(varn) .EQ. "z_i_pt" .AND. zw .NE. 1) then
+               z=1      
+               overlay(plot_z_i_wpt,plot_z_i_pt)
+               n=n+1
+               plot_ps(n) = plot_z_i_wpt           
+            end if    
+         
+            if (vNam(varn) .EQ. "wpptp0" .AND. wp .NE. 1)
+               if (wt .NE. 1) then
+                  w0=1
+                  overlay(plot_wpptp0,plot_wpptp)
+                  overlay(plot_wpptp0,plot_wpt) 
+                  n=n+1
+                  plot_ps(n) = plot_wpptp0 
+               end if            
+            end if
+            if (vNam(varn) .EQ. "wpptp" .AND. w0 .NE. 1)
+               if (wt .NE. 1) then
+                 wp=1
+                 overlay(plot_wpptp0,plot_wpptp)
+                 overlay(plot_wpptp0,plot_wpt) 
+                 n=n+1
+                 plot_ps(n) = plot_wpptp0  
+               end if           
+            end if
+            if (vNam(varn) .EQ. "wpt" .AND. wp .NE. 1)
+               if (w0 .NE. 1) then
+                  wt=1
+                  overlay(plot_wpptp0,plot_wpptp)
+                  overlay(plot_wpptp0,plot_wpt) 
+                  n=n+1
+                  plot_ps(n) = plot_wpptp0             
+               end if
+            end if
+
+            if (vNam(varn) .EQ. "pt_0_" .AND. pz .NE. 1) then
+               p=1      
+               overlay(plot_pt_0_,plot_pt_zp_)
+               n=n+1 
+               plot_ps(n) = plot_pt_0_           
+            end if 
+            if (vNam(varn) .EQ. "pt_zp_" .AND. p .NE. 1) then
+               pz=1        
+               overlay(plot_pt_0_,plot_pt_zp_)
+               n=n+1 
+               plot_ps(n) = plot_pt_0_           
+            end if
+            
+            if (vNam(varn) .EQ. "splptx" .AND. y .NE. 1)
+               if (z .NE.1 ) then
+                  x=1      
+                  overlay(plot_splptx,plot_splpty)
+                  overlay(plot_splptx,plot_splptz)
+                  n=n+1
+                  plot_ps(n) = plot_splptx           
+               end if
+            end if 
+            if (vNam(varn) .EQ. "splpty" .AND. x .NE. 1)
+               if(z .NE.1 ) then
+                  y=1         
+                  overlay(plot_splptx,plot_splpty)
+                  overlay(plot_splptx,plot_splptz)
+                  n=n+1
+                  plot_ps(n) = plot_splptx
+               end if           
+            end if 
+            if (vNam(varn) .EQ. "splptz" .AND. y .NE. 1)
+               if(x .NE.1 ) then
+                  z=1          
+                  overlay(plot_splptx,plot_splpty)
+                  overlay(plot_splptx,plot_splptz)
+                  n=n+1
+                  plot_ps(n) = plot_splptx    
+               end if       
+            end if
+
+            if(vNam(varn) .NE. "splptz" .AND. vNam(varn) .NE. "splpty" .AND. vNam(varn) .NE. "splptx" .AND. vNam(varn) .NE. "pt_zp_" .AND. vNam(varn) .NE. "pt_0_" .AND. vNam(varn) .NE. "wpt" .AND. vNam(varn) .NE. "wpptp" .AND. vNam(varn) .NE. "wpptp0" .AND. vNam(varn) .NE. "z_i_pt" .AND. vNam(varn) .NE. "z_i_wpt" .AND. vNam(varn) .NE. "wmax" .AND. vNam(varn) .NE. "vmax" .AND. vNam(varn) .NE. "umax" .AND. vNam(varn) .NE. "ws" .AND.  vNam(varn) .NE. "us" .AND. vNam(varn) .NE. "Es" .AND. vNam(varn) .NE. "E") then
+               n=n+1
+               res@xyLineColors   = (/237/)
+               res@xyLabelMode    = False
+               res@gsnLeftString  = vNam(varn)
+               res@gsnRightString = unit(varn) 
+               res@trYMaxF        = min(data(varn,:))
+               res@trYMinF        = max(data(varn,:))
+               if (min(data(varn,:)) .EQ. max(data(varn,:))) then
+                  if (min(data(varn,:)) .EQ. 0)then
+                     res@trYMaxF = min(data(varn,:))-1.
+                     res@trYMinF = max(data(varn,:))+1.
+                  end if
+                  if (min(data(varn,:)) .LT. 0)then
+                     res@trYMaxF = min(data(varn,:))+(min(data(varn,:)))/2
+                     res@trYMinF = max(data(varn,:))-(max(data(varn,:)))/2
+                  end if
+                  if (min(data(varn,:)) .GT. 0)then
+                     res@trYMaxF = min(data(varn,:))-(min(data(varn,:)))/2
+                     res@trYMinF = max(data(varn,:))+(max(data(varn,:)))/2
+                  end if
+               end if
+               plot_ps(n) = gsn_csm_xy(wks_ps,t,data(varn,:),res)  
+            end if
+       
+         else
+            
+            print("plot of " + vNam(varn)) 
+            
+            n=n+1
+            res@xyLineColors            = (/237/)
+            res@gsnLeftString       = vNam(varn)
+            res@gsnRightString      = unit(varn)
+            res@tiXAxisString        = " time [s] "
+            res@tiYAxisString        = " " 
+            res@tiXAxisFontHeightF   = 0.07
+            res@txFontHeightF        = 0.07
+            res@tiYAxisFontHeightF   = 0.07
+            res@trYMaxF        = min(data(varn,:))
+            res@trYMinF        = max(data(varn,:)) 
+            if (min(data(varn,:)) .EQ. max(data(varn,:))) then
+               if (min(data(varn,:)) .EQ. 0)then
+                     res@trYMaxF = min(data(varn,:))-1.
+                     res@trYMinF = max(data(varn,:))+1.
+                  end if
+                  if (min(data(varn,:)) .LT. 0)then
+                     res@trYMaxF = min(data(varn,:))+(min(data(varn,:)))/2
+                     res@trYMinF = max(data(varn,:))-(max(data(varn,:)))/2
+                  end if
+                  if (min(data(varn,:)) .GT. 0)then
+                     res@trYMaxF = min(data(varn,:))-(min(data(varn,:)))/2
+                     res@trYMinF = max(data(varn,:))+(max(data(varn,:)))/2
+                  end if
+            end if
+            plot_ps(n) = gsn_csm_xy(wks_ps,t,data(varn,:),res)
+           
+         end if     
+      end if
    end do
-
+  
    ; ***************************************************
    ; merge plots onto one page 
    ; ***************************************************
-
-   do np = 1,n,no_lines*no_columns
-    
-      if ( np + no_lines*no_columns .gt. n) then   
-         gsn_panel(wks_ps, plot_ps(np:n),(/no_lines,no_columns/),resP)
-      else
-         gsn_panel(wks_ps, plot_ps(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
-      end if
-     
-   end do
+  
+   if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
+      gsn_panel(wks_ps,plot_ps(1:n),(/n,1/),resP)
+   else
+      do np = 1,n,no_lines*no_columns    
+         if ( np + no_lines*no_columns .gt. n) then   
+            gsn_panel(wks_ps, plot_ps(np:n),(/no_lines,no_columns/),resP)
+         else
+            gsn_panel(wks_ps, plot_ps(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),resP)
+         end if
+      end do
+   end if
 
    print(" ")
