les2htcp.c (27096B)
1 /* Copyright (C) 2018, 2020-2023, 2025, 2025 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2018 Centre National de la Recherche Scientifique 3 * Copyright (C) 2018 Université Paul Sabatier 4 * 5 * This program is free software: you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation, either version 3 of the License, or 8 * (at your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 * GNU General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 17 18 #define _POSIX_C_SOURCE 200112L /* getopt/close support */ 19 20 #include "les2htcp.h" 21 22 #include <rsys/cstr.h> 23 #include <rsys/math.h> 24 #include <rsys/mem_allocator.h> 25 #include <rsys/rsys.h> 26 27 #include <alloca.h> 28 #include <errno.h> 29 #include <netcdf.h> 30 #include <string.h> 31 #include <unistd.h> /* getopt, sysconf */ 32 #include <fcntl.h> /* open */ 33 #include <sys/stat.h> /* S_IRUSR & S_IWUSR */ 34 35 /******************************************************************************* 36 * Command arguments 37 ******************************************************************************/ 38 struct args { 39 const char* output; 40 const char* input; 41 double fp_to_meter; 42 long pagesize; 43 int force_overwrite; 44 int check; 45 int no_output; 46 int quit; /* Quit the application */ 47 }; 48 #define ARGS_DEFAULT__ {NULL,NULL,1.0,-1,0,0,0,0} 49 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__; 50 51 static INLINE void 52 usage(const char* cmd) 53 { 54 ASSERT(cmd); 55 printf("usage: %s [-cfhqv] [-m float_to_meter] [-o output] [-p pagesize] " 56 "-i netcdf\n", cmd); 57 } 58 59 static void 60 args_release(struct args* args) 61 { 62 ASSERT(args); 63 *args = ARGS_DEFAULT; 64 } 65 66 static res_T 67 args_init(struct args* args, const int argc, char** argv) 68 { 69 long system_pagesize; 70 int opt; 71 res_T res = RES_OK; 72 ASSERT(args && argc && argv); 73 74 system_pagesize = sysconf(_SC_PAGESIZE); 75 if(system_pagesize == -1) { 76 fprintf(stderr, 77 "Error when querying the size of a system page -- %s\n", 78 strerror(errno)); 79 res = RES_BAD_OP; 80 goto error; 81 } 82 83 while((opt = getopt(argc, argv, "cfhi:m:o:p:qv")) != -1) { 84 switch(opt) { 85 case 'c': args->check = 1; break; 86 case 'f': args->force_overwrite = 1; break; 87 case 'h': 88 usage(argv[0]); 89 args_release(args); 90 args->quit = 1; 91 goto exit; 92 case 'i': args->input = optarg; break; 93 case 'm': 94 res = cstr_to_double(optarg, &args->fp_to_meter); 95 if(res == RES_OK && args->fp_to_meter <= 0) res = RES_BAD_ARG; 96 break; 97 case 'o': args->output = optarg; break; 98 case 'p': 99 res = cstr_to_long(optarg, &args->pagesize); 100 if(res == RES_OK && !IS_POW2(args->pagesize)) res = RES_BAD_ARG; 101 if(res == RES_OK && args->pagesize < system_pagesize) res = RES_BAD_ARG; 102 break; 103 case 'q': args->no_output = 1; break; 104 case 'v': 105 printf("%s %d.%d.%d\n", 106 argv[0], 107 LES2HTCP_VERSION_MAJOR, 108 LES2HTCP_VERSION_MINOR, 109 LES2HTCP_VERSION_PATCH); 110 args->quit = 1; 111 goto exit; 112 default: res = RES_BAD_ARG; break; 113 } 114 if(res != RES_OK) { 115 if(optarg) { 116 fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n", 117 argv[0], optarg, opt); 118 } 119 goto error; 120 } 121 } 122 123 if(!args->input) { 124 fprintf(stderr, "Missing input file.\n"); 125 res = RES_BAD_ARG; 126 goto error; 127 } 128 129 if(args->no_output) args->output = NULL; 130 131 /* Use the default page size if not explicitly defined by the caller. */ 132 if(args->pagesize < system_pagesize) { 133 args->pagesize = system_pagesize; 134 } 135 136 exit: 137 return res; 138 error: 139 usage(argv[0]); 140 args_release(args); 141 goto exit; 142 } 143 144 /******************************************************************************* 145 * Grid header 146 ******************************************************************************/ 147 struct grid { 148 int32_t nx, ny, nz, ntimes; /* Definition */ 149 int8_t is_z_irregular; 150 double lower[3]; /* Lower position of the grid */ 151 double vxsz_x; /* Size of the voxel is X */ 152 double vxsz_y; /* Size of the voxel in Y */ 153 double* vxsz_z; /* Size of the voxel in Z */ 154 }; 155 156 #define GRID_NULL__ {0,0,0,0,0,{0,0,0},0,0,NULL} 157 static const struct grid GRID_NULL = GRID_NULL__; 158 159 static void 160 grid_release(struct grid* grid) 161 { 162 CHK(grid); 163 if(grid->vxsz_z) mem_rm(grid->vxsz_z); 164 } 165 166 /******************************************************************************* 167 * NetCDF helper functions and macros 168 ******************************************************************************/ 169 #define INVALID_ID -1 170 171 #define NC(Func) { \ 172 const int err__ = nc_ ## Func; \ 173 if(err__ != NC_NOERR) { \ 174 fprintf(stderr, "error:%i:%s\n", __LINE__, ncerr_to_str(err__)); \ 175 abort(); \ 176 } \ 177 } (void)0 178 179 #define NDIMS_MAX 4 180 181 static INLINE const char* 182 ncerr_to_str(const int err) 183 { 184 const char* str = "NC_ERR_<UNKNOWN>"; 185 switch(err) { 186 case NC_EBADGRPID: str = "NC_EBADGRPID"; break; 187 case NC_EBADID: str = "NC_EBADID"; break; 188 case NC_EBADNAME: str = "NC_EBADNAME"; break; 189 case NC_ECHAR: str = "NC_ECHAR"; break; 190 case NC_EDIMMETA: str = "NC_EDIMMETA"; break; 191 case NC_EHDFERR: str = "NC_EHDFERR"; break; 192 case NC_EINVAL: str = "NC_EINVAL"; break; 193 case NC_ENOMEM: str = "NC_ENOMEM"; break; 194 case NC_ENOTATT: str = "NC_ENOTATT"; break; 195 case NC_ENOTVAR: str = "NC_ENOTVAR"; break; 196 case NC_ERANGE: str = "NC_ERANGE"; break; 197 case NC_NOERR: str = "NC_NOERR"; break; 198 } 199 return str; 200 } 201 202 static INLINE const char* 203 nctype_to_str(const nc_type type) 204 { 205 const char* str = "NC_TYPE_<UNKNOWN>"; 206 switch(type) { 207 case NC_NAT: str = "NC_NAT"; break; 208 case NC_BYTE: str = "NC_BYTE"; break; 209 case NC_CHAR: str = "NC_CHAR"; break; 210 case NC_SHORT: str = "NC_SHORT"; break; 211 case NC_LONG: str = "NC_LONG"; break; 212 case NC_FLOAT: str = "NC_FLOAT"; break; 213 case NC_DOUBLE: str = "NC_DOUBLE"; break; 214 case NC_UBYTE: str = "NC_UBYTE"; break; 215 case NC_USHORT: str = "NC_USHORT"; break; 216 case NC_UINT: str = "NC_UINT"; break; 217 case NC_INT64: str = "NC_INT64"; break; 218 case NC_UINT64: str = "NC_UINT64"; break; 219 case NC_STRING: str = "NC_STRING"; break; 220 default: FATAL("Unreachable code.\n"); break; 221 } 222 return str; 223 } 224 225 static INLINE size_t 226 sizeof_nctype(const nc_type type) 227 { 228 size_t sz; 229 switch(type) { 230 case NC_BYTE: 231 case NC_CHAR: 232 case NC_UBYTE: 233 sz = 1; break; 234 case NC_SHORT: 235 case NC_USHORT: 236 sz = 2; break; 237 case NC_FLOAT: 238 case NC_INT: 239 case NC_UINT: 240 sz = 4; break; 241 case NC_DOUBLE: 242 case NC_INT64: 243 case NC_UINT64: 244 sz = 8; break; 245 default: FATAL("Unreachable cde.\n"); break; 246 } 247 return sz; 248 } 249 250 struct var { 251 size_t dim_len[NDIMS_MAX]; /* Length of the dimensions */ 252 int ndims; /* #dimensions */ 253 int id; /* NetCDF identifier */ 254 int type; /* Variable type */ 255 }; 256 static const struct var VAR_NULL = {{0}, 0, 0, 0}; 257 258 static INLINE res_T 259 ncvar_real_inq 260 (int nc, 261 const char* var_name, 262 const size_t expected_ndims, 263 struct var* var) 264 { 265 int i; 266 int dimids[NDIMS_MAX]; /* NetCDF id of the data dimension */ 267 int err; 268 res_T res = RES_OK; 269 ASSERT(var_name && var && expected_ndims && expected_ndims <= NDIMS_MAX); 270 271 /* Retrieve the NetCDF id of the variable */ 272 err = nc_inq_varid(nc, var_name, &var->id); 273 if(err != NC_NOERR) { 274 fprintf(stderr, "Could not inquire the '%s' variable -- %s\n", 275 var_name, ncerr_to_str(err)); 276 res = RES_BAD_ARG; 277 goto error; 278 } 279 280 /* Inquire the dimensions of the variable */ 281 NC(inq_varndims(nc, var->id, &var->ndims)); 282 if((size_t)var->ndims != expected_ndims) { 283 fprintf(stderr, "The dimension of the '%s' variable must be %lu .\n", 284 var_name, (unsigned long)(expected_ndims)); 285 res = RES_BAD_ARG; 286 goto error; 287 } 288 NC(inq_vardimid(nc, var->id, dimids)); 289 FOR_EACH(i, 0, var->ndims) NC(inq_dimlen(nc, dimids[i], var->dim_len+i)); 290 291 /* Inquire the type of the variable */ 292 NC(inq_vartype(nc, var->id, &var->type)); 293 if(var->type != NC_DOUBLE && var->type != NC_FLOAT) { 294 fprintf(stderr, 295 "The type of the '%s' variable cannot be %s. Expecting floating point data.\n", 296 var_name, nctype_to_str(var->type)); 297 res = RES_BAD_ARG; 298 goto error; 299 } 300 301 exit: 302 return res; 303 error: 304 goto exit; 305 } 306 307 /******************************************************************************* 308 * Helper functions 309 ******************************************************************************/ 310 static res_T 311 open_output_stream(const char* path, const int force_overwrite, FILE** stream) 312 { 313 int fd = -1; 314 FILE* fp = NULL; 315 res_T res = RES_OK; 316 ASSERT(path && stream); 317 318 if(force_overwrite) { 319 fp = fopen(path, "w"); 320 if(!fp) { 321 fprintf(stderr, "Could not open the output file '%s'.\n", path); 322 goto error; 323 } 324 } else { 325 fd = open(path, O_CREAT|O_WRONLY|O_EXCL|O_TRUNC, S_IRUSR|S_IWUSR); 326 if(fd >= 0) { 327 fp = fdopen(fd, "w"); 328 if(fp == NULL) { 329 fprintf(stderr, "Could not open the output file '%s'.\n", path); 330 goto error; 331 } 332 } else if(errno == EEXIST) { 333 fprintf(stderr, 334 "The output file '%s' already exists. Use -f to overwrite it.\n", path); 335 goto error; 336 } else { 337 fprintf(stderr, 338 "Unexpected error while opening the output file `'%s'.\n", path); 339 goto error; 340 } 341 } 342 343 exit: 344 *stream = fp; 345 return res; 346 error: 347 res = RES_IO_ERR; 348 if(fp) { 349 CHK(fclose(fp) == 0); 350 fp = NULL; 351 } else if(fd >= 0) { 352 CHK(close(fd) == 0); 353 } 354 goto exit; 355 } 356 357 358 static res_T 359 setup_definition 360 (int nc, int32_t* nx, int32_t* ny, int32_t* nz, int32_t* ntimes) 361 { 362 size_t len[4]; 363 int dimids[4]; 364 int id; 365 int ndims; 366 int err = NC_NOERR; 367 res_T res = RES_OK; 368 ASSERT(nx && ny && nz && ntimes); 369 370 err = nc_inq_varid(nc, "RCT", &id); 371 if(err != NC_NOERR) { 372 fprintf(stderr, "Could not inquire the 'RCT' variable -- %s\n", 373 ncerr_to_str(err)); 374 res = RES_BAD_ARG; 375 goto error; 376 } 377 378 NC(inq_varndims(nc, id, &ndims)); 379 if(ndims != 4) { 380 fprintf(stderr, "The dimension of the 'RCT' variable must be 4.\n"); 381 res = RES_BAD_ARG; 382 goto error; 383 } 384 385 NC(inq_vardimid(nc, id, dimids)); 386 NC(inq_dimlen(nc, dimids[0], len+0)); 387 NC(inq_dimlen(nc, dimids[1], len+1)); 388 NC(inq_dimlen(nc, dimids[2], len+2)); 389 NC(inq_dimlen(nc, dimids[3], len+3)); 390 391 *nx = (int32_t)len[3]; 392 *ny = (int32_t)len[2]; 393 *nz = (int32_t)len[1]; 394 *ntimes = (int32_t)len[0]; 395 396 if(*nx < 1 || *ny < 1 || *nz < 1 || *ntimes < 1) { 397 fprintf(stderr, 398 "The spatial and time definitions cannot be null.\n" 399 " #x = %i; #y = %i; #z = %i; #time = %i\n", 400 *nx, *ny, *nz, *ntimes); 401 res = RES_BAD_ARG; 402 goto error; 403 } 404 405 exit: 406 return res; 407 error: 408 goto exit; 409 } 410 411 static res_T 412 setup_regular_dimension 413 (const int nc, 414 const char* var_name, 415 double* voxel_size, 416 double* lower_pos, 417 const int check) 418 { 419 struct var var = VAR_NULL; 420 double* mem = NULL; 421 size_t start, i; 422 double vxsz; 423 res_T res = RES_OK; 424 ASSERT(nc && var_name && voxel_size && lower_pos); 425 426 res = ncvar_real_inq(nc, var_name, 1, &var); 427 if(res != RES_OK) goto error; 428 429 if(!check) { 430 const size_t len = 1; 431 start = 0; 432 NC(get_vara_double(nc, var.id, &start, &len, &vxsz)); 433 /* Assume that the coordinate is defined at the center of the cell and that 434 * the grid starts at 0 */ 435 vxsz *= 2; 436 } else { 437 mem = mem_alloc(var.dim_len[0]*sizeof(double)); 438 if(!mem) { 439 fprintf(stderr, "Could not allocate memory for the '%s' variable.\n", var_name); 440 res = RES_MEM_ERR; 441 goto error; 442 } 443 start = 0; 444 NC(get_vara_double(nc, var.id, &start, &var.dim_len[0], mem)); 445 /* Assume that the coordinate is defined at the center of the cell and that 446 * the grid starts at 0 */ 447 vxsz = mem[0] * 2; 448 } 449 450 if(vxsz <= 0) { 451 fprintf(stderr, "The '%s' variable can't have negative or null values.\n", 452 var_name); 453 res = RES_BAD_ARG; 454 goto error; 455 } 456 457 if(check) { 458 /* Check that the dimension is regular */ 459 FOR_EACH(i, 1, var.dim_len[0]) { 460 if(!eq_eps(mem[i] - mem[i-1], vxsz, 1.e-6)) { 461 fprintf(stderr, "The voxel size of the '%s' variable must be regular.\n", 462 var_name); 463 res = RES_BAD_ARG; 464 goto error; 465 } 466 } 467 } 468 469 *voxel_size = vxsz; 470 *lower_pos = 0; 471 exit: 472 if(mem) mem_rm(mem); 473 return res; 474 error: 475 goto exit; 476 } 477 478 static res_T 479 compute_Z_voxel_size 480 (const char* var_name, 481 const double* lvls, 482 const size_t nlvls, 483 double** voxel_size, 484 double* lower_pos, 485 int8_t* is_z_irregular) 486 { 487 double vxsz; /* Voxel size */ 488 double* pvxsz = NULL; /* Pointer toward voxel sizes */ 489 size_t z; 490 int irregular; 491 res_T res = RES_OK; 492 ASSERT(var_name && lvls && nlvls && voxel_size && lower_pos && is_z_irregular); 493 494 /* Assume that the coordinate is defined at the center of the cell and that 495 * the grid starts at 0 */ 496 vxsz = lvls[0] * 2.0; 497 if(vxsz <= 0) { 498 fprintf(stderr, 499 "The '%s' variable can't have negative or null values.\n", var_name); 500 res = RES_BAD_ARG; 501 goto error; 502 } 503 504 /* Check if the Z dimension is regular */ 505 FOR_EACH(z, 1, nlvls) { 506 if(!eq_eps(lvls[z] - lvls[z-1], vxsz, 1.e-6)) break; 507 } 508 irregular = (z != nlvls); 509 510 pvxsz = mem_alloc(sizeof(double) * (irregular ? nlvls : 1)); 511 if(!pvxsz) { 512 fprintf(stderr, 513 "Could not allocate memory for the voxel size along the Z dimension.\n"); 514 res = RES_MEM_ERR; 515 goto error; 516 } 517 pvxsz[0] = vxsz; 518 519 if(irregular) { 520 FOR_EACH(z, 1, nlvls) { 521 pvxsz[z] = (lvls[z] - (lvls[z-1] + pvxsz[z-1]*0.5)) * 2.0; 522 if(pvxsz[z] <= 0) { 523 fprintf(stderr, 524 "The '%s' variable can't have negative or null values.\n", var_name); 525 res = RES_BAD_ARG; 526 goto error; 527 } 528 } 529 } 530 531 *lower_pos = 0; 532 *voxel_size = pvxsz; 533 *is_z_irregular = (int8_t)irregular; 534 535 exit: 536 return res; 537 error: 538 if(pvxsz) mem_rm(pvxsz); 539 goto exit; 540 } 541 542 static res_T 543 setup_Z_dimension_VLEV 544 (const int nc, 545 double** voxel_size, 546 double* lower_pos, 547 int8_t* is_z_irregular, 548 const int check) 549 { 550 enum { Z, Y, X, NDIMS }; /* Helper constants */ 551 struct var var = VAR_NULL; 552 double* mem = NULL; /* Memory where NetCDF data are copied */ 553 double* mem2 = NULL; /* Temporary variable used to check NetCDF data */ 554 size_t len[NDIMS]; 555 size_t start[NDIMS]; 556 res_T res = RES_OK; 557 ASSERT(nc && voxel_size && lower_pos && is_z_irregular); 558 559 /* Retrieve the VLEV variable */ 560 res = ncvar_real_inq(nc, "VLEV", NDIMS, &var); 561 if(res != RES_OK) goto error; 562 563 /* Allocate te memory space where a Z column of the VLEV variable is going to 564 * be copied */ 565 mem = mem_alloc(var.dim_len[Z]*sizeof(double)); 566 if(!mem) { 567 fprintf(stderr, "Could not allocate memory for the 'VLEV' variable.\n"); 568 res = RES_MEM_ERR; 569 goto error; 570 } 571 572 /* Get only one Z column */ 573 start[X] = 0; len[X] = 1; 574 start[Y] = 0; len[Y] = 1; 575 start[Z] = 0; len[Z] = var.dim_len[Z]; 576 NC(get_vara_double(nc, var.id, start, len, mem)); 577 578 res = compute_Z_voxel_size("VLEV", mem, var.dim_len[Z], voxel_size, 579 lower_pos, is_z_irregular); 580 if(res != RES_OK) goto error; 581 582 /* Check that the others columns have the same voxel size */ 583 if(check) { 584 size_t ix, iy, iz; 585 586 mem2 = mem_alloc(var.dim_len[Z]*sizeof(double)); 587 if(!mem2) { 588 fprintf(stderr, 589 "Could not allocate memory to check the 'VLEV' variable.\n"); 590 res = RES_MEM_ERR; 591 goto error; 592 } 593 594 FOR_EACH(iy, 0, var.dim_len[Y]) { 595 start[Y] = iy; 596 FOR_EACH(ix, 1, var.dim_len[X]) { 597 start[X] = ix; 598 NC(get_vara_double(nc, var.id, start, len, mem2)); 599 FOR_EACH(iz, 0, var.dim_len[Z]) { 600 if(!eq_eps(mem2[iz], mem[iz], 1.e-6)) { 601 fprintf(stderr, 602 "The Z columns of the 'VLEV' variable must be the same.\n"); 603 res = RES_BAD_ARG; 604 goto error; 605 } 606 } 607 } 608 } 609 } 610 611 exit: 612 if(mem2) mem_rm(mem2); 613 if(mem) mem_rm(mem); 614 return res; 615 error: 616 goto exit; 617 } 618 619 static res_T 620 setup_Z_dimension_vertical_levels 621 (const int nc, 622 double** voxel_size, 623 double* lower_pos, 624 int8_t* is_z_irregular, 625 const int check) 626 { 627 struct var var = VAR_NULL; 628 double* mem = NULL; 629 size_t start = 0; 630 res_T res = RES_OK; 631 (void)check; 632 ASSERT(nc && voxel_size && lower_pos); 633 634 res = ncvar_real_inq(nc, "vertical_levels", 1, &var); 635 if(res != RES_OK) goto error; 636 637 mem = mem_alloc(var.dim_len[0]*sizeof(double)); 638 if(!mem) { 639 fprintf(stderr, 640 "Could not allocate memory for the 'vertical_levels' variable.\n"); 641 res = RES_MEM_ERR; 642 goto error; 643 } 644 645 NC(get_vara_double(nc, var.id, &start, &var.dim_len[0], mem)); 646 647 res = compute_Z_voxel_size("vertical_levels", mem, var.dim_len[0], 648 voxel_size, lower_pos, is_z_irregular); 649 if(res != RES_OK) goto error; 650 651 exit: 652 if(mem) mem_rm(mem); 653 return res; 654 error: 655 goto exit; 656 } 657 658 static res_T 659 setup_Z_dimension 660 (const int nc, 661 double** voxel_size, 662 double* lower_pos, 663 int8_t* is_z_irregular, 664 const int check) 665 { 666 int err; 667 int id; 668 res_T res = RES_OK; 669 ASSERT(nc); 670 671 err = nc_inq_varid(nc, "VLEV", &id); 672 if(err == NC_NOERR) { 673 res = setup_Z_dimension_VLEV 674 (nc, voxel_size, lower_pos, is_z_irregular, check); 675 if(res != RES_OK) goto error; 676 } else if(err == NC_ENOTVAR) { 677 err = nc_inq_varid(nc, "vertical_levels", &id); 678 if(err != NC_NOERR) { 679 fprintf(stderr, 680 "Could not inquire neither the 'VLEV' nor the 'vertical_levels' " 681 "variable -- %s\n", ncerr_to_str(err)); 682 res = RES_BAD_ARG; 683 goto error; 684 } 685 res = setup_Z_dimension_vertical_levels 686 (nc, voxel_size, lower_pos, is_z_irregular, check); 687 if(res != RES_OK) goto error; 688 } 689 exit: 690 return res; 691 error: 692 goto exit; 693 } 694 695 /* Functor used to convert NetCDF data */ 696 typedef void (*convert_data_T) 697 (double* data, /* Input data */ 698 const size_t start[4], /* Start Index of the data in the NetCDF */ 699 const size_t len[4], /* Lenght of the loaded data */ 700 void* context); /* User defined data */ 701 702 /* Write data of the variable "var_name" into "stream". The data are loaded 703 * slice by slice, where each slice is a 2D array whose dimensions are 704 * "W_E_direction X S_N_direction". The "convert" functor is invoked on each 705 * loaded slice. */ 706 static res_T 707 write_data 708 (int nc, 709 const char* var_name, 710 FILE* stream, 711 convert_data_T convert_func, /* May be NULL <=> No conversion */ 712 void* convert_ctx) /* Sent as the last argument of "convert_func" */ 713 { 714 enum { TIME, Z, Y, X, NDIMS }; /* Helper constants */ 715 716 struct var var = VAR_NULL; 717 double* mem = NULL; /* Memory where NetCDF data are copied */ 718 size_t start[NDIMS]; 719 size_t len[NDIMS]; 720 size_t slice_len; /* #elements per slice the 3D grid */ 721 size_t t; /* Id over the time dimension */ 722 size_t z; /* Id over the Z dimension */ 723 res_T res = RES_OK; 724 CHK(var_name); 725 726 res = ncvar_real_inq(nc, var_name, NDIMS, &var); 727 if(res != RES_OK) goto error; 728 729 /* Alloc the memory space where the variable data will be copied. Note that, 730 * in order to reduce memory consumption, only one slice is read at a given 731 * time and Z position. */ 732 slice_len = var.dim_len[X]*var.dim_len[Y]; /* # elements per slice */ 733 mem = mem_alloc(slice_len*sizeof(double)); 734 735 if(!mem) { 736 fprintf(stderr, "Could not allocate memory for the '%s' variable.\n", 737 var_name); 738 res = RES_MEM_ERR; 739 goto error; 740 } 741 742 /* Read the whole slice */ 743 start[X] = 0; len[X] = var.dim_len[X]; 744 start[Y] = 0; len[Y] = var.dim_len[Y]; 745 746 FOR_EACH(t, 0, var.dim_len[TIME]) { 747 start[TIME] = t, len[TIME] = 1; 748 FOR_EACH(z, 0, var.dim_len[Z]) { 749 size_t n; 750 start[Z] = z; len[Z] = 1; 751 752 NC(get_vara_double(nc, var.id, start, len, mem)); /* Fetch the slice data */ 753 754 if(convert_func) convert_func(mem, start, len, convert_ctx); 755 756 /* Write data */ 757 n = fwrite(mem, sizeof(double), slice_len, stream); 758 if(n != slice_len) { 759 fprintf(stderr, "Error writing data of the '%s' variable -- '%s'.\n", 760 var_name, strerror(errno)); 761 res = RES_IO_ERR; 762 goto error; 763 } 764 } 765 } 766 767 exit: 768 if(mem) mem_rm(mem); 769 return res; 770 error: 771 goto exit; 772 } 773 774 struct THT_to_T_context { 775 double* mem; /* Memory where to store a PABST slice */ 776 size_t len; /* lenght of the mem */ 777 int pabst_id; /* NetCDF id of the PABST variable */ 778 int nc; /* NetCDF handler */ 779 }; 780 static const struct THT_to_T_context THT_TO_T_CONTEXT_NULL = { NULL, 0, 0, 0 }; 781 782 static void 783 THT_to_T_context_release(struct THT_to_T_context* ctx) 784 { 785 ASSERT(ctx); 786 if(ctx->mem) mem_rm(ctx->mem); 787 } 788 789 static res_T 790 THT_to_T_context_init(struct THT_to_T_context* ctx, int nc) 791 { 792 enum { TIME, Z, Y, X, NDIMS }; /* Helper constants */ 793 struct var var = VAR_NULL; 794 res_T res = RES_OK; 795 ASSERT(ctx); 796 797 *ctx = THT_TO_T_CONTEXT_NULL; 798 799 res = ncvar_real_inq(nc, "PABST", NDIMS, &var); 800 if(res != RES_OK) goto error; 801 802 /* Allocate memory to store one slice of PABST data */ 803 ctx->len = var.dim_len[X] * var.dim_len[Y]; 804 ctx->mem = mem_alloc(ctx->len*sizeof(double)); 805 if(!ctx->mem) { 806 fprintf(stderr, 807 "Could not allocate the temporary memory space of the THT_to_T_context " 808 "data structure.\n"); 809 res = RES_MEM_ERR; 810 goto error; 811 } 812 813 ctx->pabst_id = var.id; 814 ctx->nc = nc; 815 816 exit: 817 return res; 818 error: 819 THT_to_T_context_release(ctx); 820 goto exit; 821 } 822 823 /* Convert potential temperature to temperature 824 * T(i)=THT(i)*(PABST(i)/P0)^(R/(M_DA*Cp0)) */ 825 static void 826 THT_to_T 827 (double* THT, 828 const size_t start[4], /* Start id of the THT var in the NetCDF file */ 829 const size_t len[4], /* Len of the loaded THT */ 830 void* context) 831 { 832 struct THT_to_T_context* ctx = context; 833 const double P0 = 101325; /* In Pa */ 834 const double R = 8.3144598; /* In kg.m^2.s^-2.mol^-1.K */ 835 const double M_DA = 28.9644e-3; /* In kg.mol^-1 */ 836 const double Cp0 = 1.005e+3; /* In J/kg/K */ 837 const double exponent = R/(M_DA * Cp0); 838 double* T = THT; 839 double* PABST; 840 size_t i, n; 841 ASSERT(THT && start && len && context); 842 843 n = len[0] * len[1] * len[2] * len[3]; 844 CHK(n <= ctx->len); 845 PABST = ctx->mem; 846 847 /* Fetch the slice pressure data */ 848 NC(get_vara_double(ctx->nc, ctx->pabst_id, start, len, PABST)); 849 850 FOR_EACH(i, 0, n) T[i] = THT[i] * pow(PABST[i]/P0, exponent); 851 } 852 853 /******************************************************************************* 854 * Program 855 ******************************************************************************/ 856 int 857 main(int argc, char** argv) 858 { 859 FILE* stream = stdout; 860 struct grid grid = GRID_NULL; 861 struct args args = ARGS_DEFAULT; 862 struct THT_to_T_context THT_to_T_ctx = THT_TO_T_CONTEXT_NULL; 863 864 int err = 0; 865 int err_nc = NC_NOERR; 866 int nc = INVALID_ID; 867 res_T res = RES_OK; 868 869 res = args_init(&args, argc, argv); 870 if(res != RES_OK) goto error; 871 if(args.quit) goto exit; 872 873 if(args.output) { 874 res = open_output_stream(args.output, args.force_overwrite, &stream); 875 if(res != RES_OK) goto error; 876 } 877 878 err_nc = nc_open(args.input, NC_NOWRITE, &nc); 879 if(err_nc != NC_NOERR) { 880 fprintf(stderr, "Error opening file `%s' -- %s.\n", 881 args.input, ncerr_to_str(err_nc)); 882 goto exit; 883 } 884 885 #define CALL(Func) { if(RES_OK != (res = Func)) goto error; } (void)0 886 887 CALL(setup_definition(nc, &grid.nx, &grid.ny, &grid.nz, &grid.ntimes)); 888 CALL(setup_regular_dimension 889 (nc, "W_E_direction", &grid.vxsz_x, grid.lower+0, args.check)); 890 CALL(setup_regular_dimension 891 (nc, "S_N_direction", &grid.vxsz_y, grid.lower+1, args.check)); 892 CALL(setup_Z_dimension 893 (nc, &grid.vxsz_z, grid.lower+2, &grid.is_z_irregular, args.check)); 894 895 if(args.fp_to_meter != 1) { /* Convert the grid in meters */ 896 grid.lower[0] *= args.fp_to_meter; 897 grid.lower[1] *= args.fp_to_meter; 898 grid.lower[2] *= args.fp_to_meter; 899 900 grid.vxsz_x *= args.fp_to_meter; 901 grid.vxsz_y *= args.fp_to_meter; 902 if(!grid.is_z_irregular) { 903 grid.vxsz_z[0] *= args.fp_to_meter; 904 } else { 905 int32_t i; 906 FOR_EACH(i, 0, grid.nz) { 907 grid.vxsz_z[i] *= args.fp_to_meter; 908 } 909 } 910 } 911 912 #define WRITE(Var, N, Name) { \ 913 if(fwrite((Var), sizeof(*(Var)), (N), stream) != (N)) { \ 914 fprintf(stderr, "Error writing the %s.\n", (Name)); \ 915 res = RES_IO_ERR; \ 916 goto error; \ 917 } \ 918 } (void)0 919 if(!args.no_output) { 920 const int64_t pagesz = (int64_t)args.pagesize; 921 922 WRITE(&pagesz, 1, "page size"); 923 WRITE(&grid.is_z_irregular, 1, "'irregular' Z boolean"); 924 WRITE(&grid.nx, 1, "X definition"); 925 WRITE(&grid.ny, 1, "Y definition"); 926 WRITE(&grid.nz, 1, "Z definition"); 927 WRITE(&grid.ntimes, 1, "time definition"); 928 WRITE(grid.lower, 3, "lower position"); 929 WRITE(&grid.vxsz_x, 1, "X voxel size"); 930 WRITE(&grid.vxsz_y, 1, "Y voxel size"); 931 WRITE(grid.vxsz_z, grid.is_z_irregular?(size_t)grid.nz:1, "Z voxel size(s)"); 932 fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/ 933 934 CALL(write_data(nc, "RVT", stream, NULL, NULL)); 935 fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/ 936 CALL(write_data(nc, "RCT", stream, NULL, NULL)); 937 fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/ 938 CALL(write_data(nc, "PABST", stream, NULL, NULL)); 939 fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/ 940 941 CALL(THT_to_T_context_init(&THT_to_T_ctx, nc)); 942 CALL(write_data(nc, "THT", stream, THT_to_T, &THT_to_T_ctx)); 943 if(!IS_ALIGNED(ftell(stream), (size_t)pagesz)) { 944 /* Padding to ensure that the size is aligned on the page size. Note that 945 * one char is written to position the EOF indicator */ 946 const char byte = 0; 947 fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz)-1, SEEK_SET); 948 WRITE(&byte, 1, "Dummy Byte"); 949 } 950 } 951 #undef WRITE 952 #undef CALL 953 954 exit: 955 grid_release(&grid); 956 THT_to_T_context_release(&THT_to_T_ctx); 957 if(stream && stream != stdout) CHK(fclose(stream) == 0); 958 if(nc != INVALID_ID) NC(close(nc)); 959 if(mem_allocated_size() != 0) { 960 fprintf(stderr, "Memory leaks: %lu Bytes\n", 961 (unsigned long)mem_allocated_size()); 962 err = -1; 963 } 964 return err; 965 error: 966 err = -1; 967 goto exit; 968 } 969