diff --git a/matlab/generate_grid.m b/matlab/generate_grid.m
index b00f5e8..47e4842 100644
--- a/matlab/generate_grid.m
+++ b/matlab/generate_grid.m
@@ -1,8 +1,16 @@
 function [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
-          xu1,yu1,zu1,xv1,yv1,zv1,xw1,yw1,zw1,xp1,yp1,zp1,dx,dy,dz]=...
-            generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,x_periodic,y_periodic,z_periodic)
+          xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,dx,dy,dz]=...
+            generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,xperiodic,yperiodic,zperiodic)
+  % [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
+  %  xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,dx,dy,dz]=...
+  %          generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,xperiodic,yperiodic,zperiodic)
+  % Generates a staggered grid.
+  % Input
+  %   a,b,c,d,e,f       domain bounds
+  %   nxp,nyp,nzp       number of points for pressure grid
+  %   x/y/zperiodic     domain periodicity
   
-  if x_periodic==1
+  if xperiodic==1
     nxu=nxp;
     nxv=nxp;
     nxw=nxp;
@@ -11,7 +19,7 @@ function [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
     nxv=nxp;
     nxw=nxp;
   end
-  if y_periodic==1
+  if yperiodic==1
     nyu=nyp;
     nyv=nyp;
     nyw=nyp;
@@ -20,7 +28,7 @@ function [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
     nyv=nyp+1;
     nyw=nyp;
   end
-  if z_periodic==1
+  if zperiodic==1
     nzu=nzp;
     nzv=nzp;
     nzw=nzp;
@@ -30,83 +38,83 @@ function [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
     nzw=nzp+1;
   end
   %grid step:
-  if x_periodic==1
+  if xperiodic==1
     dx=(b-a)/(nxp);
   else
     dx=(b-a)/(nxp-1);
   end
-  if y_periodic==1
+  if yperiodic==1
     dy=(d-c)/(nyp);
   else
     dy=(d-c)/(nyp-1);
   end
-  if z_periodic==1
+  if zperiodic==1
     dz=(f-e)/(nzp);
   else
     dz=(f-e)/(nzp-1);
   end
   %
-  if x_periodic==1
-    xu1=a+((1:nxu)-1)*dx;
+  if xperiodic==1
+    xu=a+((1:nxu)-1)*dx;
   else
-    xu1=a+((1:nxu)-3/2)*dx;
+    xu=a+((1:nxu)-3/2)*dx;
   end
-  if y_periodic==1
-    yu1=c+((1:nyu)-1/2)*dy;
+  if yperiodic==1
+    yu=c+((1:nyu)-1/2)*dy;
   else
-    yu1=c+((1:nyu)-1)*dy;
+    yu=c+((1:nyu)-1)*dy;
   end
-  if z_periodic==1
-    zu1=e+((1:nzu)-1/2)*dz;
+  if zperiodic==1
+    zu=e+((1:nzu)-1/2)*dz;
   else
-    zu1=e+((1:nzu)-1)*dz;
+    zu=e+((1:nzu)-1)*dz;
   end
   %
-  if x_periodic==1
-    xv1=a+((1:nxv)-1/2)*dx;
+  if xperiodic==1
+    xv=a+((1:nxv)-1/2)*dx;
   else
-    xv1=a+((1:nxv)-1)*dx;
+    xv=a+((1:nxv)-1)*dx;
   end
-  if y_periodic==1
-    yv1=c+((1:nyv)-1)*dy;
+  if yperiodic==1
+    yv=c+((1:nyv)-1)*dy;
   else
-    yv1=c+((1:nyv)-3/2)*dy;
+    yv=c+((1:nyv)-3/2)*dy;
   end
-  if z_periodic==1
-    zv1=e+((1:nzv)-1/2)*dz;
+  if zperiodic==1
+    zv=e+((1:nzv)-1/2)*dz;
   else
-    zv1=e+((1:nzv)-1)*dz;
+    zv=e+((1:nzv)-1)*dz;
   end
   %
-  if x_periodic==1
-    xw1=a+((1:nxw)-1/2)*dx;
+  if xperiodic==1
+    xw=a+((1:nxw)-1/2)*dx;
   else
-    xw1=a+((1:nxw)-1)*dx;
+    xw=a+((1:nxw)-1)*dx;
   end
-  if y_periodic==1
-    yw1=c+((1:nyw)-1/2)*dy;
+  if yperiodic==1
+    yw=c+((1:nyw)-1/2)*dy;
   else
-    yw1=c+((1:nyw)-1)*dy;
+    yw=c+((1:nyw)-1)*dy;
   end
-  if z_periodic==1
-    zw1=e+((1:nzw)-1)*dz;
+  if zperiodic==1
+    zw=e+((1:nzw)-1)*dz;
   else
-    zw1=e+((1:nzw)-3/2)*dz;
+    zw=e+((1:nzw)-3/2)*dz;
   end
   %
-  if x_periodic==1
-    xp1=a+((1:nxp)-1/2)*dx;
+  if xperiodic==1
+    xp=a+((1:nxp)-1/2)*dx;
   else
-    xp1=a+((1:nxp)-1)*dx;
+    xp=a+((1:nxp)-1)*dx;
   end
-  if y_periodic==1
-    yp1=c+((1:nyp)-1/2)*dy;
+  if yperiodic==1
+    yp=c+((1:nyp)-1/2)*dy;
   else
-    yp1=c+((1:nyp)-1)*dy;
+    yp=c+((1:nyp)-1)*dy;
   end
-  if z_periodic==1
-    zp1=e+((1:nzp)-1/2)*dz;
+  if zperiodic==1
+    zp=e+((1:nzp)-1/2)*dz;
   else
-    zp1=e+((1:nzp)-1)*dz;
+    zp=e+((1:nzp)-1)*dz;
   end
 end
\ No newline at end of file
diff --git a/matlab/generate_procgrid.m b/matlab/generate_procgrid.m
index 6c8dbce..08b07a6 100644
--- a/matlab/generate_procgrid.m
+++ b/matlab/generate_procgrid.m
@@ -4,9 +4,17 @@ function [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
           ibegp,iendp,jbegp,jendp,kbegp,kendp]=generate_procgrid(...
               nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
               nxp,nyp,nzp,nxprocs,nyprocs,nzprocs)
-  %
-  % determine the 3d decomposition (pointers) of fields for u,v,w,p
-  %
+  %[ibegu,iendu,jbegu,jendu,kbegu,kendu,...
+  % ibegv,iendv,jbegv,jendv,kbegv,kendv,...
+  % ibegw,iendw,jbegw,jendw,kbegw,kendw,...
+  % ibegp,iendp,jbegp,jendp,kbegp,kendp]=generate_procgrid(...
+  %     nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
+  %     nxp,nyp,nzp,nxprocs,nyprocs,nzprocs)
+  % Generates a processore grid for a given grid.
+  % Input
+  %   nxu,nyu,nzu,...     staggered grids
+  %   nxprocs,...         number of processors
+  
   chi='_';
   % U:
   [ibegu,iendu]=mpe_decomp1d(nxu,nxprocs,chi);
diff --git a/matlab/read_field_complete_ucfmf.m b/matlab/read_field_complete_ucfmf.m
deleted file mode 100644
index ee9191d..0000000
--- a/matlab/read_field_complete_ucfmf.m
+++ /dev/null
@@ -1,61 +0,0 @@
-function [q,x,y,z] = read_field_complete_ucfmf(fdir,iseq,field,varargin)
-    % [q,x,y,z] = read_field_complete_ucfmf(fdir,iseq,field,varargin)
-    % Reads a specific field from all processor chunks (multi-file tar version)
-    % Input
-    %   fdir                directory with files
-    %   iseq                sequence index
-    %   field               field to be read 
-    %                         {'u','v','w','p','s1','s2',...}
-    % ? step                index of step to be read (default: 1)
-    % ? verbosity           verbose output? (default: 0)
-    % ? debug               debug output? (default: 0)
-    % Output
-    %   q                   complete field
-    %   x,y,z               corresponding grid
-
-    % Parse optional input arguments
-    par = inputParser;
-    addParamValue(par,'step',1,@isnumeric);
-    addParamValue(par,'verbosity',0,@isnumeric);
-    addParamValue(par,'debug',0,@isnumeric);
-    parse(par,varargin{:});
-    istep = par.Results.step;
-
-    % Parse field
-    if ischar(field)
-      switch field(1)
-      case 'u'; ifield=1; fbase='uvwp'; cmesh='u';
-      case 'v'; ifield=2; fbase='uvwp'; cmesh='v';
-      case 'w'; ifield=3; fbase='uvwp'; cmesh='w';
-      case 'p'; ifield=4; fbase='uvwp'; cmesh='p';
-      case 's'; ifield=str2double(field(2:end)); fbase='scal'; cmesh='p';
-      otherwise; error('Invalid field: %s',field);
-    end
-    
-    % Read parameters
-    fparam = sprintf('%s/parameters_%04d.asc',fdir,iseq);
-    params = read_parameters_ucf(fparam,'tarmode',0);
-    
-    % Construct an array with final mesh size
-    nx = params.mesh.(sprintf('nx%s',cmesh));
-    ny = params.mesh.(sprintf('ny%s',cmesh));
-    nz = params.mesh.(sprintf('nz%s',cmesh));
-    q = zeros(nx,ny,nz);
-    
-    % Read grid
-    fgrid = sprintf('%s/grid_%04d.bin',fdir,iseq);
-    [x.u,y.u,z.u,x.v,y.v,z.v,x.w,y.w,z.w,x.p,y.p,z.p] = read_grid_ucf(fgrid,'tarmode',0);
-    x = x.(cmesh);
-    y = y.(cmesh);
-    z = z.(cmesh);
-    
-    % Get number of processors
-    nprocs = params.parallel.nprocs;
-    
-    % Now read chunk by chunk
-    for iproc=0:nprocs-1
-      fq = sprintf('%s/%s_%04d.%05d',fdir,fbase,iseq,iproc);
-      [data,ib,jb,kb,nxl,nyl,nzl] = read_field_chunk_ucf(fq,ifield,'tarmode',0,'ghost',0,'step',istep);
-      q(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1) = data;
-    end
-end