/*! $Id: elio.c,v 1.60 2012/02/07 19:19:01 pturner Exp $
 *
 * I/O routines for new file format -  pturner 11-2002
 *
 */

/*! \mainpage ELIO - Interface to ELCIRC data files.
 **
 ** \section intro Introduction
 **
 ** The ELIO library provides a C-programming language interface to output files generated by ELCIRC. The CVS repository is 'ace/elio' - use:
<pre>
    cvs checkout ace/elio
</pre>
to get the elio source code as well as the sources for other utilities using the elio library. These codes are documented at the URL:
<pre>
    <a href="http://amb1019.ccalmr.ogi.edu/doc/codes">http://amb1019.ccalmr.ogi.edu/doc/codes</a>
</pre>
 **
 ** \section use Building your programs.
 **
 ** Makefile
<pre>
# 
# Sample Makefile for use of libelio.a.
#
# \$Revision\$
# \$Log\$
#
#

INSTALLDIR = /usr/local/ace

OBJS = testel.o

# -miee for alphas
# CFLAGS = -mieee
# for intel
CFLAGS = -D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64
LIB = /usr/local/ace/lib
CC = gcc -g

testel: testel.o 
	\$(CC) \$(CFLAGS) -g testel.c -L\$(LIB) -lelio -lm -o testel
</pre>

 ** 
 ** \section examp Example of use.
 **
 ** The following code reads the header from a data file and prints out a
 ** summary of what is seen using the ElioPrintHeader() function.
 ** 
<pre>

// ginfo.c example program - pturner.
//
// \$Revision\$
//
// \$Log\$
//

\#include <stdio.h>
\#include <stdlib.h>
\#include <math.h>
\#include "elio.h"

int main(int argc, char **argv)
{
    char iname[1024];
    FILE *fp;
    ElcircHeader h;
    int err;

    if (argc != 2) {
        fprintf(stderr, "Usage: \%s <ELCIRC *.6? file>\\n", argv[0]);
        exit(1);
    }

    strcpy(iname, argv[1]);
// check that file type is correct
    if ((fp = fopen(iname, "rb")) == NULL) {
        fprintf(stderr, "\%s: Unable to open file \%s.\\n", argv[0], iname);
        exit(1);
    }
    if (ElioGetFileType(fp) < 2) {
        fprintf(stderr, "\%s: Incorrect file type for file \%s.\\n", argv[0], iname);
        exit(1);
    }
    fclose(fp);
// Get the header
    if ((err = ElioGetHeader(iname, &h)) != ELIO_OK) {
        fprintf(stderr, "\%s: Error in call to GetElcircHeader(), error # \%d.\\n", argv[0], err);
        exit(1);
    }
    ElioPrintHeader(&h);
    printf("Number of steps in file: \%d\\n", ElioGetNStepsInFile(iname, &h));
    ElioFreeHeader(&h);
}
</pre>

Makefile for ginfo.c

<pre>
# 
# Sample Makefile for use of libelio.a.
#
# \$Revision\$
# \$Log\$
#
# -miee for alphas
# CFLAGS = -mieee
# for intel
CFLAGS = -D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64
LIB = /usr/local/ace/lib
CC = gcc -g

ginfo: ginfo.o 
	\$(CC) \$(CFLAGS) -g ginfo.c -L\$(LIB) -lelio -lm -o ginfo
</pre>


 ** 
 **/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <time.h>
#include "elio.h"

/*!
 * @def WIN32
 *
 * Windows differences, if true, fix things up for cygwin and Visual C++.
 *
 */
#ifdef WIN32
#include <string.h>
#define off_t long
#define fseeko fseek
#define ftello ftell
#endif

/*!
 * Error codes, line numbers, not used right now.
 */
int ElioErr = 0;
int ElioSysErr = 0;
int ElioLineErr = 0;

/*!
 * Free memory allocated by a previous call to ElioGetHeader()
 * Calling this on an unallocated Header will likely result 
 * in a seg fault.
 */
void ElioFreeHeader(ElcircHeader * h)
{
    free(h->zcor);
    free(h->x);
    free(h->y);
    free(h->d);
    free(h->etype);
    free(h->icon[0]);
    free(h->icon[1]);
    free(h->icon[2]);
    if (h->v == 3) {
	free(h->icon[3]);
    }
    if (h->v != 4) {
	free(h->bi);
	free(h->bisave);
    }
    free(h->no);
}

/*!
 * Free a previously allocated time step.
 */
void ElioFreeTimeStep(ElcircTimeStep * t)
{
    if (t->e != NULL) {
	free(t->e);
    }
    if (t->surfind != NULL) {
	free(t->surfind);
    }
    free(t->d);
}

/*!
 * Allocate memory for a time step given information in a valid 
 * ElcircHeader.
 *
 * Returns ELIO_ERR on failure after cleaning up any current allocations.
 * Returns ELIO_OK on success.
 */
int ElioAllocateTimeStep(ElcircHeader * h, ElcircTimeStep * t)
{
    switch (h->v) {
    case 2:
    case 3:
	t->surfind = (int *) malloc(h->np * sizeof(int));
	t->d = (float *) malloc(h->nitems * sizeof(float));
	t->e = (float *) NULL;
	if (t->surfind == NULL || t->d == NULL) {
	    return ELIO_ERR;
	}
	break;
    case 4:
    case 5:
	t->surfind = (int *) malloc(h->np * sizeof(int));
	t->e = (float *) malloc(h->np * sizeof(float));
	t->d = (float *) malloc(h->nitems * sizeof(float));
	if (t->e == NULL || t->d == NULL) {
	    return ELIO_ERR;
	}
	break;
    }
    return ELIO_OK;
}

/*!
 * Print header information from a valid ElcircHeader struct.
 *
 * \param h = Previously filled in ElcircHeader.
 * 
 */
void ElioPrintHeader(ElcircHeader * h)
{
    int i;
    printf("Magic: %s\n", h->magic);
    printf("Version: %s\n", h->version);
    printf("Start Time: %s\n", h->start_time);
    printf("Variable Name: %s\n", h->variable_nm);
    printf("Variable Dim: %s\n", h->variable_dim);
    printf("Number of Steps: %d\n", h->nsteps);
    printf("Time Step: %f\n", h->timestep);
    printf("Skip Steps: %d\n", h->skip);
    printf("Variable Type: %d\n", h->ivs);
    printf("Variable Dimension: %d\n", h->i23d);
    switch (h->v) {
    case 2:
    case 3:
	printf("Data is Z Level\n");
	printf("Vertical position: %f\n", h->vpos);
	printf("Z of MSL: %f\n", h->zmsl);
	printf("Number of Vertical Levels: %d\n", h->nvrt);
	printf("Levels:\n");
	for (i = 0; i < h->nvrt; i++) {
	    printf("    %02d: %f\n", i, h->zcor[i]);
	}
	break;
    case 4:
	printf("Data is Sigma Level\n");
	printf("ivcor: %d\n", h->ivcor);
	printf("h0: %f\n", h->h0);
	printf("hc: %f\n", h->hc);
	printf("Theta b: %f\n", h->thetab);
	printf("Theta f: %f\n", h->thetaf);
	printf("Number of Sigma Levels: %d\n", h->nvrt);
	printf("Levels:\n");
	for (i = 0; i < h->nvrt; i++) {
	    printf("    %02d: %f\n", i, h->zcor[i]);
	}
	break;
    case 5:
	printf("Data is Sigma+Z\n");
	printf("Total number of Levels: %d\n", h->nvrt);
	printf("Sigma levels ks: %d\n", h->ks);
	printf("Z levels kz: %d\n", h->kz);
	printf("h0: %f\n", h->h0);
	printf("hs: %f\n", h->hs);
	printf("hc: %f\n", h->hc);
	printf("Theta b: %f\n", h->thetab);
	printf("Theta f: %f\n", h->thetaf);
	printf("Levels:\n");
	for (i = 0; i < h->nvrt; i++) {
	    printf("    %02d: %f\n", i, h->zcor[i]);
	}
	break;
    default:
	printf("Vertical position: %f\n", h->vpos);
	printf("Z of MSL: %f\n", h->zmsl);
	printf("Number of Vertical Levels: %d\n", h->nvrt);
	printf("Levels:\n");
	for (i = 0; i < h->nvrt; i++) {
	    printf("    %02d: %f\n", i, h->zcor[i]);
	}
	break;
    }
    printf("Number of Floats/Step: %d\n", h->nitems);
    printf("Header Size (bytes): %d\n", h->hsize);
    printf("Time Step Size (bytes): %d\n", h->ssize);
    printf("Number of Nodes: %d\n", h->np);
    printf("Number of Elements: %d\n", h->ne);
    printf("Internal version number: %d\n", h->v);
    printf("First Node: 1 -> %f, %f, %f\n", h->x[0], h->y[0], h->d[0]);
    printf("First Element: 1 -> %d, %d, %d\n", h->icon[0][0], h->icon[1][0], h->icon[2][0]);
    printf("Last Node: %d -> %f, %f, %f\n", h->np, h->x[h->np - 1], h->y[h->np - 1], h->d[h->np - 1]);
    printf("Last Element: %d -> %d, %d, %d\n", h->ne, h->icon[0][h->ne - 1], h->icon[1][h->ne - 1], h->icon[2][h->ne - 1]);
}

/*!
 * ElioGetTimeStep: Read data for a single time step. 
 *
 * Input variables:
 *
 * \param fp = Open file handle.
 * \param step = Time step to read starting from 0 to nsteps - 1.
 * \param h = Valid ElcircHeader struct.
 *
 * Output Variables:
 *
 * \param t = Valid ElcircTimeStep.
 *
 * On success returns ELIO_OK.
 *
 * On failure, returns ELIO_FSEEK_ERR or ELIO_FREAD_ERR depending
 * on the nature of the error.
 */
int ElioGetTimeStep(FILE * fp, int step, ElcircHeader * h, ElcircTimeStep * t)
{
    off_t skip;
    int i, n;
    short *b;
    skip = (off_t) h->hsize + (off_t) step *(off_t) h->ssize;
    if (fseeko(fp, skip, SEEK_SET)) {
	return ELIO_FSEEK_ERR;
    }
    if (fread(&t->t, sizeof(float), 1, fp) != 1)
	return ELIO_FREAD_ERR;
    if (fread(&t->it, sizeof(int), 1, fp) != 1)
	return ELIO_FREAD_ERR;
    if (h->v == 4 || h->v == 5) {
	if (fread(t->e, sizeof(float), h->np, fp) != h->np)
	    return ELIO_FREAD_ERR;
	for (i = 0; i < h->np; i++) {
	    t->surfind[i] = h->nvrt;
	}
    } else {
	if (fread(t->surfind, sizeof(int), h->np, fp) != h->np)
	    return ELIO_FREAD_ERR;
	t->e = (float *) NULL;
    }
    if ((n = fread(t->d, sizeof(float), h->nitems, fp)) != h->nitems) {
	return ELIO_FREAD_ERR;
    }
/*
    printf("Completed Read of Time Step (%d): Time, Step = %f, %d\n", step, t->t, t->it);
*/
    return ELIO_OK;
}

/*!
 * Returns the Z positions of the vertical given the elevation, e. Also
 * the bottom and surface indexes in i1 and i2. Array Zpos1 is returned 
 * with the valid z values. Only values from zpos1[i1] through zpos[i2]
 * contain valid data.
 */
int ElioGetZPos(ElcircHeader * h, double d, double e, int *i1, int *i2, double *zpos1)
{
    int i, ilev1 = 0, ilev2 = 0;
    for (i = 0; i < h->nvrt; i++) {
	zpos1[i] = -h->zmsl + h->zcor[i];
    }

    ilev1 = 0;
    while (-d >= zpos1[ilev1] && ilev1 < h->nvrt) {
	ilev1++;
    }
    zpos1[ilev1] = (-d);
    ilev2 = h->nvrt - 1;
    while (e <= zpos1[ilev2] && ilev2 > 0) {
	ilev2--;
    }
    ilev2++;
    zpos1[ilev2] = e;
    *i1 = ilev1;
    *i2 = ilev2;
    return ELIO_OK;
}

/*!
 * Return the surface index for a node.
 * 
 * Input variables:
 *
 * \param fp = Open file handle
 * \param node = node number starting from 0.
 * \param step = Time step index starting from 0.
 * \param h = Valid ElcircHeader struct.
 */
int ElioGetNodeSurfaceIndex(FILE * fp, int node, int step, ElcircHeader * h)
{
    off_t skip;
    int retval = -1;
    int ti;
    if (node >= h->np || node < 0) {
	return retval;
    }
    skip = (off_t) h->hsize + (off_t) step *(off_t) h->ssize + (off_t) 8 + (off_t) (node * 4);
    if (fseeko(fp, skip, SEEK_SET)) {
	return ELIO_FSEEK_ERR;
    }
    if (fread(&ti, sizeof(int), 1, fp) != 1) {
	return ELIO_FREAD_ERR;
    }
    return ti;
}

/*!
 * Return data for a single node from one time step.
 * 
 * Input variables:
 *
 * \param fp = Open file handle.
 * \param step = Time step to read starting from 0.
 * \param node = Node to read starting from 0.
 * \param h = Valid ElcircHeader struct.
 *
 * Output Variables:
 *
 * \param t = Pointer to the time for this time step, float.
 * \param it = Pointer to the current computational time step, integer.
 * \param bind = Pointer to the bottom index starting from 0.
 * \param sind = Pointer to the surface index starting from 0.
 * \param d = Pointer to an array of float values, with valid data from 
 *     d[bind] through d[sind].
 *
 * Returns ELIO_OK on success.
 * Returns ELIO_FSEEK_ERR or ELIO_FREAD_ERR on error.
 *
 */
int ElioGetNode(FILE * fp, int step, int node, ElcircHeader * h, float *t, int *it, int *bind, int *sind, float *d)
{
    int i, npts, itmp;
    float e;
    off_t skip;
    *bind = h->bi[node];
    if (h->i23d == 3) {		/* 3d data */
	skip = (off_t) h->hsize + (off_t) step *(off_t) h->ssize;
	if (fseeko(fp, skip, SEEK_SET)) {	/* Read time, it, and sind */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(t, sizeof(float), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	if (fread(it, sizeof(int), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	if (h->v == 4) {
	    itmp = h->nvrt;
	    *bind = 0;
	    if (fseeko(fp, node * 4L, SEEK_CUR)) {	/* elevation */
		return ELIO_FSEEK_ERR;
	    }
	    if (fread(&e, sizeof(float), 1, fp) != 1) {
		return ELIO_FREAD_ERR;
	    }
	} else if (h->v == 5) {
	    itmp = h->nvrt;
	    if (fseeko(fp, node * 4L, SEEK_CUR)) {	/* elevation */
		return ELIO_FSEEK_ERR;
	    }
	    if (fread(&e, sizeof(float), 1, fp) != 1) {
		return ELIO_FREAD_ERR;
	    }
	} else {
	    if (fseeko(fp, node * 4L, SEEK_CUR)) {	/* sind */
		return ELIO_FSEEK_ERR;
	    }
	    if (fread(&itmp, sizeof(int), 1, fp) != 1) {
		return ELIO_FREAD_ERR;
	    }
	}
	*sind = itmp - 1;
	npts = (h->nvrt - h->bi[node]) * h->ivs;	/* number of points to read for this node */
	if (fseeko(fp, skip + (off_t) ((2 + h->np + h->no[node]) * 4), SEEK_SET)) {	/* skip time, it, and the surf indexes */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(d, sizeof(float), npts, fp) != npts) {
	    return ELIO_FREAD_ERR;
	}
    } else if (h->i23d == 2) {	/* 2d data - just one or 2 values to read */
	skip = (off_t) h->hsize + (off_t) step *(off_t) h->ssize;
	if (fseeko(fp, skip, SEEK_SET)) {	/* Read time, it, and sind */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(t, sizeof(float), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	if (fread(it, sizeof(int), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	if (h->v == 4) {
	    itmp = h->nvrt;
	    *bind = 0;
	    if (fseeko(fp, node * 4L, SEEK_CUR)) {	/* elevation */
		return ELIO_FSEEK_ERR;
	    }
	    if (fread(&e, sizeof(float), 1, fp) != 1) {
		return ELIO_FREAD_ERR;
	    }
	} else if (h->v == 5) {
	    itmp = h->nvrt;
	    if (fseeko(fp, node * 4L, SEEK_CUR)) {	/* elevation */
		return ELIO_FSEEK_ERR;
	    }
	    if (fread(&e, sizeof(float), 1, fp) != 1) {
		return ELIO_FREAD_ERR;
	    }
	} else {
	    if (fseeko(fp, node * 4L, SEEK_CUR)) {	/* Read time, it, and sind */
		return ELIO_FSEEK_ERR;
	    }
	    if (fread(&itmp, sizeof(int), 1, fp) != 1) {
		return ELIO_FREAD_ERR;
	    }
	}
	*sind = itmp - 1;
	npts = h->ivs;		/* number of points to read for this node */
	if (fseeko(fp, skip + (off_t) (2L + h->np) * 4L + (off_t) (npts * node * 4L), SEEK_SET)) {	/* skip time, it, and the surf indexes */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(d, sizeof(float), npts, fp) != npts) {
	    return ELIO_FREAD_ERR;
	}
    } else {
	return ELIO_ERR;
    }
    return ELIO_OK;
}

/*!
 * Extract data for a single node from one previously time step.
 * Works for 2d and 3d data, vector and scalar.
 * 
 * Input variables:
 *
 * \param t = Pointer to previously read time step.
 * \param h = Pointer to valid ElcircHeader struct.
 * \param node = Node number starting from 0.
 * \param bind = Bottom index starting from 0.
 * \param sind = Surface index starting from 0.
 *
 * Output Variables:
 *
 * \param d = Pointer to array extracted from timestep struct.
 *
 * Returns ELIO_OK on success.
 *
 */
int ElioExtractNode(ElcircTimeStep * t, ElcircHeader * h, int node, int bind, int sind, float *d)
{
    int i;
    if (h->i23d == 2) {
	if (h->ivs == 2) {
	    d[0] = t->d[h->no[node]];
	    d[1] = t->d[h->no[node] + 1];
	} else {
	    d[0] = t->d[h->no[node]];
	}
    } else {
	for (i = bind; i <= sind; i++) {
	    if (h->ivs == 2) {
		d[2 * i] = t->d[h->no[node] + 2 * i - 2 * bind];
		d[2 * i + 1] = t->d[h->no[node] + 2 * i + 1 - 2 * bind];
	    } else {
		d[i] = t->d[h->no[node] + i - bind];
	    }
	}
    }
    return ELIO_OK;
}

/*!
 * Get data for a single node from one time step - fills in missing data 
 * values to match the old format.
 *
 * Input variables:
 *
 * \param fp = Pointer to open file handle.
 * \param step = Time step starting from 0.
 * \param node = Node number starting from 0.
 * \param h = Pointer to valid ElcircHeader struct.
 *
 * Output Variables:
 *
 * \param t = Pointer to time, float.
 * \param it = Pointer to computational time step, integer.
 * \param d = Pointer to array of return values, h.nvrt long.
 *
 * Returns ELIO_OK on success.
 * Returns ELIO_FSEEK_ERR and ELIO_FREAD_ERR on error.
 */
int ElioGetNodeOld(FILE * fp, int step, int node, ElcircHeader * h, float *t, int *it, float *d)
{
    int i, npts, bind, sind, rval = 0;
    off_t skip;
    float *tmp;
    tmp = (float *) malloc(h->nvrt * h->ivs * sizeof(float));
    bind = h->bi[node];
    if (h->i23d == 3) {		/* 3d data */
	skip = (off_t) h->hsize + (off_t) step *(off_t) h->ssize;
	if (fseeko(fp, skip, SEEK_SET)) {	/* Read time, it, and sind */
	    rval = ELIO_FSEEK_ERR;
	    goto out;
	}
	if (fread(t, sizeof(float), 1, fp) != 1) {
	    rval = ELIO_FREAD_ERR;
	    goto out;
	}
	if (fread(it, sizeof(int), 1, fp) != 1) {
	    rval = ELIO_FREAD_ERR;
	    goto out;
	}
	if (fseeko(fp, (off_t) (node * 4L), SEEK_CUR)) {
	    rval = ELIO_FSEEK_ERR;
	    goto out;
	}			/* sind */
	if (fread(&sind, sizeof(int), 1, fp) != 1) {
	    rval = ELIO_FREAD_ERR;
	    goto out;
	}
	npts = (h->nvrt - h->bi[node]) * h->ivs;	/* number of points to read for this node */
	if (fseeko(fp, skip + (off_t) ((2 + h->np + h->no[node]) * 4), SEEK_SET)) {	/* skip time, it, and the surf indexes */
	    rval = ELIO_FSEEK_ERR;
	    goto out;
	}
	if (fread(tmp, sizeof(float), npts, fp) != npts) {
	    rval = ELIO_FREAD_ERR;
	    goto out;
	}
    } else if (h->i23d == 2) {	/* 2d data - just one or 2 values to read */
	skip = (off_t) h->hsize + (off_t) step *(off_t) h->ssize;
	if (fseeko(fp, skip, SEEK_SET)) {	/* Read time, it, and sind */
	    rval = ELIO_FSEEK_ERR;
	    goto out;
	}
	if (fread(&t, sizeof(float), 1, fp) != 1) {
	    rval = ELIO_FREAD_ERR;
	    goto out;
	}
	if (fread(&it, sizeof(int), 1, fp) != 1) {
	    rval = ELIO_FREAD_ERR;
	    goto out;
	}
	if (fseeko(fp, (off_t) (node * 4L), SEEK_CUR)) {	/* Read time, it, and sind */
	    rval = ELIO_FSEEK_ERR;
	    goto out;
	}
	if (fread(&sind, sizeof(int), 1, fp) != 1) {
	    rval = ELIO_FREAD_ERR;
	    goto out;
	}
	npts = h->ivs;		/* number of points to read for this node */
	if (fseeko(fp, skip + (off_t) ((2 + h->np) * 4), SEEK_SET)) {	/* skip time, it, and the surf indexes */
	    rval = ELIO_FSEEK_ERR;
	    goto out;
	}
	if (fread(tmp, sizeof(float), npts, fp) != npts) {
	    rval = ELIO_FREAD_ERR;
	    goto out;
	}
    } else {
	rval = 9;
	goto out;
    }
    for (i = 0; i < (sind - bind) * h->ivs; i++) {
	d[i + bind * h->ivs] = tmp[i];
    }
  out:;
    free(tmp);
    return rval;
}

/*!
 * Get data for a single node from one time step should work for both 
 * sigma and z-level grids.
 *
 * Input variables:
 *
 * \param fp = Pointer to open file handle.
 * \param step = Time step starting from 0.
 * \param elem = Element number.
 * \param h = Pointer to valid ElcircHeader struct.
 * \param hh = Pointer to shape function coefficients.
 *
 * Output Variables:
 *
 * \param t = Pointer to time, float.
 * \param it = Pointer to computational time step, integer.
 * \param bind = Pointer to bottom index, integer.
 * \param sind = Pointer to surface index, integer.
 * \param d = Pointer to array of return values, h.nvrt long.
 *
 * Returns ELIO_OK on success.
 * Returns ELIO_FSEEK_ERR and ELIO_FREAD_ERR on error.
 */
int ElioGetXYData2(FILE * fp, int step, int elem, ElcircHeader * h, double *hh, float *t, int *it, int *bind, int *sind, float *d)
{
    int i, j, nn[4], tcnt, cnt[4];
    double uret, vret;
    double *u[4], *v[4], dd[4], vv[4];
    float dtmp[4][400];
    float dt[400];

/* get the nodes for the element */
    for (i = 0; i < h->nvrt * h->ivs; i++) {
	dt[i] = d[i] = 0.0;
    }
    for (i = 0; i < h->etype[elem]; i++) {
	nn[i] = h->icon[i][elem];
	ElioGetNode(fp, step, nn[i], h, t, it, &bind[i], &sind[i], dt);
	//printf("Node %d element %d %d %d\n", nn[i], elem, bind[i], sind[i]);
	for (j = 0; j < h->nvrt; j++) {
	    if (h->ivs == 2) {
		// printf("%d %f %f\n", i, dt[2 * j], dt[2 * j + 1]);
		if (bind[i] <= j && sind[i] >= j) {
		    dtmp[i][2 * j] = dt[2 * (j - bind[i])];
		    dtmp[i][2 * j + 1] = dt[2 * (j - bind[i]) + 1];
		} else {
		    dtmp[i][2 * j] = 0.0;
		    dtmp[i][2 * j + 1] = 0.0;
		}
	    } else if (h->ivs == 1) {
		//printf("ElioGetXYData2 = %d, %d: %f\n", i, j, dt[j]);
		if (bind[i] <= j && sind[i] >= j) {
		    dtmp[i][j] = dt[j - bind[i]];
		} else {
		    dtmp[i][j] = 0.0;
		}
	    }
	}
    }
    for (j = 0; j < h->nvrt; j++) {
	uret = -99.0;
	vret = 0.0;
	tcnt = 0;
	for (i = 0; i < h->etype[elem]; i++) {
	    cnt[i] = 0;
	    if (h->ivs == 2) {
		if (bind[i] <= j && sind[i] >= j) {
		    dd[i] = dtmp[i][2 * j];
		    vv[i] = dtmp[i][2 * j + 1];
		    cnt[i]++;
		    tcnt++;
		} else {
		    dd[i] = 0.0;
		    vv[i] = 0.0;
		}
	    } else if (h->ivs == 1) {
		if (bind[i] <= j && sind[i] >= j) {
		    dd[i] = dtmp[i][j];
		    cnt[i]++;
		    tcnt++;
		} else {
		    dd[i] = 0.0;
		}
	    }
	}
	if (h->ivs == 2) {
	    if (tcnt && tcnt != h->etype[elem]) {
		uret = 0.0;
		vret = 0.0;
		for (i = 0; i < h->etype[elem]; i++) {
		    uret = uret + cnt[i] * dd[i];
		    vret = vret + cnt[i] * vv[i];
		}
		uret = uret / tcnt;
		vret = vret / tcnt;
	    } else if (tcnt == h->etype[elem]) {
		ElioEval(h->etype[elem], hh, dd, &uret);
		ElioEval(h->etype[elem], hh, vv, &vret);
	    }
	    d[2 * j] = uret;
	    d[2 * j + 1] = vret;
	} else if (h->ivs == 1) {
	    if (tcnt && tcnt != h->etype[elem]) {
		uret = 0.0;
		for (i = 0; i < h->etype[elem]; i++) {
		    uret = uret + cnt[i] * dd[i];
		}
		uret = uret / tcnt;
	    } else if (tcnt == h->etype[elem]) {
		ElioEval(h->etype[elem], hh, dd, &uret);
	    }
	    d[j] = uret;
	}
    }
    return ELIO_OK;
}

/*!
 * Get data for a single node from one time step for sigma grids
 *
 * Input variables:
 *
 * \param h = Pointer to valid ElcircHeader struct.
 * \param elem = Element number.
 * \param x = x of location to interpolate to.
 * \param y = y of location to interpolate to.
 * \param hh = Pointer to shape function coefficients.
 * \param t = Pointer to valid ElcircTimeStep struct.
 *
 * Output Variables:
 *
 * \param d = Pointer to array of return values, h.nvrt long.
 *
 * Returns ELIO_OK on success.
 * Returns ELIO_FSEEK_ERR and ELIO_FREAD_ERR on error.
 */
int ElioInterpTimeStepSigma(ElcircHeader * h, int elem, double x, double y, double *hh, ElcircTimeStep * t, float *d)
{
    int i, j, nn[4], tcnt, cnt[4];
    double uret, vret;
    double *u[4], *v[4], dd[4], vv[4];
    float dtmp[4][400];
    float dt[400];
/* get the nodes for the element */
    for (i = 0; i < h->nvrt * h->ivs; i++) {
	dt[i] = d[i] = 0.0;
    }
    for (i = 0; i < h->etype[elem]; i++) {
	nn[i] = h->icon[i][elem];
	ElioExtractNode(t, h, nn[i], 0, h->nvrt - 1, dt);
	for (j = 0; j < h->nvrt; j++) {
//printf("%d, %d: %f\n", i, j, dt[j]);
	    if (h->ivs == 2) {
		dtmp[i][2 * j] = dt[2 * j];
		dtmp[i][2 * j + 1] = dt[2 * j + 1];
	    } else if (h->ivs == 1) {
		dtmp[i][j] = dt[j];
	    }
	}
    }
    for (j = 0; j < h->nvrt; j++) {
	for (i = 0; i < h->etype[elem]; i++) {
	    if (h->ivs == 2) {
		dd[i] = dtmp[i][2 * j];
		vv[i] = dtmp[i][2 * j + 1];
	    } else if (h->ivs == 1) {
		dd[i] = dtmp[i][j];
	    }
	}
	if (h->ivs == 2) {
	    ElioEval(h->etype[elem], hh, dd, &uret);
	    ElioEval(h->etype[elem], hh, vv, &vret);
	    d[2 * j] = uret;
	    d[2 * j + 1] = vret;
	} else {
	    ElioEval(h->etype[elem], hh, dd, &uret);
	    d[j] = uret;
	}
//printf("uret = %d, %d: %f\n", i, j, d[j]);
    }
    return ELIO_OK;
}

int ElioInterpTimeStep(ElcircHeader * h, int elem, double x, double y, double *hh, ElcircTimeStep * t, int *bind, int *sind, float *d)
{
    int i, j, nn[4], tcnt, cnt[4];
    double uret, vret;
    double *u[4], *v[4], dd[4], vv[4];
    float dtmp[4][400];
    float dt[400];

/* get the nodes for the element */
    for (i = 0; i < h->nvrt * h->ivs; i++) {
	dt[i] = d[i] = 0.0;
    }
    for (i = 0; i < h->etype[elem]; i++) {
	nn[i] = h->icon[i][elem];
	bind[i] = h->bi[nn[i]];
	sind[i] = t->surfind[nn[i]] - 1;
	ElioExtractNode(t, h, nn[i], bind[i], sind[i], dt);
	for (j = 0; j < h->nvrt; j++) {
	    if (h->ivs == 2) {
		if (bind[i] <= j && sind[i] >= j) {
		    dtmp[i][2 * j] = dt[2 * j];
		    dtmp[i][2 * j + 1] = dt[2 * j + 1];
		    //printf("values = %d, %d: %f %f\n", i, j, dt[2 * j], dt[2 * j + 1]);
		} else {
		    dtmp[i][2 * j] = 0.0;
		    dtmp[i][2 * j + 1] = 0.0;
		}
	    } else if (h->ivs == 1) {
		//if (h->nvrt == 1) printf("InterpTimeStep = %d, %d: %f %d %d\n", i, j, dt[j], bind[i], sind[i]);
		if (bind[i] <= j && sind[i] >= j) {
//if (tdebug && dt[j] < 0.0) {
//    printf("InterpTimeStep = %d, %d: %f %d %d\n", i, j, dt[j], bind[i], sind[i]);
//}
		    dtmp[i][j] = dt[j];
		} else {
		    dtmp[i][j] = 0.0;
		}
	    }
	}
    }
    for (j = 0; j < h->nvrt; j++) {
	uret = -99.0;
        vret = 0.0;
	tcnt = 0;
	for (i = 0; i < h->etype[elem]; i++) {
	    cnt[i] = 0;
	    if (h->ivs == 2) {
		if (bind[i] <= j && sind[i] >= j) {
		    dd[i] = dtmp[i][2 * j];
		    vv[i] = dtmp[i][2 * j + 1];
		    cnt[i]++;
		    tcnt++;
		} else {
		    dd[i] = 0.0;
		    vv[i] = 0.0;
		}
	    } else if (h->ivs == 1) {
		if (bind[i] <= j && sind[i] >= j) {
		    dd[i] = dtmp[i][j];
		    cnt[i]++;
		    tcnt++;
		} else {
		    dd[i] = 0.0;
		}
	    }
	}
	if (h->ivs == 2) {
	    if (tcnt && tcnt != h->etype[elem]) {
		uret = 0.0;
		vret = 0.0;
		for (i = 0; i < h->etype[elem]; i++) {
		    uret = uret + cnt[i] * dd[i];
		    vret = vret + cnt[i] * vv[i];
		    //printf("ElioInterpTimeStep: %d %d %d %lf %lf\n", tcnt, i, cnt[i], dd[i], vv[i]);
		}
		uret = uret / tcnt;
		vret = vret / tcnt;
	    } else if (tcnt == h->etype[elem]) {
		ElioEval(h->etype[elem], hh, dd, &uret);
		ElioEval(h->etype[elem], hh, vv, &vret);
		//printf("ElioInterpTimeStep OK: %lf %lf\n", uret, vret);
	    }
	    d[2 * j] = uret;
	    d[2 * j + 1] = vret;
	} else if (h->ivs == 1) {
	    if (tcnt && tcnt != h->etype[elem]) {
		uret = 0.0;
		for (i = 0; i < h->etype[elem]; i++) {
		    uret = uret + cnt[i] * dd[i];
		}
		uret = uret / tcnt;
	    } else if (tcnt == h->etype[elem]) {
		ElioEval(h->etype[elem], hh, dd, &uret);
	    }
	    d[j] = uret;
	}
    }
    return ELIO_OK;
}

/*!
 * int ElioGetXYData(FILE * fp, int step, double x, double y, ElcircHeader * h, float *t, int *it, int *bind, int *sind, float *d)
 *
 * Get data at (x, y) by interpolation.
 *
 * \param fp - Open file handle for the file to read. 
 * \param step  - Step to read starting from 0 - the first step. 
 * \param x - X.
 * \param y - Y.
 * \param h - Elcirc header information previously read. 
 * \param t - Return time from the file.
 * \param it - Return time step from file.
 * \param bind - Integer array of 4 values for returning bottom indexes for the element (x, y) is in.
 * \param sind - Integer array of 4 values for returning surface indexes for the element (x, y) is in.
 * \param d - Return data values.
 *
 * Returns ELIO_OK on success.
 * Returns ELIO_ERR on unknown failures.
 * Returns ELIO_FSEEK_ERR on seek failures.
 * Returns ELIO_FREAD_ERR on read failures.
 */
int ElioGetXYData(FILE * fp, int step, double x, double y, ElcircHeader * h, float *t, int *it, int *bind, int *sind, float *d)
{
    int i, j, elem, nn[4];
    double hh[4], uret;
    double *u[4], *v[4], dd[4];
    float dtmp[4][400];

/* find the element that contains the point */
    if ((elem = ElioFindElementXY(h, x, y)) < 0) {
	printf("elio.GetXY(): Unable to locate element.\n");
	return 1;
    }
/* get the nodes for the element */
    for (i = 0; i < h->etype[elem]; i++) {
	nn[i] = h->icon[i][elem];
	ElioGetNode(fp, step, nn[i], h, t, it, &bind[i], &sind[i], &dtmp[i][0]);
    }
    ElioGetCoefficients(h, elem, x, y, hh);
    for (i = 0; i < h->nvrt * h->ivs; i++) {
	d[i] = 0.0;
    }
    for (j = 0; j < h->nvrt; j++) {
	for (i = 0; i < h->etype[elem]; i++) {
	    if (bind[i] >= j && sind[i] < h->nvrt) {
		dd[i] = dtmp[i][j - bind[i]];
	    } else {
		dd[i] = -99.0;
	    }
	}
	ElioEval(h->etype[elem], hh, dd, &uret);
	d[j] = uret;
    }
    return ELIO_OK;
}

/*!
 *  int ElioGetPoint(FILE * fp, int step, int node, int level, ElcircHeader * h, float *t, int *it, float *d)
 * Get data for a single node from one time step at a particular level (numbered from 0)
 *
 * \param fp - Open file handle for the file to read. 
 * \param step  - Step to read starting from 0 - the first step. 
 * \param node  - Node to read starting from 0.
 * \param level - Level to read starting from 0. 
 * \param h - Elcirc header information previously read. 
 * \param t         - Return time value from data file. 
 * \param it,         - Return time step from data file.
 * \param d         - Reurn data value.
 *
 * Returns ELIO_OK on success.
 * Returns ELIO_ERR on unknown failures.
 * Returns ELIO_FSEEK_ERR on seek failures.
 * Returns ELIO_FREAD_ERR on read failures.
 */
int ElioGetPoint(FILE * fp, int step, int node, int level, ElcircHeader * h, float *t, int *it, float *d)
{
    int i, npts, sind;
    off_t skip;
    if (h->i23d == 3) {		/* 3d data */
	skip = (off_t) h->hsize + (off_t) step *(off_t) h->ssize;
	if (fseeko(fp, skip, SEEK_SET)) {	/* Read time, it, and sind */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(t, sizeof(float), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	if (fread(it, sizeof(int), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	if (fseeko(fp, node * 4L, SEEK_CUR)) {	/* sind */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(&sind, sizeof(int), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	npts = h->ivs;		/* number of points to read for this node */
	if (fseeko(fp, skip + (off_t) ((2 + h->np + h->no[node] + level) * 4), SEEK_SET)) {	/* skip time, it, and the surf indexes and go to the particular level of interest */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(d, sizeof(float), npts, fp) != npts) {
	    return ELIO_FREAD_ERR;
	}
	//printf("Starting level %d ending level %d\n", h->bi[node], sind);
	//for (i = 0; i < npts; i++) {
	//    printf("Data %d = %f\n", i + h->bi[node], d[i]);
	//}
	//printf("Completed Read of Time Step (%d): Time, Step = %f, %d\n", step, *t, *it);
    } else if (h->i23d == 2) {	/* 2d data - just one or 2 values to read */
	skip = (off_t) h->hsize + (off_t) step *(off_t) h->ssize;
	if (fseeko(fp, skip, SEEK_SET)) {	/* Read time, it, and sind */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(&t, sizeof(float), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	if (fread(&it, sizeof(int), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	if (fseeko(fp, node * 4L, SEEK_CUR)) {	/* Read time, it, and sind */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(&sind, sizeof(int), 1, fp) != 1) {
	    return ELIO_FREAD_ERR;
	}
	npts = h->ivs;		/* number of points to read for this node */
	if (fseeko(fp, skip + (off_t) ((2 + h->np) * 4), SEEK_SET)) {	/* skip time, it, and the surf indexes */
	    return ELIO_FSEEK_ERR;
	}
	if (fread(d, sizeof(float), npts, fp) != npts) {
	    return ELIO_FREAD_ERR;
	}
/*
	printf("Starting level %d ending level %d\n", h->bi[node], sind);
	for (i = 0; i < npts; i++) {
	    printf("Data %d = %f\n", i, d[i]);
	}
*/
	printf("Completed Read of Time Step (%d): Time, Step = %f, %d\n", step, *t, *it);
    } else {
	return ELIO_ERR;
    }
    return ELIO_OK;
}

/*!
 * Read the ElcircHeader from a data file
 *
 * Input variables:
 *
 * \param fname = Pointer to a character array with the file name to read.
 *
 * Output variables:
 *
 * \param h = Pointer to ElcircHeader header. Memory will be allocated
 *     to accomodate the internal variables.
 *
 * Returns ELIO_OK on success.
 * Returns non-zero on failure.
 *
 * Notes: There are some transformations of header variables performed for 
 * convenience. For 2d variables such as elevations and winds, nvrt (the 
 * number of levels) is set to one, the bottom index for each node is set 
 * to 0 for all nodes. In the case of 3d variables the bottom indices 
 * ('bi') are decremented by one.
 *
 */
int ElioGetHeader(char *fname, ElcircHeader * h)
{
    FILE *fp;
    int *itmp, i, j, ss;
    if (sizeof(int) != 4 || sizeof(char) != 1 || sizeof(float) != 4) {
	fprintf(stderr, "Sizeof problems, please investigate: sizeof(char,int,float) = (%d,%d,%d)\n", sizeof(int), sizeof(char), sizeof(float));
	exit(1);
    }
    fp = fopen(fname, "rb");
    if (fp == NULL) {
	fprintf(stderr, "GetElioHeader(): Unable to open file %s\n", fname);
	return 1;
    }
    if (fread(h->magic, sizeof(char), 48, fp) != 48)
	return ELIO_FREAD_ERR;
    h->magic[48] = 0;
    /* DataFormat v3.0 */
    if (strncmp(h->magic, "DataFormat v5.0", 15) == 0) {
	h->v = 5;
	h->sigma = 2;		/* Sigma + Z-level grid */
    } else if (strncmp(h->magic, "DataFormat v4.0", 15) == 0) {
	h->v = 4;
	h->sigma = 1;		/* sigma-level grid */
    } else if (strncmp(h->magic, "DataFormat v3.0", 15) == 0) {
	h->v = 3;
	h->sigma = 0;		/* Z-level grid */
    } else if (strncmp(h->magic, "DataFormat v2.0", 15) == 0) {
	h->v = 2;
	h->sigma = 0;		/* Z-level grid */
    } else {
	fprintf(stderr, "GetElioHeader(): Unknown version identifier %s\n", h->magic);
	return (3);
    }
    if (fread(h->version, sizeof(char), 48, fp) != 48)
	return ELIO_FREAD_ERR;
    h->version[48] = 0;
    if (fread(h->start_time, sizeof(char), 48, fp) != 48)
	return ELIO_FREAD_ERR;
    h->start_time[48] = 0;
    if (fread(h->variable_nm, sizeof(char), 48, fp) != 48)
	return ELIO_FREAD_ERR;
    h->variable_nm[48] = 0;
    if (fread(h->variable_dim, sizeof(char), 48, fp) != 48)
	return ELIO_FREAD_ERR;
    h->variable_dim[48] = 0;
    if (fread(&h->nsteps, sizeof(int), 1, fp) != 1)
	return ELIO_FREAD_ERR;
    if (fread(&h->timestep, sizeof(float), 1, fp) != 1)
	return ELIO_FREAD_ERR;
    if (fread(&h->skip, sizeof(int), 1, fp) != 1)
	return ELIO_FREAD_ERR;
    if (fread(&h->ivs, sizeof(int), 1, fp) != 1)
	return ELIO_FREAD_ERR;
    if (fread(&h->i23d, sizeof(int), 1, fp) != 1)
	return ELIO_FREAD_ERR;

    if (h->v == 2 || h->v == 3) {
	if (fread(&h->vpos, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->zmsl, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->nvrt, sizeof(int), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
    } else if (h->v == 4) {
	if (fread(&h->ivcor, sizeof(int), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->h0, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->hc, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->thetab, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->thetaf, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->nvrt, sizeof(int), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
    } else if (h->v == 5) {
	if (fread(&h->nvrt, sizeof(int), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->kz, sizeof(int), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	h->ks = h->nvrt - h->kz;
	if (fread(&h->h0, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->hs, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->hc, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->thetab, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(&h->thetaf, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
    }
    h->zcor = (float *) malloc(h->nvrt * sizeof(float));
    if (fread(h->zcor, sizeof(float), h->nvrt, fp) != h->nvrt)
	return ELIO_FREAD_ERR;

    if (fread(&h->np, sizeof(int), 1, fp) != 1)
	return ELIO_FREAD_ERR;
    if (fread(&h->ne, sizeof(int), 1, fp) != 1)
	return ELIO_FREAD_ERR;
    h->x = (float *) malloc(h->np * sizeof(float));
    h->y = (float *) malloc(h->np * sizeof(float));
    h->d = (float *) malloc(h->np * sizeof(float));
    h->bi = (int *) malloc(h->np * sizeof(int));
    h->bisave = (int *) malloc(h->np * sizeof(int));
    h->no = (int *) malloc(h->np * sizeof(int));	/* offset into data for each node */
    for (i = 0; i < h->np; i++) {
	if (fread(h->x + i, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(h->y + i, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
	if (fread(h->d + i, sizeof(float), 1, fp) != 1)
	    return ELIO_FREAD_ERR;

	if (h->v != 4) {
	    if (fread(h->bi + i, sizeof(int), 1, fp) != 1)
		return ELIO_FREAD_ERR;
	} else {
	    h->bi[i] = 1;
	}
    }
    h->icon[0] = (int *) malloc(h->ne * sizeof(int));
    h->icon[1] = (int *) malloc(h->ne * sizeof(int));
    h->icon[2] = (int *) malloc(h->ne * sizeof(int));
    h->etype = (int *) malloc(h->ne * sizeof(int));
    if (h->v == 2) {		/* format version 2 */
	itmp = (int *) malloc(h->ne * 3 * sizeof(int));
	if (fread(itmp, sizeof(int), h->ne * 3, fp) != h->ne * 3)
	    return ELIO_FREAD_ERR;
	for (i = 0; i < h->ne; i++) {
	    h->etype[i] = 3;
	    h->icon[0][i] = itmp[i * 3] - 1;
	    h->icon[1][i] = itmp[i * 3 + 1] - 1;
	    h->icon[2][i] = itmp[i * 3 + 2] - 1;
	}
	free(itmp);
    } else if (h->v == 3) {	/* format version 3 support for quads */
	h->icon[3] = (int *) malloc(h->ne * sizeof(int));
	for (i = 0; i < h->ne; i++) {
	    if (fread(&h->etype[i], sizeof(int), 1, fp) != 1)
		return ELIO_FREAD_ERR;
	    for (j = 0; j < h->etype[i]; j++) {
		if (fread(&h->icon[j][i], sizeof(int), 1, fp) != 1)
		    return ELIO_FREAD_ERR;
		h->icon[j][i]--;
	    }
	}
    } else if (h->v == 4 || h->v == 5) {	/* format version 4 no support for quads */
	for (i = 0; i < h->ne; i++) {
	    if (fread(&h->etype[i], sizeof(int), 1, fp) != 1)
		return ELIO_FREAD_ERR;
	    for (j = 0; j < h->etype[i]; j++) {
		if (fread(&h->icon[j][i], sizeof(int), 1, fp) != 1)
		    return ELIO_FREAD_ERR;
		h->icon[j][i]--;
	    }
	}
    }
/* compute size of data in time step and offsets to each node */
    if (h->i23d == 3) {
	if (h->v != 4) {
	    h->nitems = 0;
	    for (i = 0; i < h->np; i++) {
		h->bisave[i] = h->bi[i];
		if (h->bi[i] > 0) {
		    h->bi[i]--;
		} else {
		}
		h->no[i] = h->nitems;
		h->nitems += (h->nvrt - h->bi[i]) * h->ivs;
//printf("%d: %d %d %d %d %d\n", i, h->nvrt, h->ivs, h->bi[i], h->no[i], h->nitems);
	    }
	} else {
	    h->nitems = h->ivs * h->nvrt * h->np;
	    for (i = 0; i < h->np; i++) {
		h->bisave[i] = h->bi[i];
		h->bi[i] = 0;
		h->no[i] = i * h->ivs * h->nvrt;
	    }
	}
    } else {
	h->nvrt = 1;		/* only 1 level */
	if (h->v != 4) {
	    h->nitems = 0;
	    for (i = 0; i < h->np; i++) {
		h->bisave[i] = h->bi[i];
		h->bi[i] = 0;	/* bottom index is always 0 */
		h->no[i] = h->nitems;
		h->nitems += h->ivs;
	    }
	} else {
	    h->nitems = h->ivs * h->np;
	    for (i = 0; i < h->np; i++) {
		h->bisave[i] = h->bi[i];
		h->bi[i] = 0;	/* bottom index is always 0 */
		h->no[i] = i * h->ivs;
	    }
	}
    }
/* get the header size in bytes */
    h->hsize = ftello(fp);
/* get time step size in bytes */
    if (h->v == 3 || h->v == 4 || h->v == 5) {	/* version 4 has elevation so size is the same */
	h->ssize = 8 + h->np * 4 + h->nitems * 4;
    } else if (h->v == 2) {
	h->ssize = 8 + h->np * 4 + h->nitems * 4;
    }
    if (fclose(fp)) {
	return ELIO_FCLOSE_ERR;
    }
    return ELIO_OK;
}

/*!
 * Allocate memory for an ElcircHeader
 *
 * Input/Output variable:
 *
 * \param h = Pointer to ElcircHeader header. Memory will be allocated
 *     to accomodate the internal variables. Assumes suitable defaults
 *     have been set for np, ne, and nvrt. Also that other variables 
 *     have proper defaults set.
 *
 * Returns ELIO_OK on success.
 * Returns non-zero on failure.
 *
 */
int ElioAllocateHeader(ElcircHeader * h)
{
    h->zcor = (float *) malloc(h->nvrt * sizeof(float));
    h->x = (float *) malloc(h->np * sizeof(float));
    h->y = (float *) malloc(h->np * sizeof(float));
    h->d = (float *) malloc(h->np * sizeof(float));
    h->bi = (int *) malloc(h->np * sizeof(int));
    h->bisave = (int *) malloc(h->np * sizeof(int));	/* save the original bi values */
    h->no = (int *) malloc(h->np * sizeof(int));	/* offset into data for each node */
    h->icon[0] = (int *) malloc(h->ne * sizeof(int));
    h->icon[1] = (int *) malloc(h->ne * sizeof(int));
    h->icon[2] = (int *) malloc(h->ne * sizeof(int));
    h->icon[3] = (int *) malloc(h->ne * sizeof(int));
    h->etype = (int *) malloc(h->ne * sizeof(int));
    return ELIO_OK;
}

int ElioMakeScalarsOld(ElcircHeader * h, float *d, float *dd)
{
    int i, j, k;
    for (i = 0; i < h->np * h->nvrt; i++) {
	dd[i] = -99.0;
    }
    for (j = 0; j < h->np; j++) {
	for (k = h->bi[j]; k < h->nvrt; k++) {
	    dd[j * h->nvrt + k] = d[h->no[j] + k - h->bi[j]];
	}
    }
    return ELIO_OK;
}

int ElioMakeVectorsOld(ElcircHeader * h, float *d, float *dd)
{
    int i, j, k;
    float tmp;
    for (i = 0; i < 2 * h->np * h->nvrt; i++) {
	dd[i] = 0.0;
    }
    for (j = 0; j < h->np; j++) {
	for (k = h->bi[j]; k < h->nvrt; k++) {
	    tmp = d[h->no[j] + (k - h->bi[j]) * 2];
	    dd[2 * (j * h->nvrt + k)] = tmp;
	    tmp = d[h->no[j] + (k - h->bi[j]) * 2 + 1];
	    dd[2 * (j * h->nvrt + k) + 1] = tmp;
/*
	    dd[2 * (j * h->nvrt + k)] = d[h->no[j] + (k - h->bi[j]) * 2];
	    dd[2 * (j * h->nvrt + k) + 1] = d[h->no[j] + (k - h->bi[j]) * 2 + 1];
*/
	}
    }
    return ELIO_OK;
}

/*! 
 * Return the actual number of time steps in a data file regardless of 
 * what the ElcircHeader says.
 */
int ElioGetNStepsInFile(char *fname, ElcircHeader * h)
{
    char id[5];
    int n;
    off_t fsize = 0;
    FILE *fp = fopen(fname, "rb");
    if (fp == NULL) {
	fprintf(stderr, "GetNStepsInFile(): Unable to open file %s\n", fname);
	return ELIO_OK;
    }
    if (fseeko(fp, 0L, SEEK_END)) {
	return ELIO_FSEEK_ERR;
    }
    fsize = ftello(fp);
    if (fclose(fp)) {
	return ELIO_FCLOSE_ERR;
    }
    fsize = fsize - (off_t) h->hsize;
    if (h->ssize > 0) {
	n = fsize / h->ssize;
    } else {
	n = -1;
    }
    return n;
}

/*! 
 * Extract a Header for the extracted grid
 * x, y, isin need to be allocated in advance.
 * invert = 1 to select outside the region
 * most == 1 to include an element if a single node is in the region
 * NOTE: most does not work, if a node is in then every node of every element that uses that node is also in
 * Need to add some additional logic for most to work.
 */
int ElioExtractGrid(ElcircHeader * h1, int nb, double *xb, double *yb, int *isin, int invert, int most, ElcircHeader * h2)
{
    int i, j;
    int n;
    int np, ne;
    int n1, n2, n3, n4;
    double x, y;

    *h2 = *h1;			/* copy the easy stuff */

/*
 * Find the 'in' nodes and assign the new node numbers for the reduced grid.
 */
    for (i = 0; i < h1->np; i++) {
	isin[i] = 0;
    }
    n = 1;			/* start from 1 for truth testing (0 is not in) */
    for (i = 0; i < h1->np; i++) {
	x = (double) h1->x[i];
	y = (double) h1->y[i];
	if (!isin[i]) {		// if most == 1 the test for already in
	    if (invert) {
		if (!ElioInPolygon(x, y, nb, xb, yb)) {
		    isin[i] = n;
		    n++;
		} else {
		    isin[i] = 0;
		}
	    } else {
		if (ElioInPolygon(x, y, nb, xb, yb)) {
		    isin[i] = n;
		    n++;
		} else {
		    isin[i] = 0;
		}
	    }
	}
	if (isin[i] && most) {
// TODO need to add logic to mark all associated nodes as in also.
//          crude version - loop over all elements and find elements using this node then
//          mark all other nodes of the element as in.
//          Better version uses an edgelist.
	    for (j = 0; j < h1->ne; j++) {
		n1 = h1->icon[0][j];
		n2 = h1->icon[1][j];
		n3 = h1->icon[2][j];
		if (n1 == i) {
		    if (!isin[n2]) {
			isin[n2] = n;
			n++;
		    }
		    if (!isin[n3]) {
			isin[n3] = n;
			n++;
		    }
		} else if (n2 == i) {
		    if (!isin[n1]) {
			isin[n1] = n;
			n++;
		    }
		    if (!isin[n3]) {
			isin[n3] = n;
			n++;
		    }
		} else if (n3 == i) {
		    if (!isin[n1]) {
			isin[n1] = n;
			n++;
		    }
		    if (!isin[n2]) {
			isin[n2] = n;
			n++;
		    }
		}
	    }
	}
    }

    h2->np = n - 1;		/* one extra */
    h2->x = (float *) malloc(h2->np * sizeof(float));
    h2->y = (float *) malloc(h2->np * sizeof(float));
    h2->d = (float *) malloc(h2->np * sizeof(float));
    h2->bi = (int *) malloc(h2->np * sizeof(int));
    h2->bisave = (int *) malloc(h2->np * sizeof(int));
    h2->no = (int *) malloc(h2->np * sizeof(int));
/*
 * count the number of elements in the new grid
 */
    n = 0;
    for (i = 0; i < h1->ne; i++) {
	n1 = h1->icon[0][i];
	n2 = h1->icon[1][i];
	n3 = h1->icon[2][i];
	if (h1->etype[i] == 3) {
	    if (isin[n1] && isin[n2] && isin[n3]) {
		++n;
	    }
	} else if (h1->etype[i] == 4) {
	    n4 = h1->icon[3][i];
	    if (isin[n1] && isin[n2] && isin[n3] && isin[n4]) {
		++n;
	    }
	}
    }
    h2->ne = n;
    h2->icon[0] = (int *) malloc(h2->ne * sizeof(int));
    h2->icon[1] = (int *) malloc(h2->ne * sizeof(int));
    h2->icon[2] = (int *) malloc(h2->ne * sizeof(int));
    h2->icon[3] = (int *) malloc(h2->ne * sizeof(int));
    h2->etype = (int *) malloc(h2->ne * sizeof(int));

/*
 * Build new header
 */
    n = 0;
    for (i = 0; i < h1->np; i++) {
	if (isin[i]) {
	    h2->x[n] = h1->x[i];
	    h2->y[n] = h1->y[i];
	    h2->d[n] = h1->d[i];
	    h2->bi[n] = h1->bi[i];
	    h2->bisave[n] = h1->bisave[i];
	    n++;
	}
    }
/*
 * use the old node numbers to hash into the list containing the new
 * node numbers when writing out the table of elements.
 */
    n = 0;
    for (i = 0; i < h1->ne; i++) {
	n1 = h1->icon[0][i];
	n2 = h1->icon[1][i];
	n3 = h1->icon[2][i];
	h2->etype[n] = h1->etype[i];
	if (h1->etype[i] == 4) {
	    n4 = h1->icon[3][i];
	    if (isin[n1] && isin[n2] && isin[n3] && isin[n4]) {
		h2->icon[0][n] = isin[n1] - 1;	/* isin starts from 1 */
		h2->icon[1][n] = isin[n2] - 1;
		h2->icon[2][n] = isin[n3] - 1;
		h2->icon[3][n] = isin[n4] - 1;
		n++;
	    }
	} else {
	    if (isin[n1] && isin[n2] && isin[n3]) {
		h2->icon[0][n] = isin[n1] - 1;	/* isin starts from 1 */
		h2->icon[1][n] = isin[n2] - 1;
		h2->icon[2][n] = isin[n3] - 1;
		n++;
	    }
	}
    }
/* get the computed elements of the header */
    if (h2->i23d == 3) {
	h2->nitems = 0;
	for (i = 0; i < h2->np; i++) {
	    h2->no[i] = h2->nitems;
	    h2->nitems += (h2->nvrt - h2->bi[i]) * h2->ivs;
	}
    } else {
	h2->nvrt = 1;		/* only 1 level */
	h2->nitems = 0;
	for (i = 0; i < h2->np; i++) {
	    h2->bi[i] = 0;	/* bottom index is always 0 */
	    h2->no[i] = h2->nitems;
	    h2->nitems += h2->ivs;
	}
    }
/* get the header size in bytes TODO */
    /*h2->hsize = ftello(fp); */
/* get time step size in bytes */
    h2->ssize = 8 + h2->np * 4 + h2->nitems * 4;
    return ELIO_OK;
}

/*!
 * Extract data from a sampled grid
 * The header must be created before this function is called 
 */
int ElioExtractData(ElcircHeader * h1, ElcircHeader * h2, int *isin, ElcircTimeStep t1, ElcircTimeStep * t2)
{
    int i, j, n, *si, npts;
    float *d;
    *t2 = t1;
    si = (int *) malloc(h2->np * sizeof(int));
    n = 0;
    for (i = 0; i < h1->np; i++) {
	if (isin[i]) {
	    si[n] = t1.surfind[i];
	    n++;
	}
    }
    t2->surfind = si;
    t2->d = (float *) malloc(h2->nitems * sizeof(float));
    n = 0;
    for (i = 0; i < h1->np; i++) {
	if (isin[i]) {
	    for (j = 0; j < h1->nvrt - h1->bi[i]; j++) {
		t2->d[n] = t1.d[h1->no[i] + j];
		n++;
	    }
	}
    }
    return ELIO_OK;
}

/*!
 * Routines to determine if a point lies in a polygon
*/
int ElioIntersectToLeft(double x, double y, double x1, double y1, double x2, double y2)
{
    double xtmp, m, b;

    /* ignore horizontal lines */
    if (y1 == y2) {
	return 0;
    }
    /* not contained vertically */
    if (((y < y1) && (y < y2)) || ((y > y1) && (y > y2))) {
	return 0;
    }
    /* none of the above, compute the intersection */
    if ((xtmp = x2 - x1) != 0.0) {
	m = (y2 - y1) / xtmp;
	b = y1 - m * x1;
	xtmp = (y - b) / m;
    } else {
	xtmp = x1;
    }
    if (xtmp <= x) {
	/* check for intersections at a vertex */
	/* if this is the max ordinate then accept */
	if (y == y1) {
	    if (y1 > y2) {
		return 1;
	    } else {
		return 0;
	    }
	}
	/* check for intersections at a vertex */
	if (y == y2) {
	    if (y2 > y1) {
		return 1;
	    } else {
		return 0;
	    }
	}
	/* no vertices intersected */
	return 1;
    }
    return 0;
}

/*!
 * Determine if (x,y) is in the polygon xlist[], ylist[].
 */
int ElioInPolygon(double x, double y, int n, double *xlist, double *ylist)
{
    int i, l = 0, ll = 0;

    for (i = 0; i < n; i++) {
	if ((y < ylist[i] && y < ylist[(i + 1) % n])
	    || (y > ylist[i] && y > ylist[(i + 1) % n])) {
	    continue;
	}
	if ((x < xlist[i] && x < xlist[(i + 1) % n])) {
	    continue;
	}
	l += ElioIntersectToLeft(x, y, xlist[i], ylist[i], xlist[(i + 1) % n], ylist[(i + 1) % n]);
    }
    return l % 2;
}

/*!
 * Locate the element in which x, y is located.
 * Returns -1 if the element cannot be located, the element number
 * starting from zero otherwise.
 *
 */
int ElioFindElementXY(ElcircHeader * h, double xp, double yp)
{
    int i, j;
    int n0, n1, n2, n3;
    double x[4], y[4];
    for (i = 0; i < h->ne; i++) {
	for (j = 0; j < h->etype[i]; j++) {
	    x[j] = h->x[h->icon[j][i]];
	    y[j] = h->y[h->icon[j][i]];
	}
	if (h->etype[i] == 3) {
	    if (ElioInsideElement(xp, yp, x[0], y[0], x[1], y[1], x[2], y[2])) {
		return i;
	    }
	} else {
	    if (ElioInsideElement4(xp, yp, x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3])) {
		return i;
	    }
	}
    }
    return -1;
}

/*!
 * Tolerance for locating point in elements.
 */
static double ElioInsideElementTol = -1e-05;

/*!
 * Set the tolerance for locating points in elements.
 */
void ElioSetInsideElementTol(double tol)
{
    ElioInsideElementTol = -tol;
}

/*!
 * Determine if (x, y) is inside the triangle given by (x1, y1), (x2, y2), (x3, y3).
 */
int ElioInsideElement(double x, double y, double x1, double y1, double x2, double y2, double x3, double y3)
{
    if ((x - x1) * (y1 - y2) + (y - y1) * (x2 - x1) < ElioInsideElementTol)
	return 0;
    if ((x - x2) * (y2 - y3) + (y - y2) * (x3 - x2) < ElioInsideElementTol)
	return 0;
    if ((x - x3) * (y3 - y1) + (y - y3) * (x1 - x3) < ElioInsideElementTol)
	return 0;
    return 1;
}

/*!
 * Determine if (x, y) is inside the convex quadrangle given by (x1, y1), (x2, y2), (x3, y3), (x4, y4).
 */
int ElioInsideElement4(double x, double y, double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4)
{
    if ((x - x1) * (y1 - y2) + (y - y1) * (x2 - x1) < ElioInsideElementTol)
	return 0;
    if ((x - x2) * (y2 - y3) + (y - y2) * (x3 - x2) < ElioInsideElementTol)
	return 0;
    if ((x - x3) * (y3 - y4) + (y - y3) * (x4 - x3) < ElioInsideElementTol)
	return 0;
    if ((x - x4) * (y4 - y1) + (y - y4) * (x1 - x4) < ElioInsideElementTol)
	return 0;
    return 1;
}

/*!
 * Given previously computed nodal coefficients in h, and flow (u, v) at each node compute the flow and return in (uret, vret).
 *
 * \param n Integer number of nodes.
 * \param h Pointer to nodal coefficients.
 * \param u Pointer to nodal values for u at the nodes.
 * \param v Pointer to nodal values for u at the nodes.
 * \param uret Pointer to return value.
 * \param vret Pointer to return value.
 */
int ElioEvalFlowXY(int n, double *h, double *u, double *v, double *uret, double *vret)
{
    int i;
    *uret = 0.0;
    *vret = 0.0;
    if (n != 3 && n != 4) {
	return 1;
    } else {
	for (i = 0; i < n; i++) {
	    *uret += h[i] * u[i];
	    *vret += h[i] * v[i];
	}
    }
    return ELIO_OK;
}

int ElioEval(int n, double *h, double *c, double *uret)
{
    int i;
    *uret = 0.0;
    if (n != 3 && n != 4) {
	return 1;
    } else {
	for (i = 0; i < n; i++) {
	    *uret += h[i] * c[i];
	}
    }
    return ELIO_OK;
}

int ElioEvalScalarXY(int n, double *h, double *c, double *uret)
{
    int i;
    *uret = 0.0;
    if (n != 3 && n != 4) {
	return 1;
    } else {
	for (i = 0; i < n; i++) {
	    *uret += h[i] * c[i];
	}
    }
    return ELIO_OK;
}

/*! Compute coefficients for interpolation for both triangles and quads using
 * information from the element number.
 *
 * \param h Header struct to write.
 * \param elem  Number of the element from 0.
 * \param xp X-coordinate of point for interpolation.
 * \param yp Y-coordinate of point for interpolation.
 * \param w Coefficients of interpolator.
 *
 * Returns ELIO_OK on success and ELIO_ERR on failure.
 */
int ElioGetCoefficients(ElcircHeader * h, int elem, double xp, double yp, double *w)
{
    double aum, ado, atr, bum, bdo, btr, cum, cdo, ctr;
    double x1, y1, x2, y2, x3, y3, x4, y4, c1, c2, c3, arei;
    double xi, eta;
    int nn[4];
    nn[0] = h->icon[0][elem];
    nn[1] = h->icon[1][elem];
    nn[2] = h->icon[2][elem];
    x1 = h->x[nn[0]];
    y1 = h->y[nn[0]];
    x2 = h->x[nn[1]];
    y2 = h->y[nn[1]];
    x3 = h->x[nn[2]];
    y3 = h->y[nn[2]];
    if (h->etype[elem] == 3) {
	aum = x2 * y3 - x3 * y2;
	bum = y2 - y3;
	cum = x3 - x2;
	ado = x3 * y1 - x1 * y3;
	bdo = y3 - y1;
	cdo = x1 - x3;
	atr = x1 * y2 - x2 * y1;
	btr = y1 - y2;
	ctr = x2 - x1;
	arei = 1.0 / (aum + ado + atr);
	w[2] = (atr + btr * xp + ctr * yp) * arei;
	w[1] = (ado + bdo * xp + cdo * yp) * arei;
	w[0] = 1.0 - w[1] - w[2];
	return ELIO_OK;
    } else {
	nn[3] = h->icon[3][elem];
	x4 = h->x[nn[3]];
	y4 = h->y[nn[3]];
	ibilinear(x1, x2, x3, x4, y1, y2, y3, y4, xp, yp, &xi, &eta, w);
	return ELIO_OK;
    }
    return ELIO_ERR;
}

/*! Compute coefficients for interpolation for both triangles and quads using
 * information from the element number.
 * 
 *  \param h Header struct to write.
 *  \param elem  Number of the element from 0.
 *  \param xp X-coordinate of point for interpolation.
 *  \param yp Y-coordinate of point for interpolation.
 *  \param w Coefficients of interpolator.
 * 
 *  Returns ELIO_OK on success and ELIO_ERR on failure.
 */

int ElioGetCoefficientsGrid(ElioGrid * g, int elem, double xp, double yp, double *w)
{
    double aum, ado, atr, bum, bdo, btr, cum, cdo, ctr;
    double x1, y1, x2, y2, x3, y3, x4, y4, c1, c2, c3, arei;
    double xi, eta;
    int nn[4];
    if (g == NULL || elem < 0 || elem >= g->ne) {
	return ELIO_ERR;
    }
    nn[0] = g->icon[0][elem];
    nn[1] = g->icon[1][elem];
    nn[2] = g->icon[2][elem];
    x1 = g->x[nn[0]];
    y1 = g->y[nn[0]];
    x2 = g->x[nn[1]];
    y2 = g->y[nn[1]];
    x3 = g->x[nn[2]];
    y3 = g->y[nn[2]];
    if (g->etype[elem] == 3) {
	aum = x2 * y3 - x3 * y2;
	bum = y2 - y3;
	cum = x3 - x2;
	ado = x3 * y1 - x1 * y3;
	bdo = y3 - y1;
	cdo = x1 - x3;
	atr = x1 * y2 - x2 * y1;
	btr = y1 - y2;
	ctr = x2 - x1;
	arei = 1.0 / (aum + ado + atr);
	w[2] = (atr + btr * xp + ctr * yp) * arei;
	w[1] = (ado + bdo * xp + cdo * yp) * arei;
	w[0] = 1.0 - w[1] - w[2];
	return ELIO_OK;
    } else {
	nn[3] = g->icon[3][elem];
	x4 = g->x[nn[3]];
	y4 = g->y[nn[3]];
	ibilinear(x1, x2, x3, x4, y1, y2, y3, y4, xp, yp, &xi, &eta, w);
	return ELIO_OK;
    }
    return ELIO_ERR;
}

/*! Compute coefficients for interpolation for both triangles and quads using
 * coordinates of the nodes rather than information from the element number.
 *
 * \param n Number of nodes in the element.
 * \param x Pointer to array with x-coordinate of nodes.
 * \param y Pointer to array with y-coordinate of nodes.
 * \param xp X-coordinate of point for interpolation.
 * \param yp Y-coordinate of point for interpolation.
 * \param w Coefficients of interpolator.
 *
 * Returns ELIO_OK on success and ELIO_ERR on failure.
 */
int ElioGetCoefficientsXY(int n, double *x, double *y, double xp, double yp, double *w)
{
    double aum, ado, atr, bum, bdo, btr, cum, cdo, ctr;
    double x1, y1, x2, y2, x3, y3, x4, y4, c1, c2, c3, arei;
    double xi, eta;
    x1 = x[0];
    y1 = y[0];
    x2 = x[1];
    y2 = y[1];
    x3 = x[2];
    y3 = y[2];
    if (n == 3) {
	aum = x2 * y3 - x3 * y2;
	bum = y2 - y3;
	cum = x3 - x2;
	ado = x3 * y1 - x1 * y3;
	bdo = y3 - y1;
	cdo = x1 - x3;
	atr = x1 * y2 - x2 * y1;
	btr = y1 - y2;
	ctr = x2 - x1;
	arei = 1.0 / (aum + ado + atr);
	w[2] = (atr + btr * xp + ctr * yp) * arei;
	w[1] = (ado + bdo * xp + cdo * yp) * arei;
	w[0] = 1.0 - w[1] - w[2];
	return ELIO_OK;
    }
    if (n == 4) {
	x4 = x[3];
	y4 = y[3];
	ibilinear(x1, x2, x3, x4, y1, y2, y3, y4, xp, yp, &xi, &eta, w);
	return ELIO_OK;
    } else {
	return ELIO_ERR;
    }
    return ELIO_ERR;
}

/*!
 * Define a header for the surface.
*/
int ElioGetSurfaceHeader(ElcircHeader * h1, ElcircHeader * h2)
{
    int i, j;
    strncpy(h2->magic, h1->magic, 48);
    strncpy(h2->version, h1->version, 48);
    strncpy(h2->start_time, h1->start_time, 48);
    strncpy(h2->variable_nm, h1->variable_nm, 48);
    strncpy(h2->variable_dim, h1->variable_dim, 48);
    h2->nsteps = h1->nsteps;
    h2->timestep = h1->timestep;
    h2->skip = h1->skip;
    h2->ivs = h1->ivs;
    h2->i23d = h1->i23d;
    h2->vpos = h1->vpos;
    h2->zmsl = h1->zmsl;
    h2->nvrt = h1->nvrt;
    h2->zcor = (float *) malloc(h1->nvrt * sizeof(float));
    for (i = 0; i < h1->nvrt; i++) {
	h2->zcor = h1->zcor;
    }
    h2->np = h1->np;
    h2->ne = h1->ne;
    h2->x = (float *) malloc(h1->np * sizeof(float));
    h2->y = (float *) malloc(h1->np * sizeof(float));
    h2->d = (float *) malloc(h1->np * sizeof(float));
    h2->bi = (int *) malloc(h1->np * sizeof(int));
    h2->bisave = (int *) malloc(h1->np * sizeof(int));
    h2->no = (int *) malloc(h1->np * sizeof(int));	/* offset into data for each node */
    for (i = 0; i < h1->np; i++) {
	h2->x[i] = h1->x[i];
	h2->y[i] = h1->y[i];
	h2->d[i] = h1->d[i];
	h2->bisave[i] = h1->bisave[i];
	h2->bi[i] = h1->bi[i];
    }
    h2->icon[0] = (int *) malloc(h1->ne * sizeof(int));
    h2->icon[1] = (int *) malloc(h1->ne * sizeof(int));
    h2->icon[2] = (int *) malloc(h1->ne * sizeof(int));
    for (i = 0; i < h1->ne; i++) {
	h2->icon[0][i] = h1->icon[0][i];
	h2->icon[1][i] = h1->icon[0][i];
	h2->icon[2][i] = h1->icon[0][i];
    }
/* compute size of data in time step and offsets to each node */
    h2->nitems = 0;
    for (i = 0; i < h1->np; i++) {
	h2->bi[i] = 1;
	h2->no[i] = h2->nitems;
	h2->nitems += 1;
    }
/* get the header size in bytes */
    /*h2->hsize = ftello(fp); */
/* get time step size in bytes */
    h2->ssize = 8 + h2->np * 4 + h2->nitems * 4;
    return ELIO_OK;
}

/*!
 * Get a time step for the surface. Not implemented.
 */
int ElioGetSurfaceStep(ElcircHeader * h1, ElcircHeader * h2, ElcircTimeStep t1, ElcircTimeStep * t2)
{
    return ELIO_OK;
}

/*!
 * Define a header for the bottom. Not implemented.
 */
int ElioGetBottomHeader(ElcircHeader * h1, ElcircHeader * h2)
{
    return ELIO_OK;
}

/*! Not implemented.
 */
int ElioGetBottomStep(ElcircHeader * h1, ElcircHeader * h2, ElcircTimeStep t1, ElcircTimeStep * t2)
{
    return ELIO_OK;
}

/*! Not implemented.
 */
int ElioGetTransectHeader(ElcircHeader * h1, int n, double *x, double *y, ElcircHeader * h2)
{
    return ELIO_OK;
}

/*! Not implemented.
 */
int ElioGetTransectStep(ElcircHeader * h1, ElcircHeader * h2, int step, int n, double *x, double *y, ElcircTimeStep t1, ElcircTimeStep * t2)
{
    return ELIO_OK;
}

/*! Not implemented.
 */
int ElioGetLevelHeader(ElcircHeader * h1, ElcircHeader * h2)
{
    return ELIO_OK;
}

/*! Not implemented.
 */
int ElioGetLevelStep(ElcircHeader * h1, ElcircHeader * h2, ElcircTimeStep t1, int level, ElcircTimeStep * t2)
{
    return ELIO_OK;
}

/*! Not implemented.
 */
int ElioGetZLevelHeader(ElcircHeader * h1, ElcircHeader * h2)
{
    return ELIO_OK;
}

/*! Not implemented.
 */
int ElioGetZLevelStep(ElcircHeader * h1, ElcircHeader * h2, ElcircTimeStep t1, float z, ElcircTimeStep * t2)
{
    return ELIO_OK;
}

/*! Write the ELCIRC header.
 *  \param fp Pointer to open file handle.
 *  \param h Header struct to write.
 *
 *  The structure to write should be prepared in advance with an allocated
 *  grid and all members filled in with appropriate values.
 */
int ElioPutHeader(FILE * fp, ElcircHeader * h)
{
    int *itmp, i, j, bitmp, etmp;
    if (fp == NULL) {
	fprintf(stderr, "Bad file pointer\n");
	return ELIO_ERR;
    }
    /*if (!(strncmp(h->magic, "DataFormat v3.0", 15) == 0)) {
       fprintf(stderr, "ElioPutHeader(): Version identifier %s, must be 'DataFormat v3.0'\n", h->magic);
       return (ELIO_ERR);
       } */
    if (fwrite(h->magic, sizeof(char), 48, fp) != 48)
	return ELIO_FWRITE_ERR;
    h->magic[48] = 0;
    if (fseeko(fp, 48L, SEEK_SET)) {
	return ELIO_FSEEK_ERR;
    }
    if (fwrite(h->version, sizeof(char), 48, fp) != 48)
	return ELIO_FWRITE_ERR;
    h->version[48] = 0;
    if (fwrite(h->start_time, sizeof(char), 48, fp) != 48)
	return ELIO_FWRITE_ERR;
    h->start_time[48] = 0;
    if (fwrite(h->variable_nm, sizeof(char), 48, fp) != 48)
	return ELIO_FWRITE_ERR;
    h->variable_nm[48] = 0;
    if (fwrite(h->variable_dim, sizeof(char), 48, fp) != 48)
	return ELIO_FWRITE_ERR;
    h->variable_dim[48] = 0;
    if (fwrite(&h->nsteps, sizeof(int), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    if (fwrite(&h->timestep, sizeof(float), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    if (fwrite(&h->skip, sizeof(int), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    if (fwrite(&h->ivs, sizeof(int), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    if (fwrite(&h->i23d, sizeof(int), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    if (h->v == 3) {
	if (fwrite(&h->vpos, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->zmsl, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->nvrt, sizeof(int), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
    } else if (h->v == 4) {
	if (fwrite(&h->ivcor, sizeof(int), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->h0, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->hc, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->thetab, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->thetaf, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->nvrt, sizeof(int), 1, fp) != 1)
	    return ELIO_FREAD_ERR;
    } else if (h->v == 5) {
	if (fwrite(&h->nvrt, sizeof(int), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->kz, sizeof(int), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->h0, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->hs, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->hc, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->thetab, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(&h->thetaf, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
    }
    if (fwrite(h->zcor, sizeof(float), h->nvrt, fp) != h->nvrt)
	return ELIO_FWRITE_ERR;
    if (fwrite(&h->np, sizeof(int), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    if (fwrite(&h->ne, sizeof(int), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    for (i = 0; i < h->np; i++) {
	if (fwrite(h->x + i, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(h->y + i, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (fwrite(h->d + i, sizeof(float), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	if (h->v != 4) {
	    bitmp = h->bisave[i];	/* subtracted 1 previously, have to add it back */
	    if (fwrite(&bitmp, sizeof(int), 1, fp) != 1)
		return ELIO_FWRITE_ERR;
	}
    }
    /* for hybrid grids */
    for (i = 0; i < h->ne; i++) {
	if (fwrite(&h->etype[i], sizeof(int), 1, fp) != 1)
	    return ELIO_FWRITE_ERR;
	for (j = 0; j < h->etype[i]; j++) {
	    etmp = h->icon[j][i] + 1;	/* node number start from 0 internally, need to add 1 */
	    if (fwrite(&etmp, sizeof(int), 1, fp) != 1)
		return ELIO_FWRITE_ERR;
	}
    }
    return ELIO_OK;
}

/*!
 * Write a time step.
 *
 * \param fp Open filehandle.
 * \param h Pointer to ELCIRC header.
 * \param t Value parameter containing time step information to write.
 */
int ElioPutTimeStep(FILE * fp, ElcircHeader * h, ElcircTimeStep t)
{
    if (fwrite(&t.t, sizeof(float), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    if (fwrite(&t.it, sizeof(int), 1, fp) != 1)
	return ELIO_FWRITE_ERR;
    if (h->v != 4) {
	if (fwrite(t.surfind, sizeof(int), h->np, fp) != h->np)
	    return ELIO_FWRITE_ERR;
    } else {
	if (fwrite(t.e, sizeof(float), h->np, fp) != h->np)
	    return ELIO_FWRITE_ERR;
    }
    if (fwrite(t.d, sizeof(float), h->nitems, fp) != h->nitems)
	return ELIO_FWRITE_ERR;
    //printf("Completed Write of Time Step: Time, Step = %f, %d\n", t.t, t.it);
    return ELIO_OK;
}

/*! Not implemented
 *
 */
int ElioPutHeaderOld(FILE * fp, ElcircHeader * h)
{
    return ELIO_OK;
}

/*! Not implemented
 *
 */
int ElioPutTimeStepOld(FILE * fp, ElcircHeader * h, ElcircTimeStep t)
{
    return ELIO_OK;
}

/*! Return the file type.
 *
 * Returns -1 on error.
*/
int ElioGetFileType(FILE * fp)
{
    char id[100];
    int retval = -1;
    if (fread(id, sizeof(char), 15, fp) != 15) {
	return retval;
    }
    id[15] = 0;
    if (fseeko(fp, 0L, SEEK_SET)) {
	return retval;
    }
    if (strcmp(id, "DataFormat v2.0") == 0) {
	retval = 2;
    } else if (strcmp(id, "DataFormat v3.0") == 0) {
	retval = 3;
    } else if (strcmp(id, "DataFormat v4.0") == 0) {
	retval = 4;
    } else if (strcmp(id, "DataFormat v5.0") == 0) {
	retval = 5;
    }
    return retval;
}

/*!
 * Copy header h1 to h2
 */
int ElioCopyHeader(ElcircHeader * h1, ElcircHeader * h2)
{
    FILE *fp;
    int *itmp, i, j, ss;
    *h2 = *h1;
    memcpy(h2->magic, h1->magic, 48);
    memcpy(h2->version, h1->version, 48);
    memcpy(h2->start_time, h1->start_time, 48);
    memcpy(h2->variable_nm, h1->variable_nm, 48);
    memcpy(h2->variable_dim, h1->variable_dim, 48);
    h2->nsteps = h1->nsteps;
    h2->timestep = h1->timestep;
    h2->skip = h1->skip;
    h2->ivs = h1->ivs;
    h2->i23d = h1->i23d;
    h2->vpos = h1->vpos;
    h2->sigma = h1->sigma;
    h2->zmsl = h1->zmsl;
    h2->nvrt = h1->nvrt;
    h2->zcor = (float *) malloc(h1->nvrt * sizeof(float));
    for (i = 0; i < h1->nvrt; i++) {
	h2->zcor[i] = h2->zcor[i];
    }
    h2->np = h1->np;
    h2->ne = h1->ne;
    h2->x = (float *) malloc(h1->np * sizeof(float));
    h2->y = (float *) malloc(h1->np * sizeof(float));
    h2->d = (float *) malloc(h1->np * sizeof(float));
    h2->bi = (int *) malloc(h1->np * sizeof(int));
    h2->bisave = (int *) malloc(h1->np * sizeof(int));
    h2->no = (int *) malloc(h1->np * sizeof(int));
    memcpy(h2->x, h1->x, h1->np * sizeof(float));
    memcpy(h2->y, h1->y, h1->np * sizeof(float));
    memcpy(h2->d, h1->d, h1->np * sizeof(float));
    memcpy(h2->bi, h1->bi, h1->np * sizeof(int));
    memcpy(h2->bisave, h1->bisave, h1->np * sizeof(int));
    h2->etype = (int *) malloc(h1->ne * sizeof(int));
    memcpy(h2->etype, h1->etype, h1->ne * sizeof(int));
    h2->icon[0] = (int *) malloc(h1->ne * sizeof(int));
    h2->icon[1] = (int *) malloc(h1->ne * sizeof(int));
    h2->icon[2] = (int *) malloc(h1->ne * sizeof(int));
    h2->icon[3] = (int *) malloc(h1->ne * sizeof(int));
    memcpy(h2->icon[0], h1->icon[0], h1->ne * sizeof(int));
    memcpy(h2->icon[1], h1->icon[1], h1->ne * sizeof(int));
    memcpy(h2->icon[2], h1->icon[2], h1->ne * sizeof(int));
    if (h1->icon[3] != NULL) {
	memcpy(h2->icon[3], h1->icon[3], h1->ne * sizeof(int));
    } else {
	h2->icon[3] = NULL;
    }
    h2->nitems = h1->nitems;
    memcpy(h2->bi, h1->bi, h1->np * sizeof(int));
    memcpy(h2->bisave, h1->bisave, h1->np * sizeof(int));
    memcpy(h2->no, h1->no, h1->np * sizeof(int));
    h2->hsize = h1->hsize;
    h2->ssize = h1->ssize;
    return 0;
}

/*!
 * Write out data to a single step .61 file
 */
int ElioPut61File(char *fn, ElcircHeader * h, double *v)
{
    int i;
    FILE *fp;
    ElcircTimeStep t;

    if ((fp = fopen(fn, "w")) == NULL) {
        return ELIO_ERR;
    }
    if (ElioPutHeader(fp, h)) {
        return ELIO_ERR;
    }
    ElioAllocateTimeStep(h, &t);
    t.t = 0.0;
    t.it = 1;
    for (i = 0; i < h->np; i++) {
        t.d[i] = v[i];
    }
    ElioPutTimeStep(fp, h, t);
    ElioFreeTimeStep(&t);
    fclose(fp);
    return ELIO_OK;
}

/*!
 * Return the area for the given element. Element numbering starts from 0.
 *
 *  \param h - Previously read header.
 *  \param elem - Element number from 0 to compute the center.
 *
 *  Returns 0.0 on error.
 */
double ElioGetElementArea(ElcircHeader * h, int elem)
{
    int n1, n2, n3, n4;
    double ar = 0.0;
    double x1, y1, x2, y2, x3, y3, x4, y4;

    if (elem >= 0 && elem < h->ne) {
	n1 = h->icon[0][elem];
	n2 = h->icon[1][elem];
	n3 = h->icon[2][elem];
	x1 = h->x[n1];
	x2 = h->x[n2];
	x3 = h->x[n3];
	y1 = h->y[n1];
	y2 = h->y[n2];
	y3 = h->y[n3];
	if (h->etype[elem] == 3) {
	    ar = 0.5 * ((x1 * y2 - x2 * y1) + (x2 * y3 - x3 * y2) + (x3 * y1 - x1 * y3));
	} else if (h->etype[elem] == 4) {
	    n4 = h->icon[3][elem];
	    x4 = h->x[n4];
	    y4 = h->y[n4];
	    ar = 0.5 * ((x1 * y2 - x2 * y1) + (x2 * y3 - x3 * y2) + (x3 * y4 - x4 * y3) + (x4 * y1 - x1 * y4));
	}
    }
    return ar;
}

/*!
 * Return the center of the element 'elem' in (cx, cy).
 *
 *  \param h - Previously read header.
 *  \param elem - Element number from 0 to compute the center.
 *  \param *cx - Returns x.
 *  \param *cy - Returns y.
 *
 *  Returns ELIO_OK on success.
 *  Returns ELIO_ERR on failure.
 */
int ElioGetElementCenter(ElcircHeader * h, int elem, double *cx, double *cy)
{
    int i, n;
    if (!h || h->ne <= 0 || elem < 0 || (h->etype[elem] != 3 && h->etype[elem] != 4)) {
	return ELIO_ERR;
    }
    *cx = *cy = 0.0;
    for (i = 0; i < h->etype[elem]; i++) {
	n = h->icon[i][elem];
	*cx += h->x[n];
	*cy += h->y[n];
    }
    *cx /= h->etype[elem];
    *cy /= h->etype[elem];
    return ELIO_OK;
}

/*!
 * Return the center of the element 'elem' in (cx, cy).
 *
 *  \param g - Previously read grid.
 *  \param elem - Element number from 0 to compute the center.
 *  \param *cx - Returns x.
 *  \param *cy - Returns y.
 *
 *  Returns ELIO_OK on success.
 *  Returns ELIO_ERR on failure.
 */
int ElioGetGridElementCenter(ElioGrid * g, int elem, double *cx, double *cy)
{
    int i, n;
    if (!g || g->ne <= 0 || elem < 0 || (g->etype[elem] != 3 && g->etype[elem] != 4)) {
	return ELIO_ERR;
    }
    *cx = *cy = 0.0;
    for (i = 0; i < g->etype[elem]; i++) {
	n = g->icon[i][elem];
	*cx += g->x[n];
	*cy += g->y[n];
    }
    *cx /= g->etype[elem];
    *cy /= g->etype[elem];
    return ELIO_OK;
}

/*!
 * Return the area for the given element. Element numbering starts from 0.
 *
 *  Returns 0.0 on error.
 */
double ElioGetGridElementArea(ElioGrid * g, int elem)
{
    int n1, n2, n3, n4;
    double ar = 0.0;
    double x1, y1, x2, y2, x3, y3, x4, y4;

    if (elem >= 0 && elem < g->ne) {
	n1 = g->icon[0][elem];
	n2 = g->icon[1][elem];
	n3 = g->icon[2][elem];
	x1 = g->x[n1];
	x2 = g->x[n2];
	x3 = g->x[n3];
	y1 = g->y[n1];
	y2 = g->y[n2];
	y3 = g->y[n3];
	if (g->etype[elem] == 3) {
	    ar = 0.5 * ((x1 * y2 - x2 * y1) + (x2 * y3 - x3 * y2) + (x3 * y1 - x1 * y3));
	} else if (g->etype[elem] == 4) {
	    n4 = g->icon[3][elem];
	    x4 = g->x[n4];
	    y4 = g->y[n4];
	    ar = 0.5 * ((x1 * y2 - x2 * y1) + (x2 * y3 - x3 * y2) + (x3 * y4 - x4 * y3) + (x4 * y1 - x1 * y4));
	}
    }
    return ar;
}

/*!
 * Read a grid
 *
 *  \param gname - Grid file name
 *  \param g - Pointer to a grid structure
 *
 *  Returns ELIO_OK on success.
 *  Returns ELIO_ERR on failure.
 */
int ElioReadGrid(char *gname, ElioGrid * g)
{
    FILE *fp;
    char buf[1024];
    int i;
    if ((fp = fopen(gname, "rb")) == NULL) {
	return ELIO_ERR;
    }
    fgets(buf, 255, fp);
    fgets(buf, 255, fp);
    sscanf(buf, "%d %d", &g->ne, &g->np);
    g->x = (double *) malloc(g->np * sizeof(double));
    g->y = (double *) malloc(g->np * sizeof(double));
    g->d = (double *) malloc(g->np * sizeof(double));
    g->etype = (int *) malloc(g->ne * sizeof(int));
    for (i = 0; i < 4; i++) {
	g->icon[i] = (int *) malloc(g->ne * sizeof(int));
    }
    for (i = 0; i < g->np; i++) {
	fgets(buf, 255, fp);
	sscanf(buf, "%*d %lf %lf %lf", &g->x[i], &g->y[i], &g->d[i]);
    }
    for (i = 0; i < g->ne; i++) {
	fgets(buf, 255, fp);
	sscanf(buf, "%*d %d %d %d %d %d", &g->etype[i], &g->icon[0][i], &g->icon[1][i], &g->icon[2][i], &g->icon[3][i]);
	g->icon[0][i]--;
	g->icon[1][i]--;
	g->icon[2][i]--;
	if (g->etype[i] == 3) {
	    g->icon[3][i] = -1;
	} else {
	    g->icon[3][i]--;
	}
    }
/* could be a boundary in this grid file, try to read it
//3 = Number of open boundaries
//91 = Total number of open boundary nodes
//85 = Number of nodes for open boundary 1

    if (fgets(buf, 255, fp) != NULL) {
	if (sscanf(buf, "%d", &g->nopen) == 1) {
            if (fgets(buf, 255, fp) != NULL) {
	        if (sscanf(buf, "%d", &g->ntotalopen) == 1) {
		    for (i=0;i<g->nopen;i++) {
			if (sscanf(buf, "%d %d", &g->nopen) == 1) {
	            	    if (sscanf(buf, "%d", &g->ntotalopen) == 1) {
	            	    }
			}
		    }
	        }
	    }
	}
    } */
    fclose(fp);
    return ELIO_OK;
}

/*!
 * Free a grid
 *
 *  \param g - Pointer to a grid structure with already allocated memory
 *
 *  Returns ELIO_OK on success.
 *
 *  No checking is done on validity of the pointers being freed
 */
int ElioFreeGrid(ElioGrid *g)
{
    int i;
    free(g->x);
    free(g->y);
    free(g->d);
    free(g->etype);
    for (i = 0; i < 4; i++) {
        free(g->icon[i]);
    }
    return ELIO_OK;
}

int ElioFindElementInGrid(ElioGrid * g, double xp, double yp)
{
    int i, j;
    int n0, n1, n2, n3;
    double x[4], y[4];
    for (i = 0; i < g->ne; i++) {
	for (j = 0; j < g->etype[i]; j++) {
	    x[j] = g->x[g->icon[j][i]];
	    y[j] = g->y[g->icon[j][i]];
	}
	//printf("ElioInsideElement %lf %lf %lf %lf %lf %lf %lf %lf\n", xp, yp, x[0], y[0], x[1], y[1], x[2], y[2]);
	if (g->etype[i] == 3) {
	    if (ElioInsideElement(xp, yp, x[0], y[0], x[1], y[1], x[2], y[2])) {
		return i;
	    }
	} else {
	    if (ElioInsideElement4(xp, yp, x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3])) {
		return i;
	    }
	}
    }
    return -1;
}

int ElioGridFindNearestNode(ElioGrid * g, double x, double y)
{
    int i, ind = -1;
    double radius;
    double tmp, *xord = g->x, *yord = g->y;

    if (g->np > 0) {
	radius = hypot((x - xord[0]), (y - yord[0]));
	ind = 0;
	for (i = 1; i < g->np; i++) {
	    tmp = hypot((x - xord[i]), (y - yord[i]));
	    if (tmp < radius) {
		radius = tmp;
		ind = i;
	    }
	}
    } else {
	ind = -1;
    }
    return ind;
}

int ElioFindNearestNode(ElcircHeader * h, double x, double y)
{
    int i, ind;
    double radius;
    double tmp;
    float *xord = h->x, *yord = h->y;

    if (h->np > 0) {
	radius = hypot((x - xord[0]), (y - yord[0]));
	ind = 0;
	for (i = 1; i < h->np; i++) {
	    tmp = hypot((x - xord[i]), (y - yord[i]));
	    if (tmp < radius) {
		radius = tmp;
		ind = i;
	    }
	}
    } else {
	ind = -1;
    }
    return ind;
}

int ElioWriteGrid(char *gname, ElioGrid * g)
{
    FILE *fp;
    char buf[1024];
    int i;
    if ((fp = fopen(gname, "wb")) == NULL) {
	return ELIO_ERR;
    }
    fprintf(fp, "Grid\n");
    fprintf(fp, "%d %d\n", g->ne, g->np);
    for (i = 0; i < g->np; i++) {
	fprintf(fp, "%d %lf %lf %lf\n", i + 1, g->x[i], g->y[i], g->d[i]);
    }
    for (i = 0; i < g->ne; i++) {
	if (g->etype[i] == 3) {
	    fprintf(fp, "%d 3 %d %d %d\n", i + 1, g->icon[0][i] + 1, g->icon[1][i] + 1, g->icon[2][i] + 1);
	} else if (g->etype[i] == 4) {
	    fprintf(fp, "%d 4 %d %d %d %d\n", i + 1, g->icon[0][i] + 1, g->icon[1][i] + 1, g->icon[2][i] + 1, g->icon[3][i] + 1);
	}
    }
    fclose(fp);
    return ELIO_OK;
}

/*!
 * Time/date functions.
 */
static int corietime = 0;
static double coriestart = 0;

/*!
 * Turn on CORIE time, should be called just once in any program.
 * Multiple calls will foul things up, so don't do it.
 */
void ElioSetCorieTime(int b)
{
    corietime = b;
    if (b) {
	putenv("TZ=PST");
	coriestart = ElioGetDay(12, 31, 1995, 0, 0, 0);
    } else {
	coriestart = 0.0;
    }
}

/*!
 * Set corietime from a string. 4 digits for year.
 */
void ElioSetStartTime(int b, char *t)
{
    int y, m, d, h, mi, s;
    corietime = b;
    if (b) {
	sscanf(t, "%4d%2d%2d%2d%2d%2d", &y, &m, &d, &h, &mi, &s);
	coriestart = ElioGetDay(m, d, y, h, mi, s);
    } else {
	coriestart = 0.0;
    }
}

/*!
 * Return the corietime constant. If non-zero then corie time is active.
 */
int ElioGetCorieDayConstant(void)
{
    return corietime;
}

/*!
 * Return corieday relative to now. The variable 'when' can be +-.
 */
double ElioGetCorieDay(int when)
{
    double jd = 1;
    time_t t;
    t = time(NULL) + when * 86400;
    jd = ElioGetDay(12, 31, 1995, 0, 0, 0);
    jd = t / 86400.0 - jd;
    return jd;
}

/*!
 *  Return the CORIE day or UNIX day + hms as a real number.
 */
double ElioGetDay(int mon, int day, int year, int h, int mi, double se)
{
    double frac = (int) se - se;
    struct tm t;
    time_t j;
    t.tm_sec = se;
    t.tm_min = mi;
    t.tm_hour = h;
    t.tm_mday = day;
    t.tm_mon = mon - 1;
    t.tm_year = year - 1900;
    j = mktime(&t);
    return ((double) j / 86400.0 - coriestart);
}

/*!
 *  Convert a CORIE or UNIX day number to m-d-y hhmmss.
 */
void ElioGetDate(double jd, int *m, int *d, int *y, int *h, int *mi, double *sec)
{
    struct tm *t;
    time_t seconds = (jd + coriestart) * 86400.0;
    t = gmtime(&seconds);
    *sec = t->tm_sec;
    *mi = t->tm_min;
    *h = t->tm_hour;
    *d = t->tm_mday;
    *m = t->tm_mon + 1;
    *y = t->tm_year + 1900;
}

/*!
 *  Convert a day number to y-yearday hhmmss.
 */
void ElioGetYearDay(double jd, int *y, int *yd, int *h, int *mi, double *sec)
{
    struct tm *t;
    time_t seconds = (jd + coriestart) * 86400.0;
    t = gmtime(&seconds);
    *sec = t->tm_sec;
    *mi = t->tm_min;
    *h = t->tm_hour;
    *yd = t->tm_yday + 1;
    *y = t->tm_year + 1900;
}

/*!
 * Inverse bilinear mapping for quadrangles from Joseph.
 * Convexity of the quad must have been checked, and (x,y) must 
 * be inside the quad.
 *
 */
int ibilinear(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4, double x, double y, double *xi, double *eta, double *shapef)
{
    double axi[2], aet[2], bxy[2], root_xi[2], root_et[2];
    double x0, y0, dxi, deta, dd, beta, delta, dgamma;
    int i, j, icount, icaseno;

    static double SMALL = 1.0e-5;

//!...  Consts.
    x0 = (x1 + x2 + x3 + x4) / 4.0;
    y0 = (y1 + y2 + y3 + y4) / 4.0;
    axi[0] = x2 - x1 + x3 - x4;
    axi[1] = y2 - y1 + y3 - y4;
    aet[0] = x3 + x4 - x1 - x2;
    aet[1] = y3 + y4 - y1 - y2;
    bxy[0] = x1 - x2 + x3 - x4;
    bxy[1] = y1 - y2 + y3 - y4;
    dxi = 2 * ((x3 - x4) * (y1 - y2) - (y3 - y4) * (x1 - x2));
    deta = 2 * ((x4 - x1) * (y3 - y2) - (y4 - y1) * (x3 - x2));

//!...  Inverse mapping
    if ((fabs(bxy[0]) < SMALL && fabs(bxy[1]) < SMALL)
	|| (fabs(dxi) < SMALL && fabs(deta) < SMALL)) {
	icaseno = 1;
//    if (dxi == 0.0 && deta == 0.0) {
	dd = axi[0] * aet[1] - axi[1] * aet[0];
	if (dd == 0.0) {
	    fprintf(stderr, "Case 1 error: %lf\n", dd);
	    return ELIO_ERR;
	}
	*xi = 4 * (aet[1] * (x - x0) - aet[0] * (y - y0)) / dd;
	*eta = 4 * (axi[0] * (y - y0) - axi[1] * (x - x0)) / dd;

    } else if (fabs(dxi) < SMALL && fabs(deta) >= SMALL) {
//    } else if (dxi == 0 && deta != 0) {
	icaseno = 2;
	*eta = 4 * (bxy[1] * (x - x0) - bxy[1] * (y - y0)) / deta;
	dd = (axi[0] + *eta * bxy[0]) * (axi[0] + *eta * bxy[0]) + (axi[1] + *eta * bxy[1]) * (axi[1] + *eta * bxy[1]);
	if (dd == 0) {
	    fprintf(stderr, "Case 2 error: &lf\n", dd);
	    return ELIO_ERR;
	}
	*xi = ((4 * (x - x0) - *eta * aet[0]) * (axi[0] + *eta * bxy[0]) + (4 * (y - y0) - *eta * aet[1]) * (axi[1] + *eta * bxy[1])) / dd;
	;
    } else if (fabs(dxi) >= SMALL && fabs(deta) < SMALL) {
	icaseno = 3;
//    } else if (dxi != 0 && deta == 0) {
	*xi = 4 * (bxy[1] * (x - x0) - bxy[0] * (y - y0)) / dxi;
	dd = (aet[0] + *xi * bxy[0]) * (aet[0] + *xi * bxy[0]) + (aet[1] + *xi * bxy[1]) * (aet[1] + *xi * bxy[1]);
	if (dd == 0) {
	    fprintf(stderr, "Case 3 error: &lf\n", dd);
	    return ELIO_ERR;
	}
	*eta = ((4 * (x - x0) - *xi * axi[0]) * (aet[0] + *xi * bxy[0]) + (4 * (y - y0) - *xi * axi[1]) * (aet[1] + *xi * bxy[1])) / dd;
	;
    } else {
	icaseno = 4;
	beta = aet[1] * axi[0] - aet[0] * axi[1] - 4 * (bxy[1] * (x - x0) - bxy[0] * (y - y0));
	dgamma = 4 * (aet[0] * (y - y0) - aet[1] * (x - x0));
	delta = beta * beta - 4 * dgamma * dxi;
	if (delta == 0) {
	    *xi = (-beta / 2.0 / dxi);
	    *eta = (4 * (bxy[1] * (x - x0) - bxy[0] * (y - y0)) - *xi * dxi) / deta;
	} else if (delta > 0.0) {
	    root_xi[0] = (-beta + sqrt(delta)) / 2 / dxi;
	    root_xi[1] = (-beta - sqrt(delta)) / 2 / dxi;
	    icount = 0;
	    for (i = 0; i < 2; i++) {
		root_et[i] = (4 * (bxy[1] * (x - x0) - bxy[0] * (y - y0)) - root_xi[i] * dxi) / deta;
		if (fabs(root_xi[i]) <= 1.1 && fabs(root_et[i]) <= 1.1) {
		    *xi = root_xi[i];
		    *eta = root_et[i];
		    icount = icount + 1;
		}
	    }
	    if (icount == 2 && fabs(root_xi[0] - root_xi[1]) < SMALL) {
	    } else if (icount != 1) {
		fprintf(stderr, "Abnormal instances %lf %lf %d\n", root_xi[0], root_et[1], icount);
		return ELIO_ERR;
	    }
	} else {
	    fprintf(stderr, "No roots %lf\n", delta);
	    return ELIO_ERR;
	}
    }

    if (fabs(*xi) > 1.1 || fabs(*eta) > 1.1) {
	fprintf(stderr, "Out of bound in ibilinear %lf %lf\n", *xi, *eta);
	return ELIO_ERR;
    }

    *xi = (*xi > -1.0) ? *xi : -1.0;
    *eta = (*eta > -1.0) ? *eta : -1.0;
    *xi = (*xi < 1.0) ? *xi : 1.0;
    *eta = (*eta < 1.0) ? *eta : 1.0;
    //xi=dmin1(1.0,dmax1(xi,-1.0));
    //eta=dmin1(1.0,dmax1(eta,-1.0));

    shapef[0] = (1 - *xi) * (1 - *eta) / 4;
    shapef[1] = (1 + *xi) * (1 - *eta) / 4;
    shapef[2] = (1 + *xi) * (1 + *eta) / 4;
    shapef[3] = (1 - *xi) * (1 + *eta) / 4;
    return ELIO_OK;
}

static double C(double sigma, ElcircHeader * h)
{
    return (1 - h->thetab) * sinh(h->thetaf * sigma) / sinh(h->thetaf) + h->thetab * (tanh(h->thetaf * (sigma + 0.5)) - tanh(h->thetaf * 0.5)) / (2.0 * tanh(h->thetaf * 0.5));
}

/*!
 *  Compute the depth at a given node and level for a sigma grid.
 *
 * Input variables:
 *
 * \param node = Node number starting from 0
 * \param level = Level number starting from 0
 * \param eh = Pointer to populated ElcircHeader struct
 *
 * Returns depth.
 *
 * Fortran code from ELCIRC
 *
 * do k=1,nvrt
 *     z(i,k)=peta(i)*(1+sigma(k))+h_c*sigma(k)+(dp(i)-h_c)*cs(k)
 * !           Abnormal cases
 *     if(ivcor==2.and.(dp(i)<=h_c.or.peta(i)<=-h_c-(dp(i)-h_c)*theta_f/dsinh(theta_f))) then
 *        iback(i)=0
 *        z(i,k)=peta(i)*(1+sigma(k))+dp(i)*sigma(k)
 *     endif
 *     if(k>=2.and.z(i,k)-z(i,k-1)<=0) then
 *        write(11,*)'Inverted z-levels at:',i,k,z(i,k)-z(i,k-1),peta(i),dp(i)
 *        stop
 *     endif
 * enddo !k
 *
 */
double ElioGetSigmaDepthAtNode(int node, int level, float eta, ElcircHeader * eh)
{
    double z, sigma, h;
    sigma = eh->zcor[level];
    h = eh->d[node];
    if (eh->ivcor == 1) {
	z = (h + eta) * sigma + eta;
    } else if (eh->ivcor == 2) {
	if ((eta + h) > eh->h0 && h > eh->hc && eta > (-eh->hc - (h - eh->hc) * eh->thetaf / sinh(eh->thetaf))) {
	    z = eta * (1.0 + sigma) + eh->hc * eta + (h - eh->hc) * C(sigma, eh);
	} else {
	    z = (h + eta) * sigma + eta;
	}
    }
    return z;
}

int ElioGetSigmaDepth(int node, float eta, ElcircHeader * eh, double *depths)
{
    double z, sigma, h;
    int i;

    for (i = 0; i < eh->nvrt; i++) {
	depths[i] = ElioGetSigmaDepthAtNode(node, i, eta, eh);
    }
    return ELIO_OK;
}

/*!
 *  Compute the depth at a given XY location and level for a sigma grid.
 *
 * Input variables:
 *
 * \param x = X co-ordinate of location
 * \param y = Y co-ordinate of location
 * \param level = Level number starting from 0
 * \param eh = Pointer to populated ElcircHeader struct
 *
 * Returns depth.
 *
 */
double ElioGetSigmaDepthAtXY(double x, double y, int level, float e, ElcircHeader * eh)
{
    double H, h, eta, sigma;
    /* Need to interpolate elevation and depth */
    //H = eh->d[node] + e;
    h = H * eh->zcor[level] + e;
    return h;
}

/*!
 *  Compute the depth at a given XYZ location for a sigma grid.
 *
 * Input variables:
 *
 * \param x = X co-ordinate of location
 * \param y = Y co-ordinate of location
 * \param z = z co-ordinate of location
 * \param eh = Pointer to populated ElcircHeader struct
 *
 * Returns depth.
 *
 */
double ElioGetSigmaDepthAtXYZ(double x, double y, double z, float e, ElcircHeader * eh)
{
    double H, h, eta, sigma;
    //H = eh->d[node] + e;
    //h = H * eh->zcor[level] + e;
    return h;
}

/*!
 *  Return the min/max value in an array of doubles, utility function.
 *
 * Input variables:
 *
 * \param n = Number of doubles.
 * \param ar = Pointer to the doubles - assumed to be a valid pointer.
 *
 * Returns minimum double found in the array.
 *
 */
int ElioMinMax(int n, double *ar, double *dmin, double *dmax)
{
    int i;
    if (ar == (double *) NULL) {
	return ELIO_ERR;
    }
    *dmax = *dmin = ar[0];
    for (i = 1; i < n; i++) {
	if (*dmin > ar[i]) {
	    *dmin = ar[i];
	}
	if (*dmax < ar[i]) {
	    *dmax = ar[i];
	}
    }
    return ELIO_OK;
}

/*!
 *  Return the minimum value in an array of ints, utility function.
 *
 * Input variables:
 *
 * \param n = Number of integers.
 * \param iar = Pointer to the integers - assumed to be a valid pointer.
 *
 * Returns minimum integer found in the array.
 *
 */
int ElioIntMin(int n, int *iar)
{
    int i, im = 0;
    im = iar[0];
    for (i = 1; i < n; i++) {
	if (im > iar[i]) {
	    im = iar[i];
	}
    }
    return im;
}

/*!
 *  Return the maximum value in an array of ints, utility function.
 *
 * Input variables:
 *
 * \param n = Number of integers.
 * \param iar = Pointer to the integers - assumed to be a valid pointer.
 *
 * Returns maximum integer found in the array.
 *
 */
int ElioIntMax(int n, int *iar)
{
    int i, im = 0;
    im = iar[0];
    for (i = 1; i < n; i++) {
	if (im < iar[i]) {
	    im = iar[i];
	}
    }
    return im;
}

/*!
 *  Return the numeric file version number 2, 3, 4
 *
 * Input variables:
 *
 * \param fname = File name.
 *
 */
int ElioGetFileVersion(char *fname)
{
    FILE *fp;
    char magic[50];
    fp = fopen(fname, "rb");
    if (fp == NULL) {
	fprintf(stderr, "ElioGetFileVersion(): Unable to open file %s\n", fname);
	return -1;
    }
    if (fread(magic, sizeof(char), 48, fp) != 48)
	return ELIO_FREAD_ERR;
    magic[48] = 0;
    /* DataFormat v3.0 */
    if (strncmp(magic, "DataFormat v4.0", 15) == 0) {
	return 4;
    } else if (strncmp(magic, "DataFormat v5.0", 15) == 0) {
	return 5;
    } else if (strncmp(magic, "DataFormat v3.0", 15) == 0) {
	return 3;
    } else if (strncmp(magic, "DataFormat v2.0", 15) == 0) {
	return 2;
    } else {
	return -1;
    }
    fclose(fp);
}

/*!
 *  Find the interval in an array containing the double d, return -1
 *  if no interval found. Data point is located in between index returned
 *  and index + 1.
 *
 * Input variables:
 *
 * \param n = Number of items in *data.
 * \param data = Array to search.
 * \param d = Item to locate.
 *
 * Returns 0 based index found or -1 if less than or n if greater than.
 *
 */
int ElioFindIndex(int n, double *data, double d)
{
    int i;
    if (d < data[0]) {
	return -1;
    }
    if (d > data[n - 1]) {
	return n;
    }
    for (i = 0; i < n - 1; i++) {
	if (d >= data[i] && d <= data[i + 1]) {
	    return i;
	}
    }
/* Not reached */
    return -1;
}

/*!
 * Given an array of depths and scalar quatities, interpolate
 * the scalar values to all the points in *zp. Clamp to values at 
 * either end if point to find is out of range. Output to *sp.
 *
 * Input variables:
 *
 * \param n = Number of items in *z.
 * \param z = Depths.
 * \param s = Scalar value.
 * \param np = number of points to interpolate.
 * \param zp = Depths to interpolate to.
 * \param sp = Interpolated scalar values.
 */
void ElioInterpolateArray(int n, double *z, double *s, int np, double *zp, double *sp)
{
    int i, j, l1, l2;
    double tmp;
    for (j = 0; j < np; j++) {
/* Use minimum value */
	if (zp[j] < z[0]) {
	    sp[j] = s[0];
/* Use maximum value */
	} else if (zp[j] > z[n - 1]) {
	    sp[j] = s[n - 1];
/* Use interpolated value */
	} else {
	    for (i = 0; i < n - 1; i++) {
		l1 = i;
		l2 = i + 1;
		if (zp[j] >= z[l1] && zp[j] <= z[l2]) {
		    break;
		}
	    }
	    tmp = (z[l2] - z[l1]);
	    if (tmp != 0.0) {
		sp[j] = (zp[j] - z[l1]) / (z[l2] - z[l1]) * (s[l2] - s[l1]) + s[l1];
	    } else {
		fprintf(stderr, "Elio: In ElioInterpolateArray(), z[l2] - z[l1] == 0.0\n");
		sp[j] = -9999.0;
	    }
	}
    }
}

/*!
 * Given an array of depths and scalar quantities, interpolate
 * a scalar value given the depth. Clamp to value at 
 * either end if point to find is out of range. returns the
 * interplated value.
 *
 * Input variables:
 *
 * \param n = Number of items in *z.
 * \param z = Depths.
 * \param s = Scalar value.
 * \param zp = Depth to interpolate to.
 */
double ElioInterpolate(int n, double *z, double *s, double zp)
{
    int i, l1, l2;
    double sp = 0.0;
    double tmp;
/* Use minimum value */
    if (zp < z[0]) {
	sp = s[0];
/* Use maximum value */
    } else if (zp > z[n - 1]) {
	sp = s[n - 1];
/* Use interpolated value */
    } else {
	for (i = 0; i < n - 1; i++) {
	    l1 = i;
	    l2 = i + 1;
	    if (zp >= z[l1] && zp <= z[l2]) {
		break;
	    }
	}
	tmp = (z[l2] - z[l1]);
	if (tmp != 0.0) {
	    sp = (zp - z[l1]) / (z[l2] - z[l1]) * (s[l2] - s[l1]) + s[l1];
	} else {
	    fprintf(stderr, "Elio: In ElioInterpolate(), z[l2] - z[l1] == 0.0\n");
	    sp = -9999.0;
	}
    }
    return sp;
}

/*!
 * Given an array of depths and scalar quantities, interpolate
 * a scalar value given the index as found from ElioFindIndex(). 
 * Clamp to value at either end if point to find is out of range.
 * 
 * Input variables:
 * 
 * \param n = Number of items in *z.
 * \param z = Depths.
 * \param s = Scalar value.
 * \param ind = Index value.
 * \param zp = Depth to use for interpolation.
 */
double ElioInterpolateAtIndex(int n, double *z, double *s, int ind, double zp)
{
    int i, l1, l2;
    double sp = 0.0;
    double tmp;
/* Use minimum value */
    if (ind < 0) {
	sp = s[0];
/* Use maximum value */
    } else if (ind == n) {
	sp = s[n - 1];
/* Use interpolated value */
    } else {
	l1 = ind;
	l2 = ind + 1;
	tmp = (z[l2] - z[l1]);
	if (tmp != 0.0) {
	    sp = (zp - z[l1]) / (z[l2] - z[l1]) * (s[l2] - s[l1]) + s[l1];
	} else {
	    //fprintf(stderr, "Elio: In ElioInterpolateAtIndex(), z[l2] - z[l1] == 0.0\n");
	    sp = -9999.0;
	}
    }
    return sp;
}

/*!
 * Given an array of depths and scalar quantities, interpolate
 * a scalar value given the index as found from ElioFindIndex(). 
 * Clamp to value at either end if point to find is out of range.
 * 
 * Input variables:
 * 
 * \param n = Number of items in *z.
 * \param z = Depths.
 * \param s = Scalar value.
 * \param ind = Index value.
 * \param zp = Depth to use for interpolation.
 */
int ElioMakeDepthsOld(ElcircHeader * h, float *d, float *dd)
{
    int i, j, k;
    for (i = 0; i < h->np * h->nvrt; i++) {
	dd[i] = -99.0;
    }
    for (j = 0; j < h->np; j++) {
// Fill in missing depth values in file - only for Z levels.
	for (k = 0; k < h->kz; k++) {
	    dd[j * h->nvrt + k] = h->zcor[k];
	}
// Fill in depth values from zcor.
	for (k = h->bi[j]; k < h->nvrt; k++) {
	    dd[j * h->nvrt + k] = d[h->no[j] + k - h->bi[j]];
	}
    }
    return ELIO_OK;
}

/*!
 * Replace the nodes in the grid h with the nodes found in grid g.
 *
 * Useful for transforming a run done in one projection to another.
 *
 */
int ElioReplaceElcircGrid(ElcircHeader * h, ElioGrid * g)
{
    int i;
    if (h->np != g->np) {
	return ELIO_ERR;
    }
    for (i = 0; i < g->np; i++) {
	h->x[i] = g->x[i];
	h->y[i] = g->y[i];
    }
    return ELIO_OK;
}

/*!
 * Get the grid from a header and stuff it in an ElioGrid struture
 *
 */
int ElioGetGrid(ElcircHeader * h, ElioGrid * g)
{
    int i;
    g->ne = h->ne;
    g->np = h->np;
    g->x = (double *) malloc(g->np * sizeof(double));
    g->y = (double *) malloc(g->np * sizeof(double));
    g->d = (double *) malloc(g->np * sizeof(double));
    g->etype = (int *) malloc(g->ne * sizeof(int));
    for (i = 0; i < 4; i++) {
	g->icon[i] = (int *) malloc(g->ne * sizeof(int));
    }
    for (i = 0; i < g->np; i++) {
	g->x[i] = h->x[i];
	g->y[i] = h->y[i];
	g->d[i] = h->d[i];
    }
    for (i = 0; i < g->ne; i++) {
	g->etype[i] = h->etype[i];
	g->icon[0][i] = h->icon[0][i];
	g->icon[1][i] = h->icon[1][i];
	g->icon[2][i] = h->icon[2][i];
	if (g->etype[i] == 4) {
	    g->icon[3][i] = h->icon[3][i];
	}
    }
    return ELIO_OK;
}