YAC 3.21.0
Yet Another Coupler
Loading...
Searching...
No Matches
grid_file_common.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#include <string.h>
6#include <math.h>
7#include <netcdf.h>
8
9#include "tests.h"
10#include "grid_file_common.h"
11#include "io_utils.h"
12#include "utils_common.h"
13#include "geometry.h"
14
15static int utest_compare_cell_coords(double * a, double * b, int count) {
16
17 for (int order = -1; order <= 1; order += 2) {
18 for (int start = 0; start < count; ++start) {
19
20 int differences = 0;
21
22 for (int i = 0; i < count; ++i) {
23
24 int j = (count + start + order * i) % count;
25
26 if (fabs(a[i] - b[j]) > 1e-6) ++differences;
27 }
28
29 if (!differences) return 0;
30 }
31 }
32 return 1;
33}
34
35static int utest_compare_ints(int * a, int * b, int count) {
36
37 for (int order = -1; order <= 1; order += 2) {
38 for (int start = 0; start < count; ++start) {
39
40 int differences = 0;
41
42 for (int i = 0; i < count; ++i) {
43
44 int j = (count + start + order * i) % count;
45
46 if (a[i] != b[j]) ++differences;
47 }
48
49 if (!differences) return 0;
50 }
51 }
52 return 1;
53}
54
55void utest_check_grid_file(
56 char const * filename, char const * grid_name,
57 size_t ref_num_cells, size_t ref_num_corners_per_cell,
58 double * ref_cla, double * ref_clo, double * ref_lat, double * ref_lon,
59 int * ref_cell_global_ids,int * ref_core_cell_mask,
60 int * ref_vertex_global_ids,int * ref_core_vertex_mask,
61 int * ref_edge_global_ids,int * ref_core_edge_mask) {
62
63 if (!yac_file_exists(filename)) {
64 PUT_ERR("error file does not exist");
65 return;
66 }
67
68 size_t grid_name_len = strlen(grid_name) + 1;
69
70 char nv_dim_name[3 + grid_name_len];
71 char nc_dim_name[3 + grid_name_len];
72
73 char cla_var_name[4 + grid_name_len];
74 char clo_var_name[4 + grid_name_len];
75 char lat_var_name[4 + grid_name_len];
76 char lon_var_name[4 + grid_name_len];
77 char gid_var_name[4 + grid_name_len];
78 char cmk_var_name[4 + grid_name_len];
79 char vgid_var_name[5 + grid_name_len];
80 char vcmk_var_name[5 + grid_name_len];
81 char egid_var_name[5 + grid_name_len];
82 char ecmk_var_name[5 + grid_name_len];
83
84 snprintf(nv_dim_name, 3 + grid_name_len, "nv_%s", grid_name);
85 snprintf(nc_dim_name, 3 + grid_name_len, "nc_%s", grid_name);
86
87 snprintf(cla_var_name, 4 + grid_name_len, "%s.cla", grid_name);
88 snprintf(clo_var_name, 4 + grid_name_len, "%s.clo", grid_name);
89 snprintf(lat_var_name, 4 + grid_name_len, "%s.lat", grid_name);
90 snprintf(lon_var_name, 4 + grid_name_len, "%s.lon", grid_name);
91 snprintf(gid_var_name, 4 + grid_name_len, "%s.gid", grid_name);
92 snprintf(cmk_var_name, 4 + grid_name_len, "%s.cmk", grid_name);
93 snprintf(vgid_var_name, 5 + grid_name_len, "%s.vgid", grid_name);
94 snprintf(vcmk_var_name, 5 + grid_name_len, "%s.vcmk", grid_name);
95 snprintf(egid_var_name, 5 + grid_name_len, "%s.egid", grid_name);
96 snprintf(ecmk_var_name, 5 + grid_name_len, "%s.ecmk", grid_name);
97
98 int ncid;
99 yac_nc_open(filename, NC_NOWRITE, &ncid);
100
101 int nv_dim_id, nc_dim_id;
102 yac_nc_inq_dimid(ncid, nv_dim_name, &nv_dim_id);
103 yac_nc_inq_dimid(ncid, nc_dim_name, &nc_dim_id);
104
105 int cla_var_id, clo_var_id, lat_var_id, lon_var_id,
106 gid_var_id, cmk_var_id, vgid_var_id, vcmk_var_id,
107 egid_var_id, ecmk_var_id;
108 yac_nc_inq_varid(ncid, cla_var_name, &cla_var_id);
109 yac_nc_inq_varid(ncid, clo_var_name, &clo_var_id);
110 if (ref_lat)
111 yac_nc_inq_varid(ncid, lat_var_name, &lat_var_id);
112 if (ref_lon)
113 yac_nc_inq_varid(ncid, lon_var_name, &lon_var_id);
114 if (ref_cell_global_ids)
115 yac_nc_inq_varid(ncid, gid_var_name, &gid_var_id);
116 if (ref_core_cell_mask)
117 yac_nc_inq_varid(ncid, cmk_var_name, &cmk_var_id);
118 if (ref_vertex_global_ids)
119 yac_nc_inq_varid(ncid, vgid_var_name, &vgid_var_id);
120 if (ref_core_vertex_mask)
121 yac_nc_inq_varid(ncid, vcmk_var_name, &vcmk_var_id);
122 if (ref_edge_global_ids)
123 yac_nc_inq_varid(ncid, egid_var_name, &egid_var_id);
124 if (ref_core_edge_mask)
125 yac_nc_inq_varid(ncid, ecmk_var_name, &ecmk_var_id);
126
127 size_t nv_dim_len, nc_dim_len;
128 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, nv_dim_id, &nv_dim_len));
129 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, nc_dim_id, &nc_dim_len));
130
131 if (nv_dim_len != ref_num_corners_per_cell)
132 PUT_ERR("wrong nv dim size");
133 if (nc_dim_len < ref_num_cells) PUT_ERR("wrong nc dim size");
134
135 nc_type type;
136 int ndims;
137 int dimids[NC_MAX_VAR_DIMS];
138
139 YAC_HANDLE_ERROR(
140 nc_inq_var(ncid, cla_var_id, NULL, &type, &ndims, dimids, NULL));
141 if (type != NC_DOUBLE) PUT_ERR("wrong type for cla");
142 if (ndims != 2) PUT_ERR("wrong ndims for cla");
143 if ((dimids[0] != nc_dim_id) ||
144 (dimids[1] != nv_dim_id)) PUT_ERR("wrong dimensions for cla");
145
146 YAC_HANDLE_ERROR(
147 nc_inq_var(ncid, clo_var_id, NULL, &type, &ndims, dimids, NULL));
148 if (type != NC_DOUBLE) PUT_ERR("wrong type for clo");
149 if (ndims != 2) PUT_ERR("wrong ndims for clo");
150 if ((dimids[0] != nc_dim_id) ||
151 (dimids[1] != nv_dim_id)) PUT_ERR("wrong dimensions for clo");
152
153 if (ref_lat) {
154 YAC_HANDLE_ERROR(
155 nc_inq_var(ncid, lat_var_id, NULL, &type, &ndims, dimids, NULL));
156 if (type != NC_DOUBLE) PUT_ERR("wrong type for lat");
157 if (ndims != 1) PUT_ERR("wrong ndims for lat");
158 if (dimids[0] != nc_dim_id) PUT_ERR("wrong dimensions for lat");
159 }
160
161 if (ref_lon) {
162 YAC_HANDLE_ERROR(
163 nc_inq_var(ncid, lon_var_id, NULL, &type, &ndims, dimids, NULL));
164 if (type != NC_DOUBLE) PUT_ERR("wrong type for lon");
165 if (ndims != 1) PUT_ERR("wrong ndims for lon");
166 if (dimids[0] != nc_dim_id) PUT_ERR("wrong dimensions for lon");
167 }
168
169 if (ref_cell_global_ids) {
170 YAC_HANDLE_ERROR(
171 nc_inq_var(ncid, gid_var_id, NULL, &type, &ndims, dimids, NULL));
172 if (type != NC_INT) PUT_ERR("wrong type for gid");
173 if (ndims != 1) PUT_ERR("wrong ndims for gid");
174 if (dimids[0] != nc_dim_id) PUT_ERR("wrong dimensions for gid");
175 }
176
177 if (ref_core_cell_mask) {
178 YAC_HANDLE_ERROR(
179 nc_inq_var(ncid, cmk_var_id, NULL, &type, &ndims, dimids, NULL));
180 if (type != NC_INT) PUT_ERR("wrong type for cmk");
181 if (ndims != 1) PUT_ERR("wrong ndims for cmk");
182 if (dimids[0] != nc_dim_id) PUT_ERR("wrong dimensions for cmk");
183 }
184
185 if (ref_vertex_global_ids) {
186 YAC_HANDLE_ERROR(
187 nc_inq_var(ncid, vgid_var_id, NULL, &type, &ndims, dimids, NULL));
188 if (type != NC_INT) PUT_ERR("wrong type for vgid");
189 if (ndims != 2) PUT_ERR("wrong ndims for vgid");
190 if ((dimids[0] != nc_dim_id) ||
191 (dimids[1] != nv_dim_id)) PUT_ERR("wrong dimensions for vgid");
192 }
193
194 if (ref_core_vertex_mask) {
195 YAC_HANDLE_ERROR(
196 nc_inq_var(ncid, vcmk_var_id, NULL, &type, &ndims, dimids, NULL));
197 if (type != NC_INT) PUT_ERR("wrong type for vcmk");
198 if (ndims != 2) PUT_ERR("wrong ndims for vcmk");
199 if ((dimids[0] != nc_dim_id) ||
200 (dimids[1] != nv_dim_id)) PUT_ERR("wrong dimensions for vcmk");
201 }
202
203 if (ref_edge_global_ids) {
204 YAC_HANDLE_ERROR(
205 nc_inq_var(ncid, egid_var_id, NULL, &type, &ndims, dimids, NULL));
206 if (type != NC_INT) PUT_ERR("wrong type for egid");
207 if (ndims != 2) PUT_ERR("wrong ndims for egid");
208 if ((dimids[0] != nc_dim_id) ||
209 (dimids[1] != nv_dim_id)) PUT_ERR("wrong dimensions for egid");
210 }
211
212 if (ref_core_edge_mask) {
213 YAC_HANDLE_ERROR(
214 nc_inq_var(ncid, ecmk_var_id, NULL, &type, &ndims, dimids, NULL));
215 if (type != NC_INT) PUT_ERR("wrong type for ecmk");
216 if (ndims != 2) PUT_ERR("wrong ndims for ecmk");
217 if ((dimids[0] != nc_dim_id) ||
218 (dimids[1] != nv_dim_id)) PUT_ERR("wrong dimensions for ecmk");
219 }
220
221 double * clo = malloc(nc_dim_len * nv_dim_len * sizeof(*clo));
222 double * cla = malloc(nc_dim_len * nv_dim_len * sizeof(*cla));
223 double * lon = ref_lon?malloc(nc_dim_len * sizeof(*lon)):NULL;
224 double * lat = ref_lat?malloc(nc_dim_len * sizeof(*lat)):NULL;
225 int * gid = ref_cell_global_ids?malloc(nc_dim_len * sizeof(*gid)):NULL;
226 int * cmk = ref_core_cell_mask?malloc(nc_dim_len * sizeof(*cmk)):NULL;
227 int * vgid = ref_vertex_global_ids?malloc(nc_dim_len * nv_dim_len * sizeof(*vgid)):NULL;
228 int * vcmk = ref_core_vertex_mask?malloc(nc_dim_len * nv_dim_len * sizeof(*vcmk)):NULL;
229 int * egid = ref_edge_global_ids?malloc(nc_dim_len * nv_dim_len * sizeof(*egid)):NULL;
230 int * ecmk = ref_core_edge_mask?malloc(nc_dim_len * nv_dim_len * sizeof(*ecmk)):NULL;
231
232 YAC_HANDLE_ERROR(nc_get_var_double(ncid, cla_var_id, cla));
233 YAC_HANDLE_ERROR(nc_get_var_double(ncid, clo_var_id, clo));
234 if (ref_lat)
235 YAC_HANDLE_ERROR(nc_get_var_double(ncid, lat_var_id, lat));
236 if (ref_lon)
237 YAC_HANDLE_ERROR(nc_get_var_double(ncid, lon_var_id, lon));
238 if (ref_cell_global_ids)
239 YAC_HANDLE_ERROR(nc_get_var_int(ncid, gid_var_id, gid));
240 if (ref_core_cell_mask)
241 YAC_HANDLE_ERROR(nc_get_var_int(ncid, cmk_var_id, cmk));
242 if (ref_vertex_global_ids)
243 YAC_HANDLE_ERROR(nc_get_var_int(ncid, vgid_var_id, vgid));
244 if (ref_core_vertex_mask)
245 YAC_HANDLE_ERROR(nc_get_var_int(ncid, vcmk_var_id, vcmk));
246 if (ref_edge_global_ids)
247 YAC_HANDLE_ERROR(nc_get_var_int(ncid, egid_var_id, egid));
248 if (ref_core_edge_mask)
249 YAC_HANDLE_ERROR(nc_get_var_int(ncid, ecmk_var_id, ecmk));
250
251 int * ref_cell_flag = calloc(ref_num_cells, sizeof(*ref_cell_flag));
252
253 // for all cells in the grid file
254 for (size_t i = 0; i < nc_dim_len; ++i) {
255
256 // get the coordinates of the current cell
257 double * curr_cla = cla + i * nv_dim_len;
258 double * curr_clo = clo + i * nv_dim_len;
259
260 // match the current cell with the reference cells
261 size_t match_idx = SIZE_MAX;
262 {
263 if (ref_cell_global_ids) {
264 for (size_t j = 0; (j < ref_num_cells) && (match_idx == SIZE_MAX); ++j)
265 if (gid[i] == ref_cell_global_ids[j])
266 match_idx = j;
267 } else if (ref_lat && ref_lon) {
268 for (size_t j = 0; (j < ref_num_cells) && (match_idx == SIZE_MAX); ++j)
269 if ((fabs(lat[i] - ref_lat[j]) < 1e-3) &&
270 (fabs(lon[i] - ref_lon[j]) < 1e-3))
271 match_idx = j;
272 } else {
273 for (size_t j = 0; (j < ref_num_cells) && (match_idx == SIZE_MAX); ++j)
274 if (!utest_compare_cell_coords(
275 curr_cla, ref_cla + j * ref_num_corners_per_cell,
276 ref_num_corners_per_cell) &&
277 !utest_compare_cell_coords(
278 curr_clo, ref_clo + j * ref_num_corners_per_cell,
279 ref_num_corners_per_cell))
280 match_idx = j;
281 }
282 }
283
284 if (match_idx == SIZE_MAX) {
285
286 PUT_ERR("error no matching cell");
287
288 } else {
289
290 ref_cell_flag[match_idx] = 1;
291
292 if (utest_compare_cell_coords(
293 curr_cla, ref_cla + match_idx * ref_num_corners_per_cell,
294 ref_num_corners_per_cell))
295 PUT_ERR("wrong cla");
296 if (utest_compare_cell_coords(
297 curr_clo, ref_clo + match_idx * ref_num_corners_per_cell,
298 ref_num_corners_per_cell))
299 PUT_ERR("wrong clo");
300 if (ref_lat)
301 if (fabs(lat[i] - ref_lat[match_idx]) > 1e-3) PUT_ERR("wrong lat");
302 if (ref_lon)
303 if (fabs(lon[i] - ref_lon[match_idx]) > 1e-3) PUT_ERR("wrong lon");
304 if (ref_cell_global_ids)
305 if (gid[i] != (int)(ref_cell_global_ids[match_idx])) PUT_ERR("wrong gid");
306 if (ref_core_cell_mask)
307 if (cmk[i] != (ref_core_cell_mask[match_idx]))
308 PUT_ERR("wrong cmk");
309 if (ref_vertex_global_ids)
310 if (utest_compare_ints(
311 vgid + i * nv_dim_len,
312 ref_vertex_global_ids + match_idx * ref_num_corners_per_cell,
313 ref_num_corners_per_cell))
314 PUT_ERR("wrong vgid");
315 if (ref_core_vertex_mask)
316 if (utest_compare_ints(
317 vcmk + i * nv_dim_len,
318 ref_core_vertex_mask + match_idx * ref_num_corners_per_cell,
319 ref_num_corners_per_cell))
320 PUT_ERR("wrong vcmk");
321 if (ref_edge_global_ids)
322 if (utest_compare_ints(
323 egid + i * nv_dim_len,
324 ref_edge_global_ids + match_idx * ref_num_corners_per_cell,
325 ref_num_corners_per_cell))
326 PUT_ERR("wrong egid");
327 if (ref_core_edge_mask)
328 if (utest_compare_ints(
329 ecmk + i * nv_dim_len,
330 ref_core_edge_mask + match_idx * ref_num_corners_per_cell,
331 ref_num_corners_per_cell))
332 PUT_ERR("wrong ecmk");
333 }
334 }
335
336 size_t match_count = 0;
337 for (size_t i = 0; i < ref_num_cells; ++i)
338 if (ref_cell_flag[i]) match_count++;
339
340 if (match_count != ref_num_cells) PUT_ERR("missing cells");
341
342 free(ref_cell_flag);
343 free(ecmk);
344 free(egid);
345 free(vcmk);
346 free(vgid);
347 free(cmk);
348 free(gid);
349 free(lat);
350 free(lon);
351 free(clo);
352 free(cla);
353
354 YAC_HANDLE_ERROR(nc_close(ncid));
355}
356
357void utest_write_dummy_exodus_grid_file(
358 char const * grid_filename, size_t num_lon, size_t num_lat,
359 double lon_range[2], double lat_range[2]) {
360
361 size_t num_nodes = num_lon * num_lat;
362 size_t num_elem = (num_lon - 1) * (num_lat - 1);
363
364 { // grid file
365
366 int ncid;
367
368 // create file
369 yac_nc_create(grid_filename, NC_CLOBBER, &ncid);
370
371 int coord_dim_id;
372 int node_dim_id;
373 int elem_dim_id;
374 int num_el_blk_dim_id;
375 int elem_blk1_dim_id;
376 int node_per_el1_dim_id;
377
378 // define dimensions
379 YAC_HANDLE_ERROR(
380 nc_def_dim(ncid, "num_dim", 3, &coord_dim_id));
381 YAC_HANDLE_ERROR(
382 nc_def_dim(ncid, "num_nodes", num_nodes, &node_dim_id));
383 YAC_HANDLE_ERROR(
384 nc_def_dim(ncid, "num_elem", num_elem, &elem_dim_id));
385 YAC_HANDLE_ERROR(
386 nc_def_dim(ncid, "num_el_blk", 1, &num_el_blk_dim_id));
387 YAC_HANDLE_ERROR(
388 nc_def_dim(ncid, "num_el_in_blk1", num_elem, &elem_blk1_dim_id));
389 YAC_HANDLE_ERROR(
390 nc_def_dim(ncid, "num_nod_per_el1", 4, &node_per_el1_dim_id));
391
392 int elem_dim_ids[2] = {elem_blk1_dim_id, node_per_el1_dim_id};
393 int coord_dim_ids[2] = {coord_dim_id, node_dim_id};
394
395 int var_elem_id;
396 int var_coord_id;
397
398 // define variable
399 YAC_HANDLE_ERROR(
400 nc_def_var(ncid, "connect1", NC_INT, 2, elem_dim_ids, &var_elem_id));
401 YAC_HANDLE_ERROR(
402 nc_def_var(ncid, "coord", NC_DOUBLE, 2, coord_dim_ids, &var_coord_id));
403
404 // end definition
405 YAC_HANDLE_ERROR(nc_enddef(ncid));
406
407 // write grid data
408 {
409 int elem[num_elem][4];
410
411 for (size_t i = 0, k = 0; i < (num_lat - 1); ++i) {
412 for (size_t j = 0; j < (num_lon - 1); ++j, ++k) {
413 elem[k][0] = (int)((i + 0) * num_lon + (j + 0) + 1);
414 elem[k][1] = (int)((i + 0) * num_lon + (j + 1) + 1);
415 elem[k][2] = (int)((i + 1) * num_lon + (j + 1) + 1);
416 elem[k][3] = (int)((i + 1) * num_lon + (j + 0) + 1);
417 }
418 }
419 YAC_HANDLE_ERROR(nc_put_var_int(ncid, var_elem_id, &elem[0][0]));
420 }
421 {
422 double coords[3][num_nodes];
423
424 for (size_t i = 0, k = 0; i < num_lat; ++i) {
425 double lat =
426 ((double)i * (lat_range[1] - lat_range[0])) / (double)(num_lat - 1) +
427 lat_range[0];
428 for (size_t j = 0; j < num_lon; ++j, ++k) {
429 double lon =
430 ((double)j * (lon_range[1] - lon_range[0])) / (double)(num_lon - 1) +
431 lon_range[0];
432 double temp_coord[3];
433 LLtoXYZ_deg(lon, lat, temp_coord);
434 coords[0][k] = temp_coord[0];
435 coords[1][k] = temp_coord[1];
436 coords[2][k] = temp_coord[2];
437 }
438 }
439 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_coord_id, &coords[0][0]));
440 }
441
442 YAC_HANDLE_ERROR(nc_close(ncid));
443 }
444}
445
446void utest_write_dummy_scrip_grid_file(
447 char * grid_name, char * grid_filename, char * mask_filename,
448 int with_corners, size_t num_lon, size_t num_lat,
449 double lon_range[2], double lat_range[2]) {
450
451 write_dummy_scrip_grid_file_ext(
452 grid_name, grid_filename, mask_filename, with_corners?4:0,
453 num_lon, num_lat, lon_range, lat_range, 0);
454}
455
456void write_dummy_scrip_grid_file_ext(
457 char * grid_name, char * grid_filename, char * mask_filename,
458 size_t num_corners, size_t num_lon, size_t num_lat,
459 double lon_range[2], double lat_range[2], size_t invalid_cell_stride) {
460
461 int with_corners = num_corners > 0;
462
463 { // grid file
464 int ncid;
465
466 // create file
467 yac_nc_create(grid_filename, NC_CLOBBER, &ncid);
468
469 char crn_dim_name[128];
470 char x_dim_name[128];
471 char y_dim_name[128];
472
473 sprintf(crn_dim_name, "crn_%s", grid_name);
474 sprintf(x_dim_name, "x_%s", grid_name);
475 sprintf(y_dim_name, "y_%s", grid_name);
476
477 int dim_crn_id = -1;
478 int dim_x_id;
479 int dim_y_id;
480
481 // define dimensions
482 if (with_corners)
483 YAC_HANDLE_ERROR(nc_def_dim(ncid, crn_dim_name, num_corners, &dim_crn_id));
484 YAC_HANDLE_ERROR(nc_def_dim(ncid, x_dim_name, num_lon, &dim_x_id));
485 YAC_HANDLE_ERROR(nc_def_dim(ncid, y_dim_name, num_lat, &dim_y_id));
486
487 char cla_var_name[128];
488 char clo_var_name[128];
489 char lat_var_name[128];
490 char lon_var_name[128];
491
492 sprintf(cla_var_name, "%s.cla", grid_name);
493 sprintf(clo_var_name, "%s.clo", grid_name);
494 sprintf(lat_var_name, "%s.lat", grid_name);
495 sprintf(lon_var_name, "%s.lon", grid_name);
496
497 int corner_dim_ids[3] = {dim_crn_id, dim_y_id, dim_x_id};
498 int cell_dim_ids[2] = {dim_y_id, dim_x_id};
499
500 int var_cla_id = -1;
501 int var_clo_id = -1;
502 int var_lat_id;
503 int var_lon_id;
504
505 char degree[] = "degree";
506 char title[] = "This is a reg lon-lat dummy grid";
507
508 // define variable
509 if (with_corners) {
510 YAC_HANDLE_ERROR(
511 nc_def_var(
512 ncid, cla_var_name, NC_DOUBLE, 3, corner_dim_ids, &var_cla_id));
513 YAC_HANDLE_ERROR(
514 nc_put_att_text(ncid, var_cla_id, "units", strlen(degree), degree));
515 YAC_HANDLE_ERROR(
516 nc_put_att_text(ncid, var_cla_id, "title", strlen(title), title));
517
518 YAC_HANDLE_ERROR(
519 nc_def_var(
520 ncid, clo_var_name, NC_DOUBLE, 3, corner_dim_ids, &var_clo_id));
521 YAC_HANDLE_ERROR(
522 nc_put_att_text(ncid, var_clo_id, "units", strlen(degree), degree));
523 YAC_HANDLE_ERROR(
524 nc_put_att_text(ncid, var_clo_id, "title", strlen(title), title));
525 }
526
527 YAC_HANDLE_ERROR(
528 nc_def_var(
529 ncid, lat_var_name, NC_DOUBLE, 2, cell_dim_ids, &var_lat_id));
530 YAC_HANDLE_ERROR(
531 nc_put_att_text(ncid, var_lat_id, "units", strlen(degree), degree));
532 YAC_HANDLE_ERROR(
533 nc_put_att_text(ncid, var_lat_id, "title", strlen(title), title));
534
535 YAC_HANDLE_ERROR(
536 nc_def_var(
537 ncid, lon_var_name, NC_DOUBLE, 2, cell_dim_ids, &var_lon_id));
538 YAC_HANDLE_ERROR(
539 nc_put_att_text(ncid, var_lon_id, "units", strlen(degree), degree));
540 YAC_HANDLE_ERROR(
541 nc_put_att_text(ncid, var_lon_id, "title", strlen(title), title));
542
543
544 // end definition
545 YAC_HANDLE_ERROR(nc_enddef(ncid));
546
547 // write grid data
548
549 double cla[num_corners?num_corners:1][num_lat][num_lon];
550 double clo[num_corners?num_corners:1][num_lat][num_lon];
551 double lat[num_lat][num_lon];
552 double lon[num_lat][num_lon];
553
554 for (size_t i = 0; i < num_lon; ++i) {
555 double vertex_lon[2] =
556 {((lon_range[1] - lon_range[0]) * (double)i)/(double)num_lon,
557 ((lon_range[1] - lon_range[0]) * (double)(i+1))/(double)num_lon};
558 for (size_t j = 0; j < num_lat; ++j) {
559 double vertex_lat[2] =
560 {((lat_range[1] - lat_range[0]) * (double)j)/(double)num_lat,
561 ((lat_range[1] - lat_range[0]) * (double)(j+1))/(double)num_lat};
562 if (with_corners) {
563 cla[0][j][i] = lat_range[0] + vertex_lat[0];
564 cla[1][j][i] = lat_range[0] + vertex_lat[0];
565 cla[2][j][i] = lat_range[0] + vertex_lat[1];
566 cla[3][j][i] = lat_range[0] + vertex_lat[1];
567 clo[0][j][i] = lon_range[0] + vertex_lon[0];
568 clo[1][j][i] = lon_range[0] + vertex_lon[1];
569 clo[2][j][i] = lon_range[0] + vertex_lon[1];
570 clo[3][j][i] = lon_range[0] + vertex_lon[0];
571 // split the remaining corners evenly along the edge that closes
572 // the cell (from corner 3 back to corner 0)
573 for (size_t n = 4; n < num_corners; ++n) {
574 cla[n][j][i] =
575 lat_range[0] + vertex_lat[1] +
576 ((vertex_lat[0] - vertex_lat[1]) * (double)(n - 3)) /
577 (double)(num_corners - 3);
578 clo[n][j][i] = lon_range[0] + vertex_lon[0];
579 }
580 }
581 lat[j][i] = lat_range[0] + (vertex_lat[0] + vertex_lat[1]) * 0.5;
582 lon[j][i] = lon_range[0] + (vertex_lon[0] + vertex_lon[1]) * 0.5;
583 }
584 }
585
586
587 if (with_corners) {
588 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_cla_id, &cla[0][0][0]));
589 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_clo_id, &clo[0][0][0]));
590 }
591 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_lat_id, &lat[0][0]));
592 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_lon_id, &lon[0][0]));
593
594 YAC_HANDLE_ERROR(nc_close(ncid));
595 }
596
597 { // mask file
598 int ncid;
599
600 // create file
601 yac_nc_create(mask_filename, NC_CLOBBER, &ncid);
602
603 char x_dim_name[128];
604 char y_dim_name[128];
605
606 sprintf(x_dim_name, "x_%s", grid_name);
607 sprintf(y_dim_name, "y_%s", grid_name);
608
609 int dim_x_id;
610 int dim_y_id;
611
612 // define dimensions
613 YAC_HANDLE_ERROR(nc_def_dim(ncid, x_dim_name, num_lon, &dim_x_id));
614 YAC_HANDLE_ERROR(nc_def_dim(ncid, y_dim_name, num_lat, &dim_y_id));
615
616 char frc_var_name[128];
617 char msk_var_name[128];
618
619 sprintf(frc_var_name, "%s.frc", grid_name);
620 sprintf(msk_var_name, "%s.msk", grid_name);
621
622 int dim_ids[2] = {dim_y_id, dim_x_id};
623
624 int var_frc_id;
625 int var_msk_id;
626
627 char adim[] = "adim";
628
629 // define variable
630 YAC_HANDLE_ERROR(
631 nc_def_var(
632 ncid, frc_var_name, NC_DOUBLE, 2, dim_ids, &var_frc_id));
633 YAC_HANDLE_ERROR(
634 nc_put_att_text(ncid, var_frc_id, "units", strlen(adim), adim));
635
636 YAC_HANDLE_ERROR(
637 nc_def_var(
638 ncid, msk_var_name, NC_INT, 2, dim_ids, &var_msk_id));
639 YAC_HANDLE_ERROR(
640 nc_put_att_text(ncid, var_msk_id, "units", strlen(adim), adim));
641
642
643 // end definition
644 YAC_HANDLE_ERROR(nc_enddef(ncid));
645
646 // write grid data
647
648 double frc[num_lat][num_lon];
649 int msk[num_lat][num_lon];
650
651 for (size_t i = 0; i < num_lon; ++i) {
652 for (size_t j = 0; j < num_lat; ++j) {
653 frc[j][i] = 1;
654 // mask out every invalid_cell_stride-th cell
655 msk[j][i] =
656 (invalid_cell_stride > 0) &&
657 (((i + j * num_lon) % invalid_cell_stride) == 0);
658 }
659 }
660
661 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_frc_id, &frc[0][0]));
662 YAC_HANDLE_ERROR(nc_put_var_int(ncid, var_msk_id, &msk[0][0]));
663
664 YAC_HANDLE_ERROR(nc_close(ncid));
665 }
666}
667
668void write_dummy_mpiom_grid_file(
669 char * grid_filename, size_t num_x, size_t num_y,
670 double lon_range[2], double lat_range[2], int mask_ndims) {
671
672 int ncid;
673
674 yac_nc_create(grid_filename, NC_CLOBBER, &ncid);
675
676 // define dimensions
677 int dim_x_id, dim_y_id, dim_nv_id, dim_time_id, dim_depth_id;
678 YAC_HANDLE_ERROR(nc_def_dim(ncid, "x", num_x, &dim_x_id));
679 YAC_HANDLE_ERROR(nc_def_dim(ncid, "y", num_y, &dim_y_id));
680 YAC_HANDLE_ERROR(nc_def_dim(ncid, "nv", 4, &dim_nv_id));
681 YAC_HANDLE_ERROR(nc_def_dim(ncid, "time", 1, &dim_time_id));
682 YAC_HANDLE_ERROR(nc_def_dim(ncid, "depth", 1, &dim_depth_id));
683
684 int cell_dim_ids[2] = {dim_y_id, dim_x_id};
685 int bounds_dim_ids[3] = {dim_y_id, dim_x_id, dim_nv_id};
686 int mask_dim_ids[4] =
687 {dim_time_id,
688 (mask_ndims == 4)?dim_depth_id:dim_y_id,
689 (mask_ndims == 4)?dim_y_id:dim_x_id,
690 dim_x_id};
691
692 // define variables
693 int var_lon_id, var_lat_id, var_lon_bounds_id, var_lat_bounds_id, var_msk_id;
694 YAC_HANDLE_ERROR(
695 nc_def_var(ncid, "lon", NC_DOUBLE, 2, cell_dim_ids, &var_lon_id));
696 YAC_HANDLE_ERROR(
697 nc_def_var(ncid, "lat", NC_DOUBLE, 2, cell_dim_ids, &var_lat_id));
698 YAC_HANDLE_ERROR(
699 nc_def_var(
700 ncid, "lon_bounds", NC_DOUBLE, 3, bounds_dim_ids, &var_lon_bounds_id));
701 YAC_HANDLE_ERROR(
702 nc_def_var(
703 ncid, "lat_bounds", NC_DOUBLE, 3, bounds_dim_ids, &var_lat_bounds_id));
704 YAC_HANDLE_ERROR(
705 nc_def_var(ncid, "var1", NC_INT, mask_ndims, mask_dim_ids, &var_msk_id));
706
707 YAC_HANDLE_ERROR(nc_enddef(ncid));
708
709 // write grid data
710 double lon[num_y][num_x];
711 double lat[num_y][num_x];
712 double lon_bounds[num_y][num_x][4];
713 double lat_bounds[num_y][num_x][4];
714 int msk[num_y][num_x];
715
716 for (size_t j = 0; j < num_y; ++j) {
717 double vertex_lat[2] =
718 {lat_range[0] +
719 ((lat_range[1] - lat_range[0]) * (double)j) / (double)num_y,
720 lat_range[0] +
721 ((lat_range[1] - lat_range[0]) * (double)(j+1)) / (double)num_y};
722 for (size_t i = 0; i < num_x; ++i) {
723 double vertex_lon[2] =
724 {lon_range[0] +
725 ((lon_range[1] - lon_range[0]) * (double)i) / (double)num_x,
726 lon_range[0] +
727 ((lon_range[1] - lon_range[0]) * (double)(i+1)) / (double)num_x};
728 lon_bounds[j][i][0] = vertex_lon[0];
729 lon_bounds[j][i][1] = vertex_lon[1];
730 lon_bounds[j][i][2] = vertex_lon[1];
731 lon_bounds[j][i][3] = vertex_lon[0];
732 lat_bounds[j][i][0] = vertex_lat[0];
733 lat_bounds[j][i][1] = vertex_lat[0];
734 lat_bounds[j][i][2] = vertex_lat[1];
735 lat_bounds[j][i][3] = vertex_lat[1];
736 lon[j][i] = (vertex_lon[0] + vertex_lon[1]) * 0.5;
737 lat[j][i] = (vertex_lat[0] + vertex_lat[1]) * 0.5;
738 // a three dimensional mask is used directly, while a four dimensional
739 // one marks valid cells with the value one
740 msk[j][i] = (int)((i + j) & 1);
741 }
742 }
743
744 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_lon_id, &lon[0][0]));
745 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_lat_id, &lat[0][0]));
746 YAC_HANDLE_ERROR(
747 nc_put_var_double(ncid, var_lon_bounds_id, &lon_bounds[0][0][0]));
748 YAC_HANDLE_ERROR(
749 nc_put_var_double(ncid, var_lat_bounds_id, &lat_bounds[0][0][0]));
750 YAC_HANDLE_ERROR(nc_put_var_int(ncid, var_msk_id, &msk[0][0]));
751
752 YAC_HANDLE_ERROR(nc_close(ncid));
753}
LLtoXYZ_deg
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:278
write_dummy_mpiom_grid_file
void write_dummy_mpiom_grid_file(char *grid_filename, size_t num_x, size_t num_y, double lon_range[2], double lat_range[2], int mask_ndims)
Definition grid_file_common.c:668
write_dummy_scrip_grid_file_ext
void write_dummy_scrip_grid_file_ext(char *grid_name, char *grid_filename, char *mask_filename, size_t num_corners, size_t num_lon, size_t num_lat, double lon_range[2], double lat_range[2], size_t invalid_cell_stride)
Definition grid_file_common.c:456
grid_file_common.h
type
enum callback_type type
Definition interp_method_check.c:142
yac_nc_create
void yac_nc_create(const char *path, int cmode, int *ncidp)
Definition io_utils.c:326
yac_nc_inq_varid
void yac_nc_inq_varid(int ncid, char const *name, int *varidp)
Definition io_utils.c:369
yac_file_exists
int yac_file_exists(const char *filename)
Check whether a file exists.
Definition io_utils.c:394
yac_nc_open
void yac_nc_open(const char *path, int omode, int *ncidp)
Definition io_utils.c:311
yac_nc_inq_dimid
void yac_nc_inq_dimid(int ncid, char const *name, int *dimidp)
Definition io_utils.c:344
io_utils.h
test_core.coords
coords
Definition test_core.py:57
test_raw_exchange.i
int i
Definition test_raw_exchange.py:89
PUT_ERR
#define PUT_ERR(string)
Definition tests.h:10
YAC_HANDLE_ERROR
#define YAC_HANDLE_ERROR(exp)
Definition toy_output.c:13
utils_common.h