function []=convert_field2nctiles(varName,dataDir,gridDir,startyr,deltaT,filAvailDiag);
%addpath('/home1/03901/atnguyen/matlab/MITprof_toolbox/');
%-----------
%function []=convert_field2nctiles(varName,dataDir,gridDir,startyr,deltaT,filAvailDiag);
% Read binary files already in tile form and write nc output in ECCO form
%
% Inputs:
%	varName		the diagnostic being writen
%	dataDir   	Contains subdirectories, one for each diagnostic, with .data and generic**.meta files
%			Each subdirectory has one file for each tile, .data binary format
%			A generic**.meta file and time_list_varName file should be in the directory
%	gridDir		the grid in binary tile form  
% 	startyr		first year of simulation
%	deltaT		Time step of the simulation
%	filAvailDiag 	full path of the file available_diagnostics.log
%
%
% NC files have:
%	5 dimensions: 
%		0 itxt  30
%		1 i1  number of time steps 
%		2 i2  vertical levels
%		3 i3  y length
%		4 i4  x length
%
%	13 Variables (dims):
%		0 i1  		i1
%		1 i2  		i2
%		2 i3  		i3
%		3 i4  		i4
%		4 varName  	i4 i3 i2 i1
%		5 lon  		i4 i3
%		6 lat  		i4 i3
%		7 dep  		i2
%		8 tim  		i1
%		9 timstep  	i1
%		10 land  	i4 i3 i2
%		11 area  	i4 i3
%		12 thic  	i2
%		
% Written by Victor Ocana 2019, following ECCO standard distributed gmcfaces/ code 
% Modified by An Nguyen Jul 2022
%-----------

%Directories, files		+---------------------------------------------------------- 
inDir=[dataDir];if(~exist(inDir));error('inDir not exist');end;
flist = dir([inDir varName '*.data']);if(length(flist)==0);error('tile data in inDir missing');end;

tmp=size(dataDir);if(dataDir(end)=='/');nend=tmp(2)-1;else;nend=tmp(2);end;
outDir = [dataDir(1:nend) '_nc/'];clear tmp
if(exist(outDir)==0);mkdir(outDir);end;

%Get information about variable +---------------------------------------------------------- 
tmp=dir([inDir '*.meta']); 
if(length(tmp)>1);	%if the generic file is also in
  j=[];
  for i=1:length(tmp);
    if(length(strfind(tmp(i).name,'generic'))==0);
      j=i;
    end;
  end;
else;
  j=1;
end;

tmp=tmp(j);
metaT = read_meta(tmp.name,inDir);
clear tmp;
if(isempty(metaT)); 
   disp(['convert_field2nctiles ERROR: no .meta file found']);
   return;
end;
if(~strcmp(varName,strcat(metaT.fldList{1}))); 
   warning(...
    'convert_field2nctiles: The given varName %s is different from the name in .meta file %s',...
                 varName,strcat(metaT.fldList{1}));
   varName=strcat(metaT.fldList{1});
end;

load([inDir '/flags.mat']);
%mannually fix available_diagnostics to add the "c" if [u,v]velmass
[avail_diag]=read_avail_diag_v2(filAvailDiag,strcat(metaT.fldList{1}));%v2: 9char instead of 8char

%since i changed nr , double check here, but only works if it's 3d field:
if(avail_diag.nr>1);
  if(avail_diag.nr ~=metaT.dimList(7));avail_diag.nr=metaT.dimList(7);end;
end;

nprec=str2num(metaT.dataprec(end-1:end))/8;
if(~(nprec==4|nprec==8));disp(['WARNING!: precission=' num2str(nprec) ' It should be either 4 or 8']);end;
prec = ['real*' num2str(nprec) ''];

%Set tile and time dimensions	+---------------------------------------------------------- 
dtilex = metaT.dimList(1);
dtiley = metaT.dimList(4);
nz =1 ; if(metaT.nDims>2); nz = metaT.dimList(7);end;

%Create time variables		+---------------------------------------------------------- 
timFile = ['time_list_' varName '.txt'];
if(~exist([inDir timFile]));error('timFile missing, need to create it in step_extract_tile');end;
timstep = textread([inDir timFile],'%d');
fid=fopen([inDir timFile]); 
nlines=0;while ~feof(fid); tmp=fgetl(fid); nlines=nlines+1;end;
fclose(fid);
if length(timstep)~=nlines;
   clear timstep;
   tmp = textread([inDir timFile],'%d'); tmp = reshape(tmp,2,length(tmp)/2)';
   timstep=tmp(:,1);
   clear tmp;
end;
dd = ts2dte(timstep,deltaT,startyr);
tim=datenum(dd)-datenum([startyr 1 0]); 
tim = tim-17;	%This centers the number of days since initial date in the middle of the month
clear dd;
timUnits=['days since ' num2str(startyr) '-1-1 0:0:0'];

%Set domain dimensions	+---------------------------------------------------------- 
ntiles = size(flist,1);
ntimes = size(timstep,1);

%Set data Attributes for netcdf +---------------------------------------------------------- 
longName = avail_diag.longNameDiag;
units    = avail_diag.units;
if avail_diag.nr~=1;
  coord='lon lat dep tim';
else;
  coord='lon lat tim';
end;
if(nprec==4); xtype ='NC_FLOAT';elseif(nprec==8);xtype ='NC_DOUBLE';end;  

%Set global attributes attributes       +---------------------------------------------------------- 
gatt_names = {'Description','Institution','References','Conventions','Software',...
              'netCDF Version','matlab Version','_FillValue','missing_value'};
gatt_value{1} = {'High resolution 3.5km run'};
gatt_value{2} = {'Oden Institute for Computational Engineering and Science, UT-Austin'};
gatt_value{3} = {'  ',...
        '   Nguyen, A.T., 2022.'};
gatt_value{4} = {'CF-1.6 (http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html)'};
gatt_value{5} = {'File created using convert_FIELD2nctiles.m'};
gatt_value{6} = {netcdf.inqLibVers};
gatt_value{7} = {version};
gatt_value{8} = {NaN};
gatt_value{9} = {NaN};

if(~(size(gatt_value)==size(gatt_names))); 
  disp('Number of global attribute names and values MUST be the same');
end;

%---------------------------------------------------------------------------------------------------------------------
%keyboard
%Loop over tiles 		+-------------------------------------------------------------------------------------
   fprintf('%s %5d %5d %5d %5d \n',varName,dtilex,dtiley,nz,ntimes);
   fprintf('Processing tile ');
for itile = 1:ntiles;
%for itile = 1;
   fprintf(' %d',itile);
   idot     = find(flist(itile).name=='.');
   tileName = flist(itile).name(idot(1)+1:idot(2)-1);
   fileTile = [outDir flist(itile).name(1:idot(end)) 'nc']; 

   %Read binary			+---------------------------------------------------------- 
   data = readbin([inDir flist(itile).name],[dtilex dtiley nz ntimes],1,prec); 
   data=squeeze(data);
   dims    = size(data);
%!!!!!!!!!!  	TO COMPLY WITH netCDF conventions CF-1.6 REVERSE ORDER OF DIMENSIONS WRT MITgcm  !!!!!!!!!!!!!!
   data    = reshape(data,fliplr(dims));
   nDim    = length(size(data));			%Number of dimensions
   dimsize = size(data);				%Size of each dimension
   
   %Create netcdf file		+---------------------------------------------------------- 
   %(by default using basic netcdf, unless data file is very large) 
   if prod(dimsize)*4/1e9<1.5;	
       cmode='clobber';	%Overwrite any existing file with the same name.
   else;		
       cmode='NETCDF4';	%To allow for large file: create a NetCDF-4/HDF5 file
   end;
   ncid = netcdf.create(fileTile,cmode);

   %Global Attributes	+------------------------------------------------------------------
   nc_global=netcdf.getConstant('NC_GLOBAL');
   
   for in=1:size(gatt_names,2);
      for iv=1:size(gatt_value{in},2);
         tmp=gatt_names{in}; 
         if (size(gatt_value{in},2)>1 & iv>1); tmp=sprintf('%2.2i',(iv-1));end;
         netcdf.putAtt(ncid,nc_global,tmp,gatt_value{in}{iv});
      end;
   end;
    
   %Dimensions (and dims as variables) +---------------------------------------------------
   %Create the dimension-variables, their names and their values
   for iDim=1:nDim;		
      if dimsize(iDim)~=1;
         dimlist{iDim}=['i' num2str(iDim)];
         dimName{iDim}=['array index ' num2str(iDim)];
         eval(['dimvec.i' num2str(iDim) '=[1:dimsize(iDim)];']);
      end;
   end;
   
   %Define dims
   netcdf.defDim(ncid,'itxt',30);   
   for dd=1:length(dimlist); 
      dimid(dd)= netcdf.defDim(ncid,dimlist{dd},dimsize(dd)); 
   end;

   %Define dimensions as variables and their attributes
   for dd=1:length(dimlist); 	
      varid = netcdf.defVar(ncid,dimlist{dd},'NC_DOUBLE',dimid(dd));
      vvarid(dd)=varid;
      netcdf.putAtt(ncid,varid,'long_name',dimName{dd});
      netcdf.putAtt(ncid,varid,'units','1');
   end;
   netcdf.endDef(ncid);

   %Fill the values of the dimension variables
   for dd=1:length(dimlist);	
      varid = netcdf.inqVarID(ncid,dimlist{dd});
      netcdf.putVar(ncid,varid,getfield(dimvec,dimlist{dd}));
   end;

   %Variables  +------------------------------------------------------------------------------

   %Data       +------------------------------------------------------------------------------
   %Define diagnostic variable and its attributes
   netcdf.reDef(ncid);
   %Notice the order of the dimensions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   varid = netcdf.defVar(ncid,varName,xtype,flipdim(dimid,2));
   if ~isempty(longName); netcdf.putAtt(ncid,varid,'long_name',longName); end;
   if ~isempty(units); netcdf.putAtt(ncid,varid,'units',units); end;
   if ~isempty(coord); netcdf.putAtt(ncid,varid,'coordinates',coord); end;

   %Write data
   netcdf.endDef(ncid);
   netcdf.putVar(ncid,varid,data); 

   %GRID	+------------------------------------------------------------------------------
   gxtype ='NC_DOUBLE';   
   [grid_diag]=set_grid_diag_v2(avail_diag,tileName,gridDir);
   
   %gvar = fieldnames(grid_diag)'; 
    %%%%gvar={grid_diag.lon grid_diag.lat grid_diag.dep tim timstep grid_diag.msk grid_diag.RAC grid_diag.dz};
   if (avail_diag.nr==1); %2D diags
      gvar      = {grid_diag.lon grid_diag.lat tim timstep grid_diag.msk grid_diag.RAC};
      gname     = {'lon','lat','tim','timstep','land','area'}; 
      glongName = {'longitude','latitude','time','final time step number','land mask','grid cell area'};
      gstdName  = {'longitude','latitude','time','final time step number','land_binary_mask','cell_area'};
      gunits    = {'degrees_east','degrees_north',timUnits,'1','1','m^2'};
      gdimid    = {dimid(2:3) dimid(2:3) dimid(1) dimid(1) dimid(2:3) dimid(2:3)};
   elseif(avail_diag.nr>1);  %3D diags
      gvar      = {grid_diag.lon grid_diag.lat grid_diag.dep tim timstep grid_diag.msk ...
                   grid_diag.RAC grid_diag.dz};
      gname     = {'lon','lat','dep','tim','timstep','land','area','thic'}; 
      glongName = {'longitude','latitude','depth','time','final time step number',...
                   'land mask','grid cell area','grid cell thickness'};
      gstdName  = {'longitude','latitude','depth','time','final time step number',...
                   'land_binary_mask','cell_area','cell_thickness'};
      gunits    = {'degrees_east','degrees_north','m',timUnits,'1','1','m^2','m'};
      gdimid    = {dimid(3:4) dimid(3:4) dimid(2) dimid(1) dimid(1) dimid(2:4) dimid(3:4) dimid(2)};
   end;

   netcdf.reDef(ncid);

   for igrid=1:length(gname);
%Notice the order of the dimensions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      gvarid(igrid) = netcdf.defVar(ncid,gname{igrid},gxtype,flipdim(gdimid{igrid},2)); 
      netcdf.putAtt(ncid,gvarid(igrid),'long_name',glongName{igrid});
      netcdf.putAtt(ncid,gvarid(igrid),'standard_name',gstdName{igrid});
      netcdf.putAtt(ncid,gvarid(igrid),'units',gunits{igrid});
   end;
   netcdf.endDef(ncid);

   %Write GRID 
   for igrid=1:length(gname);
      tmp  = gvar{igrid};
      tmpd = size(tmp); tmp = reshape(tmp,fliplr(tmpd));
      netcdf.putVar(ncid,gvarid(igrid),tmp);
      clear tmp tmpd; 
   end;
   
   %CLOSE nc file	+------------------------------------------------------------------------------
   netcdf.close(ncid); 

end;
fprintf('\n');