Commits (34)
2019-10-29 Uwe Schulzweida
2021-01-25 Uwe Schulzweida
* using CGRIBEX library version 1.9.5
* Version 1.9.10 released
2021-01-07 Uwe Schulzweida
* gridDefParamRLL: set x/y units to degrees
2020-12-03 Uwe Schulzweida
* gridDefParamLCC: write latitude_of_projection_origin to lat_0 (bug fix)
2020-11-27 Uwe Schulzweida
* struct record_t: changed data type of levelID from short to int
2020-11-26 Uwe Schulzweida
* cdi.h: removed function gaussianLatitudes() from cdi interface
2020-10-09 Uwe Schulzweida
* using CGRIBEX library version 1.9.5
* Version 1.9.9 released
2020-09-25 Uwe Schulzweida
* cdi_generate_vars: set level to level1 if zaxis bounds are available (removed level = (level1 + level2) / 2)
2020-09-21 Uwe Schulzweida
* Ignore netCDF attribute _NCProperties
......
......@@ -4,10 +4,10 @@
# autoconf 2.69
# libtool 2.4.6
AC_PREREQ([2.69])
AC_PREREQ([2.70])
LT_PREREQ([2.4.6])
AC_INIT([cdi], [1.9.9rc5], [https://mpimet.mpg.de/cdi])
AC_INIT([cdi],[1.9.10],[https://mpimet.mpg.de/cdi])
AC_DEFINE_UNQUOTED(CDI, ["$PACKAGE_VERSION"], [CDI version])
......@@ -27,7 +27,7 @@ AC_CONFIG_HEADERS([src/config.h])
AM_MAINTAINER_MODE([disable])
# Check building environment
AC_PROG_CC_C99
AC_PROG_CC
dnl The check above forgets to delete the following directory
dnl generated by the MacOS linker, which results in multiple
dnl annoying messages from the rm utility:
......@@ -439,10 +439,11 @@ AC_CONFIG_FILES([tests/test_cksum_grib \
tables/gen_tableheaderfile \
util/serialrun],[chmod a+x "$ac_file"])
AC_OUTPUT([Makefile src/Makefile interfaces/Makefile app/Makefile \
AC_CONFIG_FILES([Makefile src/Makefile interfaces/Makefile app/Makefile \
tests/Makefile examples/Makefile cdi.settings \
examples/pio/Makefile src/pkgconfig/cdi.pc src/pkgconfig/cdipio.pc \
src/pkgconfig/cdi_f2003.pc])
AC_OUTPUT
# ----------------------------------------------------------------------
# Show configuration
......
......@@ -64,7 +64,7 @@ libcdi_la_SOURCES = \
extralib.c \
file.c \
file.h \
gaussgrid.c \
gaussian_latitudes.c \
gribapi.c \
gribapi.h \
gribapi_utilities.h \
......
......@@ -7,6 +7,7 @@
UINT32 get_UINT32(unsigned char *x)
{
// clang-format off
switch (HOST_ENDIANNESS)
{
case CDI_BIGENDIAN:
......@@ -17,11 +18,13 @@ UINT32 get_UINT32(unsigned char *x)
Error("Unhandled endianness %d", HOST_ENDIANNESS);
return UINT32_C(0xFFFFFFFF);
}
// clang-format on
}
UINT32 get_SUINT32(unsigned char *x)
{
// clang-format off
switch (HOST_ENDIANNESS)
{
case CDI_BIGENDIAN:
......@@ -32,11 +35,13 @@ UINT32 get_SUINT32(unsigned char *x)
Error("Unhandled endianness %d", HOST_ENDIANNESS);
return UINT32_C(0xFFFFFFFF);
}
// clang-format on
}
UINT64 get_UINT64(unsigned char *x)
{
// clang-format off
switch (HOST_ENDIANNESS)
{
case CDI_BIGENDIAN:
......@@ -49,11 +54,13 @@ UINT64 get_UINT64(unsigned char *x)
Error("Unhandled endianness %d", HOST_ENDIANNESS);
return UINT64_C(0xFFFFFFFFFFFFFFFF);
}
// clang-format on
}
UINT64 get_SUINT64(unsigned char *x)
{
// clang-format off
switch (HOST_ENDIANNESS)
{
case CDI_BIGENDIAN:
......@@ -66,6 +73,7 @@ UINT64 get_SUINT64(unsigned char *x)
Error("Unhandled endianness %d", HOST_ENDIANNESS);
return UINT64_C(0xFFFFFFFFFFFFFFFF);
}
// clang-format on
}
......@@ -85,8 +93,8 @@ size_t binReadF77Block(int fileID, int byteswap)
void binWriteF77Block(int fileID, int byteswap, size_t blocksize)
{
static unsigned int s[4] = {0, 8, 16, 24};
unsigned long ublocksize = (unsigned long) blocksize;
static const unsigned int s[4] = {0, 8, 16, 24};
const unsigned long ublocksize = (unsigned long) blocksize;
unsigned char f77block[4];
switch (HOST_ENDIANNESS)
......@@ -142,28 +150,6 @@ int binReadInt64(int fileID, int byteswap, size_t size, INT64 *ptr)
}
int binReadFlt32(int fileID, int byteswap, size_t size, FLT32 *ptr)
{
if ( sizeof(FLT32) != 4 ) Error("Not implemented for %d byte float!", sizeof(FLT32));
fileRead(fileID, (void *) ptr, 4*size);
if ( byteswap ) swap4byte(ptr, size);
return 0;
}
int binReadFlt64(int fileID, int byteswap, size_t size, FLT64 *ptr)
{
if ( sizeof(FLT64) != 8 ) Error("Not implemented for %d byte float!", sizeof(FLT64));
fileRead(fileID, (void *) ptr, 8*size);
if ( byteswap ) swap8byte(ptr, size);
return 0;
}
int binWriteInt32(int fileID, int byteswap, size_t size, INT32 *ptr)
{
if ( sizeof(INT32) != 4 ) Error("Not implemented for %d byte integer!", sizeof(INT32));
......
......@@ -2,7 +2,7 @@
#define BINARY_H
#ifdef HAVE_CONFIG_H
# include "config.h"
#include "config.h"
#endif
#include <inttypes.h>
......@@ -34,9 +34,6 @@ int binReadInt64(int fileID, int byteswap, size_t size, INT64 *ptr);
int binWriteInt32(int fileID, int byteswap, size_t size, INT32 *ptr);
int binWriteInt64(int fileID, int byteswap, size_t size, INT64 *ptr);
int binReadFlt32(int fileID, int byteswap, size_t size, FLT32 *ptr);
int binReadFlt64(int fileID, int byteswap, size_t size, FLT64 *ptr);
int binWriteFlt32(int fileID, int byteswap, size_t size, FLT32 *ptr);
int binWriteFlt64(int fileID, int byteswap, size_t size, FLT64 *ptr);
......
......@@ -27,25 +27,12 @@ int days_per_month(int calendar, int year, int month)
{
int daysperyear = calendar_dpy(calendar);
const int *dpm;
if ( daysperyear == 360 ) dpm = month_360;
else if ( daysperyear == 365 ) dpm = month_365;
else dpm = month_366;
int dayspermonth = 0;
if ( month >= 1 && month <= 12 )
dayspermonth = dpm[month-1];
/*
else
fprintf(stderr, "days_per_month: month %d out of range\n", month);
*/
if ( daysperyear == 0 && month == 2 )
{
if ( (year%4 == 0 && year%100 != 0) || year%400 == 0 )
dayspermonth = 29;
else
dayspermonth = 28;
}
const int *dpm = (daysperyear == 360) ? month_360 : ((daysperyear == 365) ? month_365 : month_366);
int dayspermonth = ( month >= 1 && month <= 12 ) ? dpm[month-1] : 0;
if (daysperyear == 0 && month == 2)
dayspermonth = ((year%4 == 0 && year%100 != 0) || year%400 == 0) ? 29 : 28;
return dayspermonth;
}
......@@ -55,11 +42,11 @@ int days_per_year(int calendar, int year)
{
int daysperyear = calendar_dpy(calendar);
if ( daysperyear == 0 )
if (daysperyear == 0)
{
if ( year == 1582 && (calendar == CALENDAR_STANDARD || calendar == CALENDAR_GREGORIAN) )
if (year == 1582 && (calendar == CALENDAR_STANDARD || calendar == CALENDAR_GREGORIAN))
daysperyear = 355;
else if ( (year%4 == 0 && year%100 != 0) || year%400 == 0 )
else if ((year%4 == 0 && year%100 != 0) || year%400 == 0)
daysperyear = 366;
else
daysperyear = 365;
......@@ -71,8 +58,6 @@ int days_per_year(int calendar, int year)
static void decode_day(int dpy, int days, int *year, int *month, int *day)
{
int i = 0;
*year = (days-1) / dpy;
days -= (*year*dpy);
......@@ -81,6 +66,7 @@ static void decode_day(int dpy, int days, int *year, int *month, int *day)
else if ( dpy == 365 ) dpm = month_365;
else if ( dpy == 366 ) dpm = month_366;
int i = 0;
if ( dpm )
for ( i = 0; i < 12; i++ )
{
......@@ -112,7 +98,7 @@ static int64_t encode_day(int dpy, int year, int month, int day)
void encode_caldaysec(int calendar, int year, int month, int day, int hour, int minute, int second,
int64_t *julday, int *secofday)
{
int dpy = calendar_dpy(calendar);
const int dpy = calendar_dpy(calendar);
if ( dpy == 360 || dpy == 365 || dpy == 366 )
*julday = encode_day(dpy, year, month, day);
......@@ -126,7 +112,7 @@ void encode_caldaysec(int calendar, int year, int month, int day, int hour, int
void decode_caldaysec(int calendar, int64_t julday, int secofday,
int *year, int *month, int *day, int *hour, int *minute, int *second)
{
int dpy = calendar_dpy(calendar);
const int dpy = calendar_dpy(calendar);
if ( dpy == 360 || dpy == 365 || dpy == 366 )
decode_day(dpy, julday, year, month, day);
......@@ -142,7 +128,7 @@ void decode_caldaysec(int calendar, int64_t julday, int secofday,
#ifdef TEST
static int date_to_calday(int calendar, int date)
{
int dpy = calendar_dpy(calendar);
const int dpy = calendar_dpy(calendar);
int year, month, day;
cdiDecodeDate(date, &year, &month, &day);
......@@ -160,14 +146,14 @@ static int date_to_calday(int calendar, int date)
static int calday_to_date(int calendar, int calday)
{
int year, month, day;
int dpy = calendar_dpy(calendar);
const int dpy = calendar_dpy(calendar);
if ( dpy == 360 || dpy == 365 || dpy == 366 )
decode_day(dpy, calday, &year, &month, &day);
else
decode_julday(calendar, calday, &year, &month, &day);
int date = cdiEncodeDate(year, month, day);
const int date = cdiEncodeDate(year, month, day);
return date;
}
......
......@@ -103,7 +103,7 @@ static int cdfOpenFile(const char *filename, const char *mode, int *filetype)
{
case 'r':
{
int status = cdf_open(filename, readmode, &ncid);
const int status = cdf_open(filename, readmode, &ncid);
if ( status > 0 && ncid < 0 ) ncid = CDI_ESYSTEM;
#ifdef HAVE_NETCDF4
else
......
This diff is collapsed.
......@@ -10,6 +10,7 @@
#include "dmemory.h"
#include "cdi.h"
#include "cdi_int.h"
#include "stream_grb.h"
#include "stream_cdf.h"
#include "cdf_int.h"
#include "vlist.h"
......@@ -83,13 +84,13 @@ void cdfGetSlapDescription(stream_t *streamptr, int varID, size_t (*start)[4], s
Message("dim = %d start = %d count = %d", idim, start[idim], count[idim]);
}
//Scans the data array for missVals, optionally applying first a scale factor and then an offset.
//Returns the number of missing + out-of-range values encountered.
// Scans the data array for missVals, optionally applying first a scale factor and then an offset.
// Returns the number of missing + out-of-range values encountered.
static
size_t cdfDoInputDataTransformationDP(int vlistID, int varID, size_t valueCount, double *data)
{
double missVal = vlistInqVarMissval(vlistID, varID);
const bool haveMissVal = vlistInqVarMissvalUsed(vlistID, varID);
const int haveMissVal = vlistInqVarMissvalUsed(vlistID, varID);
double validRange[2];
if (!(haveMissVal && vlistInqVarValidrange(vlistID, varID, validRange)))
validRange[0] = DBL_MIN, validRange[1] = DBL_MAX;
......@@ -98,20 +99,20 @@ size_t cdfDoInputDataTransformationDP(int vlistID, int varID, size_t valueCount,
double offset = vlistInqVarAddoffset(vlistID, varID);
double scaleFactor = vlistInqVarScalefactor(vlistID, varID);
const bool haveOffset = IS_NOT_EQUAL(offset, 0);
const bool haveScaleFactor = IS_NOT_EQUAL(scaleFactor, 1);
const bool missValIsNaN = DBL_IS_NAN(missVal);
const int haveOffset = IS_NOT_EQUAL(offset, 0.0);
const int haveScaleFactor = IS_NOT_EQUAL(scaleFactor, 1.0);
size_t missValCount = 0;
if ( IS_EQUAL(validMin, VALIDMISS) ) validMin = DBL_MIN;
if ( IS_EQUAL(validMax, VALIDMISS) ) validMax = DBL_MAX;
bool haveRangeCheck = (IS_NOT_EQUAL(validMax, DBL_MAX)) | (IS_NOT_EQUAL(validMin,DBL_MIN));
const int haveRangeCheck = (IS_NOT_EQUAL(validMax, DBL_MAX)) | (IS_NOT_EQUAL(validMin, DBL_MIN));
assert(!haveRangeCheck || haveMissVal);
switch ((int)haveMissVal | ((int)haveScaleFactor << 1)
| ((int)haveOffset << 2) | ((int)haveRangeCheck << 3))
switch (haveMissVal | (haveScaleFactor << 1) | (haveOffset << 2) | (haveRangeCheck << 3))
{
case 15: /* haveRangeCheck & haveMissVal & haveScaleFactor & haveOffset */
case 15: // haveRangeCheck & haveMissVal & haveScaleFactor & haveOffset
for ( size_t i = 0; i < valueCount; i++ )
{
int outOfRange = data[i] < validMin || data[i] > validMax;
......@@ -121,7 +122,7 @@ size_t cdfDoInputDataTransformationDP(int vlistID, int varID, size_t valueCount,
: isMissVal ? data[i] : data[i] * scaleFactor + offset;
}
break;
case 13: /* haveRangeCheck & haveMissVal & haveOffset */
case 13: // haveRangeCheck & haveMissVal & haveOffset
for ( size_t i = 0; i < valueCount; i++ )
{
int outOfRange = data[i] < validMin || data[i] > validMax;
......@@ -131,7 +132,7 @@ size_t cdfDoInputDataTransformationDP(int vlistID, int varID, size_t valueCount,
: isMissVal ? data[i] : data[i] + offset;
}
break;
case 11: /* haveRangeCheck & haveMissVal & haveScaleFactor */
case 11: // haveRangeCheck & haveMissVal & haveScaleFactor
for ( size_t i = 0; i < valueCount; i++ )
{
int outOfRange = data[i] < validMin || data[i] > validMax;
......@@ -141,7 +142,7 @@ size_t cdfDoInputDataTransformationDP(int vlistID, int varID, size_t valueCount,
: isMissVal ? data[i] : data[i] * scaleFactor;
}
break;
case 9: /* haveRangeCheck & haveMissVal */
case 9: // haveRangeCheck & haveMissVal
for ( size_t i = 0; i < valueCount; i++ )
{
int outOfRange = data[i] < validMin || data[i] > validMax;
......@@ -150,42 +151,50 @@ size_t cdfDoInputDataTransformationDP(int vlistID, int varID, size_t valueCount,
data[i] = outOfRange ? missVal : data[i];
}
break;
case 7: /* haveMissVal & haveScaleFactor & haveOffset */
case 7: // haveMissVal & haveScaleFactor & haveOffset
for ( size_t i = 0; i < valueCount; i++ )
if ( DBL_IS_EQUAL(data[i], missVal) )
missValCount++;
else
data[i] = data[i] * scaleFactor + offset;
break;
case 6: /* haveOffset & haveScaleFactor */
case 6: // haveOffset & haveScaleFactor
for ( size_t i = 0; i < valueCount; i++ )
data[i] = data[i] * scaleFactor + offset;
break;
case 5: /* haveMissVal & haveOffset */
case 5: // haveMissVal & haveOffset
for ( size_t i = 0; i < valueCount; i++ )
if ( DBL_IS_EQUAL(data[i], missVal) )
missValCount++;
else
data[i] += offset;
break;
case 4: /* haveOffset */
case 4: // haveOffset
for ( size_t i = 0; i < valueCount; i++ )
data[i] += offset;
break;
case 3: /* haveMissVal & haveScaleFactor */
case 3: // haveMissVal & haveScaleFactor
for ( size_t i = 0; i < valueCount; i++ )
if ( DBL_IS_EQUAL(data[i], missVal) )
missValCount++;
else
data[i] *= scaleFactor;
break;
case 2: /* haveScaleFactor */
case 2: // haveScaleFactor
for ( size_t i = 0; i < valueCount; i++ )
data[i] *= scaleFactor;
break;
case 1: /* haveMissVal */
for ( size_t i = 0; i < valueCount; i++ )
missValCount += (unsigned)DBL_IS_EQUAL(data[i], missVal);
case 1: // haveMissVal
if (missValIsNaN)
{
for ( size_t i = 0; i < valueCount; i++ )
missValCount += (size_t)DBL_IS_NAN(data[i]);
}
else
{
for ( size_t i = 0; i < valueCount; i++ )
missValCount += (size_t)IS_EQUAL(data[i], missVal);
}
break;
}
......@@ -196,7 +205,7 @@ static
size_t cdfDoInputDataTransformationSP(int vlistID, int varID, size_t valueCount, float *data)
{
double missVal = vlistInqVarMissval(vlistID, varID);
const bool haveMissVal = vlistInqVarMissvalUsed(vlistID, varID);
const int haveMissVal = vlistInqVarMissvalUsed(vlistID, varID);
double validRange[2];
if (!(haveMissVal && vlistInqVarValidrange(vlistID, varID, validRange)))
validRange[0] = DBL_MIN, validRange[1] = DBL_MAX;
......@@ -205,20 +214,20 @@ size_t cdfDoInputDataTransformationSP(int vlistID, int varID, size_t valueCount,
double offset = vlistInqVarAddoffset(vlistID, varID);
double scaleFactor = vlistInqVarScalefactor(vlistID, varID);
const bool haveOffset = IS_NOT_EQUAL(offset, 0);
const bool haveScaleFactor = IS_NOT_EQUAL(scaleFactor, 1);
const bool missValIsNaN = DBL_IS_NAN(missVal);
const int haveOffset = IS_NOT_EQUAL(offset, 0.0);
const int haveScaleFactor = IS_NOT_EQUAL(scaleFactor, 1.0);
size_t missValCount = 0;
if ( IS_EQUAL(validMin, VALIDMISS) ) validMin = DBL_MIN;
if ( IS_EQUAL(validMax, VALIDMISS) ) validMax = DBL_MAX;
bool haveRangeCheck = (IS_NOT_EQUAL(validMax, DBL_MAX)) | (IS_NOT_EQUAL(validMin,DBL_MIN));
const int haveRangeCheck = (IS_NOT_EQUAL(validMax, DBL_MAX)) | (IS_NOT_EQUAL(validMin, DBL_MIN));
assert(!haveRangeCheck || haveMissVal);
switch ((int)haveMissVal | ((int)haveScaleFactor << 1)
| ((int)haveOffset << 2) | ((int)haveRangeCheck << 3))
switch (haveMissVal | (haveScaleFactor << 1) | (haveOffset << 2) | (haveRangeCheck << 3))
{
case 15: /* haveRangeCheck & haveMissVal & haveScaleFactor & haveOffset */
case 15: // haveRangeCheck & haveMissVal & haveScaleFactor & haveOffset
for ( size_t i = 0; i < valueCount; i++ )
{
int outOfRange = data[i] < validMin || data[i] > validMax;
......@@ -228,7 +237,7 @@ size_t cdfDoInputDataTransformationSP(int vlistID, int varID, size_t valueCount,
: isMissVal ? data[i] : (float)(data[i] * scaleFactor + offset);
}
break;
case 13: /* haveRangeCheck & haveMissVal & haveOffset */
case 13: // haveRangeCheck & haveMissVal & haveOffset
for ( size_t i = 0; i < valueCount; i++ )
{
int outOfRange = data[i] < validMin || data[i] > validMax;
......@@ -238,7 +247,7 @@ size_t cdfDoInputDataTransformationSP(int vlistID, int varID, size_t valueCount,
: isMissVal ? data[i] : (float)(data[i] + offset);
}
break;
case 11: /* haveRangeCheck & haveMissVal & haveScaleFactor */
case 11: // haveRangeCheck & haveMissVal & haveScaleFactor
for ( size_t i = 0; i < valueCount; i++ )
{
int outOfRange = data[i] < validMin || data[i] > validMax;
......@@ -248,7 +257,7 @@ size_t cdfDoInputDataTransformationSP(int vlistID, int varID, size_t valueCount,
: isMissVal ? data[i] : (float)(data[i] * scaleFactor);
}
break;
case 9: /* haveRangeCheck & haveMissVal */
case 9: // haveRangeCheck & haveMissVal
for ( size_t i = 0; i < valueCount; i++ )
{
int outOfRange = data[i] < validMin || data[i] > validMax;
......@@ -257,42 +266,50 @@ size_t cdfDoInputDataTransformationSP(int vlistID, int varID, size_t valueCount,
data[i] = outOfRange ? (float)missVal : data[i];
}
break;
case 7: /* haveMissVal & haveScaleFactor & haveOffset */
case 7: // haveMissVal & haveScaleFactor & haveOffset
for ( size_t i = 0; i < valueCount; i++ )
if ( DBL_IS_EQUAL(data[i], missVal) )
missValCount++;
else
data[i] = (float)(data[i] * scaleFactor + offset);
break;
case 6: /* haveOffset & haveScaleFactor */
case 6: // haveOffset & haveScaleFactor
for ( size_t i = 0; i < valueCount; i++ )
data[i] = (float)(data[i] * scaleFactor + offset);
break;
case 5: /* haveMissVal & haveOffset */
case 5: // haveMissVal & haveOffset
for ( size_t i = 0; i < valueCount; i++ )
if ( DBL_IS_EQUAL(data[i], missVal) )
missValCount++;
else
data[i] = (float)(data[i] + offset);
break;
case 4: /* haveOffset */
case 4: // haveOffset
for ( size_t i = 0; i < valueCount; i++ )
data[i] = (float)(data[i] + offset);
break;
case 3: /* haveMissVal & haveScaleFactor */
case 3: // haveMissVal & haveScaleFactor
for ( size_t i = 0; i < valueCount; i++ )
if ( DBL_IS_EQUAL(data[i], missVal) )
missValCount++;
else
data[i] = (float)(data[i] * scaleFactor);
break;
case 2: /* haveScaleFactor */
case 2: // haveScaleFactor
for ( size_t i = 0; i < valueCount; i++ )
data[i] = (float)(data[i] * scaleFactor);
break;
case 1: /* haveMissVal */
for ( size_t i = 0; i < valueCount; i++ )
missValCount += (unsigned)DBL_IS_EQUAL(data[i], missVal);
case 1: // haveMissVal
if (missValIsNaN)
{
for ( size_t i = 0; i < valueCount; i++ )
missValCount += (size_t)DBL_IS_NAN(data[i]);
}
else
{
for ( size_t i = 0; i < valueCount; i++ )
missValCount += (size_t)IS_EQUAL(data[i], missVal);
}
break;
}
......@@ -773,9 +790,8 @@ void cdfReadVarSliceDPPart(stream_t *streamptr, int varID, int levelID, int varT
const size_t gridsize = gridInqSize(gridId);
unsigned int position = 0;
for (int i=0 ; i<4 ; i++)
if (count[i] == gridsize)
position = i;
for (int i = 0; i < 4; i++)
if (count[i] == gridsize) position = i;
start[position] = start[position]+startpoint;
count[position] = length;
......@@ -804,9 +820,8 @@ void cdfReadVarSliceSPPart(stream_t *streamptr, int varID, int levelID, int varT
size_t gridsize = gridInqSize(gridId);
unsigned int position = 0;
for (int i=0 ; i<4 ; i++)
if (count[i] == gridsize)
position = i;
for (int i = 0; i < 4; i++)
if (count[i] == gridsize) position = i;
start[position] = start[position]+startpoint;
count[position] = length;
......@@ -836,7 +851,15 @@ int cdiStreamReadVarSlicePart(int streamID, int varID, int levelID, int varType,
// currently we only care for netcdf data
switch (filetype)
{
#if defined (HAVE_LIBNETCDF)
#ifdef HAVE_LIBGRIB
case CDI_FILETYPE_GRB:
case CDI_FILETYPE_GRB2:
{
grb_read_var_slice(streamptr, varID, levelID, memtype, data, nmiss);
break;
}
#endif
#ifdef HAVE_LIBNETCDF
case CDI_FILETYPE_NC:
case CDI_FILETYPE_NC2:
case CDI_FILETYPE_NC4:
......@@ -919,7 +942,15 @@ void cdiStreamReadVarPart(int streamID, int varID, int varType, int start, size_
// currently we only care for netcdf data
switch (filetype)
{
#if defined (HAVE_LIBNETCDF)
#ifdef HAVE_LIBGRIB
case CDI_FILETYPE_GRB:
case CDI_FILETYPE_GRB2:
{
grb_read_var(streamptr, varID, memtype, data, nmiss);
break;
}
#endif
#ifdef HAVE_LIBNETCDF
case CDI_FILETYPE_NC:
case CDI_FILETYPE_NC2:
case CDI_FILETYPE_NC4:
......@@ -955,5 +986,4 @@ void streamReadVarPart(int streamID, int varID, int varType, int start, size_t s
cdiStreamReadVarPart(streamID, varID, varType, start, size, memtype, data, nmiss);
}
#endif
#endif /* HAVE_LIBNETCDF */
......@@ -23,8 +23,8 @@ bool strStartsWith(const char *vstr, const char *cstr)
bool is_equal = false;
if (vstr && cstr)
{
size_t clen = strlen(cstr);
size_t vlen = strlen(vstr);
const size_t clen = strlen(cstr);
const size_t vlen = strlen(vstr);
if (clen <= vlen) is_equal = (memcmp(vstr, cstr, clen) == 0);
}
return is_equal;
......@@ -74,8 +74,8 @@ bool is_timeaxis_units(const char *timeunits)
{
bool status = false;
size_t len = strlen(timeunits);
char *tu = (char *) Malloc((len+1)*sizeof(char));
const size_t len = strlen(timeunits);
char *tu = (char *) malloc((len+1)*sizeof(char));
memcpy(tu, timeunits, (len+1) * sizeof(char));
char *ptu = tu;
......@@ -96,7 +96,7 @@ bool is_timeaxis_units(const char *timeunits)
}
}
Free(tu);
free(tu);
return status;
}
......
......@@ -241,11 +241,11 @@ void cdfDefineCellMethods(stream_t *streamptr, int cdiID, int varID, int fileID,
taxis_t *taxis = &streamptr->tsteps[0].taxis;
if (!taxis->has_bounds) return;
int time_varid = streamptr->basetime.ncvarid;
const int time_varid = streamptr->basetime.ncvarid;
char timeVarName[CDI_MAX_NAME];
cdf_inq_varname(fileID, time_varid, timeVarName);
int stepType = vlistInqVarTsteptype(cdiID, varID);
const int stepType = vlistInqVarTsteptype(cdiID, varID);
const char *cellMethod = NULL;
if (stepType == TSTEP_AVG) cellMethod = "mean";
......@@ -282,7 +282,7 @@ void cdfDefineAttributes(int cdiID, int varID, int fileID, int ncvarID)
if ( atttype == CDI_DATATYPE_TXT )
{
size_t attSize = (size_t)attlen*sizeof(char);
const size_t attSize = (size_t)attlen*sizeof(char);
char *atttxt = (char *)resizeBuf(&attBuf, &attBufSize, attSize);
cdiInqAttTxt(cdiID, varID, attname, attlen, atttxt);
len = (size_t)attlen;
......@@ -292,7 +292,7 @@ void cdfDefineAttributes(int cdiID, int varID, int fileID, int ncvarID)
atttype == CDI_DATATYPE_INT16 || atttype == CDI_DATATYPE_UINT16 ||
atttype == CDI_DATATYPE_INT32 || atttype == CDI_DATATYPE_UINT32 )
{
size_t attSize = (size_t)attlen*sizeof(int);
const size_t attSize = (size_t)attlen*sizeof(int);
int *attint = (int *)resizeBuf(&attBuf, &attBufSize, attSize);
cdiInqAttInt(cdiID, varID, attname, attlen, &attint[0]);
len = (size_t)attlen;
......@@ -308,7 +308,7 @@ void cdfDefineAttributes(int cdiID, int varID, int fileID, int ncvarID)
}
else if ( atttype == CDI_DATATYPE_FLT32 || atttype == CDI_DATATYPE_FLT64 )
{
size_t attSize = (size_t)attlen * sizeof(double);
const size_t attSize = (size_t)attlen * sizeof(double);
double *attflt = (double *)resizeBuf(&attBuf, &attBufSize, attSize);
cdiInqAttFlt(cdiID, varID, attname, attlen, attflt);
len = (size_t)attlen;
......@@ -372,7 +372,7 @@ void cdf_get_gmapvarname(int gridID, char *gmapvarname)
if ( !gmapname[0] )
{
int projID = gridInqProj(gridID);
const int projID = gridInqProj(gridID);
if ( projID != CDI_UNDEFID )
{
pgridID = projID;
......@@ -528,7 +528,7 @@ void cdfDefineCoordinates(const stream_t *streamptr, int ncvarID, int nczvarID,
{
char name[CDI_MAX_NAME]; name[0] = 0;
cdf_inq_varname(fileID, ncyvarID, name);
size_t len = strlen(name);
const size_t len = strlen(name);
cdf_put_att_text(fileID, ncvarID, "CDI_grid_latitudes", len, name);
}
......@@ -537,7 +537,7 @@ void cdfDefineCoordinates(const stream_t *streamptr, int ncvarID, int nczvarID,
{
char name[CDI_MAX_NAME]; name[0] = 0;
cdf_inq_varname(fileID, ncrpvarID, name);
size_t len = strlen(name);
const size_t len = strlen(name);
cdf_put_att_text(fileID, ncvarID, "CDI_grid_reduced_points", len, name);
}
}
......
......@@ -206,7 +206,7 @@ typedef struct {
#define TIME_CONSTANT 0 // Time constant
#define TIME_VARYING 1 // Time varying
#define TIME_VARIABLE 1 // obsolate, use TIME_VARYING
#define TIME_VARIABLE 1 // obsolete, use TIME_VARYING
// TSTEP types
......@@ -658,21 +658,21 @@ int vlistFindLevel(int vlistID, int fvarID, int flevelID);
int vlistMergedVar(int vlistID, int varID);
int vlistMergedLevel(int vlistID, int varID, int levelID);
// cdiClearAdditionalKeys: Clear the list of additional GRIB keys
// cdiClearAdditionalKeys: Clear the list of additional GRIB keys
void cdiClearAdditionalKeys(void);
// cdiDefAdditionalKey: Register an additional GRIB key which is read when file is opened
// cdiDefAdditionalKey: Register an additional GRIB key which is read when file is opened
void cdiDefAdditionalKey(const char *string);
// vlistDefVarIntKey: Set an arbitrary keyword/integer value pair for GRIB API
// vlistDefVarIntKey: Set an arbitrary keyword/integer value pair for GRIB API
void vlistDefVarIntKey(int vlistID, int varID, const char *name, int value);
// vlistDefVarDblKey: Set an arbitrary keyword/double value pair for GRIB API
// vlistDefVarDblKey: Set an arbitrary keyword/double value pair for GRIB API
void vlistDefVarDblKey(int vlistID, int varID, const char *name, double value);
// vlistHasVarKey: returns 1 if meta-data key was read, 0 otherwise
// vlistHasVarKey: returns 1 if meta-data key was read, 0 otherwise
int vlistHasVarKey(int vlistID, int varID, const char *name);
// vlistInqVarDblKey: raw access to GRIB meta-data
// vlistInqVarDblKey: raw access to GRIB meta-data
double vlistInqVarDblKey(int vlistID, int varID, const char *name);
// vlistInqVarIntKey: raw access to GRIB meta-data
// vlistInqVarIntKey: raw access to GRIB meta-data
int vlistInqVarIntKey(int vlistID, int varID, const char *name);
// CDI attributes
......@@ -1273,7 +1273,21 @@ int vlistInqVarSubtype(int vlistID, int varID);
void gribapiLibraryVersion(int *major_version, int *minor_version, int *revision_version);
void gaussianLatitudes(double *latitudes, double *weights, size_t nlat);
#ifdef __cplusplus
}
#endif
//FINT_OFF <--- don't change or remove this line!!!
#ifdef __cplusplus
extern "C" {
#endif
#define HAVE_CDI_PROJ_FUNCS 1
extern int (*proj_lonlat_to_lcc_func)(double, double, double, double, double, double, double, size_t, double*, double*);
extern int (*proj_lcc_to_lonlat_func)(double, double, double, double, double, double, double, double, double, size_t, double*, double*);
extern int (*proj_lonlat_to_stere_func)(double, double, double, double, double, size_t, double*, double*);
extern int (*proj_stere_to_lonlat_func)(double, double, double, double, double, double, double, size_t, double*, double*);
#ifdef __cplusplus
}
......
......@@ -4,7 +4,7 @@
!
! Author:
! -------
! Uwe Schulzweida, MPI-MET, Hamburg, June 2020
! Uwe Schulzweida, MPI-MET, Hamburg, January 2021
!
INTEGER CDI_MAX_NAME
......@@ -306,6 +306,8 @@
PARAMETER (ZAXIS_REFERENCE = 25)
INTEGER ZAXIS_CHAR
PARAMETER (ZAXIS_CHAR = 26)
INTEGER ZAXIS_TROPOPAUSE
PARAMETER (ZAXIS_TROPOPAUSE = 27)
!
! SUBTYPE types
!
......@@ -2615,9 +2617,3 @@
! INTEGER revision_version)
EXTERNAL gribapiLibraryVersion
! gaussianLatitudes
! (DOUBLEPRECISION latitudes,
! DOUBLEPRECISION weights,
! INTEGER nlat)
EXTERNAL gaussianLatitudes
/* Automatically generated by make_fint.c, don't edit! */
#ifdef HAVE_CONFIG_H
# include "config.h"
#include "config.h"
#endif
#ifndef CDI_H_
......@@ -848,11 +848,6 @@ FCALLSCFUN3 (INT, subtypeInqTile, SUBTYPEINQTILE, subtypeinqtile, INT, INT, INT)
FCALLSCFUN4 (INT, subtypeInqAttribute, SUBTYPEINQATTRIBUTE, subtypeinqattribute, INT, INT, STRING, PINT)
FCALLSCFUN2 (INT, vlistInqVarSubtype, VLISTINQVARSUBTYPE, vlistinqvarsubtype, INT, INT)
FCALLSCSUB3 (gribapiLibraryVersion, GRIBAPILIBRARYVERSION, gribapilibraryversion, PINT, PINT, PINT)
static void gaussianLatitudes_fwrap(double *latitudes, double *weights, int nlat)
{
gaussianLatitudes(latitudes, weights, (size_t)nlat);
}
FCALLSCSUB3 (gaussianLatitudes_fwrap, GAUSSIANLATITUDES, gaussianlatitudes, PDOUBLE, PDOUBLE, INT)
#if defined __clang__
# pragma GCC diagnostic pop
......
......@@ -183,7 +183,7 @@ typedef struct
int ltype;
short tsteptype;
short varID;
short levelID;
int levelID;
short used;
char varname[32]; // needed for grib decoding with GRIB_API
VarScanKeys scanKeys;
......
......@@ -182,10 +182,10 @@ int cdiInqKeyLen(int cdiID, int varID, int key, int *length)
{
int status = -1;
cdi_keys_t *keysp = cdi_get_keysp(cdiID, varID);
const cdi_keys_t *keysp = cdi_get_keysp(cdiID, varID);
xassert(keysp != NULL);
const cdi_key_t *keyp = find_key(keysp, key);
const cdi_key_t *keyp = find_key_const(keysp, key);
if ( keyp != NULL )
{
*length = keyp->length;
......@@ -260,7 +260,7 @@ int cdiCopyVarKey(const cdi_keys_t *keysp1, int key, cdi_keys_t *keysp2)
{
int status = CDI_NOERR;
const cdi_key_t *keyp = find_key_const (keysp1, key);
const cdi_key_t *keyp = find_key_const(keysp1, key);
if (keyp == NULL) return -1;
cdi_define_key(keyp, keysp2);
......@@ -352,10 +352,10 @@ int cdiInqKeyInt(int cdiID, int varID, int key, int *value)
// if ( varID != CDI_GLOBAL ) status = cdiInqKeyInt(cdiID, CDI_GLOBAL, key, value);
cdi_keys_t *keysp = cdi_get_keysp(cdiID, varID);
const cdi_keys_t *keysp = cdi_get_keysp(cdiID, varID);
xassert(keysp != NULL);
cdi_key_t *keyp = find_key(keysp, key);
const cdi_key_t *keyp = find_key_const(keysp, key);
if ( keyp != NULL ) // key in use
{
if ( keyp->type == KEY_INT )
......@@ -451,10 +451,10 @@ int cdiInqKeyFloat(int cdiID, int varID, int key, double *value)
// if ( varID != CDI_GLOBAL ) status = cdiInqKeyFloat(cdiID, CDI_GLOBAL, key, value);
cdi_keys_t *keysp = cdi_get_keysp(cdiID, varID);
const cdi_keys_t *keysp = cdi_get_keysp(cdiID, varID);
xassert(keysp != NULL);
cdi_key_t *keyp = find_key(keysp, key);
const cdi_key_t *keyp = find_key_const(keysp, key);
if ( keyp != NULL ) // key in use
{
if ( keyp->type == KEY_FLOAT )
......@@ -672,9 +672,9 @@ int cdiInqKeyString(int cdiID, int varID, int key, char *string, int *length)
return status;
}
const char *cdiInqVarKeyStringPtr(cdi_keys_t *keysp, int key)
const char *cdiInqVarKeyStringPtr(const cdi_keys_t *keysp, int key)
{
const cdi_key_t *keyp = find_key(keysp, key);
const cdi_key_t *keyp = find_key_const(keysp, key);
if ( keyp != NULL ) // key in use
{
if ( keyp->type == KEY_BYTES ) return (const char *)keyp->v.s;
......
......@@ -31,10 +31,10 @@ int cdiInqVarKeyInt(const cdi_keys_t *keysp, int key);
int cdiInqVarKeyBytes(const cdi_keys_t *keysp, int key, unsigned char *bytes, int *length);
cdi_key_t *find_key(cdi_keys_t *keysp, int key);
const char *cdiInqVarKeyStringPtr(cdi_keys_t *keysp, int key);
const char *cdiInqVarKeyStringPtr(const cdi_keys_t *keysp, int key);
static inline
const char *cdiInqVarKeyString(cdi_keys_t *keysp, int key)
const char *cdiInqVarKeyString(const cdi_keys_t *keysp, int key)
{
const char *string = cdiInqVarKeyStringPtr(keysp, key);
if (string == NULL) string = "";
......
......@@ -16,17 +16,20 @@ enum {
EXT_HEADER_LEN = 4,
};
union EXT_HEADER
{
INT32 i32[EXT_HEADER_LEN];
INT64 i64[EXT_HEADER_LEN];
};
static int initExtLib = 0;
static int extDefaultPrec = 0;
static int extDefaultNumber = EXT_REAL;
/*
* A version string.
*/
// A version string.
#undef LIBVERSION
#define LIBVERSION 1.4.1
#define LIBVERSION 1.4.2
#define XSTRING(x) #x
#define STRING(x) XSTRING(x)
static const char ext_libvers[] = STRING(LIBVERSION);
......@@ -37,13 +40,13 @@ const char *extLibraryVersion(void)
}
static int EXT_Debug = 0; /* If set to 1, debugging */
static int EXT_Debug = 0; // If set to 1, debugging
void extDebug(int debug)
{
if (debug) Message("debug level %d", debug);
EXT_Debug = debug;
if (EXT_Debug) Message("debug level %d", debug);
}
......@@ -141,12 +144,12 @@ int extCheckFiletype(int fileID, int *swap)
if ( fileRead(fileID, buffer, 4) != 4 ) return 0;
size_t blocklen = (size_t) get_UINT32(buffer);
size_t sblocklen = (size_t) get_SUINT32(buffer);
const size_t blocklen = (size_t) get_UINT32(buffer);
const size_t sblocklen = (size_t) get_SUINT32(buffer);
if ( EXT_Debug )
Message("blocklen = %d sblocklen = %d", blocklen, sblocklen);
if ( EXT_Debug ) Message("blocklen = %d sblocklen = %d", blocklen, sblocklen);
// clang-format off
if ( blocklen == 16 )
{
*swap = 0;
......@@ -179,6 +182,7 @@ int extCheckFiletype(int fileID, int *swap)
pbuf = buffer+3*fact; dimxy = (size_t) get_SUINT64(pbuf);
pbuf = buffer+sblocklen+4; data = (size_t) get_SUINT32(pbuf);
}
// clang-format on
fileRewind(fileID);
......@@ -188,7 +192,7 @@ int extCheckFiletype(int fileID, int *swap)
Message("dimxy = %lu data = %lu", dimxy, data);
}
int found = data && (dimxy*fact == data || dimxy*fact*2 == data);
const int found = data && (dimxy*fact == data || dimxy*fact*2 == data);
return found;
}
......@@ -197,9 +201,9 @@ int extInqHeader(void *ext, int *header)
{
extrec_t *extp = (extrec_t *) ext;
for (size_t i = 0; i < EXT_HEADER_LEN; i++) header[i] = extp->header[i];
for (int i = 0; i < EXT_HEADER_LEN; i++) header[i] = extp->header[i];
if ( EXT_Debug ) Message("datasize = %lu", extp->datasize);
if ( EXT_Debug ) Message("datasize = %zu", extp->datasize);
return 0;
}
......@@ -209,12 +213,12 @@ int extDefHeader(void *ext, const int *header)
{
extrec_t *extp = (extrec_t *) ext;
for (size_t i = 0; i < EXT_HEADER_LEN; i++) extp->header[i] = header[i];
for (int i = 0; i < EXT_HEADER_LEN; i++) extp->header[i] = header[i];
extp->datasize = (size_t)header[3];
if ( extp->number == EXT_COMP ) extp->datasize *= 2;
if (EXT_Debug) Message("datasize = %lu", extp->datasize);
if (EXT_Debug) Message("datasize = %zu", extp->datasize);
return 0;
}
......@@ -223,10 +227,10 @@ static
int extInqData(extrec_t *extp, int prec, void *data)
{
int ierr = 0;
int byteswap = extp->byteswap;
size_t datasize = extp->datasize;
void *buffer = extp->buffer;
int rprec = extp->prec;
const int byteswap = extp->byteswap;
const size_t datasize = extp->datasize;
void *buffer = extp->buffer;
const int rprec = extp->prec;
switch ( rprec )
{
......@@ -291,26 +295,21 @@ static int extDefData(void *ext, int prec, const void *data)
{
extrec_t *extp = (extrec_t *) ext;
int rprec = extDefaultPrec ? extDefaultPrec : extp->prec;
const int rprec = extDefaultPrec ? extDefaultPrec : extp->prec;
extp->prec = rprec ? rprec : prec;
int *header = extp->header;
size_t datasize = (size_t)header[3];
if ( extp->number == EXT_COMP ) datasize *= 2;
size_t blocklen = datasize * (size_t)rprec;
const size_t blocklen = datasize * (size_t)rprec;
extp->datasize = datasize;
void *buffer = extp->buffer;
size_t buffersize = extp->buffersize;
if ( buffersize != blocklen )
if ( extp->buffersize != blocklen )
{
buffersize = blocklen;
buffer = Realloc(buffer, buffersize);
extp->buffer = buffer;
extp->buffersize = buffersize;
extp->buffersize = blocklen;
extp->buffer = Realloc(extp->buffer, extp->buffersize);
}
switch ( rprec )
......@@ -318,20 +317,20 @@ static int extDefData(void *ext, int prec, const void *data)
case EXSE_SINGLE_PRECISION:
{
if ( rprec == prec )
memcpy(buffer, data, datasize*sizeof(FLT32));
memcpy(extp->buffer, data, datasize*sizeof(FLT32));
else
for (size_t i = 0; i < datasize; i++)
((float *) buffer)[i] = (float) ((double *) data)[i];
((float *) extp->buffer)[i] = (float) ((double *) data)[i];
break;
}
case EXSE_DOUBLE_PRECISION:
{
if ( rprec == prec )
memcpy(buffer, data, datasize*sizeof(FLT64));
memcpy(extp->buffer, data, datasize*sizeof(FLT64));
else
for (size_t i = 0; i < datasize; i++)
((double *) buffer)[i] = (double) ((float *) data)[i];
((double *) extp->buffer)[i] = (double) ((float *) data)[i];
break;
}
......@@ -369,39 +368,32 @@ int extRead(int fileID, void *ext)
extp->checked = 1;
}
int byteswap = extp->byteswap;
const int byteswap = extp->byteswap;
/* read header record */
// read header record
size_t blocklen = binReadF77Block(fileID, byteswap);
if ( fileEOF(fileID) ) return -1;
if ( EXT_Debug ) Message("blocklen = %lu", blocklen);
size_t hprec = blocklen / EXT_HEADER_LEN;
const size_t hprec = blocklen / EXT_HEADER_LEN;
extp->prec = (int)hprec;
union EXT_HEADER tempheader;
switch ( hprec )
{
case EXSE_SINGLE_PRECISION:
{
INT32 tempheader[4];
binReadInt32(fileID, byteswap, EXT_HEADER_LEN, tempheader);
for (size_t i = 0; i < EXT_HEADER_LEN; i++)
extp->header[i] = (int)tempheader[i];
binReadInt32(fileID, byteswap, EXT_HEADER_LEN, tempheader.i32);
for (int i = 0; i < EXT_HEADER_LEN; i++) extp->header[i] = (int)tempheader.i32[i];
break;
}
case EXSE_DOUBLE_PRECISION:
{
INT64 tempheader[4];
binReadInt64(fileID, byteswap, EXT_HEADER_LEN, tempheader);
for (size_t i = 0; i < EXT_HEADER_LEN; i++)
extp->header[i] = (int)tempheader[i];
binReadInt64(fileID, byteswap, EXT_HEADER_LEN, tempheader.i64);
for (int i = 0; i < EXT_HEADER_LEN; i++) extp->header[i] = (int)tempheader.i64[i];
break;
}
default:
......@@ -421,19 +413,14 @@ int extRead(int fileID, void *ext)
extp->datasize = (size_t)extp->header[3];
if ( EXT_Debug ) Message("datasize = %lu", extp->datasize);
if ( EXT_Debug ) Message("datasize = %zu", extp->datasize);
blocklen = binReadF77Block(fileID, byteswap);
void *buffer = extp->buffer;
size_t buffersize = (size_t)extp->buffersize;
if ( buffersize < blocklen )
if ( extp->buffersize < blocklen )
{
buffersize = blocklen;
buffer = Realloc(buffer, buffersize);
extp->buffer = buffer;
extp->buffersize = buffersize;
extp->buffersize = blocklen;
extp->buffer = Realloc(extp->buffer, extp->buffersize);
}
size_t dprec = blocklen / extp->datasize;
......@@ -455,7 +442,7 @@ int extRead(int fileID, void *ext)
return -1;
}
fileRead(fileID, buffer, blocklen);
fileRead(fileID, extp->buffer, blocklen);
blocklen2 = binReadF77Block(fileID, byteswap);
......@@ -472,13 +459,13 @@ int extRead(int fileID, void *ext)
int extWrite(int fileID, void *ext)
{
extrec_t *extp = (extrec_t *) ext;
union { INT32 i32[EXT_HEADER_LEN]; INT64 i64[EXT_HEADER_LEN]; } tempheader;
int byteswap = extp->byteswap;
int rprec = extp->prec;
union EXT_HEADER tempheader;
const int byteswap = extp->byteswap;
const int rprec = extp->prec;
int number = extp->number;
int *header = extp->header;
/* write header record */
// write header record
size_t blocklen = EXT_HEADER_LEN * (size_t)rprec;
binWriteF77Block(fileID, byteswap, blocklen);
......@@ -487,20 +474,14 @@ int extWrite(int fileID, void *ext)
{
case EXSE_SINGLE_PRECISION:
{
for (size_t i = 0; i < EXT_HEADER_LEN; i++)
tempheader.i32[i] = (INT32) header[i];
for (int i = 0; i < EXT_HEADER_LEN; i++) tempheader.i32[i] = (INT32) header[i];
binWriteInt32(fileID, byteswap, EXT_HEADER_LEN, tempheader.i32);
break;
}
case EXSE_DOUBLE_PRECISION:
{
for (size_t i = 0; i < EXT_HEADER_LEN; i++)
tempheader.i64[i] = (INT64) header[i];
for (int i = 0; i < EXT_HEADER_LEN; i++) tempheader.i64[i] = (INT64) header[i];
binWriteInt64(fileID, byteswap, EXT_HEADER_LEN, tempheader.i64);
break;
}
default:
......@@ -512,26 +493,22 @@ int extWrite(int fileID, void *ext)
binWriteF77Block(fileID, byteswap, blocklen);
size_t datasize = (size_t)header[3];
if ( number == EXT_COMP ) datasize *= 2;
blocklen = datasize * (size_t)rprec;
extp->datasize = (size_t)header[3];
if ( number == EXT_COMP ) extp->datasize *= 2;
blocklen = extp->datasize * (size_t)rprec;
binWriteF77Block(fileID, byteswap, blocklen);
extp->datasize = datasize;
void *buffer = extp->buffer;
switch ( rprec )
{
case EXSE_SINGLE_PRECISION:
{
binWriteFlt32(fileID, byteswap, datasize, (FLT32 *) buffer);
binWriteFlt32(fileID, byteswap, extp->datasize, (FLT32 *) extp->buffer);
break;
}
case EXSE_DOUBLE_PRECISION:
{
binWriteFlt64(fileID, byteswap, datasize, (FLT64 *) buffer);
binWriteFlt64(fileID, byteswap, extp->datasize, (FLT64 *) extp->buffer);
break;
}
default:
......
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>
......@@ -10,48 +8,42 @@
#define M_SQRT2 1.41421356237309504880168872420969808
#endif
#include "dmemory.h"
#include "cdi_int.h"
static
void cpledn(size_t kn, size_t kodd, double *pfn, double pdx, int kflag,
double *pw, double *pdxn, double *pxmod)
{
/* 1.0 Newton iteration step */
// 1.0 Newton iteration step
double zdlx = pdx;
double zdlk = 0.0;
if ( kodd == 0 ) zdlk = 0.5*pfn[0];
double zdlxn = 0.0;
double zdlldn = 0.0;
size_t ik = 1;
if ( kflag == 0 )
if (kflag == 0)
{
for ( size_t jn = 2-kodd; jn <= kn; jn += 2 )
double zdlk = 0.5*pfn[0];
for (size_t jn = 2-kodd; jn <= kn; jn += 2)
{
/* normalised ordinary Legendre polynomial == \overbar{p_n}^0 */
// normalised ordinary Legendre polynomial == \overbar{p_n}^0
zdlk = zdlk + pfn[ik]*cos((double)(jn)*zdlx);
/* normalised derivative == d/d\theta(\overbar{p_n}^0) */
// normalised derivative == d/d\theta(\overbar{p_n}^0)
zdlldn = zdlldn - pfn[ik]*(double)(jn)*sin((double)(jn)*zdlx);
ik++;
}
/* Newton method */
// Newton method
double zdlmod = -(zdlk/zdlldn);