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 #include <cdi.h>
9
10 #include "process_int.h"
11
12 void *
Complextorect(void * process)13 Complextorect(void *process)
14 {
15 cdo_initialize(process);
16
17 const auto COMPLEXTORECT = cdo_operator_add("complextorect", 0, 0, nullptr);
18 const auto COMPLEXTOPOL = cdo_operator_add("complextopol", 0, 0, nullptr);
19
20 const auto operatorID = cdo_operator_id();
21
22 operator_check_argc(0);
23
24 const auto streamID1 = cdo_open_read(0);
25
26 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
27 const auto vlistID2 = vlistDuplicate(vlistID1);
28 const auto vlistID3 = vlistDuplicate(vlistID1);
29
30 const auto nvars = vlistNvars(vlistID2);
31 for (int varID = 0; varID < nvars; ++varID)
32 {
33 auto datatype = vlistInqVarDatatype(vlistID2, varID);
34 datatype = (datatype == CDI_DATATYPE_CPX64) ? CDI_DATATYPE_FLT64 : CDI_DATATYPE_FLT32;
35 vlistDefVarDatatype(vlistID2, varID, datatype);
36 vlistDefVarDatatype(vlistID3, varID, datatype);
37 }
38
39 const auto taxisID1 = vlistInqTaxis(vlistID1);
40 const auto taxisID2 = taxisDuplicate(taxisID1);
41 const auto taxisID3 = taxisDuplicate(taxisID1);
42 vlistDefTaxis(vlistID2, taxisID2);
43 vlistDefTaxis(vlistID3, taxisID3);
44
45 const auto streamID2 = cdo_open_write(1);
46 const auto streamID3 = cdo_open_write(2);
47
48 cdo_def_vlist(streamID2, vlistID2);
49 cdo_def_vlist(streamID3, vlistID3);
50
51 VarList varList1;
52 varListInit(varList1, vlistID1);
53
54 const auto gridsizemax = vlistGridsizeMax(vlistID1);
55 Varray<double> array1(2 * gridsizemax);
56 Varray<double> array2(gridsizemax), array3(gridsizemax);
57
58 int tsID = 0;
59 while (true)
60 {
61 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
62 if (nrecs == 0) break;
63
64 cdo_taxis_copy_timestep(taxisID2, taxisID1);
65 cdo_taxis_copy_timestep(taxisID3, taxisID1);
66
67 cdo_def_timestep(streamID2, tsID);
68 cdo_def_timestep(streamID3, tsID);
69
70 for (int recID = 0; recID < nrecs; recID++)
71 {
72 int varID, levelID;
73 cdo_inq_record(streamID1, &varID, &levelID);
74 cdo_def_record(streamID2, varID, levelID);
75 cdo_def_record(streamID3, varID, levelID);
76
77 const auto gridsize = varList1[varID].gridsize;
78
79 size_t nmiss;
80 cdo_read_record(streamID1, array1.data(), &nmiss);
81
82 const auto missval1 = varList1[varID].missval;
83 const auto missval2 = missval1;
84
85 if (operatorID == COMPLEXTORECT)
86 {
87 for (size_t i = 0; i < gridsize; ++i)
88 {
89 array2[i] = array1[2 * i];
90 array3[i] = array1[2 * i + 1];
91 }
92 }
93 else if (operatorID == COMPLEXTOPOL)
94 {
95 for (size_t i = 0; i < gridsize; ++i)
96 {
97 array2[i] = SQRTMN(ADDMN(MULMN(array1[2 * i], array1[2 * i]), MULMN(array1[2 * i + 1], array1[2 * i + 1])));
98 array3[i] = (DBL_IS_EQUAL(array1[2 * i], missval1) || DBL_IS_EQUAL(array1[2 * i + 1], missval1))
99 ? missval1
100 : std::atan2(array1[2 * i + 1], array1[2 * i]);
101 }
102 }
103
104 cdo_write_record(streamID2, array2.data(), nmiss);
105 cdo_write_record(streamID3, array3.data(), nmiss);
106 }
107
108 tsID++;
109 }
110
111 cdo_stream_close(streamID3);
112 cdo_stream_close(streamID2);
113 cdo_stream_close(streamID1);
114
115 vlistDestroy(vlistID2);
116 vlistDestroy(vlistID3);
117
118 cdo_finish();
119
120 return nullptr;
121 }
122