Skip to content

GitLab

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

Derivepar: cleanup.

parent 91547b16
Hide whitespace changes
Inline Side-by-side
src/Derivepar.cc
View file @ d058ea89
......@@ -17,7 +17,8 @@
/*
This module contains the following operators:
Derivepar gheight geopotential height
Derivepar gheight Geopotential height
Derivepar sealevelpressure Sea level pressure
*/
#include <cdi.h>
......@@ -39,8 +40,6 @@ void *
Derivepar(void *process)
{
ModelMode mode(ModelMode::UNDEF);
int nrecs;
size_t offset;
int varID, levelID;
int surfaceID = -1;
int sgeopotID = -1, geopotID = -1, tempID = -1, humID = -1, psID = -1, lnpsID = -1, presID = -1, gheightID = -1;
......@@ -48,7 +47,6 @@ Derivepar(void *process)
// int clwcID = -1, ciwcID = -1;
char paramstr[32];
char varname[CDI_MAX_NAME], stdname[CDI_MAX_NAME];
double *single2;
// double *lwater = nullptr, *iwater = nullptr;
size_t nmiss, nmissout = 0;
gribcode_t gribcodes;
......@@ -83,12 +81,16 @@ Derivepar(void *process)
if (zaxisIDh == -1) cdoAbort("No 3D variable with hybrid sigma pressure coordinate found!");
VarList varList1;
varListInit(varList1, vlistID1);
varListSetUniqueMemtype(varList1);
const auto nvars = vlistNvars(vlistID1);
auto useTable = false;
for (varID = 0; varID < nvars; varID++)
{
int tableNum = tableInqNum(vlistInqVarTable(vlistID1, varID));
const auto tableNum = tableInqNum(vlistInqVarTable(vlistID1, varID));
if (tableNum > 0 && tableNum != 255)
{
useTable = true;
......@@ -100,15 +102,15 @@ Derivepar(void *process)
for (varID = 0; varID < nvars; varID++)
{
const auto gridID = vlistInqVarGrid(vlistID1, varID);
const auto zaxisID = vlistInqVarZaxis(vlistID1, varID);
const auto nlevel = zaxisInqSize(zaxisID);
const auto gridID = varList1[varID].gridID;
const auto zaxisID = varList1[varID].zaxisID;
const auto nlevels = varList1[varID].nlevels;
const auto instNum = institutInqCenter(vlistInqVarInstitut(vlistID1, varID));
const auto tableNum = tableInqNum(vlistInqVarTable(vlistID1, varID));
auto code = vlistInqVarCode(vlistID1, varID);
auto param = vlistInqVarParam(vlistID1, varID);
auto code = varList1[varID].code;
const auto param = varList1[varID].param;
cdiParamToString(param, paramstr, sizeof(paramstr));
int pnum, pcat, pdis;
cdiDecodeParam(param, &pnum, &pcat, &pdis);
......@@ -158,12 +160,12 @@ Derivepar(void *process)
{
// clang-format off
if (sgeopotID == -1 && cdo_cmpstr(varname, "geosp")) code = gribcodes.geopot;
else if (psID == -1 && cdo_cmpstr(varname, "aps")) code = gribcodes.ps;
else if (psID == -1 && cdo_cmpstr(varname, "ps")) code = gribcodes.ps;
else if (lnpsID == -1 && cdo_cmpstr(varname, "lsp")) code = gribcodes.lsp;
else if (lnpsID2 == -1 && cdo_cmpstr(varname, "lnps")) code = 777;
else if (tempID == -1 && cdo_cmpstr(varname, "t")) code = gribcodes.temp;
else if (humID == -1 && cdo_cmpstr(varname, "q")) code = gribcodes.hum;
else if (psID == -1 && cdo_cmpstr(varname, "aps")) code = gribcodes.ps;
else if (psID == -1 && cdo_cmpstr(varname, "ps")) code = gribcodes.ps;
else if (lnpsID == -1 && cdo_cmpstr(varname, "lsp")) code = gribcodes.lsp;
else if (lnpsID2 == -1 && cdo_cmpstr(varname, "lnps")) code = 777;
else if (tempID == -1 && cdo_cmpstr(varname, "t")) code = gribcodes.temp;
else if (humID == -1 && cdo_cmpstr(varname, "q")) code = gribcodes.hum;
// else if ( geopotID == -1 && cdo_cmpstr(stdname, "geopotential_full")) code = gribcodes.geopot;
// else if (cdo_cmpstr(varname, "clwc")) code = 246;
// else if (cdo_cmpstr(varname, "ciwc")) code = 247;
......@@ -172,14 +174,14 @@ Derivepar(void *process)
}
// 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.hum && nlevel == nhlevf) humID = varID;
else if (code == gribcodes.ps && nlevel == 1) psID = varID;
else if (code == gribcodes.lsp && nlevel == 1) lnpsID = varID;
else if (code == 777 && nlevel == 1) lnpsID2 = 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.hum && nlevels == nhlevf) humID = varID;
else if (code == gribcodes.ps && nlevels == 1) psID = varID;
else if (code == gribcodes.lsp && nlevels == 1) lnpsID = varID;
else if (code == 777 && nlevels == 1) lnpsID2 = varID;
else if (code == gribcodes.gheight && nlevels == nhlevf) gheightID = varID;
// else if (code == 246) clwcID = varID;
// else if (code == 247) ciwcID = varID;
// clang-format on
......@@ -195,12 +197,14 @@ Derivepar(void *process)
if (Options::cdoVerbose)
{
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));
// clang-format off
if (-1 != tempID) cdoPrint(" %s -> %s", var_stdname(air_temperature), varList1[tempID].name);
if (-1 != psID) cdoPrint(" %s -> %s", var_stdname(surface_air_pressure), varList1[psID].name);
if (-1 != lnpsID) cdoPrint(" LOG(%s) -> %s", var_stdname(surface_air_pressure), varList1[lnpsID].name);
if (-1 != sgeopotID) cdoPrint(" %s -> %s", var_stdname(surface_geopotential), varList1[sgeopotID].name);
if (-1 != geopotID) cdoPrint(" %s -> %s", var_stdname(geopotential), varList1[geopotID].name);
if (-1 != gheightID) cdoPrint(" %s -> %s", var_stdname(geopotential_height), varList1[gheightID].name);
// clang-format on
}
if (lnpsID != -1 && lnpsID2 != -1) cdoAbort("Found LOG(%s) twice: lsp and lnps!", var_stdname(surface_air_pressure));
......@@ -276,7 +280,7 @@ Derivepar(void *process)
cdoAbort("Internal problem, invalid operatorID: %d!", operatorID);
vlistDefVarParam(vlistID2, varID, cdiEncodeParam(var_echamcode(var_id), 128, 255));
cdiDefKeyString(vlistID2, varID, CDI_KEY_STDNAME, var_name(var_id));
cdiDefKeyString(vlistID2, varID, CDI_KEY_NAME, var_name(var_id));
cdiDefKeyString(vlistID2, varID, CDI_KEY_STDNAME, var_stdname(var_id));
cdiDefKeyString(vlistID2, varID, CDI_KEY_UNITS, var_units(var_id));
......@@ -288,17 +292,21 @@ Derivepar(void *process)
cdoDefVlist(streamID2, vlistID2);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
taxisCopyTimestep(taxisID2, taxisID1);
cdoDefTimestep(streamID2, tsID);
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID1, &varID, &levelID);
offset = gridsize * levelID;
cdoReadRecord(streamID1, array.data(), &nmiss);
const auto offset = gridsize * levelID;
if (zaxisIDh != -1)
{
if (varID == sgeopotID)
......@@ -337,14 +345,10 @@ Derivepar(void *process)
}
{
varID = tempID;
const auto nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for (levelID = 0; levelID < nlevel; levelID++)
for (levelID = 0; levelID < varList1[tempID].nlevels; levelID++)
{
offset = gridsize * levelID;
single2 = &temp[offset];
const auto mm = varrayMinMax(gridsize, single2);
auto single = &temp[gridsize * levelID];
const auto mm = varrayMinMax(gridsize, single);
if (mm.min < MIN_T || mm.max > MAX_T)
cdoWarning("Input temperature at level %d out of range (min=%g max=%g)!", levelID + 1, mm.min, mm.max);
}
......@@ -352,16 +356,11 @@ Derivepar(void *process)
if (humID != -1)
{
varID = humID;
const auto nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for (levelID = 0; levelID < nlevel; levelID++)
for (levelID = 0; levelID < varList1[humID].nlevels; levelID++)
{
offset = gridsize * levelID;
single2 = &hum[offset];
// corr_hum(gridsize, single2, MIN_Q);
const auto mm = varrayMinMax(gridsize, single2);
auto single = &hum[gridsize * levelID];
// corr_hum(gridsize, single, MIN_Q);
const auto mm = varrayMinMax(gridsize, single);
if (mm.min < -0.1 || mm.max > MAX_Q)
cdoWarning("Input humidity at level %d out of range (min=%g max=%g)!", levelID + 1, mm.min, mm.max);
}
......@@ -375,8 +374,8 @@ Derivepar(void *process)
nmissout = 0;
varID = 0;
const auto nlevel = nhlevf;
for (levelID = 0; levelID < nlevel; levelID++)
const auto nlevels = nhlevf;
for (levelID = 0; levelID < nlevels; levelID++)
{
cdoDefRecord(streamID2, varID, levelID);
cdoWriteRecord(streamID2, gheight.data() + levelID * gridsize, nmissout);
......
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

Imprint | Privacy Policy