Commit 10430ae1 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Use VarList.

parent d058ea89
Pipeline #4910 passed with stages
in 16 minutes and 34 seconds
......@@ -108,14 +108,8 @@ calc_adisit(size_t gridsize, size_t nlevel, const Varray<double> &pressure, cons
const auto tis_missval = tis[levelID].missval;
for (size_t i = 0; i < gridsize; ++i)
{
if (DBL_IS_EQUAL(thovec[i], tho_missval) || DBL_IS_EQUAL(saovec[i], sao_missval))
{
tisvec[i] = tis_missval;
}
else
{
tisvec[i] = adisit_1(thovec[i], saovec[i], pressure[levelID]);
}
const auto is_missing = (DBL_IS_EQUAL(thovec[i], tho_missval) || DBL_IS_EQUAL(saovec[i], sao_missval));
tisvec[i] = is_missing ? tis_missval : adisit_1(thovec[i], saovec[i], pressure[levelID]);
}
}
}
......@@ -138,14 +132,8 @@ calc_adipot(size_t gridsize, size_t nlevel, const Varray<double> &pressure, cons
const auto tpot_missval = tpot[levelID].missval;
for (size_t i = 0; i < gridsize; ++i)
{
if (DBL_IS_EQUAL(tvec[i], t_missval) || DBL_IS_EQUAL(svec[i], s_missval))
{
tpotvec[i] = tpot_missval;
}
else
{
tpotvec[i] = adipot(tvec[i], svec[i], pressure[levelID]);
}
const auto is_missing = (DBL_IS_EQUAL(tvec[i], t_missval) || DBL_IS_EQUAL(svec[i], s_missval));
tpotvec[i] = is_missing ? tpot_missval : adipot(tvec[i], svec[i], pressure[levelID]);
}
}
}
......@@ -166,7 +154,7 @@ Adisit(void *process)
const auto operatorID = cdoOperatorID();
const double pin = (operatorArgc() == 1) ? parameter2double(cdoOperatorArgv(0)) : -1;
const auto pin = (operatorArgc() == 1) ? parameter2double(cdoOperatorArgv(0)) : -1.0;
const auto streamID1 = cdoOpenRead(0);
const auto vlistID1 = cdoStreamInqVlist(streamID1);
......@@ -175,7 +163,7 @@ Adisit(void *process)
for (varID = 0; varID < nvars; varID++)
{
int code = vlistInqVarCode(vlistID1, varID);
auto code = vlistInqVarCode(vlistID1, varID);
if (code <= 0)
{
vlistInqVarName(vlistID1, varID, varname);
......@@ -190,14 +178,8 @@ Adisit(void *process)
else if (cdo_cmpstr(varname, "t"))
code = 2;
if (operatorID == ADISIT)
{
if (cdo_cmpstr(stdname, "sea_water_potential_temperature")) code = 2;
}
else
{
if (cdo_cmpstr(stdname, "sea_water_temperature")) code = 2;
}
const char *cname = (operatorID == ADISIT) ? "sea_water_potential_temperature" : "sea_water_temperature";
if (cdo_cmpstr(stdname, cname)) code = 2;
}
if (code == 2 || (code == 20 && operatorID == ADIPOT))
......@@ -213,29 +195,29 @@ Adisit(void *process)
const auto gridsize = vlist_check_gridsize(vlistID1);
auto zaxisID = vlistInqVarZaxis(vlistID1, saoID);
const auto nlevel1 = zaxisInqSize(zaxisID);
const auto nlevels1 = zaxisInqSize(zaxisID);
zaxisID = vlistInqVarZaxis(vlistID1, thoID);
const auto nlevel2 = zaxisInqSize(zaxisID);
const auto nlevels2 = zaxisInqSize(zaxisID);
if (nlevel1 != nlevel2) cdoAbort("temperature and salinity have different number of levels!");
const auto nlevel = nlevel1;
if (nlevels1 != nlevels2) cdoAbort("temperature and salinity have different number of levels!");
const auto nlevels = nlevels1;
Varray<double> pressure(nlevel);
Varray<double> pressure(nlevels);
cdoZaxisInqLevels(zaxisID, pressure.data());
if (pin >= 0)
for (i = 0; i < nlevel; ++i) pressure[i] = pin;
for (i = 0; i < nlevels; ++i) pressure[i] = pin;
else
for (i = 0; i < nlevel; ++i) pressure[i] /= 10;
for (i = 0; i < nlevels; ++i) pressure[i] /= 10;
if (Options::cdoVerbose)
{
cdoPrint("Level Pressure");
for (i = 0; i < nlevel; ++i) cdoPrint("%5d %g", i + 1, pressure[i]);
for (i = 0; i < nlevels; ++i) cdoPrint("%5d %g", i + 1, pressure[i]);
}
FieldVector tho(nlevel), sao(nlevel), tis(nlevel);
for (levelID = 0; levelID < nlevel; ++levelID)
FieldVector tho(nlevels), sao(nlevels), tis(nlevels);
for (levelID = 0; levelID < nlevels; ++levelID)
{
tho[levelID].resize(gridsize);
sao[levelID].resize(gridsize);
......@@ -300,11 +282,11 @@ Adisit(void *process)
}
if (operatorID == ADISIT)
calc_adisit(gridsize, nlevel, pressure, tho, sao, tis);
calc_adisit(gridsize, nlevels, pressure, tho, sao, tis);
else
calc_adipot(gridsize, nlevel, pressure, tho, sao, tis);
calc_adipot(gridsize, nlevels, pressure, tho, sao, tis);
for (levelID = 0; levelID < nlevel; ++levelID)
for (levelID = 0; levelID < nlevels; ++levelID)
{
cdoDefRecord(streamID2, tisID2, levelID);
cdoWriteRecord(streamID2, tis[levelID].vec_d.data(), fieldNumMiss(tis[levelID]));
......
......@@ -29,7 +29,6 @@ void *
Fldrms(void *process)
{
int lastgrid = -1;
int index;
int nrecs;
int varID, levelID;
......@@ -50,7 +49,7 @@ Fldrms(void *process)
const auto taxisID3 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID3, taxisID3);
double slon = 0, slat = 0;
double slon = 0.0, slat = 0.0;
const auto gridID3 = gridCreate(GRID_LONLAT, 1);
gridDefXsize(gridID3, 1);
gridDefYsize(gridID3, 1);
......@@ -59,19 +58,18 @@ Fldrms(void *process)
const auto ngrids = vlistNgrids(vlistID1);
int ndiffgrids = 0;
for (index = 1; index < ngrids; index++)
for (int index = 1; index < ngrids; index++)
if (vlistGrid(vlistID1, 0) != vlistGrid(vlistID1, index)) ndiffgrids++;
index = 0;
const auto gridID1 = vlistGrid(vlistID1, index);
const auto gridID2 = vlistGrid(vlistID2, index);
const auto gridID1 = vlistGrid(vlistID1, 0);
const auto gridID2 = vlistGrid(vlistID2, 0);
if (gridInqSize(gridID1) != gridInqSize(gridID2)) cdoAbort("Fields have different grid size!");
if (needWeights && gridInqType(gridID1) != GRID_LONLAT && gridInqType(gridID1) != GRID_GAUSSIAN)
cdoAbort("Unsupported gridtype: %s", gridNamePtr(gridInqType(gridID1)));
for (index = 0; index < ngrids; index++) vlistChangeGridIndex(vlistID3, index, gridID3);
for (int index = 0; index < ngrids; index++) vlistChangeGridIndex(vlistID3, index, gridID3);
if (ndiffgrids > 0) cdoAbort("Too many different grids!");
......@@ -80,6 +78,10 @@ Fldrms(void *process)
const auto gridsizemax = vlistGridsizeMax(vlistID1);
VarList varList1, varList2;
varListInit(varList1, vlistID1);
varListInit(varList2, vlistID2);
Field field1, field2, field3;
field1.resize(gridsizemax);
if (needWeights) field1.weightv.resize(gridsizemax);
......@@ -105,8 +107,8 @@ Fldrms(void *process)
cdoInqRecord(streamID2, &varID, &levelID);
cdoReadRecord(streamID2, field2.vec_d.data(), &field2.nmiss);
field1.grid = vlistInqVarGrid(vlistID1, varID);
field2.grid = vlistInqVarGrid(vlistID2, varID);
field1.grid = varList1[varID].gridID;
field2.grid = varList2[varID].gridID;
if (needWeights && field1.grid != lastgrid)
{
......@@ -116,17 +118,13 @@ Fldrms(void *process)
{
const auto wstatus = gridWeights(field1.grid, field1.weightv.data());
if (wstatus != 0 && tsID == 0 && levelID == 0)
{
char varname[CDI_MAX_NAME];
vlistInqVarName(vlistID1, varID, varname);
cdoWarning("Grid cell bounds not available, using constant grid cell area weights for variable %s!", varname);
}
cdoWarning("Grid cell bounds not available, using constant grid cell area weights for variable %s!", varList1[varID].name);
}
}
field1.missval = vlistInqVarMissval(vlistID1, varID);
field2.missval = vlistInqVarMissval(vlistID1, varID);
field3.missval = vlistInqVarMissval(vlistID1, varID);
field1.missval = varList1[varID].missval;
field2.missval = varList1[varID].missval;
field3.missval = varList1[varID].missval;
vfldrms(field1, field2, field3);
......
......@@ -28,13 +28,13 @@
#include "cdo_vlist.h"
#include <mpim_grid.h>
/* routine corr copied from PINGO */
/* correclation in space */
// routine corr copied from PINGO
// correclation in space
static double
correlation_s(const Varray<double> &in0, const Varray<double> &in1, const Varray<double> &weight, double missval1,
double missval2, size_t gridsize)
{
double sum0 = 0, sum1 = 0, sum00 = 0, sum01 = 0, sum11 = 0, wsum0 = 0;
double sum0 = 0.0, sum1 = 0.0, sum00 = 0.0, sum01 = 0.0, sum11 = 0.0, wsum0 = 0.0;
for (size_t i = 0; i < gridsize; ++i)
{
......@@ -49,19 +49,19 @@ correlation_s(const Varray<double> &in0, const Varray<double> &in1, const Varray
}
}
double out = IS_NOT_EQUAL(wsum0, 0)
double out = IS_NOT_EQUAL(wsum0, 0.0)
? DIVMN((sum01 * wsum0 - sum0 * sum1), SQRTMN((sum00 * wsum0 - sum0 * sum0) * (sum11 * wsum0 - sum1 * sum1)))
: missval1;
return out;
}
/* covariance in space */
// covariance in space
static double
covariance_s(const Varray<double> &in0, const Varray<double> &in1, const Varray<double> &weight, double missval1,
double missval2, size_t gridsize)
{
double sum0 = 0, sum1 = 0, sum01 = 0, wsum0 = 0, wsum00 = 0;
double sum0 = 0.0, sum1 = 0.0, sum01 = 0.0, wsum0 = 0.0, wsum00 = 0.0;
for (size_t i = 0; i < gridsize; ++i)
{
......@@ -75,7 +75,7 @@ covariance_s(const Varray<double> &in0, const Varray<double> &in1, const Varray<
}
}
double out = IS_NOT_EQUAL(wsum0, 0) ? (sum01 * wsum0 - sum0 * sum1) / (wsum0 * wsum0) : missval1;
double out = IS_NOT_EQUAL(wsum0, 0.0) ? (sum01 * wsum0 - sum0 * sum1) / (wsum0 * wsum0) : missval1;
return out;
}
......@@ -83,14 +83,12 @@ covariance_s(const Varray<double> &in0, const Varray<double> &in1, const Varray<
void *
Fldstat2(void *process)
{
int gridID, lastgridID = -1;
int lastgridID = -1;
int nrecs;
int varID, levelID;
size_t nmiss1, nmiss2;
bool wstatus = false;
bool needWeights = true;
double sglval = 0;
char varname[CDI_MAX_NAME];
cdoInitialize(process);
......@@ -115,7 +113,11 @@ Fldstat2(void *process)
const auto taxisID3 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID3, taxisID3);
double slon = 0, slat = 0;
VarList varList1, varList2;
varListInit(varList1, vlistID1);
varListInit(varList2, vlistID2);
double slon = 0.0, slat = 0.0;
const auto gridID3 = gridCreate(GRID_LONLAT, 1);
gridDefXsize(gridID3, 1);
gridDefYsize(gridID3, 1);
......@@ -156,22 +158,20 @@ Fldstat2(void *process)
cdoReadRecord(streamID1, array1.data(), &nmiss1);
cdoReadRecord(streamID2, array2.data(), &nmiss2);
gridID = vlistInqVarGrid(vlistID1, varID);
const auto gridsize = gridInqSize(gridID);
const auto gridID = varList1[varID].gridID;
const auto gridsize = varList1[varID].gridsize;
if (needWeights && gridID != lastgridID)
{
lastgridID = gridID;
wstatus = gridWeights(gridID, weight.data()) != 0;
}
if (wstatus && tsID == 0 && levelID == 0)
{
vlistInqVarName(vlistID1, varID, varname);
cdoWarning("Using constant grid cell area weights for variable %s!", varname);
}
cdoWarning("Using constant grid cell area weights for variable %s!", varList1[varID].name);
const auto missval1 = vlistInqVarMissval(vlistID1, varID);
const auto missval2 = vlistInqVarMissval(vlistID2, varID);
const auto missval1 = varList1[varID].missval;
const auto missval2 = varList2[varID].missval;
double sglval = 0.0;
if (operfunc == func_cor)
sglval = correlation_s(array1, array2, weight, missval1, missval2, gridsize);
else if (operfunc == func_covar)
......
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