From faef7bf4890d94e9e4c063eaff426282514a575d Mon Sep 17 00:00:00 2001 From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de> Date: Fri, 2 Jun 2006 11:19:37 +0000 Subject: [PATCH] stream_grb: update level type ISENTROPIC --- .gitattributes | 2 + ChangeLog | 1 + Makefile.am | 2 +- Makefile.in | 4 +- README | 2 +- doc/tex/formats.tex | 62 +++++++- doc/tex/intro.tex | 9 +- src/cdi.h | 2 +- src/grib.h | 1 + src/griblib.c | 373 ++++++++++++++++++++++++++++++++++++++++++-- src/stream_grb.c | 5 + 11 files changed, 438 insertions(+), 25 deletions(-) diff --git a/.gitattributes b/.gitattributes index 95bbfd180..c7003d3f9 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 29ce24251..cda598860 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 d8be992fe..189502e99 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 267b2fea0..50c24b079 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 435da850e..129807661 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 4e93ff12f..d8d58e8cb 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 3cc535641..7350b2516 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 9e3e29b45..26bee02a2 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 9ae16bf49..247d6245c 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 47824f66a..31cb69bbe 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 e77e9a6ec..259ab9d43 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; -- GitLab