diff --git a/.ci/bb/breeze-mpim/init.sh b/.ci/bb/breeze-mpim/init.sh index bc419103fb804d48a2dc0eeae1919eb90d12bc65..beb5dfdd6b2a3ebcafda4780fcaf1c9f45e1bb11 100644 --- a/.ci/bb/breeze-mpim/init.sh +++ b/.ci/bb/breeze-mpim/init.sh @@ -28,7 +28,7 @@ init_gcc () { init_env - switch_for_module gcc/12.1.0 mpich/3.4.3-gcc-12.1.0 + switch_for_module gcc/12.1.0 mpich/4.1.2-gcc-12.1.0 CC=gcc CXX=g++ @@ -38,9 +38,9 @@ init_gcc () MPI_LAUNCH="$(which mpirun)" ECCODES_ROOT='/sw/bullseye-x64/packages/gcc-12.1.0/eccodes-2.26.0' - NETCDF_ROOT='/sw/bullseye-x64/packages/gcc-12.1.0/netcdf-c-4.9.0' - PPM_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/ppm-1.0.8.1-gcc-12.1.0' - YAXT_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/yaxt-0.9.3.1-gcc-12.1.0' + NETCDF_ROOT='/sw/bullseye-x64/packages/gcc-12.1.0/netcdf-c-4.9.2' + PPM_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/ppm-1.0.8.2-mpich-4.1.2-gcc-12.1.0' + YAXT_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/yaxt-0.10.0-mpich-4.1.2-gcc-12.1.0' # The installations of NetCDF and ecCodes libraries do not provide '*.la' # files, which would trigger the RPATH injection: @@ -54,7 +54,7 @@ init_nvhpc () { init_env - switch_for_module nvhpc/22.3 mpich/3.4.3-nvhpc-22.3 + switch_for_module nvhpc/23.7 mpich/4.1.2-nvhpc-23.7 CC=nvc CXX=nvc++ @@ -63,10 +63,10 @@ init_nvhpc () MPIFC=mpif90 MPI_LAUNCH="$(which mpirun)" - ECCODES_ROOT='/sw/bullseye-x64/packages/nvhpc-22.3/eccodes-2.26.0' - NETCDF_ROOT='/sw/bullseye-x64/packages/nvhpc-22.3/netcdf-c-4.9.0' - PPM_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/ppm-1.0.8.1-nvhpc-22.3' - YAXT_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/yaxt-0.9.3.1-nvhpc-22.3' + ECCODES_ROOT='/sw/bullseye-x64/packages/nvhpc-23.7/eccodes-2.26.0' + NETCDF_ROOT='/sw/bullseye-x64/packages/nvhpc-23.7/netcdf-c-4.9.2' + PPM_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/ppm-1.0.8.2-mpich-4.1.2-nvhpc-23.7' + YAXT_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/yaxt-0.10.0-mpich-4.1.2-nvhpc-23.7' # The installations of NetCDF and ecCodes libraries do not provide '*.la' # files, which would trigger the RPATH injection: @@ -93,9 +93,9 @@ init_clang () MPI_LAUNCH="$(which mpirun) --oversubscribe" ECCODES_ROOT='/sw/bullseye-x64/packages/clang-14.0.6/eccodes-2.26.0' - NETCDF_ROOT='/sw/bullseye-x64/packages/clang-14.0.6/netcdf-c-4.9.0-openmpi-4.1.3' - PPM_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/ppm-1.0.8.1-clang-14.0.6' - YAXT_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/yaxt-0.9.3.1-clang-14.0.6' + NETCDF_ROOT='/sw/bullseye-x64/packages/clang-14.0.6/netcdf-c-4.9.2-openmpi-4.1.3' + PPM_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/ppm-1.0.8.2-openmpi-4.1.3-clang-14.0.6' + YAXT_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/yaxt-0.10.0-openmpi-4.1.3-clang-14.0.6' # The installations of NetCDF and ecCodes libraries do not provide '*.la' # files, which would trigger the RPATH injection: @@ -112,7 +112,7 @@ init_nag () init_env - switch_for_module nag/7.1 mpich/3.4.3-nag-7.1 + switch_for_module nag/7.1.7125 mpich/4.1.2-nag-7.1.7125 CC=gcc CXX=g++ @@ -121,10 +121,10 @@ init_nag () MPIFC=mpif90 MPI_LAUNCH="$(which mpirun)" - ECCODES_ROOT='/sw/bullseye-x64/packages/nag-7.1/eccodes-2.26.0' - NETCDF_ROOT='/sw/bullseye-x64/packages/nag-7.1/netcdf-c-4.9.0-mpich-3.4.3' - PPM_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/ppm-1.0.8.1-nag-7.1' - YAXT_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/yaxt-0.9.3.1-nag-7.1' + ECCODES_ROOT='/sw/bullseye-x64/packages/nag-7.1.7125/eccodes-2.26.0' + NETCDF_ROOT='/sw/bullseye-x64/packages/nag-7.1.7125/netcdf-c-4.9.2-mpich-4.1.2' + PPM_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/ppm-1.0.8.2-mpich-4.1.2-nag-7.1.7125' + YAXT_ROOT='/data/mpi/sclab/sip/m300488/libcdi-ci-sw/yaxt-0.10.0-mpich-4.1.2-nag-7.1.7125' # The installations of NetCDF and ecCodes libraries do not provide '*.la' # files, which would trigger the RPATH injection: diff --git a/.ci/bb/breeze-mpim/install_prerequisites.sh b/.ci/bb/breeze-mpim/install_prerequisites.sh index a6cc27b07f3ccb5113a8daf123c0237dc8dffe99..8b9fe12c0e326dfd90ee6fdb3340f63c0041d786 100755 --- a/.ci/bb/breeze-mpim/install_prerequisites.sh +++ b/.ci/bb/breeze-mpim/install_prerequisites.sh @@ -13,103 +13,107 @@ make_cmd='make -j22' mkdir -p "${work_dir}" && cd "${work_dir}" # Get PPM: -wget https://swprojects.dkrz.de/redmine/attachments/download/521/ppm-1.0.8.1.tar.gz -tar xvf ppm-1.0.8.1.tar.gz -ppm_src_dir="${work_dir}/ppm-1.0.8.1" -ppm_name_tag='ppm-1.0.8.1' +wget https://swprojects.dkrz.de/redmine/attachments/download/525/ppm-1.0.8.2.tar.gz +tar xvf ppm-1.0.8.2.tar.gz +ppm_src_dir="${work_dir}/ppm-1.0.8.2" +ppm_name_tag='ppm-1.0.8.2' ppm_config_args='--enable-MPI --disable-netcdf --disable-hdf5 --disable-parmetis --disable-metis --disable-crypto' # Get YAXT: -git clone -b release-0.9.3.1 https://gitlab.dkrz.de/dkrz-sw/yaxt.git -git -C yaxt cherry-pick 8c1b18bc9cfb3c7185017e37f8b39f7a61c94259 +git clone -b release-0.10.0 https://gitlab.dkrz.de/dkrz-sw/yaxt.git +git -C yaxt cherry-pick 602493aad8c6e817f32c9a4889fc2a271573f896 yaxt_src_dir="${work_dir}/yaxt" -yaxt_name_tag='yaxt-0.9.3.1' +yaxt_name_tag='yaxt-0.10.0' yaxt_config_args='' export CC='mpicc' export FC='mpif90' # Install for GCC 12.1.0: +module load mpich/4.1.2-gcc-12.1.0 +mpi_name_tag='mpich-4.1.2' compiler_name_tag='gcc-12.1.0' -module load mpich/3.4.3-gcc-12.1.0 # Install PPM: -build_dir="${ppm_name_tag}-${compiler_name_tag}" +build_dir="${ppm_name_tag}-${mpi_name_tag}-${compiler_name_tag}" mkdir "${build_dir}" ( cd "${build_dir}" - "${ppm_src_dir}/configure" ${ppm_config_args} --prefix="${install_dir}/${ppm_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -fallow-argument-mismatch' + "${ppm_src_dir}/configure" ${ppm_config_args} --prefix="${install_dir}/${ppm_name_tag}-${mpi_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -fallow-argument-mismatch' $make_cmd $make_cmd check $make_cmd install ) # Install YAXT: -build_dir="${yaxt_name_tag}-${compiler_name_tag}" +build_dir="${yaxt_name_tag}-${mpi_name_tag}-${compiler_name_tag}" mkdir "${build_dir}" ( cd "${build_dir}" - "${yaxt_src_dir}/configure" ${yaxt_config_args} --prefix="${install_dir}/${yaxt_name_tag}-${compiler_name_tag}" + "${yaxt_src_dir}/configure" ${yaxt_config_args} --prefix="${install_dir}/${yaxt_name_tag}-${mpi_name_tag}-${compiler_name_tag}" $make_cmd $make_cmd check $make_cmd install ) -module unload mpich/3.4.3-gcc-12.1.0 +module unload mpich/4.1.2-gcc-12.1.0 # Install for NVHPC 22.3: -compiler_name_tag='nvhpc-22.3' -module load mpich/3.4.3-nvhpc-22.3 +module load mpich/4.1.2-nvhpc-23.7 +mpi_name_tag='mpich-4.1.2' +compiler_name_tag='nvhpc-23.7' # Install PPM: -build_dir="${ppm_name_tag}-${compiler_name_tag}" +build_dir="${ppm_name_tag}-${mpi_name_tag}-${compiler_name_tag}" mkdir "${build_dir}" ( cd "${build_dir}" - "${ppm_src_dir}/configure" ${ppm_config_args} --prefix="${install_dir}/${ppm_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -tp sandybridge' CFLAGS='-g -O2 -tp sandybridge' + "${ppm_src_dir}/configure" ${ppm_config_args} --prefix="${install_dir}/${ppm_name_tag}-${mpi_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -tp sandybridge' CFLAGS='-g -O2 -tp sandybridge' $make_cmd $make_cmd check $make_cmd install ) # Install YAXT: -build_dir="${yaxt_name_tag}-${compiler_name_tag}" +build_dir="${yaxt_name_tag}-${mpi_name_tag}-${compiler_name_tag}" mkdir "${build_dir}" ( cd "${build_dir}" - "${yaxt_src_dir}/configure" ${yaxt_config_args} --prefix="${install_dir}/${yaxt_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -tp sandybridge' CFLAGS='-g -O2 -tp sandybridge' + "${yaxt_src_dir}/configure" ${yaxt_config_args} --prefix="${install_dir}/${yaxt_name_tag}-${mpi_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -tp sandybridge' CFLAGS='-g -O2 -tp sandybridge' $make_cmd $make_cmd check $make_cmd install ) -module unload mpich/3.4.3-nvhpc-22.3 +module unload mpich/4.1.2-nvhpc-23.7 # Install for Clang 14.0.6: -compiler_name_tag='clang-14.0.6' module load openmpi/4.1.3-clang-14.0.6 +mpi_name_tag='openmpi-4.1.3' +compiler_name_tag='clang-14.0.6' # Install PPM: -build_dir="${ppm_name_tag}-${compiler_name_tag}" +build_dir="${ppm_name_tag}-${mpi_name_tag}-${compiler_name_tag}" mkdir "${build_dir}" ( cd "${build_dir}" - "${ppm_src_dir}/configure" ${ppm_config_args} --prefix="${install_dir}/${ppm_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -no-pie' + "${ppm_src_dir}/configure" ${ppm_config_args} --prefix="${install_dir}/${ppm_name_tag}-${mpi_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -no-pie' $make_cmd $make_cmd check $make_cmd install ) # Install YAXT: -build_dir="${yaxt_name_tag}-${compiler_name_tag}" +build_dir="${yaxt_name_tag}-${mpi_name_tag}-${compiler_name_tag}" mkdir "${build_dir}" ( cd "${build_dir}" - "${yaxt_src_dir}/configure" ${yaxt_config_args} --prefix="${install_dir}/${yaxt_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -no-pie' + "${yaxt_src_dir}/configure" ${yaxt_config_args} --prefix="${install_dir}/${yaxt_name_tag}-${mpi_name_tag}-${compiler_name_tag}" FCFLAGS='-g -O2 -no-pie' $make_cmd $make_cmd check $make_cmd install ) module unload openmpi/4.1.3-clang-14.0.6 -# Install for NAG 7.1: -compiler_name_tag='nag-7.1' -module load mpich/3.4.3-nag-7.1 +# Install for NAG 7.1.7125: +module load mpich/4.1.2-nag-7.1.7125 +mpi_name_tag='mpich-4.1.2' +compiler_name_tag='nag-7.1.7125' # Install PPM: -build_dir="${ppm_name_tag}-${compiler_name_tag}" +build_dir="${ppm_name_tag}-${mpi_name_tag}-${compiler_name_tag}" mkdir "${build_dir}" ( cd "${build_dir}" - "${ppm_src_dir}/configure" ${ppm_config_args} --prefix="${install_dir}/${ppm_name_tag}-${compiler_name_tag}" FC=no + "${ppm_src_dir}/configure" ${ppm_config_args} --prefix="${install_dir}/${ppm_name_tag}-${mpi_name_tag}-${compiler_name_tag}" FC=no $make_cmd $make_cmd check $make_cmd install ) @@ -118,9 +122,9 @@ mkdir "${build_dir}" build_dir="${yaxt_name_tag}-${compiler_name_tag}" mkdir "${build_dir}" ( cd "${build_dir}" - "${yaxt_src_dir}/configure" ${yaxt_config_args} --prefix="${install_dir}/${yaxt_name_tag}-${compiler_name_tag}" + "${yaxt_src_dir}/configure" ${yaxt_config_args} --prefix="${install_dir}/${yaxt_name_tag}-${mpi_name_tag}-${compiler_name_tag}" $make_cmd $make_cmd check $make_cmd install ) -module unload mpich/3.4.3-nag-7.1 +module unload mpich/4.1.2-nag-7.1.7125 diff --git a/.gitignore b/.gitignore index a0acb6ced52d7f6ec77f0ac44576e9e6e7bea362..cdea57a8cb334212bdaf9692e9455363521629f5 100644 --- a/.gitignore +++ b/.gitignore @@ -42,6 +42,7 @@ Makefile.in Makefile # Build stage files: +*.L *.la *.lo *.mod diff --git a/ChangeLog b/ChangeLog index 7baea962161960f69e145d70cfbe8fc841d9e55a..651e6e7c53e6834a2846a61c3b0b4e24205df5b9 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,15 @@ +2023-10-18 Uwe Schulzweida + + * Version 2.3.0 released + +2023-10-06 Uwe Schulzweida + + * Add GRIB2 HEALPix support (available with eccodes-2.32.0) + +2023-09-29 Uwe Schulzweida + + * Add interface function streamDefShuffle() + 2023-08-15 Uwe Schulzweida * Version 2.2.4 released diff --git a/NEWS b/NEWS index 31ee891a61c5d1056d34609d64780238f6fe25a4..66f83c2a8f9091e081aa6159899bb738c87f91a9 100644 --- a/NEWS +++ b/NEWS @@ -1,9 +1,10 @@ CDI NEWS -------- -Version 2.3.0 (23 October 2023): +Version 2.3.0 (18 October 2023): New features: + * Add GRIB2 HEALPix support (available with eccodes-2.32.0) * Add support for NetCDF attribute type NC_INT64 * Add interface function streamInqNumSteps() to inquire number of time steps * gribapi decode: add support for single precision float interface (available since ecCodes-2.30.0) @@ -14,7 +15,7 @@ Version 2.3.0 (23 October 2023): * gribapiDefGridLCC: store DxInMetres/DyInMetres as double * NetCDF4: use chunkSize only if it is less than gridsize * cdfDefineAttributes: check filetype for unsigned int attributes - * time information missing if the stream contains only constant time fields + * time information missing if the stream contains fields constant in time only Version 2.2.0 (23 February 2023): diff --git a/config/default b/config/default index 589653e5711dc6e853c90380eb714e2ba9cef148..238a6978cb297034ac71e7e7fd4c22e3e79188d5 100755 --- a/config/default +++ b/config/default @@ -66,10 +66,10 @@ case "${HOSTNAME}" in CDILIBS="--disable-iso-c-interface \ --enable-maintainer-mode \ --with-szlib=/opt/local/lib/libaec \ - --with-eccodes=$HOME/local/eccodes-2.30.0 \ + --with-eccodes=$HOME/local/eccodes-2.32.0 \ --with-netcdf=$HOME/local/netcdf-c-4.9.1" PREFIX="--prefix=$HOME/local/cdi" - LD_ADD="-Wl,-rpath,$HOME/local/eccodes-2.30.0/lib" + LD_ADD="-Wl,-rpath,$HOME/local/eccodes-2.32.0/lib" if test "$COMP" = clang ; then ${CONFPATH}configure $CONFIG_OPTS $PREFIX $CDILIBS LDFLAGS="$LD_ADD $LDFLAGS" \ CC=clang CFLAGS="-g -pipe -D_REENTRANT -Wall -Wwrite-strings -W -Wfloat-equal -pedantic -O3" diff --git a/configure.ac b/configure.ac index 19426af95840d6d6a22f04608c632fe4d180ef27..3d986890f58e4071bd6af210709829821a666534 100644 --- a/configure.ac +++ b/configure.ac @@ -7,7 +7,7 @@ AC_PREREQ([2.69]) LT_PREREQ([2.4.6]) -AC_INIT([cdi],[2.2.4],[https://mpimet.mpg.de/cdi]) +AC_INIT([cdi],[2.3.0],[https://mpimet.mpg.de/cdi]) AC_DEFINE_UNQUOTED(CDI, ["$PACKAGE_VERSION"], [CDI version]) AC_CONFIG_AUX_DIR([config]) diff --git a/doc/tex/cdi_cman.tex b/doc/tex/cdi_cman.tex index 46b074f4edad3133c743d7ce31b6160484a256fd..308ae850222f8e911542cb8b9e9f4b9ab7db7e5a 100644 --- a/doc/tex/cdi_cman.tex +++ b/doc/tex/cdi_cman.tex @@ -136,7 +136,7 @@ \end{picture} \begin{flushright} -{\large\bfseries Climate Data Interface \\ Version 2.2.0 \\ March 2023} +{\large\bfseries Climate Data Interface \\ Version 2.3.0 \\ October 2023} \end{flushright} \vfill diff --git a/doc/tex/cdi_fman.tex b/doc/tex/cdi_fman.tex index 8bbcd5326c3fe2a3eaa2b0c78c1ac7e9e8834aa9..389fc4f2eb67ee961c7a432ccac0dee5b39511f7 100644 --- a/doc/tex/cdi_fman.tex +++ b/doc/tex/cdi_fman.tex @@ -133,7 +133,7 @@ \end{picture} \begin{flushright} -{\large\bfseries Climate Data Interface \\ Version 2.2.0 \\ March 2023} +{\large\bfseries Climate Data Interface \\ Version 2.3.0 \\ October 2023} \end{flushright} \vfill diff --git a/src/cdf_write.c b/src/cdf_write.c index 5bac6a3d9098042fdf9b7333720290a5d3ddba8f..cafd1f1f28ab8e060be398e45693b22e39722ac8 100644 --- a/src/cdf_write.c +++ b/src/cdf_write.c @@ -36,11 +36,12 @@ cdf_def_var_filter(int ncid, int ncvarID, unsigned int id, size_t nparams, const } void -cdfDefVarDeflate(int ncid, int ncvarID, int compLevel) +cdfDefVarDeflate(int ncid, int ncvarID, int shuffle, int compLevel) { #ifdef HAVE_NETCDF4 - // Set chunking, shuffle, and deflate. - int shuffle = (CDI_Shuffle > 0), deflate = 1; + int deflate = 1; + + if (CDI_Shuffle > 0 && shuffle == 0) shuffle = 1; if (compLevel < 1 || compLevel > 9) compLevel = 1; @@ -487,7 +488,7 @@ cdfDefVarCompression(const stream_t *streamptr, int ncvarID, nc_type xtype) if (streamptr->filetype == CDI_FILETYPE_NC4 || streamptr->filetype == CDI_FILETYPE_NC4C || streamptr->filetype == CDI_FILETYPE_NCZARR) { - cdfDefVarDeflate(streamptr->fileID, ncvarID, streamptr->complevel); + cdfDefVarDeflate(streamptr->fileID, ncvarID, streamptr->shuffle, streamptr->complevel); } else { diff --git a/src/cdi.h b/src/cdi.h index 257934438265e2cd82ff37e62701f62ad8149c6b..eed826b95864ee063d2c17525d1842a69ba6a9a0 100644 --- a/src/cdi.h +++ b/src/cdi.h @@ -387,6 +387,8 @@ void streamDefByteorder(int streamID, int byteorder); // streamInqByteorder: Get the byteorder int streamInqByteorder(int streamID); +void streamDefShuffle(int streamID, int shuffle); + void streamDefFilter(int streamID, int filterId, int nparams, const int *params); // streamDefCompType: Define compression type @@ -1405,6 +1407,8 @@ struct CDI_GridProjParams double y_0; // False northing (optional) double x_SP; // Longitude of southern pole double y_SP; // Latitude of southern pole + int nside; // HEALPix number of points along a side (number of data points should be = 12 * nside * nside) + int order; // HEALPix ordering convention (0:ring; 1:nested) }; void gridProjParamsInit(struct CDI_GridProjParams *gridProjParams); @@ -1417,6 +1421,10 @@ int gridInqParamsLCC(int gridID, struct CDI_GridProjParams *gridProjParams); void gridDefParamsSTERE(int gridID, struct CDI_GridProjParams gridProjParams); int gridInqParamsSTERE(int gridID, struct CDI_GridProjParams *gridProjParams); +// HEALPix grid +void gridDefParamsHEALPIX(int gridID, struct CDI_GridProjParams gridProjParams); +int gridInqParamsHEALPIX(int gridID, struct CDI_GridProjParams *gridProjParams); + #define HAVE_CDI_PROJ_FUNCS 1 extern int (*proj_lonlat_to_lcc_func)(struct CDI_GridProjParams gpp, size_t, double*, double*); extern int (*proj_lcc_to_lonlat_func)(struct CDI_GridProjParams gpp, double, double, size_t, double*, double*); diff --git a/src/cdi_int.h b/src/cdi_int.h index 7612e7dd46be958768a02a24990196cb674e2dfa..c9ba92cc406cd2b29f37cbc3be8ea20519b87e21 100644 --- a/src/cdi_int.h +++ b/src/cdi_int.h @@ -301,6 +301,7 @@ typedef struct int localatts; int unreduced; int have_missval; + int shuffle; // netcdf4/HDF5 filter unsigned int filterId; size_t numParams; diff --git a/src/grb_read.c b/src/grb_read.c index 771ad2618375fe241e652ad97422949d73508838..06d6c3458baf6e85b07eeda3ae0bccbbb1f017a9 100644 --- a/src/grb_read.c +++ b/src/grb_read.c @@ -37,15 +37,17 @@ grb_decode(int filetype, int memType, int datatype, void *cgribexp, void *gribbu #endif #ifdef HAVE_LIBGRIB_API { - void *datap = data; bool useFloatInterface = (have_gribapi_float_interface() && datatype != CDI_DATATYPE_FLT32 && datatype != CDI_DATATYPE_FLT64); int memTypeX = useFloatInterface ? memType : MEMTYPE_DOUBLE; - if (!useFloatInterface && memType == MEMTYPE_FLOAT) datap = Malloc(datasize * sizeof(double)); + void *datap = (!useFloatInterface && memType == MEMTYPE_FLOAT) ? Malloc(datasize * sizeof(double)) : data; + + // if (useFloatInterface) printf("gribapi read: useFloatInterface\n"); status = gribapiDecode(memTypeX, gribbuffer, gribsize, datap, datasize, unreduced, nmiss, missval); if (!useFloatInterface && memType == MEMTYPE_FLOAT) { + // printf("gribapi read: convert double to float\n"); float *dataf = (float *) data; double *datad = (double *) datap; for (size_t i = 0; i < datasize; ++i) dataf[i] = (float) datad[i]; diff --git a/src/grb_write.c b/src/grb_write.c index 818ac30f262fb22317a902e3fb3b93c07e3bdc41..e71004aa00c719a210c39a4cadbc6d04faf18481 100644 --- a/src/grb_write.c +++ b/src/grb_write.c @@ -24,9 +24,9 @@ #include "namespace.h" static size_t -grb_encode(int filetype, int memType, int varID, int levelID, int vlistID, int gridID, int zaxisID, CdiDateTime vDateTime, - int tsteptype, int numavg, size_t datasize, const void *data, size_t nmiss, void **gribbuffer, int comptype, - void *gribContainers) +grb_encode(int filetype, int memType, int datatype, int varID, int levelID, int vlistID, int gridID, int zaxisID, + CdiDateTime vDateTime, int tsteptype, int numavg, size_t datasize, const void *data, size_t nmiss, void **gribbuffer, + int comptype, void *gribContainers) { size_t nbytes = 0; @@ -48,10 +48,17 @@ grb_encode(int filetype, int memType, int varID, int levelID, int vlistID, int g #else void *gribContainer = (void *) &((gribContainer_t *) gribContainers)[varID]; #endif + // bool useFloatInterface = (have_gribapi_float_interface() && datatype != CDI_DATATYPE_FLT32 && datatype != + // CDI_DATATYPE_FLT64); + bool useFloatInterface = false; + int memTypeX = useFloatInterface ? memType : MEMTYPE_DOUBLE; const void *datap = data; - int memTypeX = /*have_gribapi_float_interface() ? memType :*/ MEMTYPE_DOUBLE; - if (/*!have_gribapi_float_interface() &&*/ memType == MEMTYPE_FLOAT) + + // if (useFloatInterface) printf("gribapi write: useFloatInterface\n"); + + if (!useFloatInterface && memType == MEMTYPE_FLOAT) { + // printf("gribapi write: convert float to double\n"); const float *dataf = (const float *) data; double *datad = (double *) Malloc(datasize * sizeof(double)); for (size_t i = 0; i < datasize; ++i) datad[i] = (double) dataf[i]; @@ -62,7 +69,7 @@ grb_encode(int filetype, int memType, int varID, int levelID, int vlistID, int g nbytes = gribapiEncode(memTypeX, varID, levelID, vlistID, gridID, zaxisID, vDateTime, tsteptype, numavg, (long) datasize, datap, nmiss, gribbuffer, &gribbuffersize, comptype, gribContainer); - if (/*!have_gribapi_float_interface() &&*/ memType == MEMTYPE_FLOAT) Free((void *) datap); + if (!useFloatInterface && memType == MEMTYPE_FLOAT) Free((void *) datap); } #else { @@ -379,6 +386,7 @@ grb_write_var_slice(stream_t *streamptr, int varID, int levelID, int memtype, co CdiDateTime vDateTime = streamptr->tsteps[tsID].taxis.vDateTime; int numavg = (tsteptype == TSTEP_AVG) ? streamptr->tsteps[tsID].taxis.numavg : 0; int comptype = streamptr->comptype; + int datatype = vlistInqVarDatatype(vlistID, varID); if (CDI_Debug) Message("gridID = %d zaxisID = %d", gridID, zaxisID); @@ -394,8 +402,8 @@ grb_write_var_slice(stream_t *streamptr, int varID, int levelID, int memtype, co } void *gribbuffer = NULL; - size_t nbytes = grb_encode(filetype, memtype, varID, levelID, vlistID, gridID, zaxisID, vDateTime, tsteptype, numavg, datasize, - data, nmiss, &gribbuffer, comptype, streamptr->gribContainers); + size_t nbytes = grb_encode(filetype, memtype, datatype, varID, levelID, vlistID, gridID, zaxisID, vDateTime, tsteptype, numavg, + datasize, data, nmiss, &gribbuffer, comptype, streamptr->gribContainers); if (filetype == CDI_FILETYPE_GRB && (comptype == CDI_COMPRESS_SZIP || comptype == CDI_COMPRESS_AEC)) nbytes = grbSzip(filetype, gribbuffer, nbytes); diff --git a/src/gribapi.h b/src/gribapi.h index 089f866c943cc80cf3139c0abe7b1bd83d657c53..87e9a23e30fa062572fc3314a9e18d5a71d8fdf7 100644 --- a/src/gribapi.h +++ b/src/gribapi.h @@ -56,6 +56,7 @@ #define GRIB2_GTYPE_SPECTRAL 50 // Spherical harmonic coefficients #define GRIB2_GTYPE_GME 100 // Triangular grid based on an icosahedron (GME) #define GRIB2_GTYPE_UNSTRUCTURED 101 // General Unstructured Grid +#define GRIB2_GTYPE_HEALPIX 150 // HEALPix Grid const char *gribapiLibraryVersionString(void); void gribContainersNew(stream_t *streamptr); diff --git a/src/gribapi_utilities.c b/src/gribapi_utilities.c index 56a44e3de52e18a148cc9f4ae56d6dbb67d03450..ca9f8e090a75eaebbf13dba0b18ab58622d0fd1c 100644 --- a/src/gribapi_utilities.c +++ b/src/gribapi_utilities.c @@ -584,6 +584,7 @@ gribapiGetGridType(grib_handle *gh) case GRIB2_GTYPE_SPECTRAL: return GRID_SPECTRAL; case GRIB2_GTYPE_GME: return GRID_GME; case GRIB2_GTYPE_UNSTRUCTURED: return GRID_UNSTRUCTURED; + case GRIB2_GTYPE_HEALPIX: return CDI_PROJ_HEALPIX; default: { static bool lwarn = true; @@ -780,6 +781,12 @@ gribapiGetGridProj(grib_handle *gh, grid_t *grid, size_t numberOfPoints) grid->y.flag = 2; } +static void +gribapiGetGridHealpix(grib_handle *gh, grid_t *grid, size_t numberOfPoints) +{ + grid->size = numberOfPoints; +} + static void gribapiGetGridSpectral(grib_handle *gh, grid_t *grid, size_t datasize) { @@ -859,7 +866,7 @@ gribapiGetGrid(grib_handle *gh, grid_t *grid) long editionNumber = gribEditionNumber(gh); int gridtype = gribapiGetGridType(gh); int projtype = (gridtype == GRID_PROJECTION && gribapiGetIsRotated(gh)) ? CDI_PROJ_RLL : CDI_UNDEFID; - if (gridtype == CDI_PROJ_LCC || gridtype == CDI_PROJ_STERE) + if (gridtype == CDI_PROJ_LCC || gridtype == CDI_PROJ_STERE || gridtype == CDI_PROJ_HEALPIX) { projtype = gridtype; gridtype = GRID_PROJECTION; @@ -897,6 +904,10 @@ gribapiGetGrid(grib_handle *gh, grid_t *grid) { gribapiGetGridProj(gh, grid, numberOfPoints); } + else if (projtype == CDI_PROJ_HEALPIX) + { + gribapiGetGridHealpix(gh, grid, numberOfPoints); + } else if (gridtype == GRID_SPECTRAL) { gribapiGetGridSpectral(gh, grid, datasize); @@ -918,7 +929,8 @@ gribapiGetGrid(grib_handle *gh, grid_t *grid) Error("Unsupported grid type: %s", gridNamePtr(gridtype)); } - if (gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT || projtype == CDI_PROJ_RLL || projtype == CDI_PROJ_LCC) + if (gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT || projtype == CDI_PROJ_RLL || projtype == CDI_PROJ_LCC + || projtype == CDI_PROJ_HEALPIX) { long temp = 0; GRIB_CHECK(grib_get_long(gh, "uvRelativeToGrid", &temp), 0); @@ -926,7 +938,7 @@ gribapiGetGrid(grib_handle *gh, grid_t *grid) uvRelativeToGrid = (bool) temp; } - if (gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT || gridtype == GRID_PROJECTION) + if (gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT || (gridtype == GRID_PROJECTION && projtype != CDI_PROJ_HEALPIX)) { long iScansNegatively, jScansPositively, jPointsAreConsecutive; GRIB_CHECK(grib_get_long(gh, "iScansNegatively", &iScansNegatively), 0); diff --git a/src/grid.c b/src/grid.c index 1e6324a830d78cbc0ffb300f8c34fbbe4ba86948..9d1eb1bff05ed25efad77ffc2598d409f48ea274 100644 --- a/src/grid.c +++ b/src/grid.c @@ -3656,6 +3656,8 @@ gridProjParamsInit(struct CDI_GridProjParams *gpp) gpp->y_0 = CDI_Grid_Missval; // False northing (optional) gpp->x_SP = CDI_Grid_Missval; // Longitude of southern pole gpp->y_SP = CDI_Grid_Missval; // Latitude of southern pole + gpp->nside = 0; // HEALPix number of points along a side (number of data points should be = 12 * nside * nside) + gpp->order = -1; // HEALPix ordering convention (0:ring; 1:nested) // clang-format on } @@ -3859,6 +3861,24 @@ gridVerifyProjParamsSTERE(struct CDI_GridProjParams *gpp) return 0; } +int +gridVerifyProjParamsHEALPIX(struct CDI_GridProjParams *gpp) +{ + static bool lwarn = true; + + if (lwarn) + { + lwarn = false; + const char *projection = "healpix"; + if (IS_EQUAL(gpp->nside, -1)) Error("%s mapping parameter %s missing!", projection, "nside"); + if (IS_EQUAL(gpp->order, -1)) Error("%s mapping parameter %s missing!", projection, "order"); + if (gpp->nside == 0) Error("%s mapping parameter %s unsupported!", projection, "nside", gpp->nside); + if (gpp->order != 0 && gpp->order != 1) Error("%s mapping parameter %s=%d unsupported!", projection, "order", gpp->order); + } + + return 0; +} + /* @Function gridDefParamsSTERE @Title Define the parameter of a Polar stereographic grid @@ -3892,6 +3912,26 @@ gridDefParamsSTERE(int gridID, struct CDI_GridProjParams gpp) gridVerifyProj(gridID); } +void +gridDefParamsHEALPIX(int gridID, struct CDI_GridProjParams gpp) +{ + cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_VARNAME, "healpix"); + + const char *gmapname = "healpix"; + cdiDefKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname); + cdiDefAttTxt(gridID, CDI_GLOBAL, "grid_mapping_name", (int) (strlen(gmapname)), gmapname); + + cdiDefAttInt(gridID, CDI_GLOBAL, "healpix_nside", CDI_DATATYPE_INT32, 1, &gpp.nside); + const char *orderName = (gpp.order == 1) ? "nested" : "ring"; + cdiDefAttTxt(gridID, CDI_GLOBAL, "healpix_order", (int) (strlen(orderName)), orderName); + + // gridDefParamsCommon(gridID, gpp); + + grid_t *gridptr = grid_to_pointer(gridID); + gridptr->projtype = CDI_PROJ_HEALPIX; + + // gridVerifyProj(gridID); +} /* @Function gridInqParamsSTERE @@ -3958,6 +3998,64 @@ gridInqParamsSTERE(int gridID, struct CDI_GridProjParams *gpp) return status; } +int +gridInqParamsHEALPIX(int gridID, struct CDI_GridProjParams *gpp) +{ + int status = -1; + if (gridInqType(gridID) != GRID_PROJECTION) return status; + + gridProjParamsInit(gpp); + + status = -2; + const char *projection = "healpix"; + char gmapname[CDI_MAX_NAME]; + int length = CDI_MAX_NAME; + cdiInqKeyString(gridID, CDI_GLOBAL, CDI_KEY_GRIDMAP_NAME, gmapname, &length); + if (gmapname[0] && str_is_equal(gmapname, projection)) + { + int atttype, attlen; + char attname[CDI_MAX_NAME + 1]; + + int natts; + cdiInqNatts(gridID, CDI_GLOBAL, &natts); + + if (natts) status = 0; + + for (int iatt = 0; iatt < natts; ++iatt) + { + cdiInqAtt(gridID, CDI_GLOBAL, iatt, attname, &atttype, &attlen); + + if (atttype == CDI_DATATYPE_TXT) + { + char attstring[256]; + if (cdiInqAttTxt(gridID, CDI_GLOBAL, attname, (int) sizeof(attstring), attstring) == 0) + { + attstring[attlen] = 0; + if (str_is_equal(attname, "healpix_order")) gpp->order = strStartsWith(attstring, "nest"); + } + } + else + { + if (attlen > 2) continue; + double attflt[2]; + if (cdiInqAttConvertedToFloat(gridID, atttype, attname, attlen, attflt)) + { + // clang-format off + if (str_is_equal(attname, "earth_radius")) gpp->a = attflt[0]; + else if (str_is_equal(attname, "semi_major_axis")) gpp->a = attflt[0]; + else if (str_is_equal(attname, "semi_minor_axis")) gpp->b = attflt[0]; + else if (str_is_equal(attname, "inverse_flattening")) gpp->rf = attflt[0]; + else if (str_is_equal(attname, "longitudeOfFirstGridPointInDegrees")) gpp->xval_0 = attflt[0]; + else if (str_is_equal(attname, "healpix_nside")) gpp->nside = (int) lround(attflt[0]); + // clang-format on + } + } + } + } + + return status; +} + void gridDefComplexPacking(int gridID, int lcomplex) { diff --git a/src/grid.h b/src/grid.h index 8bcfee3d299f1da772801b796be4e270c8341a1e..03b034d141c968d355afd690a56be9f3ecad40c0 100644 --- a/src/grid.h +++ b/src/grid.h @@ -181,6 +181,7 @@ struct addIfNewRes cdiVlistAddGridIfNew(int vlistID, grid_t *grid, int mode); int gridVerifyProjParamsLCC(struct CDI_GridProjParams *gpp); int gridVerifyProjParamsSTERE(struct CDI_GridProjParams *gpp); +int gridVerifyProjParamsHEALPIX(struct CDI_GridProjParams *gpp); bool isGaussianLatitudes(size_t nlats, const double *latitudes); void gaussianLatitudes(size_t nlats, double *latitudes, double *weights); diff --git a/src/julian_date.c b/src/julian_date.c index 0a9f7c92d49a2becd760a7815f40e123f530b64f..4d2c709e9e138f43bdaf90dae817baafd3efdd0e 100644 --- a/src/julian_date.c +++ b/src/julian_date.c @@ -8,17 +8,17 @@ decode_julday(int calendar, int64_t julianDay, // Julian day number to convert int *mon, // Gregorian month (1-12) (out) int *day) // Gregorian day (1-31) (out) { - const int64_t a = julianDay; + int64_t a = julianDay; - const double b = floor((a - 1867216.25) / 36524.25); + double b = floor((a - 1867216.25) / 36524.25); double c = a + b - floor(b / 4) + 1525; if (calendar == CALENDAR_STANDARD || calendar == CALENDAR_GREGORIAN) if (a < 2299161) c = a + 1524; - const double d = floor((c - 122.1) / 365.25); - const double e = floor(365.25 * d); - const double f = floor((c - e) / 30.6001); + double d = floor((c - 122.1) / 365.25); + double e = floor(365.25 * d); + double f = floor((c - e) / 30.6001); *day = (int) (c - e - floor(30.6001 * f)); *mon = (int) (f - 1 - 12 * floor(f / 14)); @@ -29,8 +29,8 @@ decode_julday(int calendar, int64_t julianDay, // Julian day number to convert static int64_t encode_julday(int calendar, int year, int month, int day) { - const int iy = (month <= 2) ? year - 1 : year; - const int im = (month <= 2) ? month + 12 : month; + int iy = (month <= 2) ? year - 1 : year; + int im = (month <= 2) ? month + 12 : month; int ib = (iy < 0) ? ((iy + 1) / 400 - (iy + 1) / 100) : (iy / 400 - iy / 100); if (calendar == CALENDAR_STANDARD || calendar == CALENDAR_GREGORIAN) @@ -75,7 +75,7 @@ time_to_sec(int time) int hour, minute, second; cdiDecodeTime(time, &hour, &minute, &second); - const int seconds = hour * 3600 + minute * 60 + second; + int seconds = hour * 3600 + minute * 60 + second; return seconds; } @@ -83,9 +83,9 @@ time_to_sec(int time) int sec_to_time(int secofday) { - const int hour = secofday / 3600; - const int minute = secofday / 60 - hour * 60; - const int second = secofday - hour * 3600 - minute * 60; + int hour = secofday / 3600; + int minute = secofday / 60 - hour * 60; + int second = secofday - hour * 3600 - minute * 60; return cdiEncodeTime(hour, minute, second); } @@ -93,9 +93,9 @@ sec_to_time(int secofday) double secofday_encode(CdiTime time) { - const int hour = time.hour; - const int minute = time.minute; - const int second = time.second; + int hour = time.hour; + int minute = time.minute; + int second = time.second; return hour * 3600 + minute * 60 + second + time.ms / 1000.0; } @@ -107,11 +107,11 @@ secofday_decode(double secondOfDay) double secondOfDayIntegral; time.ms = lround(modf(secondOfDay, &secondOfDayIntegral) * 1000); - const int fullSeconds = lrint(secondOfDayIntegral); + int fullSeconds = lrint(secondOfDayIntegral); - const int hour = fullSeconds / 3600; - const int minute = fullSeconds / 60 - hour * 60; - const int second = fullSeconds - hour * 3600 - minute * 60; + int hour = fullSeconds / 3600; + int minute = fullSeconds / 60 - hour * 60; + int second = fullSeconds - hour * 3600 - minute * 60; time.hour = hour; time.minute = minute; @@ -121,9 +121,9 @@ secofday_decode(double secondOfDay) } static int64_t -calDay_encode(int calendar, CdiDate date) +calendarDay_encode(int calendar, CdiDate date) { - const int dpy = calendar_dpy(calendar); + int dpy = calendar_dpy(calendar); if (dpy == 360 || dpy == 365 || dpy == 366) return encode_calday(dpy, date.year, date.month, date.day); @@ -132,10 +132,10 @@ calDay_encode(int calendar, CdiDate date) } static CdiDate -calDay_decode(int calendar, int64_t julday) +calendarDay_decode(int calendar, int64_t julday) { int year, month, day; - const int dpy = calendar_dpy(calendar); + int dpy = calendar_dpy(calendar); if (dpy == 360 || dpy == 365 || dpy == 366) decode_calday(dpy, julday, &year, &month, &day); @@ -155,7 +155,7 @@ julianDate_encode(int calendar, CdiDateTime dt) { JulianDate julianDate; - julianDate.julianDay = calDay_encode(calendar, dt.date); + julianDate.julianDay = calendarDay_encode(calendar, dt.date); julianDate.secondOfDay = secofday_encode(dt.time); return julianDate; @@ -166,7 +166,7 @@ julianDate_decode(int calendar, JulianDate julianDate) { CdiDateTime dt; - dt.date = calDay_decode(calendar, julianDate.julianDay); + dt.date = calendarDay_decode(calendar, julianDate.julianDay); dt.time = secofday_decode(julianDate.secondOfDay); return dt; @@ -175,7 +175,7 @@ julianDate_decode(int calendar, JulianDate julianDate) static void adjust_seconds(JulianDate *julianDate) { - const double SecondsPerDay = 86400.0; + double SecondsPerDay = 86400.0; while (julianDate->secondOfDay >= SecondsPerDay) { diff --git a/src/make_fint.c b/src/make_fint.c index 10eda72d61c4121003f37caac9af4bd9987d216d..235ac44536262fa541498aa0b7c178a4ac822151 100644 --- a/src/make_fint.c +++ b/src/make_fint.c @@ -9,7 +9,7 @@ #define _GNU_SOURCES 1 #endif #else -#define VERSION "2.2.0" +#define VERSION "2.3.0" #endif #ifndef _XOPEN_SOURCE diff --git a/src/stream.c b/src/stream.c index 09d733cb5d510d3cf98b14b14247cc1d2a0333f9..e302be77389f23fd8ce12e4132dfbb738426ae94 100644 --- a/src/stream.c +++ b/src/stream.c @@ -1007,6 +1007,7 @@ streamDefaultValue(stream_t *streamptr) streamptr->have_missval = cdiHaveMissval; streamptr->comptype = CDI_COMPRESS_NONE; streamptr->complevel = 0; + streamptr->shuffle = 0; // netcdf4/HDF5 filter streamptr->filterId = 0; streamptr->numParams = 0; @@ -1201,30 +1202,34 @@ streamDestroy(stream_t *streamptr) void (*streamCloseDelegate)(stream_t * streamptr, int recordBufIsToBeDeleted) = (void (*)(stream_t *, int)) namespaceSwitchGet(NSSWITCH_STREAM_CLOSE_BACKEND).func; - if (streamptr->filetype != -1) streamCloseDelegate(streamptr, 1); + if (streamptr->filetype != CDI_FILETYPE_UNDEF) streamCloseDelegate(streamptr, 1); if (streamptr->record) { if (streamptr->record->buffer) Free(streamptr->record->buffer); Free(streamptr->record); + streamptr->record = NULL; } - streamptr->filetype = 0; + streamptr->filetype = CDI_FILETYPE_UNDEF; if (streamptr->filename) { Free(streamptr->filename); streamptr->filename = NULL; } - for (int index = 0; index < streamptr->nvars; index++) + if (streamptr->vars) { - sleveltable_t *pslev = streamptr->vars[index].recordTable; - unsigned nsub = streamptr->vars[index].subtypeSize >= 0 ? (unsigned) streamptr->vars[index].subtypeSize : 0U; - for (size_t isub = 0; isub < nsub; isub++) deallocate_sleveltable_t(pslev + isub); - if (pslev) Free(pslev); + for (int index = 0; index < streamptr->nvars; index++) + { + sleveltable_t *pslev = streamptr->vars[index].recordTable; + unsigned nsub = streamptr->vars[index].subtypeSize >= 0 ? (unsigned) streamptr->vars[index].subtypeSize : 0U; + for (size_t isub = 0; isub < nsub; isub++) deallocate_sleveltable_t(pslev + isub); + if (pslev) Free(pslev); + } + Free(streamptr->vars); + streamptr->vars = NULL; } - Free(streamptr->vars); - streamptr->vars = NULL; if (streamptr->tsteps) { @@ -1241,11 +1246,13 @@ streamDestroy(stream_t *streamptr) } Free(streamptr->tsteps); + streamptr->tsteps = NULL; } if (vlistID != -1) { - if (streamptr->filemode != 'w' && vlistInqTaxis(vlistID) != -1) taxisDestroy(vlistInqTaxis(vlistID)); + int taxisID = (streamptr->filemode != 'w') ? vlistInqTaxis(vlistID) : -1; + if (taxisID != -1) taxisDestroy(taxisID); void (*mycdiVlistDestroy_)(int, bool) = (void (*)(int, bool)) namespaceSwitchGet(NSSWITCH_VLIST_DESTROY_).func; mycdiVlistDestroy_(vlistID, true); } @@ -1601,6 +1608,17 @@ streamInqVlist(int streamID) return s->vlistID; } +void +streamDefShuffle(int streamID, int shuffle) +{ + stream_t *s = stream_to_pointer(streamID); + if (s->shuffle != shuffle) + { + s->shuffle = shuffle; + reshSetStatus(streamID, &streamOps, RESH_DESYNC_IN_USE); + } +} + void streamDefFilter(int streamID, int filterId, int numParams, const int *params) { diff --git a/src/stream_cdf.h b/src/stream_cdf.h index 25576f3d0e74ab9c52da5e8d1b64c827af7f0b67..088348664ac6e41182ce8b0571de008e1fe389e2 100644 --- a/src/stream_cdf.h +++ b/src/stream_cdf.h @@ -47,7 +47,7 @@ void cdf_write_var_slice(stream_t *streamptr, int varID, int levelID, int memtyp void cdf_write_var_chunk(stream_t *streamptr, int varID, int memtype, const int rect[][2], const void *data, size_t nmiss); -void cdfDefVarDeflate(int ncid, int ncvarid, int deflateLevel); +void cdfDefVarDeflate(int ncid, int ncvarid, int shuffle, int deflateLevel); void cdfDefTime(stream_t *streamptr); void cdf_scale_add(size_t size, double *data, double addoffset, double scalefactor); diff --git a/src/stream_cdf_i.c b/src/stream_cdf_i.c index 229bb59bc1f26a3889c630c26abaa82a12a7931b..b1ef6f1859a3a4eea8fefb1f5493f861167207e6 100644 --- a/src/stream_cdf_i.c +++ b/src/stream_cdf_i.c @@ -33,6 +33,8 @@ enum AxisType T_AXIS = 5, }; +static int axisTypeChar[] = { '?', 'X', 'Y', 'Z', 'E', 'T' }; + typedef struct { int dimid; @@ -689,8 +691,8 @@ cdf_set_dim(ncvar_t *ncvar, int dimid, int dimtype) { if (ncvar->dimtype[dimid] != CDI_UNDEFID && ncvar->dimtype[dimid] != dimtype) { - Warning("Inconsistent dimension definition for %s! dimid = %d; type = %d; newtype = %d", ncvar->name, dimid, - ncvar->dimtype[dimid], dimtype); + Warning("Inconsistent dimension definition for %s! dimid=%d type=%c newtype=%c", ncvar->name, dimid, + axisTypeChar[ncvar->dimtype[dimid]], axisTypeChar[dimtype]); } ncvar->dimtype[dimid] = dimtype; @@ -3650,9 +3652,10 @@ cdf_define_all_vars(stream_t *streamptr, int vlistID, int instID, int modelID, i */ if (ncvar->numberOfForecastsInEnsemble != -1) { - cdiDefKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, ncvar->typeOfEnsembleForecast); cdiDefKeyInt(vlistID, varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, ncvar->numberOfForecastsInEnsemble); cdiDefKeyInt(vlistID, varID, CDI_KEY_PERTURBATIONNUMBER, ncvar->perturbationNumber); + if (ncvar->numberOfForecastsInEnsemble != -1) + cdiDefKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, ncvar->typeOfEnsembleForecast); } if (ncvar->extra[0] != 0) cdiDefKeyString(vlistID, varID, CDI_KEY_CHUNKS, ncvar->extra); diff --git a/src/stream_cdf_o.c b/src/stream_cdf_o.c index 82d1aa57c2b5efaee53e58619a8426c48b9f24e0..19e2cfb99c931be870870c8c36a1d646828de862 100644 --- a/src/stream_cdf_o.c +++ b/src/stream_cdf_o.c @@ -544,7 +544,8 @@ cdfGridCompress(int fileID, int ncvarid, size_t gridsize, int filetype, int comp && (filetype == CDI_FILETYPE_NC4 || filetype == CDI_FILETYPE_NC4C || filetype == CDI_FILETYPE_NCZARR)) { cdf_def_var_chunking(fileID, ncvarid, NC_CHUNKED, chunks); - cdfDefVarDeflate(fileID, ncvarid, 1); + int shuffle = 1, compLevel = 1; + cdfDefVarDeflate(fileID, ncvarid, shuffle, compLevel); } #endif } diff --git a/src/stream_cgribex.c b/src/stream_cgribex.c index b06652536ce9c272cc5db14f69f51d08b6461314..932a96eb02442c1a32200feb35d8b8e4db125e36 100644 --- a/src/stream_cgribex.c +++ b/src/stream_cgribex.c @@ -1850,6 +1850,11 @@ cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridID, in ISEC1_Sec2Or3Flag = 0; break; } + case CDI_PROJ_HEALPIX: + { + Error("CGRIBEX library doesn't support HEALPix grids!"); + break; + } default: { static bool lwarn = true; @@ -2023,14 +2028,14 @@ cgribexDefEnsembleVar(int *isec1, int vlistID, int varID) // Put2Byte(isec1[38]); // individual ensemble member // Put2Byte(isec1[39]); // number of forecasts in ensemble - int perturbationNumber, numberOfForecastsInEnsemble, typeOfEnsembleForecast; - int r1 = cdiInqKeyInt(vlistID, varID, CDI_KEY_PERTURBATIONNUMBER, &perturbationNumber); - int r2 = cdiInqKeyInt(vlistID, varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, &numberOfForecastsInEnsemble); - int r3 = cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, &typeOfEnsembleForecast); - - if (r1 == 0 && r2 == 0 && r3 == 0) + if (ISEC1_CenterID == 252) { - if (ISEC1_CenterID == 252) + int perturbationNumber, numberOfForecastsInEnsemble, typeOfEnsembleForecast; + int r1 = cdiInqKeyInt(vlistID, varID, CDI_KEY_PERTURBATIONNUMBER, &perturbationNumber); + int r2 = cdiInqKeyInt(vlistID, varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, &numberOfForecastsInEnsemble); + int r3 = cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, &typeOfEnsembleForecast); + + if (r1 == 0 && r2 == 0 && r3 == 0) { ISEC1_LocalFLag = 1; isec1[36] = 1; diff --git a/src/stream_grb.c b/src/stream_grb.c index 09d664f5cdead9985a670786afd0d29cd3376bf8..3ec32b16d8dae74783f8c5f92c0cfa6f9b57149e 100644 --- a/src/stream_grb.c +++ b/src/stream_grb.c @@ -499,6 +499,10 @@ grbGetGridtype(int *gridID, size_t gridsize, bool *gridIsRotated, bool *gridIsCu { gridtype = CDI_PROJ_STERE; } + else if (projtype == CDI_PROJ_HEALPIX) + { + gridtype = CDI_PROJ_HEALPIX; + } } return gridtype; diff --git a/src/stream_gribapi.c b/src/stream_gribapi.c index 80ec36f7a4b73a13e66b59b8235dc9cf201aa1d7..dafaf62e08fc9460eafc290ea86e451d19d308f1 100644 --- a/src/stream_gribapi.c +++ b/src/stream_gribapi.c @@ -40,13 +40,13 @@ typedef struct } compvar2_t; static int -gribapiGetZaxisType(long editionNumber, int grib_ltype) +gribapi_get_zaxis_type(long editionNumber, int grib_ltype) { return (editionNumber <= 1) ? grib1ltypeToZaxisType(grib_ltype) : grib2ltypeToZaxisType(grib_ltype); } static int -getTimeunits(long unitsOfTime) +get_timeunits(long unitsOfTime) { switch (unitsOfTime) { @@ -63,19 +63,19 @@ getTimeunits(long unitsOfTime) } static double -timeunit_factor(int tu1, int tu2) +timeunits_factor(int timeUnits1, int timeUnits2) { - if (tu2 == TUNIT_HOUR) + if (timeUnits2 == TUNIT_HOUR) { - switch (tu1) + switch (timeUnits1) { case TUNIT_SECOND: return 3600; case TUNIT_MINUTE: return 60; case TUNIT_HOUR: return 1; - case TUNIT_3HOURS: return 1. / 3; - case TUNIT_6HOURS: return 1. / 6; - case TUNIT_12HOURS: return 1. / 12; - case TUNIT_DAY: return 1. / 24; + case TUNIT_3HOURS: return 1.0 / 3; + case TUNIT_6HOURS: return 1.0 / 6; + case TUNIT_12HOURS: return 1.0 / 12; + case TUNIT_DAY: return 1.0 / 24; } } @@ -83,14 +83,14 @@ timeunit_factor(int tu1, int tu2) } static int -gribapiGetTimeUnits(grib_handle *gh) +gribapi_get_timeunits(grib_handle *gh) { long unitsOfTime = -1; grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime); GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0); - return getTimeunits(unitsOfTime); + return get_timeunits(unitsOfTime); } static void @@ -98,8 +98,8 @@ gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endStep) { long unitsOfTime; int status = grib_get_long(gh, "stepUnits", &unitsOfTime); - int timeunits2 = (status == 0) ? getTimeunits(unitsOfTime) : timeunits; - // timeunits2 = gribapiGetTimeUnits(gh); + int timeunits2 = (status == 0) ? get_timeunits(unitsOfTime) : timeunits; + // timeunits2 = gribapi_get_timeunits(gh); long lpar; status = grib_get_long(gh, "forecastTime", &lpar); @@ -108,13 +108,13 @@ gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endStep) else { status = grib_get_long(gh, "startStep", &lpar); - if (status == 0) *startStep = (int) (((double) lpar * timeunit_factor(timeunits, timeunits2)) + 0.5); + if (status == 0) *startStep = (int) (((double) lpar * timeunits_factor(timeunits, timeunits2)) + 0.5); } *endStep = *startStep; status = grib_get_long(gh, "endStep", &lpar); - if (status == 0) *endStep = (int) (((double) lpar * timeunit_factor(timeunits, timeunits2)) + 0.5); - // printf("%d %d %d %d %d %g\n", *startStep, *endStep, lpar, timeunits, timeunits2, timeunit_factor(timeunits, timeunits2)); + if (status == 0) *endStep = (int) (((double) lpar * timeunits_factor(timeunits, timeunits2)) + 0.5); + // printf("%d %d %d %d %d %g\n", *startStep, *endStep, lpar, timeunits, timeunits2, timeunits_factor(timeunits, timeunits2)); } static CdiDateTime @@ -152,7 +152,7 @@ gribapiSetDataDateTime(grib_handle *gh, CdiDateTime dataDateTime) } static int -gribapiGetTimeUnitFactor(int timeUnits) +gribapi_get_timeunits_factor(int timeUnits) { static bool lprint = true; switch (timeUnits) @@ -197,7 +197,7 @@ gribapiGetValidityDateTime(grib_handle *gh, CdiDateTime *sDateTime) { CdiDateTime rDateTime = gribapiGetDataDateTime(gh); - int timeUnits = gribapiGetTimeUnits(gh); + int timeUnits = gribapi_get_timeunits(gh); int startStep = 0, endStep = 0; gribapiGetSteps(gh, timeUnits, &startStep, &endStep); @@ -206,7 +206,7 @@ gribapiGetValidityDateTime(grib_handle *gh, CdiDateTime *sDateTime) extern int CGRIBEX_grib_calendar; JulianDate julianDate = julianDate_encode(CGRIBEX_grib_calendar, rDateTime); - int64_t timeUnitFactor = gribapiGetTimeUnitFactor(timeUnits); + int64_t timeUnitFactor = gribapi_get_timeunits_factor(timeUnits); // if (startStep > 0) { @@ -224,7 +224,7 @@ gribapiGetValidityDateTime(grib_handle *gh, CdiDateTime *sDateTime) } static void -grib1GetLevel(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2) +grib1_get_level(grib_handle *gh, int *leveltype, int *lbounds, int *level1, int *level2) { *leveltype = 0; *lbounds = 0; @@ -300,7 +300,7 @@ grib2ScaleFactor(long factor) } static int -calcLevel(int level_sf, long factor, long level) +calc_level(int level_sf, long factor, long level) { double result = 0; if (level != GRIB_MISSING_LONG) result = (double) level * grib2ScaleFactor(factor); @@ -309,8 +309,8 @@ calcLevel(int level_sf, long factor, long level) } static void -grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, int *level2, int *level_sf, - int *level_unit) +grib2_get_level(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, int *level2, int *level_sf, + int *level_unit) { *leveltype1 = 0; *leveltype2 = -1; @@ -353,31 +353,31 @@ grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, i long factor, llevel; GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0); // 1 byte GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0); // 4 byte - *level1 = calcLevel(*level_sf, factor, llevel); + *level1 = calc_level(*level_sf, factor, llevel); if (*lbounds) { GRIB_CHECK(grib_get_long(gh, "scaleFactorOfSecondFixedSurface", &factor), 0); // 1 byte GRIB_CHECK(grib_get_long(gh, "scaledValueOfSecondFixedSurface", &llevel), 0); // 4 byte - *level2 = calcLevel(*level_sf, factor, llevel); + *level2 = calc_level(*level_sf, factor, llevel); } } } static void -gribGetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, int *level2, int *level_sf, - int *level_unit, var_tile_t *tiles) +grib_get_level(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, int *level1, int *level2, int *level_sf, + int *level_unit, var_tile_t *tiles) { if (gribEditionNumber(gh) <= 1) { - grib1GetLevel(gh, leveltype1, lbounds, level1, level2); + grib1_get_level(gh, leveltype1, lbounds, level1, level2); *leveltype2 = -1; *level_sf = 0; *level_unit = 0; } else { - grib2GetLevel(gh, leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit); + grib2_get_level(gh, leveltype1, leveltype2, lbounds, level1, level2, level_sf, level_unit); // read in tiles attributes (if there are any) tiles->tileindex = (int) gribGetLongDefault(gh, cdiSubtypeAttributeName[SUBTYPE_ATT_TILEINDEX], 0); @@ -391,7 +391,7 @@ gribGetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lbounds, in } static void -gribapiGetString(grib_handle *gh, const char *key, char *string, size_t length) +gribapi_get_string(grib_handle *gh, const char *key, char *string, size_t length) { string[0] = 0; @@ -423,21 +423,21 @@ param_to_name(int param, char *name) } static int -gribapiGetEnsembleInfo(grib_handle *gh, long *typeOfEnsembleForecast, long *numberOfForecastsInEnsemble, long *perturbationNumber) +gribapiGetEnsembleInfo(grib_handle *gh, long *numberOfForecastsInEnsemble, long *perturbationNumber, long *typeOfEnsembleForecast) { int status = 0; - if (grib_get_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast) == 0) + if (grib_get_long(gh, "numberOfForecastsInEnsemble", numberOfForecastsInEnsemble) == 0) { - GRIB_CHECK(grib_get_long(gh, "numberOfForecastsInEnsemble", numberOfForecastsInEnsemble), 0); + if (*numberOfForecastsInEnsemble > 0) status = 1; GRIB_CHECK(grib_get_long(gh, "perturbationNumber", perturbationNumber), 0); - if (*perturbationNumber > 0) status = 1; + grib_get_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast); } if (status == 0) { - *typeOfEnsembleForecast = 0; - *perturbationNumber = 0; - *numberOfForecastsInEnsemble = 0; + *numberOfForecastsInEnsemble = -1; + *perturbationNumber = -1; + *typeOfEnsembleForecast = -1; } return status; @@ -449,8 +449,8 @@ gribapiGetScanKeys(grib_handle *gh) VarScanKeys scanKeys; varScanKeysInit(&scanKeys); - long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0; - gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber); + long numberOfForecastsInEnsemble = -1, perturbationNumber = -1, typeOfEnsembleForecast = -1; + gribapiGetEnsembleInfo(gh, &numberOfForecastsInEnsemble, &perturbationNumber, &typeOfEnsembleForecast); scanKeys.perturbationNumber = (short) perturbationNumber; long typeOfGeneratingProcess = 0; @@ -466,10 +466,10 @@ gribapiGetNameKeys(grib_handle *gh, int varID) char string[CDI_MAX_NAME]; size_t vlen = CDI_MAX_NAME; - gribapiGetString(gh, "name", string, vlen); // longname + gribapi_get_string(gh, "name", string, vlen); // longname if (string[0]) varDefKeyString(varID, CDI_KEY_LONGNAME, string); - gribapiGetString(gh, "units", string, vlen); + gribapi_get_string(gh, "units", string, vlen); if (string[0]) varDefKeyString(varID, CDI_KEY_UNITS, string); string[0] = 0; @@ -518,16 +518,16 @@ gribapiGetKeys(grib_handle *gh, int varID) */ /* - Get the ensemble Info from the grib-2 Tables and update the intermediate datastructure. + Get the ensemble information from the grib-2 Tables and update the intermediate datastructure. Further update to the "vlist" is handled in the same way as for GRIB-1 by "cdi_generate_vars" */ - long typeOfEnsembleForecast = 0, numberOfForecastsInEnsemble = 0, perturbationNumber = 0; - gribapiGetEnsembleInfo(gh, &typeOfEnsembleForecast, &numberOfForecastsInEnsemble, &perturbationNumber); - if (perturbationNumber > 0) + long numberOfForecastsInEnsemble = -1, perturbationNumber = -1, typeOfEnsembleForecast = -1; + gribapiGetEnsembleInfo(gh, &numberOfForecastsInEnsemble, &perturbationNumber, &typeOfEnsembleForecast); + if (numberOfForecastsInEnsemble > 0) { - varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast); varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, (int) numberOfForecastsInEnsemble); varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, (int) perturbationNumber); + if (typeOfEnsembleForecast != -1) varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, (int) typeOfEnsembleForecast); } long section2Length = 0; @@ -678,6 +678,24 @@ gribapiDefProjSTERE(grib_handle *gh, int gridID) gridDefParamsSTERE(gridID, gpp); } +static void +gribapiDefProjHEALPIX(grib_handle *gh, int gridID) +{ + struct CDI_GridProjParams gpp; + gridProjParamsInit(&gpp); + + decode_shapeOfTheEarth(gh, &gpp); + + long lval = -1; + grib_get_long(gh, "Nside", &lval); + gpp.nside = (int) lval; + lval = -1; + grib_get_long(gh, "ordering", &lval); + gpp.order = (int) lval; + + gridDefParamsHEALPIX(gridID, gpp); +} + static void gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh, size_t recsize, off_t position, int datatype, int comptype, const char *varname, int leveltype1, int leveltype2, int lbounds, int level1, int level2, int level_sf, @@ -704,6 +722,8 @@ gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh, size_t recsize record->tiles = tiles ? *tiles : dummy_tiles; #ifdef HAVE_LIBFDB5 record->fdbItem = fdbItem; +#else + (void) fdbItem; #endif strncpy(record->varname, varname, sizeof(record->varname) - 1); @@ -734,12 +754,13 @@ gribapiAddRecord(stream_t *streamptr, int param, grib_handle *gh, size_t recsize grid_free(grid); Free(grid); } - else if (grid->projtype == CDI_PROJ_RLL) gribapiDefProjRLL(gh, gridID); - else if (grid->projtype == CDI_PROJ_LCC) gribapiDefProjLCC(gh, gridID); - else if (grid->projtype == CDI_PROJ_STERE) gribapiDefProjSTERE(gh, gridID); + else if (grid->projtype == CDI_PROJ_RLL) gribapiDefProjRLL(gh, gridID); + else if (grid->projtype == CDI_PROJ_LCC) gribapiDefProjLCC(gh, gridID); + else if (grid->projtype == CDI_PROJ_STERE) gribapiDefProjSTERE(gh, gridID); + else if (grid->projtype == CDI_PROJ_HEALPIX) gribapiDefProjHEALPIX(gh, gridID); // clang-format on - int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1); + int zaxistype = gribapi_get_zaxis_type(gribEditionNumber(gh), leveltype1); switch (zaxistype) { @@ -1003,6 +1024,7 @@ checkTime(stream_t *streamptr, const compvar2_t *compVar, CdiDateTime verificati int fdbScanTimesteps(stream_t *streamptr) { + (void) streamptr; #ifdef HAVE_LIBFDB5 void *gribbuffer = NULL; size_t buffersize = 0; @@ -1066,10 +1088,10 @@ fdbScanTimesteps(stream_t *streamptr) int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit; var_tile_t tiles = dummy_tiles; int level1 = 0, level2 = 0; - gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles); + grib_get_level(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles); char varname[256]; - gribapiGetString(gh, "shortName", varname, sizeof(varname)); + gribapi_get_string(gh, "shortName", varname, sizeof(varname)); int param = gribapiGetParam(gh); @@ -1083,7 +1105,7 @@ fdbScanTimesteps(stream_t *streamptr) vDateTime0 = vDateTime; taxis->rDateTime = gribapiGetDataDateTime(gh); fcast = gribapiTimeIsFC(gh); - if (fcast) taxis->unit = gribapiGetTimeUnits(gh); + if (fcast) taxis->unit = gribapi_get_timeunits(gh); taxis->fDateTime = taxis->rDateTime; taxis->sDateTime = sDateTime; taxis->vDateTime = vDateTime; @@ -1227,10 +1249,10 @@ gribapiScanTimestep1(stream_t *streamptr) int leveltype1 = -1, leveltype2 = -1, lbounds, level_sf, level_unit; var_tile_t tiles = dummy_tiles; int level1 = 0, level2 = 0; - gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles); + grib_get_level(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles); char varname[256]; - gribapiGetString(gh, "shortName", varname, sizeof(varname)); + gribapi_get_string(gh, "shortName", varname, sizeof(varname)); int param = gribapiGetParam(gh); @@ -1246,7 +1268,7 @@ gribapiScanTimestep1(stream_t *streamptr) vDateTime0 = vDateTime; taxis->rDateTime = gribapiGetDataDateTime(gh); fcast = gribapiTimeIsFC(gh); - if (fcast) taxis->unit = gribapiGetTimeUnits(gh); + if (fcast) taxis->unit = gribapi_get_timeunits(gh); taxis->fDateTime = taxis->rDateTime; taxis->sDateTime = sDateTime; taxis->vDateTime = vDateTime; @@ -1383,10 +1405,10 @@ gribapiScanTimestep2(stream_t *streamptr) int level1 = 0, level2 = 0, leveltype1, leveltype2, lbounds, level_sf, level_unit; var_tile_t tiles = dummy_tiles; - gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles); + grib_get_level(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles); char varname[256]; - gribapiGetString(gh, "shortName", varname, sizeof(varname)); + gribapi_get_string(gh, "shortName", varname, sizeof(varname)); int param = gribapiGetParam(gh); @@ -1402,7 +1424,7 @@ gribapiScanTimestep2(stream_t *streamptr) if (taxisInqType(taxisID) == TAXIS_RELATIVE) { taxis->type = TAXIS_RELATIVE; - taxis->unit = gribapiGetTimeUnits(gh); + taxis->unit = gribapi_get_timeunits(gh); taxis->rDateTime = gribapiGetDataDateTime(gh); } else @@ -1562,14 +1584,14 @@ gribapiScanTimestep(stream_t *streamptr) int level1 = 0, level2 = 0, leveltype1, leveltype2 = -1, lbounds, level_sf, level_unit; var_tile_t tiles = dummy_tiles; - gribGetLevel(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles); + grib_get_level(gh, &leveltype1, &leveltype2, &lbounds, &level1, &level2, &level_sf, &level_unit, &tiles); CdiDateTime sDateTime; CdiDateTime vDateTime = gribapiGetValidityDateTime(gh, &sDateTime); if (rindex == nrecs) break; - gribapiGetString(gh, "shortName", varname, sizeof(varname)); + gribapi_get_string(gh, "shortName", varname, sizeof(varname)); int param = gribapiGetParam(gh); @@ -1582,7 +1604,7 @@ gribapiScanTimestep(stream_t *streamptr) if (taxisInqType(taxisID) == TAXIS_RELATIVE) { taxis->type = TAXIS_RELATIVE; - taxis->unit = gribapiGetTimeUnits(gh); + taxis->unit = gribapi_get_timeunits(gh); taxis->rDateTime = gribapiGetDataDateTime(gh); } else @@ -1952,6 +1974,8 @@ grib2ProDefTempHasStatisticalDef(int proDefTempNum) case 46: case 47: case 61: + case 67: + case 68: case 91: case 1001: case 1101: @@ -2377,6 +2401,33 @@ gribapiDefGridSTERE(grib_handle *gh, int gridID, int uvRelativeToGrid) GRIB_CHECK(my_grib_set_long(gh, "jScansPositively", (yinc > 0)), 0); } +static void +gribapiDefGridHEALPIX(grib_handle *gh, int gridID, int uvRelativeToGrid) +{ + struct CDI_GridProjParams gpp; + gridInqParamsHEALPIX(gridID, &gpp); + gridVerifyProjParamsHEALPIX(&gpp); + // if (gpp.xval_0 < 0.0) gpp.xval_0 += 360.0; + + static const char mesg[] = "healpix"; + size_t len = sizeof(mesg) - 1; + GRIB_CHECK(my_grib_set_string(gh, "gridType", mesg, &len), 0); + + GRIB_CHECK(my_grib_set_long(gh, "Nside", gpp.nside), 0); + GRIB_CHECK(my_grib_set_long(gh, "ordering", gpp.order), 0); + GRIB_CHECK(my_grib_set_double(gh, "longitudeOfFirstGridPointInDegrees", 45.0), 0); + // GRIB_CHECK(my_grib_set_double(gh, "longitudeOfFirstGridPointInDegrees", gpp.xval_0), 0); + /* + long shapeOfTheEarth = encode_shapeOfTheEarth(&gpp); + if (shapeOfTheEarth) GRIB_CHECK(my_grib_set_long(gh, "shapeOfTheEarth", shapeOfTheEarth), 0); + if (shapeOfTheEarth == 1) GRIB_CHECK(my_grib_set_long(gh, "radiusOfTheEarth", gpp.a), 0); + + long earthIsOblate = (shapeOfTheEarth == 2 || shapeOfTheEarth == 3 || shapeOfTheEarth == 4); + if (earthIsOblate) GRIB_CHECK(my_grib_set_long(gh, "earthIsOblate", earthIsOblate), 0); + */ + if (uvRelativeToGrid >= 0) GRIB_CHECK(my_grib_set_long(gh, "uvRelativeToGrid", uvRelativeToGrid), 0); +} + static void gribapiDefGridGME(grib_handle *gh, int gridID, long gridsize) { @@ -2540,6 +2591,12 @@ gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype, int gribapiDefGridSTERE(gh, gridID, uvRelativeToGrid); break; } + case CDI_PROJ_HEALPIX: + { + if (editionNumber <= 1) Error("HEALPix grid can't be stored in GRIB edition %d!", editionNumber); + gribapiDefGridHEALPIX(gh, gridID, uvRelativeToGrid); + break; + } case GRID_SPECTRAL: { gribapiDefGridSpectral(gh, gridID); @@ -3170,6 +3227,7 @@ static void gribapiSetExtMode(grib_handle *gh, int gridID, size_t datasize, const void *data) { #ifndef HIRLAM_EXTENSIONS + (void) gh; (void) data; (void) datasize; #endif @@ -3288,12 +3346,12 @@ gribapiEncode(int memType, int varID, int levelID, int vlistID, int gridID, int { int status, perturbationNumber, numberOfForecastsInEnsemble, typeOfEnsembleForecast; - status = cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, &typeOfEnsembleForecast); - if (status == 0) grib_set_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast); status = cdiInqKeyInt(vlistID, varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, &numberOfForecastsInEnsemble); if (status == 0) grib_set_long(gh, "numberOfForecastsInEnsemble", numberOfForecastsInEnsemble); status = cdiInqKeyInt(vlistID, varID, CDI_KEY_PERTURBATIONNUMBER, &perturbationNumber); if (status == 0) grib_set_long(gh, "perturbationNumber", perturbationNumber); + status = cdiInqKeyInt(vlistID, varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, &typeOfEnsembleForecast); + if (status == 0) grib_set_long(gh, "typeOfEnsembleForecast", typeOfEnsembleForecast); } if (!gc->init) gribapiDefInstitut(gh, vlistID, varID); @@ -3411,7 +3469,9 @@ gribapiEncode(int memType, int varID, int levelID, int vlistID, int gridID, int #endif } else - GRIB_CHECK(grib_set_double_array(gh, "values", (const double *) data, datasize), 0); + { + GRIB_CHECK(grib_set_double_array(gh, "values", (const double *) data, datasize), 0); + } if (nmiss) { @@ -3438,7 +3498,7 @@ gribapiEncode(int memType, int varID, int levelID, int vlistID, int gridID, int cdi_name[0] = 0; vlistInqVarName(vlistID, varID, cdi_name); char grb_name[256]; - gribapiGetString(gh, "shortName", grb_name, sizeof(grb_name)); + gribapi_get_string(gh, "shortName", grb_name, sizeof(grb_name)); str_to_lower(cdi_name); str_to_lower(grb_name); bool checkName = (!grb_name[0] && strncmp(cdi_name, "param", 5) == 0) ? false : true;