Philippe 02/05/2016: moved all LIBTOOLS files in LIBTOOLS directory
[MNH-git_open_source-lfs.git] / LIBTOOLS / lib / vis5d / src / v5d.c
1 /* v5d.c */
2
3 /* Vis5D version 5.1 */
4
5 /*
6 Vis5D system for visualizing five dimensional gridded data sets
7 Copyright (C) 1990 - 1997 Bill Hibbard, Johan Kellum, Brian Paul,
8 Dave Santek, and Andre Battaiola.
9
10 This program is free software; you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation; either version 1, or (at your option)
13 any later version.
14
15 This program is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with this program; if not, write to the Free Software
22 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
23 */
24
25
26 /* this should be updated when the file version changes */
27 #define FILE_VERSION "4.3"
28
29
30
31 /*
32  * New grid file format for VIS-5D:
33  *
34  * The header is a list of tagged items.  Each item has 3 parts:
35  *    1. A tag which is a 4-byte integer identifying the type of item.
36  *    2. A 4-byte integer indicating how many bytes of data follow.
37  *    3. The binary data.
38  *
39  * If we need to add new information to a file header we just create a
40  * new tag and add the code to read/write the information.
41  *
42  * If we're reading a header and find an unknown tag, we can use the
43  * length field to skip ahead to the next tag.  Therefore, the file
44  * format is forward (and backward) compatible.
45  *
46  * Grid data is stored as either:
47  *     1-byte unsigned integers  (255=missing)
48  *     2-byte unsigned integers  (65535=missing)
49  *     4-byte IEEE floats        ( >1.0e30 = missing) 
50  *
51  * All numeric values are stored in big endian order.  All floating point
52  * values are in IEEE format.
53  */
54
55
56
57 /*
58  * Updates:
59  *
60  * April 13, 1995, brianp
61  *   finished Cray support for 2-byte and 4-byte compress modes
62  */
63
64
65
66
67 #include <assert.h>
68 #include <stdio.h>
69 #include <stdlib.h>
70 #include <string.h>
71 #include <math.h>
72 #include "binio.h"
73 #include "v5d.h"
74 #include "vis5d.h"
75 #ifndef SEEK_SET
76 #  define SEEK_SET 0
77 #endif
78 #ifndef SEEK_CUR
79 #  define SEEK_CUR 1
80 #endif
81 #ifndef SEEK_END
82 #  define SEEK_END 2
83 #endif
84
85
86
87 /*
88  * Currently defined tags:
89  * Note:  the notation a[i] doesn't mean a is an array of i elements,
90  * rather it just refers to the ith element of a[].
91  *
92  * Tags marked as PHASED OUT should be readable but are no longer written.
93  * Old tag numbers can't be reused!
94  * 
95  */
96
97
98 /*      TAG NAME        VALUE       DATA (comments)                     */
99 /*----------------------------------------------------------------------*/
100 #define TAG_ID          0x5635440a  /* hex encoding of "V5D\n"          */
101
102 /* general stuff 1000+ */
103 #define TAG_VERSION     1000        /* char*10 FileVersion              */
104 #define TAG_NUMTIMES    1001        /* int*4 NumTimes                   */
105 #define TAG_NUMVARS     1002        /* int*4 NumVars                    */
106 #define TAG_VARNAME     1003        /* int*4 var; char*10 VarName[var]  */
107
108 #define TAG_NR          1004        /* int*4 Nr                         */ 
109 #define TAG_NC          1005        /* int*4 Nc                         */
110 #define TAG_NL          1006        /* int*4 Nl  (Nl for all vars)      */
111 #define TAG_NL_VAR      1007        /* int*4 var; int*4 Nl[var]         */
112 #define TAG_LOWLEV_VAR  1008        /* int*4 var; int*4 LowLev[var]     */
113
114 #define TAG_TIME        1010        /* int*4 t;  int*4 TimeStamp[t]     */
115 #define TAG_DATE        1011        /* int*4 t;  int*4 DateStamp[t]     */
116
117 #define TAG_MINVAL      1012        /* int*4 var;  real*4 MinVal[var]   */
118 #define TAG_MAXVAL      1013        /* int*4 var;  real*4 MaxVal[var]   */
119
120 #define TAG_COMPRESS    1014        /* int*4 CompressMode; (#bytes/grid)*/
121
122 #define TAG_UNITS       1015        /* int *4 var; char*20 Units[var]   */
123
124 /* vertical coordinate system 2000+ */
125 #define TAG_VERTICAL_SYSTEM 2000    /* int*4 VerticalSystem             */
126 #define TAG_VERT_ARGS    2100       /* int*4 n;  real*4 VertArgs[0..n-1]*/
127
128 #define TAG_BOTTOMBOUND  2001       /* real*4 BottomBound     (PHASED OUT)  */
129 #define TAG_LEVINC       2002       /* real*4 LevInc      (PHASED OUT)      */
130 #define TAG_HEIGHT       2003    /* int*4 l;  real*4 Height[l] (PHASED OUT) */
131
132
133 /* projection 3000+ */
134 #define TAG_PROJECTION   3000       /* int*4 projection:                    */
135                                     /*   0 = generic linear                 */
136                                     /*   1 = cylindrical equidistant        */
137                                     /*   2 = Lambert conformal/Polar Stereo */
138                                     /*   3 = rotated equidistant            */
139 #define TAG_PROJ_ARGS    3100       /* int *4 n;  real*4 ProjArgs[0..n-1]   */
140
141 #define TAG_NORTHBOUND   3001       /* real*4 NorthBound       (PHASED OUT) */
142 #define TAG_WESTBOUND    3002       /* real*4 WestBound        (PHASED OUT) */
143 #define TAG_ROWINC       3003       /* real*4 RowInc           (PHASED OUT) */
144 #define TAG_COLINC       3004       /* real*4 ColInc           (PHASED OUT) */
145 #define TAG_LAT1         3005       /* real*4 Lat1             (PHASED OUT) */
146 #define TAG_LAT2         3006       /* real*4 Lat2             (PHASED OUT) */
147 #define TAG_POLE_ROW     3007       /* real*4 PoleRow          (PHASED OUT) */
148 #define TAG_POLE_COL     3008       /* real*4 PoleCol          (PHASED OUT) */
149 #define TAG_CENTLON      3009       /* real*4 CentralLon       (PHASED OUT) */
150 #define TAG_CENTLAT      3010       /* real*4 CentralLat       (PHASED OUT) */
151 #define TAG_CENTROW      3011       /* real*4 CentralRow       (PHASED OUT) */
152 #define TAG_CENTCOL      3012       /* real*4 CentralCol       (PHASED OUT) */
153 #define TAG_ROTATION     3013       /* real*4 Rotation         (PHASED OUT) */
154
155
156 #define TAG_END                9999
157
158
159
160
161
162
163 /**********************************************************************/
164 /*****                  Miscellaneous Functions                   *****/
165 /**********************************************************************/
166
167
168 float pressure_to_height(float pressure)
169 {
170   return (float) DEFAULT_LOG_EXP * log((double) pressure / DEFAULT_LOG_SCALE);
171 }
172
173 float height_to_pressure(float height)
174 {
175   return (float) DEFAULT_LOG_SCALE * exp((double) height / DEFAULT_LOG_EXP);
176 }
177
178
179 /*
180  * Return current file position.
181  * Input:  f - file descriptor
182  */
183 static off_t ltell( int f )
184 {
185    return lseek( f, 0, SEEK_CUR );
186 }
187
188
189 /*
190  * Copy up to maxlen characters from src to dst stopping upon whitespace
191  * in src.  Terminate dst with null character.
192  * Return:  length of dst.
193  */
194 static int copy_string2( char *dst, const char *src, int maxlen )
195 {
196    int i;
197
198    for (i=0;i<maxlen;i++) dst[i] = src[i];
199    for (i=maxlen-1; i>=0; i--) {
200      if (dst[i]==' ' || i==maxlen-1) dst[i] = 0;
201      else break;
202    }
203    return strlen(dst);
204 }
205
206
207
208 /*
209  * Copy up to maxlen characters from src to dst stopping upon whitespace
210  * in src.  Terminate dst with null character.
211  * Return:  length of dst.
212  */
213 static int copy_string( char *dst, const char *src, int maxlen )
214 {
215    int i;
216
217    for (i=0;i<maxlen;i++) {
218       if (src[i]==' ' || i==maxlen-1) {
219          dst[i] = 0;
220          break;
221       }
222       else {
223          dst[i] = src[i];
224       }
225    }
226    return i;
227 }
228
229
230
231 /*
232  * Convert a date from YYDDD format to days since Jan 1, 1900.
233  */
234 int v5dYYDDDtoDays( int yyddd )
235 {
236    int iy, id, idays;
237
238    iy = yyddd / 1000;
239    id = yyddd - 1000*iy;
240    if (iy < 50) iy += 100; /* WLH 31 July 96 << 31 Dec 99 */
241    idays = 365*iy + (iy-1)/4 + id;
242
243    return idays;
244 }
245
246
247 /*
248  * Convert a time from HHMMSS format to seconds since midnight.
249  */
250 int v5dHHMMSStoSeconds( int hhmmss )
251 {
252    int h, m, s;
253
254    h = hhmmss / 10000;
255    m = (hhmmss / 100) % 100;
256    s = hhmmss % 100;
257
258    return s + m*60 + h*60*60;
259 }
260
261
262
263 /*
264  * Convert a day since Jan 1, 1900 to YYDDD format.
265  */
266 int v5dDaysToYYDDD( int days )
267 {
268    int iy, id, iyyddd;
269
270    iy = (4*days)/1461;
271    id = days-(365*iy+(iy-1)/4);
272    if (iy > 99) iy = iy - 100; /* WLH 31 July 96 << 31 Dec 99 */
273    /* iy = iy + 1900; is the right way to fix this, but requires
274       changing all places where dates are printed - procrastinate */
275    iyyddd = iy*1000+id;
276
277    return iyyddd;
278 }
279
280
281 /*
282  * Convert a time in seconds since midnight to HHMMSS format.
283  */
284 int v5dSecondsToHHMMSS( int seconds )
285 {
286    int hh, mm, ss;
287
288    hh = seconds / (60*60);
289    mm = (seconds / 60) % 60;
290    ss = seconds % 60;
291    return hh*10000 + mm * 100 + ss;
292 }
293
294
295
296
297 void v5dPrintStruct( const v5dstruct *v )
298 {
299    static char day[7][10] = { "Sunday", "Monday", "Tuesday", "Wednesday",
300                               "Thursday", "Friday", "Saturday" };
301    int time, var, i;
302    int maxnl;
303
304    maxnl = 0;
305    for (var=0;var<v->NumVars;var++) {
306       if (v->Nl[var]+v->LowLev[var]>maxnl) {
307          maxnl = v->Nl[var]+v->LowLev[var];
308       }
309    }
310
311    if (v->FileFormat==0) {
312       if (v->FileVersion[0] == 0) {
313         printf("File format: v5d  version: (4.0 or 4.1)\n");
314       }
315       else {
316         printf("File format: v5d  version: %s\n", v->FileVersion);
317       }
318    }
319    else {
320       printf("File format: comp5d  (VIS-5D 3.3 or older)\n");
321    }
322
323    if (v->CompressMode==1) {
324       printf("Compression:  1 byte per gridpoint.\n");
325    }
326    else {
327       printf("Compression:  %d bytes per gridpoint.\n", v->CompressMode);
328    }
329    printf("header size=%d\n", v->FirstGridPos);
330    printf("sizeof(v5dstruct)=%d\n", sizeof(v5dstruct) );
331    printf("\n");
332
333    printf("NumVars = %d\n", v->NumVars );
334
335    printf("Var  Name       Units      Rows  Cols  Levels LowLev  MinVal       MaxVal\n");
336    for (var=0;var<v->NumVars;var++) {
337       printf("%3d  %-10s %-10s %3d   %3d   %3d    %3d",
338              var+1, v->VarName[var], v->Units[var],
339              v->Nr, v->Nc, v->Nl[var], v->LowLev[var] );
340       if (v->MinVal[var] > v->MaxVal[var]) {
341          printf("     MISSING      MISSING\n");
342       }
343       else {
344          printf("     %-12g %-12g\n", v->MinVal[var], v->MaxVal[var] );
345       }
346    }
347
348    printf("\n");
349
350    printf("NumTimes = %d\n", v->NumTimes );
351    printf("Step    Date(YYDDD)    Time(HH:MM:SS)   Day\n");
352    for (time=0;time<v->NumTimes;time++) {
353       int i = v->TimeStamp[time];
354       printf("%3d        %05d       %5d:%02d:%02d     %s\n",
355              time+1,
356              v->DateStamp[time],
357              i/10000, (i/100)%100, i%100,
358              day[ v5dYYDDDtoDays(v->DateStamp[time]) % 7 ]);
359    }
360    printf("\n");
361
362    switch (v->VerticalSystem) {
363       case 0:
364          printf("Generic linear vertical coordinate system:\n");
365          printf("\tBottom Bound: %f\n", v->VertArgs[0] );
366          printf("\tIncrement between levels:  %f\n", v->VertArgs[1] );
367          break;
368       case 1:
369          printf("Equally spaced levels in km:\n");
370          printf("\tBottom Bound: %f\n", v->VertArgs[0] );
371          printf("\tIncrement: %f\n", v->VertArgs[1] );
372          break;
373       case 2:
374          printf("Unequally spaced levels in km:\n");
375          printf("Level\tHeight(km)\n");
376          for (i=0;i<maxnl;i++) {
377             printf("%3d     %10.3f\n", i+1, v->VertArgs[i] );
378          }
379          break;
380       case 3:
381          printf("Unequally spaced levels in mb:\n");
382          printf("Level\tPressure(mb)\n");
383          for (i=0;i<maxnl;i++) {
384             printf("%3d     %10.3f\n", i+1, height_to_pressure(v->VertArgs[i]) );
385          }
386          break;
387       default:
388          printf("Bad VerticalSystem value: %d\n", v->VerticalSystem );
389    }
390    printf("\n");
391
392    switch (v->Projection) {
393       case 0:
394          printf("Generic linear projection:\n");
395          printf("\tNorth Boundary: %f\n", v->ProjArgs[0] );
396          printf("\tWest Boundary: %f\n", v->ProjArgs[1] );
397          printf("\tRow Increment: %f\n", v->ProjArgs[2] );
398          printf("\tColumn Increment: %f\n", v->ProjArgs[3] );
399          break;
400       case 1:
401          printf("Cylindrical Equidistant projection:\n");
402          printf("\tNorth Boundary: %f degrees\n", v->ProjArgs[0] );
403          printf("\tWest Boundary: %f degrees\n", v->ProjArgs[1] );
404          printf("\tRow Increment: %f degrees\n", v->ProjArgs[2] );
405          printf("\tColumn Increment: %f degrees\n", v->ProjArgs[3] );
406 /*
407          printf("\tSouth Boundary: %f degrees\n",
408                 v->NorthBound - v->RowInc * (v->Nr-1) );
409          printf("\tEast Boundary: %f degrees\n",
410                 v->WestBound - v->ColInc * (v->Nc-1) );
411 */
412          break;
413       case 2:
414          printf("Lambert Conformal projection:\n");
415          printf("\tStandard Latitude 1: %f\n", v->ProjArgs[0] );
416          printf("\tStandard Latitude 2: %f\n", v->ProjArgs[1] );
417          printf("\tNorth/South Pole Row: %f\n", v->ProjArgs[2] );
418          printf("\tNorth/South Pole Column: %f\n", v->ProjArgs[3] );
419          printf("\tCentral Longitude: %f\n", v->ProjArgs[4] );
420          printf("\tColumn Increment: %f km\n", v->ProjArgs[5] );
421          break;
422       case 3:
423          printf("Stereographic:\n");
424          printf("\tCenter Latitude: %f\n", v->ProjArgs[0] );
425          printf("\tCenter Longitude: %f\n", v->ProjArgs[1] );
426          printf("\tCenter Row: %f\n", v->ProjArgs[2] );
427          printf("\tCenter Column: %f\n", v->ProjArgs[3] );
428          printf("\tColumn Spacing: %f\n", v->ProjArgs[4] );
429          break;
430       case 4:
431          /* WLH 4-21-95 */
432          printf("Rotated equidistant projection:\n");
433          printf("\tLatitude of grid(0,0): %f\n", v->ProjArgs[0] );
434          printf("\tLongitude of grid(0,0): %f\n", v->ProjArgs[1] );
435          printf("\tRow Increment: %f degress\n", v->ProjArgs[2] );
436          printf("\tColumn Increment: %f degrees\n", v->ProjArgs[3] );
437          printf("\tCenter Latitude: %f\n", v->ProjArgs[4] );
438          printf("\tCenter Longitude: %f\n", v->ProjArgs[5] );
439          printf("\tRotation: %f degrees\n", v->ProjArgs[6] );
440          break;
441       default:
442          printf("Bad projection number: %d\n", v->Projection );
443    }
444 }
445
446
447
448 /*
449  * Compute the location of a compressed grid within a file.
450  * Input:  v - pointer to v5dstruct describing the file header.
451  *         time, var - which timestep and variable.
452  * Return:  file offset in bytes
453  */
454 static int grid_position( const v5dstruct *v, int time, int var )
455 {
456    int pos, i;
457
458    assert( time >= 0 );
459    assert( var >= 0 );
460    assert( time < v->NumTimes );
461    assert( var < v->NumVars );
462
463    pos = v->FirstGridPos + time * v->SumGridSizes;
464    for (i=0;i<var;i++) {
465       pos += v->GridSize[i];
466    }
467
468    return pos;
469 }
470
471
472
473 /*
474  * Compute the ga and gb (de)compression values for a grid.
475  * Input:  nr, nc, nl - size of grid
476  *         data - the grid data
477  *         ga, gb - arrays to store results.
478  *         minval, maxval - pointer to floats to return min, max values
479  *         compressmode - 1, 2 or 4 bytes per grid point
480  * Output:  ga, gb - the (de)compression values
481  *          minval, maxval - the min and max grid values
482  * Side effect:  the MinVal[var] and MaxVal[var] fields in g may be
483  *               updated with new values.
484  */
485 static void compute_ga_gb( int nr, int nc, int nl, 
486                            const float data[], int compressmode,
487                            float ga[], float gb[],
488                            float *minval, float *maxval )
489 {
490 #ifdef SIMPLE_COMPRESSION
491    /*
492     * Compute ga, gb values for whole grid.
493     */
494    int i, lev, allmissing, num;
495    float min, max, a, b;
496
497    min = 1.0e30;
498    max = -1.0e30;
499    num = nr * nc * nl;
500    allmissing = 1;
501    for (i=0;i<num;i++) {
502       if (!IS_MISSING(data[i])) {
503          if (data[i]<min)  min = data[i];
504          if (data[i]>max)  max = data[i];
505          allmissing = 0;
506       }
507    }
508    if (allmissing) {
509       a = 1.0;
510       b = 0.0;
511    }
512    else {
513       a = (max-min) / 254.0;
514       b = min;
515    }
516
517    /* return results */
518    for (i=0;i<nl;i++) {
519       ga[i] = a;
520       gb[i] = b;
521    }
522
523    *minval = min;
524    *maxval = max;
525 #else
526    /*
527     * Compress grid on level-by-level basis.
528     */
529 #  define SMALLVALUE -1.0e30
530 #  define BIGVALUE 1.0e30
531 #  define ABS(x)   ( ((x) < 0.0) ? -(x) : (x) )
532    float gridmin, gridmax;
533    float levmin[MAXLEVELS], levmax[MAXLEVELS];
534    float d[MAXLEVELS], dmax;
535    float ival, mval;
536    int j, k, lev, nrnc;
537
538    nrnc = nr * nc;
539
540    /* find min and max for each layer and the whole grid */
541    gridmin = BIGVALUE;
542    gridmax = SMALLVALUE;
543    j = 0;
544
545
546    for (lev=0;lev<nl;lev++) {
547       float ave, var;
548       float min, max;
549       min = BIGVALUE;
550       max = SMALLVALUE;
551       ave = 0.0;
552       var = 0.0;
553       for (k=0;k<nrnc;k++) {
554          if (!IS_MISSING(data[j]) && data[j]<min)
555             min = data[j];
556          if (!IS_MISSING(data[j]) && data[j]>max)
557             max = data[j];
558          j++;
559       }
560
561       if (min<gridmin)
562         gridmin = min;
563       if (max>gridmax)
564         gridmax = max;
565       levmin[lev] = min;
566       levmax[lev] = max;
567    }
568
569 /* WLH 2-2-95 */
570 #ifdef KLUDGE
571    /* if the grid minimum is within delt of 0.0, fudge all values */
572    /* within delt of 0.0 to delt, and recalculate mins and maxes */
573    {
574       float delt;
575       int nrncnl = nrnc * nl;
576
577       delt = (gridmax - gridmin)/100000.0;
578       if ( ABS(gridmin) < delt && gridmin!=0.0 && compressmode != 4 ) {
579          float min, max;
580          for (j=0; j<nrncnl; j++) {
581             if (!IS_MISSING(data[j]) && data[j]<delt)
582               data[j] = delt;
583          }
584          /* re-calculate min and max for each layer and the whole grid */
585          gridmin = delt;
586          for (lev=0;lev<nl;lev++) {
587             if (ABS(levmin[lev]) < delt)
588               levmin[lev] = delt;
589             if (ABS(levmax[lev]) < delt)
590               levmax[lev] = delt;
591          }
592       }
593    }
594 #endif
595
596    /* find d[lev] and dmax = MAX( d[0], d[1], ... d[nl-1] ) */
597    dmax = 0.0;
598    for (lev=0;lev<nl;lev++) {
599       if (levmin[lev]>=BIGVALUE && levmax[lev]<=SMALLVALUE) {
600          /* all values in the layer are MISSING */
601          d[lev] = 0.0;
602       }
603       else {
604          d[lev] = levmax[lev]-levmin[lev];
605       }
606       if (d[lev]>dmax)
607          dmax = d[lev];
608    }
609
610    /*** Compute ga (scale) and gb (bias) for each grid level */
611    if (dmax==0.0) {
612       /*** Special cases ***/
613       if (gridmin==gridmax) {
614          /*** whole grid is of same value ***/
615          for (lev=0; lev<nl; lev++) {
616             ga[lev] = gridmin;
617             gb[lev] = 0.0;
618          }
619       }
620       else {
621          /*** every layer is of a single value ***/
622          for (lev=0; lev<nl; lev++) {
623             ga[lev] = levmin[lev];
624             gb[lev] = 0.0;
625          }
626       }
627    }
628    else {
629       /*** Normal cases ***/
630       if (compressmode == 1) {
631 #define ORIGINAL
632 #ifdef ORIGINAL
633          ival = dmax / 254.0;
634          mval = gridmin;
635
636          for (lev=0; lev<nl; lev++) {
637             ga[lev] = ival;
638             gb[lev] = mval + ival * (int) ( (levmin[lev]-mval) / ival ); 
639          }
640 #else
641          for (lev=0; lev<nl; lev++) {
642             if (d[lev]==0.0) {
643                ival = 1.0;
644             }
645             else {
646                ival = d[lev] / 254.0;
647             }
648             ga[lev] = ival;
649             gb[lev] = levmin[lev];
650          }
651 #endif
652       }
653       else if (compressmode == 2) {
654          ival = dmax / 65534.0;
655          mval = gridmin;
656
657          for (lev=0; lev<nl; lev++) {
658             ga[lev] = ival;
659             gb[lev] = mval + ival * (int) ( (levmin[lev]-mval) / ival );
660          }
661       }
662       else {
663          assert( compressmode==4 );
664          for (lev=0; lev<nl; lev++) {
665             ga[lev] = 1.0;
666             gb[lev] = 0.0;
667          }
668       }
669    }
670
671    /* update min, max values */
672    *minval = gridmin;
673    *maxval = gridmax;
674 #endif
675 }
676
677
678
679
680 /*
681  * Compress a 3-D grid from floats to 1-byte unsigned integers.
682  * Input: nr, nc, nl - size of grid
683  *        compressmode - 1, 2 or 4 bytes per grid point
684  *        data - array of [nr*nc*nl] floats
685  *        compdata - pointer to array of [nr*nc*nl*compressmode] bytes
686  *                   to put results into.
687  *        ga, gb - pointer to arrays to put ga and gb decompression values
688  *        minval, maxval - pointers to float to return min & max values
689  * Output:  compdata - the compressed grid data
690  *          ga, gb - the decompression values
691  *          minval, maxval - the min and max grid values
692  */
693 void v5dCompressGrid( int nr, int nc, int nl, int compressmode,
694                       const float data[],
695                       void *compdata, float ga[], float gb[],
696                       float *minval, float *maxval )
697 {
698    int nrnc = nr * nc;
699    int nrncnl = nr * nc * nl;
700    V5Dubyte *compdata1 = (V5Dubyte *) compdata;
701    V5Dushort *compdata2 = (V5Dushort *) compdata;
702
703    /* compute ga, gb values */
704    compute_ga_gb( nr, nc, nl, data, compressmode, ga, gb, minval, maxval );
705
706    /* compress the data */
707    if (compressmode==1) {
708       int i, lev, p;
709       p = 0;
710       for (lev=0;lev<nl;lev++) {
711          float one_over_a, b;
712 /* WLH 5 Nov 98
713          b = gb[lev] - 0.0001;
714 */
715          /* WLH 5 Nov 98 */
716          b = gb[lev];
717                                 /* subtract an epsilon so the int((d-b)/a) */
718                                 /* expr below doesn't get mis-truncated. */
719          if (ga[lev]==0.0) {
720             one_over_a = 1.0;
721          }
722          else {
723             one_over_a = 1.0 / ga[lev];
724          }
725          for (i=0;i<nrnc;i++,p++) {
726             if (IS_MISSING(data[p])) {
727                compdata1[p] = 255;
728             }
729             else {
730 /* MJK 1.19.99
731                compdata1[p] = (V5Dubyte) (int) ((data[p]-b) * one_over_a);
732 */
733                compdata1[p] = (V5Dubyte) rint((data[p]-b) * one_over_a);
734                if (compdata1[p] >= 255){
735                   compdata1[p] = (V5Dubyte) (int) (255.0 - .0001);
736                }
737             }
738          }
739       }
740    }
741
742    else if (compressmode == 2) {
743       int i, lev, p;
744       p = 0;
745       for (lev=0;lev<nl;lev++) {
746          float one_over_a, b;
747 /* WLH 5 Nov 98
748          b = gb[lev] - 0.0001;
749 */
750          /* WLH 5 Nov 98 */
751          b = gb[lev];
752
753          if (ga[lev]==0.0) {
754             one_over_a = 1.0;
755          }
756          else {
757             one_over_a = 1.0 / ga[lev];
758          }
759 #ifdef _CRAY
760          /* this is tricky because sizeof(V5Dushort)==8, not 2 */
761          for (i=0;i<nrnc;i++,p++) {
762             V5Dushort compvalue;
763             if (IS_MISSING(data[p])) {
764                compvalue = 65535;
765             }
766             else {
767 /* MJK 3.2.99
768                compvalue = (V5Dushort) (int) ((data[p]-b) * one_over_a);
769 */
770                compvalue = (V5Dushort) rint((data[p]-b) * one_over_a);
771             }
772             compdata1[p*2+0] = compvalue >> 8;     /* upper byte */
773             compdata1[p*2+1] = compvalue & 0xffu;  /* lower byte */
774          }
775 #else
776          for (i=0;i<nrnc;i++,p++) {
777             if (IS_MISSING(data[p])) {
778                compdata2[p] = 65535;
779             }
780             else {
781                compdata2[p] = (V5Dushort) rint((data[p]-b) * one_over_a);
782
783 /*
784                compdata2[p] = (V5Dushort) (int) ((data[p]-b) * one_over_a);
785 */
786 /* MJK 3.24.99 I put this here so if the value is close
787    to the missing value and get's rounded up it won't come out
788    as missing data */
789                if (compdata2[p] == 65535){
790                   compdata2[p] = 65534;
791                }
792             }
793          }
794          /* TODO: byte-swapping on little endian??? */
795 #endif
796       }
797    }
798
799    else {
800       /* compressmode==4 */
801 #ifdef _CRAY
802       cray_to_ieee_array( compdata, data, nrncnl );
803 #else
804       /* other machines: just copy 4-byte IEEE floats */
805       assert( sizeof(float)==4 );
806       memcpy( compdata, data, nrncnl*4 );
807       /* TODO: byte-swapping on little endian??? */
808 #endif
809    }
810 }
811
812
813
814 /*
815  * Decompress a 3-D grid from 1-byte integers to 4-byte floats.
816  * Input:  nr, nc, nl - size of grid
817  *         compdata - array of [nr*nr*nl*compressmode] bytes
818  *         ga, gb - arrays of decompression factors
819  *         compressmode - 1, 2 or 4 bytes per grid point
820  *         data - address to put decompressed values
821  * Output:  data - uncompressed floating point data values
822  */
823 void v5dDecompressGrid( int nr, int nc, int nl, int compressmode,
824                         void *compdata, float ga[], float gb[],
825                         float data[] )
826 {
827    int nrnc = nr * nc;
828    int nrncnl = nr * nc * nl;
829    V5Dubyte *compdata1 = (V5Dubyte *) compdata;
830    V5Dushort *compdata2 = (V5Dushort *) compdata;
831
832    if (compressmode == 1) {
833       int p, i, lev;
834       p = 0;
835       for (lev=0;lev<nl;lev++) {
836          float a = ga[lev];
837          float b = gb[lev];
838
839          /* WLH 2-2-95 */
840          float d, aa;
841          int id;
842          if (a > 0.0000000001) {
843            d = b / a;
844            id = floor(d);
845            d = d - id;
846            aa = a * 0.000001;
847          }
848          else {
849            id = 1;
850          }
851          if (-254 <= id && id <= 0 && d < aa) {
852            for (i=0;i<nrnc;i++,p++) {
853               if (compdata1[p]==255) {
854                  data[p] = MISSING;
855               }
856               else {
857                  data[p] = (float) (int) compdata1[p] * a + b;
858                  if (fabs(data[p]) < aa) data[p] = aa;
859               }
860            }
861          }
862          else {
863            for (i=0;i<nrnc;i++,p++) {
864               if (compdata1[p]==255) {
865                  data[p] = MISSING;
866               }
867               else {
868                  data[p] = (float) (int) compdata1[p] * a + b;
869               }
870            }
871          }
872          /* end of WLH 2-2-95 */
873       }
874    }
875
876    else if (compressmode == 2) {
877       int p, i, lev;
878       p = 0;
879       for (lev=0;lev<nl;lev++) {
880          float a = ga[lev];
881          float b = gb[lev];
882 #ifdef _CRAY
883          /* this is tricky because sizeof(V5Dushort)==8, not 2 */
884          for (i=0;i<nrnc;i++,p++) {
885             int compvalue;
886             compvalue = (compdata1[p*2] << 8) | compdata1[p*2+1];
887             if (compvalue==65535) {
888                data[p] = MISSING;
889             }
890             else {
891                data[p] = (float) compvalue * a + b;
892             }
893          }
894 #else
895          /* sizeof(V5Dushort)==2! */
896          for (i=0;i<nrnc;i++,p++) {
897             if (compdata2[p]==65535) {
898                data[p] = MISSING;
899             }
900             else {
901                data[p] = (float) (int) compdata2[p] * a + b;
902             }
903          }
904 #endif
905       }
906    }
907
908    else {
909       /* compressmode==4 */
910 #ifdef _CRAY
911       ieee_to_cray_array( data, compdata, nrncnl );
912 #else
913       /* other machines: just copy 4-byte IEEE floats */
914       assert( sizeof(float)==4 );
915       memcpy( data, compdata, nrncnl*4 );
916 #endif
917    }
918 }
919
920
921
922
923 /*
924  * Return the size (in bytes) of the 3-D grid specified by time and var.
925  * Input:  v - pointer to v5dstruct describing the file
926  *         time, var - which timestep and variable
927  * Return:  number of data points.
928  */
929 int v5dSizeofGrid( const v5dstruct *v, int time, int var )
930 {
931    return v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
932 }
933
934
935
936 /*
937  * Initialize a v5dstructure to reasonable initial values.
938  * Input:  v - pointer to v5dstruct.
939  */
940 void v5dInitStruct( v5dstruct *v )
941 {
942    int i;
943
944    /* set everything to zero */
945    memset( v, 0, sizeof(v5dstruct) );
946
947    /* special cases */
948    v->Projection = -1;
949    v->VerticalSystem = -1;
950
951    for (i=0;i<MAXVARS;i++) {
952       v->MinVal[i] = MISSING;
953       v->MaxVal[i] = -MISSING;
954       v->LowLev[i] = 0;
955    }
956
957    /* set file version */
958    strcpy(v->FileVersion, FILE_VERSION);
959
960    v->CompressMode = 1;
961    v->FileDesc = -1;
962 }
963
964
965
966 /*
967  * Return a pointer to a new, initialized v5dstruct.
968  */
969 v5dstruct *v5dNewStruct( void )
970 {
971    v5dstruct *v;
972
973    v = (v5dstruct *) malloc( sizeof(v5dstruct) );
974    if (v) {
975       v5dInitStruct(v);
976    }
977    return v;
978 }
979
980
981
982 /*
983  * Free an initialized v5dstruct. (Todd Plessel)
984  */
985 void v5dFreeStruct( v5dstruct* v )
986 {
987    /*assert( v5dVerifyStruct( v ) );*/
988    free( v );
989    v = 0;
990 }
991
992
993
994 /*
995  * Do some checking that the information in a v5dstruct is valid.
996  * Input:  v - pointer to v5dstruct
997  * Return:  1 = g is ok, 0 = g is invalid
998  */
999 int v5dVerifyStruct( const v5dstruct *v )
1000 {
1001    int var, i, invalid, maxnl;
1002
1003    invalid = 0;
1004
1005    if (!v)
1006       return 0;
1007
1008    /* Number of variables */
1009    if (v->NumVars<0) {
1010       printf("Invalid number of variables: %d\n", v->NumVars );
1011       invalid = 1;
1012    }
1013    else if (v->NumVars>MAXVARS) {
1014       printf("Too many variables: %d  (Maximum is %d)\n",
1015              v->NumVars, MAXVARS);
1016       invalid = 1;
1017    }
1018
1019    /* Variable Names */
1020    for (i=0;i<v->NumVars;i++) {
1021       if (v->VarName[i][0]==0) {
1022          printf("Missing variable name: VarName[%d]=\"\"\n", i );
1023          invalid = 1;
1024       }
1025    }
1026
1027    /* Number of timesteps */
1028    if (v->NumTimes<0) {
1029       printf("Invalid number of timesteps: %d\n", v->NumTimes );
1030       invalid = 1;
1031    }
1032    else if (v->NumTimes>MAXTIMES) {
1033       printf("Too many timesteps: %d  (Maximum is %d)\n",
1034              v->NumTimes, MAXTIMES );
1035       invalid = 1;
1036    }
1037
1038    /* Make sure timestamps are increasing */
1039    for (i=1;i<v->NumTimes;i++) {
1040       int date0 = v5dYYDDDtoDays( v->DateStamp[i-1] );
1041       int date1 = v5dYYDDDtoDays( v->DateStamp[i] );
1042       int time0 = v5dHHMMSStoSeconds( v->TimeStamp[i-1] );
1043       int time1 = v5dHHMMSStoSeconds( v->TimeStamp[i] );
1044       if (time1<=time0 && date1<=date0) {
1045          printf("Timestamp for step %d must be later than step %d\n", i, i-1);
1046          invalid = 1;
1047       }
1048    }
1049
1050    /* Rows */
1051    if (v->Nr<2) {
1052       printf("Too few rows: %d (2 is minimum)\n", v->Nr );
1053       invalid = 1;
1054    }
1055    else if (v->Nr>MAXROWS) {
1056       printf("Too many rows: %d (%d is maximum)\n", v->Nr, MAXROWS );
1057       invalid = 1;
1058    }
1059
1060    /* Columns */
1061    if (v->Nc<2) {
1062       printf("Too few columns: %d (2 is minimum)\n", v->Nc );
1063       invalid = 1;
1064    }
1065    else if (v->Nc>MAXCOLUMNS) {
1066       printf("Too many columns: %d (%d is maximum)\n", v->Nc, MAXCOLUMNS );
1067       invalid = 1;
1068    }
1069
1070    /* Levels */
1071    maxnl = 0;
1072    for (var=0;var<v->NumVars;var++) {
1073       if (v->LowLev[var] < 0) {
1074          printf("Low level cannot be negative for var %s: %d\n",
1075                  v->VarName[var], v->LowLev[var] );
1076          invalid = 1;
1077       }
1078       if (v->Nl[var]<1) {
1079          printf("Too few levels for var %s: %d (1 is minimum)\n",
1080                  v->VarName[var], v->Nl[var] );
1081          invalid = 1;
1082       }
1083       if (v->Nl[var]+v->LowLev[var]>MAXLEVELS) {
1084          printf("Too many levels for var %s: %d (%d is maximum)\n",
1085                  v->VarName[var], v->Nl[var]+v->LowLev[var], MAXLEVELS );
1086          invalid = 1;
1087       }
1088       if (v->Nl[var]+v->LowLev[var]>maxnl) {
1089          maxnl = v->Nl[var]+v->LowLev[var];
1090       }
1091    }
1092
1093    if (v->CompressMode != 1 && v->CompressMode != 2 && v->CompressMode != 4) {
1094       printf("Bad CompressMode: %d (must be 1, 2 or 4)\n", v->CompressMode );
1095       invalid = 1;
1096    }
1097
1098    switch (v->VerticalSystem) {
1099       case 0:
1100       case 1:
1101          if (v->VertArgs[1]==0.0) {
1102             printf("Vertical level increment is zero, must be non-zero\n");
1103             invalid = 1;
1104          }
1105          break;
1106       case 2:
1107          /* Check that Height values increase upward */
1108          for (i=1;i<maxnl;i++) {
1109             if (v->VertArgs[i] <= v->VertArgs[i-1]) {
1110                printf("Height[%d]=%f <= Height[%d]=%f, level heights must increase\n",
1111                       i, v->VertArgs[i], i-1, v->VertArgs[i-1] );
1112                invalid = 1;
1113                break;
1114             }
1115          }
1116          break;
1117       case 3:
1118          /* Check that Pressure values decrease upward */
1119          for (i=1;i<maxnl;i++) {
1120             if (v->VertArgs[i] <= v->VertArgs[i-1]) {
1121                printf("Pressure[%d]=%f >= Pressure[%d]=%f, level pressures must decrease\n",
1122                       i, height_to_pressure(v->VertArgs[i]),
1123                       i-1, height_to_pressure(v->VertArgs[i-1]) );
1124                invalid = 1;
1125                break;
1126             }
1127          }
1128          break;
1129       default:
1130          printf("VerticalSystem = %d, must be in 0..3\n", v->VerticalSystem );
1131          invalid = 1;
1132    }
1133
1134
1135    switch (v->Projection) {
1136       case 0:  /* Generic */
1137          if (v->ProjArgs[2]==0.0) {
1138             printf("Row Increment (ProjArgs[2]) can't be zero\n");
1139             invalid = 1;
1140          }
1141          if (v->ProjArgs[3]==0.0) {
1142             printf("Column increment (ProjArgs[3]) can't be zero\n");
1143             invalid = 1;
1144          }
1145          break;
1146       case 1:  /* Cylindrical equidistant */
1147          if (v->ProjArgs[2]<0.0) {
1148             printf("Row Increment (ProjArgs[2]) = %g  (must be >=0.0)\n",
1149                    v->ProjArgs[2] );
1150             invalid = 1;
1151          }
1152          if (v->ProjArgs[3]<=0.0) {
1153             printf("Column Increment (ProjArgs[3]) = %g  (must be >=0.0)\n",
1154                    v->ProjArgs[3] );
1155             invalid = 1;
1156          }
1157          break;
1158       case 2:  /* Lambert Conformal */
1159          if (v->ProjArgs[0]<-90.0 || v->ProjArgs[0]>90.0) {
1160             printf("Lat1 (ProjArgs[0]) out of range: %g\n", v->ProjArgs[0] );
1161             invalid = 1;
1162          }
1163          if (v->ProjArgs[1]<-90.0 || v->ProjArgs[1]>90.0) {
1164             printf("Lat2 (ProjArgs[1] out of range: %g\n", v->ProjArgs[1] );
1165             invalid = 1;
1166          }
1167          if (v->ProjArgs[5]<=0.0) {
1168             printf("ColInc (ProjArgs[5]) = %g  (must be >=0.0)\n",
1169                    v->ProjArgs[5] );
1170             invalid = 1;
1171          }
1172          break;
1173       case 3:  /* Stereographic */
1174          if (v->ProjArgs[0]<-90.0 || v->ProjArgs[0]>90.0) {
1175             printf("Central Latitude (ProjArgs[0]) out of range: ");
1176             printf("%g  (must be in +/-90)\n", v->ProjArgs[0] );
1177             invalid = 1;
1178          }
1179          if (v->ProjArgs[1]<-180.0 || v->ProjArgs[1]>180.0) {
1180             printf("Central Longitude (ProjArgs[1]) out of range: ");
1181             printf("%g  (must be in +/-180)\n", v->ProjArgs[1] );
1182             invalid = 1;
1183          }
1184          if (v->ProjArgs[4]<0) {
1185             printf("Column spacing (ProjArgs[4]) = %g  (must be positive)\n",
1186                    v->ProjArgs[4]);
1187             invalid = 1;
1188          }
1189          break;
1190       case 4:  /* Rotated */
1191          /* WLH 4-21-95 */
1192          if (v->ProjArgs[2]<=0.0) {
1193             printf("Row Increment (ProjArgs[2]) = %g  (must be >=0.0)\n",
1194                    v->ProjArgs[2] );
1195             invalid = 1;
1196          }
1197          if (v->ProjArgs[3]<=0.0) {
1198             printf("Column Increment = (ProjArgs[3]) %g  (must be >=0.0)\n",
1199                    v->ProjArgs[3] );
1200             invalid = 1;
1201          }
1202          if (v->ProjArgs[4]<-90.0 || v->ProjArgs[4]>90.0) {
1203             printf("Central Latitude (ProjArgs[4]) out of range: ");
1204             printf("%g  (must be in +/-90)\n", v->ProjArgs[4] );
1205             invalid = 1;
1206          }
1207          if (v->ProjArgs[5]<-180.0 || v->ProjArgs[5]>180.0) {
1208             printf("Central Longitude (ProjArgs[5]) out of range: ");
1209             printf("%g  (must be in +/-180)\n", v->ProjArgs[5] );
1210             invalid = 1;
1211          }
1212          if (v->ProjArgs[6]<-180.0 || v->ProjArgs[6]>180.0) {
1213             printf("Central Longitude (ProjArgs[6]) out of range: ");
1214             printf("%g  (must be in +/-180)\n", v->ProjArgs[6] );
1215             invalid = 1;
1216          }
1217          break;
1218       default:
1219          printf("Projection = %d, must be in 0..4\n", v->Projection );
1220          invalid = 1;
1221    }
1222
1223    return !invalid;
1224 }
1225
1226
1227
1228 /*
1229  * Get the McIDAS file number and grid number associated with the grid
1230  * identified by time and var.
1231  * Input:  v - v5d grid struct
1232  *         time, var - timestep and variable of grid
1233  * Output:  mcfile, mcgrid - McIDAS grid file number and grid number
1234  */
1235 int v5dGetMcIDASgrid( v5dstruct *v, int time, int var,
1236                       int *mcfile, int *mcgrid )
1237 {
1238    if (time<0 || time>=v->NumTimes) {
1239       printf("Bad time argument to v5dGetMcIDASgrid: %d\n", time );
1240       return 0;
1241    }
1242    if (var<0 || var>=v->NumVars) {
1243       printf("Bad var argument to v5dGetMcIDASgrid: %d\n", var );
1244       return 0;
1245    }
1246
1247    *mcfile = (int) v->McFile[time][var];
1248    *mcgrid = (int) v->McGrid[time][var];
1249    return 1;
1250 }
1251
1252
1253
1254 /*
1255  * Set the McIDAS file number and grid number associated with the grid
1256  * identified by time and var.
1257  * Input:  v - v5d grid struct
1258  *         time, var - timestep and variable of grid
1259  *         mcfile, mcgrid - McIDAS grid file number and grid number
1260  * Return:  1 = ok, 0 = error (bad time or var)
1261  */
1262 int v5dSetMcIDASgrid( v5dstruct *v, int time, int var,
1263                       int mcfile, int mcgrid )
1264 {
1265    if (time<0 || time>=v->NumTimes) {
1266       printf("Bad time argument to v5dSetMcIDASgrid: %d\n", time );
1267       return 0;
1268    }
1269    if (var<0 || var>=v->NumVars) {
1270       printf("Bad var argument to v5dSetMcIDASgrid: %d\n", var );
1271       return 0;
1272    }
1273
1274    v->McFile[time][var] = (short) mcfile;
1275    v->McGrid[time][var] = (short) mcgrid;
1276    return 1;
1277 }
1278
1279
1280
1281 /**********************************************************************/
1282 /*****                    Input Functions                         *****/
1283 /**********************************************************************/
1284
1285
1286
1287 /*
1288  * Read the header from a COMP* file and return results in the v5dstruct.
1289  * Input:  f - the file descriptor
1290  *         v - pointer to a v5dstruct.
1291  * Return:  1 = ok, 0 = error.
1292  */
1293 static int read_comp_header( int f, v5dstruct *v )
1294 {
1295    unsigned int id;
1296
1297    /* reset file position to start of file */
1298    lseek( f, 0, SEEK_SET );
1299
1300    /* read file ID */
1301    read_int4( f, (int *) &id );
1302
1303    if (id==0x80808080 || id==0x80808081) {
1304       /* Older COMP5D format */
1305       int gridtimes, gridparms;
1306       int i, j, it, iv, nl;
1307       int gridsize;
1308       float hgttop, hgtinc;
1309       /*char *compgrid;*/
1310
1311       if (id==0x80808080) {
1312          /* 20 vars, 300 times */
1313          gridtimes = 300;
1314          gridparms = 20;
1315       }
1316       else {
1317          /* 30 vars, 400 times */
1318          gridtimes = 400;
1319          gridparms = 30;
1320       }
1321
1322       v->FirstGridPos = 12*4 + 8*gridtimes + 4*gridparms;
1323
1324       read_int4( f, &v->NumTimes );
1325       read_int4( f, &v->NumVars );
1326       read_int4( f, &v->Nr );
1327       read_int4( f, &v->Nc );
1328       read_int4( f, &nl );
1329       for (i=0;i<v->NumVars;i++) {
1330          v->Nl[i] = nl;
1331          v->LowLev[i] = 0;
1332       }
1333       read_float4( f, &v->ProjArgs[0] );
1334       read_float4( f, &v->ProjArgs[1] );
1335       read_float4( f, &hgttop );
1336       read_float4( f, &v->ProjArgs[2] );
1337       read_float4( f, &v->ProjArgs[3] );
1338       read_float4( f, &hgtinc );
1339 /*
1340       for (i=0;i<nl;i++) {
1341          v->Height[nl-i-1] = hgttop - i * hgtinc;
1342       }
1343 */
1344       v->VerticalSystem = 1;
1345       v->VertArgs[0] = hgttop - hgtinc * (nl-1);
1346       v->VertArgs[1] = hgtinc;
1347
1348       /* read dates and times */
1349       for (i=0;i<gridtimes;i++) {
1350          read_int4( f, &j );
1351          v->DateStamp[i] = v5dDaysToYYDDD( j );
1352       }
1353       for (i=0;i<gridtimes;i++) {
1354          read_int4( f, &j );
1355          v->TimeStamp[i] = v5dSecondsToHHMMSS( j );
1356       }
1357
1358       /* read variable names */
1359       for (i=0;i<gridparms;i++) {
1360          char name[4];
1361          read_bytes( f, name, 4 );
1362          /* remove trailing spaces, if any */
1363          for (j=3;j>0;j--) {
1364             if (name[j]==' ' || name[j]==0)
1365               name[j] = 0;
1366             else
1367               break;
1368          }
1369          strncpy( v->VarName[i], name, 4 );
1370          v->VarName[i][4] = 0;
1371       }
1372
1373       gridsize = ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
1374       for (i=0;i<v->NumVars;i++) {
1375          v->GridSize[i] = 8 + gridsize;
1376       }
1377       v->SumGridSizes = (8+gridsize) * v->NumVars;
1378
1379       /* read the grids and their ga,gb values to find min and max values */
1380
1381       for (i=0;i<v->NumVars;i++) {
1382          v->MinVal[i] = 999999.9;
1383          v->MaxVal[i] = -999999.9;
1384       }
1385
1386       /*compgrid = (char *) malloc( gridsize );*/
1387
1388       for (it=0; it<v->NumTimes; it++) {
1389          for (iv=0; iv<v->NumVars; iv++) {
1390             float ga, gb;
1391             float min, max;
1392
1393             read_float4( f, &ga );
1394             read_float4( f, &gb );
1395
1396             /* skip ahead by 'gridsize' bytes */
1397             if (lseek( f, gridsize, SEEK_CUR )==-1) {
1398                printf("Error:  Unexpected end of file, ");
1399                printf("file may be corrupted.\n");
1400                return 0;
1401             }
1402             min = -(125.0+gb)/ga;
1403             max = (125.0-gb)/ga;
1404             if (min<v->MinVal[iv])  v->MinVal[iv] = min;
1405             if (max>v->MaxVal[iv])  v->MaxVal[iv] = max;
1406          }
1407       }
1408
1409       /*free( compgrid );*/
1410
1411       /* done */
1412    }
1413    else if (id==0x80808082 || id==0x80808083) {
1414       /* Newer COMP5D format */
1415       int gridtimes, gridsize;
1416       int it, iv, nl, i, j;
1417       float delta;
1418
1419       read_int4( f, &gridtimes );
1420       read_int4( f, &v->NumVars );
1421       read_int4( f, &v->NumTimes );
1422       read_int4( f, &v->Nr );
1423       read_int4( f, &v->Nc );
1424       read_int4( f, &nl );
1425       for (i=0;i<v->NumVars;i++) {
1426          v->Nl[i] = nl;
1427       }
1428
1429       read_float4( f, &v->ProjArgs[2] );
1430       read_float4( f, &v->ProjArgs[3] );
1431
1432       /* Read height and determine if equal spacing */
1433       v->VerticalSystem = 1;
1434       for (i=0;i<nl;i++) {
1435          read_float4( f, &v->VertArgs[i] );
1436          if (i==1) {
1437             delta = v->VertArgs[1] - v->VertArgs[0];
1438          }
1439          else if (i>1) {
1440             if (delta != (v->VertArgs[i] - v->VertArgs[i-1])) {
1441                v->VerticalSystem = 2;
1442             }
1443          }
1444       }
1445       if (v->VerticalSystem==1) {
1446          v->VertArgs[1] = delta;
1447       }
1448
1449       /* read variable names */
1450       for (iv=0; iv<v->NumVars; iv++) {
1451          char name[8];
1452
1453          read_bytes( f, name, 8 );
1454
1455          /* remove trailing spaces, if any */
1456          for (j=7;j>0;j--) {
1457             if (name[j]==' ' || name[j]==0)
1458               name[j] = 0;
1459             else
1460               break;
1461          }
1462          strncpy( v->VarName[iv], name, 8 );
1463          v->VarName[iv][8] = 0;
1464       }
1465
1466       for (iv=0;iv<v->NumVars;iv++) {
1467          read_float4( f, &v->MinVal[iv] );
1468       }
1469       for (iv=0;iv<v->NumVars;iv++) {
1470          read_float4( f, &v->MaxVal[iv] );
1471       }
1472       for (it=0;it<gridtimes;it++) {
1473          read_int4( f, &j );
1474          v->TimeStamp[it] = v5dSecondsToHHMMSS( j );
1475       }
1476       for (it=0;it<gridtimes;it++) {
1477          read_int4( f, &j );
1478          v->DateStamp[it] = v5dDaysToYYDDD( j );
1479       }
1480       for (it=0;it<gridtimes;it++) {
1481          float nlat;
1482          read_float4( f, &nlat );
1483          if (it==0)  v->ProjArgs[0] = nlat;
1484       }
1485       for (it=0;it<gridtimes;it++) {
1486          float wlon;
1487          read_float4( f, &wlon );
1488          if (it==0)  v->ProjArgs[1] = wlon;
1489       }
1490
1491       /* calculate grid storage sizes */
1492       if (id==0x80808082) {
1493          gridsize = nl*2*4 + ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
1494       }
1495       else {
1496          /* McIDAS grid and file numbers present */
1497          gridsize = 8 + nl*2*4 + ( (v->Nr * v->Nc * nl + 3) / 4) * 4;
1498       }
1499       for (i=0;i<v->NumVars;i++) {
1500          v->GridSize[i] = gridsize;
1501       }
1502       v->SumGridSizes = gridsize * v->NumVars;
1503
1504       /* read McIDAS numbers??? */
1505
1506       /* size (in bytes) of all header info */
1507       v->FirstGridPos = 9*4 + v->Nl[0]*4 + v->NumVars*16 + gridtimes*16;
1508
1509    }
1510
1511    v->CompressMode = 1; /* one byte per grid point */
1512    v->Projection = 1;  /* Cylindrical equidistant */
1513    v->FileVersion[0] = 0;
1514
1515    return 1;
1516 }
1517
1518
1519
1520 /*
1521  * Read a compressed grid from a COMP* file.
1522  * Return:  1 = ok, 0 = error.
1523  */
1524 static int read_comp_grid( v5dstruct *v, int time, int var,
1525                            float *ga, float *gb, void *compdata )
1526 {
1527    unsigned int pos;
1528    V5Dubyte bias;
1529    int i, n, nl;
1530    int f;
1531    V5Dubyte *compdata1 = (V5Dubyte *) compdata;
1532
1533    f = v->FileDesc;
1534
1535    /* move to position in file */
1536    pos = grid_position( v, time, var );
1537    lseek( f, pos, SEEK_SET );
1538
1539    if (v->FileFormat==0x80808083) {
1540       /* read McIDAS grid and file numbers */
1541       int mcfile, mcgrid;
1542       read_int4( f, &mcfile );
1543       read_int4( f, &mcgrid );
1544       v->McFile[time][var] = (short) mcfile;
1545       v->McGrid[time][var] = (short) mcgrid;
1546    }
1547
1548    nl = v->Nl[var];
1549
1550    if (v->FileFormat==0x80808080 || v->FileFormat==0x80808081) {
1551       /* single ga,gb pair for whole grid */
1552       float a, b;
1553       read_float4( f, &a );
1554       read_float4( f, &b );
1555       /* convert a, b to new v5d ga, gb values */
1556       for (i=0;i<nl;i++) {
1557          if (a==0.0) {
1558             ga[i] = gb[i] = 0.0;
1559          }
1560          else {
1561             gb[i] = (b+128.0) / -a;
1562             ga[i] = 1.0 / a;
1563          }
1564       }
1565       bias = 128;
1566    }
1567    else {
1568       /* read ga, gb arrays */
1569       read_float4_array( f, ga, v->Nl[var] );
1570       read_float4_array( f, gb, v->Nl[var] );
1571
1572       /* convert ga, gb values to v5d system */
1573       for (i=0;i<nl;i++) {
1574          if (ga[i]==0.0) {
1575             ga[i] = gb[i] = 0.0;
1576          }
1577          else {
1578             /*gb[i] = (gb[i]+125.0) / -ga[i];*/
1579             gb[i] = (gb[i]+128.0) / -ga[i];
1580             ga[i] = 1.0 / ga[i];
1581          }
1582       }
1583       bias = 128;  /* 125 ??? */
1584    }
1585
1586    /* read compressed grid data */
1587    n = v->Nr * v->Nc * v->Nl[var];
1588    if (read_bytes( f, compdata1, n )!=n)
1589       return 0;
1590
1591    /* convert data values to v5d system */
1592    n = v->Nr * v->Nc * v->Nl[var];
1593    for (i=0;i<n;i++) {
1594       compdata1[i] += bias;
1595    }
1596
1597    return 1;
1598 }
1599
1600
1601
1602 /*
1603  * Read a v5d file header.
1604  * Input:  f - file opened for reading.
1605  *         v - pointer to v5dstruct to store header info into.
1606  * Return:  1 = ok, 0 = error.
1607  */
1608 static int read_v5d_header( v5dstruct *v )
1609 {
1610 #define SKIP(N)   lseek( f, N, SEEK_CUR )
1611    int end_of_header = 0;
1612    unsigned int id;
1613    int idlen, var, numargs;
1614    int f;
1615
1616    f = v->FileDesc;
1617
1618    /* first try to read the header id */
1619    read_int4( f, (int*) &id );
1620    read_int4( f, &idlen );
1621    if (id==TAG_ID && idlen==0) {
1622       /* this is a v5d file */
1623       v->FileFormat = 0;
1624    }
1625    else if (id>=0x80808080 && id<=0x80808083) {
1626       /* this is an old COMP* file */
1627       v->FileFormat = id;
1628       return read_comp_header( f, v );
1629    }
1630    else {
1631       /* unknown file type */
1632       printf("Error: not a v5d file\n");
1633       return 0;
1634    }
1635
1636    v->CompressMode = 1; /* default */
1637
1638    while (!end_of_header) {
1639       int tag, length;
1640       int i, var, time, nl, lev;
1641
1642       if (read_int4(f,&tag)<1 || read_int4(f,&length)<1) {
1643          printf("Error while reading header, premature EOF\n");
1644          return 0;
1645       }
1646
1647       switch (tag) {
1648          case TAG_VERSION:
1649             assert( length==10 );
1650             read_bytes( f, v->FileVersion, 10 );
1651             /* Check if reading a file made by a future version of Vis5D */
1652             if (strcmp(v->FileVersion, FILE_VERSION)>0) {
1653                /* WLH 6 Oct 98 */
1654                printf("Warning: Trying to read a version %s file,", v->FileVersion);
1655                printf(" you should upgrade Vis5D.\n");
1656             }
1657             break;
1658          case TAG_NUMTIMES:
1659             assert( length==4 );
1660             read_int4( f, &v->NumTimes );
1661             break;
1662          case TAG_NUMVARS:
1663             assert( length==4 );
1664             read_int4( f, &v->NumVars );
1665             break;
1666          case TAG_VARNAME:
1667             assert( length==14 );   /* 1 int + 10 char */
1668             read_int4( f, &var );
1669             read_bytes( f, v->VarName[var], 10 );
1670             break;
1671          case TAG_NR:
1672             /* Number of rows for all variables */
1673             assert( length==4 );
1674             read_int4( f, &v->Nr );
1675             break;
1676          case TAG_NC:
1677             /* Number of columns for all variables */
1678             assert( length==4 );
1679             read_int4( f, &v->Nc );
1680             break;
1681          case TAG_NL:
1682             /* Number of levels for all variables */
1683             assert( length==4 );
1684             read_int4( f, &nl );
1685             for (i=0;i<v->NumVars;i++) {
1686                v->Nl[i] = nl;
1687             }
1688             break;
1689          case TAG_NL_VAR:
1690             /* Number of levels for one variable */
1691             assert( length==8 );
1692             read_int4( f, &var );
1693             read_int4( f, &v->Nl[var] );
1694             break;
1695          case TAG_LOWLEV_VAR:
1696             /* Lowest level for one variable */
1697             assert( length==8 );
1698             read_int4( f, &var );
1699             read_int4( f, &v->LowLev[var] );
1700             break;
1701
1702          case TAG_TIME:
1703             /* Time stamp for 1 timestep */
1704             assert( length==8 );
1705             read_int4( f, &time );
1706             read_int4( f, &v->TimeStamp[time] );
1707             break;
1708          case TAG_DATE:
1709             /* Date stamp for 1 timestep */
1710             assert( length==8 );
1711             read_int4( f, &time );
1712             read_int4( f, &v->DateStamp[time] );
1713             break;
1714
1715          case TAG_MINVAL:
1716             /* Minimum value for a variable */
1717             assert( length==8 );
1718             read_int4( f, &var );
1719             read_float4( f, &v->MinVal[var] );
1720             break;
1721          case TAG_MAXVAL:
1722             /* Maximum value for a variable */
1723             assert( length==8 );
1724             read_int4( f, &var );
1725             read_float4( f, &v->MaxVal[var] );
1726             break;
1727          case TAG_COMPRESS:
1728             /* Compress mode */
1729             assert( length==4 );
1730             read_int4( f, &v->CompressMode );
1731             break;
1732          case TAG_UNITS:
1733             /* physical units */
1734             assert( length==24 );
1735             read_int4( f, &var );
1736             read_bytes( f, v->Units[var], 20 );
1737             break;
1738
1739          /*
1740           * Vertical coordinate system
1741           */
1742          case TAG_VERTICAL_SYSTEM:
1743             assert( length==4 );
1744             read_int4( f, &v->VerticalSystem );
1745             if (v->VerticalSystem<0 || v->VerticalSystem>3) {
1746                printf("Error: bad vertical coordinate system: %d\n",
1747                       v->VerticalSystem );
1748             }
1749             break;
1750          case TAG_VERT_ARGS:
1751             read_int4( f, &numargs );
1752             assert( numargs <= MAXVERTARGS );
1753             read_float4_array( f, v->VertArgs, numargs );
1754             assert( length==numargs*4+4 );
1755             break;
1756          case TAG_HEIGHT:
1757             /* height of a grid level */
1758             assert( length==8 );
1759             read_int4( f, &lev );
1760             read_float4( f, &v->VertArgs[lev] );
1761             break;
1762          case TAG_BOTTOMBOUND:
1763             assert( length==4 );
1764             read_float4( f, &v->VertArgs[0] );
1765             break;
1766          case TAG_LEVINC:
1767             assert( length==4 );
1768             read_float4( f, &v->VertArgs[1] );
1769             break;
1770
1771          /*
1772           * Map projection information
1773           */
1774          case TAG_PROJECTION:
1775             assert( length==4 );
1776             read_int4( f, &v->Projection );
1777             if (v->Projection<0 || v->Projection>4) { /* WLH 4-21-95 */
1778                printf("Error while reading header, bad projection (%d)\n",
1779                        v->Projection );
1780                return 0;
1781             }
1782             break;
1783          case TAG_PROJ_ARGS:
1784             read_int4( f, &numargs );
1785             assert( numargs <= MAXPROJARGS );
1786             read_float4_array( f, v->ProjArgs, numargs );
1787             assert( length==4*numargs+4 );
1788             break;
1789          case TAG_NORTHBOUND:
1790             assert( length==4 );
1791             if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
1792                read_float4( f, &v->ProjArgs[0] );
1793             }
1794             else {
1795                SKIP( 4 );
1796             }
1797             break;
1798          case TAG_WESTBOUND:
1799             assert( length==4 );
1800             if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
1801                read_float4( f, &v->ProjArgs[1] );
1802             }
1803             else {
1804                SKIP( 4 );
1805             }
1806             break;
1807          case TAG_ROWINC:
1808             assert( length==4 );
1809             if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
1810                read_float4( f, &v->ProjArgs[2] );
1811             }
1812             else {
1813                SKIP( 4 );
1814             }
1815             break;
1816          case TAG_COLINC:
1817             assert( length==4 );
1818             if (v->Projection==0 || v->Projection==1 || v->Projection==4) {
1819                read_float4( f, &v->ProjArgs[3] );
1820             }
1821             else if (v->Projection==2) {
1822                read_float4( f, &v->ProjArgs[5] );
1823             }
1824             else if (v->Projection==3) {
1825                read_float4( f, &v->ProjArgs[4] );
1826             }
1827             else {
1828                SKIP( 4 );
1829             }
1830             break;
1831          case TAG_LAT1:
1832             assert( length==4 );
1833             if (v->Projection==2) {
1834                read_float4( f, &v->ProjArgs[0] );
1835             }
1836             else {
1837                SKIP( 4 );
1838             }
1839             break;
1840          case TAG_LAT2:
1841             assert( length==4 );
1842             if (v->Projection==2) {
1843                read_float4( f, &v->ProjArgs[1] );
1844             }
1845             else {
1846                SKIP( 4 );
1847             }
1848             break;
1849          case TAG_POLE_ROW:
1850             assert( length==4 );
1851             if (v->Projection==2) {
1852                read_float4( f, &v->ProjArgs[2] );
1853             }
1854             else {
1855                SKIP( 4 );
1856             }
1857             break;
1858          case TAG_POLE_COL:
1859             assert( length==4 );
1860             if (v->Projection==2) {
1861                read_float4( f, &v->ProjArgs[3] );
1862             }
1863             else {
1864                SKIP( 4 );
1865             }
1866             break;
1867          case TAG_CENTLON:
1868             assert( length==4 );
1869             if (v->Projection==2) {
1870                read_float4( f, &v->ProjArgs[4] );
1871             }
1872             else if (v->Projection==3) {
1873                read_float4( f, &v->ProjArgs[1] );
1874             }
1875             else if (v->Projection==4) { /* WLH 4-21-95 */
1876                read_float4( f, &v->ProjArgs[5] );
1877             }
1878             else {
1879                SKIP( 4 );
1880             }
1881             break;
1882          case TAG_CENTLAT:
1883             assert( length==4 );
1884             if (v->Projection==3) {
1885                read_float4( f, &v->ProjArgs[0] );
1886             }
1887             else if (v->Projection==4) { /* WLH 4-21-95 */
1888                read_float4( f, &v->ProjArgs[4] );
1889             }
1890             else {
1891                SKIP( 4 );
1892             }
1893             break;
1894          case TAG_CENTROW:
1895             assert( length==4 );
1896             if (v->Projection==3) {
1897                read_float4( f, &v->ProjArgs[2] );
1898             }
1899             else {
1900                SKIP( 4 );
1901             }
1902             break;
1903          case TAG_CENTCOL:
1904             assert( length==4 );
1905             if (v->Projection==3) {
1906                read_float4( f, &v->ProjArgs[3] );
1907             }
1908             else {
1909                SKIP( 4 );
1910             }
1911             break;
1912          case TAG_ROTATION:
1913             assert( length==4 );
1914             if (v->Projection==4) { /* WLH 4-21-95 */
1915                read_float4( f, &v->ProjArgs[6] );
1916             }
1917             else {
1918                SKIP( 4 );
1919             }
1920             break;
1921
1922          case TAG_END:
1923             /* end of header */
1924             end_of_header = 1;
1925             lseek( f, length, SEEK_CUR );
1926             break;
1927
1928          default:
1929             /* unknown tag, skip to next tag */
1930             printf("Unknown tag: %d  length=%d\n", tag, length );
1931             lseek( f, length, SEEK_CUR );
1932             break;
1933       }
1934
1935    }
1936
1937    v5dVerifyStruct( v );
1938
1939    /* Now we're ready to read the grid data */
1940
1941    /* Save current file pointer */
1942    v->FirstGridPos = ltell(f);
1943
1944    /* compute grid sizes */
1945    v->SumGridSizes = 0;
1946    for (var=0;var<v->NumVars;var++) {
1947       v->GridSize[var] = 8 * v->Nl[var] + v5dSizeofGrid( v, 0, var );
1948       v->SumGridSizes += v->GridSize[var];
1949    }
1950
1951    return 1;
1952 #undef SKIP
1953 }
1954
1955
1956
1957
1958 /*
1959  * Open a v5d file for reading.
1960  * Input:  filename - name of v5d file to open
1961  *         v - pointer to a v5dstruct in which to put header info or NULL
1962  *             if a struct should be dynamically allocated.
1963  * Return:  NULL if error, else v or a pointer to a new v5dstruct if v was NULL
1964  */
1965 v5dstruct *v5dOpenFile( const char *filename, v5dstruct *v )
1966 {
1967    int fd;
1968
1969    fd = open( filename, O_RDONLY );
1970    if (fd==-1) {
1971       /* error */
1972       return 0;
1973    }
1974
1975    if (v) {
1976       v5dInitStruct( v );
1977    }
1978    else {
1979       v = v5dNewStruct();
1980       if (!v) {
1981          return NULL;
1982       }
1983    }
1984
1985    v->FileDesc = fd;
1986    v->Mode = 'r';
1987    if (read_v5d_header( v )) {
1988       return v;
1989    }
1990    else {
1991       return NULL;
1992    }
1993 }
1994
1995
1996
1997
1998 /*
1999  * Read a compressed grid from a v5d file.
2000  * Input:  v - pointer to v5dstruct describing the file
2001  *         time, var - which timestep and variable
2002  *         ga, gb - arrays to store grid (de)compression values
2003  *         compdata - address of where to store compressed grid data.
2004  * Return:  1 = ok, 0 = error.
2005  */
2006 int v5dReadCompressedGrid( v5dstruct *v, int time, int var,
2007                            float *ga, float *gb, void *compdata )
2008 {
2009    int pos, n, k;
2010
2011    if (time<0 || time>=v->NumTimes) {
2012       printf("Error in v5dReadCompressedGrid: bad timestep argument (%d)\n",
2013              time);
2014       return 0;
2015    }
2016    if (var<0 || var>=v->NumVars) {
2017       printf("Error in v5dReadCompressedGrid: bad var argument (%d)\n",
2018              var);
2019       return 0;
2020    }
2021
2022    if (v->FileFormat) {
2023       /* old COMP* file */
2024       return read_comp_grid( v, time, var, ga, gb, compdata );
2025    }
2026
2027    /* move to position in file */
2028    pos = grid_position( v, time, var );
2029    lseek( v->FileDesc, pos, SEEK_SET );
2030
2031    /* read ga, gb arrays */
2032    read_float4_array( v->FileDesc, ga, v->Nl[var] );
2033    read_float4_array( v->FileDesc, gb, v->Nl[var] );
2034
2035    /* read compressed grid data */
2036    n = v->Nr * v->Nc * v->Nl[var];
2037    if (v->CompressMode==1) {
2038       k = read_block( v->FileDesc, compdata, n, 1 )==n;
2039    }
2040    else if (v->CompressMode==2) {
2041       k = read_block( v->FileDesc, compdata, n, 2 )==n;
2042    }
2043    else if (v->CompressMode==4) {
2044       k = read_block( v->FileDesc, compdata, n, 4 )==n;
2045    }
2046    if (!k) {
2047       /* error */
2048       printf("Error in v5dReadCompressedGrid: read failed, bad file?\n");
2049    }
2050    return k;
2051
2052
2053 /*
2054    n = v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
2055    if (read( v->FileDesc, compdata, n )==n)
2056       return 1;
2057    else
2058       return 0;
2059 */
2060 }
2061
2062
2063
2064
2065 /*
2066  * Read a grid from a v5d file, decompress it and return it.
2067  * Input:  v - pointer to v5dstruct describing file header
2068  *         time, var - which timestep and variable.
2069  *         data - address of buffer to put grid data
2070  * Output:  data - the grid data
2071  * Return:  1 = ok, 0 = error.
2072  */
2073 int v5dReadGrid( v5dstruct *v, int time, int var, float data[] )
2074 {
2075    float ga[MAXLEVELS], gb[MAXLEVELS];
2076    void *compdata;
2077    int bytes;
2078
2079    if (time<0 || time>=v->NumTimes) {
2080       printf("Error in v5dReadGrid: bad timestep argument (%d)\n", time);
2081       return 0;
2082    }
2083    if (var<0 || var>=v->NumVars) {
2084       printf("Error in v5dReadGrid: bad variable argument (%d)\n", var);
2085       return 0;
2086    }
2087
2088    /* allocate compdata buffer */
2089    if (v->CompressMode==1) {
2090       bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned char);
2091    }
2092    else if (v->CompressMode==2) {
2093       bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned short);
2094    }
2095    else if (v->CompressMode==4) {
2096       bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(float);
2097    }
2098    compdata = (void *) malloc( bytes );
2099    if (!compdata) {
2100       printf("Error in v5dReadGrid: out of memory (needed %d bytes)\n", bytes);
2101       return 0;
2102    }
2103
2104    /* read the compressed data */
2105    if (!v5dReadCompressedGrid( v, time, var, ga, gb, compdata )) {
2106       return 0;
2107    }
2108
2109    /* decompress the data */
2110    v5dDecompressGrid( v->Nr, v->Nc, v->Nl[var], v->CompressMode,
2111                       compdata, ga, gb, data );
2112
2113    /* free compdata */
2114    free( compdata );
2115    return 1;
2116 }
2117
2118
2119
2120
2121 /**********************************************************************/
2122 /*****                   Output Functions                         *****/
2123 /**********************************************************************/
2124
2125
2126
2127 static int write_tag( v5dstruct *v, int tag, int length, int newfile )
2128 {
2129    if (!newfile) {
2130       /* have to check that there's room in header to write this tagged item */
2131       if (v->CurPos+8+length > v->FirstGridPos) {
2132          printf("Error: out of header space!\n");
2133          /* Out of header space! */
2134          return 0;
2135       }
2136    }
2137
2138    if (write_int4( v->FileDesc, tag )==0)  return 0;
2139    if (write_int4( v->FileDesc, length )==0)  return 0;
2140    v->CurPos += 8 + length;
2141    return 1;
2142 }
2143
2144
2145
2146 /*
2147  * Write the information in the given v5dstruct as a v5d file header.
2148  * Note that the current file position is restored when this function
2149  * returns normally.
2150  * Input:  f - file already open for writing
2151  *         v - pointer to v5dstruct
2152  * Return:  1 = ok, 0 = error.
2153  */
2154 static int write_v5d_header( v5dstruct *v )
2155 {
2156    int var, time, filler, maxnl;
2157    int f;
2158    int newfile;
2159
2160    if (v->FileFormat!=0) {
2161       printf("Error: v5d library can't write comp5d format files.\n");
2162       return 0;
2163    }
2164
2165    f = v->FileDesc;
2166
2167    if (!v5dVerifyStruct( v ))
2168       return 0;
2169
2170    /* Determine if we're writing to a new file */
2171    if (v->FirstGridPos==0) {
2172       newfile = 1;
2173    }
2174    else {
2175       newfile = 0;
2176    }
2177
2178    /* compute grid sizes */
2179    v->SumGridSizes = 0;
2180    for (var=0;var<v->NumVars;var++) {
2181       v->GridSize[var] = 8 * v->Nl[var] + v5dSizeofGrid( v, 0, var );
2182       v->SumGridSizes += v->GridSize[var];
2183    }
2184
2185    /* set file pointer to start of file */
2186    lseek( f, 0, SEEK_SET );
2187    v->CurPos = 0;
2188
2189    /*
2190     * Write the tagged header info
2191     */
2192 #define WRITE_TAG( V, T, L )  if (!write_tag(V,T,L,newfile))  return 0;
2193
2194    /* ID */
2195    WRITE_TAG( v, TAG_ID, 0 );
2196
2197    /* File Version */
2198    WRITE_TAG( v, TAG_VERSION, 10 );
2199    write_bytes( f, FILE_VERSION, 10 );
2200
2201    /* Number of timesteps */
2202    WRITE_TAG( v, TAG_NUMTIMES, 4 );
2203    write_int4( f, v->NumTimes );
2204
2205    /* Number of variables */
2206    WRITE_TAG( v, TAG_NUMVARS, 4 );
2207    write_int4( f, v->NumVars );
2208
2209    /* Names of variables */
2210    for (var=0;var<v->NumVars;var++) {
2211       WRITE_TAG( v, TAG_VARNAME, 14 );
2212       write_int4( f, var );
2213       write_bytes( f, v->VarName[var], 10 );
2214    }
2215
2216    /* Physical Units */
2217    for (var=0;var<v->NumVars;var++) {
2218       WRITE_TAG( v, TAG_UNITS, 24 );
2219       write_int4( f, var );
2220       write_bytes( f, v->Units[var], 20 );
2221    }
2222
2223    /* Date and time of each timestep */
2224    for (time=0;time<v->NumTimes;time++) {
2225       WRITE_TAG( v, TAG_TIME, 8 );
2226       write_int4( f, time );
2227       write_int4( f, v->TimeStamp[time] );
2228       WRITE_TAG( v, TAG_DATE, 8 );
2229       write_int4( f, time );
2230       write_int4( f, v->DateStamp[time] );
2231    }
2232
2233    /* Number of rows */
2234    WRITE_TAG( v, TAG_NR, 4 );
2235    write_int4( f, v->Nr );
2236
2237    /* Number of columns */
2238    WRITE_TAG( v, TAG_NC, 4 );
2239    write_int4( f, v->Nc );
2240
2241    /* Number of levels, compute maxnl */
2242    maxnl = 0;
2243    for (var=0;var<v->NumVars;var++) {
2244       WRITE_TAG( v, TAG_NL_VAR, 8 );
2245       write_int4( f, var );
2246       write_int4( f, v->Nl[var] );
2247       WRITE_TAG( v, TAG_LOWLEV_VAR, 8 );
2248       write_int4( f, var );
2249       write_int4( f, v->LowLev[var] );
2250       if (v->Nl[var]+v->LowLev[var]>maxnl) {
2251          maxnl = v->Nl[var]+v->LowLev[var];
2252       }
2253    }
2254
2255    /* Min/Max values */
2256    for (var=0;var<v->NumVars;var++) {
2257       WRITE_TAG( v, TAG_MINVAL, 8 );
2258       write_int4( f, var );
2259       write_float4( f, v->MinVal[var] );
2260       WRITE_TAG( v, TAG_MAXVAL, 8 );
2261       write_int4( f, var );
2262       write_float4( f, v->MaxVal[var] );
2263    }
2264
2265    /* Compress mode */
2266    WRITE_TAG( v, TAG_COMPRESS, 4 );
2267    write_int4( f, v->CompressMode );
2268
2269    /* Vertical Coordinate System */
2270    WRITE_TAG( v, TAG_VERTICAL_SYSTEM, 4 );
2271    write_int4( f, v->VerticalSystem );
2272    WRITE_TAG( v, TAG_VERT_ARGS, 4+4*MAXVERTARGS );
2273    write_int4( f, MAXVERTARGS );
2274    write_float4_array( f, v->VertArgs, MAXVERTARGS );
2275
2276    /* Map Projection */
2277    WRITE_TAG( v, TAG_PROJECTION, 4 );
2278    write_int4( f, v->Projection );
2279    WRITE_TAG( v, TAG_PROJ_ARGS, 4+4*MAXPROJARGS );
2280    write_int4( f, MAXPROJARGS );
2281    write_float4_array( f, v->ProjArgs, MAXPROJARGS );
2282
2283    /* write END tag */
2284    if (newfile) {
2285       /* We're writing to a brand new file.  Reserve 10000 bytes */
2286       /* for future header growth. */
2287       WRITE_TAG( v, TAG_END, 10000 );
2288       lseek( f, 10000, SEEK_CUR );
2289
2290       /* Let file pointer indicate where first grid is stored */
2291       v->FirstGridPos = ltell( f );
2292    }
2293    else {
2294       /* we're rewriting a header */
2295       filler = v->FirstGridPos - ltell(f);
2296       WRITE_TAG( v, TAG_END, filler-8 );
2297    }
2298
2299 #undef WRITE_TAG
2300
2301    return 1;
2302 }
2303
2304
2305
2306 /*
2307  * Open a v5d file for writing.  If the named file already exists,
2308  * it will be deleted.
2309  * Input:  filename - name of v5d file to create.
2310  *         v - pointer to v5dstruct with the header info to write.
2311  * Return:  1 = ok, 0 = error.
2312  */
2313 int v5dCreateFile( const char *filename, v5dstruct *v )
2314 {
2315    mode_t mask;
2316    int fd;
2317
2318    mask = 0666;
2319    fd = open( filename, O_WRONLY | O_CREAT | O_TRUNC, mask );
2320    if (fd==-1) {
2321       printf("Error in v5dCreateFile: open failed\n");
2322       v->FileDesc = -1;
2323       v->Mode = 0;
2324       return 0;
2325    }
2326    else {
2327       /* ok */
2328       v->FileDesc = fd;
2329       v->Mode = 'w';
2330       /* write header and return status */
2331       return write_v5d_header(v);
2332    }
2333 }
2334
2335
2336
2337 /*
2338  * Open a v5d file for updating/appending and read the header info.
2339  * Input:  filename - name of v5d file to open for updating.
2340  *         v - pointer to v5dstruct in which the file header info will be
2341  *             put.  If v is NULL a v5dstruct will be allocated and returned.
2342  * Return:  NULL if error, else v or a pointer to a new v5dstruct if v as NULL
2343  */
2344 v5dstruct *v5dUpdateFile( const char *filename, v5dstruct *v )
2345 {
2346    int fd;
2347
2348    fd = open( filename, O_RDWR );
2349    if (fd==-1) {
2350       return NULL;
2351    }
2352
2353    if (!v) {
2354       v = v5dNewStruct();
2355       if (!v) {
2356          return NULL;
2357       }
2358    }
2359
2360    v->FileDesc = fd;
2361    v->Mode = 'w';
2362
2363    if (read_v5d_header( v )) {
2364       return v;
2365    }
2366    else {
2367       return NULL;
2368    }
2369 }
2370
2371
2372
2373 /*
2374  * Write a compressed grid to a v5d file.
2375  * Input:  v - pointer to v5dstruct describing the file
2376  *         time, var - which timestep and variable
2377  *         ga, gb - the GA and GB (de)compression value arrays
2378  *         compdata - address of array of compressed data values
2379  * Return:  1 = ok, 0 = error.
2380  */
2381 int v5dWriteCompressedGrid( const v5dstruct *v, int time, int var,
2382                             const float *ga, const float *gb,
2383                             const void *compdata )
2384 {
2385    int pos, n, k;
2386
2387    /* simple error checks */
2388    if (v->Mode!='w') {
2389       printf("Error in v5dWriteCompressedGrid: file opened for reading,");
2390       printf(" not writing.\n");
2391       return 0;
2392    }
2393    if (time<0 || time>=v->NumTimes) {
2394       printf("Error in v5dWriteCompressedGrid: bad timestep argument (%d)\n",
2395              time);
2396       return 0;
2397    }
2398    if (var<0 || var>=v->NumVars) {
2399       printf("Error in v5dWriteCompressedGrid: bad variable argument (%d)\n",
2400              var);
2401       return 0;
2402    }
2403
2404    /* move to position in file */
2405    pos = grid_position( v, time, var );
2406    if (lseek( v->FileDesc, pos, SEEK_SET )<0) {
2407       /* lseek failed, return error */
2408       printf("Error in v5dWrite[Compressed]Grid: seek failed, disk full?\n");
2409       return 0;
2410    }
2411
2412    /* write ga, gb arrays */
2413    k = 0;
2414    if (write_float4_array( v->FileDesc, ga, v->Nl[var] ) == v->Nl[var] &&
2415        write_float4_array( v->FileDesc, gb, v->Nl[var] ) == v->Nl[var]) {
2416       /* write compressed grid data (k=1=OK, k=0=Error) */
2417       n = v->Nr * v->Nc * v->Nl[var];
2418       if (v->CompressMode==1) {
2419          k = write_block( v->FileDesc, compdata, n, 1 )==n;
2420       }
2421       else if (v->CompressMode==2) {
2422          k = write_block( v->FileDesc, compdata, n, 2 )==n;
2423       }
2424       else if (v->CompressMode==4) {
2425          k = write_block( v->FileDesc, compdata, n, 4 )==n;
2426       }
2427    }
2428
2429    if (k==0) {
2430       /* Error while writing */
2431       printf("Error in v5dWrite[Compressed]Grid: write failed, disk full?\n");
2432    }
2433    return k;
2434
2435 /*
2436    n = v->Nr * v->Nc * v->Nl[var] * v->CompressMode;
2437    if (write_bytes( v->FileDesc, compdata, n )!=n) {
2438       printf("Error in v5dWrite[Compressed]Grid: write failed, disk full?\n");
2439       return 0;
2440    }
2441    else {
2442       return 1;
2443    }
2444 */
2445 }
2446
2447
2448
2449
2450 /*
2451  * Compress a grid and write it to a v5d file.
2452  * Input:  v - pointer to v5dstruct describing the file
2453  *         time, var - which timestep and variable (starting at 0)
2454  *         data - address of uncompressed grid data
2455  * Return:  1 = ok, 0 = error.
2456  */
2457 int v5dWriteGrid( v5dstruct *v, int time, int var, const float data[] )
2458 {
2459    float ga[MAXLEVELS], gb[MAXLEVELS];
2460    void *compdata;
2461    int n, bytes;
2462    float min, max;
2463
2464    if (v->Mode!='w') {
2465       printf("Error in v5dWriteGrid: file opened for reading,");
2466       printf(" not writing.\n");
2467       return 0;
2468    }
2469    if (time<0 || time>=v->NumTimes) {
2470       printf("Error in v5dWriteGrid: bad timestep argument (%d)\n", time);
2471       return 0;
2472    }
2473    if (var<0 || var>=v->NumVars) {
2474       printf("Error in v5dWriteGrid: bad variable argument (%d)\n", var);
2475       return 0;
2476    }
2477
2478    /* allocate compdata buffer */
2479    if (v->CompressMode==1) {
2480       bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned char);
2481    }
2482    else if (v->CompressMode==2) {
2483       bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(unsigned short);
2484    }
2485    else if (v->CompressMode==4) {
2486       bytes = v->Nr * v->Nc * v->Nl[var] * sizeof(float);
2487    }
2488    compdata = (void *) malloc( bytes );
2489    if (!compdata) {
2490       printf("Error in v5dWriteGrid: out of memory (needed %d bytes)\n",
2491              bytes );
2492       return 0;
2493    }
2494
2495    /* compress the grid data */
2496    v5dCompressGrid( v->Nr, v->Nc, v->Nl[var], v->CompressMode, data,
2497                     compdata, ga, gb, &min, &max );
2498
2499    /* update min and max value */
2500    if (min<v->MinVal[var])
2501       v->MinVal[var] = min;
2502    if (max>v->MaxVal[var])
2503       v->MaxVal[var] = max;
2504
2505    /* write the compressed grid */
2506    n = v5dWriteCompressedGrid( v, time, var, ga, gb, compdata );
2507
2508    /* free compdata */
2509    free( compdata );
2510
2511    return n;
2512 }
2513
2514
2515
2516 /*
2517  * Close a v5d file which was opened with open_v5d_file() or
2518  * create_v5d_file().
2519  * Input: f - file descriptor
2520  * Return:  1 = ok, 0 = error
2521  */
2522 int v5dCloseFile( v5dstruct *v )
2523 {
2524    int status = 1;
2525
2526    if (v->Mode=='w') {
2527       /* rewrite header because writing grids updates the minval and */
2528       /* maxval fields */
2529       lseek( v->FileDesc, 0, SEEK_SET );
2530       status = write_v5d_header( v );
2531       lseek( v->FileDesc, 0, SEEK_END );
2532       close( v->FileDesc );
2533    }
2534    else if (v->Mode=='r') {
2535       /* just close the file */
2536       close(v->FileDesc);
2537    }
2538    else {
2539       printf("Error in v5dCloseFile: bad v5dstruct argument\n");
2540       return 0;
2541    }
2542    v->FileDesc = -1;
2543    v->Mode = 0;
2544    return status;
2545 }
2546
2547
2548
2549
2550 /**********************************************************************/
2551 /*****           Simple v5d file writing functions.               *****/
2552 /**********************************************************************/
2553
2554
2555
2556 static v5dstruct *Simple = NULL;
2557
2558
2559
2560 /*
2561  * Create a new v5d file specifying both a map projection and vertical
2562  * coordinate system.  See README file for argument details.
2563  * Return:  1 = ok, 0 = error.
2564  */
2565 int v5dCreate( const char *name, int numtimes, int numvars,
2566                int nr, int nc, const int nl[],
2567                const char varname[MAXVARS][10],
2568                const int timestamp[], const int datestamp[],
2569                int compressmode,
2570                int projection,
2571                const FLOAT proj_args[],
2572                int vertical,
2573                const FLOAT vert_args[] )
2574 {
2575    int var, time, maxnl, i;
2576
2577    /* initialize the v5dstruct */
2578    Simple = v5dNewStruct();
2579
2580    Simple->NumTimes = numtimes;
2581    Simple->NumVars = numvars;
2582    Simple->Nr = nr;
2583    Simple->Nc = nc;
2584    maxnl = nl[0];
2585    for (var=0;var<numvars;var++) {
2586       if (nl[var]>maxnl) {
2587          maxnl = nl[var];
2588       }
2589       Simple->Nl[var] = nl[var];
2590       Simple->LowLev[var] = 0;
2591       strncpy( Simple->VarName[var], varname[var], 10 );
2592       Simple->VarName[var][9] = 0;
2593    }
2594
2595    /* time and date for each timestep */
2596    for (time=0;time<numtimes;time++) {
2597       Simple->TimeStamp[time] = timestamp[time];
2598       Simple->DateStamp[time] = datestamp[time];
2599    }
2600
2601    Simple->CompressMode = compressmode;
2602
2603    /* Map projection and vertical coordinate system */
2604    Simple->Projection = projection;
2605 #ifdef VPP
2606    { 
2607        int i;
2608        for (i=0;i<MAXPROJARGS;i++)
2609          Simple->ProjArgs[i] =  (float)proj_args[i];
2610    }
2611 #else
2612    memcpy( Simple->ProjArgs, proj_args, MAXPROJARGS*sizeof(float) );
2613 #endif
2614    Simple->VerticalSystem = vertical;
2615    if (vertical == 3) {
2616      /* convert pressures to heights */
2617      for (i=0; i<MAXVERTARGS; i++) {
2618        if (vert_args[i] > 0.000001) {
2619          Simple->VertArgs[i] = pressure_to_height((float)vert_args[i]);
2620        }
2621        else Simple->VertArgs[i] = 0.0;
2622      }
2623    }
2624    else {
2625 #ifdef VPP
2626    { 
2627        int i;
2628        for (i=0;i<MAXVERTARGS;i++)
2629            Simple->VertArgs[i] =  (float)vert_args[i];
2630    }    
2631 #else
2632      memcpy( Simple->VertArgs, vert_args, MAXVERTARGS*sizeof(float) );
2633 #endif
2634    }
2635
2636    /* create the file */
2637    if (v5dCreateFile( name, Simple )==0) {
2638      printf("Error in v5dCreateSimpleFile: unable to create %s\n", name );
2639      return 0;
2640    }
2641    else {
2642       return 1;
2643    }
2644 }
2645
2646
2647
2648 /*
2649  * Create a new v5d file using minimal information.
2650  * Return:  1 = ok, 0 = error.  See README file for argument details.
2651  */
2652 int v5dCreateSimple( const char *name, int numtimes, int numvars,
2653                      int nr, int nc, int nl,
2654                      const char varname[MAXVARS][10],
2655                      const int timestamp[], const int datestamp[],
2656                      float northlat, float latinc,
2657                      float westlon, float loninc,
2658                      float bottomhgt, float hgtinc )
2659 {
2660    int nlvar[MAXVARS];
2661    int compressmode, projection, vertical;
2662    FLOAT proj_args[100], vert_args[MAXLEVELS];
2663    int i;
2664
2665    for (i=0;i<numvars;i++) {
2666       nlvar[i] = nl;
2667    }
2668
2669    compressmode = 1;
2670
2671    projection = 1;
2672    proj_args[0] = northlat;
2673    proj_args[1] = westlon;
2674    proj_args[2] = latinc;
2675    proj_args[3] = loninc;
2676
2677    vertical = 1;
2678    vert_args[0] = bottomhgt;
2679    vert_args[1] = hgtinc;
2680
2681    return v5dCreate( name, numtimes, numvars, nr, nc, nlvar,
2682                      varname, timestamp, datestamp, compressmode,
2683                      projection, proj_args, vertical, vert_args );
2684 }
2685
2686
2687
2688 /*
2689  * Set lowest levels for each variable (other than default of 0).
2690  * Input: lowlev - array [NumVars] of ints
2691  * Return:  1 = ok, 0 = error
2692  */
2693 int v5dSetLowLev( int lowlev[] )
2694 {
2695   int var;
2696
2697   if (Simple) {
2698      for (var=0;var<Simple->NumVars;var++) {
2699         Simple->LowLev[var] = lowlev[var];
2700      }
2701      return 1;
2702   }
2703   else {
2704      printf("Error: must call v5dCreate before v5dSetLowLev\n");
2705      return 0;
2706   }
2707 }
2708
2709
2710 /*
2711  * Set the units for a variable.
2712  * Input:  var - a variable in [1,NumVars]
2713  *         units - a string
2714  * Return:  1 = ok, 0 = error
2715  */
2716 int v5dSetUnits( int var, const char *units )
2717 {
2718   if (Simple) {
2719      if (var>=1 && var<=Simple->NumVars) {
2720         strncpy( Simple->Units[var-1], units, 19 );
2721         Simple->Units[var-1][19] = 0;
2722         return 1;
2723      }
2724      else {
2725         printf("Error: bad variable number in v5dSetUnits\n");
2726         return 0;
2727      }
2728   }
2729   else {
2730      printf("Error: must call v5dCreate before v5dSetUnits\n");
2731      return 0;
2732   }
2733 }
2734
2735
2736
2737 /*
2738  * Write a grid to a v5d file.
2739  * Input:  time - timestep in [1,NumTimes]
2740  *         var - timestep in [1,NumVars]
2741  *         data - array [nr*nc*nl] of floats
2742  * Return:  1 = ok, 0 = error
2743  */
2744 int v5dWrite( int time, int var, const FLOAT data[] )
2745 {
2746    if (Simple) {
2747       if (time<1 || time>Simple->NumTimes) {
2748          printf("Error in v5dWrite: bad timestep number: %d\n", time );
2749          return 0;
2750       }
2751       if (var<1 || var>Simple->NumVars) {
2752          printf("Error in v5dWrite: bad variable number: %d\n", var );
2753       }
2754 #ifdef VPP
2755       {
2756           float *rdata;
2757           int i,irep;
2758           int size = Simple->Nr * Simple->Nc * Simple->Nl[var-1]; 
2759           rdata = (float *)malloc(size * 4);
2760           if (!rdata){
2761               printf("Error in v5dWrite: out of memory\n");
2762               return 0;
2763           }
2764           for (i=0;i<size;i++)
2765               rdata[i] = (float)data[i];
2766           irep = v5dWriteGrid( Simple, time-1, var-1, rdata );
2767           free(rdata);
2768           return irep;
2769
2770       }
2771 #else      
2772       return v5dWriteGrid( Simple, time-1, var-1, data );
2773 #endif
2774    }
2775    else {
2776       printf("Error: must call v5dCreate before v5dWrite\n");
2777       return 0;
2778    }
2779 }
2780
2781
2782
2783 /*
2784  * Close a v5d file after the last grid has been written to it.
2785  * Return:  1 = ok, 0 = error
2786  */
2787 int v5dClose( void )
2788 {
2789    if (Simple) {
2790      int ok = v5dCloseFile( Simple );
2791      v5dFreeStruct( Simple );
2792      return ok;
2793    }
2794    else {
2795      printf("Error: v5dClose: no file to close\n");
2796      return 0;
2797    }
2798 }
2799
2800
2801
2802 /**********************************************************************/
2803 /*****                FORTRAN-callable simple output              *****/
2804 /**********************************************************************/
2805
2806
2807 /*
2808  * Create a v5d file.  See README file for argument descriptions.
2809  * Return:  1 = ok, 0 = error.
2810  */
2811 #ifdef UNDERSCORE
2812    int v5dcreate_
2813 #else
2814 #  ifdef _CRAY
2815      int V5DCREATE
2816 #  else
2817      int v5dcreate
2818 #  endif
2819 #endif
2820            ( const char *name, const int *numtimes, const int *numvars,
2821              const int *nr, const int *nc, const int nl[],
2822              const char varname[][10],
2823              const int timestamp[], const int datestamp[],
2824              const int *compressmode,
2825              const int *projection,
2826              const FLOAT proj_args[],
2827              const int *vertical,
2828              const FLOAT vert_args[] )
2829 {
2830    char filename[100];
2831    char names[MAXVARS][10];
2832    int i, maxnl, args;
2833
2834    /* copy name to filename and remove trailing spaces if any */
2835    copy_string( filename, name, 100 );
2836
2837    /*
2838     * Check for uninitialized arguments
2839     */
2840    if (*numtimes<1) {
2841       printf("Error: numtimes invalid\n");
2842       return 0;
2843    }
2844    if (*numvars<1) {
2845       printf("Error: numvars invalid\n");
2846       return 0;
2847    }
2848    if (*nr<2) {
2849       printf("Error: nr invalid\n");
2850       return 0;
2851    }
2852    if (*nc<2) {
2853       printf("Error: nc invalid\n");
2854       return 0;
2855    }
2856    maxnl = 0;
2857    for (i=0;i<*numvars;i++) {
2858       if (nl[i]<1) {
2859          printf("Error: nl(%d) invalid\n", i+1);
2860          return 0;
2861       }
2862       if (nl[i]>maxnl) {
2863          maxnl = nl[i];
2864       }
2865    }
2866
2867    for (i=0;i<*numvars;i++) {
2868       if (copy_string2( names[i], varname[i], 10)==0) {
2869          printf("Error: unitialized varname(%d)\n", i+1);
2870          return 0;
2871       }
2872    }
2873
2874    for (i=0;i<*numtimes;i++) {
2875       if (timestamp[i]<0) {
2876          printf("Error: times(%d) invalid\n", i+1);
2877          return 0;
2878       }
2879       if (datestamp[i]<0) {
2880          printf("Error: dates(%d) invalid\n", i+1);
2881          return 0;
2882       }
2883    }
2884
2885    if (*compressmode != 1 && *compressmode != 2 && *compressmode != 4) {
2886       printf("Error: compressmode invalid\n");
2887       return 0;
2888    }
2889
2890    switch (*projection) {
2891       case 0:
2892          args = 4;
2893          break;
2894       case 1:
2895          args = 0;
2896          if (IS_MISSING(proj_args[0])) {
2897             printf("Error: northlat (proj_args(1)) invalid\n");
2898             return 0;
2899          }
2900          if (IS_MISSING(proj_args[1])) {
2901             printf("Error: westlon (proj_args(2)) invalid\n");
2902             return 0;
2903          }
2904          if (IS_MISSING(proj_args[2])) {
2905             printf("Error: latinc (proj_args(3)) invalid\n");
2906             return 0;
2907          }
2908          if (IS_MISSING(proj_args[3])) {
2909             printf("Error: loninc (proj_args(4)) invalid\n");
2910             return 0;
2911          }
2912          break;
2913       case 2:
2914          args = 6;
2915          break;
2916       case 3:
2917          args = 5;
2918          break;
2919       case 4:
2920          args = 7;
2921          break;
2922       default:
2923          args = 0;
2924          printf("Error: projection invalid\n");
2925          return 0;
2926    }
2927    for (i=0;i<args;i++) {
2928       if (IS_MISSING(proj_args[i])) {
2929          printf("Error: proj_args(%d) invalid\n", i+1);
2930          return 0;
2931       }
2932    }
2933
2934    switch (*vertical) {
2935       case 0:
2936 /* WLH 31 Oct 96  -  just fall through 
2937          args = 4;
2938          break;
2939 */
2940       case 1:
2941          args = 0;
2942          if (IS_MISSING(vert_args[0])) {
2943             printf("Error: bottomhgt (vert_args(1)) invalid\n");
2944             return 0;
2945          }
2946          if (IS_MISSING(vert_args[1])) {
2947             printf("Error: hgtinc (vert_args(2)) invalid\n");
2948             return 0;
2949          }
2950          break;
2951       case 2:
2952       case 3:
2953          args = maxnl;
2954          break;
2955       default:
2956          args = 0;
2957          printf("Error: vertical invalid\n");
2958          return 0;
2959    }
2960    for (i=0;i<args;i++) {
2961       if (IS_MISSING(vert_args[i])) {
2962          printf("Error: vert_args(%d) invalid\n", i+1);
2963          return 0;
2964       }
2965    }
2966
2967    return v5dCreate( filename, *numtimes, *numvars, *nr, *nc, nl,
2968                      (const char(*)[10]) names, timestamp, datestamp,
2969                      *compressmode,
2970                      *projection, proj_args, *vertical, vert_args );
2971 }
2972
2973
2974
2975
2976 /*
2977  * Create a simple v5d file.  See README file for argument descriptions.
2978  * Return:  1 = ok, 0 = error.
2979  */
2980 #ifdef UNDERSCORE
2981   int v5dcreatesimple_
2982 #else
2983 #  ifdef _CRAY
2984      int V5DCREATESIMPLE
2985 #  else
2986      int v5dcreatesimple
2987 #  endif
2988 #endif
2989            ( const char *name, const int *numtimes, const int *numvars,
2990              const int *nr, const int *nc, const int *nl,
2991              const char varname[][10],
2992              const int timestamp[], const int datestamp[],
2993              const float *northlat, const float *latinc,
2994              const float *westlon, const float *loninc,
2995              const float *bottomhgt, const float *hgtinc )
2996 {
2997    int compressmode, projection, vertical;
2998    FLOAT projarg[100], vertarg[MAXLEVELS];
2999    int varnl[MAXVARS];
3000    int i;
3001
3002    for (i=0;i<MAXVARS;i++) {
3003       varnl[i] = *nl;
3004    }
3005
3006    compressmode = 1;
3007
3008    projection = 1;
3009    projarg[0] = *northlat;
3010    projarg[1] = *westlon;
3011    projarg[2] = *latinc;
3012    projarg[3] = *loninc;
3013
3014    vertical = 1;
3015    vertarg[0] = *bottomhgt;
3016    vertarg[1] = *hgtinc;
3017
3018 #ifdef UNDERSCORE
3019    return v5dcreate_
3020 #else
3021 #  ifdef _CRAY
3022    return V5DCREATE
3023 #  else
3024
3025    return v5dcreate
3026 #  endif
3027 #endif
3028                    ( name, numtimes, numvars, nr, nc, varnl,
3029                      varname, timestamp, datestamp, &compressmode,
3030                      &projection, projarg, &vertical, vertarg );
3031 }
3032
3033
3034
3035 /*
3036  * Set lowest levels for each variable (other than default of 0).
3037  * Input: lowlev - array [NumVars] of ints
3038  * Return:  1 = ok, 0 = error
3039  */
3040 #ifdef UNDERSCORE
3041    int v5dsetlowlev_
3042 #else
3043 #  ifdef _CRAY
3044      int V5DSETLOWLEV
3045 #  else
3046      int v5dsetlowlev
3047 #  endif
3048 #endif
3049           ( int *lowlev )
3050 {
3051    return v5dSetLowLev(lowlev);
3052 }
3053
3054
3055
3056 /*
3057  * Set the units for a variable.
3058  * Input: var - variable number in [1,NumVars]
3059  *        units - a character string
3060  * Return:  1 = ok, 0 = error
3061  */
3062 #ifdef UNDERSCORE
3063    int v5dsetunits_
3064 #else
3065 #  ifdef _CRAY
3066      int V5DSETUNITS
3067 #  else
3068      int v5dsetunits
3069 #  endif
3070 #endif
3071           ( int *var, char *name )
3072 {
3073    return v5dSetUnits( *var, name );
3074 }
3075
3076
3077
3078 /*
3079  * Write a grid of data to the file.
3080  * Input:  time - timestep in [1,NumTimes]
3081  *         var - timestep in [1,NumVars]
3082  *         data - array [nr*nc*nl] of floats
3083  * Return:  1 = ok, 0 = error
3084  */
3085 #ifdef UNDERSCORE
3086    int v5dwrite_
3087 #else
3088 #  ifdef _CRAY
3089      int V5DWRITE
3090 #  else
3091      int v5dwrite
3092 #  endif
3093 #endif
3094           ( const int *time, const int *var, const FLOAT *data )
3095 {
3096    return v5dWrite( *time, *var, data );
3097 }
3098
3099
3100
3101 /*
3102  * Specify the McIDAS GR3D file number and grid number which correspond
3103  * to the grid specified by time and var.
3104  * Input:  time, var - timestep and variable of grid (starting at 1)
3105  *         mcfile, mcgrid - McIDAS grid file number and grid number
3106  * Return:  1 = ok, 0 = errror (bad time or var)
3107  */
3108 #ifdef UNDERSCORE
3109    int v5dmcfile_
3110 #else
3111 #  ifdef _CRAY
3112      int V5DMCFILE
3113 #  else
3114      int v5dmcfile
3115 #  endif
3116 #endif
3117          ( const int *time, const int *var,
3118            const int *mcfile, const int *mcgrid )
3119 {
3120    if (*time<1 || *time>Simple->NumTimes) {
3121       printf("Bad time argument to v5dSetMcIDASgrid: %d\n", *time );
3122       return 0;
3123    }
3124    if (*var<1 || *var>Simple->NumVars) {
3125       printf("Bad var argument to v5dSetMcIDASgrid: %d\n", *var );
3126       return 0;
3127    }
3128
3129    Simple->McFile[*time-1][*var-1] = (short) *mcfile;
3130    Simple->McGrid[*time-1][*var-1] = (short) *mcgrid;
3131    return 1;
3132 }
3133
3134
3135
3136 /*
3137  * Close a simple v5d file.
3138  */
3139 #ifdef UNDERSCORE
3140    int v5dclose_( void )
3141 #else
3142 #  ifdef _CRAY
3143      int V5DCLOSE( void )
3144 #  else
3145      int v5dclose( void )
3146 #  endif
3147 #endif
3148 {
3149    return v5dClose();
3150 }