function [tiles] = domain2tiles_v2(field,domain_dims,tileSize,maskIn);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%function [tiles] = domain2tiles_v2(field,domain_dims,tileSize,maskIn);
%
%Given a field in binary form, return a structure with the data in a cell
%for each tile, a mask with tile numbers, a structure with the indices
%from the binary form for each tile, and a list of tiles and wet tiles.  
%If field is empty, the list of tiles and the indices are calculated
%
%Input:
%	field		a 2D or 3D data field in binary format 
%                       (e.g. [1080 1080] or [1080 1080 83])
%	domain_dims	[nx ny nz]
%	tileSize	[dtilex dtiley] (number of grid points in each tile side)
%	maskIn	        landmask to determine where wet/dry points are
%
%Output:
%	tiles		structure with fields:
%		tot		total number of tiles in the selected domain
%		wet		number of wet tiles
%		list		vector with two columns: First column numbers all tiles in the regional 
%                               domain; second column has either 0 for dry tiles or the  wet-tile number  
%		tsize		tileSize [dtilex dtiley] (number of grid points in each tile side)
%		fsize		size of the regional compact form [nx ny nz]
%		nfx		x-size of the regional domain, can be multi-dim if domain contains many "faces"
%		nfy		y-size of the 1 regional domain, can be multi-dim if domain contains many "faces"
%		mask		domain faces cell structure with the value "tilenumber" at each location
%		index		cell structure with the indices, form the compact format, for each tile
%		field		cell structure with data in tiles, field_tiles{tilenumber} has dimensions [dtilex dtiley nz]
%
% Written by Victor Ocana 2019
% Modified by An Nguyen Jul 2022
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tiles=[];
%Dimensions	%------------------------------------------------------------
sf=size(field);

nx     = domain_dims(1);
ny     = domain_dims(2);
nz     = 1;   if(length(sf)>2); nz = sf(3); end;

nfx=[nx nx];
nfy=[nx nx];

dtilex = tileSize(1);
dtiley = tileSize(2);
ntile  = nx*ny/dtilex/dtiley;

%Check dimensions 	%------------------------------------------------------------
if (~isempty(field));
   if (sf(1)~=nx || sf(2)~=ny);
      str1=['DOMAIN2tiles_v2 - ERROR!!!: Dimensions of the data field [' num2str(sf) '] '];
      str2=['are different from domain_dims [' num2str(nx) ' ' num2str(ny) ' ' num2str(nz) ']'];
      disp([str1,str2]);
      return;
   end;
else;
   disp(['DOMAIN2tiles_v2 - WARNING!!!: field is empty. Calculating tile numbers and indices']); 
end;

if(length(domain_dims)>3); 
   if(nz~=domain_dims(4));
      str1=['DOMAIN2tiles_v2 - WARNING: vertical dimension from data field ' num2str(nz)];
      str2=[' is different from vertical dimension from domain_dims ' num2str(domain_dims(4)) ];
      disp([str1 str2]);
      disp(['           Continue with vertical dimension from data field nz=' num2str(nz)]);
   end;
end;

%Create a mask with indices = tile number;	%--------------------------------------------------
%and a list of total and wet tiles
list_tiles=zeros(ntile,2);
mski=zeros(nx,ny);

tot_tiles=0;
wet_tiles=0;
sz=size(maskIn);
nxf=sz(1)/dtilex;
nyf=sz(2)/dtiley;
for j=1:nyf
  iy=(j-1)*dtiley+1:j*dtiley;
  for i=1:nxf
    ix=(i-1)*dtilex+1:i*dtilex;
    tot_tiles=tot_tiles+1;
    tmp=maskIn(ix,iy,:);
    list_tiles(tot_tiles,1)=tot_tiles;
    if(sum(tmp(:))>0);
      wet_tiles=wet_tiles+1;
      list_tiles(tot_tiles,2)=wet_tiles;
      mski(ix,iy,:)=wet_tiles;
    end;
  end%i
end;%j

mask = mski;

%Store indices and field structures, one cell for each tile
indices_tiles = [];
field_tiles = [];
for ii=1:wet_tiles;
   tmp = find(mski==ii);
   indices_tiles{ii}=tmp;
   if(~isempty(field));
      for iz=1:nz;
         ftmp = field(:,:,iz);
         ftmp1= ftmp(tmp); ftmp1=reshape(ftmp1,dtilex,dtiley);
         ftmp2(:,:,iz) = ftmp1;
      end;
      field_tiles{ii} = ftmp2;
      clear ftmp ftmp1 ftmp2;
   end;
   clear tmp;
end;

tiles.tot=tot_tiles;
tiles.wet=wet_tiles;
tiles.list=list_tiles;
tiles.tsize = tileSize;
tiles.fsize = [nx ny nz];
tiles.nfx = nfx;
tiles.nfy = nfy;
tiles.mask = mask;
tiles.field = field_tiles;
tiles.index = indices_tiles;

return