function  oct_create_grid(L,M,grdname,title)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 	Create an empty oct_netcdf gridfile
%       L: total number of psi points in x direction  
%       M: total number of psi points in y direction  
%       grdname: name of the grid file
%       title: title in the oct_netcdf file  
% 
%  Further Information:  
%  http://www.croco-ocean.org
%  
%  This file is part of CROCOTOOLS
%
%  CROCOTOOLS is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published
%  by the Free Software Foundation; either version 2 of the License,
%  or (at your option) any later version.
%
%  CROCOTOOLS is distributed in the hope that it will be useful, but
%  WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
%  MA  02111-1307  USA
%
%  Copyright (c) 2001-2006 by Pierrick Penven 
%  e-mail:Pierrick.Penven@ird.fr  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lp=L+1;
Mp=M+1;

% Create file (clobber existing)
ncid = netcdf.create(grdname,'CLOBBER');

% Define dimensions
dim_xi_u = netcdf.defDim(ncid,'xi_u',L);
dim_eta_u = netcdf.defDim(ncid,'eta_u',Mp);
dim_xi_v = netcdf.defDim(ncid,'xi_v',Lp);
dim_eta_v = netcdf.defDim(ncid,'eta_v',M);
dim_xi_rho = netcdf.defDim(ncid,'xi_rho',Lp);
dim_eta_rho = netcdf.defDim(ncid,'eta_rho',Mp);
dim_xi_psi = netcdf.defDim(ncid,'xi_psi',L);
dim_eta_psi = netcdf.defDim(ncid,'eta_psi',M);
dim_one = netcdf.defDim(ncid,'one',1);
dim_two = netcdf.defDim(ncid,'two',2);
dim_four = netcdf.defDim(ncid,'four',4);
dim_bath = netcdf.defDim(ncid,'bath',1);

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

% Helper to define scalar char
def_char = @(name,long) netcdf.defVar(ncid,name,NC_CHAR,dim_one);
def_double = @(name,dims) netcdf.defVar(ncid,name,NC_DOUBLE,dims);

% Define variables (subset of attributes kept)
var_xl = def_double('xl',dim_one);
netcdf.putAtt(ncid,var_xl,'long_name','domain length in the XI-direction');
netcdf.putAtt(ncid,var_xl,'units','meter');

var_el = def_double('el',dim_one);
netcdf.putAtt(ncid,var_el,'long_name','domain length in the ETA-direction');
netcdf.putAtt(ncid,var_el,'units','meter');

var_depthmin = def_double('depthmin',dim_one);
netcdf.putAtt(ncid,var_depthmin,'long_name','Shallow bathymetry clipping depth');
netcdf.putAtt(ncid,var_depthmin,'units','meter');

var_depthmax = def_double('depthmax',dim_one);
netcdf.putAtt(ncid,var_depthmax,'long_name','Deep bathymetry clipping depth');
netcdf.putAtt(ncid,var_depthmax,'units','meter');

var_spherical = def_char('spherical','Grid type logical switch');
netcdf.putAtt(ncid,var_spherical,'long_name','Grid type logical switch');

var_angle = def_double('angle',[dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_angle,'long_name','angle between xi axis and east');
netcdf.putAtt(ncid,var_angle,'units','radian');

var_h = def_double('h',[dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_h,'long_name','Final bathymetry at RHO-points');
netcdf.putAtt(ncid,var_h,'units','meter');

var_hraw = def_double('hraw',[dim_bath dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_hraw,'long_name','Working bathymetry at RHO-points');
netcdf.putAtt(ncid,var_hraw,'units','meter');

var_alpha = def_double('alpha',[dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_alpha,'long_name','Weights between coarse and fine grids at RHO-points');

var_f = def_double('f',[dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_f,'long_name','Coriolis parameter at RHO-points');
netcdf.putAtt(ncid,var_f,'units','second-1');

var_pm = def_double('pm',[dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_pm,'long_name','curvilinear coordinate metric in XI');
netcdf.putAtt(ncid,var_pm,'units','meter-1');

var_pn = def_double('pn',[dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_pn,'long_name','curvilinear coordinate metric in ETA');
netcdf.putAtt(ncid,var_pn,'units','meter-1');

var_dndx = def_double('dndx',[dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_dndx,'long_name','xi derivative of inverse metric factor pn');
netcdf.putAtt(ncid,var_dndx,'units','meter');

var_dmde = def_double('dmde',[dim_eta_rho dim_xi_rho]);
netcdf.putAtt(ncid,var_dmde,'long_name','eta derivative of inverse metric factor pm');
netcdf.putAtt(ncid,var_dmde,'units','meter');

var_x_rho = def_double('x_rho',[dim_eta_rho dim_xi_rho]);
var_x_u = def_double('x_u',[dim_eta_u dim_xi_u]);
var_x_v = def_double('x_v',[dim_eta_v dim_xi_v]);
var_x_psi = def_double('x_psi',[dim_eta_psi dim_xi_psi]);

var_y_rho = def_double('y_rho',[dim_eta_rho dim_xi_rho]);
var_y_u = def_double('y_u',[dim_eta_u dim_xi_u]);
var_y_v = def_double('y_v',[dim_eta_v dim_xi_v]);
var_y_psi = def_double('y_psi',[dim_eta_psi dim_xi_psi]);

var_lon_rho = def_double('lon_rho',[dim_eta_rho dim_xi_rho]);
var_lon_u = def_double('lon_u',[dim_eta_u dim_xi_u]);
var_lon_v = def_double('lon_v',[dim_eta_v dim_xi_v]);
var_lon_psi = def_double('lon_psi',[dim_eta_psi dim_xi_psi]);

var_lat_rho = def_double('lat_rho',[dim_eta_rho dim_xi_rho]);
var_lat_u = def_double('lat_u',[dim_eta_u dim_xi_u]);
var_lat_v = def_double('lat_v',[dim_eta_v dim_xi_v]);
var_lat_psi = def_double('lat_psi',[dim_eta_psi dim_xi_psi]);

var_mask_rho = def_double('mask_rho',[dim_eta_rho dim_xi_rho]);
var_mask_u = def_double('mask_u',[dim_eta_u dim_xi_u]);
var_mask_v = def_double('mask_v',[dim_eta_v dim_xi_v]);
var_mask_psi = def_double('mask_psi',[dim_eta_psi dim_xi_psi]);

% 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,'type','CROCO grid file');

% Close file
netcdf.close(ncid);