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

Replaced malloc() by std::vector.

parent 181e31a9
......@@ -22,6 +22,7 @@
*/
#include <array>
#include <cdi.h>
#include "cdo.h"
#include "cdo_int.h"
......@@ -149,7 +150,7 @@ void *Vertintap(void *argument)
int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int gridsize = vlist_check_gridsize(vlistID1);
size_t gridsize = vlist_check_gridsize(vlistID1);
int zaxistype = (operfunc == func_hl) ? ZAXIS_HEIGHT : ZAXIS_PRESSURE;
int zaxisIDp = zaxisCreate(zaxistype, nplev);
......@@ -170,11 +171,15 @@ void *Vertintap(void *argument)
if ( cdoVerbose )
{
std::vector<std::array<char, CDI_MAX_NAME>> varNames(nvars);
for ( varID = 0; varID < nvars; varID++ )
vlistInqVarName(vlistID1, varID, &varNames[varID][0]);
cdoPrint("Found:");
if ( psID != -1 ) cdoPrint(" %s", var_stdname(surface_air_pressure));
if ( apressID != -1 ) cdoPrint(" %s", var_stdname(air_pressure));
if ( dpressID != -1 ) cdoPrint(" %s", var_stdname(pressure_thickness));
if ( tempID != -1 ) cdoPrint(" %s", var_stdname(air_temperature));
if ( psID != -1 ) cdoPrint(" %s -> %s", var_stdname(surface_air_pressure), &varNames[psID][0]);
if ( apressID != -1 ) cdoPrint(" %s -> %s", var_stdname(air_pressure), &varNames[apressID][0]);
if ( dpressID != -1 ) cdoPrint(" %s -> %s", var_stdname(pressure_thickness), &varNames[dpressID][0]);
if ( tempID != -1 ) cdoPrint(" %s -> %s", var_stdname(air_temperature), &varNames[tempID][0]);
}
if ( apressID == -1 ) cdoAbort("%s not found!", var_stdname(air_pressure));
......@@ -190,15 +195,14 @@ void *Vertintap(void *argument)
if ( is_height_axis(zaxisID, nlevel) )
{
double *level = (double *) Malloc(nlevel*sizeof(double));
cdoZaxisInqLevels(zaxisID, level);
std::vector<double> level(nlevel);
cdoZaxisInqLevels(zaxisID, &level[0]);
int l;
for ( l = 0; l < nlevel; l++ )
{
if ( (l+1) != (int) (level[l]+0.5) ) break;
}
if ( l == nlevel ) mono_level = true;
Free(level);
}
if ( is_height_axis(zaxisID, nlevel) && mono_level )
......@@ -223,7 +227,8 @@ void *Vertintap(void *argument)
int maxlev = nhlevh > nplev ? nhlevh : nplev;
size_t *pnmiss = extrapolate ? NULL : (size_t *) Malloc(nplev*sizeof(size_t));
std::vector<size_t> pnmiss;
if ( !extrapolate ) pnmiss.resize(nplev);
// check levels
if ( zaxisIDh != -1 )
......@@ -243,14 +248,14 @@ void *Vertintap(void *argument)
}
}
int *vert_index = NULL;
double *ps_prog = NULL, *full_press = NULL, *half_press = NULL;
std::vector<int> vert_index;
std::vector<double> ps_prog, full_press, half_press;
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));
half_press = (double*) Malloc(gridsize*nhlevh*sizeof(double));
vert_index.resize(gridsize*nplev);
ps_prog.resize(gridsize);
full_press.resize(gridsize*nhlevf);
half_press.resize(gridsize*nhlevh);
}
else
cdoWarning("No 3D variable with generalized height level found!");
......@@ -362,51 +367,45 @@ void *Vertintap(void *argument)
{
if ( psID != -1 )
{
memcpy(ps_prog, vardata1[psID], gridsize*sizeof(double));
memcpy(&ps_prog[0], vardata1[psID], gridsize*sizeof(double));
}
else if ( dpressID != -1 )
{
for ( int i = 0; i < gridsize; i++ ) ps_prog[i] = 0;
for ( size_t i = 0; i < gridsize; i++ ) ps_prog[i] = 0;
for ( int k = 0; k < nhlevf; ++k )
for ( int i = 0; i < gridsize; i++ )
for ( size_t i = 0; i < gridsize; i++ )
ps_prog[i] += vardata1[dpressID][k*gridsize+i];
}
else
{
memcpy(ps_prog, vardata1[apressID]+gridsize*(nhlevf-1), gridsize*sizeof(double));
//for ( int i = 0; i < gridsize; i++ ) ps_prog[i] = 110000;
memcpy(&ps_prog[0], vardata1[apressID]+gridsize*(nhlevf-1), gridsize*sizeof(double));
//for ( size_t i = 0; i < gridsize; i++ ) ps_prog[i] = 110000;
}
/* check range of ps_prog */
double minval, maxval;
minmaxval(gridsize, ps_prog, NULL, &minval, &maxval);
minmaxval(gridsize, &ps_prog[0], NULL, &minval, &maxval);
if ( minval < MIN_PS || maxval > MAX_PS )
cdoWarning("Surface pressure out of range (min=%g max=%g)!", minval, maxval);
memcpy(full_press, vardata1[apressID], gridsize*nhlevf*sizeof(double));
memcpy(&full_press[0], vardata1[apressID], gridsize*nhlevf*sizeof(double));
for ( int i = 0; i < gridsize; i++ ) half_press[i] = 0;
for ( size_t i = 0; i < gridsize; i++ ) half_press[i] = 0;
for ( int k = 1; k < nhlevf; k++ )
for ( int i = 0; i < gridsize; i++ )
for ( size_t 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];
for ( size_t i = 0; i < gridsize; i++ ) half_press[(nhlevh-1)*gridsize+i] = full_press[(nhlevf-1)*gridsize+i];
if ( opertype == type_log )
{
for ( int 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 ( int i = 0; i < gridsize; i++ )
full_press[k*gridsize+i] = log(full_press[k*gridsize+i]);
for ( size_t i = 0; i < gridsize; i++ ) ps_prog[i] = log(ps_prog[i]);
for ( size_t ki = 0; ki < nhlevh*gridsize; ki++ ) half_press[ki] = log(half_press[ki]);
for ( size_t ki = 0; ki < nhlevf*gridsize; ki++ ) full_press[ki] = log(full_press[ki]);
}
genind(vert_index, plev, full_press, gridsize, nplev, nhlevf);
genind(&vert_index[0], plev, &full_press[0], gridsize, nplev, nhlevf);
if ( !extrapolate ) genindmiss(vert_index, plev, gridsize, nplev, ps_prog, pnmiss);
if ( !extrapolate ) genindmiss(&vert_index[0], plev, gridsize, nplev, &ps_prog[0], &pnmiss[0]);
}
for ( varID = 0; varID < nvars; varID++ )
......@@ -419,8 +418,8 @@ void *Vertintap(void *argument)
if ( varinterp[varID] )
{
double *hyb_press = NULL;
if ( nlevel == nhlevf ) hyb_press = full_press;
else if ( nlevel == nhlevh ) hyb_press = half_press;
if ( nlevel == nhlevf ) hyb_press = &full_press[0];
else if ( nlevel == nhlevh ) hyb_press = &half_press[0];
else
{
vlistInqVarName(vlistID1, varID, varname);
......@@ -434,9 +433,9 @@ void *Vertintap(void *argument)
}
interp_X(vardata1[varID], vardata2[varID], hyb_press,
vert_index, plev, nplev, gridsize, nlevel, missval);
&vert_index[0], plev, nplev, gridsize, nlevel, missval);
if ( !extrapolate ) memcpy(varnmiss[varID], pnmiss, nplev*sizeof(size_t));
if ( !extrapolate ) memcpy(varnmiss[varID], &pnmiss[0], nplev*sizeof(size_t));
}
}
}
......@@ -469,12 +468,6 @@ void *Vertintap(void *argument)
if ( varinterp[varID] ) Free(vardata2[varID]);
}
if ( pnmiss ) Free(pnmiss);
if ( ps_prog ) Free(ps_prog);
if ( vert_index ) Free(vert_index);
if ( full_press ) Free(full_press);
if ( half_press ) Free(half_press);
lista_destroy(flista);
cdoFinish();
......
......@@ -73,7 +73,6 @@ void *Vertintml(void *argument)
char paramstr[32];
char varname[CDI_MAX_NAME], stdname[CDI_MAX_NAME];
double minval, maxval;
double *sgeopot = NULL;
gribcode_t gribcodes;
memset(&gribcodes, 0, sizeof(gribcode_t));
lista_t *flista = lista_new(FLT_LISTA);
......@@ -145,7 +144,7 @@ void *Vertintml(void *argument)
int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int gridsize = vlist_check_gridsize(vlistID1);
size_t gridsize = vlist_check_gridsize(vlistID1);
int zaxistype = (operfunc == func_hl) ? ZAXIS_HEIGHT : ZAXIS_PRESSURE;
int zaxisIDp = zaxisCreate(zaxistype, nplev);
......@@ -193,7 +192,8 @@ void *Vertintml(void *argument)
int maxlev = nhlevh > nplev ? nhlevh : nplev;
size_t *pnmiss = extrapolate ? NULL : (size_t*) Malloc(nplev*sizeof(size_t));
std::vector<size_t> pnmiss;
if ( !extrapolate ) pnmiss.resize(nplev);
// check levels
if ( zaxisIDh != -1 )
......@@ -213,14 +213,14 @@ void *Vertintml(void *argument)
}
}
int *vert_index = NULL;
double *ps_prog = NULL, *full_press = NULL, *half_press = NULL;
std::vector<int> vert_index;
std::vector<double> ps_prog, full_press, half_press;
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));
half_press = (double*) Malloc(gridsize*nhlevh*sizeof(double));
vert_index.resize(gridsize*nplev);
ps_prog.resize(gridsize);
full_press.resize(gridsize*nhlevf);
half_press.resize(gridsize*nhlevh);
}
else
cdoWarning("No 3D variable with hybrid sigma pressure coordinate found!");
......@@ -381,20 +381,25 @@ void *Vertintml(void *argument)
if ( cdoVerbose )
{
std::vector<std::array<char, CDI_MAX_NAME>> varNames(nvars);
for ( varID = 0; varID < nvars; varID++ )
vlistInqVarName(vlistID1, varID, &varNames[varID][0]);
cdoPrint("Found:");
if ( tempID != -1 ) cdoPrint(" %s", var_stdname(air_temperature));
if ( psID != -1 ) cdoPrint(" %s", var_stdname(surface_air_pressure));
if ( lnpsID != -1 ) cdoPrint(" LOG(%s)", var_stdname(surface_air_pressure));
if ( sgeopotID != -1 ) cdoPrint(" %s", var_stdname(surface_geopotential));
if ( geopotID != -1 ) cdoPrint(" %s", var_stdname(geopotential));
if ( gheightID != -1 ) cdoPrint(" %s", var_stdname(geopotential_height));
if ( tempID != -1 ) cdoPrint(" %s -> %s", var_stdname(air_temperature), &varNames[tempID][0]);
if ( psID != -1 ) cdoPrint(" %s -> %s", var_stdname(surface_air_pressure), &varNames[psID][0]);
if ( lnpsID != -1 ) cdoPrint(" LOG(%s) -> %s", var_stdname(surface_air_pressure), &varNames[lnpsID][0]);
if ( sgeopotID != -1 ) cdoPrint(" %s -> %s", var_stdname(surface_geopotential), &varNames[sgeopotID][0]);
if ( geopotID != -1 ) cdoPrint(" %s -> %s", var_stdname(geopotential), &varNames[geopotID][0]);
if ( gheightID != -1 ) cdoPrint(" %s -> %s", var_stdname(geopotential_height), &varNames[gheightID][0]);
}
if ( tempID != -1 || gheightID != -1 ) sgeopot_needed = true;
std::vector<double> sgeopot;
if ( zaxisIDh != -1 && sgeopot_needed )
{
sgeopot = (double*) Malloc(gridsize*sizeof(double));
sgeopot.resize(gridsize);
if ( sgeopotID == -1 )
{
if ( extrapolate )
......@@ -404,7 +409,7 @@ void *Vertintml(void *argument)
else
cdoPrint("%s not found - using bottom layer of %s!", var_stdname(surface_geopotential), var_stdname(geopotential));
}
memset(sgeopot, 0, gridsize*sizeof(double));
memset(&sgeopot[0], 0, gridsize*sizeof(double));
}
}
......@@ -482,14 +487,14 @@ void *Vertintml(void *argument)
if ( sgeopot_needed )
{
if ( sgeopotID != -1 )
memcpy(sgeopot, vardata1[sgeopotID], gridsize*sizeof(double));
memcpy(&sgeopot[0], vardata1[sgeopotID], gridsize*sizeof(double));
else if ( geopotID != -1 )
memcpy(sgeopot, vardata1[geopotID]+gridsize*(nhlevf-1), gridsize*sizeof(double));
memcpy(&sgeopot[0], vardata1[geopotID]+gridsize*(nhlevf-1), gridsize*sizeof(double));
/* check range of surface geopot */
if ( extrapolate && (sgeopotID != -1 || geopotID != -1) )
{
minmaxval(gridsize, sgeopot, NULL, &minval, &maxval);
minmaxval(gridsize, &sgeopot[0], NULL, &minval, &maxval);
if ( minval < MIN_FIS || maxval > MAX_FIS )
cdoWarning("Surface geopotential out of range (min=%g max=%g) [timestep:%d]!", minval, maxval, tsID+1);
if ( gridsize > 1 && minval >= 0 && maxval <= 9000 )
......@@ -498,34 +503,28 @@ void *Vertintml(void *argument)
}
if ( presID == lnpsID )
for ( int i = 0; i < gridsize; i++ ) ps_prog[i] = exp(vardata1[lnpsID][i]);
for ( size_t i = 0; i < gridsize; i++ ) ps_prog[i] = exp(vardata1[lnpsID][i]);
else if ( presID != -1 )
memcpy(ps_prog, vardata1[presID], gridsize*sizeof(double));
memcpy(&ps_prog[0], vardata1[presID], gridsize*sizeof(double));
/* check range of ps_prog */
minmaxval(gridsize, ps_prog, NULL, &minval, &maxval);
minmaxval(gridsize, &ps_prog[0], NULL, &minval, &maxval);
if ( minval < MIN_PS || maxval > MAX_PS )
cdoWarning("Surface pressure out of range (min=%g max=%g) [timestep:%d]!", minval, maxval, tsID+1);
presh(full_press, half_press, vct, ps_prog, nhlevf, gridsize);
presh(&full_press[0], &half_press[0], vct, &ps_prog[0], nhlevf, gridsize);
if ( opertype == type_log )
{
for ( int 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 ( int i = 0; i < gridsize; i++ )
full_press[k*gridsize+i] = log(full_press[k*gridsize+i]);
for ( size_t i = 0; i < gridsize; i++ ) ps_prog[i] = log(ps_prog[i]);
for ( size_t ki = 0; ki < nhlevh*gridsize; ki++ ) half_press[ki] = log(half_press[ki]);
for ( size_t ki = 0; ki < nhlevf*gridsize; ki++ ) full_press[ki] = log(full_press[ki]);
}
genind(vert_index, plev, full_press, gridsize, nplev, nhlevf);
genind(&vert_index[0], plev, &full_press[0], gridsize, nplev, nhlevf);
if ( !extrapolate ) genindmiss(vert_index, plev, gridsize, nplev, ps_prog, pnmiss);
if ( !extrapolate ) genindmiss(&vert_index[0], plev, gridsize, nplev, &ps_prog[0], &pnmiss[0]);
}
for ( varID = 0; varID < nvars; varID++ )
......@@ -555,8 +554,8 @@ void *Vertintml(void *argument)
}
*/
double *hyb_press = NULL;
if ( nlevel == nhlevf ) hyb_press = full_press;
else if ( nlevel == nhlevh ) hyb_press = half_press;
if ( nlevel == nhlevf ) hyb_press = &full_press[0];
else if ( nlevel == nhlevh ) hyb_press = &half_press[0];
else
{
vlistInqVarName(vlistID1, varID, varname);
......@@ -577,26 +576,26 @@ void *Vertintml(void *argument)
if ( opertype == type_log && extrapolate )
cdoAbort("Log. extrapolation of temperature unsupported!");
interp_T(sgeopot, vardata1[varID], vardata2[varID],
full_press, half_press, vert_index,
interp_T(&sgeopot[0], vardata1[varID], vardata2[varID],
&full_press[0], &half_press[0], &vert_index[0],
plev, nplev, gridsize, nlevel, missval);
}
else if ( varID == gheightID )
{
for ( int i = 0; i < gridsize; ++i )
for ( size_t i = 0; i < gridsize; ++i )
vardata1[varID][gridsize*nlevel+i] = sgeopot[i]/PlanetGrav;
interp_Z(sgeopot, vardata1[varID], vardata2[varID],
full_press, half_press, vert_index, vardata1[tempID],
interp_Z(&sgeopot[0], vardata1[varID], vardata2[varID],
&full_press[0], &half_press[0], &vert_index[0], vardata1[tempID],
plev, nplev, gridsize, nlevel, missval);
}
else
{
interp_X(vardata1[varID], vardata2[varID], hyb_press,
vert_index, plev, nplev, gridsize, nlevel, missval);
&vert_index[0], plev, nplev, gridsize, nlevel, missval);
}
if ( !extrapolate ) memcpy(varnmiss[varID], pnmiss, nplev*sizeof(size_t));
if ( !extrapolate ) memcpy(varnmiss[varID], &pnmiss[0], nplev*sizeof(size_t));
}
}
}
......@@ -630,13 +629,7 @@ void *Vertintml(void *argument)
if ( varinterp[varID] ) Free(vardata2[varID]);
}
if ( pnmiss ) Free(pnmiss);
if ( sgeopot ) Free(sgeopot);
if ( ps_prog ) Free(ps_prog);
if ( vert_index ) Free(vert_index);
if ( full_press ) Free(full_press);
if ( half_press ) Free(half_press);
if ( vct ) Free(vct);
if ( vct ) Free(vct);
lista_destroy(flista);
......
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