commit 2f1ad5f9b325a653562ba43401b168b50acb2d4a
parent fac14d4552cac3da8127d0d202bc375a085930ab
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 13 Jul 2018 16:51:24 +0200
Update the htcp fileformat
Add the pressure (PABST) and temperature (T) data. Update the les2htcp
tool to add these data to the output HTCP file
Diffstat:
4 files changed, 82 insertions(+), 34 deletions(-)
diff --git a/cmake/les2htcp/CMakeLists.txt b/cmake/les2htcp/CMakeLists.txt
@@ -54,10 +54,9 @@ rcmake_prepend_path(LES2HTCP_FILES_INC ${LES2HTCP_SOURCE_DIR})
rcmake_prepend_path(LES2HTCP_FILES_DOC ${PROJECT_SOURCE_DIR}/../../)
add_executable(les2htcp ${LES2HTCP_FILES_SRC})
-target_link_libraries(les2htcp RSys ${NETCDF_C_LIBRARIES})
+target_link_libraries(les2htcp RSys ${NETCDF_C_LIBRARIES} m)
set_target_properties(les2htcp PROPERTIES
- DEFINE_SYMBOL LES2HTCP_SHARED_BUILD
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
diff --git a/doc/htcp.txt b/doc/htcp.txt
@@ -14,6 +14,10 @@
<padding>
<RCT> # align on pagesize
<padding>
+ <PABST> # align on pagesize
+ <padding>
+ <T> # align on pagesize
+ <padding>
<header> ::= <pagesize> <is-Z-irregular>
<definition> ::= <#X> <#Y> <#Z> <#time>
@@ -21,6 +25,8 @@
<voxel-size> ::= <double3> [DOUBLE ... ] # optionnal double for irregular Z
<RVT> ::= DOUBLE [DOUBLE ... ]
<RCT> ::= DOUBLE [DOUBLE ... ]
+<PABST> ::= DOUBLE [DOUBLE ... ]
+<T> ::= DOUBLE [DOUBLE ... ] # Temperature
<double3> ::= DOUBLE DOUBLE DOUBLE
diff --git a/src/htcp.c b/src/htcp.c
@@ -177,15 +177,15 @@ load_stream(struct htcp* htcp, FILE* stream, const char* stream_name)
offset = ALIGN_SIZE(ftell(stream), htcp->pagesize);
#define MMAP(Var) { \
- htcp->Var = mmap(NULL, map_len, PROT_READ, MAP_PRIVATE|MAP_POPULATE, \
+ htcp->Var = mmap(NULL, map_len, PROT_READ, MAP_PRIVATE|MAP_POPULATE, \
fileno(stream), offset); \
- if(htcp->Var == MAP_FAILED) { \
- log_err(htcp, "%s: could not map the "STR(Var)" data -- %s.\n", \
+ if(htcp->Var == MAP_FAILED) { \
+ log_err(htcp, "%s: could not map the "STR(Var)" data -- %s.\n", \
stream_name, strerror(errno)); \
res = RES_IO_ERR; \
goto error; \
} \
- htcp->Var##_length = map_len; \
+ htcp->Var##_length = map_len; \
} (void)0
MMAP(RVT); offset += (off_t)map_len;
MMAP(RCT); offset += (off_t)map_len;
diff --git a/src/les2htcp.c b/src/les2htcp.c
@@ -13,7 +13,8 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
-#define _POSIX_C_SOURCE 200112L /* getopt/close support */
+#define _POSIX_C_SOURCE 200809L /* mmap/getopt/close support */
+#define _DEFAULT_SOURCE 1 /* MAP_POPULATE support */
#include "les2htcp.h"
@@ -23,10 +24,12 @@
#include <rsys/rsys.h>
#include <alloca.h>
+#include <errno.h>
#include <netcdf.h>
#include <string.h>
#include <unistd.h> /* getopt */
#include <fcntl.h> /* open */
+#include <sys/mman.h> /* mmap */
#include <sys/stat.h> /* S_IRUSR & S_IWUSR */
/*******************************************************************************
@@ -602,8 +605,7 @@ write_data(int nc, const char* var_name, FILE* stream)
{
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 */
+ double* mem = NULL; /* Memory where NetCDF data are copied */
size_t start[NDIMS];
size_t dim_len[NDIMS];
size_t len[NDIMS];
@@ -649,17 +651,13 @@ write_data(int nc, const char* var_name, FILE* stream)
goto error;
}
- /* Alloc the memory space where the variable data wille be copied. Note that,
+ /* Alloc the memory space where the variable data will 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));
+ mem = mem_alloc(slice_len*sizeof(double));
- if(!mem || !mem2) {
+ if(!mem) {
fprintf(stderr, "Could not allocate memory for the '%s' variable.\n",
var_name);
res = RES_MEM_ERR;
@@ -676,16 +674,10 @@ write_data(int nc, const char* var_name, FILE* stream)
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];
- }
+ NC(get_vara_double(nc, varid, start, dim_len, mem)); /* Fetch the slice data */
/* Write data */
- n = fwrite(mem2, sizeof_nctype(NC_DOUBLE), slice_len, stream);
+ n = fwrite(mem, 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));
@@ -696,13 +688,28 @@ write_data(int nc, const char* var_name, FILE* stream)
}
exit:
- if(mem2 != mem) mem_rm(mem2);
if(mem) mem_rm(mem);
return res;
error:
goto exit;
}
+/* Convert potential temperature to temperature
+ * T(i)=THT(i)*(PABST(i)/P0)^(R/(M_DA*Cp0)) */
+static void
+THT_to_T(const size_t n, const double* PABST, double* THT)
+{
+ const double P0 = 101325; /* In Pa */
+ const double R = 8.3144598; /* In kg.m^2.s^-2.mol^-1.K */
+ const double M_DA = 28.9644e-3; /* In kg.mol^-1 */
+ const double Cp0 = 1.005e+3; /* In J/kg/K */
+ const double exponent = R/(M_DA * Cp0);
+ double* T = THT;
+ size_t i;
+ ASSERT(n && PABST && THT);
+ FOR_EACH(i, 0, n) T[i] = THT[i] * pow(PABST[i]/P0, exponent);
+}
+
/*******************************************************************************
* Program
******************************************************************************/
@@ -751,9 +758,14 @@ main(int argc, char** argv)
} \
} (void)0
if(!args.no_output) {
- const int64_t pagesize = (int64_t)args.pagesize;
-
- WRITE(&pagesize, 1, "page size");
+ const int64_t pagesz = (int64_t)args.pagesize;
+ double* PABST = NULL;
+ double* THT = NULL;
+ long PABST_off, THT_off;
+ size_t PABST_len, THT_len;
+ size_t n;
+
+ WRITE(&pagesz, 1, "page size");
WRITE(&grid.is_z_irregular, 1, "'irregular' Z boolean");
WRITE(&grid.nx, 1, "X definition");
WRITE(&grid.ny, 1, "Y definition");
@@ -763,14 +775,45 @@ 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);
+ fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/
+
write_data(nc, "RVT", stream);
- /* padding */
- fseek(stream, ALIGN_SIZE(ftell(stream), (off_t)pagesize), SEEK_SET);
+ fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/
+
write_data(nc, "RCT", stream);
- /*padding */
- fseek(stream, ALIGN_SIZE(ftell(stream), (off_t)pagesize), SEEK_SET);
+ fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/
+
+ PABST_off = ftell(stream);
+ write_data(nc, "PABST", stream);
+ fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/
+ PABST_len = (size_t)(ftell(stream) - PABST_off);
+
+ THT_off = ftell(stream);
+ write_data(nc, "THT", stream);
+ fseek(stream, ALIGN_SIZE(ftell(stream),(off_t)pagesz), SEEK_SET);/*padding*/
+ THT_len = (size_t)(ftell(stream) - THT_off);
+
+ CHK(fclose(stream) == 0);
+ CHK(stream = fopen(args.output, "rw"));
+
+ #define MMAP(Var, Prot) { \
+ Var = mmap(NULL, Var##_len, Prot, MAP_PRIVATE|MAP_POPULATE, \
+ fileno(stream), Var##_off); \
+ if(Var == MAP_FAILED) { \
+ fprintf(stderr, "%s: could not map the "STR(Var)" data -- %s.\n", \
+ args.input, strerror(errno)); \
+ goto error; \
+ } \
+ } (void)0
+ MMAP(PABST, PROT_READ);
+ MMAP(THT, PROT_READ|PROT_WRITE);
+ #undef MMAP
+
+ n = (size_t)(grid.nx * grid.ny * grid.nz * grid.ntimes);
+ THT_to_T(n, PABST, THT);
+
+ CHK(munmap(PABST, PABST_len) == 0);
+ CHK(munmap(THT, THT_len) == 0);
}
#undef WRITE
