1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Copyright (C) 2007 Brockmann Consult
5
6 Author: Ralf Quast
7 Uwe Schulzweida
8
9 */
10
11 /*
12 This module contains the following operators:
13
14 Seascount seascount Seasonal counts
15 */
16
17 #include <cdi.h>
18
19 #include "cdo_season.h"
20 #include "datetime.h"
21 #include "process_int.h"
22
23 void *
Seascount(void * process)24 Seascount(void *process)
25 {
26 int64_t vdate0 = 0;
27 int vtime0 = 0;
28 int seas0 = 0;
29 int oldmon = 0;
30
31 cdo_initialize(process);
32
33 cdo_operator_add("seascount", 0, 0, nullptr);
34
35 operator_check_argc(0);
36
37 const auto season_start = get_season_start();
38
39 const auto streamID1 = cdo_open_read(0);
40
41 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
42 const auto vlistID2 = vlistDuplicate(vlistID1);
43
44 const auto taxisID1 = vlistInqTaxis(vlistID1);
45 const auto taxisID2 = taxisDuplicate(taxisID1);
46 vlistDefTaxis(vlistID2, taxisID2);
47
48 const auto streamID2 = cdo_open_write(1);
49 cdo_def_vlist(streamID2, vlistID2);
50
51 const auto maxrecs = vlistNrecs(vlistID1);
52 std::vector<RecordInfo> recList(maxrecs);
53
54 auto gridsizemax = vlistGridsizeMax(vlistID1);
55 if (vlistNumber(vlistID1) != CDI_REAL) gridsizemax *= 2;
56
57 Field field;
58 field.resize(gridsizemax);
59
60 FieldVector2D vars1;
61 fields_from_vlist(vlistID1, vars1, FIELD_VEC);
62
63 int tsID = 0;
64 int otsID = 0;
65 while (true)
66 {
67 int nrecs = 0;
68 int nsets = 0;
69 bool newseas = false;
70 while (true)
71 {
72 nrecs = cdo_stream_inq_timestep(streamID1, tsID);
73 if (nrecs == 0) break;
74
75 const auto vdate = taxisInqVdate(taxisID1);
76 const auto vtime = taxisInqVtime(taxisID1);
77
78 const auto month = decode_month(vdate);
79 auto newmon = month;
80 if (season_start == START_DEC && newmon == 12) newmon = 0;
81
82 const auto seas = month_to_season(month);
83
84 if (nsets == 0)
85 {
86 seas0 = seas;
87 oldmon = newmon;
88 }
89
90 if (newmon < oldmon) newseas = true;
91
92 if ((seas != seas0) || newseas) break;
93
94 oldmon = newmon;
95
96 for (int recID = 0; recID < nrecs; recID++)
97 {
98 int varID, levelID;
99 cdo_inq_record(streamID1, &varID, &levelID);
100
101 if (tsID == 0)
102 {
103 recList[recID].varID = varID;
104 recList[recID].levelID = levelID;
105 recList[recID].lconst = (vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT);
106 }
107
108 const auto fieldsize = vars1[varID][levelID].size;
109
110 if (nsets == 0)
111 {
112 for (size_t i = 0; i < fieldsize; i++) vars1[varID][levelID].vec_d[i] = vars1[varID][levelID].missval;
113 vars1[varID][levelID].nmiss = fieldsize;
114 }
115
116 cdo_read_record(streamID1, field.vec_d.data(), &field.nmiss);
117 field.size = vars1[varID][levelID].size;
118 field.grid = vars1[varID][levelID].grid;
119 field.missval = vars1[varID][levelID].missval;
120
121 field2_count(vars1[varID][levelID], field);
122 }
123
124 vdate0 = vdate;
125 vtime0 = vtime;
126 nsets++;
127 tsID++;
128 }
129
130 if (nrecs == 0 && nsets == 0) break;
131
132 taxisDefVdate(taxisID2, vdate0);
133 taxisDefVtime(taxisID2, vtime0);
134 cdo_def_timestep(streamID2, otsID);
135
136 for (int recID = 0; recID < maxrecs; recID++)
137 {
138 if (otsID && recList[recID].lconst) continue;
139
140 const auto varID = recList[recID].varID;
141 const auto levelID = recList[recID].levelID;
142 cdo_def_record(streamID2, varID, levelID);
143 cdo_write_record(streamID2, vars1[varID][levelID].vec_d.data(), vars1[varID][levelID].nmiss);
144 }
145
146 if (nrecs == 0) break;
147 otsID++;
148 }
149
150 cdo_stream_close(streamID2);
151 cdo_stream_close(streamID1);
152
153 cdo_finish();
154
155 return nullptr;
156 }
157