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 Monarith monadd Add monthly time series
12 Monarith monsub Subtract monthly time series
13 Monarith monmul Multiply monthly time series
14 Monarith mondiv Divide monthly time series
15 */
16
17 #include <cdi.h>
18
19 #include "cdo_options.h"
20 #include "process_int.h"
21 #include "cdo_vlist.h"
22
23 void *
Monarith(void * process)24 Monarith(void *process)
25 {
26 cdo_initialize(process);
27
28 cdo_operator_add("monadd", FieldFunc_Add, 0, nullptr);
29 cdo_operator_add("monsub", FieldFunc_Sub, 0, nullptr);
30 cdo_operator_add("monmul", FieldFunc_Mul, 0, nullptr);
31 cdo_operator_add("mondiv", FieldFunc_Div, 0, nullptr);
32
33 const auto operatorID = cdo_operator_id();
34 const auto operfunc = cdo_operator_f1(operatorID);
35
36 operator_check_argc(0);
37
38 const auto streamID1 = cdo_open_read(0);
39 const auto streamID2 = cdo_open_read(1);
40
41 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
42 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
43 const auto vlistID3 = vlistDuplicate(vlistID1);
44
45 vlist_compare(vlistID1, vlistID2, CMP_ALL);
46
47 VarList varList1;
48 varListInit(varList1, vlistID1);
49
50 const auto gridsizemax = vlistGridsizeMax(vlistID1);
51
52 Field field1, field2;
53 field1.resize(gridsizemax);
54 field2.resize(gridsizemax);
55
56 const auto taxisID1 = vlistInqTaxis(vlistID1);
57 const auto taxisID2 = vlistInqTaxis(vlistID2);
58 const auto taxisID3 = taxisDuplicate(taxisID1);
59 vlistDefTaxis(vlistID3, taxisID3);
60
61 const auto streamID3 = cdo_open_write(2);
62 cdo_def_vlist(streamID3, vlistID3);
63
64 FieldVector2D vars2;
65 fields_from_vlist(vlistID2, vars2, FIELD_VEC);
66
67 int yearmon2 = -1;
68 int tsID = 0;
69 int tsID2 = 0;
70 while (true)
71 {
72 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
73 if (nrecs == 0) break;
74
75 auto vdate = taxisInqVdate(taxisID1);
76 const int yearmon1 = vdate / 100;
77
78 if (yearmon1 != yearmon2)
79 {
80 const int year1 = yearmon1 / 100;
81 const int mon1 = yearmon1 - (yearmon1 / 100) * 100;
82 if (Options::cdoVerbose) cdo_print("Process: Year = %4d Month = %2d", year1, mon1);
83
84 const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID2);
85 if (nrecs2 == 0) cdo_abort("Missing year=%4d mon=%2d in %s!", year1, mon1, cdo_get_stream_name(1));
86
87 vdate = taxisInqVdate(taxisID2);
88 yearmon2 = vdate / 100;
89
90 if (yearmon1 != yearmon2)
91 {
92 const int year2 = yearmon2 / 100;
93 const int mon2 = yearmon2 - (yearmon2 / 100) * 100;
94 cdo_abort("Timestep %d in %s has wrong date! Current year=%4d mon=%2d, expected year=%4d mon=%2d", tsID2 + 1,
95 cdo_get_stream_name(1), year2, mon2, year1, mon1);
96 }
97
98 for (int recID = 0; recID < nrecs2; recID++)
99 {
100 size_t nmiss;
101 int varID, levelID;
102 cdo_inq_record(streamID2, &varID, &levelID);
103 cdo_read_record(streamID2, vars2[varID][levelID].vec_d.data(), &nmiss);
104 vars2[varID][levelID].nmiss = nmiss;
105 }
106
107 tsID2++;
108 }
109
110 cdo_taxis_copy_timestep(taxisID3, taxisID1);
111 cdo_def_timestep(streamID3, tsID);
112
113 for (int recID = 0; recID < nrecs; recID++)
114 {
115 int varID, levelID;
116 cdo_inq_record(streamID1, &varID, &levelID);
117 cdo_read_record(streamID1, field1.vec_d.data(), &field1.nmiss);
118 field1.grid = varList1[varID].gridID;
119 field1.missval = varList1[varID].missval;
120
121 field2_function(field1, vars2[varID][levelID], operfunc);
122
123 cdo_def_record(streamID3, varID, levelID);
124 cdo_write_record(streamID3, field1.vec_d.data(), field1.nmiss);
125 }
126
127 tsID++;
128 }
129
130 cdo_stream_close(streamID3);
131 cdo_stream_close(streamID2);
132 cdo_stream_close(streamID1);
133
134 cdo_finish();
135
136 return nullptr;
137 }
138