cgribexlib.c

/* Automatically generated by m214003 at 2012-12-07, do not edit */

/* CGRIBEXLIB_VERSION="1.5.6" */

#ifdef _ARCH_PWR6
#pragma options nostrict
#endif

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <string.h>
#include <ctype.h>
#include <stdarg.h>
#include <sys/types.h>

//#undef _GET_X86_COUNTER
//#undef _GET_IBM_COUNTER
//#undef _ARCH_PWR6

#if   defined(__GNUC__) && (__GNUC__ >= 4)
#elif defined(__ICC)    && (__ICC >= 1100)
#else
#define DISABLE_SIMD
#endif

#ifdef DISABLE_SIMD
#ifndef ENABLE_AVX
#undef __AVX__
#endif
#ifndef ENABLE_SSE2
#undef __SSE2__
#endif
#endif

#ifdef _GET_IBM_COUNTER
#include <libhpc.h>
#endif

#ifdef __AVX__
#include <float.h>
#include <stdint.h>
#include <inttypes.h>
#include <immintrin.h>
#ifdef _GET_X86_COUNTER
#include <x86intrin.h>
#endif
#else
#ifdef __SSE2__
#include <float.h>
#include <stdint.h>
#include <inttypes.h>
#include <emmintrin.h>
#ifdef _GET_X86_COUNTER
#include <x86intrin.h>
#endif
#endif
#endif

#include "file.h"
#include "dmemory.h"
#include "dtypes.h"
#include "calendar.h"
#include "timebase.h"

#ifndef _GRIB_INT_H
#define _GRIB_INT_H

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>


#if ! defined   (_CGRIBEX_H)
#  include "cgribex.h"
#endif
#if ! defined   (_ERROR_H)
#  include "error.h"
#endif
#if ! defined   (_DTYPES_H)
#  include "dtypes.h"
#endif

#if ! defined   (FALSE)
#  define  FALSE  0
#endif

#if ! defined   (TRUE)
#  define  TRUE  1
#endif

#if ! defined   (UCHAR)
#  define  UCHAR  unsigned char
#endif


#if defined (CRAY) || defined (SX) || defined (__uxpch__)
#  define VECTORCODE
#endif

#if defined (VECTORCODE)
#if  defined  (INT32)
#  define  GRIBPACK     unsigned INT32
#  define  PACK_GRIB    packInt32
#  define  UNPACK_GRIB  unpackInt32
#else
#  define  GRIBPACK     unsigned INT64
#  define  PACK_GRIB    packInt64
#  define  UNPACK_GRIB  unpackInt64
#endif
#else
#  define  GRIBPACK     unsigned char
#endif

#define  U_BYTEORDER     static union {unsigned long l; unsigned char c[sizeof(long)];} u_byteorder = {1}
#define  IS_BIGENDIAN()  (u_byteorder.c[sizeof(long) - 1])

#if defined (__xlC__) /* performance problems on IBM */
#ifndef DBL_IS_NAN
#  define DBL_IS_NAN(x)     ((x) != (x))
#endif
#else
#ifndef DBL_IS_NAN
#if  defined  (HAVE_DECL_ISNAN)
#  define DBL_IS_NAN(x)     (isnan(x))
#elif  defined  (FP_NAN)
#  define DBL_IS_NAN(x)     (fpclassify(x) == FP_NAN)
#else
#  define DBL_IS_NAN(x)     ((x) != (x))
#endif
#endif
#endif

#ifndef DBL_IS_EQUAL
/*#define DBL_IS_EQUAL(x,y) (!(x < y || y < x)) */
#  define DBL_IS_EQUAL(x,y) (DBL_IS_NAN(x)||DBL_IS_NAN(y)?(DBL_IS_NAN(x)&&DBL_IS_NAN(y)?1:0):!(x < y || y < x))
#endif

#ifndef IS_EQUAL
#  define IS_NOT_EQUAL(x,y) (x < y || y < x)
#  define IS_EQUAL(x,y)     (!IS_NOT_EQUAL(x,y))
#endif


#define  JP23SET    0x7FFFFF  /* 2**23 - 1 (---> 8388607)  */

#define  POW_2_M24  0.000000059604644775390625  /*  pow(2.0, -24.0) */

double intpow2(int x);

int gribrec_len(int b1, int b2, int b3);
int correct_bdslen(int bdslen, long recsize, long gribpos);

/* CDI converter routines */

/* param format:  DDDCCCNNN */

void    cdiDecodeParam(int param, int *dis, int *cat, int *num);
int     cdiEncodeParam(int dis, int cat, int num);

/* date format:  YYYYMMDD */
/* time format:  hhmmss   */

void    cdiDecodeDate(int date, int *year, int *month, int *day);
int     cdiEncodeDate(int year, int month, int day);

void    cdiDecodeTime(int time, int *hour, int *minute, int *second);
int     cdiEncodeTime(int hour, int minute, int second);

/* CALENDAR types */

#define  CALENDAR_STANDARD        0  /* don't change this value (used also in cgribexlib)! */
#define  CALENDAR_PROLEPTIC       1
#define  CALENDAR_360DAYS         2
#define  CALENDAR_365DAYS         3
#define  CALENDAR_366DAYS         4
#define  CALENDAR_NONE            5

extern FILE *grprsm;

extern int  CGRIBEX_Debug;

void   gprintf(const char *caller, const char *fmt, ...);

void   grsdef(void);

void   prtbin(int kin, int knbit, int *kout, int *kerr);
void   confp3(double pval, int *kexp, int *kmant, int kbits, int kround);
double decfp2(int kexp, int kmant);
void   ref2ibm(double *pref, int kbits);

void   scaleComplex(double *fpdata, int pcStart, int pcScale, int trunc, int inv);
void   scatterComplex(double *fpdata, int pcStart, int trunc, int nsp);
void   gatherComplex(double *fpdata, int pcStart, int trunc, int nsp);

void   scm0(double *pdl, double *pdr, double *pfl, double *pfr, int klg);
int    rowina2(double *p, int ko, int ki, double *pw,
	       int kcode, double msval, int *kret);
int    rowina3(double *p, int ko, int ki, double *pw,
	       int kcode, double msval, int *kret, int omisng, int operio, int oveggy);
int    qu2reg2(double *pfield, int *kpoint, int klat, int klon,
	       double *ztemp, double msval, int *kret);
int    qu2reg3(double *pfield, int *kpoint, int klat, int klon,
	       double msval, int *kret, int omisng, int operio, int oveggy);

#if  defined  (INT32)
long   packInt32(unsigned INT32 *up, unsigned char *cp, long bc, long tc);
#endif
long   packInt64(unsigned INT64 *up, unsigned char *cp, long bc, long tc);
#if  defined  (INT32)
long   unpackInt32(unsigned char *cp, unsigned INT32 *up, long bc, long tc);
#endif
long   unpackInt64(unsigned char *cp, unsigned INT64 *up, long bc, long tc);

void  gribEncode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
		 double *fsec3, int *isec4, double *fsec4, int klenp, int *kgrib,
		 int kleng, int *kword, int efunc, int *kret);

void  gribDecode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
		 double *fsec3, int *isec4, double *fsec4, int klenp, int *kgrib,
		 int kleng, int *kword, int dfunc, int *kret);

#endif  /* _GRIB_INT_H */
#ifndef _GRIBDECODE_H
#define _GRIBDECODE_H

#define  UNDEFINED          9.999e20


#define  GET_INT3(a,b,c)    ((1-(int) ((unsigned) (a & 128) >> 6)) * (int) (((a & 127) << 16)+(b<<8)+c))
#define  GET_INT2(a,b)      ((1-(int) ((unsigned) (a & 128) >> 6)) * (int) (((a & 127) << 8) + b))
#define  GET_INT1(a)        ((1-(int) ((unsigned) (a & 128) >> 6)) * (int) (a&127))

/* this requires a 32-bit default integer machine */
#define  GET_UINT4(a,b,c,d) ((int) ((a << 24) + (b << 16) + (c << 8) + (d)))
#define  GET_UINT3(a,b,c)   ((int) ((a << 16) + (b << 8)  + (c)))
#define  GET_UINT2(a,b)     ((int) ((a << 8)  + (b)))
#define  GET_UINT1(a)       ((int)  (a))

#define  BUDG_START(s)      (s[0]=='B' && s[1]=='U' && s[2]=='D' && s[3]=='G')
#define  TIDE_START(s)      (s[0]=='T' && s[1]=='I' && s[2]=='D' && s[3]=='E')
#define  GRIB_START(s)      (s[0]=='G' && s[1]=='R' && s[2]=='I' && s[3]=='B')
#define  GRIB_FIN(s)        (s[0]=='7' && s[1]=='7' && s[2]=='7' && s[3]=='7')

/* GRIB1 Section 0: Indicator Section (IS) */

#define  GRIB1_SECLEN(s)     GET_INT3(s[ 4], s[ 5], s[ 6])
#define  GRIB_EDITION(s)     GET_UINT1(s[ 7])

/* GRIB1 Section 1: Product Definition Section (PDS) */

#define  PDS_Len             GET_UINT3(pds[ 0], pds[ 1], pds[ 2])
#define  PDS_CodeTable       GET_UINT1(pds[ 3])
#define  PDS_CenterID        GET_UINT1(pds[ 4])
#define  PDS_ModelID         GET_UINT1(pds[ 5])
#define  PDS_GridDefinition  GET_UINT1(pds[ 6])
#define  PDS_Sec2Or3Flag     GET_UINT1(pds[ 7])
#define  PDS_HAS_GDS         ((pds[7] & 128) != 0)
#define  PDS_HAS_BMS         ((pds[7] &  64) != 0)
#define  PDS_Parameter       GET_UINT1(pds[ 8])
#define  PDS_LevelType       GET_UINT1(pds[ 9])
#define  PDS_Level1          (pds[10])
#define  PDS_Level2	     (pds[11])
#define  PDS_Level	     GET_UINT2(pds[10], pds[11])
#define  PDS_Year            GET_INT1(pds[12])
#define  PDS_Month           GET_UINT1(pds[13])
#define  PDS_Day             GET_UINT1(pds[14])
#define  PDS_Hour            GET_UINT1(pds[15])
#define  PDS_Minute          GET_UINT1(pds[16])
#define  PDS_Date            (PDS_Year*10000+PDS_Month*100+PDS_Day)
#define  PDS_Time            (PDS_Hour*100+PDS_Minute)
#define  PDS_TimeUnit        GET_UINT1(pds[17])
#define  PDS_TimePeriod1     GET_UINT1(pds[18])
#define  PDS_TimePeriod2     GET_UINT1(pds[19])
#define  PDS_TimeRange       GET_UINT1(pds[20])
#define  PDS_AvgNum          GET_UINT2(pds[21], pds[22])
#define  PDS_AvgMiss         GET_UINT1(pds[23])
#define  PDS_Century         GET_UINT1(pds[24])
#define  PDS_Subcenter       GET_UINT1(pds[25])
#define  PDS_DecimalScale    GET_INT2(pds[26],pds[27])


/* GRIB1 Section 2: Grid Description Section (GDS) */

#define  GDS_Len             ((gds) == NULL ? 0 : GET_UINT3(gds[ 0], gds[ 1], gds[ 2]))
#define  GDS_NV              GET_UINT1(gds[ 3])
#define  GDS_PVPL            GET_UINT1(gds[ 4])
#define  GDS_PV	             ((gds[3] ==    0) ? -1 : (int) gds[4] - 1)
#define  GDS_PL	             ((gds[4] == 0xFF) ? -1 : (int) gds[3] * 4 + (int) gds[4] - 1)
#define  GDS_GridType        GET_UINT1(gds[ 5])


/* GRIB1 Triangular grid of DWD */
#define  GDS_GME_NI2         GET_UINT2(gds[ 6], gds[ 7])
#define  GDS_GME_NI3         GET_UINT2(gds[ 8], gds[ 9])
#define  GDS_GME_ND          GET_UINT3(gds[10], gds[11], gds[12])
#define  GDS_GME_NI          GET_UINT3(gds[13], gds[14], gds[15])
#define  GDS_GME_AFlag       GET_UINT1(gds[16])
#define  GDS_GME_LatPP       GET_INT3(gds[17], gds[18], gds[19])
#define  GDS_GME_LonPP       GET_INT3(gds[20], gds[21], gds[22])
#define  GDS_GME_LonMPL      GET_INT3(gds[23], gds[24], gds[25])
#define  GDS_GME_BFlag       GET_UINT1(gds[27])

/* GRIB1 Spectral */
#define  GDS_PentaJ          GET_UINT2(gds[ 6], gds[ 7])
#define  GDS_PentaK          GET_UINT2(gds[ 8], gds[ 9])
#define  GDS_PentaM          GET_UINT2(gds[10], gds[11])
#define  GDS_RepType         GET_UINT1(gds[12])
#define  GDS_RepMode         GET_UINT1(gds[13])

/* GRIB1 Regular grid */
#define  GDS_NumLon          GET_UINT2(gds[ 6], gds[ 7])
#define  GDS_NumLat          GET_UINT2(gds[ 8], gds[ 9])
#define  GDS_FirstLat        GET_INT3(gds[10], gds[11], gds[12])
#define  GDS_FirstLon        GET_INT3(gds[13], gds[14], gds[15])
#define  GDS_ResFlag         GET_UINT1(gds[16])
#define  GDS_LastLat         GET_INT3(gds[17], gds[18], gds[19])
#define  GDS_LastLon         GET_INT3(gds[20], gds[21], gds[22])
#define  GDS_LonIncr         GET_UINT2(gds[23], gds[24])
#define  GDS_LatIncr         GET_UINT2(gds[25], gds[26])
#define  GDS_NumPar          GET_UINT2(gds[25], gds[26])
#define  GDS_ScanFlag        GET_UINT1(gds[27])
#define  GDS_LatSP           GET_INT3(gds[32], gds[33], gds[34])
#define  GDS_LonSP           GET_INT3(gds[35], gds[36], gds[37])
#define  GDS_RotAngle        GET_Real(&(gds[38]))

/* GRIB1 Lambert */
#define  GDS_Lambert_Lov     GET_INT3(gds[17], gds[18], gds[19])
#define  GDS_Lambert_dx	     GET_INT3(gds[20], gds[21], gds[22])
#define  GDS_Lambert_dy	     GET_INT3(gds[23], gds[24], gds[25])
#define  GDS_Lambert_ProjFlag GET_UINT1(gds[26])
#define  GDS_Lambert_LatS1   GET_INT3(gds[28], gds[29], gds[30])
#define  GDS_Lambert_LatS2   GET_INT3(gds[31], gds[32], gds[33])
#define  GDS_Lambert_LatSP   GET_INT3(gds[34], gds[35], gds[36])
#define  GDS_Lambert_LonSP   GET_INT3(gds[37], gds[37], gds[37])

/* GRIB1 Section 3: Bit Map Section (BMS) */

#define  BMS_Len	     ((bms) == NULL ? 0 : (int) (bms[0]<<16)+(bms[1]<<8)+bms[2])
#define  BMS_UnusedBits      (bms[3])
#define  BMS_Numeric         
#define  BMS_Bitmap	     ((bms) == NULL ? NULL : (bms)+6)
#define  BMS_BitmapSize      (((((bms[0]<<16)+(bms[1]<<8)+bms[2]) - 6)<<3) - bms[3])

/* GRIB1 Section 4: Binary Data Section (BDS) */

#define  BDS_Len	    ((int) ((bds[0]<<16)+(bds[1]<<8)+bds[2]))
#define  BDS_Flag	    (bds[3])
#define  BDS_BinScale       GET_INT2(bds[ 4], bds[ 5])
#define  BDS_RefValue       decfp2((int)bds[ 6], GET_UINT3(bds[ 7], bds[ 8], bds[ 9]))
#define  BDS_NumBits        ((int) bds[10])
#define  BDS_RealCoef       decfp2((int)bds[zoff+11], GET_UINT3(bds[zoff+12], bds[zoff+13], bds[zoff+14]))
#define  BDS_PackData       ((int) ((bds[zoff+11]<<8) + bds[zoff+12]))
#define  BDS_Power          GET_INT2(bds[zoff+13], bds[zoff+14])
#define  BDS_Z              (bds[13])

/* GRIB1 Section 5: End Section (ES) */

/* GRIB2 */

#define  GRIB2_SECLEN(section)   (GET_UINT4(section[0], section[1], section[2], section[3]))
#define  GRIB2_SECNUM(section)   (GET_UINT1(section[4]))

#endif  /* _GRIBDECODE_H */
#ifndef _GRIBENCODE_H
#define _GRIBENCODE_H


#define PutnZero(n) \
{ \
  int i; \
  for ( i = z; i < z+n; i++ ) lGrib[i] = 0; \
  z += n; \
}

#define Put1Byte(Value)  (lGrib[z++] = (Value))
#define Put2Byte(Value) ((lGrib[z++] = (Value) >>  8), \
                         (lGrib[z++] = (Value)))
#define Put3Byte(Value) ((lGrib[z++] = (Value) >> 16), \
                         (lGrib[z++] = (Value) >>  8), \
                         (lGrib[z++] = (Value)))
#define Put4Byte(Value) ((lGrib[z++] = (Value) >> 24), \
                         (lGrib[z++] = (Value) >> 16), \
                         (lGrib[z++] = (Value) >>  8), \
                         (lGrib[z++] = (Value)))

#define Put1Int(Value)  {ival = Value; if ( ival < 0 ) ival =     0x80 - ival; Put1Byte(ival);}
#define Put2Int(Value)  {ival = Value; if ( ival < 0 ) ival =   0x8000 - ival; Put2Byte(ival);}
#define Put3Int(Value)  {ival = Value; if ( ival < 0 ) ival = 0x800000 - ival; Put3Byte(ival);}

#define Put1Real(Value)          \
{                                \
  confp3(Value, &exponent, &mantissa, BitsPerInt, 1); \
  Put1Byte(exponent);            \
  Put3Byte(mantissa);            \
}

#endif  /* _GRIBENCODE_H */


double _pow2tab[158] = {
 /* pow(2.0,  0.0) */  1.0,
 /* pow(2.0,  1.0) */  2.0,
 /* pow(2.0,  2.0) */  4.0,
 /* pow(2.0,  3.0) */  8.0,
 /* pow(2.0,  4.0) */  16.0,
 /* pow(2.0,  5.0) */  32.0,
 /* pow(2.0,  6.0) */  64.0,
 /* pow(2.0,  7.0) */  128.0,
 /* pow(2.0,  8.0) */  256.0,
 /* pow(2.0,  9.0) */  512.0,
 /* pow(2.0, 10.0) */  1024.0,
 /* pow(2.0, 11.0) */  2048.0,
 /* pow(2.0, 12.0) */  4096.0,
 /* pow(2.0, 13.0) */  8192.0,
 /* pow(2.0, 14.0) */  16384.0,
 /* pow(2.0, 15.0) */  32768.0,
 /* pow(2.0, 16.0) */  65536.0,
 /* pow(2.0, 17.0) */  131072.0,
 /* pow(2.0, 18.0) */  262144.0,
 /* pow(2.0, 19.0) */  524288.0,
 /* pow(2.0, 20.0) */  1048576.0,
 /* pow(2.0, 21.0) */  2097152.0,
 /* pow(2.0, 22.0) */  4194304.0,
 /* pow(2.0, 23.0) */  8388608.0,
 /* pow(2.0, 24.0) */  16777216.0,
 /* pow(2.0, 25.0) */  33554432.0,
 /* pow(2.0, 26.0) */  67108864.0,
 /* pow(2.0, 27.0) */  134217728.0,
 /* pow(2.0, 28.0) */  268435456.0,
 /* pow(2.0, 29.0) */  536870912.0,
 /* pow(2.0, 30.0) */  1073741824.0,
 /* pow(2.0, 31.0) */  2147483648.0,
 /* pow(2.0, 32.0) */  4294967296.0,
 /* pow(2.0, 33.0) */  8589934592.0,
 /* pow(2.0, 34.0) */  17179869184.0,
 /* pow(2.0, 35.0) */  34359738368.0,
 /* pow(2.0, 36.0) */  68719476736.0,
 /* pow(2.0, 37.0) */  137438953472.0,
 /* pow(2.0, 38.0) */  274877906944.0,
 /* pow(2.0, 39.0) */  549755813888.0,
 /* pow(2.0, 40.0) */  1099511627776.0,
 /* pow(2.0, 41.0) */  2199023255552.0,
 /* pow(2.0, 42.0) */  4398046511104.0,
 /* pow(2.0, 43.0) */  8796093022208.0,
 /* pow(2.0, 44.0) */  17592186044416.0,
 /* pow(2.0, 45.0) */  35184372088832.0,
 /* pow(2.0, 46.0) */  70368744177664.0,
 /* pow(2.0, 47.0) */  140737488355328.0,
 /* pow(2.0, 48.0) */  281474976710656.0,
 /* pow(2.0, 49.0) */  562949953421312.0,
 /* pow(2.0, 50.0) */  1125899906842624.0,
 /* pow(2.0, 51.0) */  2251799813685248.0,
 /* pow(2.0, 52.0) */  4503599627370496.0,
 /* pow(2.0, 53.0) */  9007199254740992.0,
 /* pow(2.0, 54.0) */  18014398509481984.0,
 /* pow(2.0, 55.0) */  36028797018963968.0,
 /* pow(2.0, 56.0) */  72057594037927936.0,
 /* pow(2.0, 57.0) */  144115188075855872.0,
 /* pow(2.0, 58.0) */  288230376151711744.0,
 /* pow(2.0, 59.0) */  576460752303423488.0,
 /* pow(2.0, 60.0) */  1152921504606846976.0,
 /* pow(2.0, 61.0) */  2305843009213693952.0,
 /* pow(2.0, 62.0) */  4611686018427387904.0,
 /* pow(2.0, 63.0) */  9223372036854775808.0,
 /* pow(2.0, 64.0) */  18446744073709551616.0,
 /* pow(2.0, 65.0) */  36893488147419103232.0,
 /* pow(2.0, 66.0) */  73786976294838206464.0,
 /* pow(2.0, 67.0) */  147573952589676412928.0,
 /* pow(2.0, 68.0) */  295147905179352825856.0,
 /* pow(2.0, 69.0) */  590295810358705651712.0,
 /* pow(2.0, 70.0) */  1180591620717411303424.0,
 /* pow(2.0, 71.0) */  2361183241434822606848.0,
 /* pow(2.0, 72.0) */  4722366482869645213696.0,
 /* pow(2.0, 73.0) */  9444732965739290427392.0,
 /* pow(2.0, 74.0) */  18889465931478580854784.0,
 /* pow(2.0, 75.0) */  37778931862957161709568.0,
 /* pow(2.0, 76.0) */  75557863725914323419136.0,
 /* pow(2.0, 77.0) */  151115727451828646838272.0,
 /* pow(2.0, 78.0) */  302231454903657293676544.0,
 /* pow(2.0, 79.0) */  604462909807314587353088.0,
 /* pow(2.0, 80.0) */  1208925819614629174706176.0,
 /* pow(2.0, 81.0) */  2417851639229258349412352.0,
 /* pow(2.0, 82.0) */  4835703278458516698824704.0,
 /* pow(2.0, 83.0) */  9671406556917033397649408.0,
 /* pow(2.0, 84.0) */  19342813113834066795298816.0,
 /* pow(2.0, 85.0) */  38685626227668133590597632.0,
 /* pow(2.0, 86.0) */  77371252455336267181195264.0,
 /* pow(2.0, 87.0) */  154742504910672534362390528.0,
 /* pow(2.0, 88.0) */  309485009821345068724781056.0,
 /* pow(2.0, 89.0) */  618970019642690137449562112.0,
 /* pow(2.0, 90.0) */  1237940039285380274899124224.0,
 /* pow(2.0, 91.0) */  2475880078570760549798248448.0,
 /* pow(2.0, 92.0) */  4951760157141521099596496896.0,
 /* pow(2.0, 93.0) */  9903520314283042199192993792.0,
 /* pow(2.0, 94.0) */  19807040628566084398385987584.0,
 /* pow(2.0, 95.0) */  39614081257132168796771975168.0,
 /* pow(2.0, 96.0) */  79228162514264337593543950336.0,
 /* pow(2.0, 97.0) */  158456325028528675187087900672.0,
 /* pow(2.0, 98.0) */  316912650057057350374175801344.0,
 /* pow(2.0, 99.0) */  633825300114114700748351602688.0,
 /* pow(2.0, 100.0) */  1267650600228229401496703205376.0,
 /* pow(2.0, 101.0) */  2535301200456458802993406410752.0,
 /* pow(2.0, 102.0) */  5070602400912917605986812821504.0,
 /* pow(2.0, 103.0) */  10141204801825835211973625643008.0,
 /* pow(2.0, 104.0) */  20282409603651670423947251286016.0,
 /* pow(2.0, 105.0) */  40564819207303340847894502572032.0,
 /* pow(2.0, 106.0) */  81129638414606681695789005144064.0,
 /* pow(2.0, 107.0) */  162259276829213363391578010288128.0,
 /* pow(2.0, 108.0) */  324518553658426726783156020576256.0,
 /* pow(2.0, 109.0) */  649037107316853453566312041152512.0,
 /* pow(2.0, 110.0) */  1298074214633706907132624082305024.0,
 /* pow(2.0, 111.0) */  2596148429267413814265248164610048.0,
 /* pow(2.0, 112.0) */  5192296858534827628530496329220096.0,
 /* pow(2.0, 113.0) */  10384593717069655257060992658440192.0,
 /* pow(2.0, 114.0) */  20769187434139310514121985316880384.0,
 /* pow(2.0, 115.0) */  41538374868278621028243970633760768.0,
 /* pow(2.0, 116.0) */  83076749736557242056487941267521536.0,
 /* pow(2.0, 117.0) */  166153499473114484112975882535043072.0,
 /* pow(2.0, 118.0) */  332306998946228968225951765070086144.0,
 /* pow(2.0, 119.0) */  664613997892457936451903530140172288.0,
 /* pow(2.0, 120.0) */  1329227995784915872903807060280344576.0,
 /* pow(2.0, 121.0) */  2658455991569831745807614120560689152.0,
 /* pow(2.0, 122.0) */  5316911983139663491615228241121378304.0,
 /* pow(2.0, 123.0) */  10633823966279326983230456482242756608.0,
 /* pow(2.0, 124.0) */  21267647932558653966460912964485513216.0,
 /* pow(2.0, 125.0) */  42535295865117307932921825928971026432.0,
 /* pow(2.0, 126.0) */  85070591730234615865843651857942052864.0,
 /* pow(2.0, 127.0) */  170141183460469231731687303715884105728.0,
 /* pow(2.0, 128.0) */  340282366920938463463374607431768211456.0,
 /* pow(2.0, 129.0) */  680564733841876926926749214863536422912.0,
 /* pow(2.0, 130.0) */  1361129467683753853853498429727072845824.0,
 /* pow(2.0, 131.0) */  2722258935367507707706996859454145691648.0,
 /* pow(2.0, 132.0) */  5444517870735015415413993718908291383296.0,
 /* pow(2.0, 133.0) */  10889035741470030830827987437816582766592.0,
 /* pow(2.0, 134.0) */  21778071482940061661655974875633165533184.0,
 /* pow(2.0, 135.0) */  43556142965880123323311949751266331066368.0,
 /* pow(2.0, 136.0) */  87112285931760246646623899502532662132736.0,
 /* pow(2.0, 137.0) */  174224571863520493293247799005065324265472.0,
 /* pow(2.0, 138.0) */  348449143727040986586495598010130648530944.0,
 /* pow(2.0, 139.0) */  696898287454081973172991196020261297061888.0,
 /* pow(2.0, 140.0) */  1393796574908163946345982392040522594123776.0,
 /* pow(2.0, 141.0) */  2787593149816327892691964784081045188247552.0,
 /* pow(2.0, 142.0) */  5575186299632655785383929568162090376495104.0,
 /* pow(2.0, 143.0) */  11150372599265311570767859136324180752990208.0,
 /* pow(2.0, 144.0) */  22300745198530623141535718272648361505980416.0,
 /* pow(2.0, 145.0) */  44601490397061246283071436545296723011960832.0,
 /* pow(2.0, 146.0) */  89202980794122492566142873090593446023921664.0,
 /* pow(2.0, 147.0) */  178405961588244985132285746181186892047843328.0,
 /* pow(2.0, 148.0) */  356811923176489970264571492362373784095686656.0,
 /* pow(2.0, 149.0) */  713623846352979940529142984724747568191373312.0,
 /* pow(2.0, 150.0) */  1427247692705959881058285969449495136382746624.0,
 /* pow(2.0, 151.0) */  2854495385411919762116571938898990272765493248.0,
 /* pow(2.0, 152.0) */  5708990770823839524233143877797980545530986496.0,
 /* pow(2.0, 153.0) */  11417981541647679048466287755595961091061972992.0,
 /* pow(2.0, 154.0) */  22835963083295358096932575511191922182123945984.0,
 /* pow(2.0, 155.0) */  45671926166590716193865151022383844364247891968.0,
 /* pow(2.0, 156.0) */  91343852333181432387730302044767688728495783936.0,
 /* pow(2.0, 157.0) */  182687704666362864775460604089535377456991567872.0,
};


double _pow16tab[71] = {
 /* pow(16.0,  0.0) */  1.0,
 /* pow(16.0,  1.0) */  16.0,
 /* pow(16.0,  2.0) */  256.0,
 /* pow(16.0,  3.0) */  4096.0,
 /* pow(16.0,  4.0) */  65536.0,
 /* pow(16.0,  5.0) */  1048576.0,
 /* pow(16.0,  6.0) */  16777216.0,
 /* pow(16.0,  7.0) */  268435456.0,
 /* pow(16.0,  8.0) */  4294967296.0,
 /* pow(16.0,  9.0) */  68719476736.0,
 /* pow(16.0, 10.0) */  1099511627776.0,
 /* pow(16.0, 11.0) */  17592186044416.0,
 /* pow(16.0, 12.0) */  281474976710656.0,
 /* pow(16.0, 13.0) */  4503599627370496.0,
 /* pow(16.0, 14.0) */  72057594037927936.0,
 /* pow(16.0, 15.0) */  1152921504606846976.0,
 /* pow(16.0, 16.0) */  18446744073709551616.0,
 /* pow(16.0, 17.0) */  295147905179352825856.0,
 /* pow(16.0, 18.0) */  4722366482869645213696.0,
 /* pow(16.0, 19.0) */  75557863725914323419136.0,
 /* pow(16.0, 20.0) */  1208925819614629174706176.0,
 /* pow(16.0, 21.0) */  19342813113834066795298816.0,
 /* pow(16.0, 22.0) */  309485009821345068724781056.0,
 /* pow(16.0, 23.0) */  4951760157141521099596496896.0,
 /* pow(16.0, 24.0) */  79228162514264337593543950336.0,
 /* pow(16.0, 25.0) */  1267650600228229401496703205376.0,
 /* pow(16.0, 26.0) */  20282409603651670423947251286016.0,
 /* pow(16.0, 27.0) */  324518553658426726783156020576256.0,
 /* pow(16.0, 28.0) */  5192296858534827628530496329220096.0,
 /* pow(16.0, 29.0) */  83076749736557242056487941267521536.0,
 /* pow(16.0, 30.0) */  1329227995784915872903807060280344576.0,
 /* pow(16.0, 31.0) */  21267647932558653966460912964485513216.0,
 /* pow(16.0, 32.0) */  340282366920938463463374607431768211456.0,
 /* pow(16.0, 33.0) */  5444517870735015415413993718908291383296.0,
 /* pow(16.0, 34.0) */  87112285931760246646623899502532662132736.0,
 /* pow(16.0, 35.0) */  1393796574908163946345982392040522594123776.0,
 /* pow(16.0, 36.0) */  22300745198530623141535718272648361505980416.0,
 /* pow(16.0, 37.0) */  356811923176489970264571492362373784095686656.0,
 /* pow(16.0, 38.0) */  5708990770823839524233143877797980545530986496.0,
 /* pow(16.0, 39.0) */  91343852333181432387730302044767688728495783936.0,
 /* pow(16.0, 40.0) */  1461501637330902918203684832716283019655932542976.0,
 /* pow(16.0, 41.0) */  23384026197294446691258957323460528314494920687616.0,
 /* pow(16.0, 42.0) */  374144419156711147060143317175368453031918731001856.0,
 /* pow(16.0, 43.0) */  5986310706507378352962293074805895248510699696029696.0,
 /* pow(16.0, 44.0) */  95780971304118053647396689196894323976171195136475136.0,
 /* pow(16.0, 45.0) */  1532495540865888858358347027150309183618739122183602176.0,
 /* pow(16.0, 46.0) */  24519928653854221733733552434404946937899825954937634816.0,
 /* pow(16.0, 47.0) */  392318858461667547739736838950479151006397215279002157056.0,
 /* pow(16.0, 48.0) */  6277101735386680763835789423207666416102355444464034512896.0,
 /* pow(16.0, 49.0) */  100433627766186892221372630771322662657637687111424552206336.0,
 /* pow(16.0, 50.0) */  1606938044258990275541962092341162602522202993782792835301376.0,
 /* pow(16.0, 51.0) */  25711008708143844408671393477458601640355247900524685364822016.0,
 /* pow(16.0, 52.0) */  411376139330301510538742295639337626245683966408394965837152256.0,
 /* pow(16.0, 53.0) */  6582018229284824168619876730229402019930943462534319453394436096.0,
 /* pow(16.0, 54.0) */  105312291668557186697918027683670432318895095400549111254310977536.0,
 /* pow(16.0, 55.0) */  1684996666696914987166688442938726917102321526408785780068975640576.0,
 /* pow(16.0, 56.0) */  26959946667150639794667015087019630673637144422540572481103610249216.0,
 /* pow(16.0, 57.0) */  431359146674410236714672241392314090778194310760649159697657763987456.0,
 /* pow(16.0, 58.0) */  6901746346790563787434755862277025452451108972170386555162524223799296.0,
 /* pow(16.0, 59.0) */  110427941548649020598956093796432407239217743554726184882600387580788736.0,
 /* pow(16.0, 60.0) */  1766847064778384329583297500742918515827483896875618958121606201292619776.0,
 /* pow(16.0, 61.0) */  28269553036454149273332760011886696253239742350009903329945699220681916416.0,
 /* pow(16.0, 62.0) */  452312848583266388373324160190187140051835877600158453279131187530910662656.0,
 /* pow(16.0, 63.0) */  7237005577332262213973186563042994240829374041602535252466099000494570602496.0,
 /* pow(16.0, 64.0) */  115792089237316195423570985008687907853269984665640564039457584007913129639936.0,
 /* pow(16.0, 65.0) */  1852673427797059126777135760139006525652319754650249024631321344126610074238976.0,
 /* pow(16.0, 66.0) */  29642774844752946028434172162224104410437116074403984394101141506025761187823616.0,
 /* pow(16.0, 67.0) */  474284397516047136454946754595585670566993857190463750305618264096412179005177856.0,
 /* pow(16.0, 68.0) */  7588550360256754183279148073529370729071901715047420004889892225542594864082845696.0,
 /* pow(16.0, 69.0) */  121416805764108066932466369176469931665150427440758720078238275608681517825325531136.0,
 /* pow(16.0, 70.0) */  1942668892225729070919461906823518906642406839052139521251812409738904285205208498176.0,
};

static int _pow2tab_size = sizeof(_pow2tab)/sizeof(double);

void gen_pow2tab(void)
{
  int jloop;

  for ( jloop = 0; jloop < 158; jloop++ )
    printf(" /* pow(2.0, %2d.0) */  %.1f,\n", jloop,  pow(2.0, (double) jloop));
}


void gen_pow16tab(void)
{
  double pval;
  int iexp;

  for ( iexp = 0; iexp < 71; iexp++ )
    {
      pval = pow(16.0, (double)(iexp));
      printf(" /* pow(16.0, %2d.0) */  %.1f,\n", iexp, pval);
    }
}


double intpow2(int x)
{
  if ( x < _pow2tab_size )
    return (_pow2tab[x]);
  else
    return (pow(2.0, (double) x));
}
//#undef _GET_X86_COUNTER
//#undef _GET_IBM_COUNTER
//#undef _GET_MACH_COUNTER
//#undef _ARCH_PWR6

#if   defined(__GNUC__) && (__GNUC__ >= 4)
#elif defined(__ICC)    && (__ICC >= 1100)
#elif defined(__clang__)
#else
#define DISABLE_SIMD
#endif

#if defined _GET_IBM_COUNTER
#elif defined _GET_X86_COUNTER
#elif defined _GET_MACH_COUNTER
#endif

//#define DISABLE_SIMD

#ifdef DISABLE_SIMD
#ifndef ENABLE_AVX
#undef __AVX__
#endif
#ifndef ENABLE_SSE2
#undef __SSE2__
#endif
#endif

#if defined __AVX__
#elif defined __SSE2__
#endif

#if defined __AVX__

static
void avx_minmax_val(const double *restrict buf, size_t nframes, double *min, double *max)
{
  double fmin[4], fmax[4];
  __m256d current_max, current_min, work;

  // load max and min values into all four slots of the YMM registers
  current_min = _mm256_set1_pd(*min);
  current_max = _mm256_set1_pd(*max);

  // Work input until "buf" reaches 32 byte alignment
  while ( ((unsigned long)buf) % 32 != 0 && nframes > 0) {

    // Load the next double into the work buffer
    work = _mm256_set1_pd(*buf);
    current_min = _mm256_min_pd(current_min, work);
    current_max = _mm256_max_pd(current_max, work);
    buf++;
    nframes--;
  }

  while (nframes >= 16) {

    (void) _mm_prefetch(buf+8, _MM_HINT_NTA);

    work = _mm256_load_pd(buf);
    current_min = _mm256_min_pd(current_min, work);
    current_max = _mm256_max_pd(current_max, work);
    buf += 4;

    work = _mm256_load_pd(buf);
    current_min = _mm256_min_pd(current_min, work);
    current_max = _mm256_max_pd(current_max, work);
    buf += 4;

    (void) _mm_prefetch(buf+8, _MM_HINT_NTA);

    work = _mm256_load_pd(buf);
    current_min = _mm256_min_pd(current_min, work);
    current_max = _mm256_max_pd(current_max, work);
    buf += 4;

    work = _mm256_load_pd(buf);
    current_min = _mm256_min_pd(current_min, work);
    current_max = _mm256_max_pd(current_max, work);
    buf += 4;
    nframes -= 16;
  }

  // work through aligned buffers
  while (nframes >= 4) {
    work = _mm256_load_pd(buf);
    current_min = _mm256_min_pd(current_min, work);
    current_max = _mm256_max_pd(current_max, work);
    buf += 4;
    nframes -= 4;
  }

  // work through the remainung values
  while ( nframes > 0) {
    work = _mm256_set1_pd(*buf);
    current_min = _mm256_min_pd(current_min, work);
    current_max = _mm256_max_pd(current_max, work);
    buf++;
    nframes--;
  }

  // find min & max value through shuffle tricks

  work = current_min;
  work = _mm256_shuffle_pd(work, work, 5);
  work = _mm256_min_pd (work, current_min);
  current_min = work;
  work = _mm256_permute2f128_pd(work, work, 1);
  work = _mm256_min_pd (work, current_min);
  _mm256_storeu_pd(fmin, work);

  work = current_max;
  work = current_max;
  work = _mm256_shuffle_pd(work, work, 5);
  work = _mm256_max_pd (work, current_max);
  current_max = work;
  work = _mm256_permute2f128_pd(work, work, 1);
  work = _mm256_max_pd (work, current_max);
  _mm256_storeu_pd(fmax, work);

  *min = fmin[0];
  *max = fmax[0];

  return;
}

#elif defined __SSE2__

static
void sse2_minmax_val(const double *restrict buf, size_t nframes, double *min, double *max)
{
  __m128d current_max, current_min, work;
  
  // load starting max and min values into all slots of the XMM registers
  current_min = _mm_set1_pd(*min);
  current_max = _mm_set1_pd(*max);
  
  // work on input until buf reaches 16 byte alignment
  while ( ((unsigned long)buf) % 16 != 0 && nframes > 0) {
    
    // load one double and replicate
    work = _mm_set1_pd(*buf);    
    current_min = _mm_min_pd(current_min, work);
    current_max = _mm_max_pd(current_max, work);    
    buf++;
    nframes--;
  }
  
  while (nframes >= 8) {
    // use 64 byte prefetch for double octetts
    // __builtin_prefetch(buf+64,0,0); // for GCC 4.3.2 +

    work = _mm_load_pd(buf);
    current_min = _mm_min_pd(current_min, work);
    current_max = _mm_max_pd(current_max, work);
    buf += 2;
    work = _mm_load_pd(buf);
    current_min = _mm_min_pd(current_min, work);
    current_max = _mm_max_pd(current_max, work);
    buf += 2;
    work = _mm_load_pd(buf);
    current_min = _mm_min_pd(current_min, work);
    current_max = _mm_max_pd(current_max, work);
    buf += 2;
    work = _mm_load_pd(buf);
    current_min = _mm_min_pd(current_min, work);
    current_max = _mm_max_pd(current_max, work);
    buf += 2;
    nframes -= 8;
  }

  // work through smaller chunks of aligned buffers without prefetching
  while (nframes >= 2) {
    work = _mm_load_pd(buf);
    current_min = _mm_min_pd(current_min, work);
    current_max = _mm_max_pd(current_max, work);
    buf += 2;
    nframes -= 2;
  }

  // work through the remaining value
  while ( nframes > 0) {
    // load the last double and replicate
    work = _mm_set1_pd(*buf);
    current_min = _mm_min_pd(current_min, work);
    current_max = _mm_max_pd(current_max, work);
    buf++;
    nframes--;
  }

  // find final min and max value through shuffle tricks
  work = current_min;
  work = _mm_shuffle_pd(work, work, _MM_SHUFFLE2(0, 1));
  work = _mm_min_pd (work, current_min);
  _mm_store_sd(min, work);
  work = current_max;
  work = _mm_shuffle_pd(work, work, _MM_SHUFFLE2(0, 1));
  work = _mm_max_pd (work, current_max);
  _mm_store_sd(max, work);

  return;
}

#endif

static
void minmax_val(const double *restrict data, size_t datasize, double *fmin, double *fmax)
{
#if defined _GET_X86_COUNTER || defined _GET_MACH_COUNTER 
  uint64_t start_minmax, end_minmax;
#endif

  if ( datasize < 1 ) return;

#ifdef _GET_X86_COUNTER 
  start_minmax = _rdtsc();
#endif
#ifdef _GET_MACH_COUNTER 
  start_minmax = mach_absolute_time();
#endif

#if defined __AVX__

  avx_minmax_val(data, datasize, fmin, fmax);

#elif defined __SSE2__

  sse2_minmax_val(data, datasize, fmin, fmax);

#else

#ifdef _ARCH_PWR6
#define __UNROLL_DEPTH_1 6

  // to allow pipelining we have to unroll 

#ifdef _GET_IBM_COUNTER 
  hpmStart(1, "minmax fsel");
#endif
  {
    size_t i, j;
    size_t residual =  datasize % __UNROLL_DEPTH_1;
    size_t ofs = datasize - residual;
    double register dmin[__UNROLL_DEPTH_1];
    double register dmax[__UNROLL_DEPTH_1];

    for ( j = 0; j < __UNROLL_DEPTH_1; j++) 
      {
	dmin[j] = data[0];
	dmax[j] = data[0];
      }
    
    for ( i = 0; i < datasize - residual; i += __UNROLL_DEPTH_1 ) 
      {
	for (j = 0; j < __UNROLL_DEPTH_1; j++) 
	  {
	    dmin[j] = __fsel(dmin[j] - data[i+j], data[i+j], dmin[j]);
	    dmax[j] = __fsel(data[i+j] - dmax[j], data[i+j], dmax[j]);
	  }
      }

    for (j = 0; j < residual; j++) 
      {
	dmin[j] = __fsel(dmin[j] - data[ofs+j], data[ofs+j], dmin[j]);
	dmax[j] = __fsel(data[ofs+j] - dmax[j], data[ofs+j], dmax[j]);
      }

    for ( j = 0; j < __UNROLL_DEPTH_1; j++) 
      {
	*fmin = __fsel(*fmin - dmin[j], dmin[j], *fmin);
	*fmax = __fsel(dmax[j] - *fmax, dmax[j], *fmax);
      }
  }
#ifdef _GET_IBM_COUNTER 
  hpmStop(1);
#endif

#undef __UNROLL_DEPTH_1

#else // original loop

#ifdef _GET_IBM_COUNTER 
  hpmStart(1, "minmax base");
#endif
  {
    size_t i;

#if   defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
    for ( i = 0; i < datasize; ++i )
      {
	if ( *fmin > data[i] ) *fmin = data[i];
	if ( *fmax < data[i] ) *fmax = data[i];
	/*
	 *fmin = *fmin < data[i] ? *fmin : data[i];
	 *fmax = *fmax > data[i] ? *fmax : data[i];
	 */
      }
  }
#ifdef _GET_IBM_COUNTER 
  hpmStop(1);
#endif

#endif // _ARCH_PWR6 && original loop
#endif // SIMD

#if defined _GET_X86_COUNTER || defined _GET_MACH_COUNTER
#ifdef _GET_X86_COUNTER 
  end_minmax = _rdtsc();
#endif
#ifdef _GET_MACH_COUNTER 
  end_minmax = mach_absolute_time();
#endif
#if defined __AVX__
  printf("AVX minmax cycles:: %" PRIu64 "\n", 
	 end_minmax-start_minmax);
  fprintf (stderr, "AVX min: %lf max: %lf\n", *fmin, *fmax);
#elif defined __SSE2__
  printf("SSE2 minmax cycles:: %" PRIu64 "\n", 
	 end_minmax-start_minmax);
  fprintf (stderr, "SSE2 min: %lf max: %lf\n", *fmin, *fmax);
#else
  printf("loop minmax cycles:: %" PRIu64 "\n", 
	 end_minmax-start_minmax);
  fprintf (stderr, "loop min: %lf max: %lf\n", *fmin, *fmax);
#endif
#endif

  return;
}
//#undef _GET_X86_COUNTER
//#undef _GET_MACH_COUNTER
//#undef _GET_IBM_COUNTER
//#undef _ARCH_PWR6


#if   defined(__GNUC__) && (__GNUC__ >= 4)
#elif defined(__ICC)    && (__ICC >= 1100)
#elif defined(__clang__)
#else
#define DISABLE_SIMD
#endif

#if defined _GET_IBM_COUNTER
#elif defined _GET_X86_COUNTER
#elif defined _GET_MACH_COUNTER
#endif

//#define DISABLE_SIMD

#ifdef DISABLE_SIMD
#ifndef ENABLE_AVX
#undef __AVX__
#endif
#ifndef ENABLE_SSE4_1
#undef __SSE4_1__
#endif
#endif

#if defined __AVX__
#elif defined __SSE4_1__
#endif

#if defined __AVX__

static
void avx_encode_double_array_2byte(size_t datasize, 
				   unsigned char * restrict lGrib,
				   const double * restrict data, 
				   double zref, double factor, size_t *gz) 
{
  size_t i, j;
  const double *dval = data;
  __m128i *sgrib = (__m128i *) (lGrib+(*gz));

  const __m128i swap = _mm_set_epi8(14, 15, 12, 13, 10, 11, 8, 9, 6, 7, 4, 5, 2, 3, 0, 1);

  const __m256d c0 = _mm256_set1_pd(zref);
  const __m256d c1 = _mm256_set1_pd(factor);
  const __m256d c2 = _mm256_set1_pd(0.5);
  
  __m256d d0, d3, d2, d1;
  __m128i i0, i1, i2, i3;
  __m128i s0, s1;  

  for (i = 0; i < datasize; i += 16)
    {
      (void) _mm_prefetch(dval+8, _MM_HINT_NTA);
      //_____________________________________________________________________________

      d0 = _mm256_loadu_pd (dval);
      d0 = _mm256_sub_pd (d0, c0);
      d0 = _mm256_mul_pd (d0, c1);
      d0 = _mm256_add_pd (d0, c2);

      i0 = _mm256_cvttpd_epi32 (d0);
      
      //_____________________________________________________________________________
      
      d1 = _mm256_loadu_pd (dval+4);
      d1 = _mm256_sub_pd (d1, c0);
      d1 = _mm256_mul_pd (d1, c1);
      d1 = _mm256_add_pd (d1, c2);
      
      i1 = _mm256_cvttpd_epi32 (d1);

      //_____________________________________________________________________________

      s0 = _mm_packus_epi32(i0, i1);
      s0 = _mm_shuffle_epi8 (s0, swap);
      (void) _mm_storeu_si128 (sgrib, s0);

      //_____________________________________________________________________________

      (void) _mm_prefetch(dval+16, _MM_HINT_NTA);

      //_____________________________________________________________________________
      
      d2 = _mm256_loadu_pd (dval+8);
      d2 = _mm256_sub_pd (d2, c0);
      d2 = _mm256_mul_pd (d2, c1);
      d2 = _mm256_add_pd (d2, c2);
      
      i2 = _mm256_cvttpd_epi32 (d2);

      //_____________________________________________________________________________
      
      d3 = _mm256_loadu_pd (dval+12);
      d3 = _mm256_sub_pd (d3, c0);
      d3 = _mm256_mul_pd (d3, c1);
      d3 = _mm256_add_pd (d3, c2);
      
      i3 = _mm256_cvttpd_epi32 (d3);

      //_____________________________________________________________________________

      s1 = _mm_packus_epi32(i2, i3);
      s1 = _mm_shuffle_epi8 (s1, swap);
      (void) _mm_storeu_si128 (sgrib+1, s1);

      //_____________________________________________________________________________
           
      dval += 16;
      sgrib += 2;
    }

  if (i != datasize)
    {
      uint16_t ui16;
      i -= 16;
      for ( j = i; j < datasize; j++ )
	{
	  ui16 = (uint16_t) ((data[j] - zref) * factor + 0.5);
	  lGrib[*gz+2*j  ] = ui16 >>  8;
	  lGrib[*gz+2*j+1] = ui16;
	}
    }
  
  *gz += 2*datasize;

  return;
}

#elif defined __SSE4_1__

static
void sse41_encode_double_array_2byte(size_t datasize, 
				     unsigned char * restrict lGrib,
				     const double * restrict data, 
				     double zref, double factor, size_t *gz) 
{
  size_t i, j;
  const double *dval = data;
  __m128i *sgrib = (__m128i *) (lGrib+(*gz));

  const __m128i swap = _mm_set_epi8(14, 15, 12, 13, 10, 11, 8, 9, 6, 7, 4, 5, 2, 3, 0, 1);

  const __m128d c0 = _mm_set1_pd(zref);
  const __m128d c1 = _mm_set1_pd(factor);
  const __m128d c2 = _mm_set1_pd(0.5);
  
  __m128d d0, d4, d3, d2, d1;
  __m128i i0, i1, i2, i3, i4;
  __m128i s0, s1;  

  for (i = 0; i < datasize; i += 16)
    {
      (void) _mm_prefetch(dval+8, _MM_HINT_NTA);
      //_____________________________________________________________________________

      d0 = _mm_loadu_pd (dval);
      d0 = _mm_sub_pd (d0, c0);
      d0 = _mm_mul_pd (d0, c1);
      d0 = _mm_add_pd (d0, c2);
      
      d4 = _mm_loadu_pd (dval+2);
      d4 = _mm_sub_pd (d4, c0);
      d4 = _mm_mul_pd (d4, c1);
      d4 = _mm_add_pd (d4, c2);

      i0 = _mm_cvttpd_epi32 (d0);
      i4 = _mm_cvttpd_epi32 (d4);  
      i0 = _mm_unpacklo_epi64 (i0, i4);

      //_____________________________________________________________________________
      
      d1 = _mm_loadu_pd (dval+4);
      d1 = _mm_sub_pd (d1, c0);
      d1 = _mm_mul_pd (d1, c1);
      d1 = _mm_add_pd (d1, c2);
      
      d4 = _mm_loadu_pd (dval+6);
      d4 = _mm_sub_pd (d4, c0);
      d4 = _mm_mul_pd (d4, c1);
      d4 = _mm_add_pd (d4, c2);
      
      i1 = _mm_cvttpd_epi32 (d1);
      i4 = _mm_cvttpd_epi32 (d4);  
      i1 = _mm_unpacklo_epi64 (i1, i4);

      //_____________________________________________________________________________

      s0 = _mm_packus_epi32(i0, i1);
      s0 = _mm_shuffle_epi8 (s0, swap);
      (void) _mm_storeu_si128 (sgrib, s0);

      //_____________________________________________________________________________

      (void) _mm_prefetch(dval+16, _MM_HINT_NTA);

      //_____________________________________________________________________________
      
      d2 = _mm_loadu_pd (dval+8);
      d2 = _mm_sub_pd (d2, c0);
      d2 = _mm_mul_pd (d2, c1);
      d2 = _mm_add_pd (d2, c2);
      
      d4 = _mm_loadu_pd (dval+10);
      d4 = _mm_sub_pd (d4, c0);
      d4 = _mm_mul_pd (d4, c1);
      d4 = _mm_add_pd (d4, c2);
      
      i2 = _mm_cvttpd_epi32 (d2);
      i4  = _mm_cvttpd_epi32 (d4);  
      i2 = _mm_unpacklo_epi64 (i2, i4);

      //_____________________________________________________________________________
      
      d3 = _mm_loadu_pd (dval+12);
      d3 = _mm_sub_pd (d3, c0);
      d3 = _mm_mul_pd (d3, c1);
      d3 = _mm_add_pd (d3, c2);
      
      d4 = _mm_loadu_pd (dval+14);
      d4 = _mm_sub_pd (d4, c0);
      d4 = _mm_mul_pd (d4, c1);
      d4 = _mm_add_pd (d4, c2);
      
      i3 = _mm_cvttpd_epi32 (d3);
      i4 = _mm_cvttpd_epi32 (d4);  
      i3 = _mm_unpacklo_epi64 (i3, i4);

      //_____________________________________________________________________________

      s1 = _mm_packus_epi32(i2, i3);
      s1 = _mm_shuffle_epi8 (s1, swap);
      (void) _mm_storeu_si128 (sgrib+1, s1);

      //_____________________________________________________________________________
           
      dval += 16;
      sgrib += 2;
    }

  if (i != datasize) 
    {
      uint16_t ui16;
      i -= 16;  
      for ( j = i; j < datasize; j++ )
	{
	  ui16 = (uint16_t) ((data[j] - zref) * factor + 0.5);
	  lGrib[*gz+2*j  ] = ui16 >>  8;
	  lGrib[*gz+2*j+1] = ui16;
	}
    }

  *gz += 2*datasize;
  
  return;
}

#endif // SIMD variants

static
void encode_double_array_common(int numBits, size_t packStart, size_t datasize, GRIBPACK *lGrib,
				const double *data, double zref, double factor, size_t *gz)
{
  size_t i, z = *gz;
  unsigned int ival;
  int cbits, jbits;
  unsigned int c;
  static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
    
  /* code from gribw routine flist2bitstream */

  cbits = 8;
  c = 0;
  for ( i = packStart; i < datasize; i++ )
    {
      /* note float -> unsigned int .. truncate */
      ival = (unsigned int) ((data[i] - zref) * factor + 0.5);
      /*
	if ( ival > max_nbpv_pow2 ) ival = max_nbpv_pow2;
	if ( ival < 0 ) ival = 0;
      */
      jbits = numBits;
      while ( cbits <= jbits ) 
	{
	  if ( cbits == 8 )
	    {
	      jbits -= 8;
	      lGrib[z++] = (ival >> jbits) & 0xFF;
	    }
	  else
	    {
	      jbits -= cbits;
	      lGrib[z++] = (c << cbits) + ((ival >> jbits) & mask[cbits]);
	      cbits = 8;
	      c = 0;
	    }
	}
      /* now jbits < cbits */
      if ( jbits )
	{
	  c = (c << jbits) + (ival & mask[jbits]);
	  cbits -= jbits;
	}
    }
  if ( cbits != 8 ) lGrib[z++] = c << cbits;

  *gz = z;
}

static
void encode_double_array_byte(int numBits, size_t packStart, size_t datasize, 
			      GRIBPACK *restrict lGrib,
			      const double *restrict data, 
			      double zref, double factor, size_t *gz)
{
#if defined _GET_X86_COUNTER || defined _GET_MACH_COUNTER 
  uint64_t start_minmax, end_minmax;
#endif

  uint16_t ui16;
  uint32_t ui32;
  size_t i, z = *gz;
  double tmp;

  data += packStart;
  datasize -= packStart;

  if      ( numBits ==  8 )
    {
#ifdef _GET_IBM_COUNTER 
      hpmStart(2, "pack 8 bit base");
#endif

#if defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
      for ( i = 0; i < datasize; i++ )
	{
	  tmp = ((data[i] - zref) * factor + 0.5);
	  lGrib[z  ] = (uint16_t) tmp;
          z++;
	}

#ifdef _GET_IBM_COUNTER 
      hpmStop(2);
#endif
    }
  else if ( numBits == 16 )
    {
#ifdef _GET_IBM_COUNTER 
      hpmStart(3, "pack 16 bit base");
#elif defined _GET_X86_COUNTER 
  start_minmax = _rdtsc();
#elif defined _GET_MACH_COUNTER 
  start_minmax = mach_absolute_time();
#endif

#if defined __AVX__

      avx_encode_double_array_2byte(datasize, lGrib, data, zref, factor, &z);

#elif defined __SSE4_1__

      sse41_encode_double_array_2byte(datasize, lGrib, data, zref, factor, &z);

#else

#if   defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
      for ( i = 0; i < datasize; i++ )
      	{
      	  tmp = ((data[i] - zref) * factor + 0.5);
          ui16 = (uint16_t) tmp;
          lGrib[z  ] = ui16 >>  8;
          lGrib[z+1] = ui16;
          z += 2;
      	}
#endif

#if defined _GET_X86_COUNTER || defined _GET_MACH_COUNTER
#if defined _GET_X86_COUNTER 
      end_minmax = _rdtsc();
#elif defined _GET_MACH_COUNTER 
      end_minmax = mach_absolute_time();
#endif
#if defined __AVX__
      printf("AVX encoding cycles:: %" PRIu64 "\n", 
	     end_minmax-start_minmax);
#elif defined __SSE4_1__
      printf("SSE 4.1 encoding cycles:: %" PRIu64 "\n", 
	     end_minmax-start_minmax);
#else
      printf("loop encoding cycles:: %" PRIu64 "\n", 
	     end_minmax-start_minmax);
#endif  
#endif
      
#ifdef _GET_IBM_COUNTER 
      hpmStop(3);
#endif
    }
  else if ( numBits == 24 )
    {
#ifdef _GET_IBM_COUNTER 
      hpmStart(4, "pack 24 bit base");
#endif

#if   defined (CRAY)
#pragma _CRI ivdep
#elif defined (SX)
#pragma vdir nodep
#elif defined (__uxp__)
#pragma loop novrec
#endif
      for ( i = 0; i < datasize; i++ )
	{