v5d.c

/* v5d.c */

/* Vis5D version 5.1 */

/*
Vis5D system for visualizing five dimensional gridded data sets
Copyright (C) 1990 - 1997 Bill Hibbard, Johan Kellum, Brian Paul,
Dave Santek, and Andre Battaiola.

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 1, or (at your option)
any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/


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



/*
 * 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 <math.h>
#include "binio.h"
#include "v5d.h"
#include "vis5d.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 ave, var;
      float min, max;
      min = BIGVALUE;
      max = SMALLVALUE;
      ave = 0.0;
      var = 0.0;
      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;
/* WLH 5 Nov 98
         b = gb[lev] - 0.0001;
*/
         /* WLH 5 Nov 98 */
         b = gb[lev];
                                /* 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 {
/* MJK 1.19.99
               compdata1[p] = (V5Dubyte) (int) ((data[p]-b) * one_over_a);
*/
               compdata1[p] = (V5Dubyte) rint((data[p]-b) * one_over_a);
               if (compdata1[p] >= 255){
                  compdata1[p] = (V5Dubyte) (int) (255.0 - .0001);
               }
            }
         }
      }
   }

   else if (compressmode == 2) {
      int i, lev, p;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float one_over_a, b;
/* WLH 5 Nov 98
         b = gb[lev] - 0.0001;
*/
         /* WLH 5 Nov 98 */
         b = gb[lev];

         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 {
/* MJK 3.2.99
               compvalue = (V5Dushort) (int) ((data[p]-b) * one_over_a);
*/
               compvalue = (V5Dushort) rint((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) rint((data[p]-b) * one_over_a);

/*
               compdata2[p] = (V5Dushort) (int) ((data[p]-b) * one_over_a);
*/
/* MJK 3.24.99 I put this here so if the value is close
   to the missing value and get's rounded up it won't come out
   as missing data */
               if (compdata2[p] == 65535){
                  compdata2[p] = 65534;
               }
            }
         }
         /* 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
   }
}



/*
 * Decompress a 3-D grid from 1-byte integers to 4-byte floats.
 * Input:  nr, nc, nl - size of grid
 *         compdata - array of [nr*nr*nl*compressmode] bytes
 *         ga, gb - arrays of decompression factors
 *         compressmode - 1, 2 or 4 bytes per grid point
 *         data - address to put decompressed values
 * Output:  data - uncompressed floating point data values
 */
void v5dDecompressGrid( int nr, int nc, int nl, int compressmode,
                        void *compdata, float ga[], float gb[],
                        float data[] )
{
   int nrnc = nr * nc;
   int nrncnl = nr * nc * nl;
   V5Dubyte *compdata1 = (V5Dubyte *) compdata;
   V5Dushort *compdata2 = (V5Dushort *) compdata;

   if (compressmode == 1) {
      int p, i, lev;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float a = ga[lev];
         float b = gb[lev];

         /* WLH 2-2-95 */
         float d, aa;
         int id;
         if (a > 0.0000000001) {
           d = b / a;
           id = floor(d);
           d = d - id;
           aa = a * 0.000001;
         }
         else {
           id = 1;
         }
         if (-254 <= id && id <= 0 && d < aa) {
           for (i=0;i<nrnc;i++,p++) {
              if (compdata1[p]==255) {
                 data[p] = MISSING;
              }
              else {
                 data[p] = (float) (int) compdata1[p] * a + b;
                 if (fabs(data[p]) < aa) data[p] = aa;
              }
           }
         }
         else {
           for (i=0;i<nrnc;i++,p++) {
              if (compdata1[p]==255) {
                 data[p] = MISSING;
              }
              else {
                 data[p] = (float) (int) compdata1[p] * a + b;
              }
           }
         }
         /* end of WLH 2-2-95 */
      }
   }

   else if (compressmode == 2) {
      int p, i, lev;
      p = 0;
      for (lev=0;lev<nl;lev++) {
         float a = ga[lev];
         float b = gb[lev];
#ifdef _CRAY
         /* this is tricky because sizeof(V5Dushort)==8, not 2 */
         for (i=0;i<nrnc;i++,p++) {
            int compvalue;
            compvalue = (compdata1[p*2] << 8) | compdata1[p*2+1];
            if (compvalue==65535) {
               data[p] = MISSING;
            }
            else {
               data[p] = (float) compvalue * a + b;
            }
         }
#else
         /* sizeof(V5Dushort)==2! */
         for (i=0;i<nrnc;i++,p++) {
            if (compdata2[p]==65535) {
               data[p] = MISSING;
            }
            else {
               data[p] = (float) (int) compdata2[p] * a + b;
            }
         }
#endif
      }
   }

   else {
      /* compressmode==4 */
#ifdef _CRAY
      ieee_to_cray_array( data, compdata, nrncnl );
#else
      /* other machines: just copy 4-byte IEEE floats */
      assert( sizeof(float)==4 );
      memcpy( data, compdata, nrncnl*4 );
#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 )
{