diff --git a/CMakeLists.txt b/CMakeLists.txt index 7fb56d3c..3f8c42c5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,7 +16,7 @@ option(USE_HDF5 "Use HDF5?" off) option(USE_REGEX "Use Regex?" on) option(USE_TIGGE "Use tigge?" on) option(USE_MYSQL "Use MySQL?" off) -option(USE_IPOLATES "iplib=0,1,3?" 3) +option(USE_IPOLATES "Use Ipolates" on) option(USE_UDF "Use UDF?" off) option(USE_OPENMP "Use OpenMP?" on) option(USE_PROJ4 "Use Proj4?" off) @@ -83,17 +83,15 @@ endif() # If user wants to use NCEPLIBS-ip, find it and the sp library. message(STATUS "Checking if the user want to use NCEPLIBS-ip...") -if(USE_IPOLATES EQUAL 1) +if(USE_IPOLATES) find_package(ip CONFIG REQUIRED) - find_package(sp CONFIG REQUIRED) -elseif(USE_IPOLATES EQUAL 3) - find_package(ip2 CONFIG REQUIRED) + list(APPEND definitions_list -DIPOLATES_LIB="ipolates_lib_d") endif() -if(NOT USE_IPOLATES EQUAL 0) - list(APPEND definitions_list -DIPOLATES_LIB="ipolates_lib") -endif() -list(APPEND definitions_list -DUSE_IPOLATES=${USE_IPOLATES}) +# if(NOT USE_IPOLATES EQUAL 0) +# list(APPEND definitions_list -DIPOLATES_LIB="ipolates_lib") +# endif() +# list(APPEND definitions_list -DUSE_IPOLATES=${USE_IPOLATES}) message(STATUS "Checking if the user want to use NetCDF...") if(USE_NETCDF4) diff --git a/wgrib2/CMakeLists.txt b/wgrib2/CMakeLists.txt index e0508874..34734269 100644 --- a/wgrib2/CMakeLists.txt +++ b/wgrib2/CMakeLists.txt @@ -24,7 +24,7 @@ Import_grib_fs.c Import_lonlat.c Import_netcdf.c init.c int8.c intpower.c Inv.c Inv_no.c Irr_grids.c itoshort_a.c JMA.c jpeg_pk.c Last.c lat2ij.c Latlon.c Level.c Limit.c Lola.c Lvl.c Macro.c manage_inv_out.c Match.c Match_fs.c Match_inv.c Mem_buffer.c Merge.c -Misc.c missing.c mk_gdt.c mk_kgds.c Model_version_date.c Mod_grib.c +Misc.c missing.c mk_gdt.c Model_version_date.c Mod_grib.c Mysql.c Mysql_dump.c Mysql_speed.c Names.c ncep_grids.c NCEP_norm.c NCEP_uv.c Ncpu.c Ndate.c Ndates.c Netcdf.c Netcdf_sup.c New_grid.c new_grid_lambertc.c New_grid_order.c openmp_util.c parse_loop.c @@ -105,9 +105,8 @@ endif() if(USE_G2CLIB) endif() -if(USE_IPOLATES EQUAL 1) +if(USE_IPOLATES) target_link_libraries(obj_lib PUBLIC ip::ip_d) - target_link_libraries(obj_lib PUBLIC sp::sp_d) # Link to the Fortran runtime library for each compiler if using ip2. # The wgrib2 exectuable is created using the C compiler and diff --git a/wgrib2/Config.c b/wgrib2/Config.c index 7e55577b..d484fb42 100644 --- a/wgrib2/Config.c +++ b/wgrib2/Config.c @@ -89,13 +89,9 @@ int f_config(ARG0) { strcat(inv_out, "tigge package is not installed\n"); #endif -#if USE_IPOLATES > 0 +#ifdef USE_IPOLATES inv_out += strlen(inv_out); - sprintf(inv_out, "IPOLATES " IPOLATES_LIB " (option %d) is installed", USE_IPOLATES); -#if USE_SPECTRAL > 0 - strcat(inv_out, " with spectral interpolation"); -#endif - strcat(inv_out, ", default vectors:\n"); + sprintf(inv_out, "IPOLATES NCEPLIBS-ip\n"); #else strcat(inv_out, "interpolation package is not installed, default vectors:\n"); #endif diff --git a/wgrib2/New_grid.c b/wgrib2/New_grid.c index 8247dea3..1fcbbf5e 100644 --- a/wgrib2/New_grid.c +++ b/wgrib2/New_grid.c @@ -32,6 +32,8 @@ * this makes dlon only xxx.xx precision whic is less * than grib1 (xxx.xxx) and grib2 (xxx.xxxxxx) expectations * should not use float version of ip2. + * 4/2024 changed to current ip lib + * removed calls to old ip lib (1) */ const char **vectors; @@ -47,6 +49,8 @@ enum new_grid_format_type new_grid_format; static const char *no_vectors[] = { NULL }; static const char *UV_vectors[] = { "UGRD", "VGRD", NULL }; +#ifdef USE_IPOLATES + /* * HEADER:111:new_grid_vectors:misc:1:change fields to vector interpolate: X=none,default,UGRD:VGRD,(U:V list) */ @@ -122,617 +126,15 @@ int f_new_grid_vectors(ARG1) { return 0; } -#if USE_IPOLATES == 1 - -/* This is the grib1 version which is now frozen. - * The grib2 version is similar but - * (1) uses the grib2 version of ipolates - * (2) double precision - * (3) spectral support was removed - */ - -#ifndef _IPOLATES_H -#define _IPOLATES_H -#include "ipolates.h" -#endif - -extern unsigned int npnts; -extern int nx,ny; -extern double *lat, *lon; -extern int decode, latlon, file_append, flush_mode; -extern int use_scale, dec_scale, bin_scale, wanted_bits, max_bits; -extern enum output_grib_type grib_type; -extern enum output_order_type output_order; -extern int save_translation; -extern enum output_order_type output_order_wanted, output_order; - -static int interpol_type = 0; -static int ipopt[20] = {-1,-1,0, 0,0,0, 0,0,0, 0}; - -/* - * HEADER:111:new_grid_interpolation:misc:1:new_grid interpolation X=bilinear,bicubic,neighbor,budget,neighbor-budget - */ - -int f_new_grid_interpolation(ARG1) { - - if (strcmp(arg1,"bilinear") == 0) { interpol_type = 0; ipopt[0] = -1; } - else if (strcmp(arg1,"bicubic") == 0) { interpol_type = 1; ipopt[0] = 0; } - else if (strcmp(arg1,"neighbor") == 0) { interpol_type = 2; ipopt[0] = 1; } - else if (strcmp(arg1,"budget") == 0) { interpol_type = 3; ipopt[0] = -1; } -// turned off spectral -- new library rarely used interpolation option -// else if (strcmp(arg1,"spectral") == 0) { interpol_type = 4; ipopt[0] = 0; ipopt[1] = -1; } -// turned off neighbor-budget - save space for rarely used interpolation option -// else if (strcmp(arg1,"neighbor-budget") == 0) { interpol_type = 6; ipopt[0] = -1; } - else fatal_error("new_grid_interpolation: unknown type %s", arg1); - - return 0; -} - -/* - * HEADER:111:new_grid_ipopt:misc:1:new_grid ipopt values X=i1:i2..:iN N <= 20 - */ -int f_new_grid_ipopt(ARG1) { - int i, k, val, m; - - i = 0; - k = sscanf(arg1, "%d%n", &val, &m); - while (k == 1) { - if (i > 19) fatal_error("new_grid_ipopt: too many ipopt values, 20 max",""); - ipopt[i++] = val; - arg1 += m; - k = sscanf(arg1, ":%d%n", &val, &m); - } - return 0; -} - -/* - * HEADER:111:new_grid_format:misc:1:new_grid output format X=bin,ieee,grib - */ - - -int f_new_grid_format(ARG1) { - if (mode == -1) { - if (strcmp(arg1,"grib") != 0) { - fatal_error("new_grid_format: for IPOLATES=iplib, only grib is supported",""); - } - } - return 0; -} - -/* - * HEADER:111:new_grid_winds:misc:1:new_grid wind orientation: X = grid, earth (no default) - */ - -enum wind_rotation_type wind_rotation; - -int f_new_grid_winds(ARG1) { - int *save; - if (mode == -2) { - free(*local); - return 0; - } - if (mode == -1) { - if ((*local = save = (int *) malloc(sizeof(int))) == NULL) fatal_error("new_grid_winds: malloc",""); - if (strcmp(arg1,"grid") == 0) *save = 0; - else if (strcmp(arg1,"earth") == 0) *save = 1; - else fatal_error("new_grid_winds: bad arg %s", arg1); - } - save = (int *) *local; - wind_rotation = (*save) ? earth : grid; - return 0; -} - -struct local_struct { - // U data - float *u_val; - int has_u, nx, ny; - unsigned char *clone_sec[9]; - char name[NAMELEN]; - - // interpolation - int npnts_out; // must be integer .. fortran call requires int - float *rlat, *rlon, *crot, *srot; - unsigned char *sec3; - int kgds_out[200]; - double radius_major, radius_minor; - - // output file - struct seq_file out; -}; - -unsigned char blank_sec1[21] = { 0,0,0,21,1, - 255,255,255,255, // center subcenter - 2,1,255, // grib master table, local table, sig ref time - 255, 255, // year - 255, 255, 255, 255, 255, // month .. second - 255, 255}; - -/* - * HEADER:111:new_grid:output:4:bilinear interpolate: X=projection Y=x0:nx:dx Z=y0:ny:dy A=grib_file alpha - */ - -int f_new_grid(ARG4) { - struct local_struct *save; - - unsigned int i; - int is_u, is_v, ftn_npnts, ftn_nout; - int kgds[200], km; - float *data_in, *data_out; - double x0, y0, dx, dy, xn, yn; - double lov, lad, latin1, latin2; - int proj; // projection: for LC 0 = NP, 128 = SP - char name[NAMELEN]; - int j, ibi, ibo, iret, nnx, nny, n_out; - unsigned char *new_sec[8], *s, *bitmap, *bitmap_out, *p; - - /* for lambertc */ - double r_maj, r_min, ref_lon, ref_lat; - - if (mode == -1) { // initialization - decode = 1; - output_order_wanted = raw; // in raw order - - -#ifdef G95 - // initialize g95 runtime library - if (g95_runstop == 0) { g95_runtime_start(0,NULL); g95_runstop = 1; } -#endif - -// if ( (sizeof(vectors) / sizeof (vectors[0])) % 2 == 1) fatal_error("new_grid: program error in vectors[]",""); - - // allocate static variables - - *local = save = (struct local_struct *) malloc( sizeof(struct local_struct)); - if (save == NULL) fatal_error("new_grid: memory allocation error",""); - - if (fopen_file(&(save->out), arg4, file_append ? "ab" : "wb") != 0) { - fatal_error("-new_grid: could not open file %s", arg4); - } - - save->has_u = 0; - save->radius_major = save->radius_minor = 0.0; - init_sec(save->clone_sec); - s = NULL; - - // parse NCEP grids */ - ncep_grids(&arg1, &arg2, &arg3); - - // for each output grid - if (strcmp(arg1,"latlon") == 0) { - if (sscanf(arg2,"%lf:%d:%lf", &x0, &nnx, &dx) != 3) - fatal_error("new_grid: XDEF wrong:%s",arg2); - if (sscanf(arg3,"%lf:%d:%lf", &y0, &nny, &dy) != 3) - fatal_error("new_grid: YDEF wrong:%s",arg3); - - if (x0 < 0.0) x0 += 360.0; - save->nx = nnx; - save->ny = nny; - save->npnts_out = n_out = nnx*nny; - if (n_out <= 0) fatal_error("new_grid: bad nx, ny",""); - - // make a new section 3 - s = sec3_lola(nnx, x0, dx, nny, y0, dy, sec); - } - else if (strcmp(arg1,"rot-latlon") == 0) { - if (sscanf(arg2,"%lf:%d:%lf", &x0, &nnx, &dx) != 3) - fatal_error("new_grid: XDEF wrong:%s",arg2); - if (sscanf(arg3,"%lf:%d:%lf", &y0, &nny, &dy) != 3) - fatal_error("new_grid: YDEF wrong:%s",arg3); - - if (x0 < 0.0) x0 += 360.0; - save->nx = nnx; - save->ny = nny; - save->npnts_out = n_out = nnx*nny; - if (n_out <= 0) fatal_error("new_grid: bad nx, ny",""); - - // make a new section 3 - s = sec3_lola(nnx, x0, dx, nny, y0, dy, sec); - } - - - else if (strncmp(arg1,"mercator:",9) == 0) { - if (sscanf(arg1,"mercator:%lf", &lad) != 1) - fatal_error("new_grid: LaD (latitude interesection) not specified",""); - if (sscanf(arg2,"%lf:%d:%lf:%lf", &x0, &nnx, &dx, &xn) != 4) - fatal_error("new_grid: XDEF wrong:%s",arg2); - if (sscanf(arg3,"%lf:%d:%lf:%lf", &y0, &nny, &dy, &yn) != 4) - - if (x0 < 0.0) x0 += 360.0; - save->nx = nnx; - save->ny = nny; - save->npnts_out = n_out = nnx*nny; - if (n_out <= 0) fatal_error("new_grid: bad nx, ny",""); - - // make a new section 3 - s = sec3_mercator(lad, nnx, x0, dx, xn, nny, y0, dy, yn, sec); - } - else if (strcmp(arg1,"gaussian") == 0) { - if (sscanf(arg2,"%lf:%d:%lf", &x0, &nnx, &dx) != 3) - fatal_error("new_grid: XDEF wrong:%s",arg2); - if (sscanf(arg3,"%lf:%d", &y0, &nny) != 2) - fatal_error("new_grid: YDEF wrong:%s",arg3); - - if (x0 < 0.0) x0 += 360.0; - save->nx = nnx; - save->ny = nny; - save->npnts_out = n_out = nnx*nny; - if (n_out <= 0) fatal_error("new_grid: bad nx, ny",""); - // make a new section 3 - s = sec3_gaussian(nnx, x0, dx, nny, y0, sec); - } - else if (strncmp(arg1,"lambert:",8) == 0) { - i = sscanf(arg1,"lambert:%lf:%lf:%lf:%lf", &lov, &latin1, &latin2, &lad); - if (i < 2) fatal_error("new_grid: arg1 wrong:%s",arg1); - if (lov < 0.0) lov += 360.0; - if (i < 3) latin2 = latin1; - if (i < 4) lad = latin2; - proj = 0; - if (latin2 < 0.0) proj = 128; - - if (sscanf(arg2,"%lf:%d:%lf", &x0, &nnx, &dx) != 3) - fatal_error("new_grid: XDEF wrong:%s",arg2); - if (sscanf(arg3,"%lf:%d:%lf", &y0, &nny, &dy) != 3) - fatal_error("new_grid: YDEF wrong:%s",arg3); - - if (x0 < 0.0) x0 += 360.0; - save->nx = nnx; - save->ny = nny; - save->npnts_out = n_out = nnx*nny; - if (n_out <= 0) fatal_error("new_grid: bad nx, ny",""); - - // make a new section 3 - s = sec3_lc(lov, lad, latin1, latin2, proj, nnx, x0, dx, nny, y0, dy, sec); - } - - /* for lambertc, input is the lon-lat of center point */ - /* can not calc grid until radius is given, so do lambert code to check args */ - - else if (strncmp(arg1,"lambertc:",9) == 0) { - i = sscanf(arg1,"lambertc:%lf:%lf:%lf:%lf", &lov, &latin1, &latin2, &lad); - if (i < 2) fatal_error("new_grid: arg1 wrong:%s",arg1); - if (lov < 0.0) lov += 360.0; - if (i < 3) latin2 = latin1; - if (i < 4) lad = latin2; - proj = 0; - if (latin2 < 0.0) proj = 128; - - if (sscanf(arg2,"%lf:%d:%lf", &x0, &nnx, &dx) != 3) - fatal_error("new_grid: XDEF wrong:%s",arg2); - if (sscanf(arg3,"%lf:%d:%lf", &y0, &nny, &dy) != 3) - fatal_error("new_grid: YDEF wrong:%s",arg3); - - if (x0 < 0.0) x0 += 360.0; - save->nx = nnx; - save->ny = nny; - save->npnts_out = n_out = nnx*nny; - if (n_out <= 0) fatal_error("new_grid: bad nx, ny",""); - - // make a new section 3 - s = sec3_lc(lov, lad, latin1, latin2, proj, nnx, x0, dx, nny, y0, dy, sec); - } - - else if (strncmp(arg1,"nps:",4) == 0 || strncmp(arg1,"sps:",4) == 0) { - if (sscanf(arg1,"%*[ns]ps:%lf:%lf", &lov, &lad) != 2) fatal_error("new_grid: arg1 wrong:%s",arg1); - if (lad != 60.0) fatal_error("New_grid: only LatD = 60 is supported",""); - proj = 0; - if (arg1[0] == 's') proj = 128; - if (sscanf(arg2,"%lf:%d:%lf", &x0, &nnx, &dx) != 3) - fatal_error("new_grid: XDEF wrong:%s",arg2); - if (sscanf(arg3,"%lf:%d:%lf", &y0, &nny, &dy) != 3) - fatal_error("new_grid: YDEF wrong:%s",arg3); - if (lov < 0.0) lov += 360.0; - - if (x0 < 0.0) x0 += 360.0; - save->nx = nnx; - save->ny = nny; - save->npnts_out = n_out = nnx*nny; - if (n_out <= 0) fatal_error("new_grid: bad nx, ny",""); - - // make a new section 3 - s = sec3_polar_stereo(lov, lad, proj, nnx, x0, dx, nny, y0, dy, sec); - } - else fatal_error("new_grid: unsupported output grid %s", arg1); - - new_sec[1] = blank_sec1; // add center info - // save new section 3 - i = (int) uint4(s); // size of section 3 - new_sec[3] = save->sec3 = (unsigned char *) malloc(i * sizeof(unsigned char)); - for (j = 0; j < i; j++) save->sec3[j] = s[j]; - - // apply wind rotation .. change flag 3.3 - if (wind_rotation == undefined && strncmp(arg1,"grib",4) != 0) { - fprintf(stderr,"Warning: -new_grid wind orientation undefined, " - "use \"-new_grid_winds (grid|earth)\", earth used\n"); - wind_rotation = earth; - } - if ((p = flag_table_3_3_location(new_sec)) != NULL) { - if (wind_rotation == grid) *p = *p | 8; - else if (wind_rotation == earth) *p = *p & (255 - 8); - } - - if (mk_kgds(new_sec, save->kgds_out)) fatal_error("new_grid: encoding output kgds",""); - - /* some vectors need by interpolation routines */ - if ((save->rlat = (float *) malloc( sizeof(float) * (size_t) n_out)) == NULL) - fatal_error("new_grid memory allocation",""); - if ((save->rlon = (float *) malloc(sizeof(float) * (size_t) n_out)) == NULL) - fatal_error("new_grid memory allocation",""); - if ((save->crot = (float *) malloc(sizeof(float) * (size_t) n_out)) == NULL) - fatal_error("new_grid memory allocation",""); - if ((save->srot = (float *) malloc(sizeof(float) * (size_t) n_out)) == NULL) - fatal_error("new_grid memory allocation",""); - - return 0; - } - - save = (struct local_struct *) *local; - - if (mode == -2) { // cleanup -#ifdef G95 - if (g95_runstop == 1) { g95_runtime_stop(); g95_runstop = 0; } -#endif - if (save->has_u > 0) { - fprintf(stderr,"-new_grid: last field %s was not interpolated (missing V)\n", save->name); - free(save->u_val); - free_sec(save->clone_sec); - } - free(save->rlon); - free(save->rlat); - free(save->crot); - free(save->srot); - free(save->sec3); - fclose_file(&(save->out)); - free(save); - - return 0; - } - - if (mode >= 0) { // processing - - /* The kgds of some output grids will change depending on input grid */ - /* for example, radius of earth is not known grib file is read, */ - /* and mass vs wind fields */ - /* right nowm, only affects lambertc */ - - if (strncmp(arg1,"lambertc:",8) == 0) { - - // lambertc depends on the radius of the earth which is - // set by the input grib file - - /* read earth radius */ - i = axes_earth(sec, &r_maj, &r_min, NULL); - if (i) fatal_error_i("axes_earth: error code %d", i); - - if (save->radius_major != r_maj || save->radius_minor != r_min) { - - // update sec3 and kgds - - i = sscanf(arg1,"lambertc:%lf:%lf:%lf:%lf", &lov, &latin1, &latin2, &lad); - if (i < 2) fatal_error("new_grid: arg1 wrong:%s",arg1); - if (lov < 0.0) lov += 360.0; - if (i < 3) latin2 = latin1; - if (i < 4) lad = latin2; - proj = 0; - if (latin2 < 0.0) proj = 128; - - if (sscanf(arg2,"%lf:%d:%lf", &x0, &nnx, &dx) != 3) - fatal_error("new_grid: XDEF wrong:%s",arg2); - if (sscanf(arg3,"%lf:%d:%lf", &y0, &nny, &dy) != 3) - fatal_error("new_grid: YDEF wrong:%s",arg3); - - if (x0 < 0.0) x0 += 360.0; - save->nx = nnx; - save->ny = nny; - save->npnts_out = n_out = nnx*nny; - if (n_out <= 0) fatal_error("new_grid: bad nx, ny",""); - - ref_lon = x0; - ref_lat = y0; - - i = new_grid_lambertc(nnx, nny, ref_lon, ref_lat, latin1, latin2, lov, lad, r_maj, r_min, dx, dy, &x0, &y0); - if (i) fatal_error_i("new_grid_lambertc: error code %d", i); - - // make a new section 3 - s = sec3_lc(lov, lad, latin1, latin2, proj, nnx, x0, dx, nny, y0, dy, sec); - - // save new section 3 - i = (int) uint4(s); // size of section 3 - for (j = 0; j < i; j++) save->sec3[j] = s[j]; - - // make kgds - new_sec[1] = blank_sec1; // add center info - new_sec[3] = save->sec3; - if (mk_kgds(new_sec, save->kgds_out)) fatal_error("new_grid: encoding output kgds",""); - - // save radius of earth, to show sec3 and kgds has been done - save->radius_major = r_maj; - save->radius_minor = r_min; - } - } - - if (output_order != raw) fatal_error("new_grid: must be in raw output order",""); - i = getName(sec, mode, NULL, name, NULL, NULL); - is_u = is_v = 0; -// for (j = 0 ; j < sizeof(vectors) / sizeof(vectors[0]); j++) { - for (j = 0; vectors[j] != NULL; j++) { - if (strcmp(name,vectors[j]) == 0) { - if (j % 2 == 0) is_u = 1; - else is_v = 1; - break; - } - } - -// fprintf(stderr, " %s isu %d isv %d has_u %d\n", name, is_u, is_v, save->has_u); -// for (i = 0; i < 12; i++) { printf("kgds_out[%d] = %d ",i,save->kgds_out[i]); } - - // check if V matches expectation - - if (is_v && (save->has_u == 0 || (same_sec0(sec,save->clone_sec) != 1 || - same_sec1(sec,save->clone_sec) != 1 || - same_sec3(sec,save->clone_sec) != 1 || - same_sec4(sec,save->clone_sec) != 1) )) { - fprintf(stderr,"-new_grid: %s doesn't pair with previous vector field, field ignored\n", name); - return 0; - } - - // if U field - save - - if (is_u) { - if (save->has_u > 0) { - fprintf(stderr,"-new_grid: missing V, %s not interpolated\n",save->name); - free(save->u_val); - free_sec(save->clone_sec); - } - copy_sec(sec, save->clone_sec); - copy_data(data,ndata,&(save->u_val)); - GB2_ParmNum(save->clone_sec) = GB2_ParmNum(sec) + 1; - save->has_u = 1; - strncpy(save->name, name,NAMELEN-1); - save->name[NAMELEN-2]=0; - return 0; - } - - // at this point will call polates with either a scalar or vector - - n_out = save->npnts_out; - nnx = save->nx; - nny = save->ny; - km = 1; // only one field - - if (mk_kgds(sec, kgds)) fatal_error("new_grid: encoding input kgds",""); - - data_in = (float *) malloc((1 + (is_v != 0)) * sizeof(float) * (size_t) npnts); - bitmap = (unsigned char *) malloc(sizeof(unsigned char) * (size_t) npnts); - bitmap_out = (unsigned char *) malloc(sizeof(unsigned char) * (size_t) n_out); - data_out = (float *) malloc((1 + (is_v != 0)) * sizeof(float) * (size_t) n_out); - - if (data_in == NULL || data_out == NULL || bitmap == NULL || bitmap_out == NULL) - fatal_error("new_grid: memory allocation problem",""); - - if (is_v) { -#pragma omp parallel for private(i) - for (i = 0; i < npnts; i++) { - if (DEFINED_VAL(data[i]) && DEFINED_VAL(save->u_val[i])) { - data_in[i] = save->u_val[i]; - data_in[i+npnts] = data[i]; - bitmap[i] = 1; - } - else { - data_in[i] = data_in[i + npnts] = 0.0; - bitmap[i] = 0; - } - } - if (mode == 98) fprintf(stderr," UV interpolation %s , %s\n", save->name, name); - } - else { -#pragma omp parallel for private(i) - for (i = 0; i < npnts; i++) { - if (DEFINED_VAL(data[i])) { - data_in[i] = data[i]; - bitmap[i] = 1; - } - else { - data_in[i] = 0.0; - bitmap[i] = 0; - } - } - } - // check if bitmap is used - ibi = 0; // input bitmap is not used - for (i = 0; i < npnts; i++) { - if (bitmap[i] == 0) { - ibi = 1; - break; - } - } - // interpolate - -// for (i = 0; i < 12; i++) { printf("\nkgds_in[%d] = %d out=%d ",i,kgds[i],save->kgds_out[i]); } - ftn_npnts = (int) npnts; - ftn_nout = (int) n_out; - if (is_v) { - IPOLATEV(&interpol_type, ipopt,kgds,save->kgds_out, - &ftn_npnts, &n_out, &km, &ibi, bitmap, data_in, data_in+npnts, - &ftn_nout,save->rlat,save->rlon, save->crot, save->srot, - &ibo, bitmap_out, data_out, data_out + n_out, &iret); - } - else { - IPOLATES(&interpol_type, ipopt,kgds,save->kgds_out, - &ftn_npnts, &n_out, &km, &ibi, bitmap, data_in, &ftn_nout, - save->rlat,save->rlon, &ibo, bitmap_out, data_out, &iret); - } - if (iret != 0) { - for (i = 0; i < 12; i++) { - fprintf(stderr," IPOLATES error: kgds[%d] input %d output %d\n", i+1,kgds[i],save->kgds_out[i]); - } - if (iret == 2) fatal_error("IPOLATES failed, unrecognized input grid or no grid point of output grid is in input grid",""); - if (iret == 3) fatal_error("IPOLATES failed, unrecognized output grid",""); - fatal_error_i("IPOLATES failed, error %d",iret); - - } - n_out = (unsigned int) ftn_nout; - - /* use bitmap to set UNDEFINED values */ - if (ibo == 1) { // has a bitmap - if (is_v) { - for (i = 0; i < n_out; i++) { - if (bitmap_out[i] == 0) data_out[i] = data_out[i+n_out] = UNDEFINED; - } - } - else { - for (i = 0; i < n_out; i++) { - if (bitmap_out[i] == 0) data_out[i] = UNDEFINED; - } - } - } - - // now to write out the grib file - - for (i = 0; i < 8; i++) new_sec[i] = sec[i]; - new_sec[3] = save->sec3; - - if (is_v != 0) { - GB2_ParmNum(new_sec) = GB2_ParmNum(new_sec) - 1; - grib_wrt(new_sec, data_out, n_out, nnx, nny, use_scale, dec_scale, bin_scale, - wanted_bits, max_bits, grib_type, &(save->out)); - GB2_ParmNum(new_sec) = GB2_ParmNum(new_sec) + 1; - grib_wrt(new_sec, data_out+n_out, n_out, nnx, nny, use_scale, dec_scale, bin_scale, - wanted_bits, max_bits, grib_type, &(save->out)); - } - else { - grib_wrt(new_sec, data_out, n_out, nnx, nny, use_scale, dec_scale, bin_scale, - wanted_bits, max_bits, grib_type, &(save->out)); - } - if (flush_mode) fflush_file(&(save->out)); - free(data_in); - free(bitmap); - free(bitmap_out); - free(data_out); - if (is_v != 0) { - save->has_u = 0; - free(save->u_val); - free_sec(save->clone_sec); - } - } - return 0; -} -#endif - -#if USE_IPOLATES == 3 - /* * this is the double precision float and 4-byte integer version of grib2 version of ipolates * this version will not handle 2G+ grids (4 byte integers) * some optimizations relative to grib1 version and the single precision grib2 versions * of the wrappers * - * spectral restored (removed from grib1 version because I didn't want to support it) + * spectral restored */ -#ifndef _IPOLATES_H -#define _IPOLATES_h -#include "ipolates.h" -#endif extern int decode, latlon, file_append, flush_mode, header; extern int use_scale, dec_scale, bin_scale, wanted_bits, max_bits; @@ -765,17 +167,17 @@ static int ipopt[20] = {-1,-1,0, 0,0,0, 0,0,0, 0}; int f_new_grid_interpolation(ARG1) { -#ifdef USE_SPECTRAL +// #ifdef USE_SPECTRAL char type; int max_wave; -#endif +// #endif if (mode >= -1) { if (strcmp(arg1,"bilinear") == 0) { interpol_type = 0; ipopt[0] = -1; } else if (strcmp(arg1,"bicubic") == 0) { interpol_type = 1; ipopt[0] = 0; } else if (strcmp(arg1,"neighbor") == 0) { interpol_type = 2; ipopt[0] = 1; } else if (strcmp(arg1,"budget") == 0) { interpol_type = 3; ipopt[0] = -1; } -#ifdef USE_SPECTRAL +// #ifdef USE_SPECTRAL else if (sscanf(arg1,"spectral-%c%d", &type, &max_wave) == 2) { if (type == 't' || type == 'T') ipopt[0] = 0; else if (type == 'r' || type == 'R') ipopt[0] = 1; @@ -784,7 +186,7 @@ int f_new_grid_interpolation(ARG1) { ipopt[1] = max_wave; interpol_type = 4; } -#endif +// #endif else if (strcmp(arg1,"neighbor-budget") == 0) { interpol_type = 6; ipopt[0] = -1; } else fatal_error("new_grid_interpolation: unknown type %s", arg1); } @@ -1454,16 +856,16 @@ int f_new_grid(ARG4) { ftn_nout = (int) n_out; if (is_v) { - IPOLATEV(&tmp_interpol_type, ipopt, &gdtnum_in, &(gdt_in[0]), &gdt_in_size, - &(save->gdtnum_out), + ipolatev_grib2_single_field(&tmp_interpol_type, ipopt, &gdtnum_in, &(gdt_in[0]), + &gdt_in_size, &(save->gdtnum_out), &(save->gdt_out[0]), &(save->gdt_out_size), &ftn_npnts, &ftn_nout, &km, &ibi, bitmap, data_in, data_in+ndata, &n_out, save->rlat,save->rlon,save->crot,save->srot, &ibo, save->bitmap_out, save->data_out, save->data_out + n_out, &iret); } else { - IPOLATES(&tmp_interpol_type, ipopt, &gdtnum_in, &(gdt_in[0]), &gdt_in_size, - &(save->gdtnum_out), + ipolates_grib2_single_field(&tmp_interpol_type, ipopt, &gdtnum_in, &(gdt_in[0]), + &gdt_in_size, &(save->gdtnum_out), &(save->gdt_out[0]), &(save->gdt_out_size), &ftn_npnts, &ftn_nout, &km, &ibi, bitmap, data_in, &n_out, save->rlat,save->rlon, &ibo, save->bitmap_out, save->data_out, &iret); @@ -1480,7 +882,7 @@ int f_new_grid(ARG4) { fatal_error_i("IPOLATES failed, error %d",iret); } - // save datda to data_wrt[] + // save data to data_wrt[] if (ibo == 1) { // has a bitmap #pragma omp parallel for private(i) @@ -1591,7 +993,11 @@ int f_new_grid(ARG4) { #endif -#if USE_IPOLATES == 0 +#ifndef USE_IPOLATES +int f_new_grid_vectors(ARG1) { + fprintf(stderr,"IPOLATES package is not installed\n"); + return 1; +} int f_new_grid_interpolation(ARG1) { fprintf(stderr,"IPOLATES package is not installed\n"); return 1; diff --git a/wgrib2/init.c b/wgrib2/init.c index 715f8f97..667bd5e6 100644 --- a/wgrib2/init.c +++ b/wgrib2/init.c @@ -261,7 +261,7 @@ void init_globals(void) { match_extra_fn_n = 0; strncpy(ndates_fmt," %s",4); -#if USE_IPOLATES != 0 +#ifdef USE_IPOLATES wind_rotation = undefined; new_grid_format = grib; #endif diff --git a/wgrib2/ncep_grids.c b/wgrib2/ncep_grids.c index a72ba3e5..7e1cc2ae 100644 --- a/wgrib2/ncep_grids.c +++ b/wgrib2/ncep_grids.c @@ -8,7 +8,7 @@ #include "fnlist.h" -#if USE_IPOLATES > 0 +#ifdef USE_IPOLATES /* * ncep_grids: diff --git a/wgrib2/wgrib2.h b/wgrib2/wgrib2.h index 2d777c8b..4ab5495e 100644 --- a/wgrib2/wgrib2.h +++ b/wgrib2/wgrib2.h @@ -696,4 +696,16 @@ double get_unixtime(int year, int month, int day, int hour, int minute, int seco int JMA_Nb(unsigned char **sec); int JMA_Nr(unsigned char **sec); +void ipolates_grib2_single_field(int *interpol, int *ipopt, int *gdt_in, int *gdttmpl_in, int *gdttmpl_size_in, + int *gdt_out, int *gdttmpl_out, int *gdttmpl_size_out, int *mi, int *mo, int *km, + int *ibi, unsigned char *bitmap, double *data_in, int *n_out, double *rlat, double *rlon, + int *ibo, unsigned char *bitmap_out, double *data_out, int *iret); + +void ipolatev_grib2_single_field(int *interpol, int *ipopt, int *gdt_in, int *gdttmpl_in, int *gdttmpl_size_in, + int *gdt_out, int *gdttmpl_out, int *gdttmpl_size_out, int *mi, int *mo, int *km, + int *ibi, unsigned char *bitmap, double *u_in, double *v_in, int *n_out, double *rlat, double *rlon, + double *crot, double *srot, int *ibo, unsigned char *bitmap_out, + double *u_out, double *v_out, int *iret); + + #endif /* _WGRIB2_H_ */ diff --git a/wgrib2/wgrib2_api.h b/wgrib2/wgrib2_api.h index 637d4773..9cdbae76 100644 --- a/wgrib2/wgrib2_api.h +++ b/wgrib2/wgrib2_api.h @@ -12,7 +12,3 @@ int wgrib2_get_reg_data(float *data, size_t size, int reg); int wgrib2_set_reg(float *data, size_t size, int reg); int wgrib2_free_file(const char *filename); // will delete files that are used/not used -#ifndef _IPOLATES_H -#define _IPOLATES_H -#include "ipolates.h" -#endif