clear all

addpath('./mitgcm_tools/');

%directory where model outputs the GRID; This is something the users will get
%if they run the mitgcm configuration
dirgrid='./GRID/';
nx0=1080;ncut1=1260;ncut2=540;ny0=2*ncut1+nx0+ncut2;nz0=90;
nfx0=[nx0 0 nx0 ncut2 ncut1];nfy0=[ncut1 0 nx0 nx0 nx0];

%set the run directory where the mitgcm model configuration is run:
dirrun=['./'];if(~exist(dirrun));error('dirrun not exist');end;

varstr={'THETA'};ivar=1;

%========================================================
%================= EXTRACT WET    =======================
%========================================================
%step0: extract only wet points for storage; this requires reading in 
%the landmask to extract wet points., quite a large file, size 1.7GB
step_extract_wet=1;
if(step_extract_wet);	%/*{{{*/ [step_extract_wet]
  tmp=dir([dirgrid 'hFacC.data']);
  precIn='real*4';if(tmp.bytes/nx0/ny0/nz0/4==2);precIn='real*8';end;
  tic;hfc=readbin([dirgrid 'hFacC.data'],[nx0 ny0 nz0],1,precIn);toc;	%takes 13-20sec
  iwet=find(hfc(:)>0);							%111578238
  clear tmp precIn hfc	%clear memory

%now search raw data and extract only wet points:
  dirin =[dirrun 'diags/' varstr{ivar} '/'];    if(~exist(dirin));error('dir where diags are not exist');end;
  dirout=[dirrun 'diags/' varstr{ivar} '_wet/'];if(~exist(dirout));mkdir(dirout);end;

  flist=dir([dirin varstr{ivar} '.*.data']);if(length(flist)==0);error('diagnostics missing');end;
  precIn='real*4';if(flist(1).bytes/nx0/ny0/nz0/4==2);precIn='real*8';end;
  for ifile=1:length(flist);
    ain=zeros(nx0,ny0,nz0);
    tic;ain=readbin([flist(ifile).folder '/' flist(ifile).name],[nx0 ny0 nz0],1,precIn);toc;%takes 100sec per file
    clear a;a=ain(iwet);clear ain
    fout=[dirout flist(ifile).name];	%store just wet points
    tic;writebin(fout,a,1,'real*4');toc;
    fprintf('%i ',ifile);fprintf('%s\n',fout);
    clear a 
  end;
  clear iwet
  write_generic_meta(dirout,[Lwet 1 1],precIn)
end;		%/*}}}*/ [step_extract_wet]

%========================================================
%================= EXTRACT TILES  =======================
%========================================================

%step1: re-organizing to tiles and containing all time record;  This step can
%be achieved from the raw fully 3D field.  However, if for some reasons the raw 
%files were already being compacted to contain only wet points (step_extract_wet
%above), then this also allow for expanding back from wet to fully 3D.)
step_extract_tiles=1;
read_from_wet=1;	%if read from fully 3D raw files; else set this to 1
if(step_extract_tiles);	%/*{{{*/ [step_extract_tiles]
  
  varName=varstr{ivar};
  ix=1:nx0;iy=nfy0(1)+1:nfy0(1)+nx0;iz=1:83;
  nx=length(ix);
  ny=length(iy);
  nz=length(iz);
  dtilex=90;dtiley=90;
  domain_dims=[nx ny nz];tileSize=[dtilex dtiley];
  if(read_from_wet);
    dirin=[dirrun '/diags/' varName '_wet/'];
  else;
    dirin=[dirrun '/diags/' varName '/'];
  end;
  if(~exist(dirin));error('dirin not exist');end;
  flist=dir([dirin varName '.*.data']);if(length(flist)==0);error('data missing in dirin');end;
%now define dirout based on chosen tile size [dtilex dtiley]:
  dirout=[dirrun '/diags/' varName '_' num2str(dtilex) 'x' num2str(dtiley) '/'];
  if(~exist(dirout));mkdir(dirout);end;

  get_domain_tile(domain_dims,tileSize,ix,iy,iz,varName,dirin,dirout,dirgrid,read_from_wet)

end;	%/*}}}*/ [step_extract_tiles]
%--------------------------------------------------------

%--------------------------------------------------------
%now extract the tiles for the grid as well:
step_extract_GRID_tiles=1;
if(step_extract_GRID_tiles);	%/*{{{*/ [step_extract_GRID_tiles]
  dtilex=90;dtiley=90;
  ix=1:nx0;iy=nfy0(1)+1:nfy0(1)+nx0;iz=1:83;
  nx=length(ix);
  ny=length(iy);
  nz=length(iz);
  domain_dims=[nx ny nz];tileSize=[dtilex dtiley];
  dirgridout=[dirrun 'GRID_' num2str(dtilex) 'x' num2str(dtiley) '/'];
  if(~exist(dirgridout));mkdir(dirgridout);end;
  get_domain_tile_GRID(domain_dims,tileSize,ix,iy,iz,dirgrid,dirgridout);
end;				%/*}}}*/
%=================================================

%========================================================
%================= NETCDF TILES  ========================
%========================================================

%This is to write the binary tiles into netcdf files for storage instead of binary,
step_write_netcdf=1;
if(step_write_netcdf);
  dtilex=90;dtiley=90;
  dirin = [dirrun '/diags/' varstr{ivar} '_' num2str(dtilex) 'x' num2str(dtiley) '/'];
  dirgridin=[dirrun 'GRID_' num2str(dtilex) 'x' num2str(dtiley) '/'];
  if(~exist(dirin));error('dirin not exist');end;
  if(~exist(dirgridin));error('dirgridin not exist');end;
  filAvailDiag=[dirrun '/available_diagnostics.log'];
  if(~exist(filAvailDiag));error('filAvailDiag missing');end;
  deltaT=240;	%<-- info from namelists/data
  yrstart=2002; %<-- info from namelists/data.cal
%Make sure to update the header inside the script below as it describes the run's details
%If there's any citation, it should go in there.
  convert_field2nctiles(varstr{ivar},dirin,dirgridin,yrstart,deltaT,filAvailDiag);
end;

step_write_grid_netcdf=1;
if(step_write_grid_netcdf);
  dtilex=90;dtiley=90;
  iz=1:83;
  inDir=[dirrun '/GRID_' num2str(dtilex) 'x' num2str(dtiley) '/'];
    if(~exist(inDir));error('missing input grid');end;
  outDir=[dirrun 'GRID_' num2str(dtilex) 'x' num2str(dtiley) '_nc/'];
    if(~exist(outDir));mkdir(outDir);end;
  tile_dims=[dtilex dtiley length(iz)];
  convert_field2nctiles_GRID(tile_dims,inDir,outDir);
end;
