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(" ")
