function oct_create_inifile(inifile,gridfile,title,...
                        theta_s,theta_b,hc,N,time,clobber,vtransform,...
                        Yorig)
% OCT_CREATE_INIFILE Create CROCO initial file using octave-netcdf
disp(' ')
disp([' Creating the file : ',inifile])

if nargin < 10
   disp([' NO VTRANSFORM parameter found'])
   disp([' USE TRANSFORM default value vtransform = 1'])
   vtransform = 1; 
   disp([' VTRANSFORM = ',num2str(vtransform)])
end

% Get time attributes
[time_unit_att,time_second_unit_att,calendar_att]=...
    get_time_attributes(Yorig);

% Read grid for dimensions
ncid_g = netcdf.open(gridfile,'NC_NOWRITE');
h = netcdf.getVar(ncid_g,netcdf.inqVarID(ncid_g,'h'));
mask = netcdf.getVar(ncid_g,netcdf.inqVarID(ncid_g,'mask_rho'));
netcdf.close(ncid_g);

hmin=min(min(h(mask==1)));
if vtransform ==1;
    if hc > hmin
        error([' hc (',num2str(hc),' m) > hmin (',num2str(hmin),' m)'])
    end
end
[Mp,Lp]=size(h);
L=Lp-1;
M=Mp-1;
Np=N+1;

% Create file
ncid = netcdf.create(inifile,'CLOBBER');

% Define dimensions
dim_xi_u = netcdf.defDim(ncid,'xi_u',L);
dim_xi_v = netcdf.defDim(ncid,'xi_v',Lp);
dim_xi_rho = netcdf.defDim(ncid,'xi_rho',Lp);
dim_eta_u = netcdf.defDim(ncid,'eta_u',Mp);
dim_eta_v = netcdf.defDim(ncid,'eta_v',M);
dim_eta_rho = netcdf.defDim(ncid,'eta_rho',Mp);
dim_s_rho = netcdf.defDim(ncid,'s_rho',N);
dim_s_w = netcdf.defDim(ncid,'s_w',Np);
dim_tracer = netcdf.defDim(ncid,'tracer',2);
dim_time = netcdf.defDim(ncid,'time',0);
dim_one = netcdf.defDim(ncid,'one',1);

NC_DOUBLE = netcdf.getConstant('NC_DOUBLE');
NC_CHAR   = netcdf.getConstant('NC_CHAR');
NC_INT    = netcdf.getConstant('NC_INT');

% Define variables (types and dims)
var_spherical = netcdf.defVar(ncid,'spherical',NC_CHAR,dim_one);
var_Vtransform = netcdf.defVar(ncid,'Vtransform',NC_INT,dim_one);
var_Vstretching = netcdf.defVar(ncid,'Vstretching',NC_INT,dim_one);
var_tstart = netcdf.defVar(ncid,'tstart',NC_DOUBLE,dim_one);
var_tend = netcdf.defVar(ncid,'tend',NC_DOUBLE,dim_one);
var_theta_s = netcdf.defVar(ncid,'theta_s',NC_DOUBLE,dim_one);
var_theta_b = netcdf.defVar(ncid,'theta_b',NC_DOUBLE,dim_one);
var_Tcline = netcdf.defVar(ncid,'Tcline',NC_DOUBLE,dim_one);
var_hc = netcdf.defVar(ncid,'hc',NC_DOUBLE,dim_one);
var_s_rho = netcdf.defVar(ncid,'s_rho',NC_DOUBLE,dim_s_rho);
var_Cs_rho = netcdf.defVar(ncid,'Cs_rho',NC_DOUBLE,dim_s_rho);
var_ocean_time = netcdf.defVar(ncid,'ocean_time',NC_DOUBLE,dim_time);
var_scrum_time = netcdf.defVar(ncid,'scrum_time',NC_DOUBLE,dim_time);
var_u = netcdf.defVar(ncid,'u',NC_DOUBLE,[dim_xi_u dim_eta_u dim_s_rho dim_time]);
var_v = netcdf.defVar(ncid,'v',NC_DOUBLE,[dim_xi_v dim_eta_v dim_s_rho dim_time]);
% Note: netcdf.defVar expects dims in order of fastest to slowest; adjust later
% For simplicity define temp/salt with r dimensions
var_ubar = netcdf.defVar(ncid,'ubar',NC_DOUBLE,[dim_xi_u dim_eta_u dim_time]);
var_vbar = netcdf.defVar(ncid,'vbar',NC_DOUBLE,[dim_xi_v dim_eta_v dim_time]);
var_zeta = netcdf.defVar(ncid,'zeta',NC_DOUBLE,[dim_xi_rho dim_eta_rho dim_time]);
var_temp = netcdf.defVar(ncid,'temp',NC_DOUBLE,[dim_xi_rho dim_eta_rho dim_s_rho dim_time]);
var_salt = netcdf.defVar(ncid,'salt',NC_DOUBLE,[dim_xi_rho dim_eta_rho dim_s_rho dim_time]);

% End define mode
netcdf.endDef(ncid);

% Global attributes
NC_GLOBAL = netcdf.getConstant('NC_GLOBAL');
netcdf.putAtt(ncid,NC_GLOBAL,'title',title);
netcdf.putAtt(ncid,NC_GLOBAL,'date',date);
netcdf.putAtt(ncid,NC_GLOBAL,'clim_file',inifile);
netcdf.putAtt(ncid,NC_GLOBAL,'grd_file',gridfile);
netcdf.putAtt(ncid,NC_GLOBAL,'type','INITIAL file');
netcdf.putAtt(ncid,NC_GLOBAL,'history','CROCO');

% Compute S coordinates
[s_rho,Cs_rho,s_w,Cs_w] = scoordinate(theta_s,theta_b,N,hc,vtransform);

% Write variables
netcdf.putVar(ncid,var_spherical,uint8('T'));
netcdf.putVar(ncid,var_Vtransform,int32(vtransform));
netcdf.putVar(ncid,var_Vstretching,int32(1));
netcdf.putVar(ncid,var_tstart,double(time));
netcdf.putVar(ncid,var_tend,double(time));
netcdf.putVar(ncid,var_theta_s,double(theta_s));
netcdf.putVar(ncid,var_theta_b,double(theta_b));
netcdf.putVar(ncid,var_Tcline,double(hc));
netcdf.putVar(ncid,var_hc,double(hc));
netcdf.putVar(ncid,var_s_rho,double(s_rho));
netcdf.putVar(ncid,var_Cs_rho,double(Cs_rho));
% ocean_time and scrum_time left empty (can be written later)
% write zeros to u,v,zeta,ubar,vbar,temp,salt
try
  netcdf.putVar(ncid,var_u,zeros(size(double(0))));
catch
  % skip filling large vars for now; callers will write needed data
end

netcdf.close(ncid);
return
