Commit 5f6ac4a9 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Merge branch develop.

parents 20b1b815 c360b6f0
......@@ -3,6 +3,10 @@
* Using CDI library version 1.9.1
* Version 1.9.1 release
2017-09-11 Uwe Schulzweida
* ap2pl: added support for input data on half levels
2017-09-09 Uwe Schulzweida
* selindexbox: breaks uvRelativeToGrid flag [Bug #7901]
......
......@@ -212,7 +212,7 @@ void *Settime(void *argument)
}
else
{
sdate = parameter2int(datestr);
sdate = *datestr ? parameter2int(datestr) : 10101;
}
if ( operatorArgc() > 1 )
......
......@@ -48,16 +48,23 @@ bool is_height_axis(int zaxisID, int nlevel)
}
return isheight;
}
/*
static
int is_hybrid_axis(int zaxisID, int nlevel)
void change_height_zaxis(int vlistID1, int vlistID2, int zaxisID2)
{
int ishybrid = FALSE;
if ( (zaxisInqType(zaxisID) == ZAXIS_HYBRID || zaxisInqType(zaxisID) == ZAXIS_HYBRID_HALF) &&
nlevel > 1 ) ishybrid = TRUE;
return ishybrid;
int nzaxis = vlistNzaxis(vlistID1);
for ( int iz = 0; iz < nzaxis; ++iz )
{
int zaxisID = vlistZaxis(vlistID1, iz);
int nlevel = zaxisInqSize(zaxisID);
if ( is_height_axis(zaxisID, nlevel) && nlevel > 1 )
{
vlistChangeZaxisIndex(vlistID2, iz, zaxisID2);
}
}
}
*/
void *Vertintap(void *argument)
{
......@@ -66,17 +73,12 @@ void *Vertintap(void *argument)
int nrecs;
int i;
int varID, levelID;
int zaxisIDp, zaxisIDh = -1;
int zaxisIDh = -1;
int nhlev = 0, nhlevf = 0, nhlevh = 0, nlevel;
int *vert_index = NULL;
int apressID = -1, dpressID = -1;
int psID = -1, tempID = -1;
//int sortlevels = TRUE;
char paramstr[32];
char stdname[CDI_MAX_NAME];
char varname[CDI_MAX_NAME];
double *vct = NULL;
double *ps_prog = NULL, *full_press = NULL;
char varname[CDI_MAX_NAME], stdname[CDI_MAX_NAME];
bool extrapolate = false;
lista_t *flista = lista_new(FLT_LISTA);
......@@ -126,10 +128,6 @@ void *Vertintap(void *argument)
}
else
{
/*
double stdlev[] = {100000, 92500, 85000, 77500, 70000, 60000, 50000, 40000, 30000, 25000, 20000,
15000, 10000, 7000, 5000, 3000, 2000, 1000, 700, 500, 300, 200, 100, 50, 20, 10};
*/
double stdlev[] = {100000, 92500, 85000, 70000, 60000, 50000, 40000, 30000, 25000, 20000, 15000,
10000, 7000, 5000, 3000, 2000, 1000 };
nplev = sizeof(stdlev)/sizeof(*stdlev);
......@@ -154,11 +152,8 @@ void *Vertintap(void *argument)
int gridsize = vlist_check_gridsize(vlistID1);
if ( operfunc == func_hl )
zaxisIDp = zaxisCreate(ZAXIS_HEIGHT, nplev);
else
zaxisIDp = zaxisCreate(ZAXIS_PRESSURE, nplev);
int zaxistype = (operfunc == func_hl) ? ZAXIS_HEIGHT : ZAXIS_PRESSURE;
int zaxisIDp = zaxisCreate(zaxistype, nplev);
zaxisDefLevels(zaxisIDp, plev);
int nvars = vlistNvars(vlistID1);
......@@ -192,7 +187,7 @@ void *Vertintap(void *argument)
if ( zaxisID == vlistInqVarZaxis(vlistID1, apressID) )
{
bool mono_level = true;
nlevel = zaxisInqSize(zaxisID);
nlevel = zaxisInqSize(zaxisID);
if ( is_height_axis(zaxisID, nlevel) )
{
......@@ -214,17 +209,18 @@ void *Vertintap(void *argument)
nhlevf = nhlev;
nhlevh = nhlevf + 1;
vlistChangeZaxisIndex(vlistID2, i, zaxisIDp);
break;
}
}
}
change_height_zaxis(vlistID1, vlistID2, zaxisIDp);
bool *vars = (bool*) Malloc(nvars*sizeof(bool));
bool *varinterp = (bool*) Malloc(nvars*sizeof(bool));
double **vardata1 = (double**) Malloc(nvars*sizeof(double*));
double **vardata2 = (double**) Malloc(nvars*sizeof(double*));
int **varnmiss = (int**) Malloc(nvars*sizeof(int*));
std::vector<bool> vars(nvars);
std::vector<bool> varinterp(nvars);
std::vector<int *> varnmiss(nvars);
std::vector<double *> vardata1(nvars);
std::vector<double *> vardata2(nvars);
int maxlev = nhlevh > nplev ? nhlevh : nplev;
......@@ -235,8 +231,8 @@ void *Vertintap(void *argument)
{
int nlev = zaxisInqSize(zaxisIDh);
if ( nlev != nhlev ) cdoAbort("Internal error, wrong number of height level!");
double *levels = (double*) Malloc(nlev*sizeof(double));
zaxisInqLevels(zaxisIDh, levels);
std::vector<double> levels(nlev);
cdoZaxisInqLevels(zaxisIDh, levels.data());
for ( int ilev = 0; ilev < nlev; ++ilev )
{
......@@ -246,30 +242,30 @@ void *Vertintap(void *argument)
break;
}
}
Free(levels);
}
int *vert_index = NULL;
double *ps_prog = NULL, *full_press = NULL, *half_press = NULL;
if ( zaxisIDh != -1 && gridsize > 0 )
{
vert_index = (int *) Malloc(gridsize*nplev*sizeof(int));
ps_prog = (double *) Malloc(gridsize*sizeof(double));
full_press = (double *) Malloc(gridsize*nhlevf*sizeof(double));
vert_index = (int*) Malloc(gridsize*nplev*sizeof(int));
ps_prog = (double*) Malloc(gridsize*sizeof(double));
full_press = (double*) Malloc(gridsize*nhlevf*sizeof(double));
half_press = (double*) Malloc(gridsize*nhlevh*sizeof(double));
}
else
cdoWarning("No 3D variable with generalized height level found!");
if ( operfunc == func_hl )
{
double *phlev = (double*) Malloc(nplev*sizeof(double));
height2pressure(phlev, plev, nplev);
std::vector<double> phlev(nplev);
height2pressure(phlev.data(), plev, nplev);
if ( cdoVerbose )
for ( int i = 0; i < nplev; ++i )
cdoPrint("level = %d height = %g pressure = %g", i+1, plev[i], phlev[i]);
memcpy(plev, phlev, nplev*sizeof(double));
Free(phlev);
memcpy(plev, phlev.data(), nplev*sizeof(double));
}
if ( opertype == type_log )
......@@ -286,7 +282,8 @@ void *Vertintap(void *argument)
vardata1[varID] = (double *) Malloc(gridsize*nlevel*sizeof(double));
if ( zaxisID == zaxisIDh )
if ( zaxisID == zaxisIDh ||
(is_height_axis(zaxisID, nlevel) && zaxisIDh != -1 && (nlevel == nhlevh || nlevel == nhlevf)) )
{
varinterp[varID] = true;
vardata2[varID] = (double *) Malloc(gridsize*nplev*sizeof(double));
......@@ -295,20 +292,26 @@ void *Vertintap(void *argument)
}
else
{
vlistInqVarName(vlistID1, varID, varname);
if ( is_height_axis(zaxisID, nlevel) && zaxisIDh != -1 && nlevel > 1 )
{
vlistInqVarName(vlistID1, varID, varname);
cdoWarning("Parameter %d has wrong number of levels, skipped! (param=%s nlevel=%d)",
varID+1, varname, nlevel);
}
if ( is_height_axis(zaxisID, nlevel) && nlevel == nhlevh )
cdoWarning("Interpolation from half level not supported, parameter %s not interpolated!", varname);
varinterp[varID] = false;
vardata2[varID] = vardata1[varID];
varnmiss[varID] = (int *) Malloc(nlevel*sizeof(int));
}
}
if ( zaxisIDh != -1 && psID == -1 && dpressID != -1 )
cdoWarning("Surface pressure not found - set to vertical sum of %s!", var_stdname(pressure_thickness));
// cdoWarning("Surface pressure not found - set to upper level of %s!", var_stdname(air_pressure));
if ( zaxisIDh != -1 && psID == -1 )
{
if ( dpressID != -1 )
cdoWarning("Surface pressure not found - set to vertical sum of %s!", var_stdname(pressure_thickness));
else
cdoWarning("Surface pressure not found - set to lower bound of %s!", var_stdname(air_pressure));
}
for ( varID = 0; varID < nvars; ++varID )
{
......@@ -348,7 +351,7 @@ void *Vertintap(void *argument)
*/
size_t offset = gridsize*levelID;
double *single1 = vardata1[varID] + offset;
pstreamReadRecord(streamID1, single1, &varnmiss[varID][levelID]);
pstreamReadRecord(streamID1, single1, &varnmiss[varID][levelID]);
vars[varID] = true;
}
......@@ -383,10 +386,20 @@ void *Vertintap(void *argument)
memcpy(full_press, vardata1[apressID], gridsize*nhlevf*sizeof(double));
for ( int i = 0; i < gridsize; i++ ) half_press[i] = 0;
for ( int k = 1; k < nhlevf; k++ )
for ( i = 0; i < gridsize; i++ )
half_press[k*gridsize+i] = 0.5*(full_press[(k-1)*gridsize+i]+full_press[k*gridsize+i]);
for ( int i = 0; i < gridsize; i++ ) half_press[(nhlevh-1)*gridsize+i] = full_press[(nhlevf-1)*gridsize+i];
if ( opertype == type_log )
{
for ( i = 0; i < gridsize; i++ ) ps_prog[i] = log(ps_prog[i]);
for ( int k = 0; k < nhlevh; k++ )
for ( int i = 0; i < gridsize; i++ )
half_press[k*gridsize+i] = log(half_press[k*gridsize+i]);
for ( int k = 0; k < nhlevf; k++ )
for ( i = 0; i < gridsize; i++ )
full_press[k*gridsize+i] = log(full_press[k*gridsize+i]);
......@@ -407,12 +420,12 @@ void *Vertintap(void *argument)
if ( varinterp[varID] )
{
double *hyb_press = NULL;
if ( nlevel == nhlevf ) hyb_press = full_press;
if ( nlevel == nhlevf ) hyb_press = full_press;
else if ( nlevel == nhlevh ) hyb_press = half_press;
else
{
int param = vlistInqVarParam(vlistID1, varID);
cdiParamToString(param, paramstr, sizeof(paramstr));
cdoAbort("Number of generalized height level differ from full/half level (param=%s)!", paramstr);
vlistInqVarName(vlistID1, varID, varname);
cdoAbort("Number of generalized height level differ from full/half level (param=%s)!", varname);
}
for ( levelID = 0; levelID < nlevel; levelID++ )
......@@ -461,13 +474,7 @@ void *Vertintap(void *argument)
if ( ps_prog ) Free(ps_prog);
if ( vert_index ) Free(vert_index);
if ( full_press ) Free(full_press);
if ( vct ) Free(vct);
Free(vars);
Free(varinterp);
Free(vardata1);
Free(vardata2);
Free(varnmiss);
if ( half_press ) Free(half_press);
lista_destroy(flista);
......
......@@ -122,10 +122,6 @@ void *Vertintml(void *argument)
}
else
{
/*
double stdlev[] = {100000, 92500, 85000, 77500, 70000, 60000, 50000, 40000, 30000, 25000, 20000,
15000, 10000, 7000, 5000, 3000, 2000, 1000, 700, 500, 300, 200, 100, 50, 20, 10};
*/
double stdlev[] = {100000, 92500, 85000, 70000, 60000, 50000, 40000, 30000, 25000, 20000, 15000,
10000, 7000, 5000, 3000, 2000, 1000 };
nplev = sizeof(stdlev)/sizeof(*stdlev);
......@@ -150,7 +146,8 @@ void *Vertintml(void *argument)
int gridsize = vlist_check_gridsize(vlistID1);
int zaxisIDp = (operfunc == func_hl) ? zaxisCreate(ZAXIS_HEIGHT, nplev) : zaxisCreate(ZAXIS_PRESSURE, nplev);
int zaxistype = (operfunc == func_hl) ? ZAXIS_HEIGHT : ZAXIS_PRESSURE;
int zaxisIDp = zaxisCreate(zaxistype, nplev);
zaxisDefLevels(zaxisIDp, plev);
int nvct = 0;
......@@ -371,8 +368,12 @@ void *Vertintml(void *argument)
else
{
if ( zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && nlevel > 1 )
cdoWarning("Parameter %d has wrong number of levels, skipped! (param=%s nlevel=%d)",
varID+1, paramstr, nlevel);
{
vlistInqVarName(vlistID1, varID, varname);
cdoWarning("Parameter %d has wrong number of levels, skipped! (param=%s nlevel=%d)",
varID+1, varname, nlevel);
}
varinterp[varID] = false;
vardata2[varID] = vardata1[varID];
varnmiss[varID] = (int*) Malloc(nlevel*sizeof(int));
......@@ -555,13 +556,12 @@ void *Vertintml(void *argument)
}
*/
double *hyb_press = NULL;
if ( nlevel == nhlevh ) hyb_press = half_press;
else if ( nlevel == nhlevf ) hyb_press = full_press;
if ( nlevel == nhlevf ) hyb_press = full_press;
else if ( nlevel == nhlevh ) hyb_press = half_press;
else
{
int param = vlistInqVarParam(vlistID1, varID);
cdiParamToString(param, paramstr, sizeof(paramstr));
cdoAbort("Number of hybrid level differ from full/half level (param=%s)!", paramstr);
vlistInqVarName(vlistID1, varID, varname);
cdoAbort("Number of hybrid level differ from full/half level (param=%s)!", varname);
}
for ( levelID = 0; levelID < nlevel; levelID++ )
......
......@@ -4169,16 +4169,17 @@ std::vector<std::string> SpectralHelp = {
" sp2sp,trunc infile outfile",
"",
"DESCRIPTION",
" This module transforms fields on Gaussian grids to spectral coefficients and vice versa.",
" This module transforms fields on a global regular Gaussian grids to spectral coefficients and vice versa.",
" Missing values are not supported.",
"",
"OPERATORS",
" sp2gp Spectral to gridpoint",
" Convert all fields with spectral coefficients to a regular Gaussian grid. The number of ",
" Convert all fields with spectral coefficients to a global regular Gaussian grid. The number of ",
" latitudes of the resulting Gaussian grid is calculated from the triangular truncation by:",
" ",
" nlat = NINT((trunc*3 + 1.)/2.)",
" sp2gpl Spectral to gridpoint (linear)",
" Convert all fields with spectral coefficients to a regular Gaussian grid. The number of ",
" Convert all fields with spectral coefficients to a global regular Gaussian grid. The number of ",
" latitudes of the resulting Gaussian grid is calculated from the triangular truncation by:",
" ",
" nlat = NINT((trunc*2 + 1.)/2.)",
......@@ -4215,7 +4216,7 @@ std::vector<std::string> WindHelp = {
" This module converts relative divergence and vorticity to U and V wind and vice versa.",
" Divergence and vorticity are spherical harmonic coefficients in spectral space and",
" U and V are on a global regular Gaussian grid. The Gaussian latitudes need to be ordered from",
" north to south.",
" north to south. Missing values are not supported.",
"",
"OPERATORS",
" dv2uv Divergence and vorticity to U and V wind",
......
Markdown is supported
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