Commit 39087bc0 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

New operator: genlevelbounds - Generate level bounds

parent 109e23c3
......@@ -4,6 +4,10 @@
* Version 1.6.9 released
* clipping: update to YAC version 1.0.3
2015-04-21 Uwe Schulzweida
* New operator: genlevelbounds - Generate level bounds
2015-04-17 Uwe Schulzweida
* Select: added parameter date
......
......@@ -11,6 +11,7 @@ Version 1.6.9 (26 September 2015):
* aexpr: Evaluate expressions and append results
* aexprf: Evaluate expression script and append results
* selzaxisname: Select z-axes by name
* genlevelbounds: Generate level bounds
Fixed bugs:
* seltimestep: abort if none of the selected timesteps are found
......
......@@ -147,6 +147,7 @@ Operator catalog:
Setgrid setgridtype Set grid type
Setgrid setgridarea Set grid cell area
Setzaxis setzaxis Set z-axis
Setzaxis genlevelbounds Generate level bounds
Setgatt setgatt Set global attribute
Setgatt setgatts Set global attributes
Invert invertlat Invert latitudes
......
......@@ -93,7 +93,7 @@ case "${HOSTNAME}" in
# x86_64-squeeze-x64-linux
thunder*)
CDOLIBS="--with-jasper=/sw/squeeze-x64/jasper-1.900.1-static \
--with-grib_api=/sw/squeeze-x64/grib_api-1.9.9-static \
--with-grib_api=/sw/squeeze-x64/grib_api-1.13.0-static-gccsys \
--with-netcdf=/sw/squeeze-x64/netcdf-4.2-static \
--with-hdf5=/sw/squeeze-x64/hdf5-1.8.8-static \
--with-szlib=/sw/squeeze-x64/szip-2.1 \
......@@ -102,7 +102,6 @@ case "${HOSTNAME}" in
if test "$COMP" = icc ; then
${CONFPATH}configure --prefix=$HOME/local --exec_prefix=$HOME/local/thunder \
--enable-all-static \
--with-fftw3 \
$CDOLIBS \
CC=icc CFLAGS="-g -Wall -O2 -qopt-report=5 -march=native -fp-model source -fast-transcendentals"
......
No preview for this file type
@BeginModule
@NewPage
@Name = Setzaxis
@Title = Set z-axis type
@Title = Set z-axis information
@Section = Modification
@Class = Modification
@Arguments = ifile ofile
@Operators = setzaxis
@Operators = setzaxis genlevelbounds
@BeginDescription
This module modifies the metadata of the vertical grid.
@EndDescription
@EndModule
......@@ -14,8 +18,17 @@
@Parameter = zaxis
@BeginDescription
This operator sets the z-axis description of all variables with the
same number of level as the new z-axis.
This operator sets the z-axis description of all variables with the same number of level as the new z-axis.
@EndDescription
@EndOperator
@BeginOperator_genlevelbounds
@Title = Generate level bounds
@Parameter = [zbot] [ztop]
@BeginDescription
Generates the layer bounds of the z-axis.
@EndDescription
@EndOperator
......@@ -23,4 +36,8 @@ same number of level as the new z-axis.
@BeginParameter
@Item = zaxis
STRING Z-axis description file or name of the target z-axis
@Item = zbot
FLOAT Specifying the bottom of the vertical column. Must have the same units as z-axis.
@Item = ztop
FLOAT Specifying the top of the vertical column. Must have the same units as z-axis.
@EndParameter
......@@ -26,46 +26,80 @@
#include "cdo_int.h"
#include "pstream.h"
void genLayerBounds(int nlev, double *levels, double *lbounds, double *ubounds);
int getkeyval_dp(const char *keyval, const char *key, double *val)
{
int status = 0;
size_t keylen = strlen(key);
if ( strncmp(keyval, key, keylen) == 0 )
{
const char *pkv = keyval+keylen;
if ( pkv[0] == '=' && pkv[1] != 0 )
{
*val = parameter2double(&pkv[1]);
status = 1;
}
else
{
cdoAbort("Syntax error for parameter %s!", keyval);
}
}
return status;
}
void *Setzaxis(void *argument)
{
int SETZAXIS;
int operatorID;
int streamID1, streamID2 = CDI_UNDEFID;
int nrecs;
int tsID, recID, varID, levelID;
int vlistID1, vlistID2;
int taxisID1, taxisID2;
int recID, varID, levelID;
int zaxisID1, zaxisID2 = -1;
int nzaxis, index;
int gridsize;
int nmiss;
int found;
double *array = NULL;
int lztop = FALSE, lzbot = FALSE;
double ztop = 0, zbot = 0;
cdoInitialize(argument);
SETZAXIS = cdoOperatorAdd("setzaxis", 0, 0, "zaxis description file");
int SETZAXIS = cdoOperatorAdd("setzaxis", 0, 0, "zaxis description file");
int GENLEVELBOUNDS = cdoOperatorAdd("genlevelbounds", 0, 0, NULL);
operatorID = cdoOperatorID();
operatorInputArg(cdoOperatorEnter(operatorID));
int operatorID = cdoOperatorID();
if ( operatorID == SETZAXIS )
{
operatorInputArg(cdoOperatorEnter(operatorID));
operatorCheckArgc(1);
zaxisID2 = cdoDefineZaxis(operatorArgv()[0]);
}
else if ( operatorID == GENLEVELBOUNDS )
{
unsigned npar = operatorArgc();
char **parnames = operatorArgv();
for ( unsigned i = 0; i < npar; i++ )
{
if ( cdoVerbose ) cdoPrint("keyval[%d]: %s", i+1, parnames[i]);
if ( !lzbot && getkeyval_dp(parnames[i], "zbot", &zbot) ) lzbot = TRUE;
else if ( !lztop && getkeyval_dp(parnames[i], "ztop", &ztop) ) lztop = TRUE;
else cdoAbort("Parameter >%s< unsupported! Supported parameter are: zbot, ztop", parnames[i]);
}
}
streamID1 = streamOpenRead(cdoStreamName(0));
int streamID1 = streamOpenRead(cdoStreamName(0));
vlistID1 = streamInqVlist(streamID1);
vlistID2 = vlistDuplicate(vlistID1);
int vlistID1 = streamInqVlist(streamID1);
int vlistID2 = vlistDuplicate(vlistID1);
taxisID1 = vlistInqTaxis(vlistID1);
taxisID2 = taxisDuplicate(taxisID1);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
if ( operatorID == SETZAXIS )
{
......@@ -83,14 +117,44 @@ void *Setzaxis(void *argument)
}
if ( ! found ) cdoWarning("No zaxis with %d levels found!", zaxisInqSize(zaxisID2));
}
else if ( operatorID == GENLEVELBOUNDS )
{
nzaxis = vlistNzaxis(vlistID1);
for ( index = 0; index < nzaxis; index++ )
{
zaxisID1 = vlistZaxis(vlistID1, index);
int nlev = zaxisInqSize(zaxisID1);
double *levels = (double *) malloc(nlev*sizeof(double));
double *lbounds = (double *) malloc(nlev*sizeof(double));
double *ubounds = (double *) malloc(nlev*sizeof(double));
if ( nlev > 1 )
{
zaxisInqLevels(zaxisID1, levels);
zaxisID2 = zaxisDuplicate(zaxisID1);
genLayerBounds(nlev, levels, lbounds, ubounds);
if ( lzbot ) lbounds[0] = zbot;
if ( lztop ) ubounds[nlev-1] = ztop;
zaxisDefLbounds(zaxisID2, lbounds);
zaxisDefUbounds(zaxisID2, ubounds);
vlistChangeZaxisIndex(vlistID2, index, zaxisID2);
}
free(levels);
free(lbounds);
free(ubounds);
}
}
streamDefVlist(streamID2, vlistID2);
gridsize = vlistGridsizeMax(vlistID1);
int gridsize = vlistGridsizeMax(vlistID1);
if ( vlistNumber(vlistID1) != CDI_REAL ) gridsize *= 2;
array = (double*) malloc(gridsize*sizeof(double));
double *array = (double*) malloc(gridsize*sizeof(double));
tsID = 0;
int tsID = 0;
while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
{
taxisCopyTimestep(taxisID2, taxisID1);
......
......@@ -73,6 +73,27 @@ void setSurfaceID(int vlistID, int surfID)
}
}
void genLayerBounds(int nlev, double *levels, double *lbounds, double *ubounds)
{
if ( nlev == 1 )
{
lbounds[0] = 0.;
ubounds[0] = 1.;
}
else
{
lbounds[0] = levels[0];
ubounds[nlev-1] = levels[nlev-1];
for ( int i = 0; i < nlev-1; ++i )
{
double bound = 0.5*(levels[i] + levels[i+1]);
lbounds[i+1] = bound;
ubounds[i] = bound;
}
}
}
static
int getLayerThickness(int genbounds, int index, int zaxisID, int nlev, double *thickness, double *weights)
{
......@@ -86,22 +107,7 @@ int getLayerThickness(int genbounds, int index, int zaxisID, int nlev, double *t
if ( genbounds )
{
status = 2;
if ( nlev == 1 )
{
lbounds[0] = 0.;
ubounds[0] = 1.;
}
else
{
lbounds[0] = levels[0];
ubounds[nlev-1] = levels[nlev-1];
for ( i = 0; i < nlev-1; ++i )
{
double bound = 0.5*(levels[i] + levels[i+1]);
lbounds[i+1] = bound;
ubounds[i] = bound;
}
}
genLayerBounds(nlev, levels, lbounds, ubounds);
}
else if ( zaxisInqLbounds(zaxisID, NULL) && zaxisInqUbounds(zaxisID, NULL) )
{
......
......@@ -405,7 +405,7 @@ void *Maggraph(void *argument);
#define SetrcanameOperators {"setrcaname"}
#define SettimeOperators {"setyear", "setmon", "setday", "setdate", "settime", "settunits", \
"settaxis", "setreftime", "setcalendar", "shifttime"}
#define SetzaxisOperators {"setzaxis"}
#define SetzaxisOperators {"setzaxis", "genlevelbounds"}
#define ShowinfoOperators {"showyear", "showmon", "showdate", "showtime", "showtimestamp", "showcode", "showunit", \
"showparam", "showname", "showstdname", "showlevel", "showltype", "showformat"}
#define SinfoOperators {"sinfo", "sinfop", "sinfon", "sinfoc", "seinfo", "seinfop", "seinfon", "seinfoc"}
......
......@@ -1040,17 +1040,25 @@ static char *SetgridHelp[] = {
static char *SetzaxisHelp[] = {
"NAME",
" setzaxis - Set z-axis type",
" setzaxis, genlevelbounds - Set z-axis information",
"",
"SYNOPSIS",
" setzaxis,zaxis ifile ofile",
" genlevelbounds[,zbot[,ztop]] ifile ofile",
"",
"DESCRIPTION",
" This operator sets the z-axis description of all variables with the",
" same number of level as the new z-axis.",
" This module modifies the metadata of the vertical grid.",
"",
"OPERATORS",
" setzaxis Set z-axis",
" This operator sets the z-axis description of all variables with the same number of level as the new z-axis.",
" genlevelbounds Generate level bounds",
" Generates the layer bounds of the z-axis.",
"",
"PARAMETER",
" zaxis STRING Z-axis description file or name of the target z-axis",
" zbot FLOAT Specifying the bottom of the vertical column. Must have the same units as z-axis. ",
" ztop FLOAT Specifying the top of the vertical column. Must have the same units as z-axis. ",
NULL
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment