Index: /palm/trunk/SCRIPTS/NCL/.ncl_preferences
===================================================================
--- /palm/trunk/SCRIPTS/NCL/.ncl_preferences	(revision 161)
+++ /palm/trunk/SCRIPTS/NCL/.ncl_preferences	(revision 162)
@@ -11,8 +11,8 @@
 ___'file_out'________[path/name; i.e. ~/.../test]_______________default = test____________[ALL]
 test
-___start_time_step___[0;1;...]__________________________________default = 1_______________[ALL]
-1
-___end_time_step_____[1;2;...]__________________________________default = nt______________[ALL]
-nt
+___start_time_step___[0;1;...]__________________________________default = t(0)____________[ALL]
+t(0)
+___end_time_step_____[1;2;...]__________________________________default = t(end)__________[ALL]
+t(end)
 ___no_columns________[1;2;...]__________________________________default = 1_______________[ALL]
 1
@@ -35,8 +35,20 @@
 ___max_z_____________[100;1000;...]_____________________________default = max(z_u)________[PR] 
 max(z_u)
+___over______________[0] or [1]_________________________________default = 0_______________[PR+TS]
+0
+___vector____________[0] or [1]_________________________________default = 0_______________[CS]
+0
+___'vec1'____________[,u,;,v,;...]______________________________default = vec1____________[CS]
+vec1
+___'vec2'____________[,u,;,v,;...]______________________________default = vec2____________[CS]
+vec2
+__'plotvec'__________[,u,;,v,;...]______________________________default = plotvec_________[CS]
+plotvec
+___ref_mag___________[0.1;0.02;...]_____________________________default = 0.05____________[CS]
+0.05
+___'mode'____________[Fill], [Line] or [Both]___________________default = Fill____________[CS]
+Fill
 ___'sort'____________[time] or [layer]__________________________default = time____________[CS]
 time
-___'mode'____________[Fill], [Line] or [Both]___________________default = Fill____________[CS]
-Fill
 ___'fill_mode'_______[AreaFill], [RasterFill] or [CellFill]_____default = AreaFill________[CS]
 AreaFill
@@ -49,27 +61,15 @@
 ___yzc_______________[0] or [1]_________________________________default = 0_______________[CS]
 0
-___xs________________[0,1,...]__________________________________default = 0_______________[CS]
-0
-___xe________________[1;2;...]__________________________________default = xdim____________[CS]
+___xs________________[0,1,...]__________________________________default = x0______________[CS+PR]
+x0
+___xe________________[1;2;...]__________________________________default = xdim____________[CS+PR]
 xdim
-___ys________________[0;1;...]__________________________________default = 0_______________[CS]
-0
+___ys________________[0;1;...]__________________________________default = y0______________[CS]
+y0
 ___ye________________[1;2;...]__________________________________default = ydim____________[CS]
 ydim
-___zs________________[0;1;...]__________________________________default = 0_______________[CS]
-0
+___zs________________[0;1;...]__________________________________default = z0______________[CS]
+z0
 ___ze________________[1;2;...]__________________________________default = zdim____________[CS]
 zdim
-___vector____________[0] or [1]_________________________________default = 0_______________[CS]
-0
-___'vec1'____________[,u,;,v,;...]______________________________default = vec1____________[CS]
-vec1
-___'vec2'____________[,u,;,v,;...]______________________________default = vec2____________[CS]
-vec2
-__'plotvec'__________[,u,;,v,;...]______________________________default = plotvec_________[CS]
-plotvec
-___ref_mag___________[0.1;0.02;...]_____________________________default = 0.05____________[CS]
-0.05
-___over______________[0] or [1]_________________________________default = 0_______________[TS+PR]
-0
 ___end of file ".ncl_preferences"_____________________________________________________________
Index: /palm/trunk/SCRIPTS/NCL/cross_sections.ncl
===================================================================
--- /palm/trunk/SCRIPTS/NCL/cross_sections.ncl	(revision 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
 
Index: /palm/trunk/SCRIPTS/NCL/profiles.ncl
===================================================================
--- /palm/trunk/SCRIPTS/NCL/profiles.ncl	(revision 161)
+++ /palm/trunk/SCRIPTS/NCL/profiles.ncl	(revision 162)
@@ -107,7 +107,7 @@
    if ( .not. isvar("over") ) then     		        ; switches overlaying plots on
       over = 0
-      if (stringtointeger(parameter(73)) .NE. 0) then
-         over = stringtointeger(parameter(73))
-         if (stringtointeger(parameter(73)) .NE. 1) then
+      if (stringtointeger(parameter(37)) .NE. 0) then
+         over = stringtointeger(parameter(37))
+         if (stringtointeger(parameter(37)) .NE. 1) then
             print(" ")
             print("Please set 'over' to 0 or 1")
@@ -456,9 +456,14 @@
       end do
    end if
- 
-   dimz  = dimsizes(z_w)
-   dimz  = dimsizes(z_u)
+   if (isvar("z_u") ) then 
+      dimz  = dimsizes(z_u)
+   else
+      if (isvar("z_w"))then 
+         dimz  = dimsizes(z_w)
+      end if
+   end if
    t_all = f->time
    nt    = dimsizes(t_all)
+   delta_t=t_all(nt-1)/nt
  
    ; ****************************************************	
@@ -467,36 +472,53 @@
 
    if ( .not. isvar("start_time_step") ) then   	
-      start_time_step = 1
-      if (parameter(13) .NE. "1") then
-         if (parameter(13) .LE. "1")
+      start_time_step=t_all(1)/3600
+      if (parameter(13) .NE. "t(0)") then
+         if (stringtodouble(parameter(13)) .GT. t_all(nt-1)/3600)
             print(" ")
-            print("Begin at least with time step 2")
+            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)) .GE. nt)
+         if (stringtofloat(parameter(13)) .LE. t_all(0)/3600)
             print(" ")
-            print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt-1))
+            print("'start_time_step' = "+ parameter(13) +"h is equal or lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
+            print(" ")
+            print("Please start profiles at least with second time step = "+t_all(1)+"s = "+t_all(1)/3600+"h") 
             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)
+      if (start_time_step .GT. t_all(nt-1)/3600)
          print(" ")
-         print("Begin at least with time step 1")
+         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 .GE. nt)
+      if (start_time_step .LE. t_all(0)/3600)
          print(" ")
-         print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt-1))
+         print("'start_time_step' = "+ start_time_step +"h is equal or lower than first time step = " + t_all(0)+"s = "+t_all(0)/3600+"h")
+         print(" ")
+         print("Please start profiles at least with second time step = "+t_all(1)+"s = "+t_all(1)/3600+"h")
          print(" ")
          exit
       end if
-      start_time_step = start_time_step - 1 
    end if
-
+   start_time_step = start_time_step*3600
+   if (start_time_step .GT. t_all(0) .AND. start_time_step .LT. t_all(1)+delta_t/2)then
+      st=1
+   end if
+   do i=2,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
@@ -504,47 +526,57 @@
 
    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 ")
-            print(" ") 
-            exit
-         end if
-         if (stringtointeger(parameter(15)) .GE. nt)
+      end_time_step = t_all(nt-1)/3600
+      if (parameter(15) .NE. "t(end)") then
+         if (stringtodouble(parameter(15)) .GT. t_all(nt-1)/3600)
             print(" ")
-            print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt-1))
+            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)) )
+         if (stringtodouble(parameter(15)) .LT. start_time_step/3600)
             print(" ")
-            print("'end_time_step' = "+ parameter(15) +" is lower than 'start_time_step' = "+parameter(13))
+            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))  
+         end_time_step = stringtodouble(parameter(15))
       end if    
    else
-      if (end_time_step .LE. 0)
+      if (end_time_step .GT. t_all(nt-1)/3600)
          print(" ")
-         print("'end_time_step' = "+end_time_step+ " is too small; 'end_time_step' should be at least 1 ")
+         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 .GE. nt)
+      if (end_time_step .LT. start_time_step/3600)
          print(" ")
-         print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt-1))
+         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
-      if (end_time_step .LT. start_time_step)
-            print(" ")
-            print("'end_time_step' = "+end_time_step +" is lower than 'start_time_step' = "+start_time_step)
-            print(" ")
-            exit
+   end if
+   end_time_step = end_time_step*3600
+   if (end_time_step .GT. t_all(0) .AND. end_time_step .LT. t_all(1)+delta_t/2)then
+      et=1
+   end if
+   do i=2,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_time_step = end_time_step - 1
-   end if
+   end do 
+
+   delete(start_time_step)
+   start_time_step=round(st,3)
+   delete(end_time_step)
+   end_time_step=round(et,3)
 
    ; ****************************************************
@@ -557,5 +589,5 @@
    end do
    
-   np = end_time_step-start_time_step  
+   np = end_time_step-start_time_step+1  
    if ( black .eq. 0 ) then
       res=True      
@@ -887,4 +919,5 @@
  
    n_o=0
+   count_var=0
    do varn = 0,dim-1 
         
@@ -908,5 +941,4 @@
          com=isStrSubset(c_var,","+vNam(varn)+"," )      
          if (com)                       
-            check = False
             data(varn,:,:) = temp(start_time_step:end_time_step,0:dimz-1)
             unit(varn) = temp@units
@@ -919,5 +951,4 @@
             mini(n_o)=min(data(varn,:,:))
             maxi(n_o)=max(data(varn,:,:)) 
-            varn=varn+1
             n_o=n_o+1
          end if
@@ -925,5 +956,6 @@
 
       if(check) then
-               
+         count_var=count_var+1
+     
          z = f->$vNam(varn+1)$
   
@@ -935,7 +967,23 @@
             res@gsnRightString     = unit(varn)
             res@trYMinF            = min_z
-            res@trYMaxF            = max_z  
-            res@trXMinF            = min(data(varn,:,:))
-            res@trXMaxF            = max(data(varn,:,:))        
+            res@trYMaxF            = max_z 
+            if (.not. isvar("xs")) then 
+               if (parameter(63) .NE. "x0") then
+                  res@trXMinF = stringtofloat(parameter(63))
+               else
+                  res@trXMinF            = min(data(varn,:,:))
+               end if
+            else 
+               res@trXMinF            = xs
+            end if
+            if (.not. isvar("xe")) then 
+               if (parameter(65) .NE. "xdim") then
+                  res@trXMaxF = stringtofloat(parameter(65))
+               else
+                  res@trXMaxF            = max(data(varn,:,:)) 
+               end if
+            else 
+               res@trXMaxF            = xe  
+            end if                     
             plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)  
          end if       
@@ -950,6 +998,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniu
-                  res@trXMaxF            = maxiu                  
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniu
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiu 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if                
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)  
                end if
@@ -964,6 +1028,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniv
-                  res@trXMaxF            = maxiv                   
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniv
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiv 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if                
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
                end if  
@@ -978,6 +1058,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniw
-                  res@trXMaxF            = maxiw                   
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniw
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiw 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if           
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if   
@@ -993,6 +1089,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minipt
-                  res@trXMaxF            = maxipt
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minipt
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxipt 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if  
@@ -1007,6 +1119,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minivpt
-                  res@trXMaxF            = maxivpt 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minivpt
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxivpt 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
                end if 
@@ -1021,6 +1149,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minilpt
-                  res@trXMaxF            = maxilpt 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minilpt
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxilpt 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if 
@@ -1036,6 +1180,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniq
-                  res@trXMaxF            = maxiq 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minilq
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxilq 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if
@@ -1050,6 +1210,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniqv
-                  res@trXMaxF            = maxiqv
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minilqv
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxilqv 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
                end if
@@ -1064,54 +1240,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniql
-                  res@trXMaxF            = maxiql 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "rho") then                
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:))
-                  plot_q = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:)) 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "s") then                
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:))
-                  plot_s = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:)) 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "sa") then                 
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:)) 
-                  plot_sa = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:))  
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniql
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiql 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if
@@ -1127,6 +1271,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minie
-                  res@trXMaxF            = maxie 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minie
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxie 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if 
@@ -1141,6 +1301,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minies
-                  res@trXMaxF            = maxies
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minies
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxies 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
                end if
@@ -1156,6 +1332,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minikm
-                  res@trXMaxF            = maxikm
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minikm
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxikm 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if 
@@ -1170,23 +1362,23 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minikh
-                  res@trXMaxF            = maxikh 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minikh
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxikh 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "l") then                 
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:))
-                  plot_l = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:)) 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if
             end if
@@ -1201,6 +1393,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpup
-                  res@trXMaxF            = maxiwpup
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpup
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpup 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1215,6 +1423,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsus
-                  res@trXMaxF            = maxiwsus 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsus
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsus 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1229,6 +1453,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwu
-                  res@trXMaxF            = maxiwu
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwu
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwu 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1244,6 +1484,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpvp
-                  res@trXMaxF            = maxiwpvp
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpvp
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpvp 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1258,6 +1514,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsvs
-                  res@trXMaxF            = maxiwsvs 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsvs
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsvs 
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1272,6 +1544,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwv
-                  res@trXMaxF            = maxiwv 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwv
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwv
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1287,6 +1575,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpptp
-                  res@trXMaxF            = maxiwpptp 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpptp
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpptp
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1301,6 +1605,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwspts
-                  res@trXMaxF            = maxiwspts 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwspts
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwspts
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1315,6 +1635,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpt
-                  res@trXMaxF            = maxiwpt 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpt
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpt
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1330,6 +1666,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsptsBC
-                  res@trXMaxF            = maxiwsptsBC 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsptsBC
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsptsBC
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if 
@@ -1344,6 +1696,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwptBC
-                  res@trXMaxF            = maxiwptBC 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwptBC
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwptBC
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1359,6 +1727,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpvptp
-                  res@trXMaxF            = maxiwpvptp 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpvptp
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpvptp
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1373,6 +1757,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsvpts
-                  res@trXMaxF            = maxiwsvpts 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsvpts
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsvpts
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1387,6 +1787,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwvpt
-                  res@trXMaxF            = maxiwvpt 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwvpt
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwvpt
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1402,6 +1818,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpqp
-                  res@trXMaxF            = maxiwpqp 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpqp
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpqp
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1416,6 +1848,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsqs
-                  res@trXMaxF            = maxiwsqs 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsqs
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsqs
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1430,6 +1878,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwq
-                  res@trXMaxF            = maxiwq 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwq
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwq
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1445,6 +1909,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpqvp
-                  res@trXMaxF            = maxiwpqvp 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpqvp
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpqvp
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if 
@@ -1459,6 +1939,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsqvs
-                  res@trXMaxF            = maxiwsqvs 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsqvs
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsqvs
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1473,6 +1969,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwqv
-                  res@trXMaxF            = maxiwqv 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwqv
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwqv
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1488,6 +2000,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpsp
-                  res@trXMaxF            = maxiwpsp 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpsp
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpsp
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if 
@@ -1502,6 +2030,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsss
-                  res@trXMaxF            = maxiwsss 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsss
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsss
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1516,6 +2060,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniws
-                  res@trXMaxF            = maxiws 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniws
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiws
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if 
@@ -1531,6 +2091,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpsap
-                  res@trXMaxF            = maxiwpsap 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpsap
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpsap
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if 
@@ -1545,6 +2121,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwssas
-                  res@trXMaxF            = maxiwssas 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwssas
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwssas
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
                end if
@@ -1559,22 +2151,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsa
-                  res@trXMaxF            = maxiwsa 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "wses") then                
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:))
-                  plot_wses = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min(data(varn,:,:))
-                  res@trXMaxF            = max(data(varn,:,:)) 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsa
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsa
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1590,6 +2182,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minius2
-                  res@trXMaxF            = maxius2 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minius2
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxius2
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if 
@@ -1604,6 +2212,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minivs2
-                  res@trXMaxF            = maxivs2 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = minivs2
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxivs2
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
                end if
@@ -1618,96 +2242,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniws2
-                  res@trXMaxF            = maxiws2 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "pts2") then
-               minipts2=min(data(varn,:,:))
-               maxipts2=max(data(varn,:,:))                 
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minipts2
-                  res@trXMaxF            = maxipts2
-                  plot_pts2 = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = minipts2
-                  res@trXMaxF            = maxipts2 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "ws3") then
-               miniws3=min(data(varn,:,:))
-               maxiws3=max(data(varn,:,:))                 
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniws3
-                  res@trXMaxF            = maxiws3
-                  plot_ws3 = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniws3
-                  res@trXMaxF            = maxiws3 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "Sw") then
-               miniSw=min(data(varn,:,:))
-               maxiSw=max(data(varn,:,:))                 
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniSw
-                  res@trXMaxF            = maxiSw
-                  plot_Sw = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniSw
-                  res@trXMaxF            = maxiSw 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "ws2pts") then
-               miniws2pts=min(data(varn,:,:))
-               maxiws2pts=max(data(varn,:,:))                 
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniws2pts
-                  res@trXMaxF            = maxiws2pts
-                  plot_ws2pts = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniws2pts
-                  res@trXMaxF            = maxiws2pts 
-                  plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
-               end if
-            end if
-
-            if (vNam(varn) .EQ. "wspts2") then
-               miniwspts2=min(data(varn,:,:))
-               maxiwspts2=max(data(varn,:,:))                 
-               if (over .EQ. 1) then
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwspts2
-                  res@trXMaxF            = maxiwspts2
-                  plot_wspts2 = gsn_csm_xy(wks,data(varn,:,:),z_w,res)
-               else
-                  res@gsnLeftString      = vNam(varn)
-                  res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwspts2
-                  res@trXMaxF            = maxiwspts2 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniws2
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiws2
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1723,6 +2273,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwsususodz
-                  res@trXMaxF            = maxiwsususodz 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwsususodz
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwsususodz
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
                end if  
@@ -1737,6 +2303,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwspsodz
-                  res@trXMaxF            = maxiwspsodz 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwspsodz
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwspsodz
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_u,res)
                end if 
@@ -1751,6 +2333,22 @@
                   res@gsnLeftString      = vNam(varn)
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = miniwpeodz
-                  res@trXMaxF            = maxiwpeodz 
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF            = miniwpeodz
+                     end if
+                  else 
+                     res@trXMinF            = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF            = maxiwpeodz
+                     end if
+                  else 
+                     res@trXMaxF            = xe  
+                  end if 
                   plot(n) = gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
                end if
@@ -1763,5 +2361,14 @@
       delete(temp) 
    end do
-   
+
+   if (count_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 (over .EQ. 1 ) then
 
@@ -1836,10 +2443,26 @@
             z = f->$vNam(varn+1)$
            
-            ;res@gsnLeftString      = vNam(varn)
+            res@gsnLeftString      = vNam(varn)
             res@gsnRightString     = unit(varn)
             res@trYMinF            = min_z
             res@trYMaxF            = max_z  
-            res@trXMinF            = min(data(varn,:,:))
-            res@trXMaxF            = max(data(varn,:,:))
+            if (.not. isvar("xs")) then 
+               if (parameter(63) .NE. "x0") then
+                  res@trXMinF = stringtofloat(parameter(63))
+               else
+                  res@trXMinF = min(data(varn,:,:))
+               end if
+            else 
+               res@trXMinF = xs
+            end if
+            if (.not. isvar("xe")) then 
+               if (parameter(65) .NE. "xdim") then
+                  res@trXMaxF = stringtofloat(parameter(65))
+               else
+                  res@trXMaxF = max(data(varn,:,:))
+            end if
+            else 
+               res@trXMaxF = xe  
+            end if 
             plot(n) = gsn_csm_xy(wks,data(varn,:,:),z,res)  
             
@@ -1848,6 +2471,22 @@
                   res@gsnLeftString      = "u, v and w"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniu,miniv,miniw/))
-                  res@trXMaxF            = max((/maxiu,maxiv,maxiw/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniu,miniv,miniw/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiu,maxiv,maxiw/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if 
                   if (vNam(varn) .EQ. "w") then 
                      plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -1885,6 +2524,22 @@
                   res@gsnLeftString      = "pt, vpt and lpt"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/minipt,minivpt,minilpt/))
-                  res@trXMaxF            = max((/maxipt,maxivpt,maxilpt/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/minipt,minivpt,minilpt/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxipt,maxivpt,maxilpt/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
@@ -1918,6 +2573,22 @@
                   res@gsnLeftString      = "q, qv and ql"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniq,miniqv,miniql/))
-                  res@trXMaxF            = max((/maxiq,maxiqv,maxiql/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniq,miniqv,miniql/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiq,maxiqv,maxiql/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
@@ -1946,20 +2617,28 @@
                   continue   
                end if
-            end if           
-            if (vNam(varn) .EQ. "rho") then
-               plot(n) = plot_rho
-            end if
-            if (vNam(varn) .EQ. "s") then
-               plot(n) = plot_s
-            end if
-            if (vNam(varn) .EQ. "sa") then
-               plot(n) = plot_sa
-            end if
+            end if    
+           
             if (vNam(varn) .EQ. "e" .OR. vNam(varn) .EQ. "es") then
                if (e .EQ. 0) then
                   res@gsnLeftString      = "e and es"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/minie,minies/))
-                  res@trXMaxF            = max((/maxie,maxies/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/minie,minies/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxie,maxies/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
@@ -1993,6 +2672,22 @@
                   res@gsnLeftString      = "km and kh"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/minikm,minikh/))
-                  res@trXMaxF            = max((/maxikm,maxikh/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/minikm,minikh/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxikm,maxikh/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_u,res) 
@@ -2022,13 +2717,27 @@
                end if
             end if           
-            if (vNam(varn) .EQ. "l") then
-               plot(n) = plot_l
-            end if
+           
             if (vNam(varn) .EQ. "wpup" .OR. vNam(varn) .EQ. "wsus" .OR. vNam(varn) .EQ. "wu") then
                if (wpup .EQ. 0) then
                   res@gsnLeftString      = "wpup, wsus and wu"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwpup,miniwsus,miniwu/))
-                  res@trXMaxF            = max((/maxiwpup,maxiwsus,maxiwu/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwpup,miniwsus,miniwu/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwpup,maxiwsus,maxiwu/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2062,6 +2771,22 @@
                   res@gsnLeftString      = "wpvp, wsus and wv"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwpvp,miniwsvs,miniwv/))
-                  res@trXMaxF            = max((/maxiwpvp,maxiwsvs,maxiwv/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwpvp,miniwsvs,miniwv/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwpvp,maxiwsvs,maxiwv/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2095,6 +2820,22 @@
                   res@gsnLeftString      = "wpptp, wspts and wv"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwpptp,miniwspts,miniwpt/))
-                  res@trXMaxF            = max((/maxiwpptp,maxiwspts,maxiwpt/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwpptp,miniwspts,miniwpt/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwpptp,maxiwspts,maxiwpt/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2128,6 +2869,22 @@
                   res@gsnLeftString      = "wsptsBC and wptBC"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwsptsBC,miniwptBC/))
-                  res@trXMaxF            = max((/maxiwsptsBC,maxiwptBC/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwsptsBC,miniwptBC/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwsptsBC,maxiwptBC/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2161,6 +2918,22 @@
                   res@gsnLeftString      = "wpvptp, wsvpts and wv"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwpvptp,miniwsvpts,miniwvpt/))
-                  res@trXMaxF            = max((/maxiwpvptp,maxiwsvpts,maxiwvpt/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwpvptp,miniwsvpts,miniwvpt/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwpvptp,maxiwsvpts,maxiwvpt/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2194,6 +2967,22 @@
                   res@gsnLeftString      = "wpqp, wsqs and wq"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwpqp,miniwsqs,miniwq/))
-                  res@trXMaxF            = max((/maxiwpqp,maxiwsqs,maxiwq/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwpqp,miniwsqs,miniwq/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwpqp,maxiwsqs,maxiwq/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2227,6 +3016,22 @@
                   res@gsnLeftString      = "wpqvp, wsqvs and wqv"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwpqp,miniwsqvs,miniwqv/))
-                  res@trXMaxF            = max((/maxiwpqp,maxiwsqvs,maxiwqv/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwpqp,miniwsqvs,miniwqv/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwpqp,maxiwsqvs,maxiwqv/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2260,6 +3065,22 @@
                   res@gsnLeftString      = "wpsp, wsss and ws"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwpsp,miniwsss,miniws/))
-                  res@trXMaxF            = max((/maxiwpsp,maxiwsss,maxiws/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwpsp,miniwsss,miniws/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwpsp,maxiwsss,maxiws/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2293,6 +3114,22 @@
                   res@gsnLeftString      = "wpsap, wssas and wsa"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwpsap,miniwssas,miniwsa/))
-                  res@trXMaxF            = max((/maxiwpsap,maxiwssas,maxiwsa/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwpsap,miniwssas,miniwsa/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwpsap,maxiwssas,maxiwsa/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2322,13 +3159,27 @@
                end if
             end if
-            if (vNam(varn) .EQ. "wses") then
-               plot(n) = plot_wses
-            end if
+          
             if (vNam(varn) .EQ. "us2" .OR. vNam(varn) .EQ. "vs2" .OR. vNam(varn) .EQ. "ws2") then
                if (us2 .EQ. 0) then
                   res@gsnLeftString      = "us2, vs2 and ws2"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/minius2,minivs2,miniws2/))
-                  res@trXMaxF            = max((/maxius2,maxivs2,maxiws2/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/minius2,minivs2,miniws2/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxius2,maxivs2,maxiws2/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2358,25 +3209,27 @@
                end if
             end if
-            if (vNam(varn) .EQ. "pts2") then
-               plot(n) = plot_pts2 
-            end if
-            if (vNam(varn) .EQ. "ws3") then
-               plot(n) = plot_ws3
-            end if
-            if (vNam(varn) .EQ. "Sw") then
-               plot(n) = plot_Sw
-            end if
-            if (vNam(varn) .EQ. "ws2pts") then
-               plot(n) = plot_ws2pts
-            end if
-            if (vNam(varn) .EQ. "wspts2") then
-               plot(n) = plot_wspts2 
-            end if
+           
             if (vNam(varn) .EQ. "wsususodz" .OR. vNam(varn) .EQ. "wspsodz" .OR. vNam(varn) .EQ. "wpeodz" ) then
                if (wsususodz .EQ. 0) then 
                   res@gsnLeftString      = "wsususodz, wspsodz and ws2"
                   res@gsnRightString     = unit(varn)
-                  res@trXMinF            = min((/miniwsususodz,miniwspsodz,miniwpeodz/))
-                  res@trXMaxF            = max((/maxiwsususodz,maxiwspsodz,maxiwpeodz/))
+                  if (.not. isvar("xs")) then 
+                     if (parameter(63) .NE. "x0") then
+                        res@trXMinF = stringtofloat(parameter(63))
+                     else
+                        res@trXMinF = min((/miniwsususodz,miniwspsodz,miniwpeodz/))
+                     end if
+                  else 
+                     res@trXMinF = xs
+                  end if
+                  if (.not. isvar("xe")) then 
+                     if (parameter(65) .NE. "xdim") then
+                        res@trXMaxF = stringtofloat(parameter(65))
+                     else
+                        res@trXMaxF = max((/maxiwsususodz,maxiwspsodz,maxiwpeodz/))
+                     end if
+                  else 
+                     res@trXMaxF = xe  
+                  end if
 
                   plot(n) =  gsn_csm_xy(wks,data(varn,:,:),z_w,res) 
@@ -2449,6 +3302,22 @@
                      res@gsnLeftString  = c_var
                      res@gsnRightString = unit(varn)
-                     res@trXMinF        = min(mini)
-                     res@trXMaxF        = max(maxi)
+                     if (.not. isvar("xs")) then 
+                        if (parameter(63) .NE. "x0") then
+                           res@trXMinF = stringtofloat(parameter(63))
+                        else
+                           res@trXMinF = min(mini)
+                        end if
+                     else 
+                        res@trXMinF = xs
+                     end if
+                     if (.not. isvar("xe")) then 
+                        if (parameter(65) .NE. "xdim") then
+                           res@trXMaxF = stringtofloat(parameter(65))
+                        else
+                           res@trXMaxF = max(maxi)
+                        end if
+                     else 
+                        res@trXMaxF = xe  
+                     end if 
                   end if
                end if
@@ -2496,5 +3365,5 @@
    ; merge plots onto one page 
    ; ***************************************************
- print(n)
+
    do m=0,n-1
       plot_(m)=plot(n-1-m)
Index: /palm/trunk/SCRIPTS/NCL/timeseries.ncl
===================================================================
--- /palm/trunk/SCRIPTS/NCL/timeseries.ncl	(revision 161)
+++ /palm/trunk/SCRIPTS/NCL/timeseries.ncl	(revision 162)
@@ -63,7 +63,7 @@
    if ( .not. isvar("over") ) then     		        ; switches overlaying plots on
       over = 0
-      if (stringtointeger(parameter(73)) .NE. 0) then
-         over = stringtointeger(parameter(73))
-         if (stringtointeger(parameter(73)) .NE. 1) then
+      if (stringtointeger(parameter(37)) .NE. 0) then
+         over = stringtointeger(parameter(37))
+         if (stringtointeger(parameter(37)) .NE. 1) then
             print(" ")
             print("Please set 'over' to 0 or 1")
@@ -93,4 +93,5 @@
    t_all = f->time
    nt  = dimsizes(t_all)
+   delta_t=t_all(nt-1)/nt
 
    ; ****************************************************	
@@ -99,36 +100,53 @@
 
    if ( .not. isvar("start_time_step") ) then   	
-      start_time_step = 0
-      if (parameter(13) .NE. "1") then
-         if (parameter(13) .LE. "0")
-            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)) .GE. t_all(nt-1)/3600)
+            print(" ")
+            print("'start_time_step' = "+ parameter(13) +"h is equal or 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)) .GE. nt)
-            print(" ")
-            print("'start_time_step' = "+ parameter(13) +" is greater than available time steps = " + (nt-1))
+         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 .GE. t_all(nt-1)/3600)
+         print(" ")
+         print("'start_time_step' = "+ start_time_step +"h is equal or 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 .GE. nt)
-         print(" ")
-         print("'start_time_step' = "+ start_time_step +" is greater than available time steps = " + (nt-1))
+      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
-
+   start_time_step = start_time_step*3600
+   do i=0,nt-2     
+      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
+   if (start_time_step .GE. t_all(nt-1)-delta_t/2 .AND. start_time_step .LT. t_all(nt-1)) then
+      st=nt-2   
+   end if
+      
    ; ****************************************************
    ; end of time step and different types of mistakes that could be done
@@ -136,47 +154,70 @@
 
    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)) .LE. t_all(0)/3600)
+            print(" ")
+            print("'end_time_step' = "+parameter(15)+ "h is lower or equal 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)) .GE. nt)
-            print(" ")
-            print("'end_time_step' = "+ parameter(15) +" is greater than available time steps = " + (nt-1))  
+         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(13))
+         if (stringtodouble(parameter(15)) .LE. start_time_step/3600)
+            print(" ")
+            print("'end_time_step' = "+ parameter(15) +"h is equal or 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 .LE. t_all(0)/3600)
+         print(" ")
+         print("'end_time_step' = "+end_time_step+ "h is lower or equal 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 .GE. nt)
-         print(" ")
-         print("'end_time_step' = "+ end_time_step +" is greater than available time steps = "+(nt-1))
+      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 .LE. start_time_step/3600)
+         print(" ")
+         print("'end_time_step' = "+ end_time_step +"h is equal or 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)
 
    t = f->time(start_time_step:end_time_step)
@@ -266,5 +307,6 @@
       plot_splptz  = gsn_csm_xy(wks_ps,t,data_0(:),res)
    end if
-  
+
+   count_var=0
    do varn = dim-1,0,1 
 
@@ -284,5 +326,6 @@
      
       if(check) then
-       
+         count_var=count_var+1
+
          data_all = f ->$vNam(varn)$ 
          unit(varn) = data_all@units
@@ -432,4 +475,13 @@
    end do
 
+   if (count_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 
+
    do varn = dim-1,0,1 
       
@@ -457,12 +509,9 @@
             res@tiXAxisFontHeightF   = 0.07
             res@txFontHeightF        = 0.07
-            res@tiYAxisFontHeightF   = 0.07
-            res@xyLabelMode          = "Custom"
-            
+            res@tiYAxisFontHeightF   = 0.07            
 
             if (vNam(varn) .EQ. "E")
                E=0
                res@xyLineColors     = (/237/)
-               res@xyExplicitLabels =  "E"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 237
@@ -474,5 +523,4 @@
                Es=0
                res@xyLineColors            = (/144/)
-               res@xyExplicitLabels =  "Es"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 144
@@ -483,5 +531,4 @@
                us=0
                res@xyLineColors            = (/237/)
-               res@xyExplicitLabels =  "us"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 237
@@ -493,5 +540,4 @@
                ws=0
                res@xyLineColors            = (/144/)
-               res@xyExplicitLabels =  "ws"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 144
@@ -502,5 +548,4 @@
                u=0
                res@xyLineColors            = (/237/)
-               res@xyExplicitLabels =  "umax"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 237
@@ -512,5 +557,4 @@
                v=0
                res@xyLineColors            = (/144/)
-               res@xyExplicitLabels =  "vmax"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 144
@@ -520,5 +564,4 @@
                w=0
                res@xyLineColors            = (/80/)
-               res@xyExplicitLabels =  "wmax"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 80
@@ -529,5 +572,4 @@
                zw=0
                res@xyLineColors            = (/237/) 
-               res@xyExplicitLabels =  "z_i_wpt"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 237
@@ -539,5 +581,4 @@
                z=0
                res@xyLineColors            = (/144/)  
-               res@xyExplicitLabels =  "z_i_pt"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 144 
@@ -548,5 +589,4 @@
                w0=0
                res@xyLineColors            = (/237/)
-               res@xyExplicitLabels =  "wpptp0"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 237
@@ -558,5 +598,4 @@
                wp=0
                res@xyLineColors            = (/144/)
-               res@xyExplicitLabels =  "wpptp"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 144
@@ -566,5 +605,4 @@
                wt=0
                res@xyLineColors            = (/80/)
-               res@xyExplicitLabels =  "wpt"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 80
@@ -575,5 +613,4 @@
                p=0
                res@xyLineColors            = (/237/)
-               res@xyExplicitLabels =  "pt_0_"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 237
@@ -585,5 +622,4 @@
                pz=0
                res@xyLineColors            = (/144/)
-               res@xyExplicitLabels =  "pt_zp_"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 144
@@ -594,5 +630,4 @@
                x=0
                res@xyLineColors            = (/237/)
-               res@xyExplicitLabels =  "splptx"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 237
@@ -604,5 +639,4 @@
                y=0
                res@xyLineColors            = (/144/)
-               res@xyExplicitLabels =  "splpty"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 144
@@ -612,5 +646,4 @@
                z=0
                res@xyLineColors            = (/80/)
-               res@xyExplicitLabels =  "splptz"
                res@xyLineLabelFontHeightF = 0.05 
                res@xyLineLabelFontColor   = 80
@@ -645,5 +678,24 @@
                overlay(plot_E,plot_Es)
                n=n+1
-               plot_ps(n) = plot_E               
+               plot_ps(n) = plot_E    
+
+               ; ***************************************************
+               ; legend for combined plot
+               ; ***************************************************
+      
+               lgres                    = True
+               lgMonoDashIndex          = False
+               lgres@lgLabelFont        = "helvetica"    
+               lgres@lgLabelFontHeightF = .1            
+               lgres@vpWidthF           = 0.4           
+               lgres@vpHeightF          = 0.4          
+               lgres@lgDashIndexes      = (/0,0,0/)
+               lgres@lgLineColors       = (/237,144,80/)
+               lbid = gsn_create_legend(wks_ps,2,(/"E","Es"/),lgres)       
+
+               amres = True
+               amres@amParallelPosF   = 0.6		    
+               amres@amOrthogonalPosF = -0.2            
+               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)            
             end if
             if (vNam(varn) .EQ. "Es" .AND. E .NE. 1) then
@@ -651,5 +703,24 @@
                overlay(plot_E,plot_Es)
                n=n+1
-               plot_ps(n) = plot_E           
+               plot_ps(n) = plot_E  
+
+               ; ***************************************************
+               ; legend for combined plot
+               ; ***************************************************
+      
+               lgres                    = True
+               lgMonoDashIndex          = False
+               lgres@lgLabelFont        = "helvetica"    
+               lgres@lgLabelFontHeightF = .1            
+               lgres@vpWidthF           = 0.4           
+               lgres@vpHeightF          = 0.4          
+               lgres@lgDashIndexes      = (/0,0,0/)
+               lgres@lgLineColors       = (/237,144,80/)
+               lbid = gsn_create_legend(wks_ps,2,(/"E","Es"/),lgres)       
+
+               amres = True
+               amres@amParallelPosF   = 0.6		    
+               amres@amOrthogonalPosF = -0.2            
+               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)         
             end if
 
@@ -659,4 +730,23 @@
                n=n+1
                plot_ps(n) = plot_us
+
+               ; ***************************************************
+               ; legend for combined plot
+               ; ***************************************************
+      
+               lgres                    = True
+               lgMonoDashIndex          = False
+               lgres@lgLabelFont        = "helvetica"    
+               lgres@lgLabelFontHeightF = .1           
+               lgres@vpWidthF           = 0.4           
+               lgres@vpHeightF          = 0.4          
+               lgres@lgDashIndexes      = (/0,0,0/)
+               lgres@lgLineColors       = (/237,144,80/)
+               lbid = gsn_create_legend(wks_ps,2,(/"us","ws"/),lgres)       
+
+               amres = True
+               amres@amParallelPosF   = 0.6		    
+               amres@amOrthogonalPosF = -0.2            
+               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
             end if
             if (vNam(varn) .EQ. "ws" .AND. us .NE. 1) then
@@ -665,4 +755,23 @@
                n=n+1
                plot_ps(n) = plot_us
+
+               ; ***************************************************
+               ; legend for combined plot
+               ; ***************************************************
+      
+               lgres                    = True
+               lgMonoDashIndex          = False
+               lgres@lgLabelFont        = "helvetica"    
+               lgres@lgLabelFontHeightF = .1            
+               lgres@vpWidthF           = 0.4           
+               lgres@vpHeightF          = 0.4          
+               lgres@lgDashIndexes      = (/0,0,0/)
+               lgres@lgLineColors       = (/237,144,80/)
+               lbid = gsn_create_legend(wks_ps,2,(/"us","ws"/),lgres)       
+
+               amres = True
+               amres@amParallelPosF   = 0.6		    
+               amres@amOrthogonalPosF = -0.2            
+               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
             end if
          
@@ -674,4 +783,23 @@
                   n=n+1
                   plot_ps(n) = plot_umax
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"umax","vmax","wmax"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
                end if
             end if
@@ -683,4 +811,23 @@
                   n=n+1
                   plot_ps(n) = plot_umax
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1          
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"umax","vmax","wmax"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
                end if
             end if
@@ -692,4 +839,23 @@
                   n=n+1
                   plot_ps(n) = plot_umax
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"umax","vmax","wmax"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
                end if
             end if
@@ -699,5 +865,24 @@
                overlay(plot_z_i_wpt,plot_z_i_pt)
                n=n+1
-               plot_ps(n) = plot_z_i_wpt           
+               plot_ps(n) = plot_z_i_wpt 
+       
+               ; ***************************************************
+               ; legend for combined plot
+               ; ***************************************************
+      
+               lgres                    = True
+               lgMonoDashIndex          = False
+               lgres@lgLabelFont        = "helvetica"    
+               lgres@lgLabelFontHeightF = .1            
+               lgres@vpWidthF           = 0.4           
+               lgres@vpHeightF          = 0.4          
+               lgres@lgDashIndexes      = (/0,0,0/)
+               lgres@lgLineColors       = (/237,144,80/)
+               lbid = gsn_create_legend(wks_ps,2,(/"z_i_wpt","z_i_pt"/),lgres)       
+
+               amres = True
+               amres@amParallelPosF   = 0.6		    
+               amres@amOrthogonalPosF = -0.2            
+               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)          
             end if 
             if (vNam(varn) .EQ. "z_i_pt" .AND. zw .NE. 1) then
@@ -705,5 +890,24 @@
                overlay(plot_z_i_wpt,plot_z_i_pt)
                n=n+1
-               plot_ps(n) = plot_z_i_wpt           
+               plot_ps(n) = plot_z_i_wpt
+
+               ; ***************************************************
+               ; legend for combined plot
+               ; ***************************************************
+      
+               lgres                    = True
+               lgMonoDashIndex          = False
+               lgres@lgLabelFont        = "helvetica"    
+               lgres@lgLabelFontHeightF = .1            
+               lgres@vpWidthF           = 0.4           
+               lgres@vpHeightF          = 0.4          
+               lgres@lgDashIndexes      = (/0,0,0/)
+               lgres@lgLineColors       = (/237,144,80/)
+               lbid = gsn_create_legend(wks_ps,2,(/"z_i_wpt","z_i_pt"/),lgres)       
+
+               amres = True
+               amres@amParallelPosF   = 0.6		    
+               amres@amOrthogonalPosF = -0.2            
+               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)           
             end if    
          
@@ -715,4 +919,23 @@
                   n=n+1
                   plot_ps(n) = plot_wpptp0 
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
                end if            
             end if
@@ -723,5 +946,24 @@
                  overlay(plot_wpptp0,plot_wpt) 
                  n=n+1
-                 plot_ps(n) = plot_wpptp0  
+                 plot_ps(n) = plot_wpptp0
+                 
+                 ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)  
                end if           
             end if
@@ -732,5 +974,24 @@
                   overlay(plot_wpptp0,plot_wpt) 
                   n=n+1
-                  plot_ps(n) = plot_wpptp0             
+                  plot_ps(n) = plot_wpptp0
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"wpptp0","wpptp","wpt"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
                end if
             end if
@@ -740,5 +1001,24 @@
                overlay(plot_pt_0_,plot_pt_zp_)
                n=n+1 
-               plot_ps(n) = plot_pt_0_           
+               plot_ps(n) = plot_pt_0_  
+      
+               ; ***************************************************
+               ; legend for combined plot
+               ; ***************************************************
+      
+               lgres                    = True
+               lgMonoDashIndex          = False
+               lgres@lgLabelFont        = "helvetica"    
+               lgres@lgLabelFontHeightF = .1            
+               lgres@vpWidthF           = 0.4           
+               lgres@vpHeightF          = 0.4          
+               lgres@lgDashIndexes      = (/0,0,0/)
+               lgres@lgLineColors       = (/237,144,80/)
+               lbid = gsn_create_legend(wks_ps,2,(/"pt_0_","pt_zp_"/),lgres)       
+
+               amres = True
+               amres@amParallelPosF   = 0.6		    
+               amres@amOrthogonalPosF = -0.2            
+               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
             end if 
             if (vNam(varn) .EQ. "pt_zp_" .AND. p .NE. 1) then
@@ -746,5 +1026,24 @@
                overlay(plot_pt_0_,plot_pt_zp_)
                n=n+1 
-               plot_ps(n) = plot_pt_0_           
+               plot_ps(n) = plot_pt_0_   
+
+               ; ***************************************************
+               ; legend for combined plot
+               ; ***************************************************
+      
+               lgres                    = True
+               lgMonoDashIndex          = False
+               lgres@lgLabelFont        = "helvetica"    
+               lgres@lgLabelFontHeightF = .1            
+               lgres@vpWidthF           = 0.4           
+               lgres@vpHeightF          = 0.4          
+               lgres@lgDashIndexes      = (/0,0,0/)
+               lgres@lgLineColors       = (/237,144,80/)
+               lbid = gsn_create_legend(wks_ps,2,(/"pt_0_","pt_zp_"/),lgres)       
+
+               amres = True
+               amres@amParallelPosF   = 0.6		    
+               amres@amOrthogonalPosF = -0.2            
+               annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
             end if
             
@@ -755,5 +1054,24 @@
                   overlay(plot_splptx,plot_splptz)
                   n=n+1
-                  plot_ps(n) = plot_splptx           
+                  plot_ps(n) = plot_splptx    
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"splptx","splpty","splptz"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
                end if
             end if 
@@ -765,4 +1083,23 @@
                   n=n+1
                   plot_ps(n) = plot_splptx
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"splptx","splpty","splptz"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
                end if           
             end if 
@@ -773,5 +1110,24 @@
                   overlay(plot_splptx,plot_splptz)
                   n=n+1
-                  plot_ps(n) = plot_splptx    
+                  plot_ps(n) = plot_splptx  
+
+                  ; ***************************************************
+                  ; legend for combined plot
+                  ; ***************************************************
+      
+                  lgres                    = True
+                  lgMonoDashIndex          = False
+                  lgres@lgLabelFont        = "helvetica"    
+                  lgres@lgLabelFontHeightF = .1            
+                  lgres@vpWidthF           = 0.4           
+                  lgres@vpHeightF          = 0.4          
+                  lgres@lgDashIndexes      = (/0,0,0/)
+                  lgres@lgLineColors       = (/237,144,80/)
+                  lbid = gsn_create_legend(wks_ps,3,(/"splptx","splpty","splptz"/),lgres)       
+
+                  amres = True
+                  amres@amParallelPosF   = 0.6		    
+                  amres@amOrthogonalPosF = -0.2            
+                  annoid1 = gsn_add_annotation(plot_ps(n),lbid,amres)
                end if       
             end if
@@ -840,5 +1196,5 @@
    ; merge plots onto one page 
    ; ***************************************************
-  
+ 
    if (format_out .EQ. "eps" .OR. format_out .EQ. "epsi") then
       gsn_panel(wks_ps,plot_ps(1:n),(/n,1/),resP)
