easyNCDF.c

/* The functions in this file are copies of some functions in 'ncdf4' package, 
 * in order for the bugfix in .ncvar_get_inner in easyNCDF to work. When ncdf4 
 * addresses this bug, these function will be removed from this package.
 * This code belongs to David Pierce.
 */
#include <stdio.h>
#include <netcdf.h>
#include <string.h>
#include <stdlib.h>

#include <Rdefines.h>
#include <R_ext/Rdynload.h>

/* These same values are hard-coded into the R source. Don't change them! 
 * Note that they are not the same as values defined in the netcdf
 * library headers, since I don't want my R code to depend on those 
 */
#define R_NC_TYPE_SHORT  1
#define R_NC_TYPE_INT    2
#define R_NC_TYPE_FLOAT  3
#define R_NC_TYPE_DOUBLE 4
#define R_NC_TYPE_TEXT   5
#define R_NC_TYPE_BYTE   6

/* New types in netcdf version 4 */
#define R_NC_TYPE_UBYTE		7
#define R_NC_TYPE_USHORT	8
#define R_NC_TYPE_UINT		9
#define R_NC_TYPE_INT64		10
#define R_NC_TYPE_UINT64	11
#define R_NC_TYPE_STRING	12

int  R_nc4_nctype_to_Rtypecode( nc_type nct );
void R_nc4_varsize( int *ncid, int *varid, int *ndims, int *varsize, int *retval );
void R_nc4_inq_vartype( int *ncid, int *varid, int *precint, int *retval );
void R_nc4_inq_varndims( int *ncid, int *varid, int *ndims, int *retval );
void R_nc4_open( char **filename, int *cmode, int *ncid, int *retval );
void R_nc4_get_vara_text( int *ncid, int *varid, int *start, int *count, char **tempstore, char **data, int *retval );
void R_nc4_close ( int *ncid );
SEXP Rsx_nc4_get_vara_double( SEXP sx_ncid, SEXP sx_varid, SEXP sx_start, SEXP sx_count, SEXP sx_fixmiss, SEXP sx_imvstate, SEXP sx_missval );
SEXP Rsx_nc4_get_vara_int   ( SEXP sx_ncid, SEXP sx_varid, SEXP sx_start, SEXP sx_count, SEXP sx_byte_style );
SEXP R_nc4_blankstring(SEXP size);
SEXP R_nc4_get_vara_string( SEXP sx_nc, SEXP sx_varid, SEXP sx_start, SEXP sx_count );
void R_nc4_inq_varid_hier( int *ncid, char **varname, int *returned_grpid, int *returned_varid );

// Non-registered
int R_nc4_util_nslashes( char *s, int *idx_first_slash );
void R_nc4_inq_varid_hier_inner( int *ncid, char *varname, int *returned_grpid, int *returned_varid );

/* For C calls that don't use SEXP type args */
static const
R_CMethodDef cMethods[] = {
	{"R_nc4_inq_varid_hier", 	(DL_FUNC) &R_nc4_inq_varid_hier, 	4},
	{"R_nc4_nctype_to_Rtypecode", 	(DL_FUNC) &R_nc4_nctype_to_Rtypecode, 	1},
	{"R_nc4_varsize", 		(DL_FUNC) &R_nc4_varsize, 		5},
	{"R_nc4_inq_vartype", 		(DL_FUNC) &R_nc4_inq_vartype, 		4},
	{"R_nc4_inq_varndims", 		(DL_FUNC) &R_nc4_inq_varndims, 		4},
	{"R_nc4_open", 			(DL_FUNC) &R_nc4_open, 			4},
	{"R_nc4_get_vara_text", 	(DL_FUNC) &R_nc4_get_vara_text,  	7},
	{"R_nc4_close", 		(DL_FUNC) &R_nc4_close, 		1},
	NULL 
};

/* For C calls that use SEXP type args */
static const
R_CallMethodDef callMethods[] = {
	{"Rsx_nc4_get_vara_double", 	(DL_FUNC) &Rsx_nc4_get_vara_double,  	7},
	{"Rsx_nc4_get_vara_int", 	(DL_FUNC) &Rsx_nc4_get_vara_int,  	5},
	{"R_nc4_blankstring", 		(DL_FUNC) &R_nc4_blankstring,  		1},
	{"R_nc4_get_vara_string", 	(DL_FUNC) &R_nc4_get_vara_string, 	4},
	NULL
};

/*********************************************************************
 * Register our extnal routines for R
*/
void R_init_ncdf4( DllInfo *info ) 
{
	R_registerRoutines( info, cMethods, callMethods, NULL, NULL );
	R_useDynamicSymbols( info, FALSE );
}

/*********************************************************************
 * Converts from type "nc_type" to an integer as defined in the beginning 
 * of this file.  We do NOT use the raw nc_type integers because then the
 * R code would have a dependency on the arbitrary values in the netcdf 
 * header files!
 */
int R_nc4_nctype_to_Rtypecode( nc_type nct )
{
	if( nct == NC_CHAR )
		return(R_NC_TYPE_TEXT);
	else if( nct == NC_SHORT )
		return(R_NC_TYPE_SHORT);
	else if( nct == NC_INT )
		return(R_NC_TYPE_INT);
	else if( nct == NC_FLOAT )
		return(R_NC_TYPE_FLOAT);
	else if( nct == NC_DOUBLE )
		return(R_NC_TYPE_DOUBLE);
	else if( nct == NC_BYTE )
		return(R_NC_TYPE_BYTE);

	else if( nct == NC_UBYTE )
		return(R_NC_TYPE_UBYTE);
	else if( nct == NC_USHORT )
		return(R_NC_TYPE_USHORT);
	else if( nct == NC_UINT )
		return(R_NC_TYPE_UINT);
	else if( nct == NC_INT64 )
		return(R_NC_TYPE_INT64);
	else if( nct == NC_UINT64 )
		return(R_NC_TYPE_UINT64);
	else if( nct == NC_STRING )
		return(R_NC_TYPE_STRING);
	else
		return(-1);
}

/*********************************************************************/
/* Returns a vector of dim sizes for the variable.
 * 'retval' is 0 for no error, or -1 for an error.
 */
void R_nc4_varsize( int *ncid, int *varid, int *ndims, int *varsize, int *retval )
{
	int 	i, err, dimid[NC_MAX_DIMS];
	size_t	dimlen;

	*retval = 0;

	/* Get ndims */
	err = nc_inq_varndims( *ncid, *varid, ndims );
	if( err != NC_NOERR ) {
		Rprintf( "Error in R_nc4_varsize on nc_inq_varndims call: %s\n", 
			nc_strerror(err) );
		*retval = -1;
		return;
		}

	/* Get dimids */
	err = nc_inq_vardimid( *ncid, *varid, dimid );
	if( err != NC_NOERR ) {
		Rprintf( "Error in R_nc4_varsize on nc_inq_vardimid call: %s\n",
			nc_strerror(err) );
		*retval = -1;
		return;
		}

	/* Get size of each dim */
	for( i=0; i<(*ndims); i++ ) {
		err = nc_inq_dimlen( *ncid, dimid[i], &dimlen );
		if( err != NC_NOERR ) {
			Rprintf( "Error in R_nc4_varsize on nc_inq_dimlen call: %s\n",
				nc_strerror(err) );
			*retval = -1;
			return;
			}
		varsize[i] = (int)dimlen;
		}
}

/*********************************************************************/
void R_nc4_inq_vartype( int *ncid, int *varid, int *precint, int *retval )
{
	nc_type nct;

	*retval = nc_inq_vartype( *ncid, *varid, &nct );

	if( *retval != NC_NOERR ) 
		Rprintf( "Error in R_nc4_inq_var: %s\n", 
			nc_strerror(*retval) );

	*precint = R_nc4_nctype_to_Rtypecode(nct);
}

/*********************************************************************/
void R_nc4_inq_varndims( int *ncid, int *varid, int *ndims, int *retval )
{
	*retval = nc_inq_varndims(*ncid, *varid, ndims );
	if( *retval != NC_NOERR ) 
		Rprintf( "Error in R_nc4_inq_varndims: %s\n", 
			nc_strerror(*retval) );
}

/*********************************************************************/
/* Inputs:
 *	sx_fixmiss  : is 1 if we want to fix the missing values in this
 *  		routine, and 0 otherwise
 * 	sx_imvstate : 0=var has no missing value; 1=var has a NA for
 * 		the missing value; 2=var has a valid, non-NA missing value
 *
 * Returns a list with elements:
 *	$error	: 0 for success, -1 for error
 *	$data   : array of integer values read in from the netcdf file
 */
SEXP Rsx_nc4_get_vara_double( SEXP sx_ncid, SEXP sx_varid, SEXP sx_start, SEXP sx_count, 
	SEXP sx_fixmiss, SEXP sx_imvstate, SEXP sx_missval )
{
	SEXP	sx_retval, sx_retnames, sx_reterr, sx_retdata;
	int	ncid, varid, i, err, ndims, fixmiss, imvstate, scalar_var;
	double	*p_data, missval, mvtol;
	size_t	s_start[MAX_NC_DIMS], s_count[MAX_NC_DIMS], tot_size;
	char	vn[2048];

	/* Make space for our returned list, which will have
	 * two elements, named $error and $data
	 */
	PROTECT( sx_retval = allocVector( VECSXP, 2 ));	/* 2 elements in returned list */

	/* Set names for the returned list */
	PROTECT( sx_retnames = allocVector( STRSXP, 2 ));
	SET_STRING_ELT( sx_retnames, 0, mkChar("error") );
	SET_STRING_ELT( sx_retnames, 1, mkChar("data" ) );
	setAttrib( sx_retval, R_NamesSymbol, sx_retnames );
	UNPROTECT(1);     				/* done with sx_retnames */

	/* Set provisional 'no error' retval */
	PROTECT(sx_reterr = allocVector( INTSXP, 1 ));
	INTEGER(sx_reterr)[0] = 0;

	ncid  	   	= INTEGER(sx_ncid    )[0];
	varid 	   	= INTEGER(sx_varid   )[0];
	fixmiss		= INTEGER(sx_fixmiss )[0];
	imvstate	= INTEGER(sx_imvstate)[0];
	missval		= REAL   (sx_missval )[0];

	/* Get number of dimensions in this variable */
	err = nc_inq_varndims( ncid, varid, &ndims );
	if( err != NC_NOERR ) {
		Rprintf( "Error in R_nc4_get_vara_double while getting ndims: %s\n", 
			nc_strerror(err) );
		INTEGER(sx_reterr)[0] = -1;
		SET_VECTOR_ELT( sx_retval, 0, sx_reterr );
		UNPROTECT(2);
		return( sx_retval );
		}

	/* Sanity check -- number of start and count elements must
	   match number of dims in the var
	*/
	scalar_var = ((ndims==0) && (GET_LENGTH(sx_start)==1) && (INTEGER(sx_start)[0]==0) && (INTEGER(sx_count)[0]==1));
	if( (!scalar_var) && (ndims != GET_LENGTH(sx_start))) {
		Rprintf( "Error in R_nc4_get_vara_double: I think var has %d dimensions, but passed start array is length %d. They must be the same!\n",
			ndims, GET_LENGTH(sx_start) );
		INTEGER(sx_reterr)[0] = -1;
		SET_VECTOR_ELT( sx_retval, 0, sx_reterr );
		UNPROTECT(2);
		return( sx_retval );
		}

	if( (!scalar_var) && (ndims != GET_LENGTH(sx_count))) {
		Rprintf( "Error in R_nc4_get_vara_double: I think var has %d dimensions, but passed count array is length %d. They must be the same!\n",
			ndims, GET_LENGTH(sx_count) );
		INTEGER(sx_reterr)[0] = -1;
		SET_VECTOR_ELT( sx_retval, 0, sx_reterr );
		UNPROTECT(2);
		return( sx_retval );
		}

	/* Get total number of elements we will be reading so we can 
	 * allocate R space
	 */
	tot_size = 1L;
	for( i=0; i<ndims; i++ ) {
		s_start[i] = (size_t)(INTEGER(sx_start)[i]);
		s_count[i] = (size_t)(INTEGER(sx_count)[i]);
		tot_size *= s_count[i];
		}
		
	/* Make space for the returned data */
	PROTECT( sx_retdata = allocVector(REALSXP, tot_size));
	p_data = REAL( sx_retdata );

	/* Actually read in the data now */
	err = nc_get_vara_double( ncid, varid, s_start, s_count, p_data );
	if( err != NC_NOERR ) {
		nc_inq_varname( ncid, varid, vn );
		Rprintf( "Error in Rsx_nc4_get_vara_double: %s\n", 
			nc_strerror( err ) );
		Rprintf( "Var: %s  Ndims: %d   Start: ", vn, ndims );
		for( i=0; i<ndims; i++ ) {
			Rprintf( "%u", (unsigned int)s_start[i] );
			if( i < ndims-1 )
				Rprintf( "," );
			}
		Rprintf( " " );
		Rprintf( "Count: " );
		for( i=0; i<ndims; i++ ) {
			Rprintf( "%u", (unsigned int)s_count[i] );
			if( i < ndims-1 )
				Rprintf( "," );
			}
		Rprintf( "\n" );
		INTEGER(sx_reterr)[0] = -1;
		SET_VECTOR_ELT( sx_retval, 0, sx_reterr );
		UNPROTECT(3);
		return( sx_retval );
		}

	/* We only do this for imvstate==2, which is a valid, non-NA
	 * missing value. If the variable has NO missing value
	 * (imvstate==0) or a NA missing value (imvstate==1) then
	 * we skip this block.
	 */
	if( (fixmiss == 1) && (imvstate == 2)) {
		if( missval == 0.0 ) 
			mvtol = 1.e-10;	/* arbitrary -- what are you supposed to use?? */
		else
			mvtol = fabs( missval ) * 1.e-5;

		for( i=0L; i<tot_size; i++ ) {
			if( fabs( p_data[i] - missval ) < mvtol )
				p_data[i] = NA_REAL;
			}
		}

	/* OK we can successfully return now */
	SET_VECTOR_ELT( sx_retval, 0, sx_reterr  );
	SET_VECTOR_ELT( sx_retval, 1, sx_retdata );
	UNPROTECT(3);
	return( sx_retval );
}

/*********************************************************************/
/* Input value byte_style is 1 for signed, 2 for unsigned
 * 
 * Returns a list with elements:
 *	$error	: 0 for success, -1 for error
 *	$data   : array of integer values read in from the netcdf file
 */
SEXP Rsx_nc4_get_vara_int( SEXP sx_ncid, SEXP sx_varid, SEXP sx_start, 
	SEXP sx_count, SEXP sx_byte_style )
{
	SEXP	sx_retval, sx_retnames, sx_reterr, sx_retdata;
	int	ncid, varid, byte_style, i, err, ndims, *p_data, scalar_var;
	size_t	s_start[MAX_NC_DIMS], s_count[MAX_NC_DIMS], tot_size, k;
	char	vn[2048];
	nc_type	nct;

	/* Make space for our returned list, which will have
	 * two elements, named $error and $data
	 */
	PROTECT( sx_retval = allocVector( VECSXP, 2 ));	/* 2 elements in returned list */

	/* Set names for the returned list */
	PROTECT( sx_retnames = allocVector( STRSXP, 2 ));
	SET_STRING_ELT( sx_retnames, 0, mkChar("error") );
	SET_STRING_ELT( sx_retnames, 1, mkChar("data" ) );
	setAttrib( sx_retval, R_NamesSymbol, sx_retnames );
	UNPROTECT(1);     				/* done with sx_retnames */

	/* Set provisional 'no error' retval */
	PROTECT(sx_reterr = allocVector( INTSXP, 1 ));
	INTEGER(sx_reterr)[0] = 0;

	ncid  	   = INTEGER(sx_ncid      )[0];
	varid 	   = INTEGER(sx_varid     )[0];
	byte_style = INTEGER(sx_byte_style)[0];

	/* Get number of dimensions in this variable */
	err = nc_inq_varndims( ncid, varid, &ndims );
	if( err != NC_NOERR ) {
		Rprintf( "Error in R_nc4_get_vara_int while getting ndims: %s\n", 
			nc_strerror(err) );
		INTEGER(sx_reterr)[0] = -1;
		SET_VECTOR_ELT( sx_retval, 0, sx_reterr );
		UNPROTECT(2);
		return( sx_retval );
		}

	/* Sanity check -- number of start and count elements must
	   match number of dims in the var
	*/
	scalar_var = ((ndims==0) && (GET_LENGTH(sx_start)==1) && (INTEGER(sx_start)[0]==0) && (INTEGER(sx_count)[0]==1));
	if( (!scalar_var) && (ndims != GET_LENGTH(sx_start))) {
		Rprintf( "Error in R_nc4_get_vara_int: I think var has %d dimensions, but passed start array is length %d. They must be the same!\n",
			ndims, GET_LENGTH(sx_start) );
		INTEGER(sx_reterr)[0] = -1;
		SET_VECTOR_ELT( sx_retval, 0, sx_reterr );
		UNPROTECT(2);
		return( sx_retval );
		}

	if( (!scalar_var) && (ndims != GET_LENGTH(sx_count))) {
		Rprintf( "Error in R_nc4_get_vara_int: I think var has %d dimensions, but passed count array is length %d. They must be the same!\n",
			ndims, GET_LENGTH(sx_count) );
		INTEGER(sx_reterr)[0] = -1;
		SET_VECTOR_ELT( sx_retval, 0, sx_reterr );
		UNPROTECT(2);
		return( sx_retval );
		}

	/* Get total number of elements we will be reading so we can 
	 * allocate R space
	 */
	tot_size = 1L;
	for( i=0; i<ndims; i++ ) {
		s_start[i] = (size_t)(INTEGER(sx_start)[i]);
		s_count[i] = (size_t)(INTEGER(sx_count)[i]);
		tot_size *= s_count[i];
		}
		
	/* Make space for the returned data */
	PROTECT( sx_retdata = allocVector(INTSXP, tot_size));
	p_data = INTEGER( sx_retdata );

	/* Actually read in the data now */
	err = nc_get_vara_int( ncid, varid, s_start, s_count, p_data );
	if( err != NC_NOERR ) {
		nc_inq_varname( ncid, varid, vn );
		Rprintf( "Error in Rsx_nc4_get_vara_int: %s\n", 
			nc_strerror( err ) );
		Rprintf( "Var: %s  Ndims: %d   Start: ", vn, ndims );
		for( i=0; i<ndims; i++ ) {
			Rprintf( "%u", (unsigned int)s_start[i] );
			if( i < ndims-1 )
				Rprintf( "," );
			}
		Rprintf( " Count: " );
		for( i=0; i<ndims; i++ ) {
			Rprintf( "%u", (unsigned int)s_count[i] );
			if( i < ndims-1 )
				Rprintf( "," );
			}
		Rprintf( "\n" );
		INTEGER(sx_reterr)[0] = -1;
		SET_VECTOR_ELT( sx_retval, 0, sx_reterr );
		UNPROTECT(3);
		return( sx_retval );
		}

	err = nc_inq_vartype( ncid, varid, &nct );
	if( nct == NC_BYTE ) {
		/* netcdf library for reading from byte to int, as we do here,
		 * is SIGNED.  So, if user requests signed, we don't have to
		 * do anything; only adjust if user asks for unsigned.
		 */
		if( byte_style == 2 ) {	/* unsigned */
			for( k=0L; k<tot_size; k++ ) {
				if( p_data[k] < 0 )
					p_data[k] += 256;
				}
			}
		}

	/* OK we can successfully return now */
	SET_VECTOR_ELT( sx_retval, 0, sx_reterr  );
	SET_VECTOR_ELT( sx_retval, 1, sx_retdata );
	UNPROTECT(3);
	return( sx_retval );
}

/*********************************************************************/
void R_nc4_get_vara_text( int *ncid, int *varid, int *start, 
	int *count, char **tempstore, char **data, int *retval )
{
	int	i, err, ndims;
	size_t	s_start[MAX_NC_DIMS], s_count[MAX_NC_DIMS], nstr, slen;
	char	vn[2048];

	err = nc_inq_varndims(*ncid, *varid, &ndims );
	if( err != NC_NOERR ) 
		Rprintf( "Error in R_nc4_get_vara_text while getting ndims: %s\n", 
			nc_strerror(*retval) );

	nstr = 1L;
	for( i=0; i<ndims; i++ ) {
		s_start[i] = (size_t)start[i];
		s_count[i] = (size_t)count[i];
		if( i < (ndims-1) ) 
			nstr *= s_count[i];
		}
	slen = s_count[ndims-1];

	*retval = nc_get_vara_text(*ncid, *varid, s_start, s_count, tempstore[0] );

	if( *retval != NC_NOERR ) {
		nc_inq_varname( *ncid, *varid, vn );
		Rprintf( "Error in R_nc4_get_vara_text: %s\n", 
			nc_strerror(*retval) );
		Rprintf( "Var: %s  Ndims: %d   Start: ", vn, ndims );
		for( i=0; i<ndims; i++ ) {
			Rprintf( "%u", (unsigned int)s_start[i] );
			if( i < ndims-1 )
				Rprintf( "," );
			}
		Rprintf( " Count: " );
		for( i=0; i<ndims; i++ ) {
			Rprintf( "%u", (unsigned int)s_count[i] );
			if( i < ndims-1 )
				Rprintf( "," );
			}
		Rprintf( "\n" );
		}

	/* Now copy each string over to the final data array */
	for( i=0; i<nstr; i++ ) {
		strncpy( data[i], tempstore[0]+i*slen, slen );
		data[i][slen] = '\0';
		}
}

/*********************************************************************/
/* Code kindly donated by Martin Morgan of the Fred Hutchinson
 * Cancer Research Center
 */
SEXP R_nc4_blankstring(SEXP size)
{
	int	i, len;
	char	*str;
	SEXP	blank, string;

	if (!IS_INTEGER(size) || LENGTH(size)!=1)
		Rf_error("'size' must be 'integer(1)'");

	len = INTEGER(size)[0];
	if (len < 0)
		Rf_error("'size' must be >= 0");

	str = R_alloc(len + 1, sizeof(char));
	for (i = 0; i < len; ++i)
		str[i] = ' ';
	str[len] = '\0';

	blank  = PROTECT(NEW_CHARACTER(1));
	string = PROTECT(mkChar(str));
	SET_STRING_ELT(blank, 0, string);
	UNPROTECT(2);

	return(blank); 
}

/*********************************************************************************
 * Read vlen strings given the numeric varid, start, and count to use
 */
SEXP R_nc4_get_vara_string( SEXP sx_nc, SEXP sx_varid, SEXP sx_start, SEXP sx_count ) 
{
	SEXP	sx_retval, sx_retnames, sx_retstrings, sx_reterror;
	int	i, ierr, nchar, varid, ncid, ndims, count_int[MAX_NC_DIMS], start_int[MAX_NC_DIMS], len_count, len_start; 
	size_t	count[MAX_NC_DIMS], start[MAX_NC_DIMS], tot_count, isz;
	char 	**ss;

	/* Convert passed parameters (which are in R format) into C format */
	ncid  = INTEGER(sx_nc   )[0];
	varid = INTEGER(sx_varid)[0];

	len_start = length(sx_start);
	for( i=0; i<len_start; i++ ) {
		start_int[i] = INTEGER(sx_start)[i];
		start    [i] = (size_t)(start_int[i]);
		}

	len_count = length(sx_count);
	for( i=0; i<len_count; i++ ) {
		count_int[i] = INTEGER(sx_count)[i];
		count    [i] = (size_t)(count_int[i]);
		}

	PROTECT( sx_retval   = allocVector( VECSXP, 2 ));       /* 2 elements in the returned list: $error, $strings */

	/* Set the names for the returned list */
	PROTECT( sx_retnames = allocVector( STRSXP, 2 ));       /* 2 elements in the returned list */
	SET_STRING_ELT( sx_retnames, 0, mkChar("error") );
	SET_STRING_ELT( sx_retnames, 1, mkChar("data") );
	setAttrib( sx_retval, R_NamesSymbol, sx_retnames );
	UNPROTECT(1); 

	/* Set the return error code */
	PROTECT( sx_reterror = allocVector( INTSXP, 1 ));
	INTEGER( sx_reterror)[0] = 0;

	/* Get number of dims in the var */
	ierr = nc_inq_varndims( ncid, varid, &ndims );

	/*--------------------------------------------------------------
	 * At this point we have all the C values we need:
	 *  	ncid
	 *	varid (numeric)
	 *	ndims
	 *	start
	 *	count
	 *---------------------------------------------------------------*/
	tot_count = 1L;
	for( i=0; i<ndims; i++ ) 
	 	tot_count *= count[i];

	ss = (char **)malloc( sizeof( char *) * tot_count );
	if( ss == NULL ) {
		INTEGER( sx_reterror)[0] = -1;
		error("ncdf4 library: routine R_nc4_get_vara_string: Error trying to allocate space to read the vlen strings: total count of strings requested: %ld\n", tot_count );
		}

	if( (ierr = nc_get_vara_string( ncid, varid, start, count, ss )) != 0 ) {
		INTEGER( sx_reterror)[0] = -2;
		error("ncdf4 library: routine R_nc4_get_vara_string: Error reading vlen strings: %s\n",
			nc_strerror(ierr));
		}

	PROTECT( sx_retstrings = allocVector( STRSXP, tot_count ));
	for( isz=0L; isz<tot_count; isz++ ) {
		nchar = strlen( ss[isz] );
		SET_STRING_ELT( sx_retstrings, isz, mkChar(ss[isz]) );
		}

	SET_VECTOR_ELT( sx_retval, 0, sx_reterror   );
	SET_VECTOR_ELT( sx_retval, 1, sx_retstrings );

	UNPROTECT(3);	

	/* Free netcdf string storage */
	nc_free_string( tot_count, ss );

	return( sx_retval );
}

/*********************************************************************/
/* cmode is 0 for read only, 1 for write access.
 */
void R_nc4_open( char **filename, int *cmode, int *ncid, int *retval )
{
	int	nc_mode;

	if( *cmode == 0 )
		nc_mode = 0;
	else if( *cmode == 1 )
		nc_mode = NC_WRITE;
	else
		{
		Rprintf( "Error in R_nc4_open: bad mode passed.  Must be 0 (read) or 1 (write)\n");
		*retval = -1;
		return;
		}

	*retval = nc_open(R_ExpandFileName(filename[0]), nc_mode, ncid);
	if( *retval != NC_NOERR ) 
		Rprintf( "Error in R_nc4_open: %s\n", 
			nc_strerror(*retval) );
}

/*********************************************************************/
void R_nc4_close( int *ncid )
{
	int	err;
	err = nc_close(*ncid);
	if( err != NC_NOERR ) 
		Rprintf( "Error in R_nc4_close: %s\n", 
			nc_strerror(err) );
}

/************************************************************************************************/
/* Returned index_first_slash is 0-based counting 
 */
int R_nc4_util_nslashes( char *s, int *idx_first_slash ) 
{
	int	i, nslashes;

	nslashes = 0;
	*idx_first_slash = -1;
	for( i=0; i<strlen(s); i++ ) {
		if( s[i] == '/' ) {
			nslashes++;
			if( *idx_first_slash == -1 ) 
				*idx_first_slash = i;
			}
		}

	return( nslashes );
}

/****************************************************************************************/
void R_nc4_inq_varid_hier_inner( int *ncid, char *varname, int *returned_grpid, int *returned_varid )
{
	int ierr, nslashes, idx_first_slash, gid;
	char gname[MAX_NC_NAME];

/* Rprintf("R_nc4_inq_varid_hier_inner: Entering with varname=>%s<\n", varname ); */

	/* Passed var name must not start with a slash */
	if( varname[0] == '/' ) {
		Rprintf( "Error in R_nc4_varid_hier: passed varname must not start with a slash!\n" );
		*returned_varid = -1;
		*returned_grpid = -1;
		return;
		}

	/* If there are no forward slashes in the name, then it is a simple var name */
	nslashes = R_nc4_util_nslashes( varname, &idx_first_slash );
/* Rprintf("R_nc4_inq_varid_hier_inner: varname >%s< has nslashes=%d and idx first slash=%d\n", varname, nslashes, idx_first_slash ); */
	if( nslashes == 0 ) {
		*returned_grpid = *ncid;
		ierr = nc_inq_varid(*ncid, varname, returned_varid );
		if( ierr != NC_NOERR ) {
			*returned_varid = -1;
			*returned_grpid = -1;
			}
/* Rprintf("R_nc4_inq_varid_hier_inner: Found id for var >%s< varid:%d gid:%d\n", varname, *returned_varid, *returned_grpid ); */
		return;
		}

	/* If we get here then the varname must have at least one slash, which indicates
	 * a group name under the current group name. Get the ID of that group
	 */
	strncpy( gname, varname, idx_first_slash );
	gname[ idx_first_slash ] = '\0';
	ierr = nc_inq_grp_ncid( *ncid, gname, &gid );
	if( ierr != NC_NOERR ) {
		Rprintf( "Error in R_nc4_varid_hier: extracted groupname not found: >%s< (coding error?)\n", gname );
		*returned_varid = -1;
		*returned_grpid = -1;
		return;
		}

	/* Call ourselves recursively with the subgroup id and the varname stripped
	 * of its leading group name
	 */
/* Rprintf("R_nc4_inq_varid_hier_inner: calling myself recursively with varname=>%s<\n", varname+idx_first_slash+1 ); */
	R_nc4_inq_varid_hier_inner( &gid, varname+idx_first_slash+1, returned_grpid, returned_varid );
}

/****************************************************************************************/
/* This version of inq_varid works with hierarchial var names such as "model1/TS" and
 * returns both the group ID and the var ID in that group, given the root group id
 * (netcdf file ID).
 *
 * If the passed var is not found in the file, then -1 is returned.
 *
 * The passed slash-ful varname must not start or end with a slash.
 */
void R_nc4_inq_varid_hier( int *ncid, char **varname, int *returned_grpid, int *returned_varid )
{
/* Rprintf("R_nc4_inq_varid_hier: entering for var >%s<\n", varname[0] ); */
	R_nc4_inq_varid_hier_inner( ncid, varname[0], returned_grpid, returned_varid );
}