#include <stdio.h>  
#include <stdlib.h>      
#include <malloc.h>
#include <assert.h>
#include "netcdf.h"


/*-------------------------------------------------------------------------
 * Function: main
 *
 *-------------------------------------------------------------------------
 */

void handle_error(int nc_status);
int read_data( const char* fname, int ndims, int *dims, float **buf );

int main(void)
{
 int  nc_status;                    /* error nc_status */
 int  ncid;                         /* netCDF ID */
 int  lat_dim, lon_dim;             /* dimension IDs */
 int  var_id;                       /* variable ID */
 int  lat_id;                       /* variable ID */
 int  lon_id;                       /* variable ID */
 int  dimids[2];                    /* variable shape */
 float *vals=NULL;                  /* array to hold data values */
 float *latbuf=NULL;                /* array to hold the latitude values */
 float *lonbuf=NULL;                /* array to hold the longitude values */
 float fill=-99;                    /* fill value */
 int   dims[2];                     /* array to hold dimensions */
 int   latdims[1];                  /* array to hold dimensions */
 int   londims[1];                  /* array to hold dimensions */

 /* read the bathymetry data. we read the dimensions and data from an ASCII file */
 if (read_data("omex_data.txt",2,dims,&vals)<0)
  return 1;

 /* read the latitude */
 if (read_data("omex_lat.txt",1,latdims,&latbuf)<0)
  return 1;

 /* read the longitude */
 if (read_data("omex_lon.txt",1,londims,&lonbuf)<0)
  return 1;

/* create the netCDF file */
 nc_status = nc_create("omex.nc", NC_CLOBBER, &ncid);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 /* define dimensions */
 nc_status = nc_def_dim(ncid, "lat", dims[0], &lat_dim);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 nc_status = nc_def_dim(ncid, "lon", dims[1], &lon_dim);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 /* define data variable */
 dimids[0] = lat_dim;
 dimids[1] = lon_dim;
 nc_status = nc_def_var (ncid, "omex", NC_FLOAT, 2, dimids, &var_id);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 /* define fill value */
 nc_status=nc_put_att_float(ncid, var_id, "_FillValue", NC_FLOAT, 1, &fill);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 /* 
 define coordinate variables
 it is legal for a variable to have the same name as a dimension. 
 such variables have no special meaning to the netCDF library.
 however there is a convention that such variables should be treated in a 
 special way by software using this library
 */

 /* define latitude variable */
 nc_status = nc_def_var (ncid, "lat", NC_FLOAT, 1, &lat_dim, &lat_id);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 /* define longitude variable */
 nc_status = nc_def_var (ncid, "lon", NC_FLOAT, 1, &lon_dim, &lon_id);
 if (nc_status != NC_NOERR) handle_error(nc_status);


 /*leave define mode*/
 nc_status=nc_enddef(ncid);  
 if (nc_status != NC_NOERR) handle_error(nc_status);

 /* write data values into netCDF variable */
 nc_status = nc_put_var_float(ncid, var_id, vals);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 /* write latitude values into netCDF variable */
 nc_status = nc_put_var_float(ncid, lat_id, latbuf);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 /* write longitude values into netCDF variable */
 nc_status = nc_put_var_float(ncid, lon_id, lonbuf);
 if (nc_status != NC_NOERR) handle_error(nc_status);


 nc_status=ncclose(ncid);
 if (nc_status != NC_NOERR) handle_error(nc_status);

 if (vals)
  free(vals);
 if (latbuf)
  free(latbuf);
 if (lonbuf)
  free(lonbuf);
 
 return 0;

}


/*-------------------------------------------------------------------------
 * handle_error
 * print an error message
 *-------------------------------------------------------------------------
 */

void handle_error(int nc_status) {
if (nc_status != NC_NOERR) {
   fprintf(stderr, "%s\n", nc_strerror(nc_status));
   exit(-1);
   }
}



/*-------------------------------------------------------------------------
 * read_data
 * utility function to read ASCII data
 * the files have a header of the type
 *
 *   dimension i
 *   n
 * 
 * followed by the data
 *
 *-------------------------------------------------------------------------
 */

int read_data( const char* fname, int ndims, int *dims, float **buf )
{
 int      i, n;
 unsigned j;
 char     str[20];
 size_t   nelms;
 FILE     *f;
 float    val;

 f = fopen( fname, "r");
 if ( f == NULL )
 {
  printf( "Could not open file %s\n", fname );
		return -1;
	}

 for (i=0, nelms=1; i < ndims; i++)
 {
  fscanf( f, "%s %d", str, &j);
  fscanf( f, "%d",&n );
  dims[i] = n;
  nelms *= n;
 }

 *buf = (float*) malloc (nelms * sizeof( float ));

 for (j = 0; j < nelms; j++)
 {
  fscanf( f, "%f",&val );
  (*buf)[j] = val;
 }
 fclose(f);

 return 1;

}