Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
mpim-sw
cdo
Commits
941aa212
Commit
941aa212
authored
Oct 24, 2020
by
Uwe Schulzweida
Browse files
Replace malloced arrays by Varray.
parent
aae43f79
Changes
6
Hide whitespace changes
Inline
Side-by-side
src/Eof3d.cc
View file @
941aa212
...
...
@@ -68,8 +68,7 @@ EOF3d(void *process)
int
calendar
=
CALENDAR_STANDARD
;
double
sum_w
;
double
**
cov
=
nullptr
;
/* TODO: covariance matrix / eigenvectors after solving */
double
*
eigv
;
double
**
cov
=
nullptr
;
// TODO: covariance matrix / eigenvectors after solving
if
(
Options
::
Timer
)
{
...
...
@@ -94,7 +93,7 @@ EOF3d(void *process)
auto
eigen_mode
=
get_eigenmode
();
auto
weight_mode
=
get_weightmode
();
/
*
eigenvalues
*/
/
/
eigenvalues
if
(
operfunc
==
EOF3D_SPATIAL
)
cdoAbort
(
"Operator not Implemented - use eof3d or eof3dtime instead"
);
...
...
@@ -104,7 +103,7 @@ EOF3d(void *process)
VarList
varList1
;
varListInit
(
varList1
,
vlistID1
);
/
*
COUNT NUMBER OF TIMESTEPS if EOF3D_ or EOF3D_TIME
*/
/
/
COUNT NUMBER OF TIMESTEPS if EOF3D_ or EOF3D_TIME
auto
nts
=
vlistNtsteps
(
vlistID1
);
if
(
nts
==
-
1
)
{
...
...
@@ -123,8 +122,7 @@ EOF3d(void *process)
auto
taxisID1
=
vlistInqTaxis
(
vlistID1
);
/* reset the requested number of eigen-function to the maximum if neccessary
*/
// reset the requested number of eigen-function to the maximum if neccessary
if
(
n_eig
>
nts
)
{
cdoWarning
(
"Solving in time-space:"
);
...
...
@@ -143,12 +141,12 @@ EOF3d(void *process)
auto
gridID1
=
varList1
[
0
].
gridID
;
auto
gridsizemax
=
vlistGridsizeMax
(
vlistID1
);
/
*
allocation of temporary fields and output structures
*/
double
*
in
=
(
double
*
)
Malloc
(
gridsizemax
*
sizeof
(
double
)
);
int
**
datacounts
=
(
int
**
)
Malloc
(
nvars
*
sizeof
(
int
*
)
);
double
***
datafields
=
(
double
***
)
Malloc
(
nvars
*
sizeof
(
double
**
)
);
double
***
eigenvectors
=
(
double
***
)
Malloc
(
nvars
*
sizeof
(
double
**
)
);
double
***
eigenvalues
=
(
double
***
)
Malloc
(
nvars
*
sizeof
(
double
**
)
);
/
/
allocation of temporary fields and output structures
Varray
<
double
>
in
(
gridsizemax
);
Varray2D
<
int
>
datacounts
(
nvars
);
Varray3D
<
double
>
datafields
(
nvars
);
Varray3D
<
double
>
eigenvectors
(
nvars
);
Varray3D
<
double
>
eigenvalues
(
nvars
);
size_t
maxlevs
=
0
;
for
(
varID
=
0
;
varID
<
nvars
;
++
varID
)
...
...
@@ -160,30 +158,19 @@ EOF3d(void *process)
if
(
nlevs
>
maxlevs
)
maxlevs
=
nlevs
;
data
count
s
[
varID
]
=
(
int
*
)
Malloc
(
nlevs
*
sizeof
(
i
nt
)
);
datafields
[
varID
]
=
(
double
**
)
Malloc
(
nts
*
sizeof
(
double
*
)
);
data
field
s
[
varID
]
.
resize
(
nt
s
);
for
(
int
tsID
=
0
;
tsID
<
nts
;
tsID
++
)
datafields
[
varID
][
tsID
].
resize
(
temp_size
,
0.0
);
for
(
int
tsID
=
0
;
tsID
<
nts
;
tsID
++
)
{
datafields
[
varID
][
tsID
]
=
(
double
*
)
Malloc
(
temp_size
*
sizeof
(
double
));
for
(
size_t
i
=
0
;
i
<
temp_size
;
++
i
)
datafields
[
varID
][
tsID
][
i
]
=
0
;
}
datacounts
[
varID
]
=
(
int
*
)
Malloc
(
temp_size
*
sizeof
(
int
));
for
(
size_t
i
=
0
;
i
<
temp_size
;
i
++
)
datacounts
[
varID
][
i
]
=
0
;
datacounts
[
varID
].
resize
(
temp_size
,
0
);
eigenvectors
[
varID
]
=
(
double
**
)
Malloc
(
n_eig
*
sizeof
(
double
*
)
);
eigenvalues
[
varID
]
=
(
double
**
)
Malloc
(
nts
*
sizeof
(
double
*
)
);
eigenvectors
[
varID
]
.
resize
(
n_eig
);
eigenvalues
[
varID
]
.
resize
(
nts
);
for
(
int
i
=
0
;
i
<
n
;
i
++
)
{
if
(
i
<
n_eig
)
{
eigenvectors
[
varID
][
i
]
=
(
double
*
)
Malloc
(
temp_size
*
sizeof
(
double
));
for
(
size_t
i2
=
0
;
i2
<
temp_size
;
++
i2
)
eigenvectors
[
varID
][
i
][
i2
]
=
missval
;
}
if
(
i
<
n_eig
)
eigenvectors
[
varID
][
i
].
resize
(
temp_size
,
missval
);
eigenvalues
[
varID
][
i
]
=
(
double
*
)
Malloc
(
1
*
sizeof
(
double
));
eigenvalues
[
varID
][
i
][
0
]
=
missval
;
eigenvalues
[
varID
][
i
].
resize
(
1
,
missval
);
}
}
...
...
@@ -191,12 +178,11 @@ EOF3d(void *process)
cdoPrint
(
"Allocated eigenvalue/eigenvector with nts=%i, n=%i, gridsize=%zu for processing in %s"
,
nts
,
n
,
gridsizemax
,
"time_space"
);
double
*
weight
=
(
double
*
)
Malloc
(
maxlevs
*
gridsizemax
*
sizeof
(
double
));
for
(
size_t
i
=
0
;
i
<
maxlevs
*
gridsizemax
;
++
i
)
weight
[
i
]
=
1.
;
Varray
<
double
>
weight
(
maxlevs
*
gridsizemax
,
1.0
);
if
(
weight_mode
==
WEIGHT_ON
)
{
int
wstatus
=
gridWeights
(
gridID1
,
weight
);
int
wstatus
=
gridWeights
(
gridID1
,
weight
.
data
()
);
if
(
wstatus
!=
0
)
{
weight_mode
=
WEIGHT_OFF
;
...
...
@@ -211,7 +197,7 @@ EOF3d(void *process)
int
tsID
=
0
;
/
*
read the data and create covariance matrices for each var & level
*/
/
/
read the data and create covariance matrices for each var & level
while
(
true
)
{
nrecs
=
cdoStreamInqTimestep
(
streamID1
,
tsID
);
...
...
@@ -224,7 +210,7 @@ EOF3d(void *process)
const
auto
gridsize
=
varList1
[
varID
].
gridsize
;
const
auto
missval
=
varList1
[
varID
].
missval
;
cdoReadRecord
(
streamID1
,
in
,
&
nmiss
);
cdoReadRecord
(
streamID1
,
in
.
data
()
,
&
nmiss
);
const
auto
offset
=
gridsize
*
levelID
;
for
(
size_t
i
=
0
;
i
<
gridsize
;
++
i
)
...
...
@@ -252,7 +238,7 @@ EOF3d(void *process)
if
(
Options
::
cdoVerbose
)
cdoPrint
(
"Read data for %i variables"
,
nvars
);
size_t
*
pack
=
(
size_t
*
)
Malloc
(
temp_size
*
sizeof
(
size_t
)
);
// TODO
Varray
<
size_t
>
pack
(
temp_size
);
// TODO
for
(
varID
=
0
;
varID
<
nvars
;
varID
++
)
{
...
...
@@ -295,7 +281,7 @@ EOF3d(void *process)
cov
=
(
double
**
)
Malloc
(
nts
*
sizeof
(
double
*
));
for
(
int
j1
=
0
;
j1
<
nts
;
j1
++
)
cov
[
j1
]
=
(
double
*
)
Malloc
(
nts
*
sizeof
(
double
));
eigv
=
(
double
*
)
Malloc
(
n
*
sizeof
(
double
)
);
Varray
<
double
>
eigv
(
n
);
if
(
Options
::
cdoVerbose
)
{
...
...
@@ -308,11 +294,11 @@ EOF3d(void *process)
#endif
for
(
int
j1
=
0
;
j1
<
nts
;
j1
++
)
{
double
*
df1p
=
datafields
[
varID
][
j1
];
double
*
df1p
=
datafields
[
varID
][
j1
]
.
data
()
;
for
(
int
j2
=
j1
;
j2
<
nts
;
j2
++
)
{
double
*
df2p
=
datafields
[
varID
][
j2
];
double
sum
=
0
;
double
*
df2p
=
datafields
[
varID
][
j2
]
.
data
()
;
double
sum
=
0.
0
;
for
(
size_t
i
=
0
;
i
<
npack
;
i
++
)
sum
+=
weight
[
pack
[
i
]
%
gridsizemax
]
*
df1p
[
pack
[
i
]]
*
df2p
[
pack
[
i
]];
cov
[
j2
][
j1
]
=
cov
[
j1
][
j2
]
=
sum
/
sum_w
/
nts
;
}
...
...
@@ -328,10 +314,10 @@ EOF3d(void *process)
if
(
Options
::
cdoVerbose
)
cdoPrint
(
"Processed correlation matrix for var %2i | npack: %zu"
,
varID
,
n
);
if
(
eigen_mode
==
JACOBI
)
cdo
::
parallel_eigen_solution_of_symmetric_matrix
(
cov
,
eigv
,
n
,
__func__
);
cdo
::
parallel_eigen_solution_of_symmetric_matrix
(
cov
,
eigv
.
data
()
,
n
,
__func__
);
else
cdo
::
eigen_solution_of_symmetric_matrix
(
cov
,
eigv
,
n
,
__func__
);
/
*
NOW: cov contains the eigenvectors, eigv the eigenvalues
*/
cdo
::
eigen_solution_of_symmetric_matrix
(
cov
,
eigv
.
data
()
,
n
,
__func__
);
/
/
NOW: cov contains the eigenvectors, eigv the eigenvalues
if
(
Options
::
cdoVerbose
)
cdoPrint
(
"Processed SVD decomposition for var %d from %d x %d matrix"
,
varID
,
n
,
n
);
...
...
@@ -341,21 +327,21 @@ EOF3d(void *process)
for
(
int
eofID
=
0
;
eofID
<
n_eig
;
eofID
++
)
{
double
*
eigenvec
=
eigenvectors
[
varID
][
eofID
];
double
*
eigenvec
=
eigenvectors
[
varID
][
eofID
]
.
data
()
;
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(varID, nts, eofID, npack, pack, cov, datafields, eigenvec)
#endif
for
(
size_t
i
=
0
;
i
<
npack
;
i
++
)
{
double
sum
=
0
;
double
sum
=
0.
0
;
for
(
int
j
=
0
;
j
<
nts
;
j
++
)
sum
+=
datafields
[
varID
][
j
][
pack
[
i
]]
*
cov
[
eofID
][
j
];
eigenvec
[
pack
[
i
]]
=
sum
;
}
// NORMALIZING
double
sum
=
0
;
double
sum
=
0.
0
;
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(eigenvec, weight, pack, npack, gridsizemax) reduction(+ : sum)
...
...
@@ -377,16 +363,15 @@ EOF3d(void *process)
#endif
for
(
size_t
i
=
0
;
i
<
npack
;
i
++
)
eigenvec
[
pack
[
i
]]
=
missval
;
}
}
/
*
for ( eofID = 0; eofID < n_eig; eofID++ )
*/
}
/
/
for ( eofID = 0; eofID < n_eig; eofID++ )
if
(
eigv
)
Free
(
eigv
);
for
(
int
i
=
0
;
i
<
n
;
i
++
)
if
(
cov
[
i
])
Free
(
cov
[
i
]);
}
/
*
for ( varID = 0; varID < nvars; varID++ )
*/
}
/
/
for ( varID = 0; varID < nvars; varID++ )
/
*
write files with eigenvalues (ID3) and eigenvectors (ID2)
*/
/
/
write files with eigenvalues (ID3) and eigenvectors (ID2)
/
*
eigenvalues
*/
/
/
eigenvalues
const
auto
streamID2
=
cdoOpenWrite
(
1
);
const
auto
vlistID2
=
vlistDuplicate
(
vlistID1
);
...
...
@@ -463,34 +448,10 @@ EOF3d(void *process)
nmiss
=
(
DBL_IS_EQUAL
(
eigenvalues
[
varID
][
tsID
][
0
],
missval
))
?
1
:
0
;
cdoDefRecord
(
streamID2
,
varID
,
0
);
cdoWriteRecord
(
streamID2
,
eigenvalues
[
varID
][
tsID
],
nmiss
);
cdoWriteRecord
(
streamID2
,
eigenvalues
[
varID
][
tsID
]
.
data
()
,
nmiss
);
}
// for ( varID = 0; ... )
}
// for ( tsID = 0; ... )
for
(
varID
=
0
;
varID
<
nvars
;
varID
++
)
{
for
(
int
i
=
0
;
i
<
nts
;
i
++
)
{
Free
(
datafields
[
varID
][
i
]);
if
(
i
<
n_eig
)
Free
(
eigenvectors
[
varID
][
i
]);
Free
(
eigenvalues
[
varID
][
i
]);
}
Free
(
datafields
[
varID
]);
Free
(
datacounts
[
varID
]);
Free
(
eigenvectors
[
varID
]);
Free
(
eigenvalues
[
varID
]);
}
Free
(
datafields
);
Free
(
datacounts
);
Free
(
eigenvectors
);
Free
(
eigenvalues
);
Free
(
in
);
Free
(
pack
);
Free
(
weight
);
cdoStreamClose
(
streamID3
);
cdoStreamClose
(
streamID2
);
cdoStreamClose
(
streamID1
);
...
...
src/Gridsearch.cc
View file @
941aa212
...
...
@@ -43,12 +43,12 @@ static grid_type *
grid_new
(
int
gridID
,
const
char
*
txt
)
{
bool
lgrid_destroy
=
false
;
int
gridtype
=
gridInqType
(
gridID
);
const
auto
gridtype
=
gridInqType
(
gridID
);
if
(
gridtype
==
GRID_GME
)
{
lgrid_destroy
=
true
;
int
gridID_gme
=
gridToUnstructured
(
gridID
,
1
);
const
auto
gridID_gme
=
gridToUnstructured
(
gridID
,
1
);
gridCompress
(
gridID_gme
);
gridID
=
gridID_gme
;
}
...
...
@@ -93,10 +93,10 @@ void
boundbox_from_corners1r
(
long
ic
,
long
nc
,
const
double
*
restrict
corner_lon
,
const
double
*
restrict
corner_lat
,
float
*
restrict
bound_box
)
{
l
on
g
inc
=
ic
*
nc
;
c
on
st
auto
inc
=
ic
*
nc
;
float
clat
=
corner_lat
[
inc
];
float
clon
=
corner_lon
[
inc
];
auto
clat
=
corner_lat
[
inc
];
auto
clon
=
corner_lon
[
inc
];
bound_box
[
0
]
=
clat
;
bound_box
[
1
]
=
clat
;
...
...
@@ -127,20 +127,17 @@ void
boundbox_from_corners
(
long
size
,
long
nc
,
const
double
*
restrict
corner_lon
,
const
double
*
restrict
corner_lat
,
float
*
restrict
bound_box
)
{
long
i4
,
inc
,
j
;
double
clon
,
clat
;
for
(
long
i
=
0
;
i
<
size
;
++
i
)
{
i4
=
i
<<
2
;
// *4
inc
=
i
*
nc
;
clat
=
corner_lat
[
inc
];
clon
=
corner_lon
[
inc
];
const
auto
i4
=
i
<<
2
;
// *4
const
auto
inc
=
i
*
nc
;
auto
clat
=
corner_lat
[
inc
];
auto
clon
=
corner_lon
[
inc
];
bound_box
[
i4
]
=
clat
;
bound_box
[
i4
+
1
]
=
clat
;
bound_box
[
i4
+
2
]
=
clon
;
bound_box
[
i4
+
3
]
=
clon
;
for
(
j
=
1
;
j
<
nc
;
++
j
)
for
(
long
j
=
1
;
j
<
nc
;
++
j
)
{
clat
=
corner_lat
[
inc
+
j
];
clon
=
corner_lon
[
inc
+
j
];
...
...
@@ -182,22 +179,20 @@ search_cells(cellsearch_type *cellsearch, long tgt_cell_add, long *srch_add)
grid_type
*
src_grid
=
cellsearch
->
src_grid
;
grid_type
*
tgt_grid
=
cellsearch
->
tgt_grid
;
float
*
src_cell_bound_box
=
cellsearch
->
src_cell_bound_box
;
float
tgt_cell_bound_box
[
4
];
float
tgt_cell_bound_box
[
4
];
boundbox_from_corners1r
(
tgt_cell_add
,
tgt_grid
->
num_cell_corners
,
tgt_grid
->
cell_corner_lon
,
tgt_grid
->
cell_corner_lat
,
tgt_cell_bound_box
);
long
src_cell_addm4
;
float
bound_box_lat1
=
tgt_cell_bound_box
[
0
];
float
bound_box_lat2
=
tgt_cell_bound_box
[
1
];
float
bound_box_lon1
=
tgt_cell_bound_box
[
2
];
float
bound_box_lon2
=
tgt_cell_bound_box
[
3
];
const
auto
bound_box_lat1
=
tgt_cell_bound_box
[
0
];
const
auto
bound_box_lat2
=
tgt_cell_bound_box
[
1
];
const
auto
bound_box_lon1
=
tgt_cell_bound_box
[
2
];
const
auto
bound_box_lon2
=
tgt_cell_bound_box
[
3
];
long
num_srch_cells
=
0
;
for
(
long
src_cell_add
=
0
;
src_cell_add
<
src_grid
->
size
;
++
src_cell_add
)
{
src_cell_addm4
=
src_cell_add
<<
2
;
const
auto
src_cell_addm4
=
src_cell_add
<<
2
;
if
((
src_cell_bound_box
[
src_cell_addm4
+
2
]
<=
bound_box_lon2
)
&&
(
src_cell_bound_box
[
src_cell_addm4
+
3
]
>=
bound_box_lon1
))
{
if
((
src_cell_bound_box
[
src_cell_addm4
]
<=
bound_box_lat2
)
&&
(
src_cell_bound_box
[
src_cell_addm4
+
1
]
>=
bound_box_lat1
))
...
...
src/Importcmsaf.cc
View file @
941aa212
...
...
@@ -29,6 +29,7 @@
#include
"cdo_options.h"
#include
"dmemory.h"
#include
"array.h"
#include
"process_int.h"
#include
"cdo_cdi_wrapper.h"
#include
"cdo_default_values.h"
...
...
@@ -169,10 +170,8 @@ defLonLatGrid(int nx, int ny, double c0, double lts, double re)
return
-
1
;
}
Varray
<
double
>
xvals
(
nx
);
Varray
<
double
>
yvals
(
ny
);
Varray
<
double
>
xbounds
(
nx
*
2
);
Varray
<
double
>
ybounds
(
ny
*
2
);
Varray
<
double
>
xvals
(
nx
),
yvals
(
ny
);
Varray
<
double
>
xbounds
(
nx
*
2
),
ybounds
(
ny
*
2
);
for
(
i
=
0
;
i
<
nx
;
++
i
)
{
...
...
@@ -218,8 +217,7 @@ defSinusoidalGrid(int nx, int ny, double xmin, double ymax, double dx, double dy
(
void
)
(
p2
);
// CDO_UNUSED
(
void
)
(
p3
);
// CDO_UNUSED
(
void
)
(
p4
);
// CDO_UNUSED
Varray
<
double
>
xvals
(
nx
);
Varray
<
double
>
yvals
(
ny
);
Varray
<
double
>
xvals
(
nx
),
yvals
(
ny
);
for
(
int
i
=
0
;
i
<
nx
;
++
i
)
xvals
[
i
]
=
xmin
+
i
*
dx
+
dx
/
2
;
...
...
@@ -240,8 +238,7 @@ defSinusoidalGrid(int nx, int ny, double xmin, double ymax, double dx, double dy
static
int
defLaeaGrid
(
int
nx
,
int
ny
,
double
xmin
,
double
ymax
,
double
dx
,
double
dy
,
double
a
,
double
lon0
,
double
lat0
)
{
Varray
<
double
>
xvals
(
nx
);
Varray
<
double
>
yvals
(
ny
);
Varray
<
double
>
xvals
(
nx
),
yvals
(
ny
);
for
(
int
i
=
0
;
i
<
nx
;
++
i
)
xvals
[
i
]
=
xmin
+
i
*
dx
+
dx
/
2
;
for
(
int
i
=
0
;
i
<
ny
;
++
i
)
yvals
[
i
]
=
ymax
-
i
*
dy
-
dy
/
2
;
...
...
@@ -1004,7 +1001,7 @@ read_dataset(hid_t loc_id, const char *name, void *opdata)
}
else
{
std
::
vector
<
int
>
iarray
(
gridsize
*
nt
);
Varray
<
int
>
iarray
(
gridsize
*
nt
);
status
=
H5Dread
(
dset_id
,
H5T_NATIVE_INT
,
H5S_ALL
,
H5S_ALL
,
H5P_DEFAULT
,
iarray
.
data
());
if
(
status
<
0
)
cdoAbort
(
"Reading of NATIVE_INT variable %s failed!"
,
varname
);
for
(
size_t
i
=
0
;
i
<
gridsize
*
nt
;
++
i
)
array
[
i
]
=
iarray
[
i
];
...
...
@@ -1080,8 +1077,8 @@ read_dataset(hid_t loc_id, const char *name, void *opdata)
for
(
size_t
i
=
0
;
i
<
gridsize
*
nt
;
i
++
)
if
(
mask
[
i
]
==
false
)
{
if
(
array
[
i
]
<
minval
)
minval
=
array
[
i
];
if
(
array
[
i
]
>
maxval
)
maxval
=
array
[
i
];
minval
=
std
::
min
(
minval
,
array
[
i
]
)
;
maxval
=
std
::
max
(
maxval
,
array
[
i
]
)
;
}
if
(
Options
::
cdoVerbose
)
...
...
@@ -1307,13 +1304,9 @@ Importcmsaf(void *process)
size_t
nmiss
;
int
ivar
;
int
varID
,
levelID
,
tsID
;
double
*
array
;
double
missval
;
herr_t
status
;
// Generic return value
datasets_t
dsets
;
int64_t
vdate
;
int
vtime
;
int
*
vtimes
=
nullptr
;
#endif
cdoInitialize
(
process
);
...
...
@@ -1348,9 +1341,10 @@ Importcmsaf(void *process)
break
;
}
Varray
<
int
>
vtimes
;
if
(
nt
>
1
)
{
vtimes
=
(
int
*
)
Malloc
(
nt
*
sizeof
(
i
nt
)
)
;
vtimes
.
resize
(
nt
);
for
(
i
=
0
;
i
<
nt
;
++
i
)
vtimes
[
i
]
=
i
*
10000
+
45
*
100
;
...
...
@@ -1365,7 +1359,7 @@ Importcmsaf(void *process)
cdoWarning
(
"Wrong time string!"
);
break
;
}
vtimes
[
i
]
=
(
int
)
itime
;
vtimes
[
i
]
=
itime
;
}
}
}
...
...
@@ -1429,7 +1423,7 @@ Importcmsaf(void *process)
get_global_att
(
file_id
,
"/"
,
vlistID
);
if
(
dsets
.
lmetadata
)
get_global_att
(
file_id
,
"Metadata"
,
vlistID
);
vdate
=
get_vdate
(
vlistID
);
int
vdate
=
get_vdate
(
vlistID
);
if
(
vdate
==
0
)
vdate
=
10101
;
const
auto
streamID
=
cdoOpenWrite
(
1
);
...
...
@@ -1439,8 +1433,7 @@ Importcmsaf(void *process)
for
(
tsID
=
0
;
tsID
<
nt
;
++
tsID
)
{
taxisDefVdate
(
taxisID
,
vdate
);
vtime
=
0
;
if
(
vtimes
)
vtime
=
vtimes
[
tsID
];
const
int
vtime
=
(
vtimes
.
empty
())
?
0
:
vtimes
[
tsID
];
taxisDefVtime
(
taxisID
,
vtime
);
cdoDefTimestep
(
streamID
,
tsID
);
...
...
@@ -1457,7 +1450,7 @@ Importcmsaf(void *process)
{
offset
=
gridsize
*
levelID
;
if
(
nz
==
1
)
offset
=
gridsize
*
tsID
;
array
=
dsets
.
obj
[
ivar
].
array
+
offset
;
double
*
array
=
dsets
.
obj
[
ivar
].
array
+
offset
;
const
auto
mm
=
varrayMinMaxMV
(
gridsize
,
array
,
missval
);
nmiss
=
gridsize
-
mm
.
n
;
...
...
@@ -1487,8 +1480,6 @@ Importcmsaf(void *process)
for
(
ivar
=
0
;
ivar
<
dsets
.
nsets
;
++
ivar
)
Free
(
dsets
.
obj
[
ivar
].
array
);
if
(
vtimes
)
Free
(
vtimes
);
cdoFinish
();
#else
cdoAbort
(
"HDF5 support not compiled in!"
);
...
...
src/Importobs.cc
View file @
941aa212
...
...
@@ -18,14 +18,13 @@
#include
<cdi.h>
#include
"cdo_options.h"
#include
"dmemory.h"
#include
"process_int.h"
#include
"array.h"
#include
"readline.h"
#include
"cdo_cdi_wrapper.h"
#include
<mpim_grid.h>
#include
"griddes.h"
#define MAX_VARS 6
static
void
init_vars
(
int
vlistID
,
int
gridID
,
int
zaxisID
,
int
nvars
)
...
...
@@ -36,7 +35,7 @@ init_vars(int vlistID, int gridID, int zaxisID, int nvars)
for
(
int
i
=
0
;
i
<
nvars
;
++
i
)
{
int
varID
=
vlistDefVar
(
vlistID
,
gridID
,
zaxisID
,
TIME_VARYING
);
const
auto
varID
=
vlistDefVar
(
vlistID
,
gridID
,
zaxisID
,
TIME_VARYING
);
vlistDefVarParam
(
vlistID
,
varID
,
cdiEncodeParam
(
code
[
i
],
255
,
255
));
cdiDefKeyString
(
vlistID
,
varID
,
CDI_KEY_NAME
,
name
[
i
]);
cdiDefKeyString
(
vlistID
,
varID
,
CDI_KEY_UNITS
,
units
[
i
]);
...
...
@@ -45,7 +44,7 @@ init_vars(int vlistID, int gridID, int zaxisID, int nvars)
}
static
void
init_data
(
int
vlistID
,
int
nvars
,
double
*
data
[]
)
init_data
(
int
vlistID
,
int
nvars
,
Varray2D
<
double
>
&
data
)
{
for
(
int
varID
=
0
;
varID
<
nvars
;
++
varID
)
{
...
...
@@ -56,15 +55,15 @@ init_data(int vlistID, int nvars, double *data[])
}
static
void
write_data
(
CdoStreamID
streamID
,
int
vlistID
,
int
nvars
,
double
*
data
[]
)
write_data
(
CdoStreamID
streamID
,
int
vlistID
,
int
nvars
,
Varray2D
<
double
>
&
data
)
{
for
(
int
varID
=
0
;
varID
<
nvars
;
++
varID
)
{
const
auto
gridsize
=
gridInqSize
(
vlistInqVarGrid
(
vlistID
,
varID
));
const
auto
missval
=
vlistInqVarMissval
(
vlistID
,
varID
);
const
auto
nmiss
=
arrayNumMV
(
gridsize
,
data
[
varID
],
missval
);
const
auto
nmiss
=
v
arrayNumMV
(
gridsize
,
data
[
varID
],
missval
);
cdoDefRecord
(
streamID
,
varID
,
0
);
cdoWriteRecord
(
streamID
,
data
[
varID
],
nmiss
);
cdoWriteRecord
(
streamID
,
data
[
varID
]
.
data
()
,
nmiss
);
}
}
...
...
@@ -82,17 +81,14 @@ void *
Importobs
(
void
*
process
)
{
char
line
[
MAX_LINE_LEN
];
int
i
,
j
;
int
nvars
=
MAX_VARS
;
size_t
i
,
j
;
long
index
;
constexpr
int
nvars
=
6
;
int
vtime
=
0
;
double
*
data
[
MAX_VARS
];
char
dummy
[
32
],
station
[
32
],
datetime
[
32
];
float
lat
,
lon
,
height1
,
pressure
,
height2
,
value
;
double
latmin
=
90
,
latmax
=
-
90
,
lonmin
=
360
,
lonmax
=
-
360
;
int
code
;
int
index
;
double
dx
,
dy
;
char
*
pstation
;
cdoInitialize
(
process
);
...
...
@@ -107,10 +103,8 @@ Importobs(void *process)
if
(
gridInqType
(
gridID
)
!=
GRID_LONLAT
)
cdoAbort
(
"Unsupported grid type: %s"
,
gridNamePtr
(
gridInqType
(
gridID
)));
const
auto
gridsize
=
gridInqSize
(
gridID
);
const
int
xsize
=
gridInqXsize
(
gridID
);
const
int
ysize
=
gridInqYsize
(
gridID
);
// printf("gridsize=%d, xsize=%d, ysize=%d\n", gridsize, xsize, ysize);
const
auto
xsize
=
gridInqXsize
(
gridID
);
const
auto
ysize
=
gridInqYsize
(
gridID
);
Varray
<
double
>
xvals
(
gridsize
),
yvals
(
gridsize
);
gridInqXvals
(
gridID
,
xvals
.
data
());
...
...
@@ -139,7 +133,8 @@ Importobs(void *process)
const
auto
taxisID
=
cdoTaxisCreate
(
TAXIS_ABSOLUTE
);
vlistDefTaxis
(
vlistID
,
taxisID
);
for
(
i
=
0
;
i
<
nvars
;
++
i
)
data
[
i
]
=
(
double
*
)
Malloc
(
gridsize
*
sizeof
(
double
));
Varray2D
<
double
>
data
(
nvars
);
for
(
i
=
0
;
i
<
nvars
;
++
i
)
data
[
i
].
resize
(
gridsize
);
init_vars
(
vlistID
,
gridID
,
zaxisID
,
nvars
);
...
...
@@ -178,11 +173,11 @@ Importobs(void *process)
if
(
lat
<
latmin
)
latmin
=
lat
;
if
(
lat
>
latmax
)
latmax
=
lat
;
dy
=
yvals
[
1
]
-
yvals
[
0
];
const
auto
dy
=
yvals
[
1
]
-
yvals
[
0
];
for
(
j
=
0
;
j
<
ysize
;
++
j
)
if
(
lat
>=
(
yvals
[
j
]
-
dy
/
2
)
&&
lat
<
(
yvals
[
j
]
+
dy
/
2
))
break
;
dx
=
xvals
[
1
]
-
xvals
[
0
];
const
auto
dx
=
xvals
[
1
]
-
xvals
[
0
];
if
(
lon
<
(
xvals
[
0
]
-
dx
/
2
)
&&
lon
<
0
)
lon
+=
360
;
for
(
i
=
0
;
i
<
xsize
;
++
i
)
if
(
lon
>=
(
xvals
[
i
]
-
dx
/
2
)
&&
lon
<
(
xvals
[
i
]
+
dx
/
2
))
break
;
...
...
@@ -196,7 +191,7 @@ Importobs(void *process)
// printf("%d %d %d %g %g %g %g\n", i, j, index, dx, dy, lon, lat);
if
(
i
<
xsize
&&
j
<
ysize
&&
index
>=
0
)
{
pstation
=
station
;
char
*
pstation
=
station
;
while
(
isalpha
(
*
pstation
))
pstation
++
;
// printf("station %s %d\n", pstation, atoi(pstation));
data
[
index
][
j
*
xsize
+
i
]
=
value
;
...
...
@@ -212,8 +207,6 @@ Importobs(void *process)
write_data
(
streamID
,
vlistID
,
nvars
,
data
);
for
(
i
=
0
;
i
<
nvars
;
++
i
)
Free
(
data
[
i
]);
if
(
Options
::
cdoVerbose
)
printf
(
"lonmin=%g, lonmax=%g, latmin=%g, latmax=%g
\n
"
,
lonmin
,
lonmax
,
latmin
,
latmax
);
processDefVarNum
(
vlistNvars
(
vlistID
));
...
...
src/Test.cc
View file @
941aa212
...
...
@@ -17,8 +17,8 @@
#include
<cdi.h>
#include
"dmemory.h"
#include
"process_int.h"
#include
"array.h"
void
*
Test
(
void
*
process
)
...
...
@@ -83,11 +83,11 @@ Testdata(void *process)
cdoDefVlist
(
streamID2
,
vlistID2
);
auto
gridsize
=
vlistGridsizeMax
(
vlistID1
);
double
*
array
=
(
double
*
)
Malloc
(
gridsize
*
sizeof
(
double
)
);
float
*
fval
=
(
float
*
)
Malloc
(
gridsize
*
sizeof
(
float
)
);
int
*
ival
=
(
int
*
)
Malloc
(
gridsize
*
sizeof
(
int
)
);
unsigned
char
*
cval
=
(
unsigned
char
*
)
Malloc
(
gridsize
*
sizeof
(
unsigned
char
)
*
4
);
unsigned
char
*
cval2
=
(
unsigned
char
*
)
Malloc
(
gridsize
*
sizeof
(
unsigned
char
)
*
4
);
V
array
<
double
>
array
(
gridsize
);