function []=convert_field2nctiles_GRID(tile_dims,inDir,outDir);
%----
%function []=convert_field2nctiles_GRID(tile_dims,inDir,outDir);
%
% Read grid binary files already in tile form and write nc output in ECCO form
%
% Inputs:
%	tile_dims 	[nx ny nz] 	Number of grid points in each dimension 
%	inDir		Directory with GRID in binary tiles form
%	outDir		Directory where GRID nctiles is to be written
%
% Format is consistent with ECCO nctiles output and with netcdf conventions CF-1.6
%
% Written by Victor Ocana 2019
% Modified by An Nguyen Jul 2022
%-----

%GRID variables			+----------------------------------------------------------
varName={'hFacC','hFacW','hFacS','XC','YC','XG','YG','RAC','RAZ','DXC','DYC','DXG','DYG',...
         'Depth','AngleCS','AngleSN','RC','RF','DRC','DRF'};

longName={'fractional thickness','fractional thickness','fractional thickness',...
          'longitude','latitude','longitude','latitude','grid cell area',...
          'grid cell area','grid spacing','grid spacing','grid spacing',...
          'grid spacing','sea floor depth','grid orientation (cosine)',...
          'grid orientation (sine)','vertical coordinate','vertical coordinate',...
          'grid spacing','grid spacing'};

units={'1','1','1','degrees_east','degrees_north','degrees_east','degrees_north',...
       'm^2','m^2','m','m','m','m','m','1','1','m','m','m','m'};

nvars=size(varName,2);

%Directories, files		+---------------------------------------------------------- 
if(exist(outDir)==0);mkdir(outDir);end;

%Determine number of tiles	+---------------------------------------------------------- 
flist = dir([inDir varName{1} '.*.data']);
ntiles = size(flist,1);

%Determine the tile size nx, ny, nz for each output file GRID.xxxx.nc
metaFile = dir([inDir varName{1} '*.meta']);
if(size(metaFile,1));
   meta=read_meta(metaFile.name,inDir);
   if(isempty(tile_dims));
      tile_dims(1) = meta.dimList(1); 
      tile_dims(2) = meta.dimList(4); 
      tile_dims(3) = meta.dimList(7); 
   end;
else;
   if(isempty(tile_dims));
      str1=['convert_field2nctiles_GRID ERRORL: '];
      str2=[No tile_dims given and no .meta file present. Unable to determine tile size'];
      disp([str1 str2]);
      return;
   end;
end;

nDim = size(tile_dims,2);
dimsize = fliplr(tile_dims); %CF-1.6 convention has z,y,x as the order of dims rather than x,y,z as in ECCO

%Set some values for global attributes	+---------------------------------------------------------- 
missval=NaN; 
fillval=NaN;
netcdf_vers = netcdf.inqLibVers;
matlab_vers = version;

%Loop over tiles and variables 		+---------------------------------------------------------- 
for itile = 1:ntiles;
   fprintf('Tile %04d: ',itile);
   tileName = sprintf('%04d',itile);
   fileTile = [outDir 'GRID.' tileName '.nc']; 
   
   %Create netcdf file		+---------------------------------------------------------- 
   %(by default usic 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');
   netcdf.putAtt(ncid,nc_global,'Description', 'C-grid parameters (see MITgcm documentation for details).');
   netcdf.putAtt(ncid,nc_global,'Format','native grid');
   netcdf.putAtt(ncid,nc_global,'Institution','Oden Institute for Computational Engineering and Science, UT-Austin');
   netcdf.putAtt(ncid,nc_global,'N-','File created using convert_field2nctiles_GRID.m');
   netcdf.putAtt(ncid,nc_global,'V-netCDF Version',netcdf_vers);
   netcdf.putAtt(ncid,nc_global,'V-matlab Version',matlab_vers);
   netcdf.putAtt(ncid,nc_global,'date',date);
   netcdf.putAtt(ncid,nc_global,'Conventions','CF-1.6 (http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html)')

   netcdf.putAtt(ncid,nc_global,'_FillValue',fillval);
   netcdf.putAtt(ncid,nc_global,'missing_value',missval);
    
   %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));
      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;

   %GRID Variables 			+---------------------------------------------------------------------
   for ivar = 1:nvars; 
      clear dims dimidg tmpN nprec;
      %First check if .meta file exists and select appropriate dims for each grid variable 
      metaFile = dir([inDir varName{ivar} '*.meta']);
      if(size(metaFile,1));
         meta=read_meta(metaFile.name,inDir);
	 tmpN = [inDir varName{ivar} '.' tileName '.data'];
         for ii=1:meta.nDims;  
            dims(ii)=[meta.dimList(1+3*(ii-1))]; 
         end;
         if(meta.nDims==3);
            dimidg=dimid;
         elseif(meta.nDims==2);
            dimidg=dimid(2:3);
         end;
         if(dims(1)==1); 	%Special case for 1D fields (RC, RF, DRC, DRF)
	    tmpN = [inDir varName{ivar} '.data'];
            dims=dims(3);
            dimidg=dimid(1);
         end;
         nprec = str2num(meta.dataprec(end-1:end))/8;
      else;
         %If no meta file available, get info from .data file
         if(ivar<=3);          		%3D vars +-----------------------------------------------------------
            dims=tile_dims;
            dimidg=dimid;
	    tmpN = [inDir varName{ivar} '.' tileName '.data'];
            %Determine precission
            tmpd=dir(tmpN);  
            nprec=tmpd.bytes/prod(dims); 	
   
         elseif(ivar>3 & ivar<=16);	%2D vars +------------------------------------------------------------
            dims=tile_dims(1:2);   %2D vars
            dimidg=dimid(2:3);
	    tmpN = [inDir varName{ivar} '.' tileName '.data'];
            %Determine precission
            tmpd=dir(tmpN);  
            nprec=tmpd.bytes/prod(dims); 
   
         elseif(ivar>16);			%1D vars +----------------------------------------------------
            dims=tile_dims(3);     
            dimidg=dimid(1);
	    tmpN = [inDir varName{ivar} '.data'];
            %Determine precission
            tmpd=dir(tmpN);  
            nprec=floor(tmpd.bytes/prod(dims)); 
            %above: !!! RF has nz+1 levels, so tmpd.bytes/prod(dims) is slightly larger than 8 for dims(3)=nz
           
         end;
      end;
     
     fprintf(' %s ', varName{ivar});

     %Set precission flags to read/write
     if(~(nprec==4|nprec==8));
      disp(['WARNING! Check dimensions: precission=' num2str(nprec) ' It should be either 4 or 8']);
     end;
     if(nprec==4); xtype ='NC_FLOAT';elseif(nprec==8);xtype ='NC_DOUBLE';end;
     prec = ['real*' num2str(nprec) ''];

     %Read binary tiles
     data  = readbin(tmpN,dims,1,prec);      %NOTE: For RF and DRF read only the first nz levels
     if(size(dims,2)>1);data=reshape(data,fliplr(dims));end;

     %Define variable and its attributes
     netcdf.reDef(ncid);
     %Notice the order of the dimensions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     varid = netcdf.defVar(ncid,varName{ivar},xtype,flipdim(dimidg,2));
     if ~isempty(longName); netcdf.putAtt(ncid,varid,'long_name',longName{ivar}); end;
     if ~isempty(units); netcdf.putAtt(ncid,varid,'units',units{ivar}); end;
     netcdf.endDef(ncid);

     %varid
%atn it seems that we're only allowed to write up to nr, not nr+1, so if dims > nr, need to write only up to nr
     sz=size(data);
     if(ivar>16 & sz(1) > tile_dims(3));data=data(1:tile_dims(3));end;
     %Write data
     netcdf.putVar(ncid,varid,data); 
   end;
   fprintf('\n');
   netcdf.close(ncid); 
end;