Index: palm/trunk/SCRIPTS/NCL/cross_sections.ncl
===================================================================
--- palm/trunk/SCRIPTS/NCL/cross_sections.ncl	(revision 161)
+++ palm/trunk/SCRIPTS/NCL/cross_sections.ncl	(revision 162)
@@ -61,25 +61,25 @@
    if ( .not. isvar("sort") ) then   			; sort of plots
       sort = "time"
-      if (parameter(37) .NE. "time") then
-         sort = parameter(37)  
+      if (parameter(51) .NE. "time") then
+         sort = parameter(51)  
       end if
    end if 
    if ( .not. isvar("mode") ) then   			; pattern of contour plots
       mode = "Fill"
-      if (parameter(39) .NE. "Fill") then
-         mode = parameter(39)  
+      if (parameter(49) .NE. "Fill") then
+         mode = parameter(49)  
       end if
    end if
    if ( .not. isvar("fill_mode") ) then 		; pattern of filling 
       fill_mode = "AreaFill"
-      if (parameter(41) .NE. "AreaFill") then
-         mode = parameter(41)  
+      if (parameter(53) .NE. "AreaFill") then
+         mode = parameter(53)  
       end if 
    end if
    if ( .not. isvar("shape") ) then   			; keeping of aspect ratio
       shape = 1
-      if (parameter(43) .NE. "1") then
-         shape = stringtointeger(parameter(43)) 
-         if (stringtointeger(parameter(43)) .NE. 0) then
+      if (parameter(55) .NE. "1") then
+         shape = stringtointeger(parameter(55)) 
+         if (stringtointeger(parameter(55)) .NE. 0) then
             print(" ")
             print("Please set 'shape' to 0 or 1")
@@ -95,7 +95,7 @@
    if ( .not. isvar("xyc") ) then   			; turn xy cross-section on or off
       xyc = 0
-      if (parameter(45) .NE. "0") then
-         xyc = stringtointeger(parameter(45)) 
-         if (stringtointeger(parameter(45)) .NE. 1) then
+      if (parameter(57) .NE. "0") then
+         xyc = stringtointeger(parameter(57)) 
+         if (stringtointeger(parameter(57)) .NE. 1) then
             print(" ")
             print("Please set 'xyc' to 0 or 1")
@@ -107,7 +107,7 @@
    if ( .not. isvar("xzc") ) then   			; turn xz cross-section on or off 
       xzc = 0
-      if (parameter(47) .NE. "0") then
-         xzc = stringtointeger(parameter(47)) 
-         if (stringtointeger(parameter(47)) .NE. 1) then
+      if (parameter(59) .NE. "0") then
+         xzc = stringtointeger(parameter(59)) 
+         if (stringtointeger(parameter(59)) .NE. 1) then
             print(" ")
             print("Please set 'xzc' to 0 or 1")
@@ -119,7 +119,7 @@
    if ( .not. isvar("yzc") ) then   			; turn yz cross-section on or off 
       yzc = 0
-      if (parameter(49) .NE. "0") then
-         yzc = stringtointeger(parameter(49)) 
-         if (stringtointeger(parameter(49)) .NE. 1) then
+      if (parameter(61) .NE. "0") then
+         yzc = stringtointeger(parameter(61)) 
+         if (stringtointeger(parameter(61)) .NE. 1) then
             print(" ")
             print("Please set 'yzc' to 0 or 1")
@@ -129,22 +129,4 @@
       end if
    end if
-   if ( .not. isvar("xs") ) then     			; start of x-coordinate range 
-      xs = 0
-      if (parameter(51) .NE. "0") then
-         xs = stringtointeger(parameter(51))
-      end if
-   end if 
-   if ( .not. isvar("ys") ) then     			; start of y-coordinate range
-      ys = 0
-      if (parameter(55) .NE. "0") then
-         ys = stringtointeger(parameter(55))
-      end if
-   end if
-   if ( .not. isvar("zs") ) then     			; start of z-coordinate range
-      zs = 0
-      if (parameter(59) .NE. "0") then
-         zs = stringtointeger(parameter(59))
-      end if
-   end if  
 
    if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
@@ -180,7 +162,7 @@
    if ( .not. isvar("vector") ) then   			; sort of plots
       vector = 0
-      if (parameter(63) .NE. "0") then
-         vector = stringtointeger(parameter(63)) 
-         if (stringtointeger(parameter(63)) .NE. 1) then
+      if (parameter(39) .NE. "0") then
+         vector = stringtointeger(parameter(39)) 
+         if (stringtointeger(parameter(39)) .NE. 1) then
             print(" ")
             print("Please set 'vector' to 0 or 1")
@@ -192,6 +174,6 @@
    if ( .not. isvar("ref_mag") ) then   			; sort of plots
       ref_mag = 0.05
-      if (parameter(71) .NE. "0.05") then
-         ref_mag = stringtofloat(parameter(71))   
+      if (parameter(47) .NE. "0.05") then
+         ref_mag = stringtofloat(parameter(47))   
       end if
    end if 
@@ -208,8 +190,77 @@
    print(" ")
    dim   = dimsizes(vNam)
+
+   ; ***************************************************
+   ; check for kind of input file
+   ; ***************************************************
+
+   check_3d = True
+   do varn=0,dim-1
+      checkxy = isStrSubset( vNam(varn),"xy")
+      checkxz = isStrSubset( vNam(varn),"xz")
+      checkyz = isStrSubset( vNam(varn),"yz")
+      if (checkxy .OR. checkxz .OR. checkyz) then
+         check_3d = False
+         break
+      end if
+   end do
+   if (.not. check_3d) then
+      if (xyc .EQ. 1 .AND. .not. checkxy) then
+         print(" ")
+         print("Your input file doesn't have values for xy cross-sections")
+         if (checkxz)then
+            print("Select another input file or xz cross-sections")
+         else
+            print("Select another input file or yz cross-sections")
+         end if
+         print(" ")
+         exit
+      else
+         print(" ")
+         print("Your input file contains xy data")
+         print(" ") 
+      end if
+      if (xzc .EQ. 1 .AND. .not. checkxz) then
+         print(" ")
+         print("Your input file doesn't have values for xz cross-sections")
+         if (checkxy)then
+            print("Select another input file or xy cross-sections")
+         else
+            print("Select another input file or yz cross-sections")
+         end if
+         print(" ")
+         exit
+      else
+         print(" ")
+         print("Your input file contains xz data")
+         print(" ")
+      end if
+      if (yzc .EQ. 1 .AND. .not. checkyz) then
+         print(" ")
+         print("Your input file doesn't have values for yz cross-sections")
+         if (checkxy)then
+            print("Select another input file or xy cross-sections")
+         else
+            print("Select another input file or xz cross-sections")
+         end if
+         print(" ")
+         exit
+      else
+         print(" ")
+         print("Your input file contains yz data")
+         print(" ")
+      end if
+   else
+      print(" ")
+      print("Your input file: '"+file_in+"'")
+      print("contains 3d or other data")
+      print(" ")
+   end if
+   
    if (dim .EQ. 0) then
       print(" ")
-      print("There are no data on file")
+      print("There is no data on file")
       print(" ")
+      exit
    end if
 
@@ -241,8 +292,11 @@
    cs_res@lgLabelFontHeightF     = .02
    cs_res@txFontHeightF          = .02
+   cs_res@cnLevelSelectionMode    = "ManualLevels"
+
    cs_resP = True
    cs_resP@txString               = f@title
    cs_resP@txFuncCode             = "~"         	; necessary for title
    cs_resP@txFontHeightF          = .02
+   
   
    if ( mode .eq. "Fill" ) then
@@ -268,44 +322,59 @@
    ; *********************************************
 
-   t    = f->time
-   nt = dimsizes(t)
-     
-   ; *********************************************	
+   t_all = f->time
+   nt = dimsizes(t_all)
+   delta_t = t_all(nt-1)/nt
+
+   ; ****************************************************	
    ; start of time step and different types of mistakes that could be done
-   ; *********************************************
+   ; ****************************************************
 
    if ( .not. isvar("start_time_step") ) then   	
-      start_time_step = 0
-      if (parameter(13) .NE. "1") then
-         if (parameter(13) .LT. "1")
-            print(" ")
-            print("Begin with time step 1")
+      start_time_step=t_all(0)/3600
+      if (parameter(13) .NE. "t(0)") then
+         if (stringtodouble(parameter(13)) .GT. t_all(nt-1)/3600)
+            print(" ")
+            print("'start_time_step' = "+ parameter(13) +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
+            print(" ")
+            print("Please select another 'start_time_step'")
             print(" ")
             exit
          end if
-         if (stringtointeger(parameter(13)) .GT. nt)
-            print(" ")
-            print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt))
+         if (stringtofloat(parameter(13)) .LT. t_all(0)/3600)
+            print(" ")
+            print("'start_time_step' = "+ parameter(13) +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
+            print(" ")
+            print("Please select another 'start_time_step'") 
             print(" ")
             exit
          end if
-         start_time_step = stringtointeger(parameter(13))-1  
+         start_time_step=stringtodouble(parameter(13)) 
       end if
    else 
-      if (start_time_step .LE. 0)
-         print(" ")
-         print("Begin with time step 1")
+      if (start_time_step .GT. t_all(nt-1)/3600)
+         print(" ")
+         print("'start_time_step' = "+ start_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
+         print(" ")
+         print("Please select another 'start_time_step'")
          print(" ")
          exit
       end if
-      if (start_time_step .GT. nt)
-         print(" ")
-         print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt)) 
+      if (start_time_step .LT. t_all(0)/3600)
+         print(" ")
+         print("'start_time_step' = "+ start_time_step +"h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
+         print(" ")
+         print("Please select another 'start_time_step'")
          print(" ")
          exit
       end if
-      start_time_step = start_time_step - 1
-   end if
-
+   end if
+   start_time_step = start_time_step*3600
+   do i=0,nt-1     
+      if (start_time_step .GE. t_all(i)-delta_t/2 .AND. start_time_step .LT. t_all(i)+delta_t/2)then
+         st=i
+         break
+      end if
+   end do
+        
    ; ****************************************************
    ; end of time step and different types of mistakes that could be done
@@ -313,48 +382,76 @@
 
    if ( .not. isvar("end_time_step") ) then     	
-      end_time_step = nt-1
-      if (parameter(15) .NE. "nt") then
-         if (parameter(15) .LE. "0")
-            print(" ")
-            print("'end_time_step' = "+parameter(15)+ " is too small; 'end_time_step' should be at least 1 ")
+      end_time_step = t_all(nt-1)/3600
+      if (parameter(15) .NE. "t(end)") then
+         if (stringtodouble(parameter(15)) .LT. t_all(0)/3600)
+            print(" ")
+            print("'end_time_step' = "+parameter(15)+ "h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
+            print(" ")
+            print("Please select another 'end_time_step'")
             print(" ")
             exit
          end if
-         if (stringtointeger(parameter(15)) .GT. nt)
-            print(" ")
-            print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt))
+         if (stringtodouble(parameter(15)) .GT. t_all(nt-1)/3600)
+            print(" ")
+            print("'end_time_step' = "+ parameter(15) +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
+            print(" ")
+            print("Please select another 'end_time_step'")  
             print(" ")
             exit
          end if
-         if (stringtointeger(parameter(15)) .LT. stringtointeger(parameter(13)) )
-            print(" ")
-            print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(15))
+         if (stringtodouble(parameter(15)) .LT. start_time_step/3600)
+            print(" ")
+            print("'end_time_step' = "+ parameter(15) +"h is lower than 'start_time_step' = "+parameter(13)+"h")
+            print(" ")
+            print("Please select another 'start_time_step' or 'end_time_step'")
             print(" ")
             exit
          end if
-         end_time_step = stringtointeger(parameter(15))-1  
+         end_time_step = stringtodouble(parameter(15))
       end if    
    else
-      if (end_time_step .LE. 0)
-         print(" ")
-         print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
+      if (end_time_step .LT. t_all(0)/3600)
+         print(" ")
+         print("'end_time_step' = "+end_time_step+ "h is lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
+         print(" ")
+         print("Please select another 'end_time_step'")
          print(" ")
          exit
       end if
-      if (end_time_step .GT. nt)
-         print(" ")
-         print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt))
-         print(" ") 
+      if (end_time_step .GT. t_all(nt-1)/3600)
+         print(" ")
+         print("'end_time_step' = "+ end_time_step +"h is greater than last time step = " + t_all(nt-1)+"s = "+t_all(nt-1)/3600+"h")
+         print(" ")
+         print("Please select another 'end_time_step'")  
+         print(" ")
          exit
       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)
+      if (end_time_step .LT. start_time_step/3600)
+         print(" ")
+         print("'end_time_step' = "+ end_time_step +"h is lower than 'start_time_step' = "+start_time_step+"h")
+         print(" ")
+         print("Please select another 'start_time_step' or 'end_time_step'")
          print(" ")
          exit
       end if
-      end_time_step = end_time_step-1 
-   end if
-
+   end if
+   end_time_step = end_time_step*3600
+   do i=0,nt-1     
+      if (end_time_step .GE. t_all(i)-delta_t/2 .AND. end_time_step .LT. t_all(i)+delta_t/2)then
+         et=i
+         break
+      end if
+   end do 
+  
+   delete(start_time_step)
+   start_time_step=round(st,3)
+   delete(end_time_step)
+   end_time_step=round(et,3)
+
+   print(" ")
+   print("Output from "+t_all(start_time_step)/3600+"h = "+t_all(start_time_step)+"s => index = "+start_time_step)
+   print("       till "+t_all(end_time_step)/3600+"h = "+t_all(end_time_step)+"s => index = "+end_time_step)
+   print(" ")
+  
    no_time=(end_time_step-start_time_step)+1
 
@@ -365,9 +462,9 @@
    if (vector .EQ. 1) then
       if ( .not. isvar("plotvec") ) then
-         plotvec = parameter(69)
+         plotvec = parameter(45)
       end if
       if ( .not. isvar("vec1") ) then
-         vec1 = parameter(65)
-         if (parameter(65) .EQ. "vec1") then
+         vec1 = parameter(41)
+         if (parameter(41) .EQ. "vec1") then
             print(" ")
             print("Please indicate Vector 1 ('vec1') for Vector-Plot")
@@ -377,6 +474,6 @@
       end if
       if ( .not. isvar("vec2") ) then
-         vec2 = parameter(67)
-         if (parameter(67) .EQ. "vec2") then
+         vec2 = parameter(43)
+         if (parameter(43) .EQ. "vec2") then
             print(" ")
             print("Please indicate Vector 2 ('vec2') for Vector-Plot")
@@ -395,141 +492,692 @@
    ; ****************************************************
 
-   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
+   if (xyc .EQ. 1) then
+      do varn=0,dim-1
+         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x")then    
+            x_d     = f->$vNam(varn)$
+            xdim    = dimsizes(x_d)-1
+            delta_x = x_d(1)-x_d(0)
+            break
+         end if
+      end do 
+      do varn=0,dim-1
+         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then   
+            y_d     = f->$vNam(varn)$
+            ydim    = dimsizes(y_d)-1
+            delta_y = y_d(1)-y_d(0)
+            break
+         end if
+      end do
+      do varn=0,dim-1
+         if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then
+            z_d     = f->$vNam(varn)$
+            zdim    = dimsizes(z_d)-1
+            delta_z = 0
+            break
+         else
+            if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then
+               z_d     = f->$vNam(varn)$
+               zdim    = dimsizes(z_d)-1 
+               delta_z = -1.d
+               break
+            else
+               zdim    = 0
+               delta_z = -1.d
+            end if
+         end if
+      end do
+   end if
+   if (xzc .EQ. 1) then
+      do varn=0,dim-1
+         if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x") then 
+            x_d     = f->$vNam(varn)$
+            xdim    = dimsizes(x_d)-1
+            delta_x = x_d(1)-x_d(0)
+            break
+         end if
+      end do
+      do varn=0,dim-1   
+         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
+            z_d     = f->$vNam(varn)$
+            zdim    = dimsizes(z_d)-1 
+            delta_z = z_d(1)-z_d(0)
+            break
+         end if
+      end do
+      do varn=0,dim-1 
+         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
+            y_d     = f->$vNam(varn)$
+            ydim    = dimsizes(y_d)-1
+            delta_y = y_d(1)-y_d(0)
+            break
+         else
+            if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
+               y_d     = f->$vNam(varn)$
+               ydim    = dimsizes(y_d)-1 
+               delta_y = -1.d
+               break
+            else
+               ydim    = 0
+               delta_y = -1.d
+            end if
+         end if 
+      end do
+   end if
+   if (yzc .EQ. 1) then
+      do varn=0,dim-1  
+         if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
+            y_d     = f->$vNam(varn)$
+            ydim    = dimsizes(y_d)-1
+            delta_y = y_d(1)-y_d(0)
+            break
+         end if
+      end do
+      do varn=0,dim-1 
+         if (vNam(varn) .eq. "zw" .or. vNam(varn) .eq. "zu" .OR. vNam(varn) .eq. "zu_3d" .or. vNam(varn) .eq. "zw_3d") then
+            z_d     = f->$vNam(varn)$
+            zdim    = dimsizes(z_d)-1 
+            delta_z = z_d(1)-z_d(0)
+            break
+         end if
+      end do
+      do varn=0,dim-1 
+         if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
+            x_d     = f->$vNam(varn)$
+            xdim    = dimsizes(x_d)-1 
+            delta_x = x_d(1)-x_d(0)
+            break
+         else
+            if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
+               x_d     = f->$vNam(varn)$
+               xdim    = dimsizes(x_d)-1 
+               delta_x = -1.d
+               break
+            else
+               xdim    = 0
+               delta_x = -1.d
+            end if
+         end if 
+      end do
+   end if
+   
+   ; ****************************************************
+   ; set up ranges of x-, y- and z-coordinates
+   ; ****************************************************         
+                    
+   if ( .not. isvar("xs") ) then     		
+      xs     = 0.0d
+      xstart = 0
+      if (parameter(63) .NE. "x0") then
+         if (delta_x .EQ. -1) then
+            print(" ") 
+            print("You cannot choose a start value for x, there are preseted layers for x")
+            print(" ")
+            xstart=0
+         else
+            if (stringtodouble(parameter(63)) .LT. 0-delta_x/2) then
+               print(" ")
+               print("range start for x-coordinate = "+parameter(63)+"m is lower than first value x = "+0+"m or xu = "+(0-delta_x/2)+"m")
+               print(" ")
+               exit
+            end if
+            if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
+               if (stringtodouble(parameter(63)) .GE. xdim*delta_x) then
+                  print(" ")
+                  print("range start for x-coordinate = "+parameter(63)+"m is equal or greater than last value = "+xdim*delta_x+"m")
+                  print(" ")
+                  exit
+               end if 
+            else
+               if (stringtodouble(parameter(63)) .GT. xdim*delta_x) then
+                  print(" ")
+                  print("range start for x-coordinate = "+parameter(63)+"m is greater than last value = "+xdim*delta_x+"m")
+                  print(" ")
+                  exit
+               end if 
+            end if
+            xs = stringtodouble(parameter(63))
+         end if
+      end if
+   else
+      if (delta_x .EQ. -1) then
+         print(" ") 
+         print("You cannot choose a start value for x, there are preseted layers for x")
+         print(" ")
+         xstart=0
       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
-         print(vNam(varn))
-         data_all = f->$vNam(varn)$
-
-         ; ****************************************************
-         ; set up ranges of x-, y- and z-coordinates
-         ; ****************************************************         
-                    
-         xdim = dimsizes(data_all(0,0,0,:)) - 1    
-         ydim = dimsizes(data_all(0,0,:,0)) - 1        
-         zdim = dimsizes(data_all(0,:,0,0)) - 1  
-       
-         if ( .not. isvar("xe")) then     ; output x-coordinate range end (in index)
-            xe = xdim
-            if (parameter(53) .NE. "xdim") then
-               if (stringtointeger(parameter(53)) .GT. xdim) then
-                  print(" ")
-                  print("range end for x-coordinate = "+parameter(53)+" is higher than available dimensions = "+xdim)
-                  print(" ") 
+         if (xs .LT. 0-delta_x/2) then
+            print(" ")
+            print("range start for x-coordinate = "+xs+"m is lower than first value = "+(0-delta_x/2)+"m")
+            print(" ")
+            exit
+         end if
+         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
+            if (xs .GE. xdim*delta_x) then
+               print(" ")
+               print("range start for x-coordinate = "+xs+"m is equal or greater than last value = "+xdim*delta_x+"m")
+               print(" ")
+               exit
+            end if 
+         else
+            if (xs .GT. xdim*delta_x) then
+               print(" ")
+               print("range start for x-coordinate = "+xs+"m is greater than last value = "+xdim*delta_x+"m")
+               print(" ")
+               exit
+            end if 
+         end if
+      end if
+   end if 
+
+   if ( .not. isvar("ys") ) then    
+      ys = 0.0d
+      ystart = 0
+      if (parameter(67) .NE. "y0") then
+         if (delta_y .EQ. -1) then
+            print(" ") 
+            print("You cannot choose a start value for y, there are preseted layers for y")
+            print(" ")
+            ystart=0
+         else
+            if (stringtodouble(parameter(67)) .LT. 0-delta_y/2) then
+               print(" ")
+               print("range start for y-coordinate = "+parameter(67)+"m is lower than first value = "+0-delta_y/2+"m")
+               print(" ")
+               exit
+            end if
+            if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
+               if (stringtodouble(parameter(67)) .GE. ydim*delta_y) then
+                  print(" ")
+                  print("range start for y-coordinate = "+parameter(67)+"m is equal or greater than last value = "+ydim*delta_y+"m")
+                  print(" ")
                   exit
                end if
-               if (stringtointeger(parameter(53)) .LT. 0 .OR. stringtointeger(parameter(53)) .LT. xs) then
-                  print(" ")
-                  print("range end for x-coordinate = "+parameter(53)+" is too small")
+            else
+               if (stringtodouble(parameter(67)) .GT. ydim*delta_y) then
+                  print(" ")
+                  print("range start for y-coordinate = "+parameter(67)+"m is greater than last value = "+ydim*delta_y+"m")
                   print(" ")
                   exit
                end if
-               xe = stringtointeger(parameter(53)) 
-            end if 
-         else
-            if (xe .GT. xdim) then
-               print(" ")
-               print("range end for x-coordinate = "+xe+" is higher than available dimensions = "+xdim)
-               print(" ")
-               exit
-            end if
-            if (xe .LT. 0 .OR. xe .LT. xs) then
-               print(" ")
-               print("range end for x-coordinate = "+xe+" is too small")
-               print(" ")
-               exit
-            end if                  
-         end if
-
-         if ( .not. isvar("ye")) then     ; output y-coordinate range end (in index)
-            ye = ydim
-            if (parameter(57) .NE. "ydim") then
-               if (stringtointeger(parameter(57)) .GT. ydim) then
-                  print(" ")
-                  print("range end for y-coordinate = "+parameter(57)+" is higher than available dimensions = "+ydim)
+            end if
+            ys = stringtodouble(parameter(67))
+         end if
+      end if
+   else
+      if (delta_y .EQ. -1) then
+         print(" ") 
+         print("You cannot choose a start value for y, there are preseted layers for y")
+         print(" ")
+         ystart=0
+      else
+         if (ys .LT. 0-delta_y/2) then
+            print(" ")
+            print("range start for y-coordinate = "+ys+"m is lower than first value = "+0-delta_y/2+"m")
+            print(" ")
+            exit
+         end if
+         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
+            if (ys .GE. ydim*delta_y) then
+               print(" ")
+               print("range start for y-coordinate = "+ys+"m is equal or greater than last value = "+ydim*delta_y+"m")
+               print(" ")
+               exit
+            end if
+         else
+            if (ys .GT. ydim*delta_y) then
+               print(" ")
+               print("range start for y-coordinate = "+ys+"m is greater than last value = "+ydim*delta_y+"m")
+               print(" ")
+               exit
+            end if
+         end if
+      end if
+   end if
+  
+   if ( .not. isvar("zs") ) then     			; start of z-coordinate range
+      zs = 0
+      if (parameter(71) .NE. "z0") then
+         if (delta_z .EQ. -1) then
+            print(" ") 
+            print("You cannot choose a start value for z, there are preseted layers for z")
+            print(" ")
+         else
+            print(" ")
+            print("Please mind to indicate start and end ranges for the z-coordinate in")
+            print("indices not in 'meters'. Corresponding index and meter:")
+            print(" ")
+            print(" = "+z_d+" m")
+            print(" ") 
+            if (stringtointeger(parameter(71)) .LT. 0) then
+               print(" ")
+               print("range start for z-coordinate = "+parameter(71)+" is lower than first gridpoint = 0")
+               print(" ")
+               exit
+            end if
+            if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
+               if (stringtointeger(parameter(71)) .GE. zdim) then
+                  print(" ")
+                  print("range start for z-coordinate = "+parameter(71)+" is equal or greater than last gridpoint = "+zdim)
                   print(" ")
                   exit
                end if
-               if (stringtointeger(parameter(57)) .LT. 0 .OR. stringtointeger(parameter(57)) .LT. ys) then
-                  print(" ")
-                  print("range end for y-coordinate = "+parameter(57)+" is too small")
+            else
+               if (stringtodouble(parameter(71)) .GT. zdim) then
+                  print(" ")
+                  print("range start for z-coordinate = "+parameter(71)+" is greater than last gridpoint = "+zdim)
                   print(" ")
                   exit
                end if
-               ye = stringtointeger(parameter(57)) 
+            end if
+            zs = stringtointeger(parameter(71))
+         end if
+      end if
+   else
+      if (delta_z .EQ. -1) then
+         print(" ") 
+         print("You cannot choose a start value for z, there are preseted layers for z")
+         print(" ")
+      else
+         if (zs .LT. 0) then
+            print(" ")
+            print("range start for z-coordinate = "+zs+" is lower than first gridpoint = 0")
+            print(" ")
+            exit
+         end if
+         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
+            if (zs .GE. zdim) then
+               print(" ")
+               print("range start for z-coordinate = "+zs+" is equal or greater than last gridpoint = "+zdim)
+               print(" ")
+               exit
+            end if
+         else
+            if (zs .GT. zdim) then
+               print(" ")
+               print("range start for z-coordinate = "+zs+" is greater than last gridpoint = "+zdim)
+               print(" ")
+               exit
+            end if
+         end if
+      end if
+   end if  
+
+   if ( .not. isvar("xe")) then    
+      xe   = xdim*delta_x
+      xend = xdim
+      if (parameter(65) .NE. "xdim") then
+         if (delta_x .EQ. -1) then
+            print(" ") 
+            print("You cannot choose an end value for x, there are preseted layers for x")
+            print(" ") 
+            xend=xdim           
+         else
+            if (stringtodouble(parameter(65)) .GT. xdim*delta_x) then
+               print(" ")
+               print("range end for x-coordinate = "+parameter(65)+"m is greater than last value = "+xdim*delta_x+"m")
+               print(" ") 
+               exit 
             end if 
-         else
-            if (ye .GT. ydim) then
-               print(" ")
-               print("range end for y-coordinate = "+ye+" is higher than available dimensions = "+ydim)
-               print(" ")
-               exit 
-            end if
-            if (ye .LT. 0 .OR. ye .LT. ys) then
-               print(" ")
-               print("range end for y-coordinate = "+ye+" is too small")
-               print(" ")
-               exit
-            end if
-         end if
-
-         if ( .not. isvar("ze")) then     ; output z-coordinate range end (in index)
-            ze = zdim
-            if (parameter(61) .NE. "zdim") then
-               if (stringtointeger(parameter(61)) .GT. zdim) then
-                  print(" ")
-                  print("range end for z-coordinate = "+parameter(61)+" is higher than available dimensions = "+zdim)
+            if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
+               if (stringtodouble(parameter(65)) .LE. 0-delta_x/2)
+                  print(" ")
+                  print("range end for x-coordinate = "+parameter(65)+"m is equal or lower than first value = "+(0-delta_x/2)+"m")
                   print(" ")
                   exit
                end if
-               if (stringtointeger(parameter(61)) .LT. 0 .OR. stringtointeger(parameter(61)) .LT. zs) then
-                  print(" ")
-                  print("range end for z-coordinate = "+parameter(61)+" is too small")
+               if (stringtodouble(parameter(65)) .LE. xs) then
+                  print(" ")
+                  print("range end for x-coordinate = "+parameter(65)+"m is equal or lower than start range = "+xs+"m")
                   print(" ")
                   exit
                end if
-               ze = stringtointeger(parameter(61)) 
+            else
+               if (stringtodouble(parameter(65)) .LT. 0-delta_x/2)
+                  print(" ")
+                  print("range end for x-coordinate = "+parameter(65)+"m is lower than first value = "+(0-delta_x/2)+"m")
+                  print(" ")
+                  exit
+               end if
+               if (stringtodouble(parameter(65)) .LT. xs) then
+                  print(" ")
+                  print("range end for x-coordinate = "+parameter(65)+"m is lower than start range = "+xs+"m")
+                  print(" ")
+                  exit
+               end if
+            end if         
+            xe = stringtodouble(parameter(65)) 
+         end if 
+      end if
+   else
+      if (delta_x .EQ. -1) then
+         print(" ") 
+         print("You cannot choose an end value for x, there are preseted layers for x")
+         print(" ")
+         xend=xdim
+      else
+         if (xe .GT. xdim*delta_x) then
+            print(" ")
+            print("range end for x-coordinate = "+xe+"m is greater than last value = "+xdim*delta_x+"m")
+            print(" ") 
+            exit 
+         end if 
+         if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
+            if (xe .LE. 0-delta_x/2)
+               print(" ")
+               print("range end for x-coordinate = "+xe+"m is equal or lower than first value = "+(0-delta_x/2)+"m")
+               print(" ")
+               exit
+            end if
+            if (xe .LE. xs) then
+               print(" ")
+               print("range end for x-coordinate = "+xe+"m is equal or lower than start range = "+xs+"m")
+               print(" ")
+               exit
+            end if
+            if (stringtodouble(xe .EQ. xs+1)) then
+               print(" ")
+               print("range end for x-coordinate = "+xe+"m must be at least two more gridpoints greater than start range = "+xs+"m")
+               print(" ")
+               exit
+            end if
+         else
+            if (xe .LT. 0-delta_x/2)
+               print(" ")
+               print("range end for x-coordinate = "+xe+"m is lower than first value = "+(0-delta_x/2)+"m")
+               print(" ")
+               exit
+            end if
+            if (xe .LT. xs) then
+               print(" ")
+               print("range end for x-coordinate = "+xe+"m is lower than start range = "+xs+"m")
+               print(" ")
+               exit
+            end if
+         end if 
+      end if               
+   end if
+
+   if ( .not. isvar("ye")) then  
+      ye   = ydim*delta_y
+      yend = ydim
+      if (parameter(69) .NE. "ydim") then
+         if (delta_y .EQ. -1) then
+            print(" ") 
+            print("You cannot choose an end value for y, there are preseted layers for y")
+            print(" ")
+            yend=ydim
+         else
+            if (stringtodouble(parameter(69)) .GT. ydim*delta_y) then
+               print(" ")
+               print("range end for y-coordinate = "+parameter(69)+"m is greater than last value = "+ydim*delta_y+"m")
+               print(" ") 
+               exit 
             end if 
-         else
-            if (ze .GT. zdim) then
-               print(" ")
-               print("range end for z-coordinate = "+ze+" is higher than available dimensions = "+zdim)
+            if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
+               if (stringtodouble(parameter(69)) .LE. 0-delta_y/2)
+                  print(" ")
+                  print("range end for y-coordinate = "+parameter(69)+"m is equal or lower than first value = "+(0-delta_y/2)+"m")
+                  print(" ")
+                  exit
+               end if
+               if (stringtodouble(parameter(69)) .LE. ys) then
+                  print(" ")
+                  print("range end for y-coordinate = "+parameter(69)+"m is equal or lower than start range = "+ys+"m")
+                  print(" ")
+                  exit
+               end if
+               if (stringtodouble(parameter(69)) .EQ. ys+1) then
+                  print(" ")
+                  print("range end for y-coordinate = "+parameter(69)+"m must be at least two more gridpoints greater than start range = "+ys+"m")
+                  print(" ")
+                  exit
+               end if 
+            else
+               if (stringtodouble(parameter(69)) .LT. 0-delta_y/2)
+                  print(" ")
+                  print("range end for y-coordinate = "+parameter(69)+"m is lower than first value = "+(0-delta_y/2)+"m")
+                  print(" ")
+                  exit
+               end if
+               if (stringtodouble(parameter(69)) .LT. ys) then
+                  print(" ")
+                  print("range end for y-coordinate = "+parameter(69)+"m is lower than start range = "+ys+"m")
+                  print(" ")
+                  exit
+               end if
+            end if         
+            ye = stringtodouble(parameter(69)) 
+         end if 
+      end if
+   else
+      if (delta_y .EQ. -1) then
+         print(" ") 
+         print("You cannot choose an end value for y, there are preseted layers for y")
+         print(" ")
+         yend=ydim
+      else
+         if (ye .GT. ydim*delta_y) then
+            print(" ")
+            print("range end for y-coordinate = "+ye+"m is greater than last value = "+ydim*delta_y+"m")
+            print(" ") 
+            exit 
+         end if 
+         if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
+            if (ye .LE. 0-delta_y/2)
+               print(" ")
+               print("range end for y-coordinate = "+ye+"m is equal or lower than first value = "+(0-delta_y/2)+"m")
+               print(" ")
+               exit
+            end if
+            if (ye .LE. ys) then
+               print(" ")
+               print("range end for y-coordinate = "+ye+"m is equal or lower than start range = "+ys+"m")
+               print(" ")
+               exit
+            end if
+            if (ye .EQ. ys+1) then
+               print(" ")
+               print("range end for y-coordinate = "+ye+"m must be at least two more gridpoints greater than start range = "+ys+"m")
+               print(" ")
+               exit
+            end if
+         else
+            if (ye .LT. 0-delta_y/2)
+               print(" ")
+               print("range end for y-coordinate = "+ye+"m is lower than first value = "+(0-delta_y/2)+"m")
+               print(" ")
+               exit
+            end if
+            if (ye .LT. ys) then
+               print(" ")
+               print("range end for y-coordinate = "+ye+"m is lower than start range = "+ys+"m")
+               print(" ")
+               exit
+            end if
+         end if
+      end if
+   end if 
+ 
+   if ( .not. isvar("ze")) then  
+      ze = zdim
+      if (parameter(73) .NE. "zdim") then
+         if (delta_z .EQ. -1) then
+            print(" ") 
+            print("You cannot choose an end value for z, there are preseted layers for z")
+            print(" ")           
+         else
+            if (stringtointeger(parameter(73)) .GT. zdim) then
+               print(" ")
+               print("range end for z-coordinate = "+parameter(73)+" is greater than last gridpoint = "+zdim)
                print(" ") 
-               exit
-            end if
-            if (ze .LT. 0 .OR. ze .LT. zs) then
-               print(" ")
-               print("range end for z-coordinate = "+ze+" is too small")
-               print(" ")
-               exit
-            end if
-         end if   
-         delete(data_all)
-      end if  
-   end do
-
+               exit 
+            end if 
+         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
+            if (stringtointeger(parameter(73)) .LE. 0)
+               print(" ")
+               print("range end for z-coordinate = "+parameter(73)+" is equal or lower than first gridpoint = 0")
+               print(" ")
+               exit
+            end if
+            if (stringtointeger(parameter(73)) .LE. zs) then
+               print(" ")
+               print("range end for z-coordinate = "+parameter(73)+" is equal or lower than start range = "+zs)
+               print(" ")
+               exit
+            end if
+            if (stringtodouble(parameter(73)) .EQ. zs+1) then
+               print(" ")
+               print("range end for z-coordinate = "+parameter(73)+" must be at least two more gridpoints greater than start range = "+zs)
+               print(" ")
+               exit
+            end if
+         else
+            if (stringtointeger(parameter(73)) .LT. 0)
+               print(" ")
+               print("range end for z-coordinate = "+parameter(73)+" is lower than first gridpoint = 0")
+               print(" ")
+               exit
+            end if
+            if (stringtointeger(parameter(73)) .LT. zs) then
+               print(" ")
+               print("range end for z-coordinate = "+parameter(73)+" is lower than start range = "+zs)
+               print(" ")
+               exit
+            end if
+         end if         
+         ze = stringtointeger(parameter(73)) 
+      end if 
+      end if
+   else
+      if (delta_z .EQ. -1) then
+         print(" ") 
+         print("You cannot choose an end value for z, there are preseted layers for z")
+         print(" ")
+         ze = zdim
+      else
+         if (ze .GT. zdim) then
+            print(" ")
+            print("range end for z-coordinate = "+ze+" is greater than last gridpoint = "+zdim)
+            print(" ") 
+            exit 
+         end if 
+         if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
+            if (ze .LE. 0)
+               print(" ")
+               print("range end for z-coordinate = "+ze+" is equal or lower than first gridpoint = 0")
+               print(" ")
+               exit
+            end if
+            if (ze .LE. zs) then
+               print(" ")
+               print("range end for z-coordinate = "+ze+" is equal or lower than start range = "+zs)
+               print(" ")
+               exit
+            end if
+            if (ze .EQ. zs+1) then
+               print(" ")
+               print("range end for z-coordinate = "+ze+" must be at least two more gridpoints greater than start range = "+zs)
+               print(" ")
+               exit
+            end if
+         else
+            if (ze .LT. 0)
+               print(" ")
+               print("range end for z-coordinate = "+ze+" is lower than first gridpoint = 0")
+               print(" ")
+               exit
+            end if
+            if (ze .LT. zs) then
+               print(" ")
+               print("range end for z-coordinate = "+ze+" is lower than start range = "+zs)
+               print(" ")
+               exit
+            end if
+         end if
+      end if
+   end if 
+
+   if (delta_x .NE. -1) then
+      do i=0,xdim    
+         if (xs .GT. x_d(i)-delta_x/2 .AND. xs .LE. x_d(i)+delta_x/2)then
+            xstart=i
+            break
+         end if
+      end do 
+      do i=0,xdim    
+         if (xe .GT. x_d(i)-delta_x/2 .AND. xe .LE. x_d(i)+delta_x/2)then
+            xend=i
+            break
+         end if
+      end do
+   end if
+   if (delta_y .NE. -1) then
+      do i=0,xdim    
+         if (ys .GT. y_d(i)-delta_y/2 .AND. ys .LE. y_d(i)+delta_y/2)then
+            ystart=i
+            break
+         end if
+      end do 
+      do i=0,xdim    
+         if (ye .GT. y_d(i)-delta_y/2 .AND. ye .LE. y_d(i)+delta_y/2)then
+            yend=i
+            break
+         end if
+      end do 
+   end if
+   
+   delete(xs)
+   delete(xe)   
+   delete(ys)
+   delete(ye)
+
+   xs=xstart
+   xe=xend
+   ys=ystart
+   ye=yend
+
+   if (xyc .EQ. 1 .OR. xzc .EQ.1)then
+      if (xe .EQ. xs+1) then
+         print(" ")
+         print("range end for x-coordinate="+parameter(65)+"m(=>value on file="+xe*delta_x+"m) must be at least two") 
+         print("more gridpoints(="+2*delta_x+"m) greater than start range="+parameter(63)+"m(=>value on file="+xs*delta_x+"m)")
+         print(" ")
+         exit
+      end if 
+   end if
+   if (xyc .EQ. 1 .OR. yzc .EQ.1)then
+      if (ye .EQ. ys+1) then
+         print(" ")
+         print("range end for y-coordinate="+parameter(69)+"m(=>value on file="+ye*delta_y+"m) must be at least two") 
+         print("more gridpoints(="+2*delta_y+"m) greater than start range="+parameter(67)+"m(=>value on file="+ys*delta_y+"m)")
+         print(" ")
+         exit
+      end if 
+   end if
+   if (xzc .EQ. 1 .OR. yzc .EQ.1)then
+      if (ze .EQ. zs+1) then
+         print(" ")
+         print("range end for x-coordinate="+parameter(73)+"(=> value on file="+ze+") must be at least two") 
+         print("more gridpoints greater than start range="+parameter(71)+"(=>value on file="+zs+")")
+         print(" ")
+         exit
+      end if 
+   end if
+  
    if (xyc .EQ. 1) then
       no_layer = (ze-zs)+1
+      data = new((/dim,nt,zdim+1,(ye-ys)+1,(xe-xs)+1/),float)
    end if
    if (xzc .EQ. 1) then
       no_layer = (ye-ys)+1
+      data = new((/dim,nt,(ze-zs)+1,ydim+1,(xe-xs)+1/),float)
    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)
+      data = new((/dim,nt,(ze-zs)+1,(ye-ys)+1,xdim+1/),float)
+   end if
+
+   MinVal = new(dim,float)
+   MaxVal = new(dim,float) 
 
    ; ****************************************************
@@ -571,5 +1219,5 @@
       
       data_all = f->$vNam(varn)$ 
-
+    
       if (vector .EQ. 1) then    
          check_vec1 = isStrSubset( vec1,","+vNam(varn)+",")
@@ -590,5 +1238,5 @@
       end if  
 
-      if(check) then 
+      if(check) then
 
          no_var=no_var+1
@@ -603,5 +1251,14 @@
          end if 
 
-         data(varn,:,:,:,:)=data_all(0:nt-1,zs:ze,ys:ye,xs:xe) 
+         if (xyc .EQ. 1) then
+            data(varn,:,:,:,:)=data_all(:,:,ys:ye,xs:xe) 
+         end if
+         if ( xzc .eq. 1 ) then
+            data(varn,:,:,:,:)=data_all(:,zs:ze,:,xs:xe) 
+         end if
+         if ( yzc .eq. 1) then
+            data(varn,:,:,:,:)=data_all(:,zs:ze,ys:ye,:) 
+         end if
+
          if (check_vec1) then
             vect1=data_all
@@ -615,5 +1272,8 @@
    	 data!2 = "z"
    	 data!3 = "y"
-   	 data!4 = "x"   
+   	 data!4 = "x"  
+
+         MinVal(varn) = min(data(varn,:,:,:,:))
+         MaxVal(varn) = max(data(varn,:,:,:,:))
 
       end if
@@ -622,6 +1282,13 @@
 
    end do
-  
-   plot_panel=new((/no_time*no_layer*no_var/),graphic)
+
+   if (no_var .EQ. 0) then
+      print(" ")
+      print("Please select a variable 'var=' or use the default value")
+      print(" ")
+      print("Your selection '"+var+"' does not exist on the input file")
+      print(" ")
+      exit
+   end if
 
    if (vector .EQ. 1) then
@@ -648,10 +1315,10 @@
    gsn_define_colormap(wks_ps,"rainbow+white")
 
-   plot=new((/no_time*no_layer*no_var/),graphic)
-
+   plot=new((/no_time*no_layer*no_var*2/),graphic)
+ 
    page = 0
    n=0
-
-   print("plots sorted by '"+sort+"'")
+   print(" ")
+   print("Plots sorted by '"+sort+"'")
    print(" ")
 
@@ -660,6 +1327,25 @@
    ; ***************************************************
 
+   if (vector .EQ. 1 .AND. parameter(45) .EQ. "plotvec") then
+      do lo = los, loe           				
+         do li = lis, lie
+            level = "=" + data&z(lo) + "m" 	  
+            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 = "t=" + data&t(li) +"s  z"+level 
+  	    vecres@tiXAxisString    = " "                                                
+            plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+            n=n+1
+         end do
+      end do
+   end if
+
    check = True
-
+ 
    do varn=dim-1,0,1 
 
@@ -682,12 +1368,20 @@
     
       if(check) then
-    
+         print(vNam(varn))
+         space=(MaxVal(varn)-MinVal(varn))/24
+  
+         cs_res@cnMinLevelValF = MinVal(varn)
+         cs_res@cnMaxLevelValF = MaxVal(varn) 
+
+         cs_res@cnLevelSpacingF = space
+      
          ; ****************************************************
          ; loops over time and layer
          ; **************************************************** 
 
+         
          do lo = los, loe           				
-            do li = lis, lie	                        
-                  
+            do li = lis, lie
+
                ; ****************************************************
                ; xy cross section
@@ -699,5 +1393,5 @@
                   cs_res@tiYAxisString = "y [m]"
                   cs_res@gsnLeftString = "Plot of "+vNam(varn)
-      
+                  
                   if ( sort .eq. "time" ) then
                      if ( data&z(lo) .eq. -1 ) then
@@ -706,35 +1400,21 @@
                         level = "=" + data&z(lo) + "m"
                      end if
-                     cs_res@gsnCenterString = "time=" + data&t(li) +"s  z"+level
+                     cs_res@gsnCenterString = "t=" + 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@gsnRightString   = " "             ; turn off right string
-  			   vecres@gsnLeftString    = " "             ; turn off left string
-  			   vecres@tiXAxisString    = " "   
-                           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
+                     if (vector .EQ. 1 .AND. check_vecp) then
+                        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   = " "
+                        vecres@gsnLeftString    = " " 
+  			vecres@tiXAxisString    = " "   
+                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                        overlay(plot(n), plot_vec)
                      end if                         
                   end if
-
+                 
                   if ( sort .eq. "layer" ) then
                      if ( data&z(li) .eq. -1 ) then
@@ -745,30 +1425,16 @@
                      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
+                     if (vector .EQ. 1 .AND. check_vecp) then
+                        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)
                      end if
                   end if
@@ -784,5 +1450,5 @@
                   cs_res@tiYAxisString   = "z [m]"
 		  cs_res@gsnLeftString = "Plot of "+vNam(varn)
-
+               
                   if ( sort .eq. "time" ) then
                      if ( data&z(lo) .eq. -1 ) then
@@ -793,30 +1459,16 @@
                      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
-  			   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,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
+                     if (vector .EQ. 1 .AND. check_vecp) then
+                        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,lo,:,:),vect2(li,lo,:,:),vecres)
+                        overlay(plot(n), plot_vec)
                      end if
                   end if
@@ -830,30 +1482,16 @@
                      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
+                     if (vector .EQ. 1 .AND. check_vecp) then
+                        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)
                      end if
                   end if                  
@@ -878,30 +1516,16 @@
                      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
-  			   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,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
+    	             if (vector .EQ. 1 .AND. check_vecp) then
+                        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,lo,:,:),vect2(li,lo,:,:),vecres)
+                        overlay(plot(n), plot_vec)
                      end if
                   end if
@@ -915,30 +1539,16 @@
                      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
+                     if (vector .EQ. 1 .AND. check_vecp)then
+                        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)
                      end if 
                   end if
@@ -947,5 +1557,5 @@
             end do 	
          end do  
-      end if 
+      end if
    end do
 
@@ -953,18 +1563,35 @@
    ; 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
+
+   if (vector .EQ. 1 .AND. parameter(45) .EQ. "plotvec")then 
+      if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
+         gsn_panel(wks_ps,plot(0:(no_time*no_layer*(no_var+1))-1),(/no_var+1,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)-1,no_lines*no_columns   
+            if ( np + no_lines*no_columns .gt. (no_layer*no_time*(no_var+1))-1) then 
+               gsn_panel(wks_ps, plot(np:(no_layer*no_time*(no_var+1))-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
+   else          
+      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
    end if
 
