diff --git a/.gitattributes b/.gitattributes index 95bbfd180c1bde3128e58139b81612aed5a76e98..c7003d3f93b88ec2030ef9fbe12943a83445d108 100644 --- a/.gitattributes +++ b/.gitattributes @@ -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 diff --git a/ChangeLog b/ChangeLog index 29ce24251ce72a69ab23d83fa935b2d2f21b7d3e..cda59886040d6d0ca2e72f3267ec556ed1bd8c2b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -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* diff --git a/Makefile.am b/Makefile.am index d8be992feba6b9a17e1b7bfe6c1854bd7b69ca15..189502e99d1b5a701706edad427ca5b148fe5287 100644 --- a/Makefile.am +++ b/Makefile.am @@ -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 *~` diff --git a/Makefile.in b/Makefile.in index 267b2fea05613862da85d58f53fe8cc04024aba4..50c24b079aefa812c1eb800efe0b43e2b126d8ef 100644 --- a/Makefile.in +++ b/Makefile.in @@ -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,/[^/]*$$,,'`; \ diff --git a/README b/README index 435da850e4aa0cbffa063d36ed482f59914143e5..1298076617e34e495f6f8325ac94aa4bd9c27750 100644 --- a/README +++ b/README @@ -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. diff --git a/doc/tex/formats.tex b/doc/tex/formats.tex index 4e93ff12fe6b02eb40b90470aad67674cfec4531..d8d58e8cbb2c87bc1add3de0cebddc061f0f4e4a 100644 --- a/doc/tex/formats.tex +++ b/doc/tex/formats.tex @@ -1,7 +1,63 @@ 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} diff --git a/doc/tex/intro.tex b/doc/tex/intro.tex index 3cc53564171e21da59ed178159bdf42004f0fc19..7350b2516ca093b9bca209032b82e632ed7aca35 100644 --- a/doc/tex/intro.tex +++ b/doc/tex/intro.tex @@ -1,5 +1,10 @@ -% \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: diff --git a/src/cdi.h b/src/cdi.h index 9e3e29b45f448de752d301f4be65950da932adb3..26bee02a28e5ac190b34402288aad4085a554048 100644 --- a/src/cdi.h +++ b/src/cdi.h @@ -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 diff --git a/src/grib.h b/src/grib.h index 9ae16bf49cb7630453cd9b9c962a655ad31b510c..247d6245c2c2d012f5212c55aa32144caba7c13e 100644 --- a/src/grib.h +++ b/src/grib.h @@ -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 diff --git a/src/griblib.c b/src/griblib.c index 47824f66accd440df2a31e3fb7918440b5da4828..31cb69bbeaa35d19bb8284854c0532c37f28a18a 100644 --- a/src/griblib.c +++ b/src/griblib.c @@ -1,5 +1,5 @@ -/* 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"; diff --git a/src/stream_grb.c b/src/stream_grb.c index e77e9a6eca59a98a5b16966e2f77c93efc7ee012..259ab9d432bbce2ffcaaf8e2026fd9f121c86079 100644 --- a/src/stream_grb.c +++ b/src/stream_grb.c @@ -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;