Skip to content

GitLab

Commit faef7bf4 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

stream_grb: update level type ISENTROPIC

parent db2d7383
Hide whitespace changes
Inline Side-by-side
.gitattributes
View file @ faef7bf4
......@@ -18,6 +18,8 @@ config/missing -text
config/mkinstalldirs -text
/configure -text
/configure.ac -text
doc/cdi_cman.pdf -text
doc/cdi_fman.pdf -text
doc/tex/FUNCTIONS -text
doc/tex/Modules -text
doc/tex/bib.tex -text
......
ChangeLog
View file @ faef7bf4
......@@ -2,6 +2,7 @@
* using GRIB library version 1.0.0
* rename *New functions to *Create
* stream_grb: update level type ISENTROPIC (report: Ag Stephens)
* stream_ieg: multiply pressure levels with 100 (report: Holger Goettel)
* gridPrint: GME support
* cdi_error: change return type to char*
......
Makefile.am
View file @ faef7bf4
......@@ -2,6 +2,6 @@
#
SUBDIRS = src prog examples tests
#
EXTRA_DIST=config/default
EXTRA_DIST=config/default doc/cdi_cman.pdf doc/cdi_fman.pdf
#
CLEANFILES = `ls *~`
Makefile.in
View file @ faef7bf4
......@@ -84,7 +84,7 @@ install_sh = @install_sh@
#
SUBDIRS = src prog examples tests
#
EXTRA_DIST = config/default
EXTRA_DIST = config/default doc/cdi_cman.pdf doc/cdi_fman.pdf
#
CLEANFILES = `ls *~`
subdir = .
......@@ -235,7 +235,7 @@ distcleancheck_listfiles = find . -type f -print
distdir: $(DISTFILES)
$(am__remove_distdir)
mkdir $(distdir)
$(mkinstalldirs) $(distdir)/config
$(mkinstalldirs) $(distdir)/config $(distdir)/doc
@list='$(DISTFILES)'; for file in $$list; do \
if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \
dir=`echo "$$file" | sed -e 's,/[^/]*$$,,'`; \
......
README
View file @ faef7bf4
......@@ -5,7 +5,7 @@ CDI - Climate Data Interface
This package contains the source code for the CDI library from
the Max-Planck-Institute for Meteorology.
CDI is a C and Fortran Interface to access Climate model Data.
Supported file formats are GRIB, netCDF, SERVICE, EXTRA and IEG.
Supported data formats are GRIB, netCDF, SERVICE, EXTRA and IEG.
CDI is licensed under the GNU General Public License, version 2.
Read the file COPYING in the source distribution for details.
......
doc/tex/formats.tex
View file @ faef7bf4
Every input and output file is a collection of 2D or 3D variables
over an unlimited number of time steps.
\section{GRIB}
\section{GRIB edition 1}
\htmlref{GRIB}{GRIB} (GRIdded Binary) is a standard format designed by the
World Meteorological Organization (WMO) to support the efficient
transmission and storage of gridded meteorological data.
A GRIB record consists of a series of header sections, followed by
a bitstream of packed data representing one horizontal grid of data
values. The header sections are intended to fully describe the data
included in the bitstream, specifying information such as the
parameter, units, and precision of the data, the grid system and
level type on which the data is provided, and the date and time
for which the data are valid.
Non-numeric descriptors are enumerated in tables, such that a 1-byte
code in a header section refers to a unique description.
The WMO provides a standard set of enumerated parameter names and
level types, but the standard also allows for the definition
of locally used parameters and geometries. Any activity
that generates and distributes GRIB records must also make
their locally defined GRIB tables available to users.
\CDI does not support the full GRIB standard. The following
data representation and level types are implemented:
\begin{itemize}
\item Latitude/Longitude Grid
\item Gaussian Latitude/Longitude Grid
\item Spherical Harmonic Coefficients
\item Icosahedral-hexagonal GME Grid
\end{itemize}
\begin{itemize}
\item Surface level
\item Isobaric level
\item Height above ground
\item Hybrid level
%\item Layer between two hybrid levels
\item Depth below land surface
%\item Layer between two depths below land surface
\item Depth below sea level
\end{itemize}
%0 Latitude/Longitude Grid
%4 Gaussian Latitude/Longitude Grid
%50 Spherical Harmonic Coefficients
%192 Icosahedral-hexagonal GME Grid
%
%1 Surface level
%100 Isobaric level
%103 fixed height level
%105 Height above ground
%109 Hybrid level
%110 Layer between two hybrid levels
%111 Depth below land surface
%112 Layer between two depths below land surface
%160 Depth below sea level
\section{NetCDF}
......@@ -12,14 +68,14 @@ The netCDF library also defines a machine-independent format for
representing scientific data. Together, the interface, library, and
format support the creation, access, and sharing of scientific data.
\CDI supports only 2D, 3D and 4D arrays and the attributes must follow the
\CDI supports only 2D, 3D and 4D arrays and the attributes should follow the
\href{http://ftp.unidata.ucar.edu/software/netcdf/docs/conventions.html}
{GDT, COARDS or CF Conventions}.
NetCDF is an external library and not part of \CDI. To use netCDF with
\CDI the netCDF library must be installed before the configuration
of the \CDI library (see \htmlref{Build}{build}).
\subsection{ncdap}
%\subsection{ncdap}
\section{SERVICE}
......
doc/tex/intro.tex
View file @ faef7bf4
% \CDI is a C and Fortran Interface to access Climate model Data.
% Supported file formats are GRIB, netCDF, SERVICE, EXTRA and IEG.
\CDI is an Interface to access Climate model Data.
The interface is independent from a specific data format
and has a C and Fortran API.
\CDI was developed for a fast and machine independent access
to GRIB and netCDF datasets with the same interface.
The local data formats SERVICE, EXTRA and IEG are also supported.
% limitations:
......
src/cdi.h
View file @ faef7bf4
......@@ -51,7 +51,7 @@ extern "C" {
#define GRID_LONLAT 4
#define GRID_SPECTRAL 5
#define GRID_FOURIER 6
#define GRID_GME 7
#define GRID_GME 7 /* Icosahedral-hexagonal GME Grid */
#define GRID_TRAJECTORY 8
#define GRID_CELL 9
#define GRID_CURVILINEAR 10
......
src/grib.h
View file @ faef7bf4
......@@ -18,6 +18,7 @@
#define LTYPE_HYBRID_LAYER 110
#define LTYPE_LANDDEPTH 111
#define LTYPE_LANDDEPTH_LAYER 112
#define LTYPE_ISENTROPIC 113
#define LTYPE_SEADEPTH 160
#define LTYPE_99_MARGIN 1000
......
src/griblib.c
View file @ faef7bf4
/* Generated automatically from m214003 on Mon May 22 13:02:46 CEST 2006 */
/* Generated automatically from m214003 on Fri Jun 2 12:42:10 CEST 2006 */
#if defined (HAVE_CONFIG_H)
# include "config.h"
......@@ -80,8 +80,12 @@ void scatterComplex(double *fpdata, int pcStart, int truncation, int dimSP);
void scm0(double *pdl, double *pdr, double *pfl, double *pfr, int klg);
int rowina2(double *p, int ko, int ki, double *pw,
int kcode, double msval, int *kret);
int rowina3(double *p, int ko, int ki, double *pw,
int kcode, double msval, int *kret);
int qu2reg2(double *pfield, int *kpoint, int klat, int klon,
double *ztemp, double msval, int *kret);
int qu2reg3(double *pfield, int *kpoint, int klat, int klon,
double *ztemp, double msval, int *kret);
#if defined (INT32)
long packInt32(INT32 *up, char *cp, long bc, long tc);
......@@ -2772,14 +2776,14 @@ void encodePDS(GRIBPACK *lpds, int pdsLen, int *isec1)
Put1Byte(ISEC1_Parameter); /* 8 Parameter Code */
Put1Byte(ISEC1_LevelType); /* 9 Type of level */
if ( (ISEC1_LevelType != 20) &&
(ISEC1_LevelType != LTYPE_99) &&
(ISEC1_LevelType != LTYPE_ISOBARIC) &&
(ISEC1_LevelType != LTYPE_ALTITUDE) &&
(ISEC1_LevelType != LTYPE_HEIGHT) &&
(ISEC1_LevelType != LTYPE_SIGMA) &&
(ISEC1_LevelType != LTYPE_HYBRID) &&
(ISEC1_LevelType != LTYPE_LANDDEPTH) &&
(ISEC1_LevelType != 113) &&
(ISEC1_LevelType != LTYPE_99) &&
(ISEC1_LevelType != LTYPE_ISOBARIC) &&
(ISEC1_LevelType != LTYPE_ALTITUDE) &&
(ISEC1_LevelType != LTYPE_HEIGHT) &&
(ISEC1_LevelType != LTYPE_SIGMA) &&
(ISEC1_LevelType != LTYPE_HYBRID) &&
(ISEC1_LevelType != LTYPE_LANDDEPTH) &&
(ISEC1_LevelType != LTYPE_ISENTROPIC) &&
(ISEC1_LevelType != 115) &&
(ISEC1_LevelType != 117) &&
(ISEC1_LevelType != 125) &&
......@@ -3697,7 +3701,7 @@ static int decodePDS(unsigned char *pds, int *isec0, int *isec1)
(ISEC1_LevelType != LTYPE_SIGMA) &&
(ISEC1_LevelType != LTYPE_HYBRID) &&
(ISEC1_LevelType != LTYPE_LANDDEPTH) &&
(ISEC1_LevelType != 113) &&
(ISEC1_LevelType != LTYPE_ISENTROPIC) &&
(ISEC1_LevelType != 115) &&
(ISEC1_LevelType != 117) &&
(ISEC1_LevelType != 125) &&
......@@ -4561,7 +4565,8 @@ void gribDecode(int *isec0, int *isec1, int *isec2, double *fsec2, int *isec3,
ztemp = (double *) malloc(ISEC4_NumValues*sizeof(double));
if ( ztemp == NULL ) SysError(func, "No Memory!");
(void) qu2reg2(fsec4, ISEC2_RowLonPtr, nlat, nlon, ztemp, FSEC3_MissVal, iret);
/* (void) qu2reg2(fsec4, ISEC2_RowLonPtr, nlat, nlon, ztemp, FSEC3_MissVal, iret); */
(void) qu2reg3(fsec4, ISEC2_RowLonPtr, nlat, nlon, ztemp, FSEC3_MissVal, iret);
free(ztemp);
if ( bitmapSize > 0 )
......@@ -5084,6 +5089,7 @@ void scm0(double *pdl, double *pdr, double *pfl, double *pfr, int klg)
}
} /* scm0 */
int rowina2(double *p, int ko, int ki, double *pw,
int kcode, double msval, int *kret)
{
......@@ -5252,6 +5258,176 @@ L900:
return 0;
} /* rowina2 */
int rowina3(double *p, int ko, int ki, double *pw,
int kcode, double msval, int *kret)
{
/* System generated locals */
int pw_dim1, pw_offset, i_1;
/* Local variables */
double zwt1, zrdi, zpos;
int jl, ip;
double zdo, zwt;
/* Parameter adjustments */
--p;
pw_dim1 = ko + 3;
pw_offset = pw_dim1;
pw -= pw_offset;
/* **** ROWINA2 - Interpolation of row of values. */
/* Input Parameters. */
/* ----------------- */
/* P - Row of values to be interpolated. */
/* Dimension must be at least KO. */
/* KO - Number of values required. */
/* KI - Number of values in P on input. */
/* PW - Working array. */
/* Dimension must be at least (0:KO+2,3). */
/* KCODE - Interpolation required. */
/* 1 , linear. */
/* 3 , cubic. */
/* PMSVAL - Value used for missing data indicator. */
/* Output Parameters. */
/* ------------------ */
/* P - Now contains KO values. */
/* KRET - Return code */
/* 0, OK */
/* Non-zero, error */
/* Author. */
/* ------- */
/* J.D.Chambers ECMWF 22.07.94 */
/* ******************************** */
/* Section 1. Linear interpolation .. */
/* ******************************** */
*kret = 0;
if ( kcode == 1 )
{
/* Move input values to work array */
for ( jl = 1; jl <= ki; ++jl )
pw[jl + pw_dim1] = p[jl];
/* Arrange wrap-around value in work array */
pw[ki + 1 + pw_dim1] = p[1];
/* Set up constants to be used to figure out weighting for */
/* values in interpolation. */
zrdi = (double) ki;
zdo = 1.0 / (double) ko;
/* Loop through the output points */
for ( jl = 1; jl <= ko; ++jl )
{
/* Calculate weight from the start of row */
zpos = (jl - 1) * zdo;
zwt = zpos * zrdi;
/* Get the current array position(minus 1) from the weight - */
/* note the implicit truncation. */
ip = (int) zwt;
/* If the left value is missing, use the right value */
if (pw[ip + 1 + pw_dim1] == msval)
{
p[jl] = pw[ip + 2 + pw_dim1];
}
/* If the right value is missing, use the left value */
else if (pw[ip + 2 + pw_dim1] == msval)
{
p[jl] = pw[ip + 1 + pw_dim1];
}
/* If neither missing, interpolate ... */
else
{
/* Adjust the weight to range (0.0 to 1.0) */
zwt -= ip;
/* Interpolate using the weighted values on either side */
/* of the output point position */
p[jl] = (1.0 - zwt) * pw[ip + 1 + pw_dim1] +
zwt * pw[ip + 2 + pw_dim1];
}
}
/* ******************************* */
/* Section 2. Cubic interpolation .. */
/* ******************************* */
}
else if ( kcode == 3 )
{
i_1 = ki;
for ( jl = 1; jl <= i_1; ++jl )
{
if ( p[jl] == msval )
{
fprintf(stderr," ROWINA2: ");
fprintf(stderr," Cubic interpolation not supported");
fprintf(stderr," for fields containing missing data.\n");
*kret = 1;
goto L900;
}
pw[jl + pw_dim1] = p[jl];
}
pw[pw_dim1] = p[ki];
pw[ki + 1 + pw_dim1] = p[1];
pw[ki + 2 + pw_dim1] = p[2];
i_1 = ki;
for ( jl = 1; jl <= i_1; ++jl )
{
pw[jl + (pw_dim1 << 1)] =
- pw[jl - 1 + pw_dim1] / 3.0 -
pw[jl + pw_dim1] * 0.5 +
pw[jl + 1 + pw_dim1] - pw[jl + 2 + pw_dim1] / 6.0;
pw[jl + 1 + pw_dim1 * 3] =
pw[jl - 1 + pw_dim1] / 6.0 -
pw[jl + pw_dim1] +
pw[jl + 1 + pw_dim1] * 0.5 +
pw[jl + 2 + pw_dim1] / 3.0;
}
scm0(&pw[(pw_dim1 << 1) + 1], &pw[pw_dim1 * 3 + 2],
&pw[pw_dim1 + 1], &pw[pw_dim1 + 2], ki);
zrdi = (double) ki;
zdo = 1.0 / (double) ko;
for ( jl = 1; jl <= ko; ++jl )
{
zpos = (jl - 1) * zdo;
zwt = zpos * zrdi;
ip = (int) zwt + 1;
zwt = zwt + 1.0 - ip;
zwt1 = 1.0 - zwt;
p[jl] = ((3.0 - zwt1 * 2.0) * pw[ip + pw_dim1] +
zwt * pw[ip + (pw_dim1 << 1)]) * zwt1 * zwt1 +
((3.0 - zwt * 2.0) * pw[ip + 1 + pw_dim1] -
zwt1 * pw[ip + 1 + pw_dim1 * 3]) * zwt * zwt;
}
}
else
{
/* ************************************** */
/* Section 3. Invalid interpolation code .. */
/* ************************************** */
fprintf(stderr," ROWINA2:");
fprintf(stderr," Invalid interpolation code = %2d\n",kcode);
*kret = 2;
}
L900:
return 0;
} /* rowina3 */
int qu2reg2(double *pfield, int *kpoint, int klat, int klon,
double *ztemp, double msval, int *kret)
{
......@@ -5264,13 +5440,16 @@ int qu2reg2(double *pfield, int *kpoint, int klat, int klon,
/* Local variables */
int ilii, ilio, icode;
double *zline = NULL;
double zwork[1929]; /* was [643][3] */
double *zwork = NULL;
int iregno, iquano, j210, j220, j230, j240, j225;
zline = (double *) malloc(klon*sizeof(double));
zline = (double *) malloc(2*klon*sizeof(double));
if ( zline == NULL ) SysError(func, "No Memory!");
zwork = (double *) malloc(3*(2*klon+3)*sizeof(double));
if ( zwork == NULL ) SysError(func, "No Memory!");
/* Parameter adjustments */
--pfield;
--kpoint;
......@@ -5406,9 +5585,173 @@ int qu2reg2(double *pfield, int *kpoint, int klat, int klon,
L900:
free(zline);
free(zwork);
return 0;
} /* qu2reg2 */
int qu2reg3(double *pfield, int *kpoint, int klat, int klon,
double *ztemp, double msval, int *kret)
{
static char func[] = "qu2reg2";
/* System generated locals */
int i_1, i_2;
int kcode = 1;
/* Local variables */
int ilii, ilio, icode;
double *zline = NULL;
double *zwork = NULL;
int iregno, iquano, j210, j220, j230, j240, j225;
zline = (double *) malloc(2*klon*sizeof(double));
if ( zline == NULL ) SysError(func, "No Memory!");
zwork = (double *) malloc(3*(2*klon+3)*sizeof(double));
if ( zwork == NULL ) SysError(func, "No Memory!");
/* Parameter adjustments */
--pfield;
--kpoint;
/* **** QU2REG - Convert quasi-regular grid data to regular. */
/* Input Parameters. */
/* ----------------- */
/* PFIELD - Array containing quasi-regular grid */
/* data. */
/* KPOINT - Array containing list of the number of */
/* points on each latitude (or longitude) of */
/* the quasi-regular grid. */
/* KLAT - Number of latitude lines */
/* KLON - Number of longitude lines */
/* KCODE - Interpolation required. */
/* 1 , linear - data quasi-regular on */
/* latitude lines. */
/* 3 , cubic - data quasi-regular on */
/* latitude lines. */
/* 11, linear - data quasi-regular on */
/* longitude lines. */
/* 13, cubic - data quasi-regular on */
/* longitude lines. */
/* PMSVAL - Value used for missing data indicator. */
/* Output Parameters. */
/* ------------------ */
/* KRET - return code */
/* 0 = OK */
/* non-zero indicates fatal error */
/* PFIELD - Array containing regular grid data. */
/* Author. */
/* ------- */
/* J.D.Chambers ECMWF 22.07.94 */
/* J.D.Chambers ECMWF 13.09.94 */
/* Add return code KRET and remove calls to ABORT. */
/* ------------------------------ */
/* Section 1. Set initial values. */
/* ------------------------------ */
*kret = 0;
/* Check input parameters. */
if (kcode != 1 && kcode != 3 && kcode != 11 && kcode != 13) {
fprintf(stderr," QU2REG :");
fprintf(stderr," Invalid interpolation type code = %2d\n",kcode);
*kret = 1;
goto L900;
}
/* Set array indices to 0. */
ilii = 0;
ilio = 0;
/* Establish values of loop parameters. */
if (kcode > 10) {
/* Quasi-regular along longitude lines. */
iquano = klon;
iregno = klat;
icode = kcode - 10;
} else {
/* Quasi-regular along latitude lines. */
iquano = klat;
iregno = klon;
icode = kcode;
}
/* -------------------------------------------------------- */
/** Section 2. Interpolate field from quasi to regular grid. */
/* -------------------------------------------------------- */
i_1 = iquano;
for (j230 = 1; j230 <= i_1; ++j230) {
if (iregno != kpoint[j230]) {
/* Line contains less values than required,so */
/* extract quasi-regular grid values for a line */
i_2 = kpoint[j230];
for (j210 = 1; j210 <= i_2; ++j210) {
++ilii;
zline[j210 - 1] = pfield[ilii];
}
/* and interpolate this line. */
rowina3(zline, iregno, kpoint[j230], zwork, icode, msval, kret);
if (*kret != 0) goto L900;
/* Add regular grid values for this line to the
temporary array. */
i_2 = iregno;
for (j220 = 1; j220 <= i_2; ++j220) {
++ilio;
ztemp[ilio - 1] = zline[j220 - 1];
}
} else {
/* Line contains the required number of values, so add */
/* this line to the temporary array. */
i_2 = iregno;
for (j225 = 1; j225 <= i_2; ++j225) {
++ilio;
++ilii;
ztemp[ilio - 1] = pfield[ilii];
}
}
}
/* Copy temporary array to user array. */
i_1 = klon * klat;
for (j240 = 1; j240 <= i_1; ++j240) {
pfield[j240] = ztemp[j240 - 1];
}
/* -------------------------------------------------------- */
/* Section 9. Return to calling routine. Format statements. */
/* -------------------------------------------------------- */
L900:
free(zline);
free(zwork);
return 0;
} /* qu2reg2_ */
} /* qu2reg3 */
FILE *grprsm = NULL;
......@@ -7396,7 +7739,7 @@ int gribUnzip(unsigned char *dbuf, long dbufsize, unsigned char *sbuf, long sbu
return (gribLen);
}
static const char grb_libvers[] = "1.0.0" " of ""May 22 2006"" ""13:02:46";
static const char grb_libvers[] = "1.0.0" " of ""Jun 2 2006"" ""12:42:10";
......
src/stream_grb.c
View file @ faef7bf4
......@@ -315,6 +315,11 @@ int gribGetZaxisType(int grbleveltype)
leveltype = ZAXIS_DEPTH_BELOW_LAND;
break;
}
case LTYPE_ISENTROPIC:
{
leveltype = ZAXIS_ISENTROPIC;
break;
}
case LTYPE_SEADEPTH:
{
leveltype = ZAXIS_DEPTH_BELOW_SEA;
......