function [ii,jj,kk,xq,yq,zq] = find_region_ucfmf(fdir,iseq,field,xb,yb,zb,varargin)
    % [ii,jj,kk,xq,yq,zq] = find_region_ucfmf(fdir,iseq,field,xb,yb,zb,varargin)
    % Determines the global index of grid points which lie within the bounds
    % xb,yb,zb specified in the parameters, i.e. the grid points for which
    %   xb(1) <= x <= xb(2), 
    %   yb(1) <= y <= yb(2),
    %   zb(1) <= z <= zb(2).
    % If the domain is periodic, the specified region can lie "outside" of the domain,
    % i.e. it is automatically translated to coordinates in the inner domain.
    % If xb(1)==xb(2), a slice located at the nearest neighbor point will be created.
    % Input
    %   fdir                directory with files
    %   iseq                sequence index
    %   field               field to be indexed 
    %                         {'u','v','w','p','s1','s2',...}
    %   xb,yb,zb            arrays which describe interval, e.g. [0,1] for all points between 0 and 1
    % ? fuzzy               Include the first point outside of the region? (default: 0)
    % ? verbosity           verbose output? (default: 0)
    % ? debug               debug output? (default: 0)
    % Output
    %   ii,jj,kk            array with field indices
    %   xq,yq,zq            corresponding grid, monotonically increasing, even if periodic boundary is crossed
  
    % Parse input
    par = inputParser;
    addParamValue(par,'verbosity',0,@isnumeric);
    addParamValue(par,'debug',0,@isnumeric);
    addParamValue(par,'fuzzy',0,@isnumeric);
    parse(par,varargin{:});
    iverb = par.Results.verbosity;
    idebug= par.Results.debug;
    ifuzzy= par.Results.fuzzy;
    
    % Parse bounds
    xmin = xb(1);
    xmax = xb(2);
    ymin = yb(1);
    ymax = yb(2);
    zmin = zb(1);
    zmax = zb(2);
    
    % Parse 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 auxillary files
    fparam = sprintf('%s/parameters_%04d.asc',fdir,iseq);
    fgrid  = sprintf('%s/grid_%04d.bin',fdir,iseq);
    fproc  = sprintf('%s/proc_%04d.bin',fdir,iseq);
    params = read_parameters_ucf(fparam,'tarmode',0,'verbosity',iverb);
    [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,'verbosity',iverb,'debug',idebug);
    x = x.(cmesh);
    y = y.(cmesh);
    z = z.(cmesh);
    [ibeg.u,iend.u,jbeg.u,jend.u,kbeg.u,kend.u,...
     ibeg.v,iend.v,jbeg.v,jend.v,kbeg.v,kend.v,...
     ibeg.w,iend.w,jbeg.w,jend.w,kbeg.w,kend.w,...
     ibeg.p,iend.p,jbeg.p,jend.p,kbeg.p,kend.p] = read_procgrid_ucf(fproc,'tarmode',0,'verbosity',iverb,'debug',idebug);
    ibeg = ibeg.(cmesh);
    jbeg = jbeg.(cmesh);
    kbeg = kbeg.(cmesh);
    iend = iend.(cmesh);
    jend = jend.(cmesh);
    kend = kend.(cmesh);
    nx = numel(x);
    ny = numel(y);
    nz = numel(z);
    
    % Parse domain information
    a = params.geometry.a;
    b = params.geometry.b;
    c = params.geometry.c;
    d = params.geometry.d;
    e = params.geometry.e;
    f = params.geometry.f;
    x_periodic = params.geometry.xperiodic;
    y_periodic = params.geometry.yperiodic;
    z_periodic = params.geometry.zperiodic;
    
    % Check if bounds are valid
    if xmax<xmin
      error('xmax must be greater than xmin.');
    end
    if ymax<ymin
      error('ymax must be greater than ymin.');
    end
    if zmax<zmin
      error('zmax must be greater than zmin.');
    end
    
    % Check if min and max are the same point (if so, treat it as slice)
    feps = 1e-12;
    islice = false;
    jslice = false;
    kslice = false;
    if abs(xmax-xmin)<feps
      islice = true;
    end
    if abs(ymax-ymin)<feps
      jslice = true;
    end
    if abs(zmax-zmin)<feps
      kslice = true;
    end
    
    % Check values outside of domain, circular shift if periodic, otherwise error
    % Keep the number of shifts and initial offset to reconstruct grid
    if x_periodic 
      nxloop   = floor((xmax-xmin)/(b-a));
      nxoffset = floor(xmin/(b-a));
      xmin = mod(xmin-a,b-a)+a;
      xmax = mod(xmax-a,b-a)+a;
    elseif xmin<a || xmin>b || xmax<a || xmax>b
      error('Interval must be located inside domain.');
    else
      nxloop   = 0;
      nxoffset = 0;
    end
    if y_periodic 
      nyloop   = floor((ymax-ymin)/(d-c));
      nyoffset = floor(ymin/(d-c));
      ymin = mod(ymin-c,d-c)+c;
      ymax = mod(ymax-c,d-c)+c;
    elseif ymin<c || ymin>d || ymax<c || ymax>d
      error('Interval must be located inside domain.');
    else
      nyloop   = 0;
      nyoffset = 0;
    end
    if z_periodic 
      nzloop   = floor((zmax-zmin)/(f-e));
      nzoffset = floor(zmin/(f-e));
      zmin = mod(zmin-e,f-e)+e;
      zmax = mod(zmax-e,f-e)+e;
    elseif zmin<e || zmin>f || zmax<e || zmax>f
      error('Interval must be located inside domain.');
    else
      nzloop   = 0;
      nzoffset = 0;
    end
    
    % Find gridpoints closest to the bounds
    if ifuzzy
      [~,imin] = find((x-xmin)<0,1,'last');
      [~,imax] = find((x-xmax)>0,1,'first');
      [~,jmin] = find((y-ymin)<0,1,'last');
      [~,jmax] = find((y-ymax)>0,1,'first');
      [~,kmin] = find((z-zmin)<0,1,'last');
      [~,kmax] = find((z-zmax)>0,1,'first');
    else
      [~,imin] = find((x-xmin)>=0,1,'first');
      [~,imax] = find((x-xmax)<=0,1,'last');
      [~,jmin] = find((y-ymin)>=0,1,'first');
      [~,jmax] = find((y-ymax)<=0,1,'last');
      [~,kmin] = find((z-zmin)>=0,1,'first');
      [~,kmax] = find((z-zmax)<=0,1,'last');
    end
    if isempty(imax); imax=1; end
    if isempty(jmax); jmax=1; end
    if isempty(kmax); kmax=1; end
    
    % If a slice is requested, correct previous results
    if islice
      [~,imin] = min(abs(x-xmin));
      imax = imin;
    end
    if jslice
      [~,jmin] = min(abs(y-ymin));
      jmax = jmin;
    end
    if kslice
      [~,kmin] = min(abs(z-zmin));
      kmax = kmin;
    end
    
    % Construct the points which are inside of the interval and reconstruct the grid
    if imin<=imax
      ii = [repmat(mod(imin+(0:nx-1)-1,nx)+1,[1,nxloop]),imin:imax];
      shift = [reshape((repmat([zeros(1,nx-imin+1),ones(1,imin-1)],nxloop,1)+(1:nxloop)')',1,[]),...
               ones(1,imax-imin+1)+nxloop]-1;
      xq = x(ii)+(shift+nxoffset)*(b-a);
    else
      ii = [imin:nx,repmat(1:nx,[1,nxloop]),1:imax];
      shift = [ones(1,nx-imin+1),...
               reshape((repmat(ones(1,nx),nxloop,1)+(1:nxloop)')',1,[]),...
               ones(1,imax)+nxloop+1]-1;
      xq = x(ii)+(shift+nxoffset)*(b-a);
    end
    if jmin<=jmax
      jj = [repmat(mod(jmin+(0:ny-1)-1,ny)+1,[1,nyloop]),jmin:jmax];
      shift = [reshape((repmat([zeros(1,ny-jmin+1),ones(1,jmin-1)],nyloop,1)+(1:nyloop)')',1,[]),...
               ones(1,jmax-jmin+1)+nyloop]-1;
      yq = y(jj)+(shift+nyoffset)*(d-c);
    else
      jj = [jmin:ny,repmat(1:ny,[1,nyloop]),1:jmax];
      shift = [ones(1,ny-jmin+1),...
               reshape((repmat(ones(1,ny),nyloop,1)+(1:nyloop)')',1,[]),...
               ones(1,jmax)+nyloop+1]-1;
      yq = y(jj)+(shift+nyoffset)*(d-c);
    end
    if kmin<=kmax
      kk = [repmat(mod(kmin+(0:nz-1)-1,nz)+1,[1,nzloop]),kmin:kmax];
      shift = [reshape((repmat([zeros(1,nz-kmin+1),ones(1,kmin-1)],nzloop,1)+(1:nzloop)')',1,[]),...
               ones(1,kmax-kmin+1)+nzloop]-1;
      zq = z(kk)+(shift+nzoffset)*(f-e);
    else
      kk = [kmin:nz,repmat(1:nz,[1,nzloop]),1:kmax];
      shift = [ones(1,nz-kmin+1),...
               reshape((repmat(ones(1,nz),nzloop,1)+(1:nzloop)')',1,[]),...
               ones(1,kmax)+nzloop+1]-1;
      zq = z(kk)+(shift+nzoffset)*(f-e);
    end
end