/* Vis5D version 4.3 */

/* this should be updated when the file version changes */
#define FILE_VERSION "4.2"

/*
 * New grid file format for VIS-5D:
 *
 * The header is a list of tagged items.  Each item has 3 parts:
 *    1. A tag which is a 4-byte integer identifying the type of item.
 *    2. A 4-byte integer indicating how many bytes of data follow.
 *    3. The binary data.
 *
 * If we need to add new information to a file header we just create a
 * new tag and add the code to read/write the information.
 *
 * If we're reading a header and find an unknown tag, we can use the
 * length field to skip ahead to the next tag.  Therefore, the file
 * format is forward (and backward) compatible.
 *
 * Grid data is stored as either:
 *     1-byte unsigned integers  (255=missing)
 *     2-byte unsigned integers  (65535=missing)
 *     4-byte IEEE floats        ( >1.0e30 = missing)
 *
 * All numeric values are stored in big endian order.  All floating point
 * values are in IEEE format.
 */

/*
 * Updates:
 *
 * April 13, 1995, brianp
 *   finished Cray support for 2-byte and 4-byte compress modes
 */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <math.h>
#include "v5d43.h"
#ifndef SEEK_SET
#  define SEEK_SET 0
#endif
#ifndef SEEK_CUR
#  define SEEK_CUR 1
#endif
#ifndef SEEK_END
#  define SEEK_END 2
#endif

/*
 * Currently defined tags:
 * Note:  the notation a[i] doesn't mean a is an array of i elements,
 * rather it just refers to the ith element of a[].
 *
 * Tags marked as PHASED OUT should be readable but are no longer written.
 * Old tag numbers can't be reused!
 *
 */

/*      TAG NAME        VALUE       DATA (comments)                     */
/*----------------------------------------------------------------------*/
#define TAG_ID          0x5635440a  /* hex encoding of "V5D\n"          */

/* general stuff 1000+ */
#define TAG_VERSION     1000        /* char*10 FileVersion              */
#define TAG_NUMTIMES    1001        /* int*4 NumTimes                   */
#define TAG_NUMVARS     1002        /* int*4 NumVars                    */
#define TAG_VARNAME     1003        /* int*4 var; char*10 VarName[var]  */

#define TAG_NR          1004        /* int*4 Nr                         */ 
#define TAG_NC          1005        /* int*4 Nc                         */
#define TAG_NL          1006        /* int*4 Nl  (Nl for all vars)      */
#define TAG_NL_VAR      1007        /* int*4 var; int*4 Nl[var]         */
#define TAG_LOWLEV_VAR  1008        /* int*4 var; int*4 LowLev[var]     */

#define TAG_TIME        1010        /* int*4 t;  int*4 TimeStamp[t]     */
#define TAG_DATE        1011        /* int*4 t;  int*4 DateStamp[t]     */

#define TAG_MINVAL      1012        /* int*4 var;  real*4 MinVal[var]   */
#define TAG_MAXVAL      1013        /* int*4 var;  real*4 MaxVal[var]   */

#define TAG_COMPRESS    1014        /* int*4 CompressMode; (#bytes/grid)*/

#define TAG_UNITS       1015        /* int *4 var; char*20 Units[var]   */

/* vertical coordinate system 2000+ */
#define TAG_VERTICAL_SYSTEM 2000    /* int*4 VerticalSystem             */
#define TAG_VERT_ARGS    2100       /* int*4 n;  real*4 VertArgs[0..n-1]*/

#define TAG_BOTTOMBOUND  2001       /* real*4 BottomBound     (PHASED OUT)  */
#define TAG_LEVINC       2002       /* real*4 LevInc      (PHASED OUT)      */
#define TAG_HEIGHT       2003    /* int*4 l;  real*4 Height[l] (PHASED OUT) */

/* projection 3000+ */
#define TAG_PROJECTION   3000       /* int*4 projection:                    */
                                    /*   0 = generic linear                 */
                                    /*   1 = cylindrical equidistant        */
                                    /*   2 = Lambert conformal/Polar Stereo */
                                    /*   3 = rotated equidistant            */
#define TAG_PROJ_ARGS    3100       /* int *4 n;  real*4 ProjArgs[0..n-1]   */

#define TAG_NORTHBOUND   3001       /* real*4 NorthBound       (PHASED OUT) */
#define TAG_WESTBOUND    3002       /* real*4 WestBound        (PHASED OUT) */
#define TAG_ROWINC       3003       /* real*4 RowInc           (PHASED OUT) */
#define TAG_COLINC       3004       /* real*4 ColInc           (PHASED OUT) */
#define TAG_LAT1         3005       /* real*4 Lat1             (PHASED OUT) */
#define TAG_LAT2         3006       /* real*4 Lat2             (PHASED OUT) */
#define TAG_POLE_ROW     3007       /* real*4 PoleRow          (PHASED OUT) */
#define TAG_POLE_COL     3008       /* real*4 PoleCol          (PHASED OUT) */
#define TAG_CENTLON      3009       /* real*4 CentralLon       (PHASED OUT) */
#define TAG_CENTLAT      3010       /* real*4 CentralLat       (PHASED OUT) */
#define TAG_CENTROW      3011       /* real*4 CentralRow       (PHASED OUT) */
#define TAG_CENTCOL      3012       /* real*4 CentralCol       (PHASED OUT) */
#define TAG_ROTATION     3013       /* real*4 Rotation         (PHASED OUT) */

#define TAG_END                9999

/**********************************************************************/
/*****                  Miscellaneous Functions                   *****/
/**********************************************************************/

float pressure_to_height(float pressure)
{
  return (float) DEFAULT_LOG_EXP * log((double) pressure / DEFAULT_LOG_SCALE);
}

float height_to_pressure(float height)
{
  return (float) DEFAULT_LOG_SCALE * exp((double) height / DEFAULT_LOG_EXP);
}

/*
 * Return current file position.
 * Input:  f - file descriptor
 */
static off_t ltell( int f )
{
   return lseek( f, 0, SEEK_CUR );
}

/* ******************************************************************
 * Copy up to maxlen characters from src to dst stopping upon whitespace
 * in src.  Terminate dst with null character.
 * Return:  length of dst.
 */
static int copy_string2( char *dst, const char *src, int maxlen )
{
   int i;

   for (i=0;i<maxlen;i++) dst[i] = src[i];
   for (i=maxlen-1; i>=0; i--) {
     if (dst[i]==' ' || i==maxlen-1) dst[i] = 0;
     else break;
   }
   return strlen(dst);
}

/* ******************************************************************
 * Copy up to maxlen characters from src to dst stopping upon whitespace
 * in src.  Terminate dst with null character.
 * Return:  length of dst.
 */
static int copy_string( char *dst, const char *src, int maxlen )
{
   int i;

   for (i=0;i<maxlen;i++) {
      if (src[i]==' ' || i==maxlen-1) {
         dst[i] = 0;
         break;
      }
      else {
         dst[i] = src[i];
      }
   }
   return i;
}

/* ******************************************************************
 * Convert a date from YYDDD format to days since Jan 1, 1900.
 */
int v5dYYDDDtoDays( int yyddd )
{
   int iy, id, idays;

   iy = yyddd / 1000;
   id = yyddd - 1000*iy;
   if (iy < 50) iy += 100; /* WLH 31 July 96 << 31 Dec 99 */
   idays = 365*iy + (iy-1)/4 + id;

   return idays;
}

/* ******************************************************************
 * Convert a time from HHMMSS format to seconds since midnight.
 */
int v5dHHMMSStoSeconds( int hhmmss )
{
   int h, m, s;

   h = hhmmss / 10000;
   m = (hhmmss / 100) % 100;
   s = hhmmss % 100;

   return s + m*60 + h*60*60;
}

/* ******************************************************************
 * Convert a day since Jan 1, 1900 to YYDDD format.
 */
int v5dDaysToYYDDD( int days )
{
   int iy, id, iyyddd;

   iy = (4*days)/1461;
   id = days-(365*iy+(iy-1)/4);
   if (iy > 99) iy = iy - 100; /* WLH 31 July 96 << 31 Dec 99 */
   /* iy = iy + 1900; is the right way to fix this, but requires
      changing all places where dates are printed - procrastinate */
   iyyddd = iy*1000+id;

   return iyyddd;
}

/* ******************************************************************
 * Convert a time in seconds since midnight to HHMMSS format.
 */
int v5dSecondsToHHMMSS( int seconds )
{
   int hh, mm, ss;

   hh = seconds / (60*60);
   mm = (seconds / 60) % 60;
   ss = seconds % 60;
   return hh*10000 + mm * 100 + ss;
}

void v5dPrintStruct( const v5dstruct *v )
{
   static char day[7][10] = { "Sunday", "Monday", "Tuesday", "Wednesday",
                              "Thursday", "Friday", "Saturday" };
   int time, var, i;
   int maxnl;

   maxnl = 0;
   for (var=0;var<v->NumVars;var++) {
      if (v->Nl[var]+v->LowLev[var]>maxnl) {
         maxnl = v->Nl[var]+v->LowLev[var];
      }
   }

   if (v->FileFormat==0) {
      if (v->FileVersion[0] == 0) {
        printf("File format: v5d  version: (4.0 or 4.1)\n");
      }
      else {
        printf("File format: v5d  version: %s\n", v->FileVersion);
      }
   }
   else {
      printf("File format: comp5d  (VIS-5D 3.3 or older)\n");
   }

   if (v->CompressMode==1) {
      printf("Compression:  1 byte per gridpoint.\n");
   }
   else {
      printf("Compression:  %d bytes per gridpoint.\n", v->CompressMode);
   }
   printf("header size=%d\n", v->FirstGridPos);
   printf("sizeof(v5dstruct)=%d\n", sizeof(v5dstruct) );
   printf("\n");

   printf("NumVars = %d\n", v->NumVars );

   printf("Var  Name  Units  Rows  Cols  Levels LowLev  MinVal  MaxVal\n");
   for (var=0;var<v->NumVars;var++) {
      printf("%3d  %-10s %-10s %3d   %3d   %3d    %3d",
             var+1, v->VarName[var], v->Units[var],
             v->Nr, v->Nc, v->Nl[var], v->LowLev[var] );
      if (v->MinVal[var] > v->MaxVal[var]) {
         printf("     MISSING      MISSING\n");
      }
      else {
         printf("     %-12g %-12g\n", v->MinVal[var], v->MaxVal[var] );
      }
   }

   printf("\n");

   printf("NumTimes = %d\n", v->NumTimes );
   printf("Step    Date(YYDDD)    Time(HH:MM:SS)   Day\n");
   for (time=0;time<v->NumTimes;time++) {
      int i = v->TimeStamp[time];
      printf("%3d        %05d       %5d:%02d:%02d     %s\n",
             time+1,
             v->DateStamp[time],
             i/10000, (i/100)%100, i%100,
             day[ v5dYYDDDtoDays(v->DateStamp[time]) % 7 ]);
   }
   printf("\n");

   switch (v->VerticalSystem) {
      case 0:
         printf("Generic linear vertical coordinate system:\n");
         printf("\tBottom Bound: %f\n", v->VertArgs[0] );
         printf("\tIncrement between levels:  %f\n", v->VertArgs[1] );
         break;
      case 1:
         printf("Equally spaced levels in km:\n");
         printf("\tBottom Bound: %f\n", v->VertArgs[0] );
         printf("\tIncrement: %f\n", v->VertArgs[1] );
         break;
      case 2:
         printf("Unequally spaced levels in km:\n");
         printf("Level\tHeight(km)\n");
         for (i=0;i<maxnl;i++) {
            printf("%3d     %10.3f\n", i+1, v->VertArgs[i] );
         }
         break;
      case 3:
         printf("Unequally spaced levels in mb:\n");
         printf("Level\tPressure(mb)\n");
         for (i=0;i<maxnl;i++) {
            printf("%3d     %10.3f\n", i+1, height_to_pressure(v->VertArgs[i]) );
         }
         break;
      default:
         printf("Bad VerticalSystem value: %d\n", v->VerticalSystem );
   }
   printf("\n");

   switch (v->Projection) {
      case 0:
         printf("Generic linear projection:\n");
         printf("\tNorth Boundary: %f\n", v->ProjArgs[0] );
         printf("\tWest Boundary: %f\n", v->ProjArgs[1] );
         printf("\tRow Increment: %f\n", v->ProjArgs[2] );
         printf("\tColumn Increment: %f\n", v->ProjArgs[3] );
         break;
      case 1:
         printf("Cylindrical Equidistant projection:\n");
         printf("\tNorth Boundary: %f degrees\n", v->ProjArgs[0] );
         printf("\tWest Boundary: %f degrees\n", v->ProjArgs[1] );
         printf("\tRow Increment: %f degrees\n", v->ProjArgs[2] );
         printf("\tColumn Increment: %f degrees\n", v->ProjArgs[3] );
/*
         printf("\tSouth Boundary: %f degrees\n",
                v->NorthBound - v->RowInc * (v->Nr-1) );
         printf("\tEast Boundary: %f degrees\n",
                v->WestBound - v->ColInc * (v->Nc-1) );
*/
         break;
      case 2:
         printf("Lambert Conformal projection:\n");
         printf("\tStandard Latitude 1: %f\n", v->ProjArgs[0] );
         printf("\tStandard Latitude 2: %f\n", v->ProjArgs[1] );
         printf("\tNorth/South Pole Row: %f\n", v->ProjArgs[2] );
         printf("\tNorth/South Pole Column: %f\n", v->ProjArgs[3] );
         printf("\tCentral Longitude: %f\n", v->ProjArgs[4] );
         printf("\tColumn Increment: %f km\n", v->ProjArgs[5] );
         break;
      case 3:
         printf("Stereographic:\n");
         printf("\tCenter Latitude: %f\n", v->ProjArgs[0] );
         printf("\tCenter Longitude: %f\n", v->ProjArgs[1] );
         printf("\tCenter Row: %f\n", v->ProjArgs[2] );
         printf("\tCenter Column: %f\n", v->ProjArgs[3] );
         printf("\tColumn Spacing: %f\n", v->ProjArgs[4] );
         break;
      case 4:
         /* WLH 4-21-95 */
         printf("Rotated equidistant projection:\n");
         printf("\tLatitude of grid(0,0): %f\n", v->ProjArgs[0] );
         printf("\tLongitude of grid(0,0): %f\n", v->ProjArgs[1] );
         printf("\tRow Increment: %f degress\n", v->ProjArgs[2] );
         printf("\tColumn Increment: %f degrees\n", v->ProjArgs[3] );
         printf("\tCenter Latitude: %f\n", v->ProjArgs[4] );
         printf("\tCenter Longitude: %f\n", v->ProjArgs[5] );
         printf("\tRotation: %f degrees\n", v->ProjArgs[6] );
         break;
      default:
         printf("Bad projection number: %d\n", v->Projection );
   }
}

/* ******************************************************************
 * Compute the location of a compressed grid within a file.
 * Input:  v - pointer to v5dstruct describing the file header.
 *         time, var - which timestep and variable.
 * Return:  file offset in bytes
 */
static int grid_position( const v5dstruct *v, int time, int var )
{
   int pos, i;

   assert( time >= 0 );
   assert( var >= 0 );
   assert( time < v->NumTimes );
   assert( var < v->NumVars );

   pos = v->FirstGridPos + time * v->SumGridSizes;
   for (i=0;i<var;i++) {
      pos += v->GridSize[i];
   }

   return pos;
}

/* ******************************************************************
 * Compute the ga and gb (de)compression values for a grid.
 * Input:  nr, nc, nl - size of grid
 *         data - the grid data
 *         ga, gb - arrays to store results.
 *         minval, maxval - pointer to floats to return min, max values
 *         compressmode - 1, 2 or 4 bytes per grid point
 * Output:  ga, gb - the (de)compression values
 *          minval, maxval - the min and max grid values
 * Side effect:  the MinVal[var] and MaxVal[var] fields in g may be
 *               updated with new values.
 */
static void compute_ga_gb( int nr, int nc, int nl,
                           const float data[], int compressmode,
                           float ga[], float gb[],
                           float *minval, float *maxval )
{
#ifdef SIMPLE_COMPRESSION
   /*
    * Compute ga, gb values for whole grid.
    */
   int i, lev, allmissing, num;
   float min, max, a, b;

   min = 1.0e30;
   max = -1.0e30;
   num = nr * nc * nl;
   allmissing = 1;
   for (i=0;i<num;i++) {
      if (!IS_MISSING(data[i])) {
         if (data[i]<min)  min = data[i];
         if (data[i]>max)  max = data[i];
         allmissing = 0;
      }
   }
   if (allmissing) {
      a = 1.0;
      b = 0.0;
   }
   else {
      a = (max-min) / 254.0;
      b = min;
   }

   /* return results */
   for (i=0;i<nl;i++) {
      ga[i] = a;
      gb[i] = b;
   }

   *minval = min;
   *maxval = max;
#else
   /*
    * Compress grid on level-by-level basis.
    */
#  define SMALLVALUE -1.0e30
#  define BIGVALUE 1.0e30
#  define ABS(x)   ( ((x) < 0.0) ? -(x) : (x) )
   float gridmin, gridmax;
   float levmin[MAXLEVELS], levmax[MAXLEVELS];
   float d[MAXLEVELS], dmax;
   float ival, mval;
   int j, k, lev, nrnc;

   nrnc = nr * nc;

   /* find min and max for each layer and the whole grid */
   gridmin = BIGVALUE;
   gridmax = SMALLVALUE;
   j = 0;
   for (lev=0;lev<nl;lev++) {
      float min, max;
      min = BIGVALUE;
      max = SMALLVALUE;
      for (k=0;k<nrnc;k++) {
         if (!IS_MISSING(data[j]) && data[j]<min)
            min = data[j];
         if (!IS_MISSING(data[j]) && data[j]>max)
            max = data[j];
         j++;
      }
      if (min<gridmin)
        gridmin = min;
      if (max>gridmax)
        gridmax = max;
      levmin[lev] = min;
      levmax[lev] = max;
   }

/* WLH 2-2-95 */
#ifdef KLUDGE
   /* if the grid minimum is within delt of 0.0, fudge all values */
   /* within delt of 0.0 to delt, and recalculate mins and maxes */
   {
      float delt;
      int nrncnl = nrnc * nl;

      delt = (gridmax - gridmin)/100000.0;
      if ( ABS(gridmin) < delt && gridmin!=0.0 && compressmode != 4 ) {
         float min, max;
         for (j=0; j<nrncnl; j++) {
            if (!IS_MISSING(data[j]) && data[j]<delt)
              data[j] = delt;
         }
         /* re-calculate min and max for each layer and the whole grid */
         gridmin = delt;
         for (lev=0;lev<nl;lev++) {
            if (ABS(levmin[lev]) < delt)
              levmin[lev] = delt;
            if (ABS(levmax[lev]) < delt)
              levmax[lev] = delt;
         }
      }
   }
#endif

   /* find d[lev] and dmax = MAX( d[0], d[1], ... d[nl-1] ) */
   dmax = 0.0;
   for (lev=0;lev<nl;lev++) {
      if (levmin[lev]>=BIGVALUE && levmax[lev]<=SMALLVALUE) {
         /* all values in the layer are MISSING */
         d[lev] = 0.0;
      }
      else {
         d[lev] = levmax[lev]-levmin[lev];
      }
      if (d[lev]>dmax)
         dmax = d[lev];
   }

   /*** Compute ga (scale) and gb (bias) for each grid level */
   if (dmax==0.0) {
      /*** Special cases ***/
      if (gridmin==gridmax) {
         /*** whole grid is of same value ***/
         for (lev=0; lev<nl; lev++) {
            ga[lev] = gridmin;
            gb[lev] = 0.0;
         }
      }
      else {
         /*** every layer is of a single value ***/
         for (lev=0; lev<nl; lev++) {
            ga[lev] = levmin[lev];
            gb[lev] = 0.0;
         }
      }
   }
   else {
      /*** Normal cases ***/
      if (compressmode == 1) {
#define ORIGINAL
#ifdef ORIGINAL
         ival = dmax / 254.0;
         mval = gridmin;
         for (lev=0; lev<nl; lev++) {
            ga[lev] = ival;
            gb[lev] = mval + ival * (int) ( (levmin[lev]-mval) / ival ); 
         }
#else
         for (lev=0; lev<nl; lev++) {
            if (d[lev]==0.0) {
               ival = 1.0;
            }
            else {
               ival = d[lev] / 254.0;
            }
            ga[lev] = ival;
            gb[lev] = levmin[lev];
         }
#endif
      }
      else if (compressmode == 2) {
         ival = dmax / 65534.0;
         mval = gridmin;
         for (lev=0; lev<nl; lev++) {
            ga[lev] = ival;
            gb[lev] = mval + ival * (int) ( (levmin[lev]-mval) / ival );
         }
      }
      else {
         assert( compressmode==4 );
         for (lev=0; lev<nl; lev++) {
            ga[lev] = 1.0;
            gb[lev] = 0.0;
         }
      }
   }

   /* update min, max values */
   *minval = gridmin;
   *maxval = gridmax;
#endif
}

/* ******************************************************************
 * Compress a 3-D grid from floats to 1-byte unsigned integers.
 * Input: nr, nc, nl - size of grid
 *        compressmode - 1, 2 or 4 bytes per grid point
 *        data - array of [nr*nc*nl] floats
 *        compdata - pointer to array of [nr*nc*nl*compressmode] bytes
 *                   to put results into.
 *        ga, gb - pointer to arrays to put ga and gb decompression values
 *        minval, maxval - pointers to float to return min & max values
 * Output:  compdata - the compressed grid data
 *          ga, gb - the decompression values
 *          minval, maxval - the min and max grid values
 */
void v5dCompressGrid( int nr, int nc, int nl, int compressmode,
                      const float data[],
                      void *compdata, float ga[], float gb[],
                      float *minval, float *maxval )
{
   int nrnc = nr * nc;
   int nrncnl = nr * nc * nl;
   V5Dubyte *compdata1 = (V5Dubyte *) compdata;
   V5Dushort *compdata2 = (V5Dushort *) compdata;

   /* compute ga, gb values */
   compute_ga_gb( nr, nc, nl, data, compressmode, ga, gb, minval, maxval );

   /* compress the data */
   if (compressmode==1) {
      int i, lev, p;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float one_over_a, b;
         b = gb[lev] - 0.0001;  /* subtract an epsilon so the int((d-b)/a) */
                                /* expr below doesn't get mis-truncated. */
         if (ga[lev]==0.0) {
            one_over_a = 1.0;
         }
         else {
            one_over_a = 1.0 / ga[lev];
         }
         for (i=0;i<nrnc;i++,p++) {
            if (IS_MISSING(data[p])) {
               compdata1[p] = 255;
            }
            else {
               compdata1[p] = (V5Dubyte) (int) ((data[p]-b) * one_over_a);
            }
         }
      }
   }

   else if (compressmode == 2) {
      int i, lev, p;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float one_over_a, b;
         b = gb[lev] - 0.0001;
         if (ga[lev]==0.0) {
            one_over_a = 1.0;
         }
         else {
            one_over_a = 1.0 / ga[lev];
         }
#ifdef _CRAY
         /* this is tricky because sizeof(V5Dushort)==8, not 2 */
         for (i=0;i<nrnc;i++,p++) {
            V5Dushort compvalue;
            if (IS_MISSING(data[p])) {
               compvalue = 65535;
            }
            else {
               compvalue = (V5Dushort) (int) ((data[p]-b) * one_over_a);
            }
            compdata1[p*2+0] = compvalue >> 8;     /* upper byte */
            compdata1[p*2+1] = compvalue & 0xffu;  /* lower byte */
         }
#else
         for (i=0;i<nrnc;i++,p++) {
            if (IS_MISSING(data[p])) {
               compdata2[p] = 65535;
            }
            else {
               compdata2[p] = (V5Dushort) (int) ((data[p]-b) * one_over_a);
            }
         }
         /* TODO: byte-swapping on little endian??? */
#endif
      }
   }

   else {
      /* compressmode==4 */
#ifdef _CRAY
      cray_to_ieee_array( compdata, data, nrncnl );
#else
      /* other machines: just copy 4-byte IEEE floats */
      assert( sizeof(float)==4 );
      memcpy( compdata, data, nrncnl*4 );
      /* TODO: byte-swapping on little endian??? */
#endif
   }
}

/* ******************************************************************
 * Return the size (in bytes) of the 3-D grid specified by time and var.
 * Input:  v - pointer to v5dstruct describing the file
 *         time, var - which timestep and variable
 * Return:  number of data points.
 */
int v5dSizeofGrid( const v5dstruct *v, int time, int var )
{
   return v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
}

/* ******************************************************************
 * Initialize a v5dstructure to reasonable initial values.
 * Input:  v - pointer to v5dstruct.
 */
void v5dInitStruct( v5dstruct *v )
{
   int i;

   /* set everything to zero */
   memset( v, 0, sizeof(v5dstruct) );

   /* special cases */
   v->Projection = -1;
   v->VerticalSystem = -1;

   for (i=0;i<MAXVARS;i++) {
      v->MinVal[i] = MISSING;
      v->MaxVal[i] = -MISSING;
      v->LowLev[i] = 0;
   }

   /* set file version */
   strcpy(v->FileVersion, FILE_VERSION);

   v->CompressMode = 1;
   v->FileDesc = -1;
}

/* ******************************************************************
 * Return a pointer to a new, initialized v5dstruct.
 */
v5dstruct *v5dNewStruct( void )
{
   v5dstruct *v;

   v = (v5dstruct *) malloc( sizeof(v5dstruct) );
   if (v) {
      v5dInitStruct(v);
   }
   return v;
}

/* ******************************************************************
 * Free an initialized v5dstruct. (Todd Plessel)
 */
void v5dFreeStruct( v5dstruct* v )
{
   /*assert( v5dVerifyStruct( v ) );*/
   free( v );
   v = 0;
}

/* ******************************************************************
 * Do some checking that the information in a v5dstruct is valid.
 * Input:  v - pointer to v5dstruct
 * Return:  1 = g is ok, 0 = g is invalid
 */
int v5dVerifyStruct( const v5dstruct *v )
{
   int var, i, invalid, maxnl;

   invalid = 0;

   if (!v)
      return 0;

   /* Number of variables */
   if (v->NumVars<0) {
      printf("Invalid number of variables: %d\n", v->NumVars );
      invalid = 1;
   }
   else if (v->NumVars>MAXVARS) {
      printf("Too many variables: %d  (Maximum is %d)\n",
             v->NumVars, MAXVARS);
      invalid = 1;
   }

   /* Variable Names */
   for (i=0;i<v->NumVars;i++) {
      if (v->VarName[i][0]==0) {
         printf("Missing variable name: VarName[%d]=\"\"\n", i );
         invalid = 1;
      }
   }

   /* Number of timesteps */
   if (v->NumTimes<0) {
      printf("Invalid number of timesteps: %d\n", v->NumTimes );
      invalid = 1;
   }
   else if (v->NumTimes>MAXTIMES) {
      printf("Too many timesteps: %d  (Maximum is %d)\n",
             v->NumTimes, MAXTIMES );
      invalid = 1;
   }

   /* Make sure timestamps are increasing */
   for (i=1;i<v->NumTimes;i++) {
      int date0 = v5dYYDDDtoDays( v->DateStamp[i-1] );
      int date1 = v5dYYDDDtoDays( v->DateStamp[i] );
      int time0 = v5dHHMMSStoSeconds( v->TimeStamp[i-1] );
      int time1 = v5dHHMMSStoSeconds( v->TimeStamp[i] );
      if (time1<=time0 && date1<=date0) {
         printf("Timestamp for step %d must be later than step %d\n", i, i-1);
         invalid = 1;
      }
   }

   /* Rows */
   if (v->Nr<2) {
      printf("Too few rows: %d (2 is minimum)\n", v->Nr );
      invalid = 1;
   }
   else if (v->Nr>MAXROWS) {
      printf("Too many rows: %d (%d is maximum)\n", v->Nr, MAXROWS );
      invalid = 1;
   }

   /* Columns */
   if (v->Nc<2) {
      printf("Too few columns: %d (2 is minimum)\n", v->Nc );
      invalid = 1;
   }
   else if (v->Nc>MAXCOLUMNS) {
      printf("Too many columns: %d (%d is maximum)\n", v->Nc, MAXCOLUMNS );
      invalid = 1;
   }

   /* Levels */
   maxnl = 0;
   for (var=0;var<v->NumVars;var++) {
      if (v->LowLev[var] < 0) {
         printf("Low level cannot be negative for var %s: %d\n",
                 v->VarName[var], v->LowLev[var] );
         invalid = 1;
      }
      if (v->Nl[var]<1) {
         printf("Too few levels for var %s: %d (1 is minimum)\n",
                 v->VarName[var], v->Nl[var] );
         invalid = 1;
      }
      if (v->Nl[var]+v->LowLev[var]>MAXLEVELS) {
         printf("Too many levels for var %s: %d (%d is maximum)\n",
                 v->VarName[var], v->Nl[var]+v->LowLev[var], MAXLEVELS );
         invalid = 1;
      }
      if (v->Nl[var]+v->LowLev[var]>maxnl) {
         maxnl = v->Nl[var]+v->LowLev[var];
      }
   }

   if (v->CompressMode != 1 && v->CompressMode != 2 && v->CompressMode != 4) {
      printf("Bad CompressMode: %d (must be 1, 2 or 4)\n", v->CompressMode );
      invalid = 1;
   }

   switch (v->VerticalSystem) {
      case 0:
      case 1:
         if (v->VertArgs[1]==0.0) {
            printf("Vertical level increment is zero, must be non-zero\n");
            invalid = 1;
         }
         break;
      case 2:
         /* Check that Height values increase upward */
         for (i=1;i<maxnl;i++) {
            if (v->VertArgs[i] <= v->VertArgs[i-1]) {
      printf("Height[%d]=%f <= Height[%d]=%f, level heights must increase\n",
                      i, v->VertArgs[i], i-1, v->VertArgs[i-1] );
               invalid = 1;
               break;
            }
         }
         break;
      case 3:
         /* Check that Pressure values decrease upward */
         for (i=1;i<maxnl;i++) {
            if (v->VertArgs[i] <= v->VertArgs[i-1]) {
 printf("Pressure[%d]=%f >= Pressure[%d]=%f, level pressures must decrease\n",
                      i, height_to_pressure(v->VertArgs[i]),
                      i-1, height_to_pressure(v->VertArgs[i-1]) );
               invalid = 1;
               break;
            }
         }
         break;
      default:
         printf("VerticalSystem = %d, must be in 0..3\n", v->VerticalSystem );
         invalid = 1;
   }

   switch (v->Projection) {
      case 0:  /* Generic */
         if (v->ProjArgs[2]==0.0) {
            printf("Row Increment (ProjArgs[2]) can't be zero\n");
            invalid = 1;
         }
         if (v->ProjArgs[3]==0.0) {
            printf("Column increment (ProjArgs[3]) can't be zero\n");
            invalid = 1;
         }
         break;
      case 1:  /* Cylindrical equidistant */
         if (v->ProjArgs[2]<0.0) {
            printf("Row Increment (ProjArgs[2]) = %g  (must be >=0.0)\n",
                   v->ProjArgs[2] );
            invalid = 1;
         }
         if (v->ProjArgs[3]<=0.0) {
            printf("Column Increment (ProjArgs[3]) = %g  (must be >=0.0)\n",
                   v->ProjArgs[3] );
            invalid = 1;
         }
         break;
      case 2:  /* Lambert Conformal */
         if (v->ProjArgs[0]<-90.0 || v->ProjArgs[0]>90.0) {
            printf("Lat1 (ProjArgs[0]) out of range: %g\n", v->ProjArgs[0] );
            invalid = 1;
         }
         if (v->ProjArgs[1]<-90.0 || v->ProjArgs[1]>90.0) {
            printf("Lat2 (ProjArgs[1] out of range: %g\n", v->ProjArgs[1] );
            invalid = 1;
         }
         if (v->ProjArgs[5]<=0.0) {
            printf("ColInc (ProjArgs[5]) = %g  (must be >=0.0)\n",
                   v->ProjArgs[5] );
            invalid = 1;
         }
         break;
      case 3:  /* Stereographic */
         if (v->ProjArgs[0]<-90.0 || v->ProjArgs[0]>90.0) {
            printf("Central Latitude (ProjArgs[0]) out of range: ");
            printf("%g  (must be in +/-90)\n", v->ProjArgs[0] );
            invalid = 1;
         }
         if (v->ProjArgs[1]<-180.0 || v->ProjArgs[1]>180.0) {
            printf("Central Longitude (ProjArgs[1]) out of range: ");
            printf("%g  (must be in +/-180)\n", v->ProjArgs[1] );
            invalid = 1;
         }
         if (v->ProjArgs[4]<0) {
            printf("Column spacing (ProjArgs[4]) = %g  (must be positive)\n",
                   v->ProjArgs[4]);
            invalid = 1;
         }
         break;
      case 4:  /* Rotated */
         /* WLH 4-21-95 */
         if (v->ProjArgs[2]<=0.0) {
            printf("Row Increment (ProjArgs[2]) = %g  (must be >=0.0)\n",
                   v->ProjArgs[2] );
            invalid = 1;
         }
         if (v->ProjArgs[3]<=0.0) {
            printf("Column Increment = (ProjArgs[3]) %g  (must be >=0.0)\n",
                   v->ProjArgs[3] );
            invalid = 1;
         }
         if (v->ProjArgs[4]<-90.0 || v->ProjArgs[4]>90.0) {
            printf("Central Latitude (ProjArgs[4]) out of range: ");
            printf("%g  (must be in +/-90)\n", v->ProjArgs[4] );
            invalid = 1;
         }
         if (v->ProjArgs[5]<-180.0 || v->ProjArgs[5]>180.0) {
            printf("Central Longitude (ProjArgs[5]) out of range: ");
            printf("%g  (must be in +/-180)\n", v->ProjArgs[5] );
            invalid = 1;
         }
         if (v->ProjArgs[6]<-180.0 || v->ProjArgs[6]>180.0) {
            printf("Central Longitude (ProjArgs[6]) out of range: ");
            printf("%g  (must be in +/-180)\n", v->ProjArgs[6] );
            invalid = 1;
         }
         break;
      default:
         printf("Projection = %d, must be in 0..4\n", v->Projection );
         invalid = 1;
   }

   return !invalid;
}

/**********************************************************************/
/*****                   Output Functions                         *****/
/**********************************************************************/

static int write_tag( v5dstruct *v, int tag, int length, int newfile )
{
   if (!newfile) {
      /* have to check that there's room in header to write this tagged item */
      if (v->CurPos+8+length > v->FirstGridPos) {
         printf("Error: out of header space!\n");
         /* Out of header space! */
         return 0;
      }
   }

   if (write_int4( v->FileDesc, tag )==0)  return 0;
   if (write_int4( v->FileDesc, length )==0)  return 0;
   v->CurPos += 8 + length;
   return 1;
}

/* ******************************************************************
 * Write the information in the given v5dstruct as a v5d file header.
 * Note that the current file position is restored when this function
 * returns normally.
 * Input:  f - file already open for writing
 *         v - pointer to v5dstruct
 * Return:  1 = ok, 0 = error.
 */
static int write_v5d_header( v5dstruct *v )
{
   int var, time, filler, maxnl;
   int f;
   int newfile;

   if (v->FileFormat!=0) {
      printf("Error: v5d library can't write comp5d format files.\n");
      return 0;
   }

   f = v->FileDesc;

   if (!v5dVerifyStruct( v ))
      return 0;

   /* Determine if we're writing to a new file */
   if (v->FirstGridPos==0) {
      newfile = 1;
   }
   else {
      newfile = 0;
   }

   /* compute grid sizes */
   v->SumGridSizes = 0;
   for (var=0;var<v->NumVars;var++) {
      v->GridSize[var] = 8 * v->Nl[var] + v5dSizeofGrid( v, 0, var );
      v->SumGridSizes += v->GridSize[var];
   }

   /* set file pointer to start of file */
   lseek( f, 0, SEEK_SET );
   v->CurPos = 0;

   /*
    * Write the tagged header info
    */
#define WRITE_TAG( V, T, L )  if (!write_tag(V,T,L,newfile))  return 0;

   /* ID */
   WRITE_TAG( v, TAG_ID, 0 );

   /* File Version */
   WRITE_TAG( v, TAG_VERSION, 10 );
   write_bytes( f, FILE_VERSION, 10 );

   /* Number of timesteps */
   WRITE_TAG( v, TAG_NUMTIMES, 4 );
   write_int4( f, v->NumTimes );

   /* Number of variables */
   WRITE_TAG( v, TAG_NUMVARS, 4 );
   write_int4( f, v->NumVars );

   /* Names of variables */
   for (var=0;var<v->NumVars;var++) {
      WRITE_TAG( v, TAG_VARNAME, 14 );
      write_int4( f, var );
      write_bytes( f, v->VarName[var], 10 );
   }

   /* Physical Units */
   for (var=0;var<v->NumVars;var++) {
      WRITE_TAG( v, TAG_UNITS, 24 );
      write_int4( f, var );
      write_bytes( f, v->Units[var], 20 );
   }

   /* Date and time of each timestep */
   for (time=0;time<v->NumTimes;time++) {
      WRITE_TAG( v, TAG_TIME, 8 );
      write_int4( f, time );
      write_int4( f, v->TimeStamp[time] );
      WRITE_TAG( v, TAG_DATE, 8 );
      write_int4( f, time );
      write_int4( f, v->DateStamp[time] );
   }

   /* Number of rows */
   WRITE_TAG( v, TAG_NR, 4 );
   write_int4( f, v->Nr );

   /* Number of columns */
   WRITE_TAG( v, TAG_NC, 4 );
   write_int4( f, v->Nc );

   /* Number of levels, compute maxnl */
   maxnl = 0;
   for (var=0;var<v->NumVars;var++) {
      WRITE_TAG( v, TAG_NL_VAR, 8 );
      write_int4( f, var );
      write_int4( f, v->Nl[var] );
      WRITE_TAG( v, TAG_LOWLEV_VAR, 8 );
      write_int4( f, var );
      write_int4( f, v->LowLev[var] );
      if (v->Nl[var]+v->LowLev[var]>maxnl) {
         maxnl = v->Nl[var]+v->LowLev[var];
      }
   }

   /* Min/Max values */
   for (var=0;var<v->NumVars;var++) {
      WRITE_TAG( v, TAG_MINVAL, 8 );
      write_int4( f, var );
      write_float4( f, v->MinVal[var] );
      WRITE_TAG( v, TAG_MAXVAL, 8 );
      write_int4( f, var );
      write_float4( f, v->MaxVal[var] );
   }

   /* Compress mode */
   WRITE_TAG( v, TAG_COMPRESS, 4 );
   write_int4( f, v->CompressMode );

   /* Vertical Coordinate System */
   WRITE_TAG( v, TAG_VERTICAL_SYSTEM, 4 );
   write_int4( f, v->VerticalSystem );
   WRITE_TAG( v, TAG_VERT_ARGS, 4+4*MAXVERTARGS );
   write_int4( f, MAXVERTARGS );
   write_float4_array( f, v->VertArgs, MAXVERTARGS );

   /* Map Projection */
   WRITE_TAG( v, TAG_PROJECTION, 4 );
   write_int4( f, v->Projection );
   WRITE_TAG( v, TAG_PROJ_ARGS, 4+4*MAXPROJARGS );
   write_int4( f, MAXPROJARGS );
   write_float4_array( f, v->ProjArgs, MAXPROJARGS );

   /* write END tag */
   if (newfile) {
      /* We're writing to a brand new file.  Reserve 10000 bytes */
      /* for future header growth. */
      WRITE_TAG( v, TAG_END, 10000 );
      lseek( f, 10000, SEEK_CUR );

      /* Let file pointer indicate where first grid is stored */
      v->FirstGridPos = ltell( f );
   }
   else {
      /* we're rewriting a header */
      filler = v->FirstGridPos - ltell(f);
      WRITE_TAG( v, TAG_END, filler-8 );
   }

#undef WRITE_TAG

   return 1;
}

/* **************************************************************
 * Open a v5d file for writing.  If the named file already exists,
 * it will be deleted.
 * Input:  filename - name of v5d file to create.
 *         v - pointer to v5dstruct with the header info to write.
 * Return:  1 = ok, 0 = error.
 */
int v5dCreateFile( const char *filename, v5dstruct *v )
{
   mode_t mask;
   int fd;

   mask = 0666;
   fd = open( filename, O_WRONLY | O_CREAT | O_TRUNC, mask );
   if (fd==-1) {
      printf("Error in v5dCreateFile: open failed\n");
      v->FileDesc = -1;
      v->Mode = 0;
      return 0;
   }
   else {
      /* ok */
      v->FileDesc = fd;
      v->Mode = 'w';
      /* write header and return status */
      return write_v5d_header(v);
   }
}

/* **************************************************************
 * Write a compressed grid to a v5d file.
 * Input:  v - pointer to v5dstruct describing the file
 *         time, var - which timestep and variable
 *         ga, gb - the GA and GB (de)compression value arrays
 *         compdata - address of array of compressed data values
 * Return:  1 = ok, 0 = error.
 */
int v5dWriteCompressedGrid( const v5dstruct *v, int time, int var,
                            const float *ga, const float *gb,
                            const void *compdata )
{
   int pos, n, k;

   /* simple error checks */
   if (v->Mode!='w') {
      printf("Error in v5dWriteCompressedGrid: file opened for reading,");
      printf(" not writing.\n");
      return 0;
   }
   if (time<0 || time>=v->NumTimes) {
      printf("Error in v5dWriteCompressedGrid: bad timestep argument (%d)\n",
             time);
      return 0;
   }
   if (var<0 || var>=v->NumVars) {
      printf("Error in v5dWriteCompressedGrid: bad variable argument (%d)\n",
             var);
      return 0;
   }

   /* move to position in file */
   pos = grid_position( v, time, var );
   if (lseek( v->FileDesc, pos, SEEK_SET )<0) {
      /* lseek failed, return error */
      printf("Error in v5dWrite[Compressed]Grid: seek failed, disk full?\n");
      return 0;
   }

   /* write ga, gb arrays */
   k = 0;
   if (write_float4_array( v->FileDesc, ga, v->Nl[var] ) == v->Nl[var] &&
       write_float4_array( v->FileDesc, gb, v->Nl[var] ) == v->Nl[var]) {
      /* write compressed grid data (k=1=OK, k=0=Error) */
      n = v->Nr * v->Nc * v->Nl[var];
      if (v->CompressMode==1) {
         k = write_block( v->FileDesc, compdata, n, 1 )==n;
      }
      else if (v->CompressMode==2) {
         k = write_block( v->FileDesc, compdata, n, 2 )==n;
      }
      else if (v->CompressMode==4) {
         k = write_block( v->FileDesc, compdata, n, 4 )==n;
      }
   }

   if (k==0) {
      /* Error while writing */
      printf("Error in v5dWrite[Compressed]Grid: write failed, disk full?\n");
   }
   return k;

/*
   n = v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
   if (write_bytes( v->FileDesc, compdata, n )!=n) {
      printf("Error in v5dWrite[Compressed]Grid: write failed, disk full?\n");
      return 0;
   }
   else {
      return 1;
   }
*/
}

/* **************************************************************
 * Compress a grid and write it to a v5d file.
 * Input:  v - pointer to v5dstruct describing the file
 *         time, var - which timestep and variable (starting at 0)
 *         data - address of uncompressed grid data
 * Return:  1 = ok, 0 = error.
 */
int v5dWriteGrid( v5dstruct *v, int time, int var, const float data[] )
{
   float ga[MAXLEVELS], gb[MAXLEVELS];
   void *compdata;
   int n, bytes;
   float min, max;

   if (v->Mode!='w') {
      printf("Error in v5dWriteGrid: file opened for reading,");
      printf(" not writing.\n");
      return 0;
   }
   if (time<0 || time>=v->NumTimes) {
      printf("Error in v5dWriteGrid: bad timestep argument (%d)\n", time);
      return 0;
   }
   if (var<0 || var>=v->NumVars) {
      printf("Error in v5dWriteGrid: bad variable argument (%d)\n", var);
      return 0;
   }

   /* allocate compdata buffer */
   if (v->CompressMode==1) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned char);
   }
   else if (v->CompressMode==2) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned short);
   }
   else if (v->CompressMode==4) {
      bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(float);
   }
   compdata = (void *) malloc( bytes );
   if (!compdata) {
      printf("Error in v5dWriteGrid: out of memory (needed %d bytes)\n",
             bytes );
      return 0;
   }

   /* compress the grid data */
   v5dCompressGrid( v->Nr, v->Nc, v->Nl[var], v->CompressMode, data,
                    compdata, ga, gb, &min, &max );

   /* update min and max value */
   if (min<v->MinVal[var])
      v->MinVal[var] = min;
   if (max>v->MaxVal[var])
      v->MaxVal[var] = max;

   /* write the compressed grid */
   n = v5dWriteCompressedGrid( v, time, var, ga, gb, compdata );

   /* free compdata */
   free( compdata );

   return n;
}

/* **************************************************************
 * Close a v5d file which was opened with open_v5d_file() or
 * create_v5d_file().
 * Input: f - file descriptor
 * Return:  1 = ok, 0 = error
 */
int v5dCloseFile( v5dstruct *v )
{
   int status = 1;

   if (v->Mode=='w') {
      /* rewrite header because writing grids updates the minval and */
      /* maxval fields */
      lseek( v->FileDesc, 0, SEEK_SET );
      status = write_v5d_header( v );
      lseek( v->FileDesc, 0, SEEK_END );
      close( v->FileDesc );
   }
   else if (v->Mode=='r') {
      /* just close the file */
      close(v->FileDesc);
   }
   else {
      printf("Error in v5dCloseFile: bad v5dstruct argument\n");
      return 0;
   }
   v->FileDesc = -1;
   v->Mode = 0;
   return status;
}

/**********************************************************************/
/*****           Simple v5d file writing functions.               *****/
/**********************************************************************/

static v5dstruct *Simple = NULL;

/* **************************************************************
 * Create a new v5d file specifying both a map projection and vertical
 * coordinate system.  See README file for argument details.
 * Return:  1 = ok, 0 = error.
 */
int v5dCreate( const char *name, int numtimes, int numvars,
               int nr, int nc, const int nl[],
               const char varname[MAXVARS][10],
               const int timestamp[], const int datestamp[],
               int compressmode,
               int projection,
               const float proj_args[],
               int vertical,
               const float vert_args[] )
{
   int var, time, maxnl, i;

   /* initialize the v5dstruct */
   Simple = v5dNewStruct();

   Simple->NumTimes = numtimes;
   Simple->NumVars = numvars;
   Simple->Nr = nr;
   Simple->Nc = nc;
   maxnl = nl[0];
   for (var=0;var<numvars;var++) {
      if (nl[var]>maxnl) {
         maxnl = nl[var];
      }
      Simple->Nl[var] = nl[var];
      Simple->LowLev[var] = 0;
      strncpy( Simple->VarName[var], varname[var], 10 );
      Simple->VarName[var][9] = 0;
   }

   /* time and date for each timestep */
   for (time=0;time<numtimes;time++) {
      Simple->TimeStamp[time] = timestamp[time];
      Simple->DateStamp[time] = datestamp[time];
   }

   Simple->CompressMode = compressmode;

   /* Map projection and vertical coordinate system */
   Simple->Projection = projection;
   memcpy( Simple->ProjArgs, proj_args, MAXPROJARGS*sizeof(float) );

   Simple->VerticalSystem = vertical;
   if (vertical == 3) {
     /* convert pressures to heights */
     for (i=0; i<MAXVERTARGS; i++) {
       if (vert_args[i] > 0.000001) {
         Simple->VertArgs[i] = pressure_to_height(vert_args[i]);
       }
       else Simple->VertArgs[i] = 0.0;
     }
   }
   else {
     memcpy( Simple->VertArgs, vert_args, MAXVERTARGS*sizeof(float) );
   }

   /* create the file */
   if (v5dCreateFile( name, Simple )==0) {
     printf("Error in v5dCreateSimpleFile: unable to create %s\n", name );
     return 0;
   }
   else {
      return 1;
   }
}

/* **************************************************************
 * Set lowest levels for each variable (other than default of 0).
 * Input: lowlev - array [NumVars] of ints
 * Return:  1 = ok, 0 = error
 */
int v5dSetLowLev( int lowlev[] )
{
  int var;

  if (Simple) {
     for (var=0;var<Simple->NumVars;var++) {
        Simple->LowLev[var] = lowlev[var];
     }
     return 1;
  }
  else {
     printf("Error: must call v5dCreate before v5dSetLowLev\n");
     return 0;
  }
}

/* **************************************************************
 * Set the units for a variable.
 * Input:  var - a variable in [1,NumVars]
 *         units - a string
 * Return:  1 = ok, 0 = error
 */
int v5dSetUnits( int var, const char *units )
{
  if (Simple) {
     if (var>=1 && var<=Simple->NumVars) {
        strncpy( Simple->Units[var-1], units, 19 );
        Simple->Units[var-1][19] = 0;
        return 1;
     }
     else {
        printf("Error: bad variable number in v5dSetUnits\n");
        return 0;
     }
  }
  else {
     printf("Error: must call v5dCreate before v5dSetUnits\n");
     return 0;
  }
}

/* **************************************************************
 * Write a grid to a v5d file.
 * Input:  time - timestep in [1,NumTimes]
 *         var - timestep in [1,NumVars]
 *         data - array [nr*nc*nl] of floats
 * Return:  1 = ok, 0 = error
 */
int v5dWrite( int time, int var, const float data[] )
{
   if (Simple) {
      if (time<1 || time>Simple->NumTimes) {
         printf("Error in v5dWrite: bad timestep number: %d\n", time );
         return 0;
      }
      if (var<1 || var>Simple->NumVars) {
         printf("Error in v5dWrite: bad variable number: %d\n", var );
      }
      return v5dWriteGrid( Simple, time-1, var-1, data );
   }
   else {
      printf("Error: must call v5dCreate before v5dWrite\n");
      return 0;
   }
}

/* **************************************************************
 * Close a v5d file after the last grid has been written to it.
 * Return:  1 = ok, 0 = error
 */
int v5dClose( void )
{
   if (Simple) {
     int ok = v5dCloseFile( Simple );
     v5dFreeStruct( Simple );
     return ok;
   }
   else {
     printf("Error: v5dClose: no file to close\n");
     return 0;
   }
}

/**********************************************************************/
/*****                FORTRAN-callable simple output              *****/
/**********************************************************************/

/* **************************************************************
 * Create a v5d file.  See README file for argument descriptions.
 * Return:  1 = ok, 0 = error.
 */
#ifdef UNDERSCORE
   int v5dcreate_
#else
#  ifdef _CRAY
     int V5DCREATE
#  else
     int v5dcreate
#  endif
#endif
           ( const char *name, const int *numtimes, const int *numvars,
             const int *nr, const int *nc, const int nl[],
             const char varname[][10],
             const int timestamp[], const int datestamp[],
             const int *compressmode,
             const int *projection,
             const float proj_args[],
             const int *vertical,
             const float vert_args[] )
{
   char filename[13];
   char names[MAXVARS][10];
   int i, maxnl, args;

   /* copy name to filename and remove trailing spaces if any */
   copy_string( filename, name, 14 );

   /*
    * Check for uninitialized arguments
    */
   if (*numtimes<1) {
      printf("Error: numtimes invalid\n");
      return 0;
   }
   if (*numvars<1) {
      printf("Error: numvars invalid\n");
      return 0;
   }
   if (*nr<2) {
      printf("Error: nr invalid\n");
      return 0;
   }
   if (*nc<2) {
      printf("Error: nc invalid\n");
      return 0;
   }
   maxnl = 0;
   for (i=0;i<*numvars;i++) {
      if (nl[i]<1) {
         printf("Error: nl(%d) invalid\n", i+1);
         return 0;
      }
      if (nl[i]>maxnl) {
         maxnl = nl[i];
      }
   }

   for (i=0;i<*numvars;i++) {
      if (copy_string2( names[i], varname[i], 10)==0) {
         printf("Error: unitialized varname(%d)\n", i+1);
         return 0;
      }
   }

   for (i=0;i<*numtimes;i++) {
      if (timestamp[i]<0) {
         printf("Error: times(%d) invalid\n", i+1);
         return 0;
      }
      if (datestamp[i]<0) {
         printf("Error: dates(%d) invalid\n", i+1);
         return 0;
      }
   }

   if (*compressmode != 1 && *compressmode != 2 && *compressmode != 4) {
      printf("Error: compressmode invalid\n");
      return 0;
   }

   switch (*projection) {
      case 0:
         args = 4;
         break;
      case 1:
         args = 0;
         if (IS_MISSING(proj_args[0])) {
            printf("Error: northlat (proj_args(1)) invalid\n");
            return 0;
         }
         if (IS_MISSING(proj_args[1])) {
            printf("Error: westlon (proj_args(2)) invalid\n");
            return 0;
         }
         if (IS_MISSING(proj_args[2])) {
            printf("Error: latinc (proj_args(3)) invalid\n");
            return 0;
         }
         if (IS_MISSING(proj_args[3])) {
            printf("Error: loninc (proj_args(4)) invalid\n");
            return 0;
         }
         break;
      case 2:
         args = 6;
         break;
      case 3:
         args = 5;
         break;
      case 4:
         args = 7;
         break;
      default:
         args = 0;
         printf("Error: projection invalid\n");
         return 0;
   }
   for (i=0;i<args;i++) {
      if (IS_MISSING(proj_args[i])) {
         printf("Error: proj_args(%d) invalid\n", i+1);
         return 0;
      }
   }

   switch (*vertical) {
      case 0:
/* WLH 31 Oct 96  -  just fall through
         args = 4;
         break;
*/
      case 1:
         args = 0;
         if (IS_MISSING(vert_args[0])) {
            printf("Error: bottomhgt (vert_args(1)) invalid\n");
            return 0;
         }
         if (IS_MISSING(vert_args[1])) {
            printf("Error: hgtinc (vert_args(2)) invalid\n");
            return 0;
         }
         break;
      case 2:
      case 3:
         args = maxnl;
         break;
      default:
         args = 0;
         printf("Error: vertical invalid\n");
         return 0;
   }
   for (i=0;i<args;i++) {
      if (IS_MISSING(vert_args[i])) {
         printf("Error: vert_args(%d) invalid\n", i+1);
         return 0;
      }
   }

   return v5dCreate( filename, *numtimes, *numvars, *nr, *nc, nl,
                     (const char(*)[10]) names, timestamp, datestamp,
                     *compressmode,
                     *projection, proj_args, *vertical, vert_args );
}

/* **************************************************************
 * Set lowest levels for each variable (other than default of 0).
 * Input: lowlev - array [NumVars] of ints
 * Return:  1 = ok, 0 = error
 */
#ifdef UNDERSCORE
   int v5dsetlowlev_
#else
#  ifdef _CRAY
     int V5DSETLOWLEV
#  else
     int v5dsetlowlev
#  endif
#endif
          ( int *lowlev )
{
   return v5dSetLowLev(lowlev);
}

/* **************************************************************
 * Set the units for a variable.
 * Input: var - variable number in [1,NumVars]
 *        units - a character string
 * Return:  1 = ok, 0 = error
 */
#ifdef UNDERSCORE
   int v5dsetunits_
#else
#  ifdef _CRAY
     int V5DSETUNITS
#  else
     int v5dsetunits
#  endif
#endif
          ( int *var, char *name )
{
   return v5dSetUnits( *var, name );
}

/* **************************************************************
 * Write a grid of data to the file.
 * Input:  time - timestep in [1,NumTimes]
 *         var - timestep in [1,NumVars]
 *         data - array [nr*nc*nl] of floats
 * Return:  1 = ok, 0 = error
 */
#ifdef UNDERSCORE
   int v5dwrite_
#else
#  ifdef _CRAY
     int V5DWRITE
#  else
     int v5dwrite
#  endif
#endif
          ( const int *time, const int *var, const float *data )
{
   return v5dWrite( *time, *var, data );
}

/* **************************************************************
 * Close a simple v5d file.
 */
#ifdef UNDERSCORE
   int v5dclose_( void )
#else
#  ifdef _CRAY
     int V5DCLOSE( void )
#  else
     int v5dclose( void )
#  endif
#endif
{
   return v5dClose();
}

/*
 * Functions to do binary I/O of floats, ints.
 *
 * >>>> These functions are built on top of Unix I/O functions, not stdio! <<<<
 *
 * The file format is assumed to be BIG-ENDIAN.
 * If this code is compiled with -DLITTLE and executes on a little endian
 * CPU then byte-swapping will be done.
 *
 * If an ANSI compiler is used prototypes and ANSI function declarations
 * are used.  Otherwise use K&R conventions.
 *
 * If we're running on a CRAY (8-byte ints and floats), conversions will
 * be done as needed.
 */

/*
 * Updates:
 *
 * April 13, 1995, brianp
 *   added cray_to_ieee and iee_to_cray array conversion functions.
 *   fixed potential cray bug in write_float4_array function.
 *
 */

/**********************************************************************/
/******                     Byte Flipping                         *****/
/**********************************************************************/

#define FLIP4( n )  (  (n & 0xff000000) >> 24     \
                     | (n & 0x00ff0000) >> 8      \
                     | (n & 0x0000ff00) << 8      \
                     | (n & 0x000000ff) << 24  )

#define FLIP2( n )  (((unsigned short) (n & 0xff00)) >> 8  |  (n & 0x00ff) << 8)

/*
 * Flip the order of the 4 bytes in an array of 4-byte words.
 */
void flip4( const unsigned int *src, unsigned int *dest, int n )
{
   int i;

   for (i=0;i<n;i++) {
      unsigned int tmp = src[i];
      dest[i] = FLIP4( tmp );
   }
}

/*
 * Flip the order of the 2 bytes in an array of 2-byte words.
 */
void flip2( const unsigned short *src, unsigned short *dest, int n )
{
   int i;

   for (i=0;i<n;i++) {
      unsigned short tmp = src[i];
      dest[i] = FLIP2( tmp );
   }
}

#ifdef _CRAY
/*****************************************************************************
*
* The following source code is in the public domain.
* Specifically, we give to the public domain all rights for future licensing
* of the source code, all resale rights, and all publishing rights.
*
* We ask, but do not require, that the following message be included in all
* derived works:
*
* Portions developed at the National Center for Supercomputing Applications at
* the University of Illinois at Urbana-Champaign.
*
* THE UNIVERSITY OF ILLINOIS GIVES NO WARRANTY, EXPRESSED OR IMPLIED, FOR THE
* SOFTWARE AND/OR DOCUMENTATION PROVIDED, INCLUDING, WITHOUT LIMITATION,
* WARRANTY OF MERCHANTABILITY AND WARRANTY OF FITNESS FOR A PARTICULAR PURPOSE
*
****************************************************************************/

/** THESE ROUTINES MUST BE COMPILED ON THE CRAY ONLY SINCE THEY **/
/** REQUIRE 8-BYTES PER C-TYPE LONG                             **/

/* Cray to IEEE single precision */
static void c_to_if( long *t, const long *f)
{
    if (*f != 0){
        *t = (((*f & 0x8000000000000000) |      /* sign bit */
                 ((((*f & 0x7fff000000000000) >> 48)-16258) << 55)) + /* exp */
                 (((*f & 0x00007fffff000000) +
                    ((*f & 0x0000000000800000) << 1)) << 8));  /* mantissa */
    }
    else *t = *f;
}

#define C_TO_IF( T, F )                                                 \
        if (F != 0) {                                                   \
                T = (((F & 0x8000000000000000) |                        \
                ((((F & 0x7fff000000000000) >> 48)-16258) << 55)) +     \
                (((F & 0x00007fffff000000) +                            \
                ((F & 0x0000000000800000) << 1)) << 8));                \
        }                                                               \
        else {                                                          \
                T = F;                                                  \
        }

/* Cray to IEEE single precision */
static void c_to_if1( long *t, const long *f)
{
  /*
   * Clamp values to [-1e35, -1e-35] U {0} U [1e-35, 1e+35]
   * to prevent overflows and underflows that can occur when converting
   * from 8 byte to 4 byte floats.
   */

    const float* fp   = (const float*) f;
    float        x    = *fp;

    if ( x != 0 )
    {
#define SIGN(x) ((x) < 0 ? -1 : 1)
#define ABST(x)  ((x) < 0 ? -(x) : (x))

      const float  huge = 1e35;
      const float  tiny = 1e-35;

      printf("ABST( x )\n");
      if (      ABST( x ) < tiny ) x = tiny * SIGN( x );
      else if ( ABST( x ) > huge ) x = huge * SIGN( x );
      printf("huge\n");

      f = (const long*) &x;
      printf("f = \n");

#undef SIGN
#undef ABST
    }

    if (*f != 0){
    printf("8000000000000000 \n");
        *t = (((*f & 0x8000000000000000) |      /* sign bit */
                 ((((*f & 0x7fff000000000000) >> 48)-16258) << 55)) + /* exp */
                 (((*f & 0x00007fffff000000) +
                    ((*f & 0x0000000000800000) << 1)) << 8));  /* mantissa */
    printf("0x00000 \n");
    }
    else *t = *f;
}


#define C_TO_IF( T, F )							\
	if (F != 0) {							\
		T = (((F & 0x8000000000000000) |			\
		((((F & 0x7fff000000000000) >> 48)-16258) << 55)) +	\
		(((F & 0x00007fffff000000) +				\
		((F & 0x0000000000800000) << 1)) << 8));		\
	}								\
	else {								\
		T = F;							\
	}

/* IEEE single precison to Cray */
static void if_to_c( long *t, const long *f)
{
    if (*f != 0) {
        *t = (((*f & 0x8000000000000000) |
                ((*f & 0x7f80000000000000) >> 7) +
                (16258 << 48)) |
                (((*f & 0x007fffff00000000) >> 8) | (0x0000800000000000)));
        if ((*f << 1) == 0) *t = 0;
    }
    else *t = *f;
}

/* T and F must be longs! */
#define IF_TO_C( T, F )							\
	if (F != 0) {							\
		T = (((F & 0x8000000000000000) |			\
		((F & 0x7f80000000000000) >> 7) +			\
		(16258 << 48)) |					\
		(((F & 0x007fffff00000000) >> 8) | (0x0000800000000000)));  \
		if ((F << 1) == 0) T = 0;				\
	}								\
	else {								\
		T = F;							\
	}

/*
 * Convert an array of Cray 8-byte floats to an array of IEEE 4-byte floats.
 */
void cray_to_ieee_array( long *dest, const float *source, int n )
{
   long *dst;
   const long *src;
   long tmp1, tmp2;
   int i;

   dst = dest;
   src = (const long *) source;

   for (i=0;i<n;i+=2) {       /* add 1 in case n is odd */
      c_to_if( &tmp1, &src[i] );
      if ((i+1) >= n) {
         tmp2 = 0;
      }
      else {
         c_to_if( &tmp2, &src[i+1] );
      }
      *dst = (tmp1 & 0xffffffff00000000) | (tmp2 >> 32);
      dst++;
   }
}

/*
 * Convert an array of IEEE 4-byte floats to an array of 8-byte Cray floats.
 */
void ieee_to_cray_array( float *dest, const long *source, int n )
{
   long *dst;
   const long *src;
   int i;
   long ieee;

   src = source;
   dst = (long *) dest;

   for (i=0;i<n;i++) {
      /* most significant 4-bytes of ieee contain bit pattern to convert */
      if ((i&1)==0) {
         /* get upper half */
         ieee = src[i/2] & 0xffffffff00000000;
      }
      else {
         /* get lower half */
         ieee = src[i/2] << 32;
      }
      if_to_c( dst, &ieee );
      dst++;
   }
}

#endif /*_CRAY*/

/**********************************************************************/
/*****                         Write Functions                    *****/
/**********************************************************************/

/*
 * Write a block of bytes.
 * Input:  f - the file descriptor to write to.
 *         b - address of buffer to write.
 *         n - number of bytes to write.
 * Return:  number of bytes written, 0 if error.
 */
int write_bytes( int f, const void *b, int n )
{
   return write( f, b, n );
}

/*
 * Write an array of 2-byte integers.
 * Input:  f - file descriptor
 *         iarray - address to put integers
 *         n - number of integers to write.
 * Return:  number of integers written
 */
int write_int2_array( int f, const short *iarray, int n )
{
#ifdef _CRAY
   printf("write_int2_array not implemented!\n");
   exit(1);
#else
   int nwritten;
#ifdef LITTLE
   flip2( (const unsigned short *) iarray, (unsigned short *) iarray, n );
#endif
   nwritten = write( f, iarray, 2*n );
#ifdef LITTLE
   flip2( (const unsigned short *) iarray, (unsigned short *) iarray, n );
#endif
   if (nwritten<=0)
      return 0;
   return nwritten/2;
#endif
}

/*
 * Write an array of 2-byte unsigned integers.
 * Input:  f - file descriptor
 *         iarray - address to put integers
 *         n - number of integers to write.
 * Return:  number of integers written
 */
int write_uint2_array( int f, const unsigned short *iarray, int n )
{
#ifdef _CRAY
   int i, nwritten;
   unsigned char *buffer;
   buffer = (unsigned char *) malloc( 2*n );
   if (!buffer)  return 0;
   for (i=0;i<n;i++) {
      buffer[i*2] = (iarray[i] >> 8) & 0xff;
      buffer[i*2+1] = iarray[i] & 0xff;
   }
   nwritten = write( f, buffer, 2*n );
   free( buffer );
   if (nwritten<=0)
      return 0;
   else
      return nwritten/2;
#else
   int nwritten;
#ifdef LITTLE
   flip2( iarray, (unsigned short *) iarray, n );
#endif
   nwritten = write( f, iarray, 2*n );
#ifdef LITTLE
   flip2( iarray, (unsigned short *) iarray, n );
#endif
   if (nwritten<=0)
      return 0;
   else
      return nwritten/2;
#endif
}

/*
 * Write a 4-byte integer.
 *Input:  f - the file descriptor
 *         i - the integer
 * Return:  1 = ok, 0 = error
 */
int write_int4( int f, int i )
{
#ifdef _CRAY
   i = i << 32;
   return write( f, &i, 4 ) > 0;
#else
#  ifdef LITTLE
     i = FLIP4( i );
#  endif
   return write( f, &i, 4 ) > 0;
#endif
}

/*
 * Write an array of 4-byte integers.
 * Input:  f - the file descriptor
 *         i - the array of ints
 *           n - the number of ints in array
 *  Return:  number of integers written.
 */
int write_int4_array( int f, const int *i, int n )
{
#ifdef _CRAY
   int j, nwritten;
   char *buf, *b, *ptr;

   b = buf = (char *) malloc( n*4 + 8 );
   if (!b)
      return 0;
   ptr = (char *) i;
   for (j=0;j<n;j++) {
      ptr += 4;      /* skip upper 4 bytes */
      *b++ = *ptr++;
      *b++ = *ptr++;
      *b++ = *ptr++;
      *b++ = *ptr++;
   }
   nwritten = write( f, buf, 4*n );
   free( buf );
   if (nwritten<=0)
      return 0;
   else
      return nwritten / 4;
#else
#  ifdef LITTLE
      int nwritten;
      flip4( (const unsigned int *) i, (unsigned int *) i, n );
      nwritten = write( f, i, 4*n );
      flip4( (const unsigned int *) i, (unsigned int *) i, n );
      if (nwritten<=0)
         return 0;
      else
        return nwritten / 4;
#  else
      return write( f, i, 4*n ) / 4;
#  endif
#endif
}

/*
 * Write a 4-byte IEEE float.
 * Input:  f - the file descriptor
 *         x - the float
 * Return:  1 = ok, 0 = error
 */
int write_float4( int f, float x )
{
#ifdef _CRAY
   char buffer[8];
   c_to_if( (long *) buffer, (const long *) &x );
   return write( f, buffer, 4 ) > 0;
#else
#  ifdef LITTLE
      float y;
      unsigned int *iptr = (unsigned int *) &y, temp;
      y = (float) x;
      temp = FLIP4( *iptr );
      return write( f, &temp, 4 ) > 0;
#  else
      float y;
      y = (float) x;
      return write( f, &y, 4 ) > 0;
#  endif
#endif
}

/*
 * Write an array of 4-byte IEEE floating point numbers.
 * Input:  f - the file descriptor
 *         x - the array of floats
 *         n - number of floats in array
 * Return:  number of float written.
 */
int write_float4_array( int f, const float *x, int n )
{
#ifdef _CRAY
   /* convert cray floats to IEEE and put into buffer */
   int nwritten;
   long *buffer;
   buffer = (long *) malloc( n*4 + 8 );
   if (!buffer)
      return 0;
   cray_to_ieee_array( buffer, x, n );
   nwritten = write( f, buffer, 4*n );
   free( buffer );
   if (nwritten<=0)
      return 0;
   else
      return nwritten / 4;
#else
#  ifdef LITTLE
      int nwritten;
      flip4( (const unsigned int *) x, (unsigned int *) x, n );
      nwritten = write( f, x, 4*n );
      flip4( (const unsigned int *) x, (unsigned int *) x, n );
      if (nwritten<=0)
         return 0;
      else
         return nwritten / 4;
#  else
      return write( f, x, 4*n ) / 4;
#  endif
#endif
}

/*
 * Write a block of memory.
 * Input:  f - file descriptor
 *         data - address of first byte
 *         elements - number of elements to write
 *         elsize - size of each element to write (1, 2 or 4)
 * Return: number of elements written
 */
int write_block( int f, const void *data, int elements, int elsize )
{
   if (elsize==1) {
      return write( f, data, elements );
   }
   else if (elsize==2) {
#ifdef LITTLE
      int n;
      flip2( (const unsigned short *) data, (unsigned short *) data, elements);
      n = write( f, data, elements*2 ) / 2;
      flip2( (const unsigned short *) data, (unsigned short *) data, elements);
      return n;
#else
      return write( f, data, elements*2 ) / 2;
#endif
   }
   else if (elsize==4) {
#ifdef LITTLE
      int n;
      flip4( (const unsigned int *) data, (unsigned int *) data, elements );
      n = write( f, data, elements*4 ) / 4;
      flip4( (const unsigned int *) data, (unsigned int *) data, elements );
      return n;
#else
      return write( f, data, elements*4 ) / 4;
#endif
   }
   else {
      printf("Fatal error in write_block(): bad elsize (%d)\n", elsize );
      abort();
   }
   return 0;
}