function oct_create_bulk(frcname,grdname,title,bulkt,bulkc,Yorig)
% oct_create_bulk: create an empty bulk forcing file using octave-netcdf

% Read grid to get sizes
ng = netcdf.open(grdname,'NC_NOWRITE');
L = netcdf.inqDim(ncIdToUse=netcdf.inqDimID(ng,'xi_psi'));
% Above is a no-op used to keep compatibility with older Octave; instead query dims:
L = netcdf.inqDim(ng, netcdf.inqDimID(ng,'xi_psi'));
M = netcdf.inqDim(ng, netcdf.inqDimID(ng,'eta_psi'));
% netcdf.inqDim returns [name,len], some Octave versions return two outputs
try
  [~,Llen] = netcdf.inqDim(ng, netcdf.inqDimID(ng,'xi_psi'));
  [~,Mlen] = netcdf.inqDim(ng, netcdf.inqDimID(ng,'eta_psi'));
  L = Llen; M = Mlen;
catch
  % fallback: assume previous values
end
netcdf.close(ng);
Lp = L+1;
Mp = M+1;

% time attribute helpers
[insecond,time_unit_att,time_second_unit_att,calendar_att] = deal(0,'days since 0001-01-01','seconds since 0001-01-01','gregorian');
if ~exist('Yorig','var') , Yorig = [] ; end
[time_unit_att,time_second_unit_att,calendar_att] = get_time_attributes(Yorig);

% Create file
ncid = netcdf.create(frcname,'NC_CLOBBER');
% Define dimensions
netcdf.defDim(ncid,'xi_rho',Lp);
netcdf.defDim(ncid,'eta_rho',Mp);
netcdf.defDim(ncid,'xi_psi',L);
netcdf.defDim(ncid,'eta_psi',M);
netcdf.defDim(ncid,'xi_u',L);
netcdf.defDim(ncid,'eta_u',Mp);
netcdf.defDim(ncid,'xi_v',Lp);
netcdf.defDim(ncid,'eta_v',M);
netcdf.defDim(ncid,'bulk_time',0);

% Define variables (double)
var_bulk_time = netcdf.defVar(ncid,'bulk_time','NC_DOUBLE','bulk_time');
netcdf.putAtt(ncid,var_bulk_time,'long_name','bulk formulation execution time');
% write other var definitions
netcdf.defVar(ncid,'tair','NC_DOUBLE', ['bulk_time','eta_rho','xi_rho']);
netcdf.putAtt(ncid,netcdf.inqVarID(ncid,'tair'),'long_name','surface air temperature');
netcdf.putAtt(ncid,netcdf.inqVarID(ncid,'tair'),'units','Celsius');
% For simplicity define remaining vars using dimension IDs
% Get dim IDs
d_btime = netcdf.inqDimID(ncid,'bulk_time');
d_eta_rho = netcdf.inqDimID(ncid,'eta_rho');
d_xi_rho = netcdf.inqDimID(ncid,'xi_rho');
d_eta_u = netcdf.inqDimID(ncid,'eta_u');
d_xi_u = netcdf.inqDimID(ncid,'xi_u');
d_eta_v = netcdf.inqDimID(ncid,'eta_v');
d_xi_v = netcdf.inqDimID(ncid,'xi_v');

% use arrays of dim IDs for multi-dim vars
netcdf.defVar(ncid,'rhum','NC_DOUBLE',[d_btime d_eta_rho d_xi_rho]);
netcdf.putAtt(ncid,netcdf.inqVarID(ncid,'rhum'),'long_name','relative humidity');
netcdf.putAtt(ncid,netcdf.inqVarID(ncid,'rhum'),'units','fraction');

netcdf.defVar(ncid,'prate','NC_DOUBLE',[d_btime d_eta_rho d_xi_rho]);
netcdf.putAtt(ncid,netcdf.inqVarID(ncid,'prate'),'long_name','precipitation rate');
netcdf.putAtt(ncid,netcdf.inqVarID(ncid,'prate'),'units','cm day-1');

netcdf.defVar(ncid,'wspd','NC_DOUBLE',[d_btime d_eta_rho d_xi_rho]);
netcdf.defVar(ncid,'radlw','NC_DOUBLE',[d_btime d_eta_rho d_xi_rho]);
netcdf.defVar(ncid,'radlw_in','NC_DOUBLE',[d_btime d_eta_rho d_xi_rho]);
netcdf.defVar(ncid,'radsw','NC_DOUBLE',[d_btime d_eta_rho d_xi_rho]);
netcdf.defVar(ncid,'sustr','NC_DOUBLE',[d_btime d_eta_u d_xi_u]);
netcdf.defVar(ncid,'svstr','NC_DOUBLE',[d_btime d_eta_v d_xi_v]);
netcdf.defVar(ncid,'uwnd','NC_DOUBLE',[d_btime d_eta_u d_xi_u]);
netcdf.defVar(ncid,'vwnd','NC_DOUBLE',[d_btime d_eta_v d_xi_v]);

% Global attributes
netcdf.putAtt(ncid, netcdf.getConstant('NC_GLOBAL'), 'title', title);
netcdf.putAtt(ncid, netcdf.getConstant('NC_GLOBAL'), 'date', date);
netcdf.putAtt(ncid, netcdf.getConstant('NC_GLOBAL'), 'grd_file', grdname);
netcdf.putAtt(ncid, netcdf.getConstant('NC_GLOBAL'), 'type', 'CROCO heat flux bulk forcing file');

netcdf.endDef(ncid);

% Write time variable bulkt
varid = netcdf.inqVarID(ncid,'bulk_time');
netcdf.putVar(ncid,varid,0,length(bulkt),bulkt);

netcdf.close(ncid);
return
