function get_domain_tile(domain_dims,tileSize,ix,iy,iz,varName,dirin,dirout,dirgrid,read_wet);
%--------------------------------------
%function get_domain_tile(domain_dims,tileSize,ix,iy,iz,varName,dirin,dirout,dirgrid,read_wet);
%
%Read the full domain model 3D outputs, divide into tiles and produce binary files, one per tile
%
%Input:
%       domain_dims     [nx ny nz]
%       tileSize        [ntilex ntiley]
%       ix,iy,iz        indices, e.g., for the Arctic region within a larger domain
%       varName         variable name, e.g., 'THETA', a string 
%       dirin           Directory (full path) where the model original or wet diags lives 
%                       This is either, e.g., 'diags/THETA/' or 'diags/THETA_wet/';
%       dirout          Directory (full path) where tiled output will be written 
%                       (a directory, e.g., /.../THETA_90x90/'] will be created
%       dirgrid         Directory (full path) where the model original grid lives
%       read_wet        1 for reading wet points only, 0 for reading original raw 3D files
%Output:
%       None.  But files will be written out to dirout, one file per tile, named e.g., THETA.0001.data
%
% written by Victor Ocana 2019
% modified by An Nguyen Jul 2022
%
%Example of values of inputs:
%%input:
%  varstr={'THETA'};ivar=1;
%  varName=varstr{ivar};
%  ix=1:nx0;iy=nfy0(1)+1:nfy0(1)+nx0;iz=1:83;
%  dtilex=90;dtiley=90;
%  nx=length(ix);ny=length(iy);nz=length(iz);
%  domain_dims=[nx ny nz];tileSize=[dtilex dtiley];
%%now define dirout based on chosen tile size [dtilex dtiley]:
%  dirgrid=[dirrun 'GRID/'];
%  if(read_wet);
%    dirin=[dirrun '/diags/' varName '_wet/'];
%  else;
%    dirin =[dirrun '/diags/' varName '/'];
%  end;
%  dirout=[dirrun '/diags/' varName '_' num2str(dtilex) 'x' num2str(dtiley) '/'];
%  read_wet=1;
%--------------------------------------
  
%begin function below:

%check dirin and diagnostics existence:
%=====================================================================
if(~exist(dirin));error('dirin where diags are not exist');end;
flist=dir([dirin varName '.*.data']);L=length(flist);if(L==0);error('diagnostic files missing');end;

%Directories and dimensions
%=====================================================================
  nx=domain_dims(1);
  nz=1; if(length(domain_dims)==3); nz=domain_dims(3);end;
  ny=domain_dims(2);

  ntilex=tileSize(1);
  ntiley=tileSize(2);
  ntile=nx*ny/ntilex/ntiley;

%Set output directory   %------------------------------------------------------------------
if(0);
if(~exist(dirout));
  mkdir(dirout);
else;
   disp(['The directory ' dirout ' already exists. If Continue = Y or Return Key, it will be overwritten']);
   prompt = 'Continue? Y/N [Y]: ';
   str = input(prompt,'s');
   if isempty(str);
      str = 'Y';
   end;
   if(~strcmp(str,'Y')); error('Script interrupted');end;
end;
end;

%Get list of tiles with wet points
%=====================================================================
meta=read_meta('hFacC.meta',dirgrid);
%original size of raw data:
nx0=domain_dims(1);ny0=meta.dimList(1)*meta.dimList(4)/nx0;nz0=meta.dimList(7);
precIn='real*4';if(str2num(meta.dataprec(end-1:end))==64);precIn='real*8';end;
mskC=readbin([dirgrid 'hFacC.data'],[nx0 ny0],1,precIn);
mskC=mskC(ix,iy);mskC(find(mskC(:)~=0))=1;
 
%if read_wet, we need to also load in landmask to find wet points:
%first read in the landmask, this takes a LONG TIME
if(read_wet);
  tic;hfc=readbin([dirgrid 'hFacC.data'],[nx0 ny0 nz0],1,precIn);toc;	%~10-30sec
  ind=find(hfc(:)>0);clear hfc;			%111578238
else;
  ind=1:nx*ny*nz;
end;
Lmax=length(ind);

%define tile size , 90 x 90 is good; smaller tiles = too many files; larger = too large of file size
flags=[];
flags.name = varName;
flags.dDir = dirin;
flags.freq = 2635200;%monthly; 3600; %hourly %freq; %<-- choose based on data.diagnostics
flags.prec = 'R';
flags.vert = ' ';
flags.cumu = ' ';

%read_meta; note if read from wet, the meta is written out to a "2d" meta file:
if(read_wet);
  fmeta=['generic_2d_real' precIn(end) '.meta'];
else;
  fmeta=[flist(1).name(1:end-5)];
end;
clear meta;[meta]  =read_meta(fmeta,dirin);
if(read_wet);meta.fldList={varName};end;

%Set precision for read/write binaries
precIn='real*4';if(strcmp(meta.dataprec,'float64'));precIn='real*8';end;

%Set vertical dimension
if(meta.nDims==2); nz=1; elseif(meta.nDims==3); nz=meta.dimList(7); end;
if(meta.nDims==3 & nz>iz(end));nz=iz(end);end;
if(read_wet); nz=iz(end);end;

%Summary of the parameters for this routine
fprintf('Data file: %s\n',varName);
fprintf('Reading data from: %s\n',dirin);
fprintf('Writing tiles to : %s\n',dirout);
fprintf('Tile size: %2d x %2d x %2d\n',ntilex,ntiley,nz);
fprintf('Data precission %s\n',precIn);
fprintf('Number of files (timesteps) to be read: %i\n',L);

%Get list of tiles with wet points and indices for each tile
%=====================================================================
%first get the mask of land/ocean:

  [tiles] = domain2tiles_v2([],domain_dims,tileSize,mskC);
        %      tot: 144			%      nfx: [1080 1080]
        %      wet: 129                 %      nfy: [1080 1080]
        %     list: [144x2 double]      %     mask: [1080x1080 double]
        %    tsize: [90 90]             %    field: []
        %    fsize: [1080 1080 1]       %    index: {1x129 cell}
        
  ntile_wet = tiles.wet;
  fprintf('Number of wet tiles: %2d\n',ntile_wet);fprintf('\n');

  write_generic_meta(dirout,[tileSize nz],precIn);	%write a file generic_3d_real4.meta

  mout = [dirout strcat(meta.fldList{1}) '_' num2str(ntilex) 'x' num2str(ntiley) '.meta'];
  write_meta(mout,[tileSize nz],precIn,{meta.fldList{1}},L); %write a file THETA_90x90.meta 

%=================================================
  %Open files to write tiles
  %%% **** IMPORTANT **** we're opening ALL tile files at once! LOTS OF FILES ARE BEING OPENED
  %%% **** make sure to run fclose all after this portion of code is done ********************
  fclose all;
  for jtile=1:ntile_wet;
     fname{jtile}=[dirout strcat(meta.fldList{1}) '.' sprintf('%4.4i',jtile) '.data'];
     %<-- critical line: opening ntile_wet files at once -->%
     eval(['fid' num2str(jtile) '=fopen(fname{jtile},''w'',''b'');']);
  end;

%--------------------------------------------
%now reading diags:
%looping through all available model output files, either raw 3D or wet:
  idot=find(flist(1).name=='.');idot=idot(1)+1:idot(2)-1;
  tt=zeros(L,1);
  for itime=1:L
    t0=clock;
    ts=flist(itime).name(idot);
    tt(itime)=str2num(flist(itime).name(idot));

%actual reading of field, can take quite long depending on reading raw 3D (~90sec) or wet (~40sec):
    tic;ain=readbin([flist(itime).folder '/' flist(itime).name],[1 Lmax],1,precIn);toc;
    if(read_wet);
      tic;a=zeros(nx,ny0,nz0);toc;                                   %0.02sec
      tic;a(ind)=ain;toc;                                            %6.2sec
    else;
      a=ain;
    end;
    clear ain
    a=reshape(a,nx0,ny0,nz0);
%now extract just the Arctic domain:
    a=a(ix,iy,iz);
%now split entire Arctic domain into tiles:
    [field_tiles] = field2tiles_v2(a,tiles.index,tileSize);
    for jtile=1:ntile_wet;
       eval(['fid=fid' num2str(jtile) ';']);
       fwrite(fid,field_tiles{jtile},meta.dataprec);
    end;
    fprintf('%i %f ',[itime,etime(clock,t0)]);
  end;
  fprintf('\n');
  clear a
%--------------------------------------------

  %close tile files
  for jtile=1:ntile_wet;
     eval(['fclose(fid' num2str(jtile) ');']);
  end;
  fclose all;
  ftime=[dirout 'time_list_' strcat(meta.fldList{1}) '.txt'];

  fid=fopen(ftime,'w');
  for it=1:L;fprintf(fid,'%10.0i\n',tt(it));end;
  fclose(fid);
  fprintf('Wrote %s\n',ftime);
  save([dirout 'flags.mat'],'flags');

  fprintf('\n');
  fclose all;

return