%
% This is the name of the data file we will create.
ncfile='pres_temp_4D.nc';

%
% We are writing 4D data, a 2 x 6 x 12 lvl-lat-lon grid, with 2
% timesteps of data. 
NDIMS=4;
NLAT=6;
NLON=12;
NREC=2;
NLVL=2;

% 
% These are used to construct some example data. 
SAMPLE_PRESSURE = 900;
SAMPLE_TEMP = 9.0;
START_LAT = 25.0;
START_LON = -125.0;

%
% Create some pretend data. If this wasn't an example program, we
% would have some real data to write, for example, model
% output.
lats = [25:0.5:27.5];
NLAT = length(lats);
lons = [-125:0.5:-119.5];
NLON = length(lons);

pres_out = ones(NLVL,NLAT,NLON);
temp_out = ones(NLVL,NLAT,NLON);

%
% This is a little inefficient from a matlab point of view, but 
% accomplishes the same thing as the corresponding C program.
inc = 0;
for lvl = 1:NLVL
    for lat = 1:NLAT
        for lon = 1:NLON
            pres_out(lvl,lat,lon) = SAMPLE_PRESSURE + inc;
            temp_out(lvl,lat,lon) = SAMPLE_TEMP + inc;
            inc = inc + 1;
        end
    end
end

%
% Create the file.
[ncid, status] = mexnc ( 'CREATE', ncfile, nc_clobber_mode );
if status, error(mexnc('STRERROR',status)), end

%
% Define the dimensions. The record dimension is defined to have
% unlimited length - it can grow as needed. In this example it is
% the time dimension.
[lvl_dimid, status] = mexnc ( 'DEF_DIM', ncid, 'level', NLVL );
if status, error(mexnc('STRERROR',status)), end
[lat_dimid, status] = mexnc ( 'DEF_DIM', ncid, 'latitude', NLAT );
if status, error(mexnc('STRERROR',status)), end
[lon_dimid, status] = mexnc ( 'DEF_DIM', ncid, 'longitude', NLON );
if status, error(mexnc('STRERROR',status)), end
[rec_dimid, status] = mexnc ( 'DEF_DIM', ncid, 'time', 0 );
if status, error(mexnc('STRERROR',status)), end

%
% Define the coordinate variables. We will only define coordinate
% variables for lat and lon.  
[lat_varid, status] = mexnc ( 'DEF_VAR', ncid, 'latitude', nc_float, 1, lat_dimid );
if status, error(mexnc('STRERROR',status)), end
[lon_varid, status] = mexnc ( 'DEF_VAR', ncid, 'longitude', nc_float, 1, lon_dimid );
if status, error(mexnc('STRERROR',status)), end


%
% Assign units attributes to coordinate variables.  This could also be
% done with "nc_attput" in SNCTOOLS after we close the file.  
status = mexnc ( 'PUT_ATT_TEXT', ncid, lat_varid, 'units', nc_char, length('degrees_north'), 'degrees_north' );
if status, error(mexnc('STRERROR',status)), end
status = mexnc ( 'PUT_ATT_TEXT', ncid, lon_varid, 'units', nc_char, length('degrees_east'), 'degrees_east' );
if status, error(mexnc('STRERROR',status)), end

% 
% The dimids array is used to pass the dimids of the dimensions of
% the netCDF variables. Both of the netCDF variables we are
% creating share the same four dimensions. In both C and MATLAB, the
% unlimited dimension must come first on the list of dimids. 
dimids = [ rec_dimid lvl_dimid lat_dimid lon_dimid ];

   
% 
% Define the netCDF variables for the pressure and temperature
% data. 
[pres_varid, status] = mexnc ( 'DEF_VAR', ncid, 'pressure', nc_float, 4, dimids );
if status, error(mexnc('STRERROR',status)), end
[temp_varid, status] = mexnc ( 'DEF_VAR', ncid, 'temperature', nc_float, 4, dimids );
if status, error(mexnc('STRERROR',status)), end

%
% Assign units attributes to the netCDF variables.
status = mexnc ( 'PUT_ATT_TEXT', ncid, pres_varid, 'units', nc_char, length('hPa'), 'hPa' );
if status, error(mexnc('STRERROR',status)), end
status = mexnc ( 'PUT_ATT_TEXT', ncid, temp_varid, 'units', nc_char, length('celcius'), 'celcius' );
if status, error(mexnc('STRERROR',status)), end

%
% Close the file.  We then use SNCTOOLS to write the data.
status = mexnc ( 'CLOSE', ncid );
if status, error(mexnc('STRERROR',status)), end


%
% Write the coordinate variable data. This will put the latitudes
% and longitudes of our data grid into the netCDF file. 
nc_varput ( ncfile, 'latitude', lats );
nc_varput ( ncfile, 'longitude', lons );

%
% These settings tell netcdf to write one timestep of data. (The
% setting of start(1) inside the loop below tells netCDF which
% timestep to write.) 
count = [1 NLVL NLAT NLON ];
start = [0    0    0    0 ];

%
% Write the pretend data. This will write our surface pressure and
% surface temperature data. The arrays only hold one timestep worth
% of data. We will just rewrite the same data for each timestep. In
% a real application, the data would change between timesteps. 
%
% One little gotcha here is that nc_varput expects the rank of the
% input data to match the rank of the netcdf variable.  The input data
% here has three dimensions, but the netcdf variables each have four.
% We get around this by stuffing the input data into a 4D matrix with
% a leading dimension extent of one.
for rec = 0:NREC-1

    start(1) = rec;

    out_data(1,:,:,:) = pres_out;
    nc_varput ( ncfile, 'pressure', out_data, start, count );

    out_data(1,:,:,:) = temp_out;
    nc_varput ( ncfile, 'temperature', out_data, start, count );

end


fprintf(1, '*** SUCCESS writing example file %s!\n', ncfile);