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 <mpim_grid.h>
17 #include "cdi_lockedIO.h"
18
19 static void
varrms(const Varray<double> & w,const FieldVector & field1,const FieldVector & field2,Field & field3)20 varrms(const Varray<double> &w, const FieldVector &field1, const FieldVector &field2, Field &field3)
21 {
22 const auto grid1 = field1[0].grid;
23 const auto grid2 = field2[0].grid;
24 const auto missval1 = field1[0].missval;
25 const auto missval2 = field2[0].missval;
26 double rsum = 0, rsumw = 0;
27
28 const auto nlev = field1.size();
29 const auto len = gridInqSize(grid1);
30 if (len != gridInqSize(grid2)) cdo_abort("fields have different size!");
31
32 // if ( nmiss1 )
33 {
34 for (size_t k = 0; k < nlev; k++)
35 {
36 const auto array1 = field1[k].vec_d;
37 const auto array2 = field2[k].vec_d;
38 for (size_t i = 0; i < len; i++) /* if ( !DBL_IS_EQUAL(w[i], missval1) ) */
39 {
40 rsum = ADDMN(rsum, MULMN(w[i], MULMN(SUBMN(array2[i], array1[i]), SUBMN(array2[i], array1[i]))));
41 rsumw = ADDMN(rsumw, w[i]);
42 }
43 }
44 }
45 /*
46 else
47 {
48 for ( i = 0; i < len; i++ )
49 {
50 rsum += w[i] * array1[i];
51 rsumw += w[i];
52 }
53 }
54 */
55
56 const double ravg = SQRTMN(DIVMN(rsum, rsumw));
57
58 size_t rnmiss = 0;
59 if (DBL_IS_EQUAL(ravg, missval1)) rnmiss++;
60
61 field3.vec_d[0] = ravg;
62 field3.nmiss = rnmiss;
63 }
64
65 void *
Varrms(void * process)66 Varrms(void *process)
67 {
68 int lastgrid = -1;
69 int oldcode = 0;
70
71 cdo_initialize(process);
72
73 operator_check_argc(0);
74
75 const auto needWeights = true;
76
77 const auto streamID1 = cdo_open_read(0);
78 const auto streamID2 = cdo_open_read(1);
79
80 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
81 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
82
83 double slon = 0.0, slat = 0.0;
84 const auto gridID3 = gridCreate(GRID_LONLAT, 1);
85 gridDefXsize(gridID3, 1);
86 gridDefYsize(gridID3, 1);
87 gridDefXvals(gridID3, &slon);
88 gridDefYvals(gridID3, &slat);
89
90 vlistClearFlag(vlistID1);
91 const auto nvars = vlistNvars(vlistID1);
92 for (int varID = 0; varID < nvars; varID++) vlistDefFlag(vlistID1, varID, 0, true);
93
94 auto vlistID3 = vlistCreate();
95 cdo_vlist_copy_flag(vlistID3, vlistID1);
96
97 const auto taxisID1 = vlistInqTaxis(vlistID1);
98 const auto taxisID3 = taxisDuplicate(taxisID1);
99 vlistDefTaxis(vlistID3, taxisID3);
100
101 const auto ngrids = vlistNgrids(vlistID1);
102 int index = 0;
103 const auto gridID1 = vlistGrid(vlistID1, index);
104
105 if (needWeights && gridInqType(gridID1) != GRID_LONLAT && gridInqType(gridID1) != GRID_GAUSSIAN)
106 cdo_abort("Unsupported gridtype: %s", gridNamePtr(gridInqType(gridID1)));
107
108 vlistChangeGridIndex(vlistID3, index, gridID3);
109 if (ngrids > 1) cdo_abort("Too many different grids!");
110
111 const auto streamID3 = cdo_open_write(2);
112 cdo_def_vlist(streamID3, vlistID3);
113
114 FieldVector2D vars1, vars2;
115 fields_from_vlist(vlistID1, vars1, FIELD_VEC);
116 fields_from_vlist(vlistID2, vars2, FIELD_VEC);
117
118 const auto gridsizemax = vlistGridsizeMax(vlistID1);
119 Varray<double> weights;
120 if (needWeights) weights.resize(gridsizemax);
121
122 Field field3;
123 field3.resize(1);
124 field3.grid = gridID3;
125
126 int tsID = 0;
127 while (true)
128 {
129 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
130 if (nrecs == 0) break;
131
132 const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
133 if (nrecs2 == 0) cdo_abort("Input streams have different number of timesteps!");
134
135 cdo_taxis_copy_timestep(taxisID3, taxisID1);
136 cdo_def_timestep(streamID3, tsID);
137
138 for (int recID = 0; recID < nrecs; recID++)
139 {
140 size_t nmiss;
141 int varID, levelID;
142 cdo_inq_record(streamID1, &varID, &levelID);
143 cdo_read_record(streamID1, vars1[varID][levelID].vec_d.data(), &nmiss);
144 if (nmiss) cdo_abort("Missing values unsupported for this operator!");
145
146 cdo_inq_record(streamID2, &varID, &levelID);
147 cdo_read_record(streamID2, vars2[varID][levelID].vec_d.data(), &nmiss);
148 if (nmiss) cdo_abort("Missing values unsupported for this operator!");
149 }
150
151 for (int varID = 0; varID < nvars; varID++)
152 {
153 auto wstatus = false;
154 const auto gridID = vars1[varID][0].grid;
155 if (needWeights && gridID != lastgrid)
156 {
157 lastgrid = gridID;
158 wstatus = gridWeights(gridID, weights.data());
159 }
160 const auto code = vlistInqVarCode(vlistID1, varID);
161 if (wstatus != 0 && tsID == 0 && code != oldcode) cdo_warning("Using constant area weights for code %d!", oldcode = code);
162
163 field3.missval = vars1[varID][0].missval;
164 varrms(weights, vars1[varID], vars2[varID], field3);
165
166 cdo_def_record(streamID3, varID, 0);
167 cdo_write_record(streamID3, field3.vec_d.data(), field3.nmiss);
168 }
169
170 tsID++;
171 }
172
173 cdo_stream_close(streamID3);
174 cdo_stream_close(streamID2);
175 cdo_stream_close(streamID1);
176
177 vlistDestroy(vlistID3);
178
179 cdo_finish();
180
181 return nullptr;
182 }
183