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 <algorithm>
9 #include <climits>
10
11 #include <cdi.h>
12
13 #include "cdo_options.h"
14 #include "process_int.h"
15 #include "param_conversion.h"
16 #include "libncl.h"
17 #include "pmlist.h"
18
19 static void
uv2dv_cfd_W(double missval,double * u,double * v,double * lon,double * lat,size_t nlon,size_t nlat,size_t nlev,int boundOpt,double * div)20 uv2dv_cfd_W(double missval, double *u, double *v, double *lon, double *lat, size_t nlon, size_t nlat, size_t nlev, int boundOpt,
21 double *div)
22 {
23 int ierror;
24
25 // Test dimension sizes.
26 if ((nlon > INT_MAX) || (nlat > INT_MAX)) cdo_abort("nlat and/or nlon is greater than INT_MAX!");
27
28 int inlon = (int) nlon;
29 int inlat = (int) nlat;
30
31 size_t gridsize_uv = nlat * nlon;
32
33 for (size_t k = 0; k < nlev; ++k)
34 {
35 double *tmp_u = u + k * gridsize_uv;
36 double *tmp_v = v + k * gridsize_uv;
37 double *tmp_div = div + k * gridsize_uv;
38 // Init output array.
39 varray_fill(gridsize_uv, tmp_div, 0.0);
40 // Call the Fortran routine.
41 #ifdef HAVE_CF_INTERFACE
42 DDVFIDF(tmp_u, tmp_v, lat, lon, inlon, inlat, missval, boundOpt, tmp_div, ierror);
43 #else
44 cdo_abort("Fortran support not compiled in!");
45 #endif
46 }
47 }
48
49 static void
uv2vr_cfd_W(double missval,double * u,double * v,double * lon,double * lat,size_t nlon,size_t nlat,size_t nlev,int boundOpt,double * vort)50 uv2vr_cfd_W(double missval, double *u, double *v, double *lon, double *lat, size_t nlon, size_t nlat, size_t nlev, int boundOpt,
51 double *vort)
52 {
53 int ierror;
54
55 // Test dimension sizes.
56 if ((nlon > INT_MAX) || (nlat > INT_MAX)) cdo_abort("nlat and/or nlon is greater than INT_MAX!");
57
58 int inlon = (int) nlon;
59 int inlat = (int) nlat;
60
61 size_t gridsize_uv = nlat * nlon;
62
63 for (size_t k = 0; k < nlev; ++k)
64 {
65 double *tmp_u = u + k * gridsize_uv;
66 double *tmp_v = v + k * gridsize_uv;
67 double *tmp_vort = vort + k * gridsize_uv;
68 // Init output array.
69 varray_fill(gridsize_uv, tmp_vort, 0.0);
70 // Call the Fortran routine.
71 #ifdef HAVE_CF_INTERFACE
72 DVRFIDF(tmp_u, tmp_v, lat, lon, inlon, inlat, missval, boundOpt, tmp_vort, ierror);
73 #else
74 cdo_abort("Fortran support not compiled in!");
75 #endif
76 }
77 }
78
79 static int
find_name(int vlistID,char * name)80 find_name(int vlistID, char *name)
81 {
82 char varname[CDI_MAX_NAME];
83
84 const auto nvars = vlistNvars(vlistID);
85 for (int varID = 0; varID < nvars; ++varID)
86 {
87 vlistInqVarName(vlistID, varID, varname);
88 if (cdo_cmpstr(name, varname)) return varID;
89 }
90
91 return CDI_UNDEFID;
92 }
93
94 enum struct OutMode
95 {
96 NEW,
97 APPEND,
98 REPLACE
99 };
100
101 // Parameter
102 OutMode outMode(OutMode::NEW);
103 int boundOpt = -1;
104 char name_u[CDI_MAX_NAME], name_v[CDI_MAX_NAME];
105
106 static void
print_parameter(void)107 print_parameter(void)
108 {
109 cdo_print("u=%s, v=%s, boundOpt=%d, outMode=%s", name_u, name_v, boundOpt,
110 outMode == OutMode::NEW ? "new" : outMode == OutMode::APPEND ? "append" : "replace");
111 }
112
113 static void
set_parameter(void)114 set_parameter(void)
115 {
116 strcpy(name_u, "u");
117 strcpy(name_v, "v");
118
119 const auto pargc = cdo_operator_argc();
120 if (pargc)
121 {
122 auto pargv = cdo_get_oper_argv();
123
124 KVList kvlist;
125 kvlist.name = "PARAMETER";
126 if (kvlist.parse_arguments(pargc, pargv) != 0) cdo_abort("Parse error!");
127 if (Options::cdoVerbose) kvlist.print();
128
129 for (const auto &kv : kvlist)
130 {
131 const auto &key = kv.key;
132 if (kv.nvalues > 1) cdo_abort("Too many values for parameter key >%s<!", key.c_str());
133 if (kv.nvalues < 1) cdo_abort("Missing value for parameter key >%s<!", key.c_str());
134 const auto &value = kv.values[0];
135
136 // clang-format off
137 if (key == "u") strcpy(name_u, value.c_str());
138 else if (key == "v") strcpy(name_v, value.c_str());
139 else if (key == "boundOpt") boundOpt = parameter_to_int(value);
140 else if (key == "outMode")
141 {
142 if (value == "new") outMode = OutMode::NEW;
143 else if (value == "append") outMode = OutMode::APPEND;
144 else if (value == "replace") outMode = OutMode::REPLACE;
145 else cdo_abort("Invalid parameter key value: outMode=%s (valid are: new/append/replace)", value.c_str());
146 }
147 else cdo_abort("Invalid parameter key >%s<!", key.c_str());
148 // clang-format on
149 }
150 }
151
152 if (Options::cdoVerbose) print_parameter();
153 }
154
155 void *
NCL_wind(void * process)156 NCL_wind(void *process)
157 {
158 cdo_initialize(process);
159
160 const auto UV2DV_CFD = cdo_operator_add("uv2dv_cfd", 0, 0, "[u, v, boundsOpt, outMode]");
161 const auto UV2VR_CFD = cdo_operator_add("uv2vr_cfd", 0, 0, "[u, v, boundsOpt, outMode]");
162
163 const auto operatorID = cdo_operator_id();
164
165 set_parameter();
166
167 const auto streamID1 = cdo_open_read(0);
168
169 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
170 int vlistID2 = CDI_UNDEFID;
171 if (outMode == OutMode::NEW)
172 vlistID2 = vlistCreate();
173 else if (outMode == OutMode::APPEND)
174 vlistID2 = vlistDuplicate(vlistID1);
175 else
176 cdo_abort("outMode=%d unsupported!", outMode);
177
178 const auto varIDu = find_name(vlistID1, name_u);
179 const auto varIDv = find_name(vlistID1, name_v);
180
181 if (varIDu == CDI_UNDEFID) cdo_abort("%s not found!", name_u);
182 if (varIDv == CDI_UNDEFID) cdo_abort("%s not found!", name_v);
183
184 const auto gridIDu = vlistInqVarGrid(vlistID1, varIDu);
185 const auto gridIDv = vlistInqVarGrid(vlistID1, varIDv);
186 const auto gridtype = gridInqType(gridIDu);
187 const auto gridsizeuv = gridInqSize(gridIDu);
188
189 if (!((gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN) && gridtype == gridInqType(gridIDv)))
190 cdo_abort("u and v must be on a regular lonlat or Gaussian grid!");
191
192 if (gridsizeuv != gridInqSize(gridIDv)) cdo_abort("u and v must have the same grid size!");
193
194 if (boundOpt == -1) boundOpt = gridIsCircular(gridIDu) ? 1 : 0;
195 if (Options::cdoVerbose) print_parameter();
196 if (boundOpt < 0 || boundOpt > 3) cdo_abort("Parameter boundOpt=%d out of bounds (0-3)!", boundOpt);
197
198 const auto nlon = gridInqXsize(gridIDu);
199 const auto nlat = gridInqYsize(gridIDu);
200
201 const auto zaxisIDu = vlistInqVarZaxis(vlistID1, varIDu);
202 const auto nlev = zaxisInqSize(zaxisIDu);
203
204 if (nlev != zaxisInqSize(vlistInqVarZaxis(vlistID1, varIDv))) cdo_abort("u and v must have the same number of level!");
205
206 const auto missvalu = vlistInqVarMissval(vlistID1, varIDu);
207 const auto missvalv = vlistInqVarMissval(vlistID1, varIDv);
208
209 const auto timetype = vlistInqVarTimetype(vlistID1, varIDu);
210 const auto varIDo = vlistDefVar(vlistID2, gridIDu, zaxisIDu, timetype);
211 if (operatorID == UV2DV_CFD)
212 {
213 cdiDefKeyString(vlistID2, varIDo, CDI_KEY_NAME, "d");
214 cdiDefKeyString(vlistID2, varIDo, CDI_KEY_LONGNAME, "divergence");
215 cdiDefKeyString(vlistID2, varIDo, CDI_KEY_UNITS, "1/s");
216 }
217 else if (operatorID == UV2VR_CFD)
218 {
219 cdiDefKeyString(vlistID2, varIDo, CDI_KEY_NAME, "vo");
220 cdiDefKeyString(vlistID2, varIDo, CDI_KEY_LONGNAME, "vorticity");
221 cdiDefKeyString(vlistID2, varIDo, CDI_KEY_UNITS, "1/s");
222 }
223
224 vlistDefVarMissval(vlistID2, varIDo, missvalu);
225
226 Varray<double> lon(nlon), lat(nlat);
227 gridInqXvals(gridIDu, &lon[0]);
228 gridInqYvals(gridIDu, &lat[0]);
229
230 const auto taxisID1 = vlistInqTaxis(vlistID1);
231 const auto taxisID2 = taxisDuplicate(taxisID1);
232 vlistDefTaxis(vlistID2, taxisID2);
233
234 const auto streamID2 = cdo_open_write(1);
235
236 cdo_def_vlist(streamID2, vlistID2);
237
238 const auto gridsizemax = vlistGridsizeMax(vlistID1);
239 Varray<double> array(gridsizemax);
240 Varray<double> arrayu(nlev * gridsizeuv);
241 Varray<double> arrayv(nlev * gridsizeuv);
242 Varray<double> arrayo(nlev * gridsizeuv);
243
244 int tsID = 0;
245 while (true)
246 {
247 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
248 if (nrecs == 0) break;
249
250 size_t nmissu = 0, nmissv = 0;
251
252 cdo_taxis_copy_timestep(taxisID2, taxisID1);
253 cdo_def_timestep(streamID2, tsID);
254
255 for (int recID = 0; recID < nrecs; recID++)
256 {
257 int varID, levelID;
258 cdo_inq_record(streamID1, &varID, &levelID);
259 size_t nmiss;
260 cdo_read_record(streamID1, &array[0], &nmiss);
261
262 if (varID == varIDu || varID == varIDv)
263 {
264 if (varID == varIDu)
265 {
266 std::copy_n(&array[0], gridsizeuv, &arrayu[levelID * gridsizeuv]);
267 nmissu += nmiss;
268 }
269 if (varID == varIDv)
270 {
271 std::copy_n(&array[0], gridsizeuv, &arrayv[levelID * gridsizeuv]);
272 nmissv += nmiss;
273 }
274 }
275
276 if (outMode == OutMode::APPEND)
277 {
278 cdo_def_record(streamID2, varID, levelID);
279 cdo_write_record(streamID2, &array[0], nmiss);
280 }
281 }
282
283 if (nmissu != nmissv)
284 {
285 cdo_abort("u and v have different number of missing values!");
286 if (nmissu && !DBL_IS_EQUAL(missvalu, missvalv))
287 {
288 for (int levelID = 0; levelID < nlev; ++levelID)
289 {
290 auto parray = &arrayv[levelID * gridsizeuv];
291 for (size_t i = 0; i < gridsizeuv; ++i)
292 if (DBL_IS_EQUAL(parray[i], missvalv)) parray[i] = missvalu;
293 }
294 }
295 }
296
297 if (operatorID == UV2DV_CFD)
298 uv2dv_cfd_W(missvalu, &arrayu[0], &arrayv[0], &lon[0], &lat[0], nlon, nlat, nlev, boundOpt, &arrayo[0]);
299 else if (operatorID == UV2VR_CFD)
300 uv2vr_cfd_W(missvalu, &arrayu[0], &arrayv[0], &lon[0], &lat[0], nlon, nlat, nlev, boundOpt, &arrayo[0]);
301
302 for (int levelID = 0; levelID < nlev; ++levelID)
303 {
304 auto parray = &arrayo[levelID * gridsizeuv];
305 auto nmiss = array_num_mv(gridsizeuv, parray, missvalu);
306 cdo_def_record(streamID2, varIDo, levelID);
307 cdo_write_record(streamID2, parray, nmiss);
308 }
309
310 tsID++;
311 }
312
313 cdo_stream_close(streamID1);
314 cdo_stream_close(streamID2);
315
316 vlistDestroy(vlistID2);
317
318 cdo_finish();
319
320 return nullptr;
321 }
322