Commit 38e94139 authored by Fabian Wachsmann's avatar Fabian Wachsmann
Browse files

Solved merge conflicts

parents 3f7f0362 19fd4f6b
......@@ -3,4 +3,6 @@ Ralf Mueller, <ralf.mueller AT mpimet.mpg.de>
Oliver Heidmann, <oliver.heidmann AT mpimet.mpg.de>
Cedrick Ansorge, <cedrick.ansorge AT mpimet.mpg.de>
Luis Kornblueh, <luis.kornblueh AT mpimet.mpg.de>
Fabian Wachsmann <wachsmann@dkrz.de>
Modali Kameswarrao <kameswarrao.modali@mpimet.mpg.de>
Ralf Quast, <ralf.quast AT brockmann-consult.de>
Uwe Schulzweida, Max Planck Institute for Meteorologie, (2018),
Climate Data Operators (CDO) User Guide, Version 1.9.3,
https://code.mpimet.mpg.de/projects/cdo/embedded/cdo.pdf.
If you use the Climate Data Operators (CDO) to process data for an article in a scientific publication, please cite:
Schulzweida, Uwe. (2019, February 6). CDO User Guide (Version 1.9.6). Zenodo. http://doi.org/10.5281/zenodo.2558193
BibTeX:
@misc{schulzweida_uwe_2019_2558193,
author = {Schulzweida, Uwe},
title = {CDO User Guide},
month = feb,
year = 2019,
doi = {10.5281/zenodo.2558193},
url = {https://doi.org/10.5281/zenodo.2558193}
}
2019-01-24 Uwe Schulzweida
2019-05-23 Uwe Schulzweida
* Using CDI library version 1.9.7
* Version 1.9.7 release
2019-03-25 Uwe Schulzweida
* expr:zonSTAT: wrong result (bug fix)
* expr::vertmean: fix wrong warning message about layer bounds
2019-03-25 Uwe Schulzweida
* gridarea: use radius from grid description if available
2019-03-15 Uwe Schulzweida
* Fix compile error: EXIT_FAILURE not declared in cdoDebugOutput.h [Bug #8899]
2019-03-07 Uwe Schulzweida
* adipot: use code 20 as input
2019-03-05 Uwe Schulzweida
* inttime, intntime: handling of missing values is incorrect (bug fix)
2019-03-01 Uwe Schulzweida
* fldmean: added support for zonal mean data without longitude information
2019-02-26 Uwe Schulzweida
* varsavg, varsmean, varsstd, varsvar: wrong result if first record contains missing values (bug fix)
2019-02-20 Uwe Schulzweida
* uvRelativeToGrid: changed flag from grid to variable
2019-02-18 Uwe Schulzweida
* select: combination of some parameter (var, grid, zaxis) doesn't work (bug fix)
2019-02-09 Uwe Schulzweida
* griddes: print text attributes containing double quotes, in single quotes
2019-02-07 Uwe Schulzweida
* Using CDI library version 1.9.6
* Version 1.9.6 release
2019-02-01 Uwe Schulzweida
* Relative time axis (-r) returns wrong first timestep in operator chain for NetCDF
2019-01-17 Uwe Schulzweida
* setgridtype,regular: set nx=4*N+16 for octahedral reduced Gaussian grids (bug fix)
......@@ -14,7 +64,7 @@
2019-01-10 Uwe Schulzweida
* Wrong result with fldmean on zonal mean data (bug introduce in 1.9.6) [Bug #8834]
* Wrong result with fldmean on zonal mean data (bug introduce in 1.9.5) [Bug #8834]
2019-01-09 Uwe Schulzweida
......
......@@ -3,7 +3,7 @@ CDO NEWS
Improvement
Version 1.9.6 (24 January 2019):
Version 1.9.6 (7 February 2019):
New features:
* Added support for polar stereographic projection
......@@ -27,7 +27,8 @@ Version 1.9.6 (24 January 2019):
* Operator masklonlatbox: wrong result if lon1 > first lon || lon2 < last lon (bug introduce in 1.9.4) [Bug #8695]
* Operator maskindexbox: wrong result if idx1 > 1 || idx2 < nlon (bug introduce in 1.9.4) [Bug #8695]
* Absolute time axis (-a) returns wrong units in operator chain for NetCDF [Bug #8777]
* Wrong result with fldmean on zonal mean data (bug introduce in 1.9.6) [Bug #8834]
* Relative time axis (-r) returns wrong first timestep in operator chain for NetCDF
* Wrong result with fldmean on zonal mean data (bug introduce in 1.9.5) [Bug #8834]
Version 1.9.5 (9 August 2018):
......
......@@ -52,16 +52,23 @@ case "${HOSTNAME}" in
CC=icc CFLAGS="-g -Wall -O2 -qopt-report=5 -march=native -fp-model strict"
elif test "$COMP" = pgi ; then
${CONFPATH}configure --disable-openmp \
$CDOLIBS \
$CDOLIBS --without-eccodes \
CXX=pgc++ CXXFLAGS="-g -O1" \
CC=pgcc CFLAGS="-g -O1"
elif test "$COMP" = clang0 ; then
${CONFPATH}configure --prefix=$HOME/local \
--enable-maintainer-mode \
$CDOLIBS \
LDFLAGS="-Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
CXX=clang++ CXXFLAGS="-g -Wall -Wextra -O0" \
CC=clang CFLAGS="-g -Wall -Wextra -O0"
elif test "$COMP" = clang ; then
${CONFPATH}configure --prefix=$HOME/local \
--enable-maintainer-mode \
$CDOLIBS \
LDFLAGS="-Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
CXX=clang++ CXXFLAGS="-g -Wall -Wextra -O3" \
CC=clang CFLAGS="-g -Wall -Wextra -O3"
CXX=clang++ CXXFLAGS="-g -Wall -Wextra -O3 -march=native" \
CC=clang CFLAGS="-g -Wall -Wextra -O3 -march=native"
elif test "$COMP" = gnu4.9 ; then
${CONFPATH}configure --prefix=$HOME/local \
$CDOLIBS \
......@@ -77,13 +84,19 @@ case "${HOSTNAME}" in
elif test "$COMP" = gnu_pic ; then
${CONFPATH}configure \
$CDOLIBS \
LDFLAGS="/opt/local/lib/gcc6/libstdc++.6.dylib -Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
LDFLAGS="/opt/local/lib/gcc8/libstdc++.6.dylib -Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
CXX=g++ CXXFLAGS="-g -fPIC" \
CC=gcc CFLAGS="-g -fPIC"
elif test "$COMP" = gnu_stack ; then
${CONFPATH}configure \
$CDOLIBS \
LDFLAGS="/opt/local/lib/gcc8/libstdc++.6.dylib -Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
CXX=g++ CXXFLAGS="-g -O2 -fstack-protector-strong" \
CC=gcc CFLAGS="-g -O2 -fstack-protector-strong"
elif test "$COMP" = gribapi ; then
${CONFPATH}configure --disable-cgribex \
$CDOLIBS \
LDFLAGS="/opt/local/lib/gcc6/libstdc++.6.dylib -Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
LDFLAGS="/opt/local/lib/gcc8/libstdc++.6.dylib -Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
CXX=g++ CXXFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native" \
CC=gcc CFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native"
elif test "$COMP" = sanitize_address ; then
......@@ -104,7 +117,7 @@ case "${HOSTNAME}" in
${CONFPATH}configure --prefix=$HOME/local \
--enable-maintainer-mode \
$CDOLIBS \
LDFLAGS="/opt/local/lib/gcc6/libstdc++.6.dylib -Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
LDFLAGS="/opt/local/lib/gcc8/libstdc++.6.dylib -Wl,-rpath,$HOME/local/eccodes-2.6.0/lib" \
CXX=g++ CXXFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native" \
CC=gcc CFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native"
fi
......@@ -215,9 +228,9 @@ case "${HOSTNAME}" in
fi
;;
# jessie workstation x64
breeze*)
breeze1)
CDOLIBS="--with-eecodes=/sw/jessie-x64/eccodes/eccodes-2.4.1-gccsys \
--with-netcdf=/sw/jessie-x64/netcdf-4.3.3.1-gccsys \
--with-netcdf=/sw/jessie-x64/netcdf_c-4.4.1.1 \
--with-udunits2=/sw/jessie-x64/udunits-2.2.20-gccsys \
--with-proj=/sw/jessie-x64/proj4-4.9.3-gccsys LIBS=-lz"
if test "$COMP" = intel ; then
......@@ -230,10 +243,37 @@ case "${HOSTNAME}" in
$CDOLIBS \
CC=pgcc CXX=pgc++ CFLAGS="-g"
else
${CONFPATH}configure \
--with-fftw3 \
$CDOLIBS \
CXX=g++ CXXFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native" \
CC=gcc CFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native"
fi
;;
# stretch workstation x64
breeze2)
CDOLIBS="--with-netcdf=/sw/stretch-x64/netcdf/netcdf_c-4.6.1 \
--with-udunits2=/sw/stretch-x64/sys/udunits-2.2.24-gccsys"
if test "$COMP" = intel ; then
${CONFPATH}configure --prefix=$HOME/local --exec_prefix=$HOME/local/$COMP \
--with-fftw3 \
$CDOLIBS \
CC=icc CXX=icpc CFLAGS="-g -Wall -O2 -qopt-report=5 -march=native"
elif test "$COMP" = pgi ; then
${CONFPATH}configure --prefix=$HOME/local --exec_prefix=$HOME/local/$COMP \
$CDOLIBS \
CC=pgcc CXX=pgc++ CFLAGS="-g"
elif test "$COMP" = gnu_stack ; then
${CONFPATH}configure \
$CDOLIBS \
CXX=g++ CXXFLAGS="-g -O2 -fstack-protector-strong" \
CC=gcc CFLAGS="-g -O2 -fstack-protector-strong"
else
${CONFPATH}configure \
--with-fftw3 \
$CDOLIBS \
CFLAGS='-g -Wall -O3'
CXX=g++ CXXFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native" \
CC=gcc CFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native"
fi
;;
# mistral
......
......@@ -5,7 +5,7 @@
# libtool 2.4.2
AC_PREREQ([2.68])
AC_INIT([cdo], [1.9.6rc5], [http://mpimet.mpg.de/cdo])
AC_INIT([cdo], [1.9.7rc1], [http://mpimet.mpg.de/cdo])
AC_DEFINE_UNQUOTED(CDO, ["$PACKAGE_VERSION"], [CDO version])
......@@ -287,11 +287,11 @@ AC_PROG_AWK
AC_CONFIG_FILES([test/File.test test/Read_grib.test test/Read_netcdf.test test/Copy_netcdf.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Comp.test test/Compc.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Cat.test test/Gridarea.test test/Genweights.test test/Remap.test test/Remap2.test test/Remap3.test test/Remapgrid.test test/Remapeta.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/EOF.test test/Select.test test/Spectral.test test/Vertint.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Timstat.test test/Timstat2.test test/Timselstat.test test/Seasstat.test test/Timpctl.test test/Runstat.test test/Runpctl.test test/Multiyearstat.test test/Ydrunstat.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Gridboxstat.test test/Vertstat.test test/Fldstat.test test/Fldpctl.test test/Ensstat.test test/Enspctl.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/EOF.test test/Select.test test/Spectral.test test/Inttime.test test/Vertint.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Timstat.test test/Timstat2.test test/Yearmonstat.test test/Timselstat.test test/Seasstat.test test/Timpctl.test test/Runstat.test test/Runpctl.test test/Multiyearstat.test test/Ydrunstat.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Gridboxstat.test test/Vertstat.test test/Varsstat.test test/Fldstat.test test/Fldpctl.test test/Ensstat.test test/Enspctl.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Merstat.test test/Zonstat.test test/Mergetime.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Afterburner.test test/Detrend.test test/Arithc.test test/Arith.test test/Ymonarith.test test/Expr.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Afterburner.test test/Detrend.test test/Arithc.test test/Arith.test test/Monarith.test test/Ymonarith.test test/Expr.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Gradsdes.test test/Collgrid.test test/threads.test test/tsformat.test test/wildcard.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Setmiss.test test/Smooth.test test/MapReduce.test test/Ninfo.test],[chmod a+x "$ac_file"])
AC_CONFIG_FILES([test/Filter.test ],[chmod a+x "$ac_file"])
......
......@@ -198,7 +198,7 @@ keepaspectratio]{cdo_libdep.pdf}}%
\end{picture}
\begin{flushright}
\large \textbf{Climate Data Operator \\ Version 1.9.6 \\ January 2019}
\large \textbf{Climate Data Operator \\ Version 1.9.7 \\ May 2019}
\end{flushright}
\vfill
......
......@@ -14,7 +14,7 @@
\put(0,0.0){\line(1,0){3.95}}
\end{picture}
\begin{flushright}
{\small{Climate Data Operator \\ Version 1.9.6 \\ January 2019}}
{\small{Climate Data Operator \\ Version 1.9.7 \\ May 2019}}
\end{flushright}
\vspace*{0mm}
......
......@@ -321,6 +321,55 @@ attributes are available:
\item \textbf{polar\_stereographic}
\end{itemize}
It is recommend to set the attribute \textbf{proj4\_params} also for
the above projections to make sure all \textbf{proj4} parameter are set correctly.
\vspace{2mm}
%\begin{minipage}[t]{\textwidth}
Here is an example of a {\CDO} grid description using the attribute
\textbf{proj4\_params} to define the proj4 parameter of a polar
stereographic projection:
\begin{lstlisting}[frame=single, backgroundcolor=\color{pcolor1}, basicstyle=\footnotesize]
gridtype = projection
xsize = 11
ysize = 11
xunits = "meter"
yunits = "meter"
xfirst = -638000
xinc = 150
yfirst = -3349350
yinc = 150
grid_mapping = crs
grid_mapping_name = polar_stereographic
proj4_params = "+proj=stere +lon_0=-45 +lat_ts=70 +lat_0=90 +x_0=0 +y_0=0"
\end{lstlisting}
%\end{minipage}
\vspace{2mm}
%\begin{minipage}[t]{\textwidth}
The result is the same as using the CF conform Grid Mapping Attributes:
\begin{lstlisting}[frame=single, backgroundcolor=\color{pcolor1}, basicstyle=\footnotesize]
gridtype = projection
xsize = 11
ysize = 11
xunits = "meter"
yunits = "meter"
xfirst = -638000
xinc = 150
yfirst = -3349350
yinc = 150
grid_mapping = crs
grid_mapping_name = polar_stereographic
straight_vertical_longitude_from_pole = -45.
standard_parallel = 70.
latitude_of_projection_origin = 90.
false_easting = 0.
false_northing = 0.
\end{lstlisting}
%\end{minipage}
\vspace{2mm}
%\begin{minipage}[t]{\textwidth}
......
......@@ -314,6 +314,7 @@ while (<MOFILE>) {
print HELPFILE "\n";
print HELPFILE "static const std::vector<std::string> ${mname}Help = {\n";
# print HELPFILE "static const std::string ${mname}Help[] = {\n";
@hkeys = split(" ", $moperators);
# print "$#hkeys @hkeys \n";
......@@ -1039,6 +1040,7 @@ while (<MOFILE>) {
# print HELPFILE " \"\",\n";
# print HELPFILE " \"\@End_${operator}\",\n";
# print HELPFILE " \"ENDHELP\"\n";
print HELPFILE "};\n";
print TCFILE "\\end{tabular*}\n";
......
......@@ -5,18 +5,16 @@
[ ! -f cdo_eca.pdf ] && ./makepdfeca
[ ! -f cdo_magics.pdf ] && ./makepdfmagics
if [ ! -d ../html ];then mkdir ../html ;fi
htlatex cdo.tex "html" "" -d../html/
cd ../html
#htlatex cdo.tex "html" "" -d../html/
htlatex cdo.tex 'xhtml,charset="us-ascii"'
#check for html cleanup tool
which tidy
if [[ 0 = $? ]]; then
tidy -c -i -wrap 200 cdo.html > index.html
tidy -latin0 -c -i -wrap 200 cdo.html > index.html
else
cp cdo.html index.html
fi
cp ../tex/*pdf .
mkdir grids
mv cell.png grids
mv curv.png grids
zip cdo-doc.zip *
zip cdo-doc.zip *.pdf *png *html grids/* figures/* logo/*
......@@ -25,6 +25,8 @@ The "afterburner" is the standard post processor for @cite{ECHAM} data which pro
This operator reads selection parameters as namelist from stdin.
Use the UNIX redirection "<namelistfile" to read the namelist from file.
The input files can't be combined with other CDO operators because of an optimized reader for this operator.
@EndDescription
@EndOperator
......
libcdi @ 4f669581
Subproject commit e8d41cde3c56bb0fb7c9ac0e25d20f65bd69e9bb
Subproject commit 4f6695814c2be9f5c78d81e8c945a6b51fb0edb3
......@@ -27,7 +27,7 @@
#include "cdo_options.h"
#include "cdo_int.h"
#include "param_conversion.h"
#include "pstream_int.h"
#include "util_string.h"
/*
......@@ -44,25 +44,25 @@
*/
/* compute insitu temperature from potential temperature */
static double
adisit_1(double tpot, double sal, double p)
static inline double
adisit_1(const double tpot, const double sal, const double p)
{
double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
double qc = p * (a_a1 + p * (a_c1 - a_e1 * p));
double qv = p * (a_b1 - a_d * p);
double dc = 1. + p * (-a_a2 + p * (a_c2 - a_e2 * p));
double dv = a_b2 * p;
double qnq = -p * (-a_a3 + p * a_c3);
double qn3 = -p * a_a4;
double tpo = tpot;
double qvs = qv * (sal - 35.) + qc;
double dvs = dv * (sal - 35.) + dc;
constexpr double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
const double qc = p * (a_a1 + p * (a_c1 - a_e1 * p));
const double qv = p * (a_b1 - a_d * p);
const double dc = 1. + p * (-a_a2 + p * (a_c2 - a_e2 * p));
const double dv = a_b2 * p;
const double qnq = -p * (-a_a3 + p * a_c3);
const double qn3 = -p * a_a4;
const double tpo = tpot;
const double qvs = qv * (sal - 35.) + qc;
const double dvs = dv * (sal - 35.) + dc;
double t = (tpo + qvs) / dvs;
double fne = -qvs + t * (dvs + t * (qnq + t * qn3)) - tpo;
double fst = dvs + t * (2. * qnq + 3. * qn3 * t);
const double fne = -qvs + t * (dvs + t * (qnq + t * qn3)) - tpo;
const double fst = dvs + t * (2. * qnq + 3. * qn3 * t);
t = t - fne / fst;
return t;
......@@ -70,76 +70,80 @@ adisit_1(double tpot, double sal, double p)
/* compute potential temperature from insitu temperature */
/* Ref: Gill, p. 602, Section A3.5:Potential Temperature */
static double
adipot(double t, double s, double p)
static inline double
adipot(const double t, const double s, const double p)
{
double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
constexpr double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
double s_rel = s - 35.0;
const double s_rel = s - 35.0;
double aa = (a_a1 + t * (a_a2 - t * (a_a3 - a_a4 * t)));
double bb = s_rel * (a_b1 - a_b2 * t);
double cc = (a_c1 + t * (-a_c2 + a_c3 * t));
double cc1 = a_d * s_rel;
double dd = (-a_e1 + a_e2 * t);
const double aa = (a_a1 + t * (a_a2 - t * (a_a3 - a_a4 * t)));
const double bb = s_rel * (a_b1 - a_b2 * t);
const double cc = (a_c1 + t * (-a_c2 + a_c3 * t));
const double cc1 = a_d * s_rel;
const double dd = (-a_e1 + a_e2 * t);
double tpot = t - p * (aa + bb + p * (cc - cc1 + p * dd));
const double tpot = t - p * (aa + bb + p * (cc - cc1 + p * dd));
return tpot;
}
static void
calc_adisit(long gridsize, long nlevel, double *pressure, Field tho, Field sao, Field tis)
calc_adisit(size_t gridsize, size_t nlevel, const std::vector<double> &pressure, const FieldVector &tho,
const FieldVector &sao, FieldVector &tis)
{
// pressure units: hPa
// tho units: Celsius
// sao units: psu
for (long levelID = 0; levelID < nlevel; ++levelID)
for (size_t levelID = 0; levelID < nlevel; ++levelID)
{
long offset = gridsize * levelID;
double *thoptr = tho.ptr + offset;
double *saoptr = sao.ptr + offset;
double *tisptr = tis.ptr + offset;
for (long i = 0; i < gridsize; ++i)
const std::vector<double> &thovec = tho[levelID].vec;
const std::vector<double> &saovec = sao[levelID].vec;
std::vector<double> &tisvec = tis[levelID].vec;
const double tho_missval = tho[levelID].missval;
const double sao_missval = sao[levelID].missval;
const double tis_missval = tis[levelID].missval;
for (size_t i = 0; i < gridsize; ++i)
{
if (DBL_IS_EQUAL(thoptr[i], tho.missval) || DBL_IS_EQUAL(saoptr[i], sao.missval))
if (DBL_IS_EQUAL(thovec[i], tho_missval) || DBL_IS_EQUAL(saovec[i], sao_missval))
{
tisptr[i] = tis.missval;
tisvec[i] = tis_missval;
}
else
{
tisptr[i] = adisit_1(thoptr[i], saoptr[i], pressure[levelID]);
tisvec[i] = adisit_1(thovec[i], saovec[i], pressure[levelID]);
}
}
}
}
static void
calc_adipot(long gridsize, long nlevel, double *pressure, Field t, Field s, Field tpot)
calc_adipot(size_t gridsize, size_t nlevel, const std::vector<double> &pressure, const std::vector<Field> &t,
const std::vector<Field> &s, std::vector<Field> &tpot)
{
// pressure units: hPa
// t units: Celsius
// s units: psu
for (long levelID = 0; levelID < nlevel; ++levelID)
for (size_t levelID = 0; levelID < nlevel; ++levelID)
{
long offset = gridsize * levelID;
double *tptr = t.ptr + offset;
double *sptr = s.ptr + offset;
double *tpotptr = tpot.ptr + offset;
for (long i = 0; i < gridsize; ++i)
const std::vector<double> &tvec = t[levelID].vec;
const std::vector<double> &svec = s[levelID].vec;
std::vector<double> &tpotvec = tpot[levelID].vec;
const double t_missval = t[levelID].missval;
const double s_missval = s[levelID].missval;
const double tpot_missval = tpot[levelID].missval;
for (size_t i = 0; i < gridsize; ++i)
{
if (DBL_IS_EQUAL(tptr[i], t.missval) || DBL_IS_EQUAL(sptr[i], s.missval))
if (DBL_IS_EQUAL(tvec[i], t_missval) || DBL_IS_EQUAL(svec[i], s_missval))
{
tpotptr[i] = tpot.missval;
tpotvec[i] = tpot_missval;
}
else
{
tpotptr[i] = adipot(tptr[i], sptr[i], pressure[levelID]);
tpotvec[i] = adipot(tvec[i], svec[i], pressure[levelID]);
}
}
}
......@@ -150,45 +154,36 @@ Adisit(void *process)
{
int nrecs;
int varID, levelID;
size_t offset;
int i;
size_t nmiss;
int thoID = -1, saoID = -1;
char varname[CDI_MAX_NAME], stdname[CDI_MAX_NAME];
double pin = -1;
double *single;
cdoInitialize(process);
int ADISIT = cdoOperatorAdd("adisit", 1, 1, "");
int ADIPOT = cdoOperatorAdd("adipot", 1, 1, "");
UNUSED(ADIPOT);
const int ADISIT = cdoOperatorAdd("adisit", 1, 1, "");
const int ADIPOT = cdoOperatorAdd("adipot", 1, 1, "");
int operatorID = cdoOperatorID();
const int operatorID = cdoOperatorID();
if (operatorArgc() == 1) pin = parameter2double(operatorArgv()[0]);
const double pin = (operatorArgc() == 1) ? parameter2double(operatorArgv()[0]) : -1;
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int vlistID1 = cdoStreamInqVlist(streamID1);
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
const int vlistID1 = cdoStreamInqVlist(streamID1);
int nvars = vlistNvars(vlistID1);
const int nvars = vlistNvars(vlistID1);
for (varID = 0; varID < nvars; varID++)
{
int code = vlistInqVarCode(vlistID1, varID);
if (code <= 0)
{
vlistInqVarName(vlistID1, varID, varname);
vlistInqVarStdname(vlistID1, varID, stdname);
strtolower(varname);
if (cstrIsEqual(varname, "s"))
if (cstrIsEqual(varname, "s") || cstrIsEqual(stdname, "sea_water_salinity"))
code = 5;
else if (cstrIsEqual(varname, "t"))
code = 2;
else if (cstrIsEqual(stdname, "sea_water_salinity"))
code = 5;
if (operatorID == ADISIT)
{
......@@ -200,25 +195,25 @@ Adisit(void *process)
}
}