1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 /*
9 This module contains the following operators:
10
11 */
12
13 #include <cdi.h>
14
15 #include "process_int.h"
16 #include "cdo_vlist.h"
17 #include <mpim_grid.h>
18
19 void genGridIndex(int gridID1, int gridID2, std::vector<long> &index);
20
21 void *
Mergegrid(void * process)22 Mergegrid(void *process)
23 {
24 cdo_initialize(process);
25
26 operator_check_argc(0);
27
28 const auto streamID1 = cdo_open_read(0);
29
30 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
31 const auto taxisID1 = vlistInqTaxis(vlistID1);
32 const auto taxisID3 = taxisDuplicate(taxisID1);
33
34 const auto streamID2 = cdo_open_read(1);
35 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
36
37 vlist_compare(vlistID1, vlistID2, CMP_NAME | CMP_NLEVEL);
38
39 int ndiffgrids = 0;
40 for (int index = 1; index < vlistNgrids(vlistID1); index++)
41 if (vlistGrid(vlistID1, 0) != vlistGrid(vlistID1, index)) ndiffgrids++;
42
43 if (ndiffgrids > 0) cdo_abort("Too many different grids in %s!", cdo_get_stream_name(0));
44
45 ndiffgrids = 0;
46 for (int index = 1; index < vlistNgrids(vlistID2); index++)
47 if (vlistGrid(vlistID2, 0) != vlistGrid(vlistID2, index)) ndiffgrids++;
48
49 if (ndiffgrids > 0) cdo_abort("Too many different grids in %s!", cdo_get_stream_name(1));
50
51 const auto gridID1 = vlistGrid(vlistID1, 0);
52 const auto gridID2 = vlistGrid(vlistID2, 0);
53
54 const auto gridsize1 = gridInqSize(gridID1);
55 const auto gridsize2 = gridInqSize(gridID2);
56
57 Varray<double> array1(gridsize1), array2(gridsize2);
58 std::vector<long> gindex(gridsize2);
59
60 genGridIndex(gridID1, gridID2, gindex);
61
62 const auto vlistID3 = vlistDuplicate(vlistID1);
63 const auto streamID3 = cdo_open_write(2);
64
65 vlistDefTaxis(vlistID3, taxisID3);
66 cdo_def_vlist(streamID3, vlistID3);
67
68 int tsID = 0;
69 while (true)
70 {
71 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
72 if (nrecs == 0) break;
73
74 cdo_taxis_copy_timestep(taxisID3, taxisID1);
75
76 const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
77 if (nrecs2 == 0) cdo_abort("Input streams have different number of timesteps!");
78
79 if (nrecs != nrecs2) cdo_abort("Input streams have different number of records!");
80
81 cdo_def_timestep(streamID3, tsID);
82
83 for (int recID = 0; recID < nrecs; recID++)
84 {
85 int varID, levelID;
86 cdo_inq_record(streamID2, &varID, &levelID);
87 size_t nmiss2;
88 cdo_read_record(streamID2, array2.data(), &nmiss2);
89
90 const auto missval2 = vlistInqVarMissval(vlistID2, varID);
91
92 cdo_inq_record(streamID1, &varID, &levelID);
93 size_t nmiss1;
94 cdo_read_record(streamID1, array1.data(), &nmiss1);
95
96 const auto missval1 = vlistInqVarMissval(vlistID1, varID);
97
98 for (size_t i = 0; i < gridsize2; i++)
99 {
100 if (gindex[i] >= 0 && !DBL_IS_EQUAL(array2[i], missval2))
101 {
102 array1[gindex[i]] = array2[i];
103 }
104 }
105
106 if (nmiss1)
107 {
108 nmiss1 = 0;
109 for (size_t i = 0; i < gridsize1; i++)
110 if (DBL_IS_EQUAL(array1[i], missval1)) nmiss1++;
111 }
112
113 cdo_def_record(streamID3, varID, levelID);
114 cdo_write_record(streamID3, array1.data(), nmiss1);
115 }
116
117 tsID++;
118 }
119
120 cdo_stream_close(streamID3);
121 cdo_stream_close(streamID2);
122 cdo_stream_close(streamID1);
123
124 cdo_finish();
125
126 return nullptr;
127 }
128