commit 04fc6e940e71212626127f68503f70da3a225ed2
parent 7b1c65b4507b212597de99c3d54565fa07c91f12
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 28 Sep 2018 15:13:10 +0200
Handle LES with no VLEV
If the "VLEV" variable is not found in the LES, use the "vertical_levels"
variable to compute the voxel sizes along the Z dimension.
Diffstat:
2 files changed, 151 insertions(+), 52 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -71,7 +71,8 @@ if(NOT NO_TEST)
CEN2D.1.ARMCU.008.diaKCL
CERK4.1.ARMCu.008
DZVAR.1.ARMCU.008.diaKCL
- L25km.1.ARMCU.008.diaKCL)
+ L25km.1.ARMCU.008.diaKCL
+ L51km.1.ARMCu.008)
set(VAR_LISTS RVT RCT PABST THT)
diff --git a/src/les2htcp.c b/src/les2htcp.c
@@ -492,9 +492,72 @@ error:
goto exit;
}
-/* Return 1 if Z is irregular */
-static int
-setup_Z_dimension
+static res_T
+compute_Z_voxel_size
+ (const char* var_name,
+ const double* lvls,
+ const size_t nlvls,
+ double** voxel_size,
+ double* lower_pos,
+ int8_t* is_z_irregular)
+{
+ double vxsz; /* Voxel size */
+ double* pvxsz = NULL; /* Pointer toward voxel sizes */
+ size_t z;
+ int irregular;
+ res_T res = RES_OK;
+ ASSERT(var_name && lvls && nlvls && voxel_size && lower_pos && is_z_irregular);
+
+ /* Assume that the coordinate is defined at the center of the cell and that
+ * the grid starts at 0 */
+ vxsz = lvls[0] * 2.0;
+ if(vxsz <= 0) {
+ fprintf(stderr,
+ "The '%s' variable can't have negative or null values.\n", var_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Check if the Z dimension is regular */
+ FOR_EACH(z, 1, nlvls) {
+ if(!eq_eps(lvls[z] - lvls[z-1], vxsz, 1.e-6)) break;
+ }
+ irregular = (z != nlvls);
+
+ pvxsz = mem_alloc(sizeof(double) * (irregular ? nlvls : 1));
+ if(!pvxsz) {
+ fprintf(stderr,
+ "Could not allocate memory for the voxel size along the Z dimension.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ pvxsz[0] = vxsz;
+
+ if(irregular) {
+ FOR_EACH(z, 1, nlvls) {
+ pvxsz[z] = (lvls[z] - (lvls[z-1] + pvxsz[z-1]*0.5)) * 2.0;
+ if(pvxsz[z] <= 0) {
+ fprintf(stderr,
+ "The '%s' variable can't have negative or null values.\n", var_name);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ }
+ }
+
+ *lower_pos = 0;
+ *voxel_size = pvxsz;
+ *is_z_irregular = (int8_t)irregular;
+
+exit:
+ return res;
+error:
+ if(pvxsz) mem_rm(pvxsz);
+ goto exit;
+}
+
+static res_T
+setup_Z_dimension_VLEV
(const int nc,
double** voxel_size,
double* lower_pos,
@@ -507,10 +570,6 @@ setup_Z_dimension
double* mem2 = NULL; /* Temporary variable used to check NetCDF data */
size_t len[NDIMS];
size_t start[NDIMS];
- size_t z; /* Id over the Z dimension */
- double vxsz; /* Voxel size */
- double* pvxsz = NULL; /* Pointer toward voxel sizes */
- int irregular; /* Boolean defining if the Z dimension is irregular */
res_T res = RES_OK;
ASSERT(nc && voxel_size && lower_pos && is_z_irregular);
@@ -533,46 +592,13 @@ setup_Z_dimension
start[Z] = 0; len[Z] = var.dim_len[Z];
NC(get_vara_double(nc, var.id, start, len, mem));
- /* Assume that the coordinate is defined at the center of the cell and that
- * the grid starts at 0 */
- vxsz = mem[0] * 2.0;
- if(vxsz <= 0) {
- fprintf(stderr, "The 'VLEV' variable can't have negative or null values.\n");
- res = RES_BAD_ARG;
- goto error;
- }
-
- /* Check if the Z dimension is regular */
- FOR_EACH(z, 1, var.dim_len[Z]) {
- if(!eq_eps(mem[z] - mem[z-1], vxsz, 1.e-6)) break;
- }
- irregular = (z != var.dim_len[Z]);
-
- /* Setup the size of the voxel */
- pvxsz = mem_alloc(sizeof(double) * (irregular ? var.dim_len[Z] : 1));
- if(!pvxsz) {
- fprintf(stderr,
- "Could not allocate memory for the voxel size along the Z dimension.\n");
- res = RES_MEM_ERR;
- goto error;
- }
- pvxsz[0] = vxsz;
-
- if(irregular) {
- FOR_EACH(z, 1, var.dim_len[Z]) {
- pvxsz[z] = (mem[z] - (mem[z-1] + pvxsz[z-1]*0.5)) * 2.0;
- if(pvxsz[z] <= 0) {
- fprintf(stderr,
- "The 'VLEV' variable can't have negative or null values.\n");
- res = RES_BAD_ARG;
- goto error;
- }
- }
- }
+ res = compute_Z_voxel_size("VLEV", mem, var.dim_len[Z], voxel_size,
+ lower_pos, is_z_irregular);
+ if(res != RES_OK) goto error;
/* Check that the others columns have the same voxel size */
if(check) {
- size_t ix, iy;
+ size_t ix, iy, iz;
mem2 = mem_alloc(var.dim_len[Z]*sizeof(double));
if(!mem2) {
@@ -587,8 +613,8 @@ setup_Z_dimension
FOR_EACH(ix, 1, var.dim_len[X]) {
start[X] = ix;
NC(get_vara_double(nc, var.id, start, len, mem2));
- FOR_EACH(z, 0, var.dim_len[Z]) {
- if(!eq_eps(mem2[z], mem[z], 1.e-6)) {
+ FOR_EACH(iz, 0, var.dim_len[Z]) {
+ if(!eq_eps(mem2[iz], mem[iz], 1.e-6)) {
fprintf(stderr,
"The Z columns of the 'VLEV' variable must be the same.\n");
res = RES_BAD_ARG;
@@ -599,10 +625,6 @@ setup_Z_dimension
}
}
- *lower_pos = 0;
- *voxel_size = pvxsz;
- *is_z_irregular = (int8_t)irregular;
-
exit:
if(mem2) mem_rm(mem2);
if(mem) mem_rm(mem);
@@ -611,6 +633,82 @@ error:
goto exit;
}
+static res_T
+setup_Z_dimension_vertical_levels
+ (const int nc,
+ double** voxel_size,
+ double* lower_pos,
+ int8_t* is_z_irregular,
+ const int check)
+{
+ struct var var = VAR_NULL;
+ double* mem = NULL;
+ size_t start = 0;
+ res_T res = RES_OK;
+ (void)check;
+ ASSERT(nc && voxel_size && lower_pos);
+
+ res = ncvar_real_inq(nc, "vertical_levels", 1, &var);
+ if(res != RES_OK) goto error;
+
+ mem = mem_alloc(var.dim_len[0]*sizeof(double));
+ if(!mem) {
+ fprintf(stderr,
+ "Could not allocate memory for the 'vertical_levels' variable.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ NC(get_vara_double(nc, var.id, &start, &var.dim_len[0], mem));
+
+ res = compute_Z_voxel_size("vertical_levels", mem, var.dim_len[0],
+ voxel_size, lower_pos, is_z_irregular);
+ if(res != RES_OK) goto error;
+
+exit:
+ if(mem) mem_rm(mem);
+ return res;
+error:
+ goto exit;
+}
+
+static res_T
+setup_Z_dimension
+ (const int nc,
+ double** voxel_size,
+ double* lower_pos,
+ int8_t* is_z_irregular,
+ const int check)
+{
+ int err;
+ int id;
+ res_T res = RES_OK;
+ ASSERT(nc);
+
+ err = nc_inq_varid(nc, "VLEV", &id);
+ if(err == NC_NOERR) {
+ res = setup_Z_dimension_VLEV
+ (nc, voxel_size, lower_pos, is_z_irregular, check);
+ if(res != RES_OK) goto error;
+ } else if(err == NC_ENOTVAR) {
+ err = nc_inq_varid(nc, "vertical_levels", &id);
+ if(err != NC_NOERR) {
+ fprintf(stderr,
+ "Could not inquire neither the 'VLEV' nor the 'vertical_levels' "
+ "variable -- %s\n", ncerr_to_str(err));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ res = setup_Z_dimension_vertical_levels
+ (nc, voxel_size, lower_pos, is_z_irregular, check);
+ if(res != RES_OK) goto error;
+ }
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/* Functor used to convert NetCDF data */
typedef void (*convert_data_T)
(double* data, /* Input data */
@@ -794,7 +892,7 @@ main(int argc, char** argv)
if(res != RES_OK) goto error;
}
- err_nc = nc_open(args.input, NC_WRITE|NC_MMAP, &nc);
+ err_nc = nc_open(args.input, NC_NOWRITE, &nc);
if(err_nc != NC_NOERR) {
fprintf(stderr, "Error opening file `%s' -- %s.\n",
args.input, ncerr_to_str(err_nc));
