Index: palm/trunk/SCRIPTS/NCL/cross_sections.ncl
===================================================================
--- palm/trunk/SCRIPTS/NCL/cross_sections.ncl	(revision 195)
+++ palm/trunk/SCRIPTS/NCL/cross_sections.ncl	(revision 218)
@@ -37,5 +37,5 @@
    if (file_1 .EQ. "File in") then
       print(" ")
-      print("Please provide 1st input file 'file_1=' in 'ncl_preferences.ncl'")
+      print("Declare input file 'file_1=' in 'ncl_preferences.ncl' or prompt")
       print(" ")
       exit 
@@ -43,10 +43,4 @@
       file_in = file_1   
    end if
-   if (.not. isfilepresent(file_in)) then
-      print(" ")
-      print("1st input file: '"+file_in+"' does not exist")
-      print(" ")
-      exit
-   end if 
 
    if (format_out .NE. "x11" .AND. format_out .NE. "pdf" .AND. format_out .NE. "eps" .AND. format_out .NE. "ps" .AND. format_out .NE. "epsi" .AND. format_out .NE. "ncgm")then
@@ -73,5 +67,5 @@
    if (fill_mode .NE. "AreaFill" .AND. fill_mode .NE. "RasterFill" .AND. fill_mode .NE. "CellFill") then
       print(" ")
-      print("Your 'fill_mode'= "+fill_mode+" is invalid and set to 'AreaFill'")
+      print("'fill_mode'= "+fill_mode+" is invalid and set to 'AreaFill'")
       print(" ")
       fill_mode = "AreaFill"
@@ -108,5 +102,5 @@
    if (xyc .EQ. 0 .AND. xzc .EQ. 0 .AND. yzc .EQ. 0) then
       print(" ")
-      print("Please select one crossection (xyc=1 or xzc=1 or yzc=1)")
+      print("Select one crossection (xyc=1 or xzc=1 or yzc=1)")
       print(" ")
       exit
@@ -115,5 +109,5 @@
       if (xzc .EQ. 1 .OR. yzc .EQ. 1) then
          print(" ")
-         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
+         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
          print(" ") 
          exit
@@ -123,5 +117,5 @@
       if (xyc .EQ. 1 .OR. yzc .EQ. 1) then
          print(" ")
-         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
+         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
          print(" ")
          exit
@@ -131,5 +125,5 @@
       if (xyc .EQ. 1 .OR. xzc .EQ. 1) then
          print(" ")
-         print("Please select just one crossection (xyc=1 or xzc=1 or yzc=1)")
+         print("Select just one crossection (xyc=1 or xzc=1 or yzc=1)")
          print(" ")
          exit
@@ -139,16 +133,102 @@
    if (vector .NE. 0 .AND. vector .NE. 1) then
       print(" ")
-      print("Please set 'vector' to 0 or 1")
+      print("Set 'vector' to 0 or 1")
       print(" ")
       exit
-   end if   
+   end if  
+   
+   if (norm_x .EQ. 0) then 
+      print(" ")
+      print("Normalising with 0 is not allowed, 'norm_x' is set to 1.0")
+      print(" ")
+      norm_x = 1.0
+   end if
+   if (norm_y .EQ. 0) then 
+      print(" ")
+      print("Normalising with 0 is not allowed, 'norm_y' is set to 1.0")
+      print(" ")
+      norm_y = 1.0
+   end if 
+   if (norm_z .EQ. 0) then 
+      print(" ")
+      print("Normalising with 0 is not allowed, 'norm_z' is set to 1.0")
+      print(" ")
+      norm_z = 1.0
+   end if 
  
    ; ***************************************************
    ; open input file
    ; ***************************************************
-
-   f  = addfile( file_in, "r" )
-
-   vNam  = getfilevarnames(f) 
+   
+   file_in_1 = False
+   if (isStrSubset(file_in, ".nc"))then
+      start_f = -2
+      end_f = -2 
+      file_in_1 = True     
+   end if
+
+   if (start_f .EQ. -1)then
+      print(" ")
+      print("'start_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
+      print(" ")  
+      exit
+   end if 
+   if (end_f .EQ. -1)then
+      print(" ")
+      print("'end_f' must be one of the cyclic numbers (at least 0) of your input file(s)")
+      print(" ")  
+      exit
+   end if 
+
+   files=new(end_f-start_f+1,string)
+
+   if (file_in_1)then
+      if (isfilepresent(file_in))then
+         files(0)=file_in
+      else
+         print(" ")
+         print("1st input file: '"+file_in+"' does not exist")
+         print(" ")
+         exit
+      end if
+   else   
+      if (start_f .EQ. 0)then
+         if (isfilepresent(file_in+".nc"))then
+            files(0)=file_in+".nc"
+            do i=1,end_f
+	       if (isfilepresent(file_in+"."+i+".nc"))then    
+                  files(i)=file_in+"."+i+".nc"
+	       else
+	          print(" ")
+                  print("Input file: '"+file_in+"."+i+".nc' does not exist")
+                  print(" ")
+                  exit  
+	       end if	    
+            end do         
+         else
+            print(" ")
+            print("Input file: '"+file_in+".nc' does not exist")
+            print(" ")
+            exit
+         end if
+      else
+         do i=start_f,end_f
+            if (isfilepresent(file_in+"."+i+".nc"))then    
+               files(i-start_f)=file_in+"."+i+".nc"
+	    else
+	       print(" ")
+               print("Input file: '"+file_in+"."+i+".nc' does not exist")
+               print(" ")
+               exit  
+	    end if
+         end do
+      end if
+   end if
+   
+   f=addfiles(files,"r")
+   f_att=addfile(files(0),"r")
+   ListSetType(f,"cat")
+
+   vNam  = getfilevarnames(f_att) 
    print(" ")
    print("Variables in input file:")
@@ -225,6 +305,5 @@
    else
       print(" ")
-      print("Your input file: '"+file_in+"'")
-      print("contains 3d or other data")
+      print("Your input file contains 3d or other data")
       print(" ")
    end if
@@ -242,4 +321,5 @@
 
    cs_res                          = True
+   vecres                          = True
    cs_res@gsnYAxisIrregular2Linear = True  
   
@@ -247,30 +327,26 @@
    cs_res@gsnFrame                = False
    cs_res@gsnMaximize             = True 
-   cs_res@gsnPaperOrientation     = "portrait"
-   cs_res@gsnPaperWidth           = 8.27
-   cs_res@gsnPaperHeight          = 11.69
-   cs_res@gsnPaperMargin          = 0.79
-   cs_res@tmXBLabelFontHeightF   = .03
-   cs_res@tmYLLabelFontHeightF   = .03
-   cs_res@tiXAxisFontHeightF     = .03
-   cs_res@tiYAxisFontHeightF     = .03
+   
+   cs_res@tmXBLabelFontHeightF   = font_size
+   cs_res@tmYLLabelFontHeightF   = font_size
+   cs_res@tiXAxisFontHeightF     = font_size
+   cs_res@tiYAxisFontHeightF     = font_size
    cs_res@tmXBMode                ="Automatic"
    cs_res@tmYLMode                ="Automatic"
-   cs_res@lgTitleFontHeightF     = .03 
-   cs_res@lgLabelFontHeightF     = .03
-   cs_res@txFontHeightF          = .03
+  
    cs_res@cnLevelSelectionMode    = "ManualLevels"
+   cs_res@lbLabelFontHeightF = font_size_legend
+   cs_res@lbLabelStride = legend_label_stride
+   
 
    cs_resP = True
-   cs_resP@txString               = f@title
+   cs_resP@txString               = f_att@title
    cs_resP@txFuncCode             = "~"         	
-   cs_resP@txFontHeightF          = .02
-   
+   cs_resP@txFontHeightF          = 0.015        
   
    if ( mode .eq. "Fill" ) then
       cs_res@cnFillOn             = True
       cs_res@gsnSpreadColors      = True
-      cs_res@cnFillMode           = fill_mode
-      cs_res@lbOrientation        = "Vertical"  	
+      cs_res@cnFillMode           = fill_mode	
       cs_res@cnLinesOn            = False	
       cs_res@cnLineLabelsOn       = False
@@ -280,6 +356,5 @@
       cs_res@cnFillOn             = True
       cs_res@gsnSpreadColors      = True
-      cs_res@cnFillMode           = fill_mode
-      cs_res@lbOrientation        = "Vertical"  
+      cs_res@cnFillMode           = fill_mode 
       cs_res@cnLinesOn            = True
       cs_res@cnLineLabelsOn       = True
@@ -290,5 +365,5 @@
    ; *********************************************
 
-   t_all = f->time
+   t_all = f[:]->time
    nt = dimsizes(t_all)
    delta_t = t_all(nt-1)/nt
@@ -305,5 +380,5 @@
          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("Select another 'start_time_step'")
          print(" ")
          exit
@@ -313,10 +388,10 @@
          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("Select another 'start_time_step'")
          print(" ")
          exit
       end if
    end if
-   do i=0,nt-1     
+   do i=0,nt-1    
       if (start_time_step .GE. (t_all(i)-delta_t/2)/3600 .AND. start_time_step .LT. (t_all(i)+delta_t/2)/3600)then
          st=i
@@ -329,5 +404,5 @@
       print("'start_time_step' = "+ start_time_step +"h is invalid")
       print(" ")
-      print("Please select another 'start_time_step'")
+      print("Select another 'start_time_step'")
       print(" ")
       exit
@@ -345,5 +420,5 @@
          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("Select another 'end_time_step'")
          print(" ")
          exit
@@ -353,5 +428,5 @@
          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("Select another 'end_time_step'")  
          print(" ")
          exit
@@ -361,5 +436,5 @@
          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("Select another 'start_time_step' or 'end_time_step'")
          print(" ")
          exit
@@ -377,5 +452,5 @@
       print("'end_time_step' = "+ end_time_step +"h is invalid")
       print(" ")
-      print("Please select another 'end_time_step'")
+      print("Select another 'end_time_step'")
       print(" ")
       exit
@@ -401,5 +476,5 @@
       if (vec1 .EQ. "component1") then
          print(" ")
-         print("Please indicate Vector 1 ('vec1') for Vector-Plot or set 'vector' to 0")
+         print("Indicate Vector 1 ('vec1') for Vector-Plot or set 'vector' to 0")
          print(" ")
          exit
@@ -407,5 +482,5 @@
       if (vec2 .EQ. "component2") then
          print(" ")
-         print("Please indicate Vector 2 ('vec2') for Vector-Plot or set 'vector' to 0")
+         print("Indicate Vector 2 ('vec2') for Vector-Plot or set 'vector' to 0")
          print(" ")
          exit
@@ -424,5 +499,5 @@
       do varn=0,dim-1
          if (vNam(varn) .eq. "xu" .OR. vNam(varn) .eq. "x")then    
-            x_d     = f->$vNam(varn)$
+            x_d     = f_att->$vNam(varn)$
             xdim    = dimsizes(x_d)-1
             delta_x = x_d(1)-x_d(0)
@@ -432,5 +507,5 @@
       do varn=0,dim-1
          if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y")then   
-            y_d     = f->$vNam(varn)$
+            y_d     = f_att->$vNam(varn)$
             ydim    = dimsizes(y_d)-1
             delta_y = y_d(1)-y_d(0)
@@ -440,5 +515,5 @@
       do varn=0,dim-1
          if (vNam(varn) .eq. "zu_3d" .OR. vNam(varn) .eq. "zw_3d")then
-            z_d     = f->$vNam(varn)$
+            z_d     = f_att->$vNam(varn)$
             zdim    = dimsizes(z_d)-1
             delta_z = 0
@@ -446,5 +521,5 @@
          else
             if (vNam(varn) .eq. "zu_xy" .OR. vNam(varn) .eq. "zw_xy") then
-               z_d     = f->$vNam(varn)$
+               z_d     = f_att->$vNam(varn)$
                zdim    = dimsizes(z_d)-1 
                delta_z = -1.d
@@ -460,5 +535,5 @@
       do varn=0,dim-1
          if (vNam(varn) .eq. "xu"  .OR. vNam(varn) .eq. "x") then 
-            x_d     = f->$vNam(varn)$
+            x_d     = f_att->$vNam(varn)$
             xdim    = dimsizes(x_d)-1
             delta_x = x_d(1)-x_d(0)
@@ -468,5 +543,5 @@
       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)$
+            z_d     = f_att->$vNam(varn)$
             zdim    = dimsizes(z_d)-1 
             delta_z = z_d(1)-z_d(0)
@@ -476,5 +551,5 @@
       do varn=0,dim-1 
          if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
-            y_d     = f->$vNam(varn)$
+            y_d     = f_att->$vNam(varn)$
             ydim    = dimsizes(y_d)-1
             delta_y = y_d(1)-y_d(0)
@@ -482,5 +557,5 @@
          else
             if (vNam(varn) .eq. "y_xz" .or. vNam(varn) .eq. "yv_xz") then
-               y_d     = f->$vNam(varn)$
+               y_d     = f_att->$vNam(varn)$
                ydim    = dimsizes(y_d)-1 
                delta_y = -1.d
@@ -496,5 +571,5 @@
       do varn=0,dim-1  
          if (vNam(varn) .eq. "yv" .or. vNam(varn) .eq. "y") then
-            y_d     = f->$vNam(varn)$
+            y_d     = f_att->$vNam(varn)$
             ydim    = dimsizes(y_d)-1
             delta_y = y_d(1)-y_d(0)
@@ -504,5 +579,5 @@
       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)$
+            z_d     = f_att->$vNam(varn)$
             zdim    = dimsizes(z_d)-1 
             delta_z = z_d(1)-z_d(0)
@@ -512,5 +587,5 @@
       do varn=0,dim-1 
          if (vNam(varn) .eq. "xu" .or. vNam(varn) .eq. "x")
-            x_d     = f->$vNam(varn)$
+            x_d     = f_att->$vNam(varn)$
             xdim    = dimsizes(x_d)-1 
             delta_x = x_d(1)-x_d(0)
@@ -518,5 +593,5 @@
          else
             if (vNam(varn) .eq. "xu_yz" .or. vNam(varn) .eq. "x_yz" ) then
-               x_d     = f->$vNam(varn)$
+               x_d     = f_att->$vNam(varn)$
                xdim    = dimsizes(x_d)-1 
                delta_x = -1.d
@@ -848,20 +923,53 @@
       end do 
    end if
-   
+  
    if( shape .eq. 1 ) then 
       if (xyc .EQ. 1)then
          cs_res@vpWidthF    = (xe-xs)/(ye-ys)        
          cs_res@vpHeightF   = 1
+	 vecres@vpWidthF    = (xe-xs)/(ye-ys)        
+         vecres@vpHeightF   = 1
+	 if (xe-xs .GT. ye-ys)then
+	    cs_res@gsnPaperOrientation     = "landscape"
+	    vecres@gsnPaperOrientation     = "landscape"
+	    cs_res@lbOrientation = "Horizontal"   
+	 else
+	    cs_res@gsnPaperOrientation     = "portrait"
+	    vecres@gsnPaperOrientation     = "portrait"
+	    cs_res@lbOrientation = "Vertical" 
+	 end if
       end if
       if (xzc .EQ. 1)then
          cs_res@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))        
          cs_res@vpHeightF   = 1
+	 vecres@vpWidthF    = (xe-xs)/(delta_x*(ze-zs))        
+         vecres@vpHeightF   = 1
+	 if (xe-xs .GT. (delta_x*(ze-zs)))then
+	    cs_res@gsnPaperOrientation     = "landscape"
+	    vecres@gsnPaperOrientation     = "landscape"
+	    cs_res@lbOrientation = "Horizontal"   
+	 else
+	    cs_res@gsnPaperOrientation     = "portrait"
+	    vecres@gsnPaperOrientation     = "portrait"
+	    cs_res@lbOrientation = "Vertical" 
+	 end if
       end if
       if (yzc .EQ. 1)then
          cs_res@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))        
          cs_res@vpHeightF   = 1
-      end if
-   end if
-   
+	 vecres@vpWidthF    = (ye-ys)/(delta_y*(ze-zs))        
+         vecres@vpHeightF   = 1
+	 if (ye-ys .GT. (delta_y*(ze-zs)))then
+            cs_res@gsnPaperOrientation     = "landscape"
+	    vecres@gsnPaperOrientation     = "landscape"
+	    cs_res@lbOrientation = "Horizontal"   
+	 else
+	    cs_res@gsnPaperOrientation     = "portrait"
+	    vecres@gsnPaperOrientation     = "portrait"
+	    cs_res@lbOrientation = "Vertical" 
+	 end if
+      end if
+   end if
+
    delete(xs)
    delete(xe)   
@@ -874,4 +982,178 @@
    ye=yend
 
+   if (xyc .EQ. 1) then
+      d=(ye-ys+1)/(major_ticks_y-1)
+      e=(xe-xs+1)/(major_ticks_x-1)
+      array_yl       =new(major_ticks_y,integer)
+      array_yl_labels=new(major_ticks_y,double)
+      array_minor_yl =new((major_ticks_y-1)*4,double)
+      array_xb       =new(major_ticks_x,integer)
+      array_xb_labels=new(major_ticks_x,double)
+      array_minor_xb =new((major_ticks_x-1)*4,double)
+      array_yl(0)=ys
+      array_xb(0)=xs
+      array_yl_labels(0)=y_d(ys)
+      array_xb_labels(0)=x_d(xs)
+      do ar=1,major_ticks_y-1
+         array_yl(ar)=d*(ar-1)+d-1
+	 array_yl_labels(ar) = y_d(array_yl(ar))
+	 if (ar .GT. 0)
+	    do min_ar=0,3
+	       array_minor_yl(4*(ar-1)+min_ar)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-array_yl(ar-1))/5*(min_ar+1)
+	    end do 
+	 end if	 
+      end do 
+      do br=1,major_ticks_x-1
+         array_xb(br)=e*(br-1)+e-1 
+	 array_xb_labels(br) = x_d(array_xb(br)) 
+	 if (br .GT. 0)
+	    do min_br=0,3
+	       array_minor_xb(4*(br-1)+min_br)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-array_xb(br-1))/5*(min_br+1)
+	    end do 
+	 end if	
+      end do 
+   end if
+ 
+   if (xzc .EQ. 1) then
+      d=(ze-zs+1)/(major_ticks_y-1)
+      e=(xe-xs+1)/(major_ticks_x-1)
+      array_yl       =new(major_ticks_y,integer)
+      array_yl_labels=new(major_ticks_y,double)
+      array_minor_yl =new((major_ticks_y-1)*4,double)
+      array_xb       =new(major_ticks_x,integer)
+      array_xb_labels=new(major_ticks_x,double)
+      array_minor_xb =new((major_ticks_x-1)*4,double)
+      array_yl(0)=zs
+      array_xb(0)=xs
+      array_yl_labels(0)=z_d(zs)
+      array_xb_labels(0)=x_d(xs)
+      do ar=1,major_ticks_y-1
+         array_yl(ar)=d*(ar-1)+d-1
+	 array_yl_labels(ar) = z_d(array_yl(ar))
+	 if (ar .GT. 0)
+	    do min_ar=0,3
+	       array_minor_yl(4*(ar-1)+min_ar)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-array_yl(ar-1))/5*(min_ar+1)
+	    end do 
+	 end if	 
+      end do 
+      do br=1,major_ticks_x-1
+         array_xb(br)=e*(br-1)+e-1 
+	 array_xb_labels(br) = x_d(array_xb(br)) 
+	 if (br .GT. 0)
+	    do min_br=0,3
+	       array_minor_xb(4*(br-1)+min_br)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-array_xb(br-1))/5*(min_br+1)
+	    end do 
+	 end if	
+      end do 
+   end if
+
+   if (yzc .EQ. 1) then
+      d=(ze-zs+1)/(major_ticks_y-1)
+      e=(ye-ys+1)/(major_ticks_x-1)
+      array_yl       =new(major_ticks_y,integer)
+      array_yl_labels=new(major_ticks_y,double)
+      array_minor_yl =new((major_ticks_y-1)*4,double)
+      array_xb       =new(major_ticks_x,integer)
+      array_xb_labels=new(major_ticks_x,double)
+      array_minor_xb =new((major_ticks_x-1)*4,double)
+      array_yl(0)=zs
+      array_xb(0)=ys
+      array_yl_labels(0)=z_d(zs)
+      array_xb_labels(0)=y_d(ys)
+      do ar=1,major_ticks_y-1
+         array_yl(ar)=d*(ar-1)+d-1
+	 array_yl_labels(ar) = y_d(array_yl(ar))
+	 if (ar .GT. 0)
+	    do min_ar=0,3
+	       array_minor_yl(4*(ar-1)+min_ar)= int2dble(array_yl(ar-1))+int2dble(array_yl(ar)-array_yl(ar-1))/5*(min_ar+1)
+	    end do 
+	 end if	 
+      end do 
+      do br=1,major_ticks_x-1
+         array_xb(br)=e*(br-1)+e-1 
+	 array_xb_labels(br) = x_d(array_xb(br)) 
+	 if (br .GT. 0)
+	    do min_br=0,3
+	       array_minor_xb(4*(br-1)+min_br)= int2dble(array_xb(br-1))+int2dble(array_xb(br)-array_xb(br-1))/5*(min_br+1)
+	    end do 
+	 end if	
+      end do
+   end if
+
+   if (axes_explicit .EQ. 1)then
+      cs_res@tmYLMode = "Explicit"
+      cs_res@tmXBMode = "Explicit"
+      cs_res@tmYLValues = array_yl
+      if (xyc .EQ. 1)then
+         cs_res@tmYLLabels = array_yl_labels/norm_y
+	 cs_res@tmXBLabels = array_xb_labels/norm_x
+	 if (norm_x .NE. 1.)then
+	    cs_res@tiXAxisString = "x ["+unit_x+"]"
+	 else
+	    cs_res@tiXAxisString = "x [m]"
+	 end if 
+	 if (norm_y .NE. 1.)then
+            cs_res@tiYAxisString = "y ["+unit_y+"]"
+	 else
+	    cs_res@tiYAxisString = "y [m]"
+	 end if   
+      end if
+      if (xzc .EQ. 1)then
+         cs_res@tmYLLabels = array_yl_labels/norm_z
+	 cs_res@tmXBLabels = array_xb_labels/norm_x
+	 if (norm_x .NE. 1.)then
+	    cs_res@tiXAxisString = "x ["+unit_x+"]"
+	 else
+	    cs_res@tiXAxisString = "x [m]"
+	 end if 
+	 if (norm_z .NE. 1.)then
+            cs_res@tiYAxisString = "z ["+unit_z+"]"
+	 else
+	    cs_res@tiYAxisString = "z [m]"
+	 end if
+      end if
+      if (yzc .EQ. 1)then
+         cs_res@tmYLLabels = array_yl_labels/norm_z
+	 cs_res@tmXBLabels = array_xb_labels/norm_y
+	 if (norm_y .NE. 1.)then
+	    cs_res@tiXAxisString = "y ["+unit_y+"]"
+	 else
+	    cs_res@tiXAxisString = "y [m]"
+	 end if 
+	 if (norm_z .NE. 1.)then
+            cs_res@tiYAxisString = "z ["+unit_z+"]"
+	 else
+	    cs_res@tiYAxisString = "z [m]"
+	 end if
+      end if
+      cs_res@tmXBValues = array_xb      
+      cs_res@tmYLMinorValues = array_minor_yl
+      cs_res@tmXBMinorValues = array_minor_xb
+   else 
+      if (axes_explicit .EQ. 0)then  
+         cs_res@tmYLMinorPerMajor = 4
+         cs_res@tmXBMinorPerMajor = 4
+      else
+         print(" ")
+	 print("'axes_explicit' is invalid and set to 0")
+	 print(" ")
+	 axes_explicit = 0
+	 cs_res@tmYLMinorPerMajor = 4
+         cs_res@tmXBMinorPerMajor = 4
+      end if
+      if (xyc .EQ. 1)then
+         cs_res@tiXAxisString = "x [gridpoints]"
+         cs_res@tiYAxisString = "y [gridpoints]"
+      end if
+      if (xzc .EQ. 1)then
+         cs_res@tiXAxisString = "x [gridpoints]"
+         cs_res@tiYAxisString = "z [gridpoints]"
+      end if
+      if (yzc .EQ. 1)then
+         cs_res@tiXAxisString = "y [gridpoints]"
+         cs_res@tiYAxisString = "z [gridpoints]"
+      end if
+   end if
+
    if (xyc .EQ. 1 .OR. xzc .EQ.1)then
       if (xe .EQ. xs+1) then
@@ -916,5 +1198,6 @@
 
    MinVal = new(dim,float)
-   MaxVal = new(dim,float) 
+   MaxVal = new(dim,float)
+   unit   = new(dim,string) 
 
    ; ****************************************************
@@ -980,11 +1263,17 @@
 
          if (xyc .EQ. 1) then
-            data(varn,:,:,:,:)=f->$vNam(varn)$(:,zs:ze,ys:ye,xs:xe) 	
+	    temp = f[:]->$vNam(varn)$
+	    data_att = f_att->$vNam(varn)$
+            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe) 	    	
          end if
          if ( xzc .eq. 1 ) then
-            data(varn,:,:,:,:)=f->$vNam(varn)$(:,zs:ze,ys:ye,xs:xe) 
+	    temp = f[:]->$vNam(varn)$
+	    data_att = f_att->$vNam(varn)$
+            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe) 
          end if
          if ( yzc .eq. 1) then
-            data(varn,:,:,:,:)=f->$vNam(varn)$(:,zs:ze,ys:ye,xs:xe) 
+	    temp = f[:]->$vNam(varn)$
+	    data_att = f_att->$vNam(varn)$
+            data(varn,:,:,:,:)=temp(:,zs:ze,ys:ye,xs:xe) 
          end if
          if (check_vec1) then
@@ -1003,4 +1292,6 @@
          MinVal(varn) = min(data(varn,:,:,:,:))
          MaxVal(varn) = max(data(varn,:,:,:,:))
+	 
+	 unit(varn) = data_att@units
 
       end if
@@ -1075,6 +1366,10 @@
    gsn_define_colormap(wks_ps,"rainbow+white")
 
-   plot=new((/no_time*no_layer*no_var/),graphic)
- 
+   if (vector .EQ. 1 .AND. plotvec .EQ. "plotvec") then
+      plot=new((/no_time*no_layer*no_var+no_time*no_layer/),graphic)
+   else
+      plot=new((/no_time*no_layer*no_var/),graphic)
+   end if
+
    page = 0
    n=0
@@ -1090,5 +1385,25 @@
       do lo = los, loe           				
          do li = lis, lie
-            level = "=" + data&z(lo) + "m" 	  
+	    if (xyc .EQ. 1)then
+	       if (sort .EQ. "time")then
+                  level = "z=" + z_d(lo) + "m"
+	       else
+	          level = "z=" + z_d(li) + "m"
+	       end if
+	    end if
+	    if (xzc .EQ. 1)then
+	       if (sort .EQ. "time")then
+                  level = "y=" + y_d(lo) + "m" 
+	       else
+	          level = "y=" + y_d(li) + "m"
+	       end if
+	    end if
+	    if (yzc .EQ. 1)then
+	       if (sort .EQ. "time")then
+                  level = "x=" + x_d(lo) + "m" 
+	       else
+	          level = "x=" + x_d(li) + "m" 
+	       end if
+	    end if	    	  
             vecres                  = True            ; vector only resources
             vecres@gsnDraw          = False           ; don't draw
@@ -1098,7 +1413,35 @@
   	    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)
+	    vecres@tmXBLabelFontHeightF   = font_size
+            vecres@tmYLLabelFontHeightF   = font_size
+            vecres@tiXAxisFontHeightF     = font_size
+            vecres@tiYAxisFontHeightF     = font_size
+	    if (sort .EQ. "time")then
+               vecres@gsnLeftString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h   "+level	
+	    else
+	       vecres@gsnLeftString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) +"h   "+level	
+	    end if
+	    vecres@tiXAxisString    = " " 
+	    if (xyc .EQ. 1)then  
+	       if (sort .EQ. "time")then                                          
+                  plot(n) = gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+	       else
+	          plot(n) = gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+	       end if
+	    end if
+	    if (xzc .EQ. 1) then
+	       if (sort .EQ. "time")then 
+	          plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
+	       else
+	          plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
+	       end if
+	    end if
+	    if (yzc .EQ. 1) then
+	       if (sort .EQ. "time")then 
+	          plot(n) = gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
+	       else
+	          plot(n) = gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
+	       end if
+	    end if
             n=n+1
          end do
@@ -1123,5 +1466,5 @@
     
       if(check) then
-   
+
          space=(MaxVal(varn)-MinVal(varn))/24
   
@@ -1144,18 +1487,17 @@
 
                if ( xyc .eq. 1 ) then 
-                
-                  cs_res@tiXAxisString = "x [m]"
-                  cs_res@tiYAxisString = "y [m]"
+          
                   cs_res@gsnLeftString = "Plot of "+vNam(varn)
+		  cs_res@gsnRightString = unit(varn)
                   
                   if ( sort .eq. "time" ) then
-                     if ( data&z(lo) .eq. -1)then
+                     if ( z_d(lo) .eq. -1)then
 		        if (delta_z .EQ. -1) then
                            level = "-average"
 			else
-			   level = "=" + data&z(lo) + "m"
+			   level = "=" + z_d(lo) + "m"
 			end if
                      else 
-                        level = "=" + data&z(lo) + "m"
+                        level = "=" + z_d(lo) + "m"
                      end if
                      cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) +"h  z"+level		     
@@ -1177,12 +1519,12 @@
                  
                   if ( sort .eq. "layer" ) then
-		     if (data&z(li) .eq. -1) then
+		     if (z_d(li) .eq. -1) then
 		        if (delta_z .EQ. -1) then
                            level = "-average"
 			else
-			   level = "=" + data&z(li) + "m"
+			   level = "=" + z_d(li) + "m"
 			end if
                      else 
-                        level = "=" + data&z(li) + "m"
+                        level = "=" + z_d(li) + "m"
                      end if 
                      cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  z"+ level		
@@ -1209,16 +1551,14 @@
 
                if ( xzc .eq. 1 ) then
-        
-                  cs_res@tiXAxisString   = "x [m]"
-                  cs_res@tiYAxisString   = "z [m]"
+                  
 		  cs_res@gsnLeftString = "Plot of "+vNam(varn)
                
                   if ( sort .eq. "time" ) then
-                     if ( data&y(lo) .eq. -1 ) then
+                     if ( y_d(lo) .eq. -1 ) then
                         level = "-average"
                      else 
-                        level = "=" + data&y(lo) + "m"
+                        level = "=" + y_d(lo) + "m"
                      end if
-                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "s  y"+ level
+                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  y"+ level
                      plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,lo,:),cs_res)
                      if (vector .EQ. 1 .AND. check_vecp) then
@@ -1232,5 +1572,5 @@
   			vecres@gsnLeftString    = " "             ; turn off left string
   			vecres@tiXAxisString    = " "   
-                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,lo,:),vect2(li,:,lo,:),vecres)
                         overlay(plot(n), plot_vec)
                      end if
@@ -1238,10 +1578,10 @@
 
                   if ( sort .eq. "layer" ) then
-                     if ( data&y(li) .eq. -1 ) then
+                     if ( y_d(li) .eq. -1 ) then
                         level = "-average"
                      else 
-                        level = "=" + data&y(li) + "m"
+                        level = "=" + y_d(li) + "m"
                      end if
-                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "s  y"+ level
+                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  y"+ level
                      plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,li,:),cs_res)
                      if (vector .EQ. 1 .AND. check_vecp) then
@@ -1255,5 +1595,5 @@
   			vecres@gsnLeftString    = " "             ; turn off left string
   			vecres@tiXAxisString    = " "   
-                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,li,:),vect2(lo,:,li,:),vecres)
                         overlay(plot(n), plot_vec)
                      end if
@@ -1266,16 +1606,14 @@
 
                if ( yzc .eq. 1 ) then
-               
-                  cs_res@tiXAxisString   = "y [m]"
-                  cs_res@tiYAxisString   = "z [m]"
+                                  
 		  cs_res@gsnLeftString = "Plot of "+vNam(varn)
 
                   if ( sort .eq. "time" ) then
-                     if ( data&x(lo) .eq. -1 ) then
+                     if ( x_d(lo) .eq. -1 ) then
                         level = "-average"
                      else 
-                        level = "=" + data&x(lo) + "m"
+                        level = "=" + x_d(lo) + "m"
                      end if
-                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "s  x"+ level
+                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(li)/3600,2,True) + "h  x"+ level
                      plot(n) = gsn_csm_contour(wks_ps,data(varn,li,:,:,lo),cs_res)
     	             if (vector .EQ. 1 .AND. check_vecp) then
@@ -1289,5 +1627,5 @@
   			vecres@gsnLeftString    = " "             ; turn off left string
   			vecres@tiXAxisString    = " "   
-                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,lo,:,:),vect2(li,lo,:,:),vecres)
+                        plot_vec=gsn_csm_vector(wks_ps,vect1(li,:,:,lo),vect2(li,:,:,lo),vecres)
                         overlay(plot(n), plot_vec)
                      end if
@@ -1295,10 +1633,10 @@
 
                   if ( sort .eq. "layer" ) then
-                     if ( data&x(li) .eq. -1 ) then
+                     if ( x_d(li) .eq. -1 ) then
                         level = "-average"
                      else 
-                        level = "=" + data&x(li) + "m"
+                        level = "=" + x_d(li) + "m"
                      end if
-                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "s  x"+ level
+                     cs_res@gsnCenterString = "t=" + decimalPlaces(t_all(lo)/3600,2,True) + "h  x"+ level
                      plot(n) = gsn_csm_contour(wks_ps,data(varn,lo,:,:,li),cs_res)
                      if (vector .EQ. 1 .AND. check_vecp)then
@@ -1312,10 +1650,10 @@
   			vecres@gsnLeftString    = " "             ; turn off left string
   			vecres@tiXAxisString    = " "   
-                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,li,:,:),vect2(lo,li,:,:),vecres)
+                        plot_vec=gsn_csm_vector(wks_ps,vect1(lo,:,:,li),vect2(lo,:,:,li),vecres)
                         overlay(plot(n), plot_vec)
                      end if 
                   end if
                end if	      
-               n=n+1   
+               n=n+1  
             end do 	
          end do  
@@ -1334,13 +1672,13 @@
          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)
+         do np = 0,no_layer*no_time*(no_var+1)-1,no_rows*no_columns   
+            if ( np + no_rows*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_rows,no_columns/),cs_resP)
             else
-               gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
+               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
             end if
          end do
       end if
-   else          
+   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)
@@ -1349,9 +1687,9 @@
          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)
+         do np = 0,no_layer*no_time*no_var-1,no_rows*no_columns   
+            if ( np + no_rows*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_rows,no_columns/),cs_resP)
             else
-               gsn_panel(wks_ps, plot(np:np+no_lines*no_columns-1),(/no_lines,no_columns/),cs_resP)
+               gsn_panel(wks_ps, plot(np:np+no_rows*no_columns-1),(/no_rows,no_columns/),cs_resP)
             end if
          end do
