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

Afterburner: cleanup.

parent 7453790c
Pipeline #4819 failed with stages
in 9 minutes and 19 seconds
...@@ -414,20 +414,16 @@ after_readTimestep(void *arg) ...@@ -414,20 +414,16 @@ after_readTimestep(void *arg)
const auto code = vlistInqVarCode(globs->ivlistID, varID); const auto code = vlistInqVarCode(globs->ivlistID, varID);
if (code <= 0 || code >= MaxCodes) continue; if (code <= 0 || code >= MaxCodes) continue;
/* Skip records containing unneeded codes */ // Skip records containing unneeded codes
if (!vars[code].needed0) continue; if (!vars[code].needed0) continue;
vlistInqVar(globs->ivlistID, varID, &gridID, &zaxisID, &timeID); vlistInqVar(globs->ivlistID, varID, &gridID, &zaxisID, &timeID);
const auto leveltype = zaxisInqType(zaxisID); const auto leveltype = zaxisInqType(zaxisID);
/* Skip records with unselected levels */ // Skip records with unselected levels
int levelOffset = -1; int levelOffset = -1;
/* // if ( vars[code].ozaxisID != vars[code].izaxisID && ! Lhybrid2pressure )
if ( vars[code].ozaxisID != vars[code].izaxisID && ! Lhybrid2pressure )
*/
if ((vars[code].ozaxisID != vars[code].izaxisID) && (leveltype == ZAXIS_PRESSURE)) if ((vars[code].ozaxisID != vars[code].izaxisID) && (leveltype == ZAXIS_PRESSURE))
{ {
const auto level = (int) zaxisInqLevel(zaxisID, levelID); const auto level = (int) zaxisInqLevel(zaxisID, levelID);
...@@ -527,25 +523,23 @@ after_setEndOfInterval(AfterControl &globs, int nrecs) ...@@ -527,25 +523,23 @@ after_setEndOfInterval(AfterControl &globs, int nrecs)
static void static void
after_moveTimestep(struct Variable *vars) after_moveTimestep(struct Variable *vars)
{ {
int code; for (int code = 0; code < MaxCodes; code++) vars[code].nmiss = vars[code].nmiss0;
for (code = 0; code < MaxCodes; code++) vars[code].nmiss = vars[code].nmiss0;
for (code = 0; code < MaxCodes; code++) for (int code = 0; code < MaxCodes; code++)
if (vars[code].hybrid0) if (vars[code].hybrid0)
{ {
vars[code].hybrid = vars[code].hybrid0; vars[code].hybrid = vars[code].hybrid0;
vars[code].hybrid0 = nullptr; vars[code].hybrid0 = nullptr;
} }
for (code = 0; code < MaxCodes; code++) for (int code = 0; code < MaxCodes; code++)
if (vars[code].spectral0) if (vars[code].spectral0)
{ {
vars[code].spectral = vars[code].spectral0; vars[code].spectral = vars[code].spectral0;
vars[code].spectral0 = nullptr; vars[code].spectral0 = nullptr;
} }
for (code = 0; code < MaxCodes; code++) for (int code = 0; code < MaxCodes; code++)
if (vars[code].grid0) if (vars[code].grid0)
{ {
vars[code].grid = vars[code].grid0; vars[code].grid = vars[code].grid0;
...@@ -558,7 +552,7 @@ after_check_content(struct Variable *vars, int timestep) ...@@ -558,7 +552,7 @@ after_check_content(struct Variable *vars, int timestep)
{ {
for (int code = 0; code < 272; code++) for (int code = 0; code < 272; code++)
{ {
/* if ( code == GEOPOTENTIAL ) continue; */ // if ( code == GEOPOTENTIAL ) continue;
if (code == SLP) continue; if (code == SLP) continue;
if (code == GEOPOTHEIGHT) continue; if (code == GEOPOTHEIGHT) continue;
if (code == STREAM) continue; if (code == STREAM) continue;
...@@ -1845,7 +1839,7 @@ after_processing(AfterControl &globs, struct Variable *vars) ...@@ -1845,7 +1839,7 @@ after_processing(AfterControl &globs, struct Variable *vars)
after_defineGrid(globs, vars); after_defineGrid(globs, vars);
after_postcntl(globs, vars); /* define output variables */ after_postcntl(globs, vars); // define output variables
after_control(globs, vars); after_control(globs, vars);
...@@ -1885,7 +1879,7 @@ Afterburner(void *process) ...@@ -1885,7 +1879,7 @@ Afterburner(void *process)
struct Variable vars[MaxCodes + 5]; struct Variable vars[MaxCodes + 5];
for (int code = 0; code < MaxCodes + 5; code++) after_variable_init(&vars[code]); for (int code = 0; code < MaxCodes + 5; code++) after_variable_init(&vars[code]);
after_parini(globs, vars); /* read namelist parameter */ after_parini(globs, vars); // read namelist parameter
if (CdoDefault::FileType != CDI_UNDEFID) ofiletype = CdoDefault::FileType; if (CdoDefault::FileType != CDI_UNDEFID) ofiletype = CdoDefault::FileType;
......
...@@ -19,22 +19,9 @@ ...@@ -19,22 +19,9 @@
#include "cdi.h" #include "cdi.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#ifdef CDO
#include "cimdOmp.h" #include "cimdOmp.h"
#define streamDefRecord cdoDefRecord #define streamDefRecord cdoDefRecord
#define streamWriteRecord cdoWriteRecord #define streamWriteRecord cdoWriteRecord
#else
#define OPENMP4 201307
#ifdef _OPENMP &&defined(OPENMP4) && _OPENMP >= OPENMP4
#define HAVE_OPENMP4 1
#endif
#endif
#endif
#include "afterburner.h" #include "afterburner.h"
#include "constants.h" #include "constants.h"
...@@ -132,10 +119,7 @@ VarQuaSum(double *Variance, const double *restrict Sum, int len, int n) ...@@ -132,10 +119,7 @@ VarQuaSum(double *Variance, const double *restrict Sum, int len, int n)
for (int i = 0; i < len; i++) Variance[i] = (Variance[i] - Sum[i] * Sum[i] * n) * rn1; for (int i = 0; i < len; i++) Variance[i] = (Variance[i] - Sum[i] * Sum[i] * n) * rn1;
for (int i = 0; i < len; i++) for (int i = 0; i < len; i++)
if (Variance[i] > 0.0) Variance[i] = (Variance[i] > 0.0) ? std::sqrt(Variance[i]) : 0.0;
Variance[i] = std::sqrt(Variance[i]);
else
Variance[i] = 0.0;
} }
static void static void
...@@ -171,12 +155,7 @@ MultVectorScalar(double *dest, const double *restrict src, double factor, size_t ...@@ -171,12 +155,7 @@ MultVectorScalar(double *dest, const double *restrict src, double factor, size_t
if (nmiss) if (nmiss)
{ {
for (size_t i = 0; i < len; i++) for (size_t i = 0; i < len; i++)
{ dest[i] = (IS_EQUAL(src[i], missval)) ? missval : src[i] * factor;
if (IS_EQUAL(src[i], missval))
dest[i] = missval;
else
dest[i] = src[i] * factor;
}
} }
else else
{ {
...@@ -288,7 +267,7 @@ after_FC2GP(double *fc, double *gp, long nlat, long nlon, long nlev, long nfc) ...@@ -288,7 +267,7 @@ after_FC2GP(double *fc, double *gp, long nlat, long nlon, long nlev, long nfc)
fc2gp(trig, ifax, fc, gp, nlat, nlon, nlev, nfc); fc2gp(trig, ifax, fc, gp, nlat, nlon, nlev, nfc);
} }
/* HUMTEST */ // HUMTEST
static void static void
sh2rh(int AnalysisData, double *sphum, double *rhum, double *t, int lev, int dimgpout, const double *level, double *fullpresshybrid) sh2rh(int AnalysisData, double *sphum, double *rhum, double *t, int lev, int dimgpout, const double *level, double *fullpresshybrid)
...@@ -349,14 +328,13 @@ rh2sh(double *sphum, double *rhum, double *t, int lev, int dimgpout, const doubl ...@@ -349,14 +328,13 @@ rh2sh(double *sphum, double *rhum, double *t, int lev, int dimgpout, const doubl
int lp, i; int lp, i;
int lpi; int lpi;
double es, qsat; double es, qsat;
double RALP, RBET, RGAM;
/* ***************************************************** */ /* ***************************************************** */
/* Define constants for calculation in presence of water */ /* Define constants for calculation in presence of water */
/* ***************************************************** */ /* ***************************************************** */
const double RGAMW = (C_RCW - C_RCPV) / C_RV; constexpr double RGAMW = (C_RCW - C_RCPV) / C_RV;
const double RBETW = C_RLVTT / C_RV + RGAMW * C_RTT; constexpr double RBETW = C_RLVTT / C_RV + RGAMW * C_RTT;
const double RALPW = std::log(C_RESTT) + RBETW / C_RTT + RGAMW * std::log(C_RTT); constexpr double RALPW = std::log(C_RESTT) + RBETW / C_RTT + RGAMW * std::log(C_RTT);
/* ***************************************************** */ /* ***************************************************** */
/* Define constants for calculation in presence of ice */ /* Define constants for calculation in presence of ice */
...@@ -371,9 +349,9 @@ rh2sh(double *sphum, double *rhum, double *t, int lev, int dimgpout, const doubl ...@@ -371,9 +349,9 @@ rh2sh(double *sphum, double *rhum, double *t, int lev, int dimgpout, const doubl
/* Hint of Michael Ponater 08.10.97 */ /* Hint of Michael Ponater 08.10.97 */
/***************************************************/ /***************************************************/
RGAM = RGAMW; constexpr double RGAM = RGAMW;
RBET = RBETW; constexpr double RBET = RBETW;
RALP = RALPW; constexpr double RALP = RALPW;
for (lp = 0; lp < lev; lp++) for (lp = 0; lp < lev; lp++)
{ {
for (i = 0; i < dimgpout; i++) for (i = 0; i < dimgpout; i++)
...@@ -908,16 +886,15 @@ static void ...@@ -908,16 +886,15 @@ static void
Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, double *halfpress, double *fullpress, double *dpsdx, Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, double *halfpress, double *fullpress, double *dpsdx,
double *dpsdy, double *vct, int dimgp, int nlev) double *dpsdy, double *vct, int dimgp, int nlev)
{ {
int i, j;
double DeltaHybrid, Cterm, Pterm; double DeltaHybrid, Cterm, Pterm;
double *diver, *halfp, *fullp, *uwind, *vwind; double *diver, *halfp, *fullp, *uwind, *vwind;
double *omega = omega_in; double *omega = omega_in;
/* Compute half level part of vertical velocity */ // Compute half level part of vertical velocity
for (i = 0; i < dimgp; i++) omega[i] = 0.0; for (int i = 0; i < dimgp; i++) omega[i] = 0.0;
for (j = 0; j < nlev; j++) for (int j = 0; j < nlev; j++)
{ {
omega = omega_in + j * dimgp; omega = omega_in + j * dimgp;
halfp = halfpress + j * dimgp; halfp = halfpress + j * dimgp;
...@@ -929,7 +906,7 @@ Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, doub ...@@ -929,7 +906,7 @@ Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, doub
#ifdef HAVE_OPENMP4 #ifdef HAVE_OPENMP4
#pragma omp parallel for simd #pragma omp parallel for simd
#endif #endif
for (i = 0; i < dimgp; i++) for (int i = 0; i < dimgp; i++)
{ {
omega[i + dimgp] omega[i + dimgp]
= omega[i] - diver[i] * (halfp[i + dimgp] - halfp[i]) - DeltaHybrid * (uwind[i] * dpsdx[i] + vwind[i] * dpsdy[i]); = omega[i] - diver[i] * (halfp[i + dimgp] - halfp[i]) - DeltaHybrid * (uwind[i] * dpsdx[i] + vwind[i] * dpsdy[i]);
...@@ -938,21 +915,21 @@ Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, doub ...@@ -938,21 +915,21 @@ Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, doub
/* interpolate to full levels */ /* interpolate to full levels */
for (j = 0; j < nlev; j++) for (int j = 0; j < nlev; j++)
{ {
omega = omega_in + j * dimgp; omega = omega_in + j * dimgp;
#ifdef HAVE_OPENMP4 #ifdef HAVE_OPENMP4
#pragma omp parallel for simd #pragma omp parallel for simd
#endif #endif
for (i = 0; i < dimgp; i++) omega[i] = 0.5 * (omega[i] + omega[i + dimgp]); for (int i = 0; i < dimgp; i++) omega[i] = 0.5 * (omega[i] + omega[i + dimgp]);
} }
/* compute full level part of vertical velocity */ /* compute full level part of vertical velocity */
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel for default(shared) private(i, omega, halfp, fullp, uwind, vwind, DeltaHybrid, Cterm, Pterm) #pragma omp parallel for default(shared) private(omega, halfp, fullp, uwind, vwind, DeltaHybrid, Cterm, Pterm)
#endif #endif
for (j = 0; j < nlev; j++) for (int j = 0; j < nlev; j++)
{ {
omega = omega_in + j * dimgp; omega = omega_in + j * dimgp;
halfp = halfpress + j * dimgp; halfp = halfpress + j * dimgp;
...@@ -967,7 +944,7 @@ Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, doub ...@@ -967,7 +944,7 @@ Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, doub
#ifdef HAVE_OPENMP4 #ifdef HAVE_OPENMP4
#pragma omp simd #pragma omp simd
#endif #endif
for (i = 0; i < dimgp; i++) for (int i = 0; i < dimgp; i++)
{ {
if (IS_NOT_EQUAL(Cterm, 0.0)) if (IS_NOT_EQUAL(Cterm, 0.0))
Pterm = Cterm / (halfp[i + dimgp] - halfp[i]) * std::log(halfp[i + dimgp] / halfp[i]); Pterm = Cterm / (halfp[i + dimgp] - halfp[i]) * std::log(halfp[i + dimgp] / halfp[i]);
...@@ -984,57 +961,51 @@ Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, doub ...@@ -984,57 +961,51 @@ Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, doub
void void
MakeGeopotHeight(double *geop, double *gt, double *gq, double *ph, int nhor, int nlev) MakeGeopotHeight(double *geop, double *gt, double *gq, double *ph, int nhor, int nlev)
{ {
int i, j; constexpr double z2log2 = 2.0 * std::log(2.0);
double vtmp; const double vtmp = (C_RV / PlanetRD) - 1.0;
double zrg; const double zrg = 1.0 / PlanetGrav;
double z2log2;
double *restrict geopl, *restrict gtl, *restrict gql, *restrict phl;
z2log2 = 2.0 * std::log(2.0);
vtmp = (C_RV / PlanetRD) - 1.0;
zrg = 1.0 / PlanetGrav;
if (gq) /* Humidity is present */ if (gq) /* Humidity is present */
{ {
for (j = nlev; j > 1; j--) for (int j = nlev; j > 1; j--)
{ {
geopl = geop + nhor * (j - 1); double *geopl = geop + nhor * (j - 1);
gtl = gt + nhor * (j - 1); const double *restrict gtl = gt + nhor * (j - 1);
gql = gq + nhor * (j - 1); const double *restrict gql = gq + nhor * (j - 1);
phl = ph + nhor * (j - 1); const double *restrict phl = ph + nhor * (j - 1);
#ifdef HAVE_OPENMP4 #ifdef HAVE_OPENMP4
#pragma omp parallel for simd #pragma omp parallel for simd
#endif #endif
for (i = 0; i < nhor; i++) for (int i = 0; i < nhor; i++)
geopl[i] = geopl[i + nhor] + PlanetRD * gtl[i] * (1.0 + vtmp * gql[i]) * std::log(phl[i + nhor] / phl[i]); geopl[i] = geopl[i + nhor] + PlanetRD * gtl[i] * (1.0 + vtmp * gql[i]) * std::log(phl[i + nhor] / phl[i]);
} }
#ifdef HAVE_OPENMP4 #ifdef HAVE_OPENMP4
#pragma omp parallel for simd #pragma omp parallel for simd
#endif #endif
for (i = 0; i < nhor; i++) geop[i] = geop[i + nhor] + PlanetRD * gt[i] * (1.0 + vtmp * gq[i]) * z2log2; for (int i = 0; i < nhor; i++) geop[i] = geop[i + nhor] + PlanetRD * gt[i] * (1.0 + vtmp * gq[i]) * z2log2;
} }
else /* No humidity */ else // No humidity
{ {
geopl = geop + nhor; double *geopl = geop + nhor;
phl = ph + nhor; const double *phl = ph + nhor;
for (j = nlev; j > 1; j--) for (int j = nlev; j > 1; j--)
#ifdef HAVE_OPENMP4 #ifdef HAVE_OPENMP4
#pragma omp simd #pragma omp simd
#endif #endif
for (i = nhor * (j - 1); i < nhor * j; i++) geop[i] = geopl[i] + PlanetRD * gt[i] * std::log(phl[i] / ph[i]); for (int i = nhor * (j - 1); i < nhor * j; i++) geop[i] = geopl[i] + PlanetRD * gt[i] * std::log(phl[i] / ph[i]);
#ifdef HAVE_OPENMP4 #ifdef HAVE_OPENMP4
#pragma omp simd #pragma omp simd
#endif #endif
for (i = 0; i < nhor; i++) geop[i] = geopl[i] + PlanetRD * gt[i] * z2log2; for (int i = 0; i < nhor; i++) geop[i] = geopl[i] + PlanetRD * gt[i] * z2log2;
} }
#ifdef HAVE_OPENMP4 #ifdef HAVE_OPENMP4
#pragma omp parallel for simd #pragma omp parallel for simd
#endif #endif
for (i = 0; i < nhor * (nlev + 1); i++) geop[i] *= zrg; for (int i = 0; i < nhor * (nlev + 1); i++) geop[i] *= zrg;
} }
constexpr double SCALESLP = 101325.0; constexpr double SCALESLP = 101325.0;
...@@ -1046,25 +1017,24 @@ constexpr double SCALESLP = 101325.0; ...@@ -1046,25 +1017,24 @@ constexpr double SCALESLP = 101325.0;
void void
LayerWater(double *ww, double *ll, double pmax, double pmin, int DimGP, int HalfLevels, double *vct) LayerWater(double *ww, double *ll, double pmax, double pmin, int DimGP, int HalfLevels, double *vct)
{ {
int i, k;
int MaxLev, MinLev;
double pph[MaxLevel]; double pph[MaxLevel];
int k;
for (k = 0; k < HalfLevels; k++) pph[k] = vct[k] + vct[k + HalfLevels] * SCALESLP; for (k = 0; k < HalfLevels; k++) pph[k] = vct[k] + vct[k + HalfLevels] * SCALESLP;
for (k = 0; k < HalfLevels; k++) for (k = 0; k < HalfLevels; k++)
if (pph[k] > pmax) break; if (pph[k] > pmax) break;
MaxLev = k - 1; const auto MaxLev = k - 1;
for (k = HalfLevels - 1; k >= 0; k--) for (k = HalfLevels - 1; k >= 0; k--)
if (pph[k] < pmin) break; if (pph[k] < pmin) break;
MinLev = k; const auto MinLev = k;
varrayFill(DimGP, ll, 0.0); varrayFill(DimGP, ll, 0.0);
for (k = MaxLev; k <= MinLev; k++) for (k = MaxLev; k <= MinLev; k++)
{ {
for (i = 0; i < DimGP; i++) ll[i] += ww[i + k * DimGP] * (pph[k + 1] - pph[k]); for (int i = 0; i < DimGP; i++) ll[i] += ww[i + k * DimGP] * (pph[k + 1] - pph[k]);
} }
for (i = 0; i < DimGP; i++) ll[i] /= PlanetGrav; for (int i = 0; i < DimGP; i++) ll[i] /= PlanetGrav;
} }
/* ================================================= */ /* ================================================= */
...@@ -1074,28 +1044,27 @@ LayerWater(double *ww, double *ll, double pmax, double pmin, int DimGP, int Half ...@@ -1074,28 +1044,27 @@ LayerWater(double *ww, double *ll, double pmax, double pmin, int DimGP, int Half
void void
LayerCloud(double *cc, double *ll, double pmax, double pmin, int DimGP, int HalfLevels, double *vct) LayerCloud(double *cc, double *ll, double pmax, double pmin, int DimGP, int HalfLevels, double *vct)
{ {
int i, k;
int MaxLev, MinLev;
double pph[MaxLevel]; double pph[MaxLevel];
double ZEPSEC = 1.0e-12; constexpr double ZEPSEC = 1.0e-12;
int k;
for (k = 0; k < HalfLevels; k++) pph[k] = vct[k] + vct[k + HalfLevels] * SCALESLP; for (k = 0; k < HalfLevels; k++) pph[k] = vct[k] + vct[k + HalfLevels] * SCALESLP;
for (k = 0; k < HalfLevels; k++) for (k = 0; k < HalfLevels; k++)
if (pph[k] > pmax) break; if (pph[k] > pmax) break;
MaxLev = k - 1; const auto MaxLev = k - 1;
for (k = HalfLevels - 1; k >= 0; k--) for (k = HalfLevels - 1; k >= 0; k--)
if (pph[k] < pmin) break; if (pph[k] < pmin) break;
MinLev = k; const auto MinLev = k;
for (i = 0; i < DimGP; i++) ll[i] = 1. - cc[i + MaxLev * DimGP]; for (int i = 0; i < DimGP; i++) ll[i] = 1. - cc[i + MaxLev * DimGP];
for (k = MaxLev + 1; k <= MinLev; k++) for (k = MaxLev + 1; k <= MinLev; k++)
{ {
for (i = 0; i < DimGP; i++) for (int i = 0; i < DimGP; i++)
ll[i] ll[i]
*= (1. - std::max(cc[i + (k - 1) * DimGP], cc[i + k * DimGP])) / (1. - std::min(cc[i + (k - 1) * DimGP], 1. - ZEPSEC)); *= (1. - std::max(cc[i + (k - 1) * DimGP], cc[i + k * DimGP])) / (1. - std::min(cc[i + (k - 1) * DimGP], 1. - ZEPSEC));
} }
for (i = 0; i < DimGP; i++) ll[i] = 1. - ll[i]; for (int i = 0; i < DimGP; i++) ll[i] = 1. - ll[i];
} }
// Grid Point Computations // Grid Point Computations
......
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