Skip to content

GitLab

Commit 71e226d2 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Vertintml: Added memory support for 32-bit float data.

parent 30c526b4
Hide whitespace changes
Inline Side-by-side
ChangeLog
View file @ 71e226d2
2020-03-24 Uwe Schulzweida
* Vertintml: Added memory support for 32-bit float data.
2020-03-21 Uwe Schulzweida
* Diff: Added memory support for 32-bit float data.
......
src/Derivepar.cc
View file @ 71e226d2
......@@ -370,7 +370,7 @@ Derivepar(void *process)
if (operatorID == GHEIGHT)
{
presh(nullptr, half_press.data(), vct.data(), ps.data(), nhlevf, gridsize);
presh((double*)nullptr, half_press.data(), vct.data(), ps.data(), nhlevf, gridsize);
arrayCopy(gridsize, sgeopot.data(), gheight.data() + gridsize * nhlevf);
MakeGeopotHeight(gheight.data(), temp.data(), hum.data(), half_press.data(), gridsize, nhlevf);
......
src/Makefile.am
View file @ 71e226d2
......@@ -91,6 +91,8 @@ libcdo_la_SOURCES = after_dvtrans.cc \
fieldmem.cc \
fieldmer.cc \
fieldzon.cc \
field_vinterp.cc \
field_vinterp.h \
fileStream.cc \
fileStream.h \
functs.h \
......
src/Remapeta.cc
View file @ 71e226d2
......@@ -704,7 +704,7 @@ Remapeta(void *process)
}
else if (operatorID == REMAPETAZ)
{
presh(nullptr, half_press1.data(), vct1.data(), ps1.data(), nhlevf1, gridsize);
presh((double*)nullptr, half_press1.data(), vct1.data(), ps1.data(), nhlevf1, gridsize);
for (int k = 0; k < nhlevf1; ++k)
for (size_t i = 0; i < gridsize; ++i)
{
......@@ -713,7 +713,7 @@ Remapeta(void *process)
}
vertSumw(sum1, vars1[iv], gridsize, nhlevf1, deltap1);
presh(nullptr, half_press2.data(), vct2.data(), ps1.data(), nhlevf2, gridsize);
presh((double*)nullptr, half_press2.data(), vct2.data(), ps1.data(), nhlevf2, gridsize);
for (int k = 0; k < nhlevf2; ++k)
for (size_t i = 0; i < gridsize; ++i)
{
......
src/Verifygrid.cc
View file @ 71e226d2
......@@ -769,8 +769,16 @@ Verifygrid(void *argument)
cdiInqKeyInt(gridID, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, &number);
if (number > 0)
{
gridID = referenceToGrid(gridID);
if (gridID == -1) cdoAbort("Reference to source grid not found!");
auto newgridID = referenceToGrid(gridID);
if (newgridID == -1)
{
lgeo = false;
cdoPrint("Reference to source grid not found!");
}
else
{
gridID = newgridID;
}
}
}
}
......
src/Vertintap.cc
View file @ 71e226d2
......@@ -26,7 +26,7 @@
#include "cdo_options.h"
#include "process_int.h"
#include "cdo_vlist.h"
#include "vertical_interp.h"
#include "field_vinterp.h"
#include "stdnametable.h"
#include "util_string.h"
#include "const.h"
......@@ -90,46 +90,6 @@ calc_half_press(const Field3D &full_press, Field3D &half_press)
calc_half_press(full_press.gridsize, full_press.nlevels, full_press.vec_d, half_press.nlevels, half_press.vec_d);
}
static void
genind(std::vector<int> &vert_index, Varray<double> &plev, Field3D &full_press, size_t gridsize)
{
const auto nplev = plev.size();
const auto nhlevf = full_press.nlevels;
if (full_press.memType == MemType::Float)
genind(vert_index.data(), plev.data(), full_press.vec_f.data(), gridsize, nplev, nhlevf);
else
genind(vert_index.data(), plev.data(), full_press.vec_d.data(), gridsize, nplev, nhlevf);
}
static void
genindmiss(std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize, Field &ps_prog, std::vector<size_t> &pnmiss)
{
const auto nplev = plev.size();
if (ps_prog.memType == MemType::Float)
genindmiss(vert_index.data(), plev.data(), gridsize, nplev, ps_prog.vec_f.data(), pnmiss.data());
else
genindmiss(vert_index.data(), plev.data(), gridsize, nplev, ps_prog.vec_d.data(), pnmiss.data());
}
static void
vertical_interp_X(size_t nlevels, Field3D &full_press, Field3D &half_press, Field3D &field1, Field3D &field2, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize)
{
const auto nplev = plev.size();
const auto missval = field1.missval;
if (field1.memType == MemType::Float)
{
const auto &hyb_press = (nlevels == full_press.nlevels) ? full_press.vec_f : half_press.vec_f;
vertical_interp_X(field1.vec_f.data(), field2.vec_f.data(), hyb_press.data(), &vert_index[0], plev.data(), nplev,
gridsize, nlevels, missval);
}
else
{
const auto &hyb_press = (nlevels == full_press.nlevels) ? full_press.vec_d : half_press.vec_d;
vertical_interp_X(field1.vec_d.data(), field2.vec_d.data(), hyb_press.data(), &vert_index[0], plev.data(), nplev,
gridsize, nlevels, missval);
}
}
void *
Vertintap(void *process)
{
......@@ -145,7 +105,7 @@ Vertintap(void *process)
};
int nrecs;
int varID, levelID;
int nhlev = 0, nhlevf = 0, nhlevh = 0, nlevel;
int nhlev = 0, nhlevf = 0, nhlevh = 0;
int apressID = -1, dpressID = -1;
int psID = -1, tempID = -1;
char stdname[CDI_MAX_NAME];
......@@ -251,24 +211,24 @@ Vertintap(void *process)
if (zaxisID == varList1[apressID].zaxisID)
{
bool mono_level = true;
nlevel = zaxisInqSize(zaxisID);
int nlevels = zaxisInqSize(zaxisID);
if (is_height_axis(zaxisID, nlevel))
if (is_height_axis(zaxisID, nlevels))
{
Varray<double> level(nlevel);
Varray<double> level(nlevels);
cdoZaxisInqLevels(zaxisID, &level[0]);
int l;
for (l = 0; l < nlevel; l++)
for (l = 0; l < nlevels; l++)
{
if ((l + 1) != (int) (level[l] + 0.5)) break;
}
if (l == nlevel) mono_level = true;
if (l == nlevels) mono_level = true;
}
if (is_height_axis(zaxisID, nlevel) && mono_level)
if (is_height_axis(zaxisID, nlevels) && mono_level)
{
zaxisIDh = zaxisID;
nhlev = nlevel;
nhlev = nlevels;
nhlevf = nhlev;
nhlevh = nhlevf + 1;
......@@ -311,7 +271,7 @@ Vertintap(void *process)
full_press.init(var3Df);
var3Dh.gridsize = gridsize;
var3Dh.nlevels = nhlevf;
var3Dh.nlevels = nhlevh;
var3Dh.memType = memtype;
half_press.init(var3Dh);
}
......@@ -336,14 +296,14 @@ Vertintap(void *process)
{
const auto gridID = varList1[varID].gridID;
const auto zaxisID = varList1[varID].zaxisID;
const auto nlevel = varList1[varID].nlevels;
const auto nlevels = varList1[varID].nlevels;
if (gridInqType(gridID) == GRID_SPECTRAL) cdoAbort("Spectral data unsupported!");
vardata1[varID].init(varList1[varID]);
varinterp[varID]
= (zaxisID == zaxisIDh || (is_height_axis(zaxisID, nlevel) && zaxisIDh != -1 && (nlevel == nhlevh || nlevel == nhlevf)));
= (zaxisID == zaxisIDh || (is_height_axis(zaxisID, nlevels) && zaxisIDh != -1 && (nlevels == nhlevh || nlevels == nhlevf)));
if (varinterp[varID])
{
......@@ -352,10 +312,10 @@ Vertintap(void *process)
}
else
{
if (is_height_axis(zaxisID, nlevel) && zaxisIDh != -1 && nlevel > 1)
cdoWarning("Parameter %d has wrong number of levels, skipped! (param=%s nlevel=%d)", varID + 1, varList1[varID].name, nlevel);
if (is_height_axis(zaxisID, nlevels) && zaxisIDh != -1 && nlevels > 1)
cdoWarning("Parameter %d has wrong number of levels, skipped! (param=%s nlevel=%d)", varID + 1, varList1[varID].name, nlevels);
varnmiss[varID].resize(nlevel);
varnmiss[varID].resize(nlevels);
}
}
......@@ -382,8 +342,8 @@ Vertintap(void *process)
for (varID = 0; varID < nvars; ++varID)
{
vars[varID] = false;
nlevel = varList1[varID].nlevels;
for (levelID = 0; levelID < nlevel; levelID++) varnmiss[varID][levelID] = 0;
int nlevels = varList1[varID].nlevels;
for (levelID = 0; levelID < nlevels; levelID++) varnmiss[varID][levelID] = 0;
}
taxisCopyTimestep(taxisID2, taxisID1);
......@@ -391,9 +351,9 @@ Vertintap(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
vars[varID] = true;
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, vardata1[varID], levelID, &varnmiss[varID][levelID]);
vars[varID] = true;
}
for (varID = 0; varID < nvars; varID++)
......
src/Vertintml.cc
View file @ 71e226d2
......@@ -27,7 +27,7 @@
#include "cdo_options.h"
#include "process_int.h"
#include "cdo_vlist.h"
#include "vertical_interp.h"
#include "field_vinterp.h"
#include "stdnametable.h"
#include "constants.h"
#include "util_string.h"
......@@ -183,12 +183,14 @@ Vertintml(void *process)
VarList varList1, varList2;
varListInit(varList1, vlistID1);
varListInit(varList2, vlistID2);
varListSetUniqueMemtype(varList1);
const auto memtype = varList1[0].memType;
const auto nvars = vlistNvars(vlistID1);
std::vector<bool> vars(nvars), varinterp(nvars);
std::vector<std::vector<size_t>> varnmiss(nvars);
Varray2D<double> vardata1(nvars), vardata2(nvars);
Field3DVector vardata1(nvars), vardata2(nvars);
const int maxlev = nhlevh > nplev ? nhlevh : nplev;
......@@ -203,13 +205,23 @@ Vertintml(void *process)
}
std::vector<int> vert_index;
Varray<double> ps_prog, full_press, half_press;
Field3D full_press, half_press;
if (zaxisIDh != -1 && gridsize > 0)
{
vert_index.resize(gridsize * nplev);
ps_prog.resize(gridsize);
full_press.resize(gridsize * nhlevf);
half_press.resize(gridsize * nhlevh);
CdoVar var3Df, var3Dh;
var3Df.gridsize = gridsize;
var3Df.nlevels = nhlevf;
var3Df.memType = memtype;
full_press.init(var3Df);
var3Dh.gridsize = gridsize;
var3Dh.nlevels = nhlevh;
var3Dh.memType = memtype;
half_press.init(var3Dh);
}
else
cdoWarning("No 3D variable with hybrid sigma pressure coordinate found!");
......@@ -246,8 +258,7 @@ Vertintml(void *process)
const auto gridID = varList1[varID].gridID;
const auto zaxisID = varList1[varID].zaxisID;
const auto zaxistype = zaxisInqType(zaxisID);
// gridsize = gridInqSize(gridID);
const auto nlevel = zaxisInqSize(zaxisID);
const auto nlevels = varList1[varID].nlevels;
const auto instNum = institutInqCenter(vlistInqVarInstitut(vlistID1, varID));
const auto tableNum = tableInqNum(vlistInqVarTable(vlistID1, varID));
......@@ -272,9 +283,8 @@ Vertintml(void *process)
echam_gribcodes(&gribcodes);
}
// KNMI: HIRLAM model version 7.2 uses tableNum=1 (LAMH_D11*)
// KNMI: HARMONIE model version 36 uses tableNum=1 (grib*)
// (opreational NWP version) KNMI: HARMONIE model version 38 uses
// tableNum=253 (grib,grib_md) and tableNum=1 (grib_sfx) (research version)
// KNMI: HARMONIE model version 36 uses tableNum=1 (grib*) (opreational NWP version)
// KNMI: HARMONIE model version 38 uses tableNum=253 (grib,grib_md) and tableNum=1 (grib_sfx) (research version)
else if (tableNum == 1 || tableNum == 253)
{
mode = ModelMode::HIRLAM;
......@@ -320,21 +330,21 @@ Vertintml(void *process)
if (mode == ModelMode::ECHAM)
{
// clang-format off
if (code == gribcodes.geopot && nlevel == 1) sgeopotID = varID;
else if (code == gribcodes.geopot && nlevel == nhlevf) geopotID = varID;
else if (code == gribcodes.temp && nlevel == nhlevf) tempID = varID;
else if (code == gribcodes.ps && nlevel == 1) psID = varID;
else if (code == gribcodes.lsp && nlevel == 1) lnpsID = varID;
else if (code == gribcodes.gheight && nlevel == nhlevf) gheightID = varID;
if (code == gribcodes.geopot && nlevels == 1) sgeopotID = varID;
else if (code == gribcodes.geopot && nlevels == nhlevf) geopotID = varID;
else if (code == gribcodes.temp && nlevels == nhlevf) tempID = varID;
else if (code == gribcodes.ps && nlevels == 1) psID = varID;
else if (code == gribcodes.lsp && nlevels == 1) lnpsID = varID;
else if (code == gribcodes.gheight && nlevels == nhlevf) gheightID = varID;
// clang-format on
}
else if (mode == ModelMode::WMO || mode == ModelMode::HIRLAM)
{
// clang-format off
if (code == gribcodes.geopot && nlevel == 1) sgeopotID = varID;
else if (code == gribcodes.geopot && nlevel == nhlevf) geopotID = varID;
else if (code == gribcodes.temp && nlevel == nhlevf) tempID = varID;
else if (code == gribcodes.ps && nlevel == 1) psID = varID;
if (code == gribcodes.geopot && nlevels == 1) sgeopotID = varID;
else if (code == gribcodes.geopot && nlevels == nhlevf) geopotID = varID;
else if (code == gribcodes.temp && nlevels == nhlevf) tempID = varID;
else if (code == gribcodes.ps && nlevels == 1) psID = varID;
// clang-format on
}
......@@ -342,27 +352,27 @@ Vertintml(void *process)
if (gridInqType(gridID) == GRID_SPECTRAL) cdoAbort("Spectral data unsupported!");
const size_t maxlevel = (varID == gheightID) ? (nlevel + 1) : nlevel;
vardata1[varID].resize(gridsize * maxlevel);
if (varID == gheightID) varList1[varID].nlevels = nlevels + 1;
vardata1[varID].init(varList1[varID]);
// varinterp[varID] = ( zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && nlevel == nhlev );
// varinterp[varID] = ( zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && nlevels == nhlev );
varinterp[varID]
= (zaxisID == zaxisIDh || (zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && (nlevel == nhlevh || nlevel == nhlevf)));
= (zaxisID == zaxisIDh || (zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && (nlevels == nhlevh || nlevels == nhlevf)));
if (varinterp[varID])
{
vardata2[varID].resize(gridsize * nplev);
varnmiss[varID].resize(maxlev, 0);
vardata2[varID].init(varList2[varID]);
}
else
{
if (zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && nlevel > 1)
if (zaxis_is_hybrid(zaxistype) && zaxisIDh != -1 && nlevels > 1)
{
vlistInqVarName(vlistID1, varID, varname);
cdoWarning("Parameter %d has wrong number of levels, skipped! (param=%s nlevel=%d)", varID + 1, varname, nlevel);
cdoWarning("Parameter %d has wrong number of levels, skipped! (param=%s nlevel=%d)", varID + 1, varname, nlevels);
}
varnmiss[varID].resize(nlevel);
varnmiss[varID].resize(nlevels);
}
}
......@@ -382,23 +392,6 @@ Vertintml(void *process)
if (tempID != -1 || gheightID != -1) sgeopot_needed = true;
Varray<double> sgeopot;
if (zaxisIDh != -1 && sgeopot_needed)
{
sgeopot.resize(gridsize);
if (sgeopotID == -1)
{
if (extrapolate)
{
if (geopotID == -1)
cdoWarning("%s not found - set to zero!", var_stdname(surface_geopotential));
else
cdoPrint("%s not found - using bottom layer of %s!", var_stdname(surface_geopotential), var_stdname(geopotential));
}
varrayFill(sgeopot, 0.0);
}
}
if (zaxisIDh != -1 && gheightID != -1 && tempID == -1)
cdoAbort("Temperature not found, needed for vertical interpolation of geopotheight!");
......@@ -415,11 +408,30 @@ Vertintml(void *process)
if (Options::cdoVerbose)
{
vlistInqVarName(vlistID1, presID, varname);
if (presID == lnpsID)
cdoPrint("using LOG(%s) from %s", var_stdname(surface_air_pressure), varname);
cdoPrint("using LOG(%s) from %s", var_stdname(surface_air_pressure), varList1[presID].name);
else
cdoPrint("using %s from %s", var_stdname(surface_air_pressure), varname);
cdoPrint("using %s from %s", var_stdname(surface_air_pressure), varList1[presID].name);
}
Field ps_prog;
ps_prog.init(varList1[presID]);
Field sgeopot;
if (zaxisIDh != -1 && sgeopot_needed)
{
sgeopot.init(varList1[presID]);
if (sgeopotID == -1)
{
if (extrapolate)
{
if (geopotID == -1)
cdoWarning("%s not found - set to zero!", var_stdname(surface_geopotential));
else
cdoPrint("%s not found - using bottom layer of %s!", var_stdname(surface_geopotential), var_stdname(geopotential));
}
fieldFill(sgeopot, 0.0);
}
}
// check VCT
......@@ -445,15 +457,13 @@ Vertintml(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID1, &varID, &levelID);
// gridsize = varList1[varID].gridsize;
const auto zaxisID = varList1[varID].zaxisID;
const auto nlevel = varList1[varID].nlevels;
const auto nlevels = varList1[varID].nlevels;
if (linvertvct && zaxisIDh != -1 && zaxisID == zaxisIDh) levelID = nlevels - 1 - levelID;
if (linvertvct && zaxisIDh != -1 && zaxisID == zaxisIDh) levelID = nlevel - 1 - levelID;
cdoReadRecord(streamID1, vardata1[varID], levelID, &varnmiss[varID][levelID]);
const auto offset = gridsize * levelID;
auto single = &vardata1[varID][offset];
cdoReadRecord(streamID1, single, &varnmiss[varID][levelID]);
vars[varID] = true;
}
......@@ -462,14 +472,14 @@ Vertintml(void *process)
if (sgeopot_needed)
{
if (sgeopotID != -1)
varrayCopy(gridsize, vardata1[sgeopotID], sgeopot);
fieldCopy(vardata1[sgeopotID], sgeopot);
else if (geopotID != -1)
arrayCopy(gridsize, &vardata1[geopotID][gridsize * (nhlevf - 1)], &sgeopot[0]);
fieldCopy(vardata1[geopotID], nhlevf - 1, sgeopot);
// check range of surface geopot
if (extrapolate && (sgeopotID != -1 || geopotID != -1))
{
const auto minmax = varrayMinMax(sgeopot);
const auto minmax = fieldMinMax(sgeopot);
if (minmax.first < MIN_FIS || minmax.second > MAX_FIS)
cdoWarning("Surface geopotential out of range (min=%g max=%g) [timestep:%d]!", minmax.first, minmax.second,
tsID + 1);
......@@ -480,102 +490,93 @@ Vertintml(void *process)
}
if (presID == lnpsID)
for (size_t i = 0; i < gridsize; i++) ps_prog[i] = std::exp(vardata1[lnpsID][i]);
{
if (memtype == MemType::Float)
for (size_t i = 0; i < gridsize; i++) ps_prog.vec_f[i] = std::exp(vardata1[lnpsID].vec_f[i]);
else
for (size_t i = 0; i < gridsize; i++) ps_prog.vec_d[i] = std::exp(vardata1[lnpsID].vec_d[i]);
}
else if (presID != -1)
varrayCopy(gridsize, vardata1[presID], ps_prog);
fieldCopy(vardata1[presID], ps_prog);
// check range of ps_prog
const auto minmax = varrayMinMax(ps_prog);
const auto minmax = fieldMinMax(ps_prog);
if (minmax.first < MIN_PS || minmax.second > MAX_PS)
cdoWarning("Surface pressure out of range (min=%g max=%g) [timestep:%d]!", minmax.first, minmax.second, tsID + 1);
presh(full_press.data(), half_press.data(), vct.data(), ps_prog.data(), nhlevf, gridsize);
if (memtype == MemType::Float)
presh(full_press.vec_f.data(), half_press.vec_f.data(), vct.data(), ps_prog.vec_f.data(), nhlevf, gridsize);
else
presh(full_press.vec_d.data(), half_press.vec_d.data(), vct.data(), ps_prog.vec_d.data(), nhlevf, gridsize);
if (useLogType)
{
for (size_t i = 0; i < gridsize; i++) ps_prog[i] = std::log(ps_prog[i]);
for (size_t ki = 0; ki < nhlevh * gridsize; ki++) half_press[ki] = std::log(half_press[ki]);
for (size_t ki = 0; ki < nhlevf * gridsize; ki++) full_press[ki] = std::log(full_press[ki]);
if (memtype == MemType::Float)
{
for (size_t i = 0; i < gridsize; i++) ps_prog.vec_f[i] = std::log(ps_prog.vec_f[i]);
for (size_t ki = 0; ki < nhlevh * gridsize; ki++) half_press.vec_f[ki] = std::log(half_press.vec_f[ki]);
for (size_t ki = 0; ki < nhlevf * gridsize; ki++) full_press.vec_f[ki] = std::log(full_press.vec_f[ki]);
}
else
{
for (size_t i = 0; i < gridsize; i++) ps_prog.vec_d[i] = std::log(ps_prog.vec_d[i]);
for (size_t ki = 0; ki < nhlevh * gridsize; ki++) half_press.vec_d[ki] = std::log(half_press.vec_d[ki]);
for (size_t ki = 0; ki < nhlevf * gridsize; ki++) full_press.vec_d[ki] = std::log(full_press.vec_d[ki]);
}
}
genind(vert_index.data(), plev.data(), full_press.data(), gridsize, nplev, nhlevf);
genind(vert_index, plev, full_press, gridsize);
if (!extrapolate) genindmiss(vert_index.data(), plev.data(), gridsize, nplev, ps_prog.data(), pnmiss.data());
if (!extrapolate) genindmiss(vert_index, plev, gridsize, ps_prog, pnmiss);
}
for (varID = 0; varID < nvars; varID++)
{
if (vars[varID])
{
const auto nlevel = varList1[varID].nlevels;
const auto missval = varList1[varID].missval;
if (varinterp[varID])
{
/*
if ( nlevel == nhlevh )
{
for ( int k = 1; k < nlevel; k++ )
{
double *vl1 = &vardata1[varID][gridsize*(k-1)];
double *vl2 = &vardata1[varID][gridsize*(k)];
for ( size_t i = 0; i < gridsize; i++ )
vl1[i] = 0.5*(vl1[i] + vl2[i]);
}
nlevel = nhlevf;
}
*/
double *hyb_press = (nlevel == nhlevf) ? &full_press[0] : (nlevel == nhlevh) ? &half_press[0] : nullptr;
if (hyb_press == nullptr)
{
vlistInqVarName(vlistID1, varID, varname);
cdoAbort("Number of hybrid level differ from full/half level (param=%s)!", varname);
}
const auto nlevels = varList1[varID].nlevels;
if (nlevels != nhlevf && nlevels != nhlevh)
cdoAbort("Number of hybrid level differ from full/half level (param=%s)!", varList1[varID].name);
for (levelID = 0; levelID < nlevel; levelID++)
for (levelID = 0; levelID < nlevels; levelID++)
{
if (varnmiss[varID][levelID]) cdoAbort("Missing values unsupported for this operator!");
}
if (varID == tempID)
{
if (nlevel == nhlevh) cdoAbort("Temperature on half level unsupported!");
if (nlevels == nhlevh) cdoAbort("Temperature on half level unsupported!");
if (useLogType && extrapolate) cdoAbort("Log. extrapolation of temperature unsupported!");
vertical_interp_T(&sgeopot[0], &vardata1[varID][0], &vardata2[varID][0], &full_press[0], &half_press[0],
&vert_index[0], plev.data(), nplev, gridsize, nlevel, missval);
vertical_interp_T(nlevels, full_press, half_press, vardata1[varID], vardata2[varID], sgeopot,
vert_index, plev, gridsize);
}
else if (varID == gheightID)
{
for (size_t i = 0; i < gridsize; ++i) vardata1[varID][gridsize * nlevel + i] = sgeopot[i] / PlanetGrav;
if (memtype == MemType::Float)
for (size_t i = 0; i < gridsize; ++i) vardata1[varID].vec_f[gridsize * nlevels + i] = sgeopot.vec_f[i] / PlanetGrav;
else
for (size_t i = 0; i < gridsize; ++i) vardata1[varID].vec_d[gridsize * nlevels + i] = sgeopot.vec_d[i] / PlanetGrav;
vertical_interp_Z(&sgeopot[0], &vardata1[varID][0], &vardata2[varID][0], &full_press[0], &half_press[0],
&vert_index[0], &vardata1[tempID][0], plev.data(), nplev, gridsize, nlevel, missval);
vertical_interp_Z(nlevels, full_press, half_press, vardata1[varID], vardata2[varID], vardata1[tempID], sgeopot,
vert_index, plev, gridsize);
}
else
{
vertical_interp_X(&vardata1[varID][0], &vardata2[varID][0], hyb_press, &vert_index[0], plev.data(), nplev,
gridsize, nlevel, missval);
vertical_interp_X(nlevels, full_press, half_press, vardata1[varID], vardata2[varID], vert_index, plev, gridsize);
}
if (!extrapolate) varrayCopy(nplev, pnmiss, varnmiss[varID]);
}
}
}
for (varID = 0; varID < nvars; varID++)
{
if (vars[varID])
{
const auto gridsize = varList2[varID].gridsize;
const auto nlevel = varList2[varID].nlevels;