Commit 1c8ee1da authored by Ralf Mueller's avatar Ralf Mueller
Browse files

Added operators 'stdatm' and 'intlevel3d' (closes #1192)

* stdatm creates vertical single gridpoint values for pressure and temperature for a given list of zlevels
* intlevel3d interpolates 3d variables from/to 3d vertical coordinates like height fields
parent 5bd8d23d
......@@ -144,6 +144,7 @@ doc/tex/mod/Info -text
doc/tex/mod/Input -text
doc/tex/mod/Intgrid -text
doc/tex/mod/Intlevel -text
doc/tex/mod/Intlevel3d -text
doc/tex/mod/Inttime -text
doc/tex/mod/Intyear -text
doc/tex/mod/Invert -text
......@@ -315,6 +316,7 @@ src/Input.c -text
src/Intgrid.c -text
src/Intgridtraj.c -text
src/Intlevel.c -text
src/Intlevel3d.c -text
src/Intntime.c -text
src/Inttime.c -text
src/Intyear.c -text
......
......@@ -3,6 +3,11 @@
* using CDI library version 1.5.1
* Version 1.5.1 released
2011-07-11 Ralf Mueller <ralf.mueller@zmaw.de>
* Added operator stdatm: standard atmosphere values for T,PS
* Added operator intlevel3d: vertical interpolation to/from 3d vertical coordinates
2011-07-04 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* Exprf: wrong result for missing values != (double) -9.e33 (bug fix) [report: Michael Boettinger]
......
......@@ -398,6 +398,8 @@ Operator catalog:
Vertint ml2pl Model to pressure level interpolation
Vertint ml2hl Model to height level interpolation
Intlevel intlevel Linear level interpolation
Intlevel3d intlevel3d Linear level interpolation onto a 3d vertical coordinate
Intlevel3d intlevelx3d like intlevel3d but with extrapolation
Inttime inttime Interpolation between time steps
Inttime intntime Interpolation between time steps
Intyear intyear Interpolation between two years
......
......@@ -13,8 +13,8 @@
m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.66],,
[m4_warning([this file was generated for autoconf 2.66.
m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.68],,
[m4_warning([this file was generated for autoconf 2.68.
You have another version of autoconf. It may work, but is not guaranteed to.
If you have problems, you may need to regenerate the build system entirely.
To do so, use the procedure documented by the package, typically `autoreconf'.])])
......
This diff is collapsed.
......@@ -198,7 +198,9 @@ inputext -inputext \
inputsrv -inputsrv \
int -int \
intlevel -intlevel \
intlevel3d -intlevel3d \
intlevelx -intlevelx \
intlevelx3d -intlevelx3d \
intntime -intntime \
inttime -inttime \
intyear -intyear \
......
......@@ -198,7 +198,9 @@ inputext \
inputsrv \
int \
intlevel \
intlevel3d \
intlevelx \
intlevelx3d \
intntime \
inttime \
intyear \
......
......@@ -198,7 +198,9 @@ inputext -inputext \
inputsrv -inputsrv \
int -int \
intlevel -intlevel \
intlevel3d -intlevel3d \
intlevelx -intlevelx \
intlevelx3d -intlevelx3d \
intntime -intntime \
inttime -inttime \
intyear -intyear \
......
......@@ -87,6 +87,7 @@ Remap Interpolation
Remapeta Interpolation
Vertint Interpolation
Intlevel Interpolation
Intlevel3d Interpolation
Inttime Interpolation
Intyear Interpolation
Spectral Transformation
......
@BeginModule
@NewPage
@Name = Intlevel3d
@Title = Linear level interpolation from/to 3d vertical coordinates
@Section = Interpolation
@Class = Interpolation
@Arguments = ifile1 ifile2 ofile
@Operators = intlevel3d intlevelx3d
@BeginDescription
This operator performs a linear vertical interpolation of 3D variables fields
with given 3D vertical coordinates.
@EndDescription
@EndModule
@BeginOperator_intlevel3d
@Title = Linear level interpolation onto a 3d vertical coordinate
@Parameter = icoordinate
@BeginDescription
@IfMan
To interpolate 3D variables from one set of 3d height levels
into another one where
- icoordinate contains a single 3d variable,
which represents the input 3d vert. coordinate
- ifile1 contains the source data,
which the vertical coordinate from icoordinate belongs to
- ifile2 only contains the target 3d height levels
call
cdo intlevel3,icoordinate ifile1 ifile2 ofile
@EndifMan
@EndDescription
@EndOperator
@BeginOperator_intlevelx3d
@Title = like intlevel3d but with extrapolation
@Parameter = icoordinate
@EndOperator
@BeginParameter
@Item = icoordinate
STRING filename for vertical source coordinates variable
@Item = ifile2
STRING target vertical coordinate field (intlevel3d only)
@EndParameter
@IfDoc
@BeginExample
To interpolate 3D variables from one set of 3d height levels into another one where
@BeginItemize
@Item
@file{icoordinate} contains a single 3d variable, which represents the input 3d vert. coordinate
@Item
@file{ifile1} contains the source data, which the vertical coordinate from icoordinate belongs to
@Item
@file{ifile2} only contains the target 3d height levels
@EndItemize
@BeginVerbatim
cdo intlevel3d,icoordinate ifile1 ifile2 ofile
@EndVerbatim
@EndExample
@EndifDoc
/*
This file is part of CDO. CDO is a collection of Operators to
manipulate and analyse Climate model Data.
Copyright (C) 2003-2011 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
See COPYING file for copying and redistribution conditions.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
*/
/*
This module contains the following operators:
Intlevel intlevel3d Linear level interpolation on a 3d vertical coordinates variable
Intlevel intlevelx3d Linear level interpolation on a 3d vertical coordinates variable with extrapolation
*/
#include <ctype.h>
#include <cdi.h>
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "list.h"
/*
* 3d vertical interpolation routine (see interp_lev() in src/Intlevel.c)
*/
static
void interp_lev3d(int gridsize, double missval, double *vardata1, double *vardata2,
int nlev2, int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int i, ilev;
int idx1, idx2;
double wgt1, wgt2;
double w1, w2;
double var1L1, var1L2, *var2;
for ( ilev = 0; ilev < nlev2; ilev++ )
{
var2 = vardata2+(ilev*gridsize);
for ( i = 0; i < gridsize; i++ )
{
idx1 = lev_idx1[ilev*gridsize+i];
idx2 = lev_idx2[ilev*gridsize+i];
wgt1 = lev_wgt1[ilev*gridsize+i];
wgt2 = lev_wgt2[ilev*gridsize+i];
/* upper/lower values from input field */
var1L1 = *(vardata1+idx1);
var1L2 = *(vardata1+idx2);
/* if (cdoVerbose) printf("i:%d level %d: idx1=%d idx2=%d (ilev*gridsize+i:%d) wgt1=%g wgt2=%g var1L1:%g var1L2:%g ",
* i, ilev, idx1, idx2, ilev*gridsize+i, wgt1, wgt2, var1L1, var1L2);
*/
w1 = wgt1;
w2 = wgt2;
if ( DBL_IS_EQUAL(var1L1, missval) ) w1 = 0;
if ( DBL_IS_EQUAL(var1L2, missval) ) w2 = 0;
if ( IS_EQUAL(w1, 0) && IS_EQUAL(w2, 0) )
{
var2[i] = missval;
}
else if ( IS_EQUAL(w1, 0) )
{
if ( w2 >= 0.5 )
var2[i] = var1L2;
else
var2[i] = missval;
}
else if ( IS_EQUAL(w2, 0) )
{
if ( w1 >= 0.5 )
var2[i] = var1L1;
else
var2[i] = missval;
}
else
{
var2[i] = var1L1*w1 + var1L2*w2;
}
}
}
}
/*
* Create weights for the 3d vertical coordinate
*
* The resulting index sets lev_idx1 and lev_idx2 contain absolute numbers,i.e.
* wrt. the given gridsize. They can directly be used to read values from 3d
* data fields.
*
* 3d version of gen_weights() (src/Intlevel.c)
*/
static
void gen_weights3d(int expol, int nlev1, int gridsize, double *lev1, int nlev2, double *lev2,
int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
int i,i1, i2;
double val1, val2 = 0;
int idx1 = 0, idx2 = 0;
for ( i = 0; i < gridsize; i++ )
{
for ( i2 = 0; i2 < nlev2; i2++ )
{
/* Because 2 levels were added to the source vertical coordinate (one on
* top, one at the bottom), its loop starts at 1 */
for ( i1 = 1; i1 < nlev1; i1++ )
{
if ( lev1[(i1-1)*gridsize+i] < lev1[i1*gridsize+i] )
{
idx1 = (i1 - 1)*gridsize+i;
idx2 = i1*gridsize+i;
}
else
{
idx1 = i1*gridsize+i;
idx2 = (i1-1)*gridsize+i;
}
val1 = lev1[idx1];
val2 = lev1[idx2];
if ( lev2[i2*gridsize+i] > val1 && lev2[i2*gridsize+i] <= val2 ) break;
}
if ( i1 == nlev1 )
{
if ( expol )
cdoAbort("Level %g at index %d not found! Use extrapolation", lev2[i2*gridsize],i2);
else
cdoAbort("Level %g at index %d not found!");
}
if ( i1-1 == 0 ) /* destination levels ios not covert by the first two input z levels */
{
lev_idx1[i2*gridsize+i] = gridsize+i;
lev_idx2[i2*gridsize+i] = gridsize+i;
lev_wgt1[i2*gridsize+i] = 0;
if ( expol || IS_EQUAL(lev2[i2*gridsize+i], val2) )
lev_wgt2[i2*gridsize+i] = 1;
else
lev_wgt2[i2*gridsize+i] = 0;
}
else if ( i1 == nlev1-1 ) /* destination level is beyond the last value of the input z field */
{
lev_idx1[i2*gridsize+i] = (nlev1-2)*gridsize+i;
lev_idx2[i2*gridsize+i] = (nlev1-2)*gridsize+i;
if ( expol || IS_EQUAL(lev2[i2*gridsize+i], val2) )
lev_wgt1[i2*gridsize+i] = 1;
else
lev_wgt1[i2*gridsize+i] = 0;
lev_wgt2[i2*gridsize+i] = 0;
}
else /* target z values has two bounday values in input z field */
{
lev_idx1[i2*gridsize+i] = idx1;
lev_idx2[i2*gridsize+i] = idx2;
lev_wgt1[i2*gridsize+i] = (lev1[idx2] - lev2[i2*gridsize+i]) / (lev1[idx2] - lev1[idx1]);
lev_wgt2[i2*gridsize+i] = (lev2[i2*gridsize+i] - lev1[idx1]) / (lev1[idx2] - lev1[idx1]);
}
/* if (cdoVerbose)
* {
* printf("i:%d i2:%d\ti2*gridsize+i:%d\tlev2[i2*gridsize+i]:%g\tidx1:%d\tidx2:%d\tlev1[idx1]:%g\tlev1[idx2]:%g\t",
* i, i2, i2*gridsize+i, lev2[i2*gridsize+i], idx1, idx2, lev1[idx1], lev1[idx2]);
* printf("\tlev_wgt1:%g\tlev_wgt2:%g\n", lev_wgt1[i2*gridsize+i], lev_wgt2[i2*gridsize+i]);
* }
*/
/* backshift of the indices because if the two additional levels in input vertical coordinate */
lev_idx1[i2*gridsize+i] -= gridsize;
lev_idx2[i2*gridsize+i] -= gridsize;
}
}
}
void *Intlevel3d(void *argument)
{
int INTLEVEL, INTLEVELX, INTLEVEL3D, INTLEVELX3D;
int operatorID;
int streamID0, streamID1, streamID2,streamID3;
int vlistID0, vlistID1, vlistID2, vlistID3;
int gridsize,gridSize,gridsizei,gridsizeo;
int recID, nrecs;
int i, offset;
int tsID, varID, levelID;
int nvars,nvct;
int nzaxis;
int nmiss;
int zaxisID1 = -1, zaxisID3;
int gridID, zaxisID;
int nlevi, nlevo, nlevel = 0, maxlev;
int lup, ldown;
int **varnmiss = NULL;
int *varinterp = NULL;
int *vars = NULL;
int expol = FALSE;
double missval;
double *lev1 = NULL, *lev2 = NULL;
double *single1, *single2;
double **vardata1 = NULL, **vardata2 = NULL;
int taxisID1, taxisID2, taxisID3;
int *lev_idx1, *lev_idx2;
double *lev_wgt1, *lev_wgt2;
double *zlevels_in, *zlevels_out;
int zlevels_in_miss, zlevels_out_miss;
char varname[10];
cdoInitialize(argument);
INTLEVEL3D = cdoOperatorAdd("intlevel3d", 0, 0, NULL);
INTLEVELX3D = cdoOperatorAdd("intlevelx3d", 0, 0, NULL);
operatorID = cdoOperatorID();
if ( cdoVerbose )
cdoPrint("Called '%s'\n",__func__);
if ( operatorID == INTLEVELX3D ) expol = TRUE;
operatorInputArg("icoordinate");
/* Read filename from Parameter */
operatorInputArg("grid description file or name, remap file (SCRIP netCDF)");
operatorCheckArgc(1);
streamID0 = streamOpenRead(operatorArgv()[0]); /* 3d vertical input coordinate */
streamID1 = streamOpenRead(cdoStreamName(0)); /* input data */
streamID2 = streamOpenRead(cdoStreamName(1)); /* 3d target vertical coordinate */
streamID3 = streamOpenWrite(cdoStreamName(2),cdoFiletype()); /* output stream */
vlistID0 = streamInqVlist(streamID0);
vlistID1 = streamInqVlist(streamID1); taxisID1 = vlistInqTaxis(vlistID1);
vlistID2 = streamInqVlist(streamID2);
vlistID3 = vlistDuplicate(vlistID1); taxisID3 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID3, taxisID1);
/*
* Read 3d input coordinate (streamID0)
* * two additional levels are added (top + bottom) for later extrapolation checking
*/
{
nvars = vlistNvars(vlistID0);
if (nvars != 1)
cdoAbort("Only one single variable is allowed!");
gridID = vlistInqVarGrid(vlistID0, 0);
zaxisID = vlistInqVarZaxis(vlistID0, 0);
gridsize = gridInqSize(gridID);
nlevel = zaxisInqSize(zaxisID);
zlevels_in = (double *) malloc(gridsize*(nlevel+2)*sizeof(double*));
nlevi = nlevel; /* number of input levels for later use */
gridsizei = gridsize; /* horizontal gridsize of input z coordinate */
nrecs = streamInqTimestep(streamID0, 0);
if (cdoVerbose)
cdoPrint("%d records input 3d vertical height\n",nrecs);
{
for ( recID = 0; recID < nrecs; recID++ )
{
streamInqRecord(streamID0, &varID, &levelID);
gridsize = gridInqSize(vlistInqVarGrid(vlistID0, varID));
offset = gridsize + gridsize*levelID;
single1 = zlevels_in + offset;
streamReadRecord(streamID0, single1, &zlevels_in_miss);
}
}
}
/*
* Read 3d output coordinate (streamID2)
*/
{
nvars = vlistNvars(vlistID2);
if (nvars != 1) cdoAbort("Only one single variable is allowed!");
gridID = vlistInqVarGrid(vlistID2, varID);
zaxisID = vlistInqVarZaxis(vlistID2, varID);
gridsize = gridInqSize(gridID);
nlevel = zaxisInqSize(zaxisID);
zlevels_out = (double *) malloc(gridsize*nlevel*sizeof(double));
nlevo = nlevel; /* number of output levels for later use */
gridsizeo = gridsize;/* horizontal gridsize of output z coordinate */
nrecs = streamInqTimestep(streamID2, 0);
if (cdoVerbose)
cdoPrint("%d records target 3d vertical height and gridsize %d\n",nrecs,gridsize);
for ( recID = 0; recID < nrecs; recID++ )
{
streamInqRecord(streamID2, &varID, &levelID);
gridsize = gridInqSize(vlistInqVarGrid(vlistID2, varID));
offset = gridsize*levelID;
single1 = zlevels_out + offset;
streamReadRecord(streamID2, single1, &zlevels_out_miss);
}
}
/* Missing values are not allowed for coordinate variables */
if ( 0 != zlevels_in_miss )
cdoAbort("Input vertical coordinate variables are not allowd to contian missing values.");
else
cdoPrint("Input vertical coordinate has no missing values.");
if ( 0 != zlevels_out_miss )
cdoAbort("Output vertical coordinate variables are not allowd to contian missing values.");
else
cdoPrint("Output vertical coordinate has no missing values.");
/*
* gridsize of input and output vertical coordinate must be equal
* (later use of gridsizeo ONLY)
*/
if ( gridsizei != gridsizeo )
cdoAbort("Input and output vertical coordinate must have the same gridsize");
gridSize = gridsizeo;
/*
* Check for the correct vertical axis in the input: Variables with the same
* number of levels as the input vertical levels from operators parameter
* (streamID0). Variables with a different z-axis should be copied into output.
*/
nzaxis = vlistNzaxis(vlistID1);
for ( i = 0; i < nzaxis; i++ )
{
zaxisID = vlistZaxis(vlistID1, i);
nlevel = zaxisInqSize(zaxisID);
if ( nlevel == nlevi )
{
zaxisID1 = zaxisID;
break;
}
}
if ( i == nzaxis ) cdoAbort("No processable variable found!");
/*
* Check monotony of vertical levels
*/
lup = FALSE;
ldown = FALSE;
lev1 = zlevels_in + gridSize;
for ( i = 0; i < nlevi-1; i++ )
{
if ( i == 1 )
{
if ( lev1[(i+1)*gridSize] > lev1[i*gridSize] )
lup = TRUE;
else if ( lev1[(i+1)*gridSize] < lev1[i*gridSize] )
ldown = TRUE;
}
else
{
if ( lup )
{
if ( !(lev1[(i+1)*gridSize] > lev1[i*gridSize]) ) lup = FALSE;
}
else if ( ldown )
{
if ( !(lev1[(i+1)*gridSize] < lev1[i*gridSize]) ) ldown = FALSE;
}
}
}
/* Add artificial values for intication of extrapolation areas (lowermost + upmost levels) */
if ( lup )
{
for ( i = 0; i < gridSize ;i++)
{
zlevels_in[i] = -1.e33;
zlevels_in[(nlevi+1)*gridSize + i] = 1.e33;
}
}
else if ( ldown )
{
for ( i = 0; i < gridSize ;i++)
{
zlevels_in[i] = 1.e33;
zlevels_in[(nlevi+1)*gridSize + i] = -1.e33;
}
}
else
cdoWarning("Non monotonic zaxis!");
/*
* Create weights for later interpolation - assumption: input vertical correct is constant in time
*/
lev_idx1 = (int *) malloc(nlevo*gridSize*sizeof(int));
lev_idx2 = (int *) malloc(nlevo*gridSize*sizeof(int));
lev_wgt1 = (double *) malloc(nlevo*gridSize*sizeof(double));
lev_wgt2 = (double *) malloc(nlevo*gridSize*sizeof(double));
gen_weights3d(expol, nlevi+2, gridSize, zlevels_in, nlevo, zlevels_out, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
/*
* Copy z-axis information to output z-axis
*/
zaxisID3 = zaxisCreate(zaxisInqType(zaxisID1), nlevo);
lev2 = (double *) malloc(nlevo*sizeof(double));
/* fill values with its indices */
for (i=0;i<nlevo;i++)
lev2[i] = (double) i;
zaxisDefLevels(zaxisID3, lev2);
zaxisDefName(zaxisID3, "lev");
/* copy VCT from input vlistID1 to output vlistID3 if there is one */
nvct = zaxisInqVctSize(zaxisID1);
if ( nvct > 0 ) zaxisDefVct(zaxisID3,zaxisInqVctSize(zaxisID1), zaxisInqVctPtr(zaxisID1));
for ( i = 0; i < nzaxis; i++ )
if ( zaxisID1 == vlistZaxis(vlistID1, i) )
vlistChangeZaxisIndex(vlistID3, i, zaxisID3);
/* add the vertical output field to the output stream */
int oz3dvarID = vlistDefVar(vlistID3,0,zaxisID3,TIME_VARIABLE);
{
char str[256];
str[0] = 0;
vlistInqVarName(vlistID2,0,str);
vlistDefVarName(vlistID3,oz3dvarID,str);
str[0] = 0;
vlistInqVarLongname(vlistID2,0,str);
if ( str[0] ) vlistDefVarLongname(vlistID3,zaxisID1, str);
str[0] = 0;
vlistInqVarUnits(vlistID2,0, str);
if ( str[0] ) vlistDefVarUnits(vlistID3,oz3dvarID, str);
}
streamDefVlist(streamID3, vlistID3);
maxlev = nlevi > nlevo ? nlevi : nlevo;
nvars = vlistNvars(vlistID1);
vars = (int *) malloc(nvars*sizeof(int));
vardata1 = (double **) malloc(nvars*sizeof(double*)); /* input */
vardata2 = (double **) malloc(nvars*sizeof(double*)); /* output */
varnmiss = (int **) malloc(nvars*sizeof(int*)); /* can for missing values of arbitrary variables */
varinterp = (int *) malloc(nvars*sizeof(int)); /* marker for variables to be interpolated */
/* by default no variable should be interpolated */
for (i = 0; i < nvars;i++)
varinterp[varID] = FALSE;
for ( varID = 0; varID < nvars; varID++ )
{
gridID = vlistInqVarGrid(vlistID1, varID);
zaxisID = vlistInqVarZaxis(vlistID1, varID);
gridsize = gridInqSize(gridID);
nlevel = zaxisInqSize(zaxisID);
vlistInqVarName(vlistID1, varID, varname);
vardata1[varID] = (double *) malloc(gridsize*nlevel*sizeof(double));
/* variabls for interpolation:
* * have the required vertical axis, i.e. the correct number of levels (nlevi)
* * have the same number of horizontal grid points (i.e. same gridSize) like the two vertical coordinates
* * are NOT the output vertical coordinates itself
*/
if ( zaxisID == zaxisID1 && varID != oz3dvarID && gridsize ==