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

cdiwrite update

parent ef13de48
......@@ -15,256 +15,91 @@
GNU General Public License for more details.
*/
/*
This module contains the following operators:
Vargen const Create a constant field
Vargen random Field with random values
Vargen stdatm Field values for pressure and temperature for
the standard atmosphere
*/
#include <cdi.h>
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
#include "list.h"
#if defined (__GNUC__)
#if __GNUC__ > 2
# define WITH_DATA 1
#endif
#endif
#if defined(WITH_DATA)
static double etopo_scale = 3;
static double etopo_offset = 11000;
static const unsigned short etopo[] = {
#include "etopo.h"
};
static double temp_scale = 500;
static double temp_offset = -220;
static const unsigned short temp[] = {
#include "temp.h"
};
static double mask_scale = 1;
static double mask_offset = 0;
static const unsigned short mask[] = {
#include "mask.h"
};
#endif
/* some Constants for creating temperatur and pressure for the standard atmosphere */
#define T_ZERO (213.0)
#define T_DELTA (75.0)
#define SCALEHEIGHT (10000.0) /* [m] */
#define P_ZERO (1013.25) /* surface pressure [hPa] */
#define C_EARTH_GRAV (9.80665)
#define C_R (287.05) /* specific gas constant for air */
static double TMP4PRESSURE = (C_EARTH_GRAV*SCALEHEIGHT)/(C_R*T_ZERO);
static double
std_atm_temperatur(double height)
{
/*
Compute the temperatur for the given height (in meters) according to the
solution of the hydrostatic atmosphere
*/
return (T_ZERO + T_DELTA * exp((-1)*(height/SCALEHEIGHT)));
}
static double
std_atm_pressure(double height)
{
/*
Compute the pressure for the given height (in meters) according to the
solution of the hydrostatic atmosphere
*/
return (P_ZERO * exp((-1)*TMP4PRESSURE*log((exp(height/SCALEHEIGHT)*T_ZERO + T_DELTA)/(T_ZERO + T_DELTA))));
}
void *CDIwrite(void *argument)
{
int RANDOM, CONST, FOR, TOPO, TEMP, MASK, STDATM;
int nvars = 1, nlevs = 0, ntimesteps = 0;
int operatorID;
int streamID;
int nrecs,nvars, ntimesteps, nlevels = 1;
int tsID, recID, varID, varID2, levelID;
int tsID, varID, levelID;
int gridsize, i;
int rval, rstart, rinc;
int vlistID;
int gridID = -1, zaxisID, taxisID;
int vdate, vtime, julday;
unsigned int seed = 1;
const char *gridfile;
double rval, rstart = 0, rstop = 0, rinc = 0;
double rconst = 0;
double *array, *levels;
LIST *flist = listNew(FLT_LIST);
double *levels = NULL;
double ***vars = NULL;
cdoInitialize(argument);
RANDOM = cdoOperatorAdd("random", 0, 0, "grid description file or name, <seed>");
CONST = cdoOperatorAdd("const", 0, 0, "constant value, grid description file or name");
FOR = cdoOperatorAdd("for", 0, 0, "start, end, <increment>");
TOPO = cdoOperatorAdd("topo", 0, 0, NULL);
TEMP = cdoOperatorAdd("temp", 0, 0, NULL);
MASK = cdoOperatorAdd("mask", 0, 0, NULL);
STDATM = cdoOperatorAdd("stdatm", 0, 0, "levels");
operatorInputArg("grid, <nlevs, <ntimesteps, <nvars>>>");
if ( operatorArgc() < 1 ) cdoAbort("Too few arguments!");
if ( operatorArgc() > 4 ) cdoAbort("Too many arguments!");
operatorID = cdoOperatorID();
if ( operatorID == RANDOM )
{
unsigned int seed = 1;
operatorInputArg(cdoOperatorEnter(operatorID));
if ( operatorArgc() < 1 ) cdoAbort("Too few arguments!");
if ( operatorArgc() > 2 ) cdoAbort("Too many arguments!");
gridfile = operatorArgv()[0];
gridID = cdoDefineGrid(gridfile);
if ( operatorArgc() == 2 )
{
long idum;
idum = atol(operatorArgv()[1]);
if ( idum >= 0 && idum < 0x7FFFFFFF )
seed = idum;
}
srand(seed);
}
else if ( operatorID == CONST )
{
operatorInputArg(cdoOperatorEnter(operatorID));
operatorCheckArgc(2);
rconst = atof(operatorArgv()[0]);
gridfile = operatorArgv()[1];
gridID = cdoDefineGrid(gridfile);
}
else if ( operatorID == TOPO || operatorID == TEMP || operatorID == MASK )
{
int nlon, nlat, i;
double lon[720], lat[360];
nlon = 720;
nlat = 360;
gridID = gridCreate(GRID_LONLAT, nlon*nlat);
gridDefXsize(gridID, nlon);
gridDefYsize(gridID, nlat);
gridfile = operatorArgv()[0];
gridID = cdoDefineGrid(gridfile);
if ( operatorArgc() >= 2 ) nlevs = atol(operatorArgv()[1]);
if ( operatorArgc() >= 3 ) ntimesteps = atol(operatorArgv()[2]);
if ( operatorArgc() >= 4 ) nvars = atol(operatorArgv()[3]);
for ( i = 0; i < nlon; i++ ) lon[i] = -179.75 + i*0.5;
for ( i = 0; i < nlat; i++ ) lat[i] = -89.75 + i*0.5;
srand(seed);
gridDefXvals(gridID, lon);
gridDefYvals(gridID, lat);
}
else if ( operatorID == FOR )
{
double lon = 0, lat = 0;
operatorInputArg(cdoOperatorEnter(operatorID));
if ( operatorArgc() < 2 ) cdoAbort("Too few arguments!");
if ( operatorArgc() > 3 ) cdoAbort("Too many arguments!");
rstart = atof(operatorArgv()[0]);
rstop = atof(operatorArgv()[1]);
if ( operatorArgc() == 3 )
rinc = atof(operatorArgv()[2]);
else
rinc = 1;
if ( DBL_IS_EQUAL(rinc, 0.0) ) cdoAbort("Increment is zero!");
gridID = gridCreate(GRID_LONLAT, 1);
gridDefXsize(gridID, 1);
gridDefYsize(gridID, 1);
gridDefXvals(gridID, &lon);
gridDefYvals(gridID, &lat);
}
else if ( operatorID == STDATM )
{
double lon = 0, lat = 0;
operatorInputArg("levels");
nlevels = args2fltlist(operatorArgc(), operatorArgv(), flist);
levels = (double *) listArrayPtr(flist);
if ( cdoVerbose ) for ( i = 0; i < nlevels; ++i ) printf("levels %d: %g\n", i, levels[i]);
gridsize = gridInqSize(gridID);
gridID = gridCreate(GRID_LONLAT, 1);
gridDefXsize(gridID, 1);
gridDefYsize(gridID, 1);
gridDefXvals(gridID, &lon);
gridDefYvals(gridID, &lat);
}
zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
if ( operatorID == STDATM )
if ( cdoVerbose )
{
zaxisID = zaxisCreate(ZAXIS_HEIGHT, nlevels);
zaxisDefLevels(zaxisID , levels);
zaxisDefName(zaxisID , "level");
zaxisDefLongname(zaxisID, "Level");
zaxisDefUnits(zaxisID , "m");
}
else
cdoPrint("gridsize : %d", gridInqSize);
cdoPrint("nlevs : %d", nlevs);
cdoPrint("ntimesteps : %d", ntimesteps);
cdoPrint("nvars : %d", nvars);
}
if ( nlevs <= 0 ) nlevs = 1;
if ( ntimesteps <= 0 ) ntimesteps = 1;
if ( nvars <= 0 ) nvars = 1;
vars = (double ***) malloc(nvars*sizeof(double **));
for ( varID = 0; varID < nvars; varID++ )
{
zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
nlevels = 1;
vars[varID] = (double **) malloc(nlevs*sizeof(double *));
for ( levelID = 0; levelID < nlevs; levelID++ )
{
vars[varID][levelID] = (double *) malloc(gridsize*sizeof(double));
for ( i = 0; i < gridsize; ++i )
vars[varID][levelID][i] = varID + rand()/(RAND_MAX+1.0);
}
}
vlistID = vlistCreate();
if ( operatorID == FOR )
varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARIABLE);
else
varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_CONSTANT);
/*
For the standard atmosphere two output variables are generated: pressure and
temperatur. The first (varID) is pressure, second (varID2) is temperatur.
Add an additional variable for the standard atmosphere.
*/
if ( operatorID == STDATM )
varID2 = vlistDefVar(vlistID, gridID, zaxisID, TIME_CONSTANT);
if ( operatorID == MASK )
vlistDefVarDatatype(vlistID, varID, DATATYPE_INT8);
if ( operatorID != STDATM )
vlistDefVarName(vlistID, varID, cdoOperatorName(operatorID));
else
for ( i = 0; i < nvars; ++i )
{
vlistDefVarName(vlistID , varID , "P");
vlistDefVarStdname(vlistID , varID , "air_pressure");
vlistDefVarLongname(vlistID, varID , "pressure");
vlistDefVarUnits(vlistID , varID , "hPa");
vlistDefVarName(vlistID , varID2, "T");
vlistDefVarStdname(vlistID , varID2, "air_temperature");
vlistDefVarLongname(vlistID, varID2, "temperature");
vlistDefVarUnits(vlistID , varID2, "K");
varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARIABLE);
vlistDefVarCode(vlistID, varID, varID+1);
// vlistDefVarName(vlistID, varID, );
}
taxisID = taxisCreate(TAXIS_RELATIVE);
vlistDefTaxis(vlistID, taxisID);
if ( operatorID == RANDOM || operatorID == CONST || operatorID == TOPO ||
operatorID == TEMP || operatorID == MASK || operatorID == STDATM )
vlistDefNtsteps(vlistID, 1);
// vlistDefNtsteps(vlistID, 1);
streamID = streamOpenWrite(cdoStreamName(0), cdoFiletype());
streamDefVlist(streamID, vlistID);
gridsize = gridInqSize(gridID);
array = (double *) malloc(gridsize*sizeof(double));
if ( operatorID == FOR )
ntimesteps = 1.001 + ((rstop-rstart)/rinc);
else
ntimesteps = 1;
julday = date_to_julday(CALENDAR_PROLEPTIC, 10101);
nvars = vlistNvars(vlistID);
nrecs = vlistNrecs(vlistID);
for ( tsID = 0; tsID < ntimesteps; tsID++ )
{
rval = rstart + rinc*tsID;
......@@ -276,58 +111,10 @@ void *CDIwrite(void *argument)
for ( varID = 0; varID < nvars; varID++ )
{
nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID, varID));
for ( levelID = 0; levelID < nlevels; levelID++ )
for ( levelID = 0; levelID < nlevs; levelID++ )
{
streamDefRecord(streamID, varID, levelID);
if ( operatorID == RANDOM )
{
for ( i = 0; i < gridsize; i++ )
array[i] = rand()/(RAND_MAX+1.0);
}
else if ( operatorID == CONST )
{
for ( i = 0; i < gridsize; i++ )
array[i] = rconst;
}
else if ( operatorID == TOPO )
{
#if defined(WITH_DATA)
for ( i = 0; i < gridsize; i++ )
array[i] = etopo[i]/etopo_scale - etopo_offset;
#else
cdoAbort("Operator support disabled!");
#endif
}
else if ( operatorID == TEMP )
{
#if defined(WITH_DATA)
for ( i = 0; i < gridsize; i++ )
array[i] = temp[i]/temp_scale - temp_offset;
#else
cdoAbort("Operator support disabled!");
#endif
}
else if ( operatorID == MASK )
{
#if defined(WITH_DATA)
for ( i = 0; i < gridsize; i++ )
array[i] = mask[i]/mask_scale - mask_offset;
#else
cdoAbort("Operator support disabled!");
#endif
}
else if ( operatorID == FOR )
{
array[0] = rval;
}
else if ( operatorID == STDATM )
{
array[0] = (varID == varID2) ? std_atm_temperatur(levels[levelID]) : std_atm_pressure(levels[levelID]);
}
streamWriteRecord(streamID, array, 0);
streamWriteRecord(streamID, vars[varID][levelID], 0);
}
}
}
......@@ -336,7 +123,12 @@ void *CDIwrite(void *argument)
vlistDestroy(vlistID);
if ( array ) free(array);
for ( varID = 0; varID < nvars; varID++ )
{
free(vars[varID]);
for ( levelID = 0; levelID < nlevs; levelID++ ) free(vars[varID][levelID]);
}
free(vars);
cdoFinish();
......
......@@ -873,6 +873,9 @@ int main(int argc, char *argv[])
#endif
#if defined (__SSE2__)
fprintf(stderr, "Predefined: __SSE2__\n");
#endif
#if defined (__GNUC__)
fprintf(stderr, "Predefined: __GNUC__\n");
#endif
fprintf(stderr, "\n");
......
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