commit 09a981163886fdfb7c5c402f2ecc1a604d00d8b2
parent 4ec8ec7374795728edb3621c12b27c4f1ac12d5c
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 30 Apr 2018 15:28:36 +0200
Fix the les2htcop program
LES data saved as float are now corectly converted in double in the
htcop output file.
Diffstat:
| M | src/les2htcop.c | | | 125 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------- |
1 file changed, 85 insertions(+), 40 deletions(-)
diff --git a/src/les2htcop.c b/src/les2htcop.c
@@ -455,19 +455,26 @@ setup_Z_dimension
int8_t* is_z_irregular,
const int check)
{
- enum { Z, X, Y, NDIMS };
- char* mem = NULL;
- char* mem2 = NULL;
- double* pvxsz = NULL;
- size_t len[NDIMS], start[NDIMS], dim_len[NDIMS], iz;
- size_t i;
- double vxsz;
- nc_type type;
- int varid, dimids[NDIMS], irregular, ndims;
+ enum { Z, X, Y, NDIMS }; /* Helper constants */
+ char* mem = NULL; /* Memory where NetCDF data are copied */
+ char* mem2 = NULL; /* Temporary memory used to check the VLEV data */
+ size_t len[NDIMS];
+ size_t start[NDIMS];
+ size_t dim_len[NDIMS];
+ size_t z; /* Id over the Z dimension */
+ size_t i; /* Common id */
+ double vxsz; /* Voxel size */
+ double* pvxsz = NULL; /* Pointer toward voxel sizes */
+ nc_type type; /* Type of the NetCDF VLEV variable */
+ int varid; /* NetCDF id of the VLEV variable */
+ int dimids[NDIMS]; /* NetCDF id of the VLEV dimensions */
+ int irregular; /* Boolean defining if the Z dimension is irregular */
+ int ndims; /* #dims of the VLEV variable */
int err = NC_NOERR;
res_T res = RES_OK;
ASSERT(nc && voxel_size && lower_pos && is_z_irregular);
+ /* Retrieve the VLEV variable */
err = nc_inq_varid(nc, "VLEV", &varid);
if(err != NC_NOERR) {
fprintf(stderr, "Could not inquire the 'VLEV' variable -- %s\n",
@@ -476,6 +483,7 @@ setup_Z_dimension
goto error;
}
+ /* Inquire the VLEV dimensionality */
NC(inq_varndims(nc, varid, &ndims));
if(ndims != NDIMS) {
fprintf(stderr, "The dimension of the 'VLEV' variable must be %i .\n",
@@ -483,10 +491,10 @@ setup_Z_dimension
res = RES_BAD_ARG;
goto error;
}
-
NC(inq_vardimid(nc, varid, dimids));
FOR_EACH(i, 0, NDIMS) NC(inq_dimlen(nc, dimids[i], len+i));
+ /* Inquire the type of the VLEV variable */
NC(inq_vartype(nc, varid, &type));
if(type != NC_DOUBLE && type != NC_FLOAT) {
fprintf(stderr, "The type of 'VLEV' variable must be float or double.\n");
@@ -494,6 +502,8 @@ setup_Z_dimension
goto error;
}
+ /* Allocate te memory space where a Z column of the VLEV variable is going to
+ * be copied */
mem = mem_alloc(len[Z]*sizeof_nctype(type));
if(!mem) {
fprintf(stderr, "Could not allocate memory for the 'VLEV' variable.\n");
@@ -501,10 +511,10 @@ setup_Z_dimension
goto error;
}
+ /* Get only one Z column */
start[X] = 0; dim_len[X] = 1;
start[Y] = 0; dim_len[Y] = 1;
start[Z] = 0; dim_len[Z] = len[Z];
-
NC(get_vara(nc, varid, start, dim_len, mem));
#define AT(Mem, Id) (type == NC_FLOAT ? ((float*)Mem)[Id] : ((double*)Mem)[Id])
@@ -518,10 +528,10 @@ setup_Z_dimension
}
/* Check if the Z dimension is regular */
- FOR_EACH(iz, 1, len[Z]) {
- if(!eq_eps(AT(mem, iz) - AT(mem, iz-1), vxsz, 1.e-6)) break;
+ FOR_EACH(z, 1, len[Z]) {
+ if(!eq_eps(AT(mem, z) - AT(mem, z-1), vxsz, 1.e-6)) break;
}
- irregular = (iz != len[Z]);
+ irregular = (z != len[Z]);
/* Setup the size of the voxel */
pvxsz = mem_alloc(sizeof(double) * (irregular ? len[Z] : 1));
@@ -534,9 +544,9 @@ setup_Z_dimension
pvxsz[0] = vxsz;
if(irregular) {
- FOR_EACH(iz, 1, len[Z]) {
- pvxsz[iz] = (AT(mem, iz) - (AT(mem, iz-1) + pvxsz[iz-1]*0.5)) * 2.0;
- if(pvxsz[iz] <= 0) {
+ FOR_EACH(z, 1, len[Z]) {
+ pvxsz[z] = (AT(mem, z) - (AT(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;
@@ -562,8 +572,8 @@ setup_Z_dimension
FOR_EACH(ix, 1, len[X]) {
start[X] = ix;
NC(get_vara(nc, varid, start, dim_len, mem2));
- FOR_EACH(iz, 0, len[Z]) {
- if(!eq_eps(AT(mem2, iz), AT(mem, iz), 1.e-6)) {
+ FOR_EACH(z, 0, len[Z]) {
+ if(!eq_eps(AT(mem2, z), AT(mem, z), 1.e-6)) {
fprintf(stderr,
"The Z columns of the 'VLEV' variable must be the same.\n");
res = RES_BAD_ARG;
@@ -590,15 +600,26 @@ error:
static res_T
write_data(int nc, const char* var_name, FILE* stream)
{
- enum { TIME, Z, X, Y, NDIMS };
- char* mem = NULL;
- size_t start[NDIMS], dim_len[NDIMS], len[NDIMS], i, grid_len;
+ enum { TIME, Z, X, Y, NDIMS }; /* Helper constants */
+
+ char* mem = NULL; /* Memory where NetCDF data are copied */
+ char* mem2 = NULL; /* Tmp memory where NetCDF data are converted */
+ size_t start[NDIMS];
+ size_t dim_len[NDIMS];
+ size_t len[NDIMS];
+ size_t slice_len; /* #elements per slice the 3D grid */
+ size_t i; /* Common id */
+ size_t t; /* Id over the time dimension */
+ size_t z; /* Id over the Z dimension */
nc_type type;
- int varid, ndims, dimids[NDIMS];
+ int varid; /* NetCDF id of the data to write */
+ int ndims; /* #dimensions of the NetCDF data to write */
+ int dimids[NDIMS]; /* NetCDF id of the data dimension */
int err = NC_NOERR;
res_T res = RES_OK;
CHK(var_name);
+ /* Retrieve the NetCDF id of the variable */
err = nc_inq_varid(nc, var_name, &varid);
if(err != NC_NOERR) {
fprintf(stderr, "Could not inquire the '%s' variable -- %s\n",
@@ -607,6 +628,7 @@ write_data(int nc, const char* var_name, FILE* stream)
goto error;
}
+ /* Inquire the dimensions of the variable */
NC(inq_varndims(nc, varid, &ndims));
if(ndims != NDIMS) {
fprintf(stderr, "The dimension of the '%s' variable must be %i .\n",
@@ -614,10 +636,10 @@ write_data(int nc, const char* var_name, FILE* stream)
res = RES_BAD_ARG;
goto error;
}
-
NC(inq_vardimid(nc, varid, dimids));
FOR_EACH(i, 0, NDIMS) NC(inq_dimlen(nc, dimids[i], len+i));
+ /* Inquire the type of the variable */
NC(inq_vartype(nc, varid, &type));
if(type != NC_DOUBLE && type != NC_FLOAT) {
fprintf(stderr,
@@ -627,31 +649,54 @@ write_data(int nc, const char* var_name, FILE* stream)
goto error;
}
- grid_len = len[X]*len[Y]*len[Z];
- mem = mem_alloc(grid_len*sizeof_nctype(type));
- if(!mem) {
+ /* Alloc the memory space where the variable data wille be copied. Note that,
+ * in order to reduce memory consumption, only one slice is read at a given
+ * time and Z position. */
+ slice_len = len[X]*len[Y]; /* # elements per slice */
+ mem = mem_alloc(slice_len*sizeof_nctype(type));
+
+ /* Allocate a temporary memory space to store the variable data converted in
+ * double precision if they are stored in single precision. */
+ mem2 = type==NC_DOUBLE ? mem : mem_alloc(slice_len*sizeof_nctype(NC_DOUBLE));
+
+ if(!mem || !mem2) {
fprintf(stderr, "Could not allocate memory for the '%s' variable.\n",
var_name);
res = RES_MEM_ERR;
goto error;
}
+ /* Read the whole slice */
start[X] = 0; dim_len[X] = len[X];
start[Y] = 0; dim_len[Y] = len[Y];
- start[Z] = 0; dim_len[Z] = len[Z];
- FOR_EACH(i, 0, len[TIME]) {
- start[TIME] = i;
- dim_len[TIME] = 1;
- NC(get_vara(nc, varid, start, dim_len, mem));
- if(fwrite(mem, sizeof_nctype(type), grid_len, stream) != grid_len) {
- fprintf(stderr, "Error writing data of the '%s' variable -- '%s'.\n",
- var_name, strerror(errno));
- res = RES_IO_ERR;
- goto error;
+ FOR_EACH(t, 0, len[TIME]) {
+ start[TIME] = t, dim_len[TIME] = 1;
+ FOR_EACH(z, 0, len[Z]) {
+ size_t n;
+ start[Z] = z; dim_len[Z] = 1;
+
+ NC(get_vara(nc, varid, start, dim_len, mem)); /* Fetch the slice data */
+
+ /* In order to be compliant with the htgop fileformat, single precision
+ * data must be converted in double precision */
+ if(type == NC_FLOAT) {
+ FOR_EACH(i, 0, slice_len) ((double*)mem2)[i] = ((float*)mem)[i];
+ }
+
+ /* Write data */
+ n = fwrite(mem2, sizeof_nctype(NC_DOUBLE), slice_len, stream);
+ if(n != slice_len) {
+ fprintf(stderr, "Error writing data of the '%s' variable -- '%s'.\n",
+ var_name, strerror(errno));
+ res = RES_IO_ERR;
+ goto error;
+ }
}
}
+
exit:
+ if(mem2 != mem) mem_rm(mem2);
if(mem) mem_rm(mem);
return res;
error:
@@ -718,13 +763,13 @@ main(int argc, char** argv)
WRITE(&grid.vxsz_x, 1, "X voxel size");
WRITE(&grid.vxsz_y, 1, "Y voxel size");
WRITE(grid.vxsz_z, grid.is_z_irregular?(size_t)grid.nz:1, "Z voxel size(s)");
-
+ /* padding */
fseek(stream, ALIGN_SIZE(ftell(stream), (off_t)pagesize), SEEK_SET);
write_data(nc, "RVT", stream);
-
+ /* padding */
fseek(stream, ALIGN_SIZE(ftell(stream), (off_t)pagesize), SEEK_SET);
write_data(nc, "RCT", stream);
-
+ /*padding */
fseek(stream, ALIGN_SIZE(ftell(stream), (off_t)pagesize), SEEK_SET);
}
#undef WRITE
