xref: /dports/graphics/gmt/gmt-6.3.0/src/gmt_nc.c

/*--------------------------------------------------------------------
 *
 *	Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
 *	See LICENSE.TXT file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU Lesser General Public License as published by
 *	the Free Software Foundation; version 3 or any later version.
 *
 *	This program 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 Lesser General Public License for more details.
 *
 *	Contact info: www.generic-mapping-tools.org
 *--------------------------------------------------------------------*/
/*
 *
 *	G M T _ N C . C   R O U T I N E S
 *
 * Takes care of all grd input/output built on NCAR's NetCDF routines.
 * This version is intended to provide more general support for reading
 * NetCDF files that were not generated by GMT. At the same time, the grids
 * written by these routines are intended to be more conform COARDS conventions.
 * These routines are to eventually replace the older gmt_cdf_ routines.
 *
 * Most functions will return with error message if an internal error is returned.
 * There functions are only called indirectly via the GMT_* grdio functions.
 *
 * Author:  Remko Scharroo
 * Date:    04-AUG-2005
 * Version: 1
 *
 * Added support for chunked I/O, Florian Wobbe, June 2012.
 *
 * Functions include:
 *
 *  gmt_nc_read_grd_info:   Read header from file
 *  gmt_nc_read_grd:        Read data set from file
 *  gmt_nc_update_grd_info: Update header in existing file
 *  gmt_nc_write_grd_info:  Write header to new file
 *  gmt_nc_write_grd:       Write header and data set to new file
 *  gmt_nc_read_cube_info:  Read information from cube file
 *  gmt_nc_write_cube:      rite header and cube to new file(s)
 *  gmtlib_is_nc_grid:	    Determine if we have a nc grid
 *
 * Private functions:
 *  gmtnc_setup_chunk_cache:      Change the default HDF5 chunk cache settings
 *  gmtnc_pad_grid:               Add padding to a grid
 *  gmtnc_unpad_grid:             Remove padding from a grid
 *  gmtnc_padding_copy:           Fill padding by replicating the border cells
 *  gmtnc_padding_zero:           Fill padding with zeros
 *  gmtnc_n_chunked_rows_in_cache Determines how many chunks to read at once
 *  gmtnc_io_nc_grid              Does the actual netcdf I/O
 *  gmtnc_netcdf_libvers          returns the netCDF library version
 *  gmtnc_right_shift_grid
 *  gmtnc_set_optimal_chunksize   Determines the optimal chunksize
 *  gmtnc_get_units
 *  gmtnc_put_units
 *  gmtnc_check_step
 *  gmtnc_grd_info
 *  gmtnc_grid_fix_repeat_col
 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/

#include "gmt_dev.h"
#include "gmt_internals.h"
#include <netcdf.h>

/* Declaration modifier for netcdf DLL support
 * annoying: why can't netcdf.h do this on its own? */
#if defined WIN32 && ! defined NETCDF_STATIC
#define DLL_NETCDF
#endif

/* HDF5 chunk cache: reasonable defaults assuming min. chunk size of 128x128 and type byte */
#define NC_CACHE_SIZE       33554432 /* 32MiB */
#define NC_CACHE_NELEMS     2053     /* prime > NC_CACHE_SIZE / (128*128*1byte) */
#define NC_CACHE_PREEMPTION 0.75f

int gmt_cdf_grd_info (struct GMT_CTRL *GMT, int ncid, struct GMT_GRID_HEADER *header, char job);
int gmt_cdf_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode);

static int nc_libvers[] = {-1, -1, -1, -1}; /* holds the version of the netCDF library */

static char *regtype[2] = {"gridline", "pixel"};

GMT_LOCAL const int * gmtnc_netcdf_libvers (void) {
	static bool inquired = false;

	if (!inquired) {
		const char *vers_string = nc_inq_libvers();
		sscanf (vers_string, "%d.%d.%d.%d",
				nc_libvers, nc_libvers+1, nc_libvers+2, nc_libvers+3);
		inquired = true;
	}

	return nc_libvers; /* return pointer to version array */
}

/* netcdf I/O mode */
enum Netcdf_io_mode {
	k_put_netcdf = 0,
	k_get_netcdf
};

/* Wrapper around gmt_nc_put_vara_grdfloat and gmt_nc_get_vara_grdfloat */
static inline int gmtnc_vara_grdfloat (int ncid, int varid, const size_t *startp, const size_t *countp, gmt_grdfloat *fp, unsigned io_mode) {
	if (io_mode == k_put_netcdf)	/* write netcdf */
		return gmt_nc_put_vara_grdfloat (ncid, varid, startp, countp, fp);
	/* read netcdf */
	return gmt_nc_get_vara_grdfloat (ncid, varid, startp, countp, fp);
}

/* Wrapper around gmt_nc_put_varm_grdfloat and gmt_nc_get_varm_grdfloat */
static inline int gmtnc_varm_grdfloat (int ncid, int varid, const size_t *startp, const size_t *countp, const ptrdiff_t *stridep, const ptrdiff_t *imapp, gmt_grdfloat *fp, unsigned io_mode) {
	if (io_mode == k_put_netcdf)	/* write netcdf */
		return gmt_nc_put_varm_grdfloat (ncid, varid, startp, countp, stridep, imapp, fp);
	/* read netcdf */
	return gmt_nc_get_varm_grdfloat (ncid, varid, startp, countp, stridep, imapp, fp);
}

/* Get number of chunked rows that fit into cache (32MiB) */
GMT_LOCAL int gmtnc_n_chunked_rows_in_cache (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, unsigned width, unsigned height, size_t *n_contiguous_chunk_rows, size_t *chunksize) {
	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
	nc_type z_type;		/* type of z variable */
	size_t z_size;		/* size of z variable */
	size_t z_bytes;		/* Number of bytes */
	size_t width_t = (size_t)width;
	size_t height_t = (size_t)height;
	unsigned yx_dim[2] = {HH->xy_dim[1], HH->xy_dim[0]}; /* because xy_dim not row major */
	int err, storage_in;

	gmt_M_err_trap (nc_inq_vartype(HH->ncid, HH->z_id, &z_type)); /* type of z */
	gmt_M_err_trap (nc_inq_type(HH->ncid, z_type, NULL, &z_size)); /* size of z elements in bytes */
	gmt_M_err_trap (nc_inq_var_chunking (HH->ncid, HH->z_id, &storage_in, chunksize));
	if (storage_in != NC_CHUNKED) {
		/* default if NC_CONTIGUOUS */
		chunksize[yx_dim[0]] = 128;   /* 128 rows */
		chunksize[yx_dim[1]] = width_t; /* all columns */
	}

	z_bytes = height_t * height_t * z_size;
	if (z_bytes > NC_CACHE_SIZE) {
		/* memory needed for subset exceeds the cache size */
		unsigned int level;
		size_t chunks_per_row = (size_t) ceil ((double)width / chunksize[yx_dim[1]]);
		*n_contiguous_chunk_rows = NC_CACHE_SIZE / (width_t * z_size) / chunksize[yx_dim[0]];
#ifdef NC4_DEBUG
		level = GMT_MSG_WARNING;
#else
		level = GMT_MSG_DEBUG;
#endif
		GMT_Report (GMT->parent, level,
				"processing at most %" PRIuS " (%" PRIuS "x%" PRIuS ") chunks at a time (%.1lf MiB)...\n",
				*n_contiguous_chunk_rows * chunks_per_row,
				*n_contiguous_chunk_rows, chunks_per_row,
				*n_contiguous_chunk_rows * z_size * width * chunksize[yx_dim[0]] / 1048576.0);
	}
	else
		*n_contiguous_chunk_rows = 0; /* all chunks fit into cache */

	return GMT_NOERROR;
}

/* Read and write classic or chunked netcdf files */
GMT_LOCAL int gmtnc_io_nc_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, unsigned dim[], unsigned origin[], size_t stride, unsigned io_mode, gmt_grdfloat* grid, bool cube) {
	/* io_mode = k_get_netcdf: read a netcdf file to grid
	 * io_mode = k_put_netcdf: write a grid to netcdf */
	int status = NC_NOERR;
	unsigned int d_off = (cube) ? 1 : 0;
	unsigned width = dim[1+d_off], height = dim[0+d_off];
	unsigned yx_dim[2];  /* because xy_dim is not row major! */
	size_t width_t = (size_t)width;
	size_t height_t = (size_t)height;
	size_t chunksize[5]; /* chunksize of z */
	size_t start[5] = {0,0,0,0,0}, count[5] = {1,1,1,1,1};
	size_t n_contiguous_chunk_rows = 0;  /* that are processed at once, 0 = all */
	ptrdiff_t imap[5] = {1,1,1,1,1}; /* mapping between dims of netCDF and in-memory grid */
	const ptrdiff_t onestride[5] = {1,1,1,1,1};	/* Passing this instead of NULL bypasses netCDF bug in 4.6.2 */
	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);

	/* catch illegal io_mode in debug */
	assert (io_mode == k_put_netcdf || io_mode == k_get_netcdf);

#ifdef NC4_DEBUG
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "%s n_columns:%u n_rows:%u x0:%u y0:%u y-order:%s\n",
			io_mode == k_put_netcdf ? "writing," : "reading,",
			dim[1+d_off], dim[0+d_off], origin[1], origin[0],
			HH->row_order == k_nc_start_south ? "S->N" : "N->S");
#endif

	/* set index of input origin */
	yx_dim[0] = HH->xy_dim[1], yx_dim[1] = HH->xy_dim[0]; /* xy_dim not row major */
	memcpy (start, HH->t_index, 3 * sizeof(size_t)); /* set lower dimensions first (e.g. layer) */
	start[yx_dim[0]] = origin[0]; /* first row */
	start[yx_dim[1]] = origin[1]; /* first col */

	/* set mapping of complex grids or if reading a part of a grid */
	imap[yx_dim[0]] = (stride == 0 ? width : stride); /* distance between each row */
	imap[yx_dim[1]] = 1;                              /* distance between values in a row */
	if (cube) imap[0] = header->nm;                   /* distance between values in a layer (?) */

	/* determine how many chunks to process at once */
	gmtnc_n_chunked_rows_in_cache (GMT, header, width, height, &n_contiguous_chunk_rows, chunksize);

	if (n_contiguous_chunk_rows) {
		/* read/write grid in chunks to keep memory footprint low */
		size_t remainder;
#ifdef NC4_DEBUG
		unsigned row_num = 0;
			GMT_Report (GMT->parent, GMT_MSG_WARNING, "stride: %u width: %u\n",
					stride, width);
#endif

		/* adjust row count, so that it ends on the bottom of a chunk */
		count[yx_dim[0]] = chunksize[yx_dim[0]] * n_contiguous_chunk_rows;
		remainder = start[yx_dim[0]] % chunksize[yx_dim[0]];
		count[yx_dim[0]] -= remainder;

		count[yx_dim[1]] = width_t;
		while ( (start[yx_dim[0]] + count[yx_dim[0]]) <= height_t && status == NC_NOERR) {
#ifdef NC4_DEBUG
			GMT_Report (GMT->parent, GMT_MSG_WARNING, "chunked row #%u start-y:%" PRIuS " height:%" PRIuS "\n",
					++row_num, start[yx_dim[0]], count[yx_dim[0]]);
#endif
			/* get/put chunked rows */
			if (stride)
				status = gmtnc_varm_grdfloat (HH->ncid, HH->z_id, start, count, onestride, imap, grid, io_mode);
			else
				status = gmtnc_vara_grdfloat (HH->ncid, HH->z_id, start, count, grid, io_mode);

			/* advance grid location and set new origin */
			grid += count[yx_dim[0]] * ((stride == 0 ? width_t : stride));
			start[yx_dim[0]] += count[yx_dim[0]];
			if (remainder) {
				/* reset count to full chunk height */
				count[yx_dim[0]] += remainder;
				remainder = 0;
			}
		}
		if ( start[yx_dim[0]] != height_t && status == NC_NOERR ) {
			/* get/put last chunked row */
			count[yx_dim[0]] = height_t - start[yx_dim[0]] + origin[0];
#ifdef NC4_DEBUG
			GMT_Report (GMT->parent, GMT_MSG_WARNING, "chunked row #%u start-y:%" PRIuS " height:%" PRIuS "\n",
					++row_num, start[yx_dim[0]], count[yx_dim[0]]);
#endif
			if (stride)
				status = gmtnc_varm_grdfloat (HH->ncid, HH->z_id, start, count, onestride, imap, grid, io_mode);
			else
				status = gmtnc_vara_grdfloat (HH->ncid, HH->z_id, start, count, grid, io_mode);
		}
	}
	else {
		/* get/put whole grid contiguous */
		count[yx_dim[0]] = height_t;
		count[yx_dim[1]] = width_t;
		if (stride)
			status = gmtnc_varm_grdfloat (HH->ncid, HH->z_id, start, count, onestride, imap, grid, io_mode);
		else
			status = gmtnc_vara_grdfloat (HH->ncid, HH->z_id, start, count, grid, io_mode);
	}
	return status;
}

GMT_LOCAL void gmtnc_get_units (struct GMT_CTRL *GMT, int ncid, int varid, char *name_units) {
	/* Get attributes long_name and units for given variable ID
	 * and assign variable name if attributes are not available.
	 * ncid, varid		: as in nc_get_att_text
	 * nameunit		: long_name and units in form "long_name [units]"
	 */
	char name[GMT_GRID_UNIT_LEN80], units[GMT_GRID_UNIT_LEN80];
	if (gmtlib_nc_get_att_text (GMT, ncid, varid, "long_name", name, GMT_GRID_UNIT_LEN80))
		nc_inq_varname (ncid, varid, name);
	if (!gmtlib_nc_get_att_text (GMT, ncid, varid, "units", units, GMT_GRID_UNIT_LEN80) && units[0])
		snprintf (name_units, GMT_GRID_UNIT_LEN80, "%s [%s]", name, units);
	else
		strncpy (name_units, name, GMT_GRID_UNIT_LEN80-1);
}

GMT_LOCAL void gmtnc_put_units (int ncid, int varid, char *name_units) {
	/* Put attributes long_name and units for given variable ID based on
	 * string name_unit in the form "long_name [units]".
	 * ncid, varid		: as is nc_put_att_text
	 * name_units		: string in form "long_name [units]"
	 */
	bool copy = false;
	int i = 0, j = 0;
	char name[GMT_GRID_UNIT_LEN80], units[GMT_GRID_UNIT_LEN80];

	strncpy (name, name_units, GMT_GRID_UNIT_LEN80-1);
	units[0] = '\0';
	for (i = 0; i < GMT_GRID_UNIT_LEN80 && name[i]; i++) {
		if (name[i] == ']') copy = false, units[j] = '\0';
		if (copy) units[j++] = name[i];
		if (name[i] == '[') {
			name[i] = '\0';
			if (i > 0 && name[i-1] == ' ') name[i-1] = '\0';
			copy = true;
		}
	}
	if (name[0]) nc_put_att_text (ncid, varid, "long_name", strlen(name), name);
	if (units[0]) nc_put_att_text (ncid, varid, "units", strlen(units), units);
	if (strstr(units, "degrees_east")) {
		nc_put_att_text (ncid, varid, "standard_name", 9, "longitude");
		nc_put_att_text (ncid, varid, "axis", 1, "X");
	}
	else if (strstr(units, "degrees_north")) {
		nc_put_att_text (ncid, varid, "standard_name", 8, "latitude");
		nc_put_att_text (ncid, varid, "axis", 1, "Y");
	}
}

GMT_LOCAL int gmtnc_check_step (struct GMT_CTRL *GMT, uint32_t n, double *x, char *varname, char *file, bool save_xy_array) {
	/* Check if all steps in range are the same (within 0.1%). Returns 0 if OK and 1 if variable spacing */
	double step, step_min, step_max;
	uint32_t i;
	unsigned int wlevel = (save_xy_array) ? GMT_MSG_INFORMATION : GMT_MSG_WARNING;
	if (n < 2) return 0;
	step_min = step_max = x[1]-x[0];
	for (i = 2; i < n; i++) {
		step = x[i]-x[i-1];
		if (step < step_min) step_min = step;
		if (step > step_max) step_max = step;
	}
	if (fabs (step_min-step_max)/(fabs (step_min)*0.5 + fabs (step_max)*0.5) > 0.001) {
		GMT_Report (GMT->parent, wlevel,
			"The step size of coordinate (%s) in grid %s is not constant.\n", varname, file);
		GMT_Report (GMT->parent, wlevel,
			"GMT will use a constant step size of %g; the original ranges from %g to %g.\n",
			(x[n-1]-x[0])/(n-1), step_min, step_max);
		return 1;
	}
	return 0;
}

GMT_LOCAL void gmtnc_set_optimal_chunksize (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
	/* For optimal performance, set the number of elements in a given chunk
	 * dimension (n) to be the ceiling of the number of elements in that
	 * dimension of the array variable (d) divided by a natural number N>1.
	 * That is, set n = ceil (d / N).  Using a chunk size slightly larger than
	 * this value is also acceptable.  For example: 129 = ceil (257 / 2).
	 * Do NOT set n = floor (d / N), for example 128 = floor (257 / 2). */

	double chunksize[2] = {128, 128};                            /* default min chunksize */
	const size_t min_chunk_pixels = (size_t)(chunksize[0] * chunksize[1]); /* min pixel count per chunk */

	if (GMT->current.setting.io_nc4_chunksize[0] == k_netcdf_io_classic)
		/* no chunking with classic model */
		return;

	if (GMT->current.setting.io_nc4_chunksize[0] != k_netcdf_io_chunked_auto && (
			header->n_rows >= GMT->current.setting.io_nc4_chunksize[0] ||
			header->n_columns >= GMT->current.setting.io_nc4_chunksize[1])) {
		/* if chunk size is smaller than grid size */
		return;
	}

	/* here, chunk size is either k_netcdf_io_chunked_auto or the chunk size is
	 * larger than grid size */

	if ( header->nm < min_chunk_pixels ) {
		/* the grid dimension is too small for chunking to make sense. switch to
		 * classic model */
		GMT->current.setting.io_nc4_chunksize[0] = k_netcdf_io_classic;
		return;
	}

	/* adjust default chunk sizes for grids that have more than min_chunk_pixels
	 * cells but less than chunksize (default 128) cells in one dimension */
	if (header->n_rows < chunksize[0]) {
		chunksize[0] = header->n_rows;
		chunksize[1] = floor (min_chunk_pixels / chunksize[0]);
	}
	else if (header->n_columns < chunksize[1]) {
		chunksize[1] = header->n_columns;
		chunksize[0] = floor (min_chunk_pixels / chunksize[1]);
	}

	/* determine optimal chunk size in the range [chunksize,2*chunksize) */
	GMT->current.setting.io_nc4_chunksize[0] = (size_t) ceil (header->n_rows / floor (header->n_rows / chunksize[0]));
	GMT->current.setting.io_nc4_chunksize[1] = (size_t) ceil (header->n_columns / floor (header->n_columns / chunksize[1]));
}

GMT_LOCAL bool gmtnc_not_obviously_global (double *we) {
	/* If range is not 360 and boundaries are not 0, 180, 360 then we pass true */
	if (!gmt_M_360_range (we[0], we[1])) return true;
	if (!(doubleAlmostEqualZero (we[0], 0.0) || doubleAlmostEqual (we[0], -180.0))) return true;
	return false;
}

GMT_LOCAL bool gmtnc_not_obviously_polar (double *se) {
	/* If range is not 180 and boundaries are not -90/90 then we pass true */
	if (!gmt_M_180_range (se[0], se[1])) return true;
	if (!(doubleAlmostEqual (se[0], -90.0) && doubleAlmostEqual (se[1], 90.0))) return true;
	return false;
}

GMT_LOCAL int gmtnc_put_xy_vectors (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
	/* Write the x and y vectors to the netCDF file */
	int err = GMT_NOERROR;
	/* Create enough memory to store the x- and y-coordinate values */
	double *xy = gmt_M_memory (GMT, NULL, MAX (header->n_columns,header->n_rows), double);
	unsigned int col, row;
	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);

	for (col = 0; col < header->n_columns; col++) xy[col] = gmt_M_grd_col_to_x (GMT, col, header);
	gmt_M_err_trap (nc_put_var_double (HH->ncid, HH->xyz_id[GMT_X], xy));

	/* Depending on row_order, write y-coordinate array bottom-to-top or top-to-bottom */
	if (HH->row_order == k_nc_start_south) {
		for (row = 0; row < header->n_rows; row++)
			xy[row] = (double) gmt_M_col_to_x (GMT, row, header->wesn[YLO], header->wesn[YHI], header->inc[GMT_Y], 0.5 * header->registration, header->n_rows);
	}
	else {
		for (row = 0; row < header->n_rows; row++)
			xy[row] = (double) gmt_M_row_to_y (GMT, row, header->wesn[YLO], header->wesn[YHI], header->inc[GMT_Y], 0.5 * header->registration, header->n_rows);
	}
	gmt_M_err_trap (nc_put_var_double (HH->ncid, HH->xyz_id[GMT_Y], xy));
	gmt_M_free (GMT, xy);
	return GMT_NOERROR;
}

GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, char job, bool cube, uint64_t n_layers) {
	/* Used when reading and writing 2-D grids and when writing 3-D data cubes */
	int j, err, has_vector, has_range, registration, var_spacing = 0;
	int old_fill_mode, status;
	unsigned int d_off = (cube) ? 1 : 0;
	size_t len;
	double dummy[2], min, max;
	char dimname[GMT_GRID_UNIT_LEN80], coord[GMT_LEN8];
	nc_type z_type;
	bool save_xy_array = !strncmp (GMT->init.module_name, "grd2xyz", 7U);
	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);

	/* Dimension ids, variable ids, etc.. */
	int i, ncid, z_id = GMT_NOTSET, ids[5] = {-1,-1,-1,-1,-1}, gm_id = GMT_NOTSET, dims[5], nvars, ndims = 0;
	size_t lens[5], item[2];

	/* If not yet determined, attempt to get the layer IDs from the variable name */
	int n_t_index_by_value = 0;
	double t_value[3];

	if (HH->varname[0]) {
		char *c = strpbrk (HH->varname, "(["); /* find first occurrence of ( or [ */
		if (c != NULL && *c == '(') {
			n_t_index_by_value = sscanf (c+1, "%lf,%lf,%lf", &t_value[0], &t_value[1], &t_value[2]);
			*c = '\0';
			/* t_index will be determined later from t_value when the nc-file is opened */
		}
		else if (c != NULL && *c == '[') {
			sscanf (c+1, "%" SCNuS ",%" SCNuS ",%" SCNuS "", &HH->t_index[0], &HH->t_index[1], &HH->t_index[2]);
			*c = '\0';
		}
	}

	/* Open NetCDF file */
	if (!strcmp (HH->name,"=")) return (GMT_GRDIO_NC_NO_PIPE);
	switch (job) {
		case 'r':
			gmt_M_err_trap (gmt_nc_open (GMT, HH->name, NC_NOWRITE, &ncid));
			break;
		case 'u':
			gmt_M_err_trap (gmt_nc_open (GMT, HH->name, NC_WRITE, &ncid));
			gmt_M_err_trap (nc_set_fill (ncid, NC_NOFILL, &old_fill_mode));
			break;
		default:
			/* create new nc-file */
			gmtnc_set_optimal_chunksize (GMT, header);
			if (GMT->current.setting.io_nc4_chunksize[0] != k_netcdf_io_classic) {
				/* create chunked nc4 file */
				gmt_M_err_trap (gmt_nc_create (GMT, HH->name, NC_NETCDF4, &ncid));
				HH->is_netcdf4 = true;
			}
			else {
				/* create nc using classic data model */
				gmt_M_err_trap (gmt_nc_create (GMT, HH->name, NC_CLOBBER, &ncid));
				HH->is_netcdf4 = false;
			}
			gmt_M_err_trap (nc_set_fill (ncid, NC_NOFILL, &old_fill_mode));
			break;
	}

	/* Retrieve or define dimensions and variables */

	if (job == 'r' || job == 'u') {
		int kind;
		/* determine netCDF data model */
		gmt_M_err_trap (nc_inq_format(ncid, &kind));
		HH->is_netcdf4 = (kind == NC_FORMAT_NETCDF4 || kind == NC_FORMAT_NETCDF4_CLASSIC);

		/* First see if this is an old NetCDF formatted file */
		if (!nc_inq_dimid (ncid, "xysize", &i)) return (gmt_cdf_grd_info (GMT, ncid, header, job));

		/* Find first 2-dimensional (z) variable or specified variable */
		if (!HH->varname[0]) {	/* No specific named variable was requested; search for suitable matrix */
			int ID = GMT_NOTSET, dim = 0;
			gmt_M_err_trap (nc_inq_nvars (ncid, &nvars));
			i = 0;
			while (i < nvars && z_id < 0) {	/* Look for first 2D grid, with fallback to first higher-dimension (3-4D) grid if 2D not found */
				gmt_M_err_trap (nc_inq_varndims (ncid, i, &ndims));
				if (ndims == 2)	/* Found the first 2-D grid */
					z_id = i;
				else if (ID == GMT_NOTSET && ndims > 2 && ndims < 5) {	/* Also look for higher-dim grid in case no 2-D */
					ID = i;
					dim = ndims;
				}
				i++;
			}
			if (z_id < 0) {	/* No 2-D grid found, check if we found a higher dimension cube */
				if (ID == GMT_NOTSET) return (GMT_GRDIO_NO_2DVAR);	/* No we didn't */
				z_id = ID;		/* Pick the higher dimensioned cube instead, get its name, and warn */
				nc_inq_varname (ncid, z_id, HH->varname);
				GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No 2-D array in file %s.  Selecting first 2-D slice in the %d-D array %s\n", HH->name, dim, HH->varname);
			}
		}
		else if (nc_inq_varid (ncid, HH->varname, &z_id) == NC_NOERR) {	/* Gave a named variable that exist, is it 2-4 D? */
			gmt_M_err_trap (nc_inq_varndims (ncid, z_id, &ndims));
			if (ndims < 2 || ndims > 5) return (GMT_GRDIO_BAD_DIM);
		}
		else	/* No such named variable in the grid */
			return (GMT_GRDIO_NO_VAR);

		if (z_id < 0) return (GMT_GRDIO_NO_2DVAR);

		/* Get the z data type and determine its dimensions */
		gmt_M_err_trap (nc_inq_vartype (ncid, z_id, &z_type));
		gmt_M_err_trap (nc_inq_vardimid (ncid, z_id, dims));
		switch (z_type) {
			case NC_UBYTE:  header->type = GMT_GRID_IS_NB; HH->orig_datatype = GMT_UCHAR; break;
			case NC_BYTE:   header->type = GMT_GRID_IS_NB; HH->orig_datatype = GMT_CHAR; break;
			case NC_SHORT:  header->type = GMT_GRID_IS_NS; HH->orig_datatype = GMT_SHORT; break;
			case NC_INT:    header->type = GMT_GRID_IS_NI; HH->orig_datatype = GMT_INT; break;
			case NC_FLOAT:  header->type = GMT_GRID_IS_NF; HH->orig_datatype = GMT_FLOAT; break;
			case NC_DOUBLE: header->type = GMT_GRID_IS_ND; HH->orig_datatype = GMT_DOUBLE; break;
			default:        header->type = k_grd_unknown_fmt; break;
		}

		/* Get the ids of the x and y (and depth and time) coordinate variables */
		for (i = 0; i < ndims; i++) {
			gmt_M_err_trap (nc_inq_dim (ncid, dims[i], dimname, &lens[i]));
			//if (nc_inq_varid (ncid, dimname, &ids[i])) return (GMT_GRDIO_NC_NOT_COARDS);
			if ((status = nc_inq_varid (ncid, dimname, &ids[i])) != NC_NOERR)
				GMT_Report (GMT->parent, GMT_MSG_WARNING, "\"%s\", %s\n\tIf something bad happens later, try importing via GDAL.\n",
				                                         dimname, nc_strerror(status));
		}
		HH->xy_dim[0] = ndims-1;
		HH->xy_dim[1] = ndims-2;
		HH->xyz_id[GMT_X] = ids[1];
		HH->xyz_id[GMT_Y] = ids[0];

		/* Check if LatLon variable exists, then we may need to flip x and y */
		if (nc_inq_varid (ncid, "LatLon", &i) == NC_NOERR) nc_get_var_int (ncid, i, HH->xy_dim);
		header->n_columns = (uint32_t) lens[HH->xy_dim[0]];
		header->n_rows = (uint32_t) lens[HH->xy_dim[1]];

		/* Check if the spatial_ref attribute exists, first via the z_id variable, and if not then via a global attribute */
		if (nc_inq_attlen (ncid, z_id, "grid_mapping", &len) == NC_NOERR) {	/* The data layer has an attribute with the name of the variable that has the spatial_ref */
			char *v_name = gmt_M_memory (GMT, NULL, len+1, char);           /* Allocate the needed space for the variable name */
			if (nc_get_att_text (ncid, z_id, "grid_mapping", v_name) == NC_NOERR) {	/* Get the name of the variable */
				if (nc_inq_varid (ncid, v_name, &i) == NC_NOERR)	/* Got the id of this variable */
					gm_id = i;
			}
			gmt_M_free (GMT, v_name);
		}
		else if (nc_inq_varid (ncid, "grid_mapping", &i) == NC_NOERR)	/* A global variable has it */
			gm_id = i;
	} /* if (job == 'r' || job == 'u') */
	else {
		/* Define dimensions of z variable */
		ndims = (cube) ? 3 : 2;
		HH->xy_dim[0] = 1 + d_off;
		HH->xy_dim[1] = 0 + d_off;

		strcpy (coord, (gmt_M_x_is_lon (GMT, GMT_OUT)) ? "lon" : (gmt_M_type (GMT, GMT_OUT, GMT_X) & GMT_IS_RATIME) ? "time" : "x");
		gmt_M_err_trap (nc_def_dim (ncid, coord, (size_t) header->n_columns, &dims[1+d_off]));
		gmt_M_err_trap (nc_def_var (ncid, coord, NC_DOUBLE, 1, &dims[1+d_off], &ids[1+d_off]));
		HH->xyz_id[GMT_X] = ids[1+d_off];

		strcpy (coord, (gmt_M_y_is_lat (GMT, GMT_OUT)) ? "lat" : (gmt_M_type (GMT, GMT_OUT, GMT_Y) & GMT_IS_RATIME) ? "time" : "y");
		gmt_M_err_trap (nc_def_dim (ncid, coord, (size_t) header->n_rows, &dims[0+d_off]));
		gmt_M_err_trap (nc_def_var (ncid, coord, NC_DOUBLE, 1, &dims[0+d_off], &ids[0+d_off]));
		HH->xyz_id[GMT_Y] = ids[0+d_off];

		if (cube) {	/* Allow for cube expansion by setting layer dimension to 0 (NC_UNLIMITED) */
			strcpy (coord, (gmt_M_type (GMT, GMT_OUT, GMT_Z) & GMT_IS_RATIME) ? "time" : "z");
			gmt_M_err_trap (nc_def_dim (ncid, coord, (size_t) n_layers, &dims[0]));
			gmt_M_err_trap (nc_def_var (ncid, coord, NC_DOUBLE, 1, &dims[0], &HH->xyz_id[GMT_Z]));
		}

		switch (header->type) {
			case GMT_GRID_IS_NB: z_type = NC_BYTE; break;
			case GMT_GRID_IS_NS: z_type = NC_SHORT; break;
			case GMT_GRID_IS_NI: z_type = NC_INT; break;
			case GMT_GRID_IS_NF: z_type = NC_FLOAT; break;
			case GMT_GRID_IS_ND: z_type = NC_DOUBLE; break;
			default: z_type = NC_NAT;
		}

		/* Variable name is given, or defaults to "z" (or "cube" under cube) */
		if (!HH->varname[0])
			strcpy (HH->varname, (cube) ? "cube" : "z");
		gmt_M_err_trap (nc_def_var (ncid, HH->varname, z_type, (cube) ? 3 : 2, dims, &z_id));

		/* set deflation and chunking */
		if (GMT->current.setting.io_nc4_chunksize[0] != k_netcdf_io_classic) {
			/* set chunk size */
			size_t nc4_chunksize[3] = {0, 0, 0};
			gmt_M_memcpy (&nc4_chunksize[cube], GMT->current.setting.io_nc4_chunksize, 2U, size_t);
			gmt_M_err_trap (nc_def_var_chunking (ncid, z_id, NC_CHUNKED, nc4_chunksize));
			/* set deflation level and shuffle for x, y[, z], and z[|w] variable */
			if (GMT->current.setting.io_nc4_deflation_level) {
				gmt_M_err_trap (nc_def_var_deflate (ncid, HH->xyz_id[GMT_X], true, true, GMT->current.setting.io_nc4_deflation_level));
				gmt_M_err_trap (nc_def_var_deflate (ncid, HH->xyz_id[GMT_Y], true, true, GMT->current.setting.io_nc4_deflation_level));
				if (cube && n_layers == 0) gmt_M_err_trap (nc_def_var_deflate (ncid, HH->xyz_id[GMT_Z], true, true, GMT->current.setting.io_nc4_deflation_level));
				gmt_M_err_trap (nc_def_var_deflate (ncid, z_id, true, true, GMT->current.setting.io_nc4_deflation_level));
			}
		} /* GMT->current.setting.io_nc4_chunksize[0] != k_netcdf_io_classic */
	} /* if (job == 'r' || job == 'u') */
	HH->z_id = z_id;
	HH->ncid = ncid;

	/* Query or assign attributes */

	if (job == 'u') gmt_M_err_trap (nc_redef (ncid));

	if (job == 'r') {
		bool set_reg = true;
		double dx = 0, dy = 0, threshold = 0.0;
		/* Create enough memory to store the x- and y-coordinate values */
		double *xy = gmt_M_memory (GMT, NULL, MAX (header->n_columns,header->n_rows), double);
		/* Get global information */
		if (gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "title", header->title, GMT_GRID_TITLE_LEN80))
			gmtlib_nc_get_att_text (GMT, ncid, z_id, "long_name", header->title, GMT_GRID_TITLE_LEN80);
		if (gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "history", header->command, GMT_GRID_COMMAND_LEN320))
			gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "source", header->command, GMT_GRID_COMMAND_LEN320);
		gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "description", header->remark, GMT_GRID_REMARK_LEN160);
		header->registration = GMT_GRID_NODE_REG;
		if (!nc_get_att_int (ncid, NC_GLOBAL, "node_offset", &i)) {	/* GMT wrote the registration in the grid */
			header->registration = i;
			set_reg = false;	/* Do not update it below since we know the registration */
		}

		if (gm_id != GMT_NOTSET && nc_inq_attlen (ncid, gm_id, "spatial_ref", &len) == NC_NOERR) {	/* Good to go */
			char *attrib = NULL;
			gmt_M_str_free (header->ProjRefWKT);	/* Make sure we didn't have a previously allocated one */
			attrib = gmt_M_memory (GMT, NULL, len+1, char);		/* and allocate the needed space */
			gmt_M_err_trap (nc_get_att_text (ncid, gm_id, "spatial_ref", attrib));
			header->ProjRefWKT = strdup (attrib);	/* Turn it into a strdup allocation to be compatible with other instances elsewhere */
			gmt_M_free (GMT, attrib);
		}

		/* Explanation for the logic below: Not all netCDF grids are proper COARDS grids and hence we sometime must guess
		 * regarding the settings.  The x and y coordinates may be written as arrays, which reflect the positions of the
		 * nodes.  There may also be the attributes actual_range which specifies the x range.  These will differ if the
		 * grid is pixel registered, hence when both are present we use this difference to detect a pixel grid. However,
		 * some external tools such as xarray may slice a grid but not update the attributes.  In this case the actual_range
		 * may have an initial range that is no longer the case.  We have added a check if these differ by more than a
		 * half grid increment.  If not then we can trust it.  If actual_range is missing then we have to guess the registration
		 * which we do by checking if the start coordinate is an integer multiple of the increment.  If not, we guess pixel registration
		 * but cannot know if this is the case unless the adjusted coordinates in x has a range of 360 and in y a range of 180.
		 * Finally, if there is no array just the actual_range, then we cannot tell the registration from the range but try
		 * and leave it as gridline registration. */

		/* Get information about x variable */
		gmtnc_get_units (GMT, ncid, ids[HH->xy_dim[0]], header->x_units);
		/* Set default range to number of nodes in case nothing is found further down */
		dummy[0] = 0.0, dummy[1] = (double) header->n_columns-1; /* Default */
		registration = GMT_GRID_NODE_REG;

		/* Look for the x-coordinate vector */
		if ((has_vector = !nc_get_var_double (ncid, ids[HH->xy_dim[0]], xy)) && header->n_columns > 1) {
			var_spacing = gmtnc_check_step (GMT, header->n_columns, xy, header->x_units, HH->name, save_xy_array);
			if (save_xy_array && var_spacing) {
				if (GMT->current.io.nc_xarray) gmt_M_free (GMT, GMT->current.io.nc_xarray);
				GMT->current.io.nc_xarray = gmt_M_memory (GMT, NULL, header->n_columns, double);
				gmt_M_memcpy (GMT->current.io.nc_xarray, xy, header->n_columns, double);
				HH->var_spacing[GMT_X] = var_spacing;
			}
			dx = fabs (xy[1] - xy[0]);	/* Grid spacing in x */
		}
		if (has_vector && header->n_columns == 1) has_vector = false;	/* One coordinate does not a vector make */

		/* Look for the x-coordinate range attributes */
		has_range = (!nc_get_att_double (ncid, ids[HH->xy_dim[0]], "actual_range", dummy) ||
			!nc_get_att_double (ncid, ids[HH->xy_dim[0]], "valid_range", dummy) ||
			!(nc_get_att_double (ncid, ids[HH->xy_dim[0]], "valid_min", &dummy[0]) +
			nc_get_att_double (ncid, ids[HH->xy_dim[0]], "valid_max", &dummy[1])));

		if (has_vector && has_range) {	/* Has both so we can do a basic sanity check */
			threshold = (0.5+GMT_CONV5_LIMIT) * dx;
			min = xy[0];	max = xy[header->n_columns-1];
			if (min > max) gmt_M_double_swap (min, max);	/* Got it backwards in the array */
			if (fabs (dummy[0] - min) > threshold || fabs (dummy[1] - max) > threshold) {
				if (gmtnc_not_obviously_global (dummy)) {
					GMT_Report (GMT->parent, GMT_MSG_WARNING, "The x-coordinates and range attribute are in conflict; must rely on coordinates only\n");
					dummy[0] = xy[0], dummy[1] = xy[header->n_columns-1];
					has_range = false;	/* Since useless information */
					/* For registration, we have to assume that the actual range is an integer multiple of increment.
					 * If so, then if the coordinates are off by 0.5*dx then we assume we have pixel registration */
					threshold = (0.5-GMT_CONV5_LIMIT) * dx;
					if (set_reg && fabs (fmod (dummy[0], dx)) > threshold) {	/* Pixel registration */
						registration = header->registration = GMT_GRID_PIXEL_REG;
						dummy[0] -= 0.5 * dx;	dummy[1] += 0.5 * dx;
						if (gmt_M_360_range (dummy[0], dummy[1]))
							GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid x-coordinates after pixel registration adjustment have exactly 360 range\n");
					}
				}
				else {	/* Got clean global longitude range, stick with that information and ignore xy */
					GMT_Report (GMT->parent, GMT_MSG_WARNING, "The x-coordinates and range attribute are in conflict but range is exactly 360; we rely on this range\n");
					if (set_reg && (header->n_columns%2) == 0) {	/* Pixel registration */
						registration = header->registration = GMT_GRID_PIXEL_REG;
						GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Global longitudes, guessing registration to be %s since nx is odd\n", regtype[header->registration]);
					}
				}
			}
			else {	/* Data seems OK; determine registration */
				dummy[0] = xy[0], dummy[1] = xy[header->n_columns-1];
				if (set_reg && (fabs(dummy[1] - dummy[0]) / fabs(xy[header->n_columns-1] - xy[0]) - 1.0 > 0.5 / (header->n_columns - 1)))
					registration = header->registration = GMT_GRID_PIXEL_REG;
			}
		}
		else if (has_vector) {	/* No attribute for range, use coordinates */
			dummy[0] = xy[0], dummy[1] = xy[header->n_columns-1];
			threshold = (0.5-GMT_CONV5_LIMIT) * dx;
			if (set_reg && fabs (fmod (dummy[0], dx)) > threshold) {	/* Most likely pixel registration since off by dx/2 */
				registration = header->registration = GMT_GRID_PIXEL_REG;
				dummy[0] -= 0.5 * dx;	dummy[1] += 0.5 * dx;
				if (gmt_M_360_range (dummy[0], dummy[1]))
					GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid x-coordinates after pixel registration adjustment have exactly 360 range\n");
			}
			GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No range attribute, guessing registration to be %s\n", regtype[header->registration]);
		}
		else {	/* Only has the valid_range settings.  If no registration set, and no dx available, guess based on nx */
			if (set_reg && (header->n_columns%2) == 0) {	/* Pixel registration */
				registration = header->registration = GMT_GRID_PIXEL_REG;
				GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No x-coordinates, guessing registration to be %s since nx is odd\n", regtype[header->registration]);
			}
		}

		/* Determine grid step */
		header->inc[GMT_X] = gmt_M_get_inc (GMT, dummy[0], dummy[1], header->n_columns, registration);
		if (gmt_M_is_dnan(header->inc[GMT_X]) || gmt_M_is_zero (header->inc[GMT_X])) header->inc[GMT_X] = 1.0;
		if (header->n_columns == 1) registration = GMT_GRID_PIXEL_REG;	/* The only way to have a grid like that */

#ifdef NC4_DEBUG
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "x registration: %u\n", header->registration);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "x dummy: %g %g\n", dummy[0], dummy[1]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "x[0] x[nx-1]: %g %g\n", xy[0], xy[header->n_columns-1]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "xinc: %g %g\n", header->inc[GMT_X]);
#endif

		/* Extend x boundaries by half if we found pixel registration */
		if (registration == GMT_GRID_NODE_REG && header->registration == GMT_GRID_PIXEL_REG)
			header->wesn[XLO] = dummy[0] - header->inc[GMT_X] / 2.0,
			header->wesn[XHI] = dummy[1] + header->inc[GMT_X] / 2.0;
		else	/* Use as is */
			header->wesn[XLO] = dummy[0], header->wesn[XHI] = dummy[1];

		/* Get information about y variable */
		gmtnc_get_units (GMT, ncid, ids[HH->xy_dim[1]], header->y_units);
		/* Set default range to number of nodes in case nothing is found further down */
		dummy[0] = 0.0, dummy[1] = (double) header->n_rows-1; /* Default */
		registration = GMT_GRID_NODE_REG;

		/* Read the y-coordinate vector (if available), otherwise just look for range attributes */
		if ((has_vector = !nc_get_var_double (ncid, ids[HH->xy_dim[1]], xy)) && header->n_rows > 1) {
			var_spacing = gmtnc_check_step (GMT, header->n_rows, xy, header->y_units, HH->name, save_xy_array);
			if (save_xy_array && var_spacing) {
				if (GMT->current.io.nc_yarray) gmt_M_free (GMT, GMT->current.io.nc_yarray);
				GMT->current.io.nc_yarray = gmt_M_memory (GMT, NULL, header->n_rows, double);
				gmt_M_memcpy (GMT->current.io.nc_yarray, xy, header->n_rows, double);
				HH->var_spacing[GMT_Y] = var_spacing;
				/* Flip-ud y-array since row = 0 is the last value */
				gmt_grd_flip_vertical (GMT->current.io.nc_yarray, 1, header->n_rows, 0, sizeof(double));
			}
			dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1];
			dy = fabs (xy[1] - xy[0]);	/* Grid spacing in y */
		}
		if (has_vector && header->n_rows == 1) has_vector = false;	/* One coordinate does not a vector make */
		has_range = (!nc_get_att_double (ncid, ids[HH->xy_dim[1]], "actual_range", dummy) ||
			!nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_range", dummy) ||
			!(nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_min", &dummy[0]) +
			nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_max", &dummy[1])));

		if (has_vector && has_range) {	/* Has both so we can do a basic sanity check */
			threshold = (0.5+GMT_CONV5_LIMIT) * dy;
			min = xy[0];	max = xy[header->n_rows-1];
			if (min > max) gmt_M_double_swap (min, max);	/* Got it backwards in the array */
			if (fabs (dummy[0] - min) > threshold || fabs (dummy[1] - max) > threshold) {
				if (gmtnc_not_obviously_polar (dummy)) {
					GMT_Report (GMT->parent, GMT_MSG_WARNING, "The y-coordinates and range attribute are in conflict; must rely on coordinates only\n");
					dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1];
					has_range = false;	/* Since useless information */
					/* Registration was set using x values, so here we just check that we get the same result.
					 * If the coordinates are off by 0.5*dy then we assume we have pixel registration */
					threshold = (0.5-GMT_CONV5_LIMIT) * dy;
					if (fabs (fmod (dummy[0], dy)) > threshold) {	/* Pixel registration? */
						if (header->registration == GMT_GRID_NODE_REG)	/* No, somehow messed up now */
							GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]);
						else {	/* Pixel registration confirmed */
							dummy[0] -= 0.5 * dy;	dummy[1] += 0.5 * dy;
							registration = GMT_GRID_PIXEL_REG;
							if (gmt_M_180_range (dummy[0], dummy[1]))
								GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid y-coordinates after pixel registration adjustment have exactly 180 range\n");
						}
					}
				}
				else {
					GMT_Report (GMT->parent, GMT_MSG_WARNING, "The y-coordinates and range attribute are in conflict but range is exactly 180; we rely on this range\n");
					if ((header->n_rows%2) == 1 && header->registration == GMT_GRID_NODE_REG)	/* Pixel registration? */
						GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]);
					else
						registration = header->registration;
				}
			}
			else {	/* Data seems OK; determine registration and set dummy from data coordinates */
				dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1];
				if ((fabs(dummy[1] - dummy[0]) / fabs(xy[header->n_rows-1] - xy[0]) - 1.0 > 0.5 / (header->n_rows - 1)) && header->registration == GMT_GRID_NODE_REG)
					GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]);
			}
		}
		else if (has_vector) {	/* No attribute for range, use coordinates */
			threshold = (0.5-GMT_CONV5_LIMIT) * dy;
			dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1];
			if (fabs (fmod (dummy[0], dy)) > threshold) {	/* Most likely pixel registration since off by dy/2 */
				if (header->registration == GMT_GRID_NODE_REG)	/* No, somehow messed up now */
					GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]);
				else {	/* Pixel registration confirmed */
					if (dummy[0] > dummy[1]) {		/* Check for reverse order of y-coordinate */
						gmt_M_double_swap (dummy[0], dummy[1]);
						HH->row_order = k_nc_start_north;
					}
					dummy[0] -= 0.5 * dy;	dummy[1] += 0.5 * dy;
					registration = GMT_GRID_PIXEL_REG;
					if (gmt_M_180_range (dummy[0], dummy[1]))
						GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Grid y-coordinates after pixel registration adjustment have exactly 180 range\n");
				}
			}
		}
		else {	/* Only has the valid_range settings.  If no registration set, and no dy available, guess based on ny */
			if ((header->n_rows%2) == 1 && header->registration == GMT_GRID_NODE_REG)	/* Pixel registration? */
				GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]);
		}

		/* Check for reverse order of y-coordinate */
		if (dummy[0] > dummy[1]) {
			HH->row_order = k_nc_start_north;
			gmt_M_double_swap (dummy[0], dummy[1]);
		}
		else if (has_vector && (xy[0] > xy[header->n_rows-1]) && (dummy[0] > dummy[1])) 	/* Here the lat vector is top-down but range is bottum-up */
			HH->row_order = k_nc_start_north;
		else if (HH->row_order != k_nc_start_north)		/* If not already set to north in an above test */
			HH->row_order = k_nc_start_south;

		/* Determine grid step */
		header->inc[GMT_Y] = gmt_M_get_inc (GMT, dummy[0], dummy[1], header->n_rows, registration);
		if (gmt_M_is_dnan(header->inc[GMT_Y]) || gmt_M_is_zero (header->inc[GMT_Y])) header->inc[GMT_Y] = 1.0;
		if (header->n_rows == 1) registration = GMT_GRID_PIXEL_REG;	/* The only way to have a grid like that */

#ifdef NC4_DEBUG
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "y registration: %u\n", header->registration);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "y dummy: %g %g\n", dummy[0], dummy[1]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "y[0] y[ny-1]: %g %g\n", xy[0], xy[header->n_rows-1]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "xinc: %g %g\n", header->inc[GMT_Y]);
#endif

		/* Extend y boundaries by half if we found pixel registration */
		if (registration == GMT_GRID_NODE_REG && header->registration == GMT_GRID_PIXEL_REG)
			header->wesn[YLO] = dummy[0] - header->inc[GMT_Y] / 2.0,
			header->wesn[YHI] = dummy[1] + header->inc[GMT_Y] / 2.0;
		else	/* Use as is */
			header->wesn[YLO] = dummy[0], header->wesn[YHI] = dummy[1];

		gmt_M_free (GMT, xy);	/* Done with the array */

		/* Get information about z variable */
		gmtnc_get_units (GMT, ncid, z_id, header->z_units);
		if (nc_get_att_double (ncid, z_id, "scale_factor", &header->z_scale_factor)) header->z_scale_factor = 1.0;
		if (nc_get_att_double (ncid, z_id, "add_offset", &header->z_add_offset)) header->z_add_offset = 0.0;
#ifdef DOUBLE_PRECISION_GRID
		if (nc_get_att_double (ncid, z_id, "_FillValue", &header->nan_value))
		    nc_get_att_double (ncid, z_id, "missing_value", &header->nan_value);
#else
		if (nc_get_att_float (ncid, z_id, "_FillValue", &header->nan_value))
		    nc_get_att_float (ncid, z_id, "missing_value", &header->nan_value);
#endif
		if (!nc_get_att_double (ncid, z_id, "actual_range", dummy)) {
			/* z-limits need to be converted from actual to internal grid units. */
			header->z_min = (dummy[0] - header->z_add_offset) / header->z_scale_factor;
			header->z_max = (dummy[1] - header->z_add_offset) / header->z_scale_factor;
		}
		else if (!nc_get_att_double (ncid, z_id, "valid_range", dummy) ||
			!(nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_min", &dummy[0]) +
			nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_max", &dummy[1]))) {
			/* Valid range is already in packed units, so do not convert */
			header->z_min = dummy[0], header->z_max = dummy[1];
		}
		if (gmt_M_is_dnan (header->z_min) && gmt_M_is_dnan (header->z_max)) {
			GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "netCDF grid %s information has zmin = zmax = NaN. Reset to 0/0.\n", HH->name);
			header->z_min = header->z_max = 0.0;
		}
		{	/* Get deflation and chunking info */
			int storage_mode, shuffle = 0, deflate = 0, deflate_level = 0;
			size_t chunksize[5]; /* chunksize of z */
			gmt_M_err_trap (nc_inq_var_chunking (ncid, z_id, &storage_mode, chunksize));
			if (storage_mode == NC_CHUNKED) {
				HH->z_chunksize[0] = chunksize[dims[0]]; /* chunk size in vertical dimension */
				HH->z_chunksize[1] = chunksize[dims[1]]; /* chunk size in horizontal dimension */
			}
			else { /* NC_CONTIGUOUS */
				HH->z_chunksize[0] = HH->z_chunksize[1] = 0;
			}
			if (HH->is_netcdf4) gmt_M_err_trap (nc_inq_var_deflate (ncid, z_id, &shuffle, &deflate, &deflate_level));
			HH->z_shuffle = shuffle ? true : false; /* if shuffle filter is turned on */
			HH->z_deflate_level = deflate ? deflate_level : 0; /* if deflate filter is in use */
		}

		/* Determine t_index from t_value */
		item[0] = 0;
		for (i = 0; i < n_t_index_by_value; i++) {
			item[1] = lens[i]-1;
			if (nc_get_att_double (ncid, ids[i], "actual_range", dummy)) {
				gmt_M_err_trap (nc_get_var1_double (ncid, ids[i], &item[0], &dummy[0]));
				gmt_M_err_trap (nc_get_var1_double (ncid, ids[i], &item[1], &dummy[1]));
			}
			if (item[1] != 0 && dummy[0] != dummy[1]) { /* avoid dvision by 0 */
				double index = (t_value[i] - dummy[0]) / (dummy[1] - dummy[0]) * item[1];
				if (index > 0)
					HH->t_index[i] = (size_t)index;
			}
		}

#ifdef NC4_DEBUG
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->wesn: %g %g %g %g\n",
		            header->wesn[XLO], header->wesn[XHI], header->wesn[YLO], header->wesn[YHI]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->registration:%u\n", header->registration);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->row_order: %s\n",
		            HH->row_order == k_nc_start_south ? "S->N" : "N->S");
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->n_columns: %3d   head->n_rows:%3d\n", header->n_columns, header->n_rows);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->inc: %g %g\n", header->inc[GMT_X], header->inc[GMT_Y]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->mx: %3d   head->my:%3d\n", header->mx, header->my);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->nm: %3d head->size:%3d\n", header->nm, header->size);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->t-index %d,%d,%d\n",
		            HH->t_index[0], HH->t_index[1], HH->t_index[2]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->pad xlo:%u xhi:%u ylo:%u yhi:%u\n",
		            header->pad[XLO], header->pad[XHI], header->pad[YLO], header->pad[YHI]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->BC  xlo:%u xhi:%u ylo:%u yhi:%u\n",
		            HH->BC[XLO], HH->BC[XHI], HH->BC[YLO], HH->BC[YHI]);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->grdtype:%u %u\n", HH->grdtype, GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT);
#endif
	}
	else {	/* Here we are writing */
		/* Store global attributes */
		const int *nc_vers = gmtnc_netcdf_libvers();
		GMT_Report (GMT->parent, GMT_MSG_DEBUG, "netCDF Library version: %d\n", *nc_vers);
		gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "Conventions", strlen(GMT_NC_CONVENTION), GMT_NC_CONVENTION));
		if (header->title[0]) gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "title", strlen(header->title), header->title));
		gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "history", strlen(header->command), header->command));
		if (header->remark[0]) gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "description", strlen(header->remark), header->remark));
		gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "GMT_version", strlen(GMT_VERSION), (const char *) GMT_VERSION));
		if (header->registration == GMT_GRID_PIXEL_REG) {
			int reg = header->registration;
			gmt_M_err_trap (nc_put_att_int (ncid, NC_GLOBAL, "node_offset", NC_LONG, 1U, &reg));
		}
		else
			nc_del_att (ncid, NC_GLOBAL, "node_offset");

#ifdef HAVE_GDAL
		/* If we have projection information create a container variable named "grid_mapping" with an attribute
		   "spatial_ref" that will hold the projection info in WKT format. GDAL and Mirone know use this info */
		if ((header->ProjRefWKT != NULL) || (header->ProjRefPROJ4 != NULL)) {
			int id[1], dim[1];

			if (header->ProjRefWKT == NULL || !header->ProjRefWKT[0]) {				/* Must convert from proj4 string to WKT */
				OGRSpatialReferenceH hSRS = OSRNewSpatialReference(NULL);

				if (header->ProjRefPROJ4 && (!strncmp(header->ProjRefPROJ4, "+unavailable", 4) || strlen(header->ProjRefPROJ4) <= 5)) {	/* Silently jump out of here */
					OSRDestroySpatialReference(hSRS);
					goto L100;
				}
				GMT_Report(GMT->parent, GMT_MSG_INFORMATION, "Proj4 string to be converted to WKT:\n\t%s\n", header->ProjRefPROJ4);
				if (OSRImportFromProj4(hSRS, header->ProjRefPROJ4) == CE_None) {
					char *pszPrettyWkt = NULL;
					OSRExportToPrettyWkt(hSRS, &pszPrettyWkt, false);
					header->ProjRefWKT = strdup(pszPrettyWkt);
					CPLFree(pszPrettyWkt);
					GMT_Report(GMT->parent, GMT_MSG_INFORMATION, "WKT converted from proj4 string:\n%s\n", header->ProjRefWKT);
				}
				else {
					header->ProjRefWKT = NULL;
					GMT_Report(GMT->parent, GMT_MSG_WARNING, "gmtnc_grd_info failed to convert the proj4 string\n%s\n to WKT\n",
							header->ProjRefPROJ4);
				}
				OSRDestroySpatialReference(hSRS);
			}

			if (header->ProjRefWKT != NULL) {			/* It may be NULL if the above conversion failed */
				if (nc_inq_varid(ncid, "grid_mapping", &id[0]) != NC_NOERR) {
					gmt_M_err_trap(nc_def_dim(ncid, "grid_mapping", 12U, &dim[0]));
					gmt_M_err_trap(nc_def_var(ncid, "grid_mapping", NC_CHAR,  1, dim, &id[0]));
				}
				gmt_M_err_trap(nc_put_att_text(ncid, id[0], "spatial_ref", strlen(header->ProjRefWKT), header->ProjRefWKT));
				gmt_M_err_trap(nc_put_att_text(ncid, z_id, "grid_mapping", 12U, "grid_mapping"));	/* Create attrib in z variable */
			}
		}
L100:
#endif

		/* Avoid NaN increments */
		if (gmt_M_is_dnan(header->inc[GMT_X])) header->inc[GMT_X] = 1.0;
		if (gmt_M_is_dnan(header->inc[GMT_Y])) header->inc[GMT_Y] = 1.0;

		/* Define x variable */
		gmtnc_put_units (ncid, ids[HH->xy_dim[0]], header->x_units);
		dummy[0] = header->wesn[XLO], dummy[1] = header->wesn[XHI];
		gmt_M_err_trap (nc_put_att_double (ncid, ids[HH->xy_dim[0]], "actual_range", NC_DOUBLE, 2U, dummy));
		gmt_M_err_trap (nc_put_att_text (ncid, ids[HH->xy_dim[0]], "axis", 1U, "X"));

		/* Define y variable */
		gmtnc_put_units (ncid, ids[HH->xy_dim[1]], header->y_units);
		dummy[0] = header->wesn[YLO], dummy[1] = header->wesn[YHI];
		gmt_M_err_trap (nc_put_att_double (ncid, ids[HH->xy_dim[1]], "actual_range", NC_DOUBLE, 2U, dummy));
		gmt_M_err_trap (nc_put_att_text (ncid, ids[HH->xy_dim[1]], "axis", 1U, "Y"));

		/* When varname is given, and z_units is default, overrule z_units with varname */
		if (HH->varname[0] && !strcmp (header->z_units, cube ? "cube" : "z"))
			strncpy (header->z_units, HH->varname, GMT_GRID_UNIT_LEN80-1);

		/* Define z variable. Attempt to remove "scale_factor" or "add_offset" when no longer needed */
		gmtnc_put_units (ncid, z_id, header->z_units);

		if (header->z_scale_factor != 1.0) {
			gmt_M_err_trap (nc_put_att_double (ncid, z_id, "scale_factor", NC_DOUBLE, 1U, &header->z_scale_factor));
		}
		else if (job == 'u')
			nc_del_att (ncid, z_id, "scale_factor");

		if (header->z_add_offset != 0.0) {
			gmt_M_err_trap (nc_put_att_double (ncid, z_id, "add_offset", NC_DOUBLE, 1U, &header->z_add_offset));
		}
		else if (job == 'u')
			nc_del_att (ncid, z_id, "add_offset");

		if (z_type != NC_FLOAT && z_type != NC_DOUBLE)
			header->nan_value = rintf (header->nan_value); /* round to integer */
		if (job == 'u' && HH->is_netcdf4) {
			/* netCDF-4 has a bug and crash when * rewriting the _FillValue attribute in netCDF-4 files
			   https://bugtracking.unidata.ucar.edu/browse/NCF-187
			   To work-around it we implement the renaming trick advised on NCF-133
         Edit: This work-around should eventually be removed because the bug was fixed as of 2015-04-02 (any version >4.3.3.1 should be fine).
			*/
			gmt_M_err_trap (nc_rename_att (ncid, z_id, "_FillValue", "_fillValue"));
#ifdef DOUBLE_PRECISION_GRID
			gmt_M_err_trap (nc_put_att_double (ncid, z_id, "_fillValue", z_type, 1U, &header->nan_value));
#else
			gmt_M_err_trap (nc_put_att_float (ncid, z_id, "_fillValue", z_type, 1U, &header->nan_value));
#endif
			gmt_M_err_trap (nc_rename_att (ncid, z_id, "_fillValue", "_FillValue"));
		}
		else {
#ifdef DOUBLE_PRECISION_GRID
			gmt_M_err_trap (nc_put_att_double (ncid, z_id, "_FillValue", z_type, 1U, &header->nan_value));
#else
			gmt_M_err_trap (nc_put_att_float (ncid, z_id, "_FillValue", z_type, 1U, &header->nan_value));
#endif
		}

		/* Limits need to be stored in actual, not internal grid, units */
		if (header->z_min <= header->z_max) {
			dummy[0] = header->z_min * header->z_scale_factor + header->z_add_offset;
			dummy[1] = header->z_max * header->z_scale_factor + header->z_add_offset;
		}
		else
			dummy[0] = 0.0, dummy[1] = 0.0;
		gmt_M_err_trap (nc_put_att_double (ncid, z_id, "actual_range", NC_DOUBLE, 2U, dummy));

		/* Store values along x and y axes */
		nc_inq_nvars (ncid, &nvars);
		for (j = 0; j < nvars; j++) {
			gmt_M_err_trap (nc_inq_varndims (ncid, j, &ndims));
		}
		if (!cube) {	/* Hold off on this until we also have the 3rd dimension */
			gmt_M_err_trap (nc_enddef (ncid));	/* End definition mode */
			gmtnc_put_xy_vectors (GMT, header);	/* Place x and y vectors */
		}
	}

	/* Close NetCDF file, unless job == 'W' */

	if (job != 'W') gmt_M_err_trap (gmt_nc_close (GMT, ncid));
	return (GMT_NOERROR);
}

/* Shift columns in a grid to the right (n_shift < 0) or to the left (n_shift < 0) */
GMT_LOCAL void gmtnc_right_shift_grid (void *gridp, const unsigned n_cols, const unsigned n_rows, int n_shift, size_t cell_size) {
	char *tmp = NULL, *grid = (char*)gridp;
	size_t row, n_shift_abs = abs (n_shift), nm;
	size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows;

	assert (n_shift_abs != 0 && n_cols_t > n_shift_abs && n_cols_t > 0 && n_rows_t > 0);

	tmp = calloc (n_shift_abs, cell_size);

	if (n_shift > 0) { /* right shift */
		for (row = 0; row < n_rows_t; ++row) {
			nm = row * n_cols_t;
			/* copy last n_shift_abs cols into tmp buffer */
			memcpy (tmp, grid + (nm + n_cols_t - n_shift_abs) * cell_size, n_shift_abs * cell_size);
			/* right shift row */
			memmove (grid + (nm + n_shift_abs) * cell_size,
							 grid + nm * cell_size,
							 (n_cols_t - n_shift_abs) * cell_size);
			/* prepend tmp buffer */
			memcpy (grid + nm * cell_size, tmp, n_shift_abs * cell_size);
		}
	}
	else { /* n_shift_abs < 0 */
		for (row = 0; row < n_rows_t; ++row) {
			nm = row * n_cols_t;
			/* copy first n_shift_abs cols into tmp buffer */
			memcpy (tmp, grid + nm * cell_size, n_shift_abs * cell_size);
			/* left shift row */
			memmove (grid + nm * cell_size,
							 grid + (nm + n_shift_abs) * cell_size,
							 (n_cols_t - n_shift_abs) * cell_size);
			/* append tmp buffer */
			memcpy (grid + (nm + n_cols_t - n_shift_abs) * cell_size, tmp, n_shift_abs * cell_size);
		}
	}
	gmt_M_str_free (tmp);
}

/* Fill padding by replicating the border cells or wrapping around a
 * row if columns are periodic */
GMT_LOCAL void gmtnc_padding_copy (void *gridp, const unsigned n_cols, const unsigned n_rows, const unsigned *n_pad, size_t cell_size, bool periodic_cols) {
	/* n_cols and n_rows are dimensions of the padded grid */
	char *grid = (char*)gridp;
	size_t row, cell, nm, nm2;
	size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows;
	size_t pad_t[4];

	assert (n_cols > n_pad[XLO] + n_pad[XHI] && n_rows > n_pad[YLO] + n_pad[YHI] &&
		n_pad[XLO] + n_pad[XHI] + n_pad[YLO] + n_pad[YHI] > 0 && cell_size > 0);

	for (row = 0; row < 4; row++) pad_t[row] = (size_t)n_pad[row];

	if (periodic_cols) {
		/* A periodic grid wraps around */
		for (row = pad_t[YHI]; (row + pad_t[YLO]) < n_rows_t; ++row) {
			nm = row * n_cols_t;
			/* Iterate over rows that contain data */
			for (cell = 0; cell < pad_t[XLO]; ++cell) {
				/* Copy end of this row into first pad_t[XLO] columns:
				 * X X 0 1 2 3 4 5 X X -> 4 5 0 1 2 3 4 5 X X */
				memcpy (grid + (nm + cell) * cell_size,
								grid + (nm + n_cols_t + cell - pad_t[XLO] - pad_t[XHI]) * cell_size,
								cell_size);
			}
			for (cell = 0; cell < pad_t[XHI]; ++cell) {
				/* Copy start of this row into last n_pad[XHI] columns:
				 * 4 5 0 1 2 3 4 5 X X -> 4 5 0 1 2 3 4 5 0 1 */
				memcpy (grid + (nm + n_cols_t - cell - 1) * cell_size,
								grid + (nm + pad_t[XLO] + pad_t[XHI] - cell - 1) * cell_size,
								cell_size);
			}
		}
	}
	else { /* !periodic_cols */
		for (row = pad_t[YHI]; (row + pad_t[YLO]) < n_rows_t; ++row) {
			nm = row * n_cols_t;
			/* Iterate over rows that contain data */
			for (cell = 0; cell < pad_t[XLO]; ++cell) {
				/* Duplicate first pad_t[XLO] columns in this row:
				 * 4 5 0 1 2 3 4 5 X X -> 0 0 0 1 2 3 4 5 X X */
				memcpy (grid + (nm + cell) * cell_size,
								grid + (nm + pad_t[XLO]) * cell_size,
								cell_size);
			}
			for (cell = 0; cell < pad_t[XHI]; ++cell) {
				/* Duplicate last pad_t[XHI] columns in this row:
				 * 0 0 0 1 2 3 4 5 X X -> 0 0 0 1 2 3 4 5 5 5 */
				memcpy (grid + (nm + n_cols_t - cell - 1) * cell_size,
								grid + (nm + n_cols_t - pad_t[XHI] - 1) * cell_size,
								cell_size);
			}
		}
	}

	for (cell = 0; cell < pad_t[YHI]; ++cell) {
		nm = cell * n_cols_t;
		/* Duplicate pad_t[YHI] rows in the beginning */
		memcpy (grid + nm * cell_size,
					 grid + pad_t[YHI] * n_cols_t * cell_size,
					 n_cols_t * cell_size);
	}
	nm2 = (n_rows_t - pad_t[YLO] - 1) * n_cols_t;
	for (cell = 0; cell < pad_t[YLO]; ++cell) {
		nm = (n_rows_t - cell - 1) * n_cols_t;
		/* Duplicate last pad_t[YLO] rows */
		memcpy (grid + nm * cell_size,
					 grid + nm2 * cell_size,
					 n_cols_t * cell_size);
	}
}

/* Fill padding with zeros */
GMT_LOCAL void gmtnc_padding_zero (void *gridp, const unsigned n_cols, const unsigned n_rows, const unsigned *n_pad, size_t cell_size) {
	/* n_cols and n_rows are dimensions of the padded grid */
	char *grid = (char*)gridp;
	size_t row, nm;
	size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows;
	size_t pad_t[4];

	assert (n_cols > n_pad[XLO] + n_pad[XHI] && n_rows > n_pad[YLO] + n_pad[YHI] &&
		n_pad[XLO] + n_pad[XHI] + n_pad[YLO] + n_pad[YHI] > 0 && cell_size > 0);

	for (row = 0; row < 4; row++) pad_t[row] = (size_t)n_pad[row];

	/* Iterate over rows that contain data */
	for (row = pad_t[YHI]; (row + pad_t[YLO]) < n_rows_t; ++row) {
		nm = row * n_cols_t;
		/* Zero n cells at beginning of row */
		memset (grid + nm * cell_size, 0, pad_t[XLO] * cell_size);
		/* Zero n cells at end of row */
		memset (grid + (nm + n_cols_t - pad_t[XHI]) * cell_size, 0, pad_t[XHI] * cell_size);
	}
	/* Zero pad_t[YHI] rows in the beginning */
	memset (grid, 0, pad_t[YHI] * n_cols_t * cell_size);
	/* Zero last pad_t[YLO] rows */
	memset(grid + (n_rows_t-pad_t[YLO]) * n_cols_t * cell_size, 0, pad_t[YLO] * n_cols_t * cell_size);
}

/* Fill mode for grid padding */
enum Grid_padding_mode {
	k_pad_fill_none = 0, /* Leave padded cells untouched */
	k_pad_fill_zero,     /* Fill padded grid cells with zeros */
	k_pad_fill_copy,     /* Padded cells get the value of their nearest neighbor */
	k_pad_fill_copy_wrap /* Padded cells get wrapped values from the other side of the row (gridswith periodic columns) */
};

/* Add padding to a matrix/grid and reshape data */
GMT_LOCAL void gmtnc_pad_grid (void *gridp, const unsigned n_cols, const unsigned n_rows, const unsigned *n_pad, size_t cell_size, unsigned filltype) {
	/* n_cols and n_rows are dimensions of the grid without padding
	 * cell_size is the size in bytes of each element in grid
	 * n_pad[4] contains the number of cols/rows to pad on each side {W,E,S,N}
	 *
	 * Note: when grid is complex, we pass 2x n_rows */
	char *grid = (char*)gridp;
	size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows;
	size_t new_row;
	size_t old_row = n_rows_t-1;
	size_t n_new_cols = n_cols_t + n_pad[XLO] + n_pad[XHI];
	size_t n_new_rows = n_rows_t + n_pad[YLO] + n_pad[YHI];

#ifdef NC4_DEBUG
	fprintf (stderr, "pad grid w:%u e:%u s:%u n:%u\n",
			n_pad[XLO], n_pad[XHI], n_pad[YLO], n_pad[YHI]);
#endif
	if (n_pad[XLO] + n_pad[XHI] + n_pad[YLO] + n_pad[YHI] == 0)
		return; /* nothing to pad */

	assert (n_cols > 0 && n_rows > 0 && cell_size > 0);

	/* Reshape matrix */
	if (n_pad[XLO] + n_pad[XHI] + n_pad[YHI] != 0) {
		/* When padding W, E, and N (not necessary when padding S only). */
		for (new_row = n_new_rows - n_pad[YLO] - 1; new_row + 1 > n_pad[YHI]; --new_row, --old_row) {
			/* Copy original row to new row, bottom upwards */
			void *from = grid + old_row * n_cols_t * cell_size;
			void *to   = grid + (new_row * n_new_cols + ((size_t)n_pad[XLO])) * cell_size;
			if (n_pad[YHI] == 0) /* rows overlap! */
				memmove (to, from, n_cols_t * cell_size);
			else /* no overlap, memcpy is safe */
				memcpy  (to, from, n_cols_t * cell_size);
		}
	}

	/* Fill padded grid cells */
	switch (filltype) {
		case k_pad_fill_zero:
			gmtnc_padding_zero (grid, (const unsigned)n_new_cols, (const unsigned)n_new_rows, n_pad, cell_size);
			break;
		case k_pad_fill_copy:
			gmtnc_padding_copy (grid, (const unsigned)n_new_cols, (const unsigned)n_new_rows, n_pad, cell_size, false);
			break;
		case k_pad_fill_copy_wrap:
			gmtnc_padding_copy (grid, (const unsigned)n_new_cols, (const unsigned)n_new_rows, n_pad, cell_size, true);
			break;
	}
}

/* Remove padding from a matrix/grid and reshape data */
GMT_LOCAL void gmtnc_unpad_grid (void *gridp, const unsigned n_cols, const unsigned n_rows, const unsigned *n_pad, size_t cell_size) {
	/* n_cols and n_rows are dimensions of the grid without padding
	 * cell_size is the size in bytes of each element in grid
	 * n_pad[4] contains the number of cols/rows to pad on each side {W,E,S,N}
	 *
	 * Note: when grid is complex, we pass 2x n_rows */
	char *grid = (char*)gridp;
	size_t n_old_cols = n_cols + n_pad[XLO] + n_pad[XHI];
	size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows;
	size_t row;

#ifdef NC4_DEBUG
	fprintf (stderr, "unpad grid w:%u e:%u s:%u n:%u\n",
			n_pad[XLO], n_pad[XHI], n_pad[YLO], n_pad[YHI]);
#endif
	if (n_pad[XLO] + n_pad[XHI] + n_pad[YHI] == 0)
		return; /* nothing to unpad (we just ignore n_pad[YLO]) */

	assert (n_cols > 0 && n_rows > 0 && cell_size > 0);

	/* Reshape matrix */
	for (row = 0; row < n_rows_t; ++row) {
		size_t old_row = row + (size_t)n_pad[YHI];
		void *from = grid + (old_row * n_old_cols + ((size_t)n_pad[XLO])) * cell_size;
		void *to   = grid + row * n_cols_t * cell_size;
		/* Copy original row to new row */
		if (n_pad[YHI] == 0) /* rows overlap! */
			memmove (to, from, n_cols_t * cell_size);
		else /* no overlap, memcpy is safe */
			memcpy  (to, from, n_cols_t * cell_size);
	}
}

/* Ensure that repeating columns in geographic gridline registered grids
 * do not contain conflicting information */
GMT_LOCAL void gmtnc_grid_fix_repeat_col (struct GMT_CTRL *GMT, void *gridp, const unsigned n_cols, const unsigned n_rows, size_t cell_size) {
	/* Note: when grid is complex, pass 2x n_rows */
	char *grid = (char*)gridp;
	unsigned n_conflicts = 0;
	size_t row;
	size_t n_cols_t = (size_t)n_cols, n_rows_t = (size_t)n_rows;

	for (row = 0; row < n_rows_t; ++row) {
		char *first = grid + row * n_cols_t * cell_size;                /* first element in row */
		char *last =  grid + (row * n_cols_t + n_cols_t - 1) * cell_size; /* last element in row */
		if ( memcmp(last, first, cell_size) ) {
			/* elements differ: replace value of last element in row with value of first */
			memcpy (last, first, cell_size);
			++n_conflicts;
		}
	}

	if (n_conflicts)
		GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Detected %u inconsistent values along east boundary of grid. Values fixed by duplicating west boundary.\n", n_conflicts);
}

/* Change the default chunk cache settings in the HDF5 library for all variables
 * in the nc-file. The settings apply for subsequent file opens/creates. */
static inline void gmtnc_setup_chunk_cache (void) {
	static bool already_setup = false;
	if (!already_setup) {
		nc_set_chunk_cache (NC_CACHE_SIZE, NC_CACHE_NELEMS, NC_CACHE_PREEMPTION);
		already_setup = true;
	}
}


GMT_LOCAL int gmtnc_grd_prep_io (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, double wesn[4], unsigned int *width, unsigned int *height, int *n_shift, unsigned origin[2], unsigned dim[2], unsigned origin2[2], unsigned dim2[2]) {
	/* Determines which rows and columns to extract from a grid, based on
	 * w,e,s,n.  This routine first rounds the w,e,s,n boundaries to the nearest
	 * gridlines or pixels, then determines the first and last columns and rows,
	 * and the width and height of the subset (in cells).
	 */
	bool is_gridline_reg, is_global, is_global_repeat;
	unsigned last_row, first_col, last_col, first_row;
	unsigned first_col2, last_col2;
	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);

	*n_shift = 0;
	memset (origin2, 0, 2 * sizeof(unsigned));
	memset (dim2, 0, 2 * sizeof(unsigned));

	is_global = gmt_grd_is_global (GMT, header);
	is_global_repeat = HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT;
	is_gridline_reg = header->registration != GMT_GRID_PIXEL_REG;

#ifdef NC4_DEBUG
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "  x-region: %g %g, grid: %g %g\n", wesn[XLO], wesn[XHI], header->wesn[XLO], header->wesn[XHI]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "  y-region: %g %g, grid: %g %g\n", wesn[YLO], wesn[YHI], header->wesn[YLO], header->wesn[YHI]);
#endif

	if (wesn[XLO] == 0 && wesn[XHI] == 0 && wesn[YLO] == 0 && wesn[YHI] == 0) {
		/* When -R... matches the exact dimensions of the grid: read whole grid */
		*width  = header->n_columns;
		*height = header->n_rows;
		first_col = first_row = 0;
		last_col  = header->n_columns - 1;
		last_row  = header->n_rows - 1;
		gmt_M_memcpy (wesn, header->wesn, 4, double);
	}
	else {
		/* Must deal with a subregion */
		double x;
		x = fabs (header->wesn[YLO] - wesn[YLO]);	/* if |x| < GMT_CONV4_LIMIT * header->inc[GMT_Y] we set wesn to the grid limit */
		if (x > 0.0 && x < GMT_CONV4_LIMIT * header->inc[GMT_Y]) wesn[YLO] = header->wesn[YLO];	/* Avoid snafu */
		x = fabs (header->wesn[YHI] - wesn[YHI]);	/* if |x| < GMT_CONV4_LIMIT * header->inc[GMT_Y] we set wesn to the grid limit */
		if (x > 0.0 && x < GMT_CONV4_LIMIT * header->inc[GMT_Y]) wesn[YHI] = header->wesn[YHI];	/* Avoid snafu */
		if (wesn[YLO] < header->wesn[YLO] || wesn[YHI] > header->wesn[YHI])
			return (GMT_GRDIO_DOMAIN_VIOLATION);	/* Calling program goofed... */

		/* Make sure w,e,s,n are proper multiples of x_inc,y_inc away from x_min,y_min */
		gmt_M_err_pass (GMT, gmt_adjust_loose_wesn (GMT, wesn, header), HH->name);

		/* Global grids: ensure that wesn >= header->wesn (w+e only) */
		if ( is_global ) {	/* Deal with wrapping issues */
			while (wesn[XLO] > header->wesn[XHI]) {
				wesn[XLO] -= 360;
				wesn[XHI] -= 360;
			}
			while (wesn[XLO] < header->wesn[XLO]) {
				wesn[XLO] += 360;
				wesn[XHI] += 360;
			}
		}
		else /* Apply a 10^-4 times inc sloppiness in the test [it was zero which was too harsh for sloppy grids] */
			assert ((wesn[XLO]+GMT_CONV4_LIMIT*header->inc[GMT_X]) >= header->wesn[XLO] && (wesn[XHI]-GMT_CONV4_LIMIT*header->inc[GMT_X]) <= header->wesn[XHI]);

		/* Get dimension of subregion */
		*width  = urint ((wesn[XHI] - wesn[XLO]) * HH->r_inc[GMT_X]) + is_gridline_reg;
		*height = urint ((wesn[YHI] - wesn[YLO]) * HH->r_inc[GMT_Y]) + is_gridline_reg;

		/* Get first and last row and column numbers */
		first_col = urint ((wesn[XLO] - header->wesn[XLO]) * HH->r_inc[GMT_X]);
		last_col  = urint ((wesn[XHI] - header->wesn[XLO]) * HH->r_inc[GMT_X]) + is_gridline_reg - 1;
		first_row = urint ((header->wesn[YHI] - wesn[YHI]) * HH->r_inc[GMT_Y]);
		last_row  = urint ((header->wesn[YHI] - wesn[YLO]) * HH->r_inc[GMT_Y]) + is_gridline_reg - 1;

		/* Adjust first_row */
		if (HH->row_order == k_nc_start_south)
			first_row = header->n_rows - 1 - last_row;

		/* Global grids: if -R + padding is equal or exceeds grid bounds */
		if (is_global && 1 + is_global_repeat + last_col - first_col >= header->n_columns) {
			/* Number of requested cols >= n_columns: read whole grid and shift */
			*n_shift  = -(int)first_col;
			first_col = 0;
			last_col  = header->n_columns - 1;
		}
		else if (is_global && last_col + 1 > header->n_columns) {
			/* Subset of a global grid that wraps around east boundary. This means
			 * we have to read the grid in two parts and stich them together. */
			first_col2 = 0;
			last_col2 = is_global_repeat + last_col - header->n_columns;
			last_col = header->n_columns - 1;
#ifdef NC4_DEBUG
			GMT_Report (GMT->parent, GMT_MSG_WARNING, "col2: %u %u\n", first_col2, last_col2);
#endif

			origin2[0] = first_row;
			origin2[1] = first_col2;
			dim2[0] = *height;
			assert (first_col2 <= 1 + last_col2); /* check for unsigned overflow */
			dim2[1] = 1 + last_col2 - first_col2;
		}
		assert (last_col + 1 <= header->n_columns);
	}

	/* Do not read last column from global repeating grids */
	if (is_global_repeat && last_col + 1 == header->n_columns)
		--last_col;

	origin[0] = first_row;
	origin[1] = first_col;
	dim[0] = *height;
	assert (first_col <= 1 + last_col); /* check for unsigned overflow */
	dim[1] = 1 + last_col - first_col;

#ifdef NC4_DEBUG
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "-> x-region: %g %g, grid: %g %g\n", wesn[XLO], wesn[XHI], header->wesn[XLO], header->wesn[XHI]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "-> y-region: %g %g, grid: %g %g\n", wesn[YLO], wesn[YHI], header->wesn[YLO], header->wesn[YHI]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "row: %u %u  col: %u %u  r_shift: %d\n", first_row, last_row, first_col, last_col, *n_shift);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "origin : %u,%u  dim : %u,%u\n", origin[0], origin[1], dim[0], dim[1]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "origin2: %u,%u  dim2: %u,%u\n", origin2[0], origin2[1], dim2[0], dim2[1]);
#endif

	return GMT_NOERROR;
}

int gmt_nc_create (struct GMT_CTRL *GMT, char *path, int mode, int *id) {
	int err = nc_create (path, mode, id);
	GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Calling nc_create on %s, ncid = %d, err = %d\n", path, *id, err);
	return (err);
}

int gmt_nc_open (struct GMT_CTRL *GMT, char *path, int mode, int *id) {
	int err = nc_open (path, mode, id);	/* Open the file, return error code */
	GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Calling nc_open on %s, ncid = %d, err = %d\n", path, *id, err);
	return err;	/* Open the file, return error code */
}

int gmt_nc_close (struct GMT_CTRL *GMT, int id) {
	int err = nc_close (id);
	GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Calling nc_close on ncid %d, err = %d\n", id, err);
	return err;
}

int gmtlib_is_nc_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
	/* Returns GMT_NOERROR if NetCDF grid */
	int ncid, z_id = GMT_NOTSET, j = 0, nvars, ndims, err, old = false;
	nc_type z_type;
	char varname[GMT_GRID_VARNAME_LEN80];
	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);

	/* Extract levels name from variable name */
	strncpy (varname, HH->varname, GMT_GRID_VARNAME_LEN80-1);
	if (varname[0]) {
		j = 0;
		while (varname[j] && varname[j] != '[' && varname[j] != '(') j++;
		if (varname[j])
			varname[j] = '\0';
	}
	if (!strcmp (HH->name, "="))
		return (GMT_GRDIO_NC_NO_PIPE);

	/* Open the file and look for the required variable */
	if (gmt_access (GMT, HH->name, F_OK))
		return (GMT_GRDIO_FILE_NOT_FOUND);
	if ((err = gmt_nc_open (GMT, HH->name, NC_NOWRITE, &ncid)))
		return (GMT_GRDIO_OPEN_FAILED);
	if (!nc_inq_dimid (ncid, "xysize", &z_id)) {
		/* Old style GMT netCDF grid */
		old = true;
		if (nc_inq_varid (ncid, "z", &z_id))
			return (GMT_GRDIO_NO_VAR);
	}
	else if (varname[0]) {
		/* ?<varname> used */
		if (nc_inq_varid (ncid, varname, &z_id))
			return (GMT_GRDIO_NO_VAR);
	}
	else {	/* Look for first 2D grid, with fallback to first higher-dimension (3-4D) grid if 2D not found */
		int ID = GMT_NOTSET, dim = 0;
		nc_inq_nvars (ncid, &nvars);
		for (j = 0; j < nvars && z_id < 0; j++) {
			gmt_M_err_trap (nc_inq_varndims (ncid, j, &ndims));
			if (ndims == 2)	/* Found the first 2-D grid */
				z_id = j;
			else if (ID == GMT_NOTSET && ndims > 2 && ndims < 5) {	/* Also look for higher-dim grid in case no 2-D */
				ID = j;
				dim = ndims;
			}
		}
		if (z_id < 0) {	/* No 2-D grid found, check if we found a higher dimension cube */
			if (ID == GMT_NOTSET) return (GMT_GRDIO_NO_2DVAR);	/* No we didn't */
			z_id = ID;	/* Pick the higher dimensioned cube instead, get its name, and warn */
			nc_inq_varname (ncid, z_id, varname);
			GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No 2-D array in file %s.  Selecting first 2-D slice in the %d-D array %s\n", HH->name, dim, varname);
		}
	}

	gmt_M_err_trap (nc_inq_vartype (ncid, z_id, &z_type));
	switch (z_type) {
		case NC_UBYTE:  header->type = old ? GMT_GRID_IS_CB : GMT_GRID_IS_NB; HH->orig_datatype = GMT_UCHAR;  break;
		case NC_BYTE:   header->type = old ? GMT_GRID_IS_CB : GMT_GRID_IS_NB; HH->orig_datatype = GMT_CHAR;   break;
		case NC_SHORT:  header->type = old ? GMT_GRID_IS_CS : GMT_GRID_IS_NS; HH->orig_datatype = GMT_SHORT;  break;
		case NC_INT:    header->type = old ? GMT_GRID_IS_CI : GMT_GRID_IS_NI; HH->orig_datatype = GMT_INT;    break;
		case NC_FLOAT:  header->type = old ? GMT_GRID_IS_CF : GMT_GRID_IS_NF; HH->orig_datatype = GMT_FLOAT;  break;
		case NC_DOUBLE: header->type = old ? GMT_GRID_IS_CD : GMT_GRID_IS_ND; HH->orig_datatype = GMT_DOUBLE; break;
		default:        header->type = k_grd_unknown_fmt; break;
	}
	gmt_nc_close (GMT, ncid);
	return GMT_NOERROR;
}

int gmt_nc_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
	return (gmtnc_grd_info (GMT, header, 'r', false, 0));
}

int gmt_nc_update_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
	return (gmtnc_grd_info (GMT, header, 'u', false, 0));
}

int gmt_nc_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
	return (gmtnc_grd_info (GMT, header, 'w', false, 0));
}

int gmt_nc_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
	/* header:       grid structure header
	 * grid:         array with final grid
	 * wesn:         Sub-region to extract  [Use entire file if 0,0,0,0]
	 * padding:      # of empty rows/columns to add on w, e, s, n of grid, respectively
	 * complex_mode:	&4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only)
	 *		Note: The file has only real values, we simply allow space in the complex array
	 *		for real and imaginary parts when processed by grdfft etc.
	 *
	 * Reads a subset of a grdfile and optionally pads the array with extra rows and columns
	 * header values for n_columns and n_rows are reset to reflect the dimensions of the logical array,
	 * not the physical size (i.e., the padding is not counted in n_columns and n_rows)
	 */

	bool adj_nan_value; /* if we need to change the fill value */
	int err;            /* netcdf errors */
	int n_shift;
	int error;
	unsigned dim[2], dim2[2], origin[2], origin2[2]; /* dimension and origin {y,x} of subset to read from netcdf */
	unsigned width = 0, height = 0;
	size_t width_t, height_t, row, stride_t;
	uint64_t imag_offset;
	gmt_grdfloat *pgrid = NULL;
	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);

	/* Check type: is file in old NetCDF format or not at all? */
	if (GMT->session.grdformat[header->type][0] == 'c')
		return (gmt_cdf_read_grd (GMT, header, grid, wesn, pad, complex_mode));
	else if (GMT->session.grdformat[header->type][0] != 'n')
		return (NC_ENOTNC);

	if ((error = gmt_M_err_fail (GMT, gmtnc_grd_prep_io (GMT, header, wesn, &width, &height, &n_shift, origin, dim, origin2, dim2), HH->name)))
		return (error);

	/* Set stride and offset if complex */
	(void)gmtlib_init_complex (header, complex_mode, &imag_offset);	/* Set offset for imaginary complex component */
	pgrid = grid + imag_offset;	/* Start of this complex component (or start of non-complex grid) */

#ifdef NC4_DEBUG
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "      wesn: %g %g %g %g\n", wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->wesn: %g %g %g %g\n", header->wesn[XLO], header->wesn[XHI], header->wesn[YLO], header->wesn[YHI]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->registration:%u\n", header->registration);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->row_order: %s\n", HH->row_order == k_nc_start_south ? "S->N" : "N->S");
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "width:    %3d     height:%3d\n", width, height);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->n_columns: %3d   head->n_rows:%3d\n", header->n_columns, header->n_rows);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->mx: %3d   head->my:%3d\n", header->mx, header->my);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->nm: %3d head->size:%3d\n", header->nm, header->size);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->t-index %d,%d,%d\n", HH->t_index[0], HH->t_index[1], HH->t_index[2]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "      pad xlo:%u xhi:%u ylo:%u yhi:%u\n", pad[XLO], pad[XHI], pad[YLO], pad[YHI]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->pad xlo:%u xhi:%u ylo:%u yhi:%u\n", header->pad[XLO], header->pad[XHI], header->pad[YLO], header->pad[YHI]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->BC  xlo:%u xhi:%u ylo:%u yhi:%u\n", HH->BC[XLO], HH->BC[XHI], HH->BC[YLO], HH->BC[YHI]);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "head->grdtype:%u %u\n", HH->grdtype, GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT);
	GMT_Report (GMT->parent, GMT_MSG_WARNING, "imag_offset: %" PRIu64 "\n", imag_offset);
#endif

	/* Open the NetCDF file */
	if (!strcmp (HH->name,"="))
		return (GMT_GRDIO_NC_NO_PIPE);
	gmtnc_setup_chunk_cache();
	gmt_M_err_trap (gmt_nc_open (GMT, HH->name, NC_NOWRITE, &HH->ncid));

	/* Read grid */
	if (dim2[1] == 0)
		gmtnc_io_nc_grid (GMT, header, dim, origin, HH->stride, k_get_netcdf, pgrid + HH->data_offset, false);
	else {
		/* Read grid in two parts */
		unsigned int stride_or_width = HH->stride != 0 ? HH->stride : width;
		gmtnc_io_nc_grid (GMT, header, dim, origin, stride_or_width, k_get_netcdf, pgrid + HH->data_offset, false);
		gmtnc_io_nc_grid (GMT, header, dim2, origin2, stride_or_width, k_get_netcdf, pgrid + HH->data_offset + dim[1], false);
	}

	/* If we need to shift grid */
	if (n_shift) {
#ifdef NC4_DEBUG
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtnc_right_shift_grid: %d\n", n_shift);
#endif
		gmtnc_right_shift_grid (pgrid, dim[1], dim[0], n_shift, sizeof(grid[0]));
	}

	/* If dim[1] + dim2[1] was < requested width: wrap-pad east border */
	if (gmt_grd_is_global(GMT, header) && width > dim[1] + dim2[1]) {
		unsigned int fix_pad[4] = {0,0,0,0};
		fix_pad[XHI] = width - dim[1] - dim2[1];
		gmtnc_pad_grid (pgrid, width - fix_pad[XHI], height, fix_pad, sizeof(grid[0]), k_pad_fill_copy_wrap);
	}
	else
		assert (width == dim[1] + dim2[1]);

	/* Get stats */
	width_t = (size_t)width;
	height_t = (size_t)height;
	stride_t = (size_t)HH->stride;
	header->z_min = DBL_MAX;
	header->z_max = -DBL_MAX;
	adj_nan_value = !isnan (header->nan_value);
	HH->has_NaNs = GMT_GRID_NO_NANS;	/* We are about to check for NaNs and if none are found we retain 1, else 2 */
	for (row = 0; row < height_t; ++row) {
		gmt_grdfloat *p_data = pgrid + row * (HH->stride ? stride_t : width_t);
		unsigned col;
		for (col = 0; col < width; col ++) {
			if (adj_nan_value && p_data[col] == header->nan_value) {
				p_data[col] = (gmt_grdfloat)NAN;
				HH->has_NaNs = GMT_GRID_HAS_NANS;
				continue;
			}
			else if (!isnan (p_data[col])) {
				header->z_min = MIN (header->z_min, p_data[col]);
				header->z_max = MAX (header->z_max, p_data[col]);
			}
			else
				HH->has_NaNs = GMT_GRID_HAS_NANS;
		}
	}
	/* Check limits */
	if (header->z_min > header->z_max) {
		header->z_min = NAN;
		header->z_max = NAN;
		GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No valid values in grid [%s]\n", HH->name);
	}
	else {
		/* Report z-range of grid (with scale and offset applied): */
		unsigned int level;
#ifdef NC4_DEBUG
		level = GMT_MSG_WARNING;
#else
		level = GMT_MSG_DEBUG;
#endif
		GMT_Report (GMT->parent, level,
				"packed z-range: [%g,%g]\n", header->z_min, header->z_max);
	}

	/* Flip grid upside down */
	if (HH->row_order == k_nc_start_south)
		gmt_grd_flip_vertical (pgrid + HH->data_offset, width, height, HH->stride, sizeof(grid[0]));

	/* Add padding with border replication */
	gmtnc_pad_grid (pgrid, width, height, pad, sizeof(grid[0]), k_pad_fill_zero);

#ifdef NC4_DEBUG
	if (header->size < 160) {
		unsigned n;
		unsigned pad_x = pad[XLO] + pad[XHI];
		unsigned stride = HH->stride ? HH->stride : width;
		gmt_grdfloat *p_data = pgrid + HH->data_offset;
		for (n = 0; n < (stride + pad_x) * (height + pad[YLO] + pad[YHI]); n++) {
			if (n % (stride + pad_x) == 0)
				fprintf (stderr, "\n");
			fprintf (stderr, "%4.0f", p_data[n]);
		}
		fprintf (stderr, "\n");
	}
#endif

	/* Adjust header */
	gmt_M_memcpy (header->wesn, wesn, 4, double);
	header->n_columns = width;
	header->n_rows = height;

	gmt_M_err_trap (gmt_nc_close (GMT, HH->ncid));

	return GMT_NOERROR;
}

int gmt_nc_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode)
{ /* header:       grid structure header
	 * grid:         array with final grid
	 * wesn:         Sub-region to write out  [Use entire file if 0,0,0,0]
	 * padding:      # of empty rows/columns to add on w, e, s, n of grid, respectively
	 * complex_mode:	&4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only)
	 *		Note: The file has only real values, we simply allow space in the complex array
	 *		for real and imaginary parts when processed by grdfft etc.
	 */

	int status = NC_NOERR;
	bool adj_nan_value;   /* if we need to change the fill value */
	bool do_round = true; /* if we need to round to integral */
	unsigned int width, height, *actual_col = NULL;
	unsigned int dim[2], origin[2]; /* dimension and origin {y,x} of subset to write to netcdf */
	int first_col, last_col, first_row, last_row;
	uint64_t imag_offset;
	size_t n, nm;
	size_t width_t, height_t;
	size_t was_chunk_setting = GMT->current.setting.io_nc4_chunksize[0];	/* Remember current setting */
	double limit[2];      /* minmax of z variable */
	gmt_grdfloat *pgrid = NULL;
	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);

	/* Determine the value to be assigned to missing data, if not already done so */
	switch (header->type) {
		case GMT_GRID_IS_NB:
			if (isnan (header->nan_value))
				header->nan_value = NC_MIN_BYTE;
			break;
		case GMT_GRID_IS_NS:
			if (isnan (header->nan_value))
				header->nan_value = NC_MIN_SHORT;
			break;
		case GMT_GRID_IS_NI:
			if (isnan (header->nan_value))
				header->nan_value = NC_MIN_INT;
			break;
		case GMT_GRID_IS_ND:
#ifndef DOUBLE_PRECISION_GRID
			GMT_Report (GMT->parent, GMT_MSG_WARNING, "Precision loss! GMT's internal grid representation is 32-bit float.\n");
#endif
			/* Intentionally fall through */
		default: /* don't round float */
			do_round = false;
	}

	gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width, &height, &first_col, &last_col, &first_row, &last_row, &actual_col), HH->name);
	gmt_M_free (GMT, actual_col);

	/* Adjust header */
	gmt_M_memcpy (header->wesn, wesn, 4, double);
	header->n_columns = width;
	header->n_rows = height;

	/* Adjust first_row */
	if (HH->row_order == k_nc_start_south)
		first_row = header->n_rows - 1 - last_row;

	/* Write grid header without closing file afterwards */
	gmtnc_setup_chunk_cache();
	status = gmtnc_grd_info (GMT, header, 'W', false, 0);
	if (status != NC_NOERR)
		goto nc_err;

	/* Set stride and offset if complex */
	(void)gmtlib_init_complex (header, complex_mode, &imag_offset);	/* Set offset for imaginary complex component */
	pgrid = grid + imag_offset;	/* Beginning of this complex component (or the regular non-complex grid) */

	/* Remove padding from grid */
	gmtnc_unpad_grid (pgrid, width, height, pad, sizeof(grid[0]));

	/* Check that repeating columns do not contain conflicting information */
	if (HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT)
		gmtnc_grid_fix_repeat_col (GMT, pgrid, width, height, sizeof(grid[0]));

	/* Flip grid upside down */
	if (HH->row_order == k_nc_start_south)
		gmt_grd_flip_vertical (pgrid, width, height, 0, sizeof(grid[0]));

	/* Get stats */
	header->z_min = DBL_MAX;
	header->z_max = -DBL_MAX;
	adj_nan_value = !isnan (header->nan_value);
	n = 0;
	width_t  = (size_t)width;
	height_t = (size_t)height;
	nm = width_t * height_t;
	while (n < nm) {
		if (adj_nan_value && isnan (pgrid[n]))
			pgrid[n] = header->nan_value;
		else if (!isnan (pgrid[n])) {
			if (do_round)
				pgrid[n] = rintf (pgrid[n]); /* round to int */
			header->z_min = MIN (header->z_min, pgrid[n]);
			header->z_max = MAX (header->z_max, pgrid[n]);
		}
		n++;
	}

	/* Write grid */
	dim[0]    = height,    dim[1]    = width;
	origin[0] = first_row, origin[1] = first_col;
	status = gmtnc_io_nc_grid (GMT, header, dim, origin, 0, k_put_netcdf, pgrid, false);
	if (status != NC_NOERR)
		goto nc_err;

	if (header->z_min <= header->z_max) {
		/* Warn if z-range exceeds the precision of a single precision float: */
		static const uint32_t exp2_24 = 0x1000000; /* exp2 (24) */
		unsigned int level;
		if (fabs(header->z_min) >= exp2_24 || fabs(header->z_max) >= exp2_24)
			GMT_Report (GMT->parent, GMT_MSG_DEBUG, "The z-range, [%g,%g], might exceed the significand's precision of 24 bits; round-off errors may occur.\n", header->z_min, header->z_max);

		/* Report z-range of grid (with scale and offset applied): */
#ifdef NC4_DEBUG
		level = GMT_MSG_WARNING;
#else
		level = GMT_MSG_DEBUG;
#endif
		GMT_Report (GMT->parent, level,
				"packed z-range: [%g,%g]\n", header->z_min, header->z_max);

		/* Limits need to be written in actual, not internal grid, units: */
		limit[0] = header->z_min * header->z_scale_factor + header->z_add_offset;
		limit[1] = header->z_max * header->z_scale_factor + header->z_add_offset;
	}
	else {
		GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No valid values in grid [%s]\n", HH->name);
		limit[0] = limit[1] = NAN; /* Set limit to NaN */
	}
	status = nc_put_att_double (HH->ncid, HH->z_id, "actual_range", NC_DOUBLE, 2, limit);
	if (status != NC_NOERR)
		goto nc_err;

	/* Close grid */
	status = gmt_nc_close (GMT, HH->ncid);
	if (status != NC_NOERR)
		goto nc_err;

	GMT->current.setting.io_nc4_chunksize[0] = was_chunk_setting;	/* Restore to previous status */

	return GMT_NOERROR;

nc_err:
	/* exit gracefully */
	gmt_nc_close (GMT, HH->ncid); /* close nc-file */
	unlink (HH->name);  /* remove nc-file */
	if (status == NC_ERANGE) {
		/* report out of range z variable */
		GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot write format %s.\n", GMT->session.grdformat[header->type]);
		GMT_Report (GMT->parent, GMT_MSG_ERROR, "The packed z-range, [%g,%g], exceeds the maximum representable size. Adjust scale and offset parameters or remove out-of-range values.\n", header->z_min, header->z_max);
	}
	return status;
}

bool gmt_nc_is_cube (struct GMTAPI_CTRL *API, char *file) {
	/* Return true if the file is a 3-D netCDF cube */
	bool is_cube = false;
	int i, ID = GMT_NOTSET, ncid, z_id = GMT_NOTSET, nvars, ndims = 0;
	char varname[GMT_GRID_NAME_LEN256] = {""}, *c = NULL;
	gmt_M_unused (API);

	if (file == NULL || file[0] == '\0') return false;	/* Got no file */
	if (!strcmp (file,"=")) return (false);	/* NetCDF is not pipeable */
	if (gmt_M_file_is_memory (file)) return ((file[GMTAPI_OBJECT_FAMILY_START] == 'U') ? true : false);	/* Memory cube or not */
	if ((c = strchr (file, '?'))) {	/* Specific variable */
		strcpy (varname, &c[1]);
		c[0] = '\0';	/* Chop off for now */
	}
	if (gmt_nc_open (API->GMT, file, NC_NOWRITE, &ncid)) {
		if (c) c[0] = '?';	/* Restore */
		return false;	/* Unable to open, so probably not netCDF file */
	}
	/* Here file is open so must close before returning */
	if (c) c[0] = '?';	/* Restore */
	if (nc_inq_nvars (ncid, &nvars)) goto nc_cube_return;	/* No variables, no good */
	if (c) {
		if (nc_inq_varid (ncid, varname, &z_id) == NC_NOERR) {	/* Gave a named variable that exist, is it 2-4 D? */
			if (nc_inq_varndims (ncid, z_id, &ndims)) goto nc_cube_return;	/* Not god */
			if (ndims != 3) z_id = GMT_NOTSET;	/* But not 3-D */
			if (ndims > 2 && ndims <= 5) ID = z_id;
		}
	}
	else {	/* Must search */
		i = 0;
		while (i < nvars && z_id < 0) {	/* Look for first 3D grid, with fallback to first higher-dimension (4-%D) grid if 3D not found */
			if (nc_inq_varndims (ncid, i, &ndims)) goto nc_cube_return;	/* Not good */
			if (ndims == 3)	/* Found the first 3-D grid */
				z_id = i;
			else if (ID == GMT_NOTSET && ndims > 3 && ndims < 5)	/* Also look for higher-dim grid in case no 3-D */
				ID = i;
			i++;
		}
	}
	if (ID != GMT_NOTSET || z_id >= 0) is_cube = true;	/* Found  3-D grid */
nc_cube_return:
	gmt_nc_close (API->GMT, ncid); /* close nc-file */
	return (is_cube);
}

/* Examine the netCDF data cube and determine if it is a 3-D cube and return the knots */

int gmt_nc_read_cube_info (struct GMT_CTRL *GMT, char *file, double *w_range, uint64_t *nz, double **zarray, char *z_unit) {
	int i, err, ID = GMT_NOTSET, dim = 0, ncid, z_id = GMT_NOTSET, ids[5] = {-1,-1,-1,-1,-1}, dims[5], nvars;
	int has_vector, ndims = 0, z_dim, status;
	uint64_t n_layers = 0;
	size_t lens[5];
	char varname[GMT_GRID_VARNAME_LEN80], dimname[GMT_GRID_UNIT_LEN80], z_units[GMT_GRID_UNIT_LEN80];
	double *z = NULL, dummy[2] = {0.0, 0.0};

	gmt_M_err_trap (gmt_nc_open (GMT, file, NC_NOWRITE, &ncid));

	gmt_M_err_trap (nc_inq_nvars (ncid, &nvars));
	i = 0;
	while (i < nvars && z_id < 0) {	/* Look for first 3D grid, with fallback to first higher-dimension (4-5D) grid if 3D not found */
		gmt_M_err_trap (nc_inq_varndims (ncid, i, &ndims));
		if (ndims == 3)	/* Found the first 3-D grid */
			z_id = i;
		else if (ID == GMT_NOTSET && ndims > 3 && ndims < 5) {	/* Also look for higher-dim grid in case no 3-D */
			ID = i;
			dim = ndims;
		}
		i++;
	}
	if (z_id < 0) {	/* No 3-D grid found, check if we found a higher dimension cube */
		if (ID == GMT_NOTSET) {	/* No we didn't */
			GMT_Report (GMT->parent, GMT_MSG_ERROR, "No 3-D data cube found in file %s.\n", file);
			return GMT_GRDIO_NO_2DVAR;
		}
		z_id = ID;	/* Pick the higher dimensioned cube instead, get its name, and warn */
		nc_inq_varname (ncid, z_id, varname);
		GMT_Report (GMT->parent, GMT_MSG_WARNING, "No 3-D array in file %s.  Selecting first 3-D slice in the %d-D array %s\n", file, dim, varname);
	}
	gmt_M_err_trap (nc_inq_vardimid (ncid, z_id, dims));

	/* Get the ids of the x and y (and depth and time) coordinate variables */
	for (i = 0; i < ndims; i++) {
		gmt_M_err_trap (nc_inq_dim (ncid, dims[i], dimname, &lens[i]));
		if ((status = nc_inq_varid (ncid, dimname, &ids[i])) != NC_NOERR)
			GMT_Report (GMT->parent, GMT_MSG_WARNING, "\"%s\", %s\n\tIf something bad happens later, try importing via GDAL.\n",
				dimname, nc_strerror(status));
	}
	z_dim = ndims-3;
	n_layers = lens[z_dim];

	/* Create enough memory to store the level-coordinate values */
	z = gmt_M_memory (GMT, NULL, n_layers, double);

	/* Get information about z variable */
	gmtnc_get_units (GMT, ncid, ids[z_dim], z_units);
	if (strstr (z_units, "seconds since 1970-01-01"))
		gmt_set_column_type (GMT, GMT_IN, GMT_Z, GMT_IS_ABSTIME);
	else if (z_units[0] == '\0')
		strcpy (z_units, "z");	/* Set default z-unit name */
	/* Look for the z-coordinate vector */
	if ((has_vector = nc_get_var_double (ncid, ids[z_dim], z))) {
		GMT_Report (GMT->parent, GMT_MSG_ERROR, "No 3rd-dimension coordinate vector found in %s\n", file);
		return GMT_GRDIO_NO_2DVAR;
	}

	if (w_range) gmt_M_memset (w_range, 2U, double);
	if (w_range && (!nc_get_att_double (ncid, z_id, "actual_range", dummy) ||
			!nc_get_att_double (ncid, z_id, "valid_range", dummy) ||
			!(nc_get_att_double (ncid, z_id, "valid_min", &dummy[0]) +
			nc_get_att_double (ncid, z_id, "valid_max", &dummy[1]))))
		gmt_M_memcpy (w_range, dummy, 2U, double);

	*zarray = z;
	*nz = n_layers;
	if (z_unit)	/* Passed pointer to storage for z-unit in cubes */
		strncpy (z_unit, z_units, GMT_GRID_UNIT_LEN80);

	return GMT_NOERROR;
}

/* Write a 3-D cube to file; cube is represented internally by a stack of 2-D grids and a layer z-array */
int gmt_nc_write_cube (struct GMT_CTRL *GMT, struct GMT_CUBE *C, double wesn[], const char *file) {
	/* Depending on mode, we either write individual layer grid files or a single 3-D data cube */
	int status = NC_NOERR, err = 0;
	uint64_t k, k0, k1, n_layers_used, save_n_bands, here = 0;
	struct GMTAPI_CTRL *API = GMT->parent;

	/* Determine which layers we want to write */
	k0 = 0;	k1 = C->header->n_bands - 1;	/* Default is all */
	if (wesn && wesn[ZHI] > wesn[ZLO]) {	/* Want to write a subset of layers */
		if (gmt_get_active_layers (GMT, C, &(wesn[ZLO]), &k0, &k1) == 0) {
			gmtlib_report_error (API, GMT_RUNTIME_ERROR);
		}
	}
	n_layers_used = k1 - k0 + 1;	/* Total number of layers actually to be written */
	if (n_layers_used == 0) {
		GMT_Report (API, GMT_MSG_ERROR, "gmt_nc_write_cube: No layers selected from GMT_IS_CUBE.\n");
		return (gmtlib_report_error (API, GMT_DIM_TOO_SMALL));
	}
	save_n_bands = C->header->n_bands;	/* Remember how many bands */
	here = k0 * C->header->size;	/* Start position in the cube for layer k0 */

	if (strchr (file, '%') || n_layers_used == 1) {	/* Format specifier found, do individual layer grids, or a single layer so no format needed */
		char gfile[PATH_MAX] = {""};
		struct GMT_GRID *G = NULL;

		/* !! Remember that the gmt_nc.c codes will unpad and mess up C->data */
		G = gmt_get_grid (GMT);	/* Need a dummy grid for writing */
		G->header = C->header;	/* Use this pointer for now */
		G->header->n_bands = 1;	/* Grids only have one band [must undo this below] */
		/* Write the layers via a grid */
		for (k = k0; k <= k1; k++) {	/* For all selected output levels */
			if (n_layers_used > 1 || strchr (file, '%'))	/* Create the k'th layer file name from template */
				sprintf (gfile, file, C->z[k]);
			else	/* Just this one layer grid */
				sprintf (gfile, "%s", file);
			G->data = &C->data[here];	/* Point to start of this layer */
			GMT_Report (API, GMT_MSG_DEBUG, "gmt_nc_write_cube: Layer %" PRIu64 ", offset = %" PRIu64 ".\n", k, here);
			if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, wesn, gfile, G) != GMT_NOERROR) {
				return (API->error);
			}
			here += C->header->size;	/* Move to next level */
		}
		/* Wipe and free the temporary grid structure */
		G->header->n_bands = save_n_bands;	/* Restore bands */
		G->data = NULL;
		G->header = NULL;
		gmt_free_grid (GMT, &G, true);
		return (GMT_NOERROR);
	}
	else {	/* Here we must write a 3-D netcdf data cube NOT COMPLETED SO NOT WORKING */
		bool adj_nan_value;   /* if we need to change the fill value */
		bool do_round = true; /* if we need to round to integral */
		unsigned int width, height;
		unsigned int dim[3], origin[2]; /* dimension and origin {y,x} of subset to write to netcdf */
		int first_col, first_row, last_row;
		size_t n, nm;
		size_t width_t, height_t;
		double level_min, level_max;      /* minmax of level variable */
		gmt_grdfloat *pgrid = NULL;
		struct GMT_GRID_HEADER *header = C->header;
		struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);

		width = header->n_columns;	height = header->n_rows;
		header->n_bands = n_layers_used;	/* Number of layers to actually write */
		strncpy (HH->name, file, GMT_GRID_NAME_LEN256);
		/* Determine the value to be assigned to missing data, if not already done so */
		switch (header->type) {
			case GMT_GRID_IS_NB:
				if (isnan (header->nan_value))
					header->nan_value = NC_MIN_BYTE;
				break;
			case GMT_GRID_IS_NS:
				if (isnan (header->nan_value))
					header->nan_value = NC_MIN_SHORT;
				break;
			case GMT_GRID_IS_NI:
				if (isnan (header->nan_value))
					header->nan_value = NC_MIN_INT;
				break;
			case GMT_GRID_IS_ND:
#ifndef DOUBLE_PRECISION_GRID
				GMT_Report (GMT->parent, GMT_MSG_WARNING, "Precision loss! GMT's internal grid representation is 32-bit float.\n");
#endif
				/* Intentionally fall through */
			default: /* don't round float */
				do_round = false;
		}

		first_col = first_row = 0;
		last_row  = header->n_rows - 1;
		level_min = DBL_MAX;
		level_max = -DBL_MAX;

		/* Get stats */
		adj_nan_value = !isnan (header->nan_value);
		width_t  = (size_t)width;
		height_t = (size_t)height;
		nm = width_t * height_t;
		for (k = k0; k <= k1; k++) {	/* For all selected output levels */
			pgrid = &C->data[here];
			/* Update stats and set NaN-values */
			n = 0;
			while (n < nm) {
				if (adj_nan_value && isnan (pgrid[n]))
					pgrid[n] = header->nan_value;
				else if (!isnan (pgrid[n])) {
					if (do_round)
						pgrid[n] = rintf (pgrid[n]); /* round to int */
					level_min = MIN (level_min, pgrid[n]);
					level_max = MAX (level_max, pgrid[n]);
				}
				n++;
			}
			here += C->header->size;
		}
		/* The min/max of the cube */
		header->z_min = level_min;
		header->z_max = level_max;

		/* Adjust first_row */
		if (HH->row_order == k_nc_start_south)
			first_row = header->n_rows - 1 - last_row;

		/* Write grid header without closing file afterwards so more items can be added */
		gmtnc_setup_chunk_cache();
		status = gmtnc_grd_info (GMT, header, 'W', true, n_layers_used);
		if (status != NC_NOERR)
			goto nc_err;

		/* Add cube parameters, etc not done by gmtnc_grd_info */

		/* Define 3-D z variable and write the attributes */

		gmtnc_put_units (HH->ncid, HH->xyz_id[GMT_Z], C->units);

	//	if (C->name[0])
	//		nc_put_att_text (HH->ncid, HH->z_id, "long_name", strlen(C->name), C->name);
	//	else
	//		nc_put_att_text (HH->ncid, HH->z_id, "long_name", 4U, "cube");
	//	if (C->units[0]) {	/* This is the 3rd dimension name and units */
	//		int i = 0, j = 0;
	//		bool copy = false;
	//		char name[GMT_GRID_UNIT_LEN80], units[GMT_GRID_UNIT_LEN80];
	//		strncpy (name, C->units, GMT_GRID_UNIT_LEN80-1);
	//		units[0] = '\0';
	//		for (i = 0; i < GMT_GRID_UNIT_LEN80 && name[i]; i++) {
	//			if (name[i] == ']') copy = false, units[j] = '\0';
	//			if (copy) units[j++] = name[i];
	//			if (name[i] == '[') {
	//				name[i] = '\0';
	//				if (i > 0 && name[i-1] == ' ') name[i-1] = '\0';
	//				copy = true;
	//			}
	//		}
	//		if (name[0]) nc_put_att_text (HH->ncid, HH->xyz_id[GMT_Z], "long_name", strlen(name), name);
	//		if (units[0]) nc_put_att_text (HH->ncid, HH->xyz_id[GMT_Z], "units", strlen(units), units);
	//		if (name[0] && strstr (name, "time"))
	//			nc_put_att_text (HH->ncid, HH->xyz_id[GMT_Z], "standard_name", 4U, "time");
	//		else if (name[0] == '\0')
	//			nc_put_att_text (HH->ncid, HH->xyz_id[GMT_Z], "long_name", 1U, "z");
	//	}
	//	else	/* No information about z, so generic z */
	//		nc_put_att_text (HH->ncid, HH->xyz_id[GMT_Z], "long_name", 1U, "z");
		if (gmt_M_is_geographic (GMT, GMT_OUT))
			nc_put_att_text (HH->ncid, HH->xyz_id[GMT_Z], "axis", 1, "Z");
		gmt_M_err_trap (nc_put_att_double (HH->ncid, HH->xyz_id[GMT_Z], "actual_range", NC_DOUBLE, 2U, C->z_range));

		gmt_M_err_trap (nc_enddef (HH->ncid));	/* Belatedly end definition mode */
		gmtnc_put_xy_vectors (GMT, header);	/* Place x and y vectors */
		{
			size_t z_zero = 0, z_count = n_layers_used;
			gmt_M_err_trap (nc_put_vara_double (HH->ncid, HH->xyz_id[GMT_Z], &z_zero, &z_count, &C->z[k0]));
		}

		/* Get stats */
		//dim[0]    = height,    dim[1]    = width;    dim[2]    = n_layers_used;
		dim[1]    = height,    dim[2]    = width;    dim[0]    = n_layers_used;
		origin[0] = 0, origin[1] = 0;
		here = k0 * C->header->size;	/* Start position in the cube for layer k0 */

		for (k = k0; k <= k1; k++) {	/* For all selected output levels */
			pgrid = &C->data[here];
			HH->t_index[0] = k - k0;	/* Layer index */
			/* Remove padding from grid */
			gmtnc_unpad_grid (pgrid, width, height, header->pad, sizeof(pgrid[0]));

			/* Check that repeating columns do not contain conflicting information */
			if (HH->grdtype == GMT_GRID_GEOGRAPHIC_EXACT360_REPEAT)
				gmtnc_grid_fix_repeat_col (GMT, pgrid, width, height, sizeof(pgrid[0]));

			/* Flip grid upside down */
			if (HH->row_order == k_nc_start_south)
				gmt_grd_flip_vertical (pgrid, width, height, 0, sizeof(pgrid[0]));

			/* Write grid layer */
			status = gmtnc_io_nc_grid (GMT, header, dim, origin, 0, k_put_netcdf, pgrid, true);
			if (status != NC_NOERR)
				goto nc_err;
			here += C->header->size;
		}
		if (level_min <= level_max) {
			/* Warn if level-range exceeds the precision of a single precision float: */
			static const uint32_t exp2_24 = 0x1000000; /* exp2 (24) */
			unsigned int level;
			if (fabs(level_min) >= exp2_24 || fabs(level_max) >= exp2_24)
				GMT_Report (GMT->parent, GMT_MSG_WARNING, "The level-range, [%g,%g], might exceed the significand's precision of 24 bits; round-off errors may occur.\n", level_min, level_max);

			/* Report level-range of grid layer (with scale and offset applied): */
#ifdef NC4_DEBUG
			level = GMT_MSG_WARNING;
#else
			level = GMT_MSG_DEBUG;
#endif
			GMT_Report (GMT->parent, level,
					"packed z-range: [%g,%g]\n", level_min, level_max);

		}
		else
			GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No valid values in grid [%s]\n", HH->name);

		if (status != NC_NOERR)
			goto nc_err;
		/* Close grid */
		status = gmt_nc_close (GMT, HH->ncid);
		if (status != NC_NOERR)
			goto nc_err;

		return GMT_NOERROR;

nc_err:
		/* exit gracefully */
		gmt_nc_close (GMT, HH->ncid); /* close nc-file */
		unlink (HH->name);  /* remove nc-file */
		if (status == NC_ERANGE) {
			/* report out of range z variable */
			GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot write format %s.\n", GMT->session.grdformat[header->type]);
			GMT_Report (GMT->parent, GMT_MSG_ERROR, "The packed z-range, [%g,%g], exceeds the maximum representable size. Adjust scale and offset parameters or remove out-of-range values.\n", level_min, level_max);
		}
		return status;
	}
}