Commit 8ab4d09a authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Merge declaration and initialization.

parent 94148dad
...@@ -534,9 +534,7 @@ int set_remapgrids(int filetype, int vlistID, int ngrids, std::vector<bool>& rem ...@@ -534,9 +534,7 @@ int set_remapgrids(int filetype, int vlistID, int ngrids, std::vector<bool>& rem
} }
} }
for ( index = 0; index < ngrids; index++ ) for ( index = 0; index < ngrids; index++ ) if ( remapgrids[index] ) break;
if ( remapgrids[index] ) break;
if ( index == ngrids ) cdoAbort("No remappable grid found!"); if ( index == ngrids ) cdoAbort("No remappable grid found!");
return index; return index;
...@@ -776,10 +774,8 @@ void *Remap(void *argument) ...@@ -776,10 +774,8 @@ void *Remap(void *argument)
int streamID2 = -1; int streamID2 = -1;
int nrecs; int nrecs;
int varID, levelID; int varID, levelID;
size_t gridsize, gridsize2; size_t gridsize2;
int gridID1 = -1, gridID2;
size_t nmiss1, nmiss2; size_t nmiss1, nmiss2;
size_t j;
int r = -1; int r = -1;
int nremaps = 0; int nremaps = 0;
int norm_opt = NORM_OPT_NONE; int norm_opt = NORM_OPT_NONE;
...@@ -835,7 +831,6 @@ void *Remap(void *argument) ...@@ -835,7 +831,6 @@ void *Remap(void *argument)
{ {
operatorInputArg("grid description file or name, remap weights file (SCRIP NetCDF)"); operatorInputArg("grid description file or name, remap weights file (SCRIP NetCDF)");
operatorCheckArgc(2); operatorCheckArgc(2);
gridID2 = cdoDefineGrid(operatorArgv()[0]);
remap_file = operatorArgv()[1]; remap_file = operatorArgv()[1];
} }
else else
...@@ -843,7 +838,6 @@ void *Remap(void *argument) ...@@ -843,7 +838,6 @@ void *Remap(void *argument)
operatorInputArg("grid description file or name"); operatorInputArg("grid description file or name");
if ( operfunc == REMAPDIS && operatorArgc() == 2 ) if ( operfunc == REMAPDIS && operatorArgc() == 2 )
{ {
gridID2 = cdoDefineGrid(operatorArgv()[0]);
int inum = parameter2int(operatorArgv()[1]); int inum = parameter2int(operatorArgv()[1]);
// if ( inum < 1 || inum > 9 ) cdoAbort("Number of nearest neighbors out of range (1-9)!", inum); // if ( inum < 1 || inum > 9 ) cdoAbort("Number of nearest neighbors out of range (1-9)!", inum);
num_neighbors = inum; num_neighbors = inum;
...@@ -851,9 +845,10 @@ void *Remap(void *argument) ...@@ -851,9 +845,10 @@ void *Remap(void *argument)
else else
{ {
operatorCheckArgc(1); operatorCheckArgc(1);
gridID2 = cdoDefineGrid(operatorArgv()[0]);
} }
} }
int gridID2 = cdoDefineGrid(operatorArgv()[0]);
if ( gridInqType(gridID2) == GRID_GENERIC ) cdoAbort("Unsupported target grid type (generic)!"); if ( gridInqType(gridID2) == GRID_GENERIC ) cdoAbort("Unsupported target grid type (generic)!");
...@@ -871,7 +866,7 @@ void *Remap(void *argument) ...@@ -871,7 +866,7 @@ void *Remap(void *argument)
int ngrids = vlistNgrids(vlistID1); int ngrids = vlistNgrids(vlistID1);
std::vector<bool> remapgrids(ngrids); std::vector<bool> remapgrids(ngrids);
int index = set_remapgrids(filetype, vlistID1, ngrids, remapgrids); int index = set_remapgrids(filetype, vlistID1, ngrids, remapgrids);
gridID1 = vlistGrid(vlistID1, index); int gridID1 = vlistGrid(vlistID1, index);
for ( index = 0; index < ngrids; index++ ) for ( index = 0; index < ngrids; index++ )
if ( remapgrids[index] ) if ( remapgrids[index] )
...@@ -979,7 +974,7 @@ void *Remap(void *argument) ...@@ -979,7 +974,7 @@ void *Remap(void *argument)
double *array1 = (double*) Malloc(grid1sizemax*sizeof(double)); double *array1 = (double*) Malloc(grid1sizemax*sizeof(double));
int *imask = (int*) Malloc(grid1sizemax*sizeof(int)); int *imask = (int*) Malloc(grid1sizemax*sizeof(int));
gridsize = gridInqSize(gridID2); size_t gridsize = gridInqSize(gridID2);
double *array2 = (double*) Malloc(gridsize*sizeof(double)); double *array2 = (double*) Malloc(gridsize*sizeof(double));
if ( ! lwrite_remap ) if ( ! lwrite_remap )
...@@ -1119,7 +1114,7 @@ void *Remap(void *argument) ...@@ -1119,7 +1114,7 @@ void *Remap(void *argument)
(gridInqType(gridID1) == GRID_GME || gridInqType(gridID1) == GRID_UNSTRUCTURED) ) (gridInqType(gridID1) == GRID_GME || gridInqType(gridID1) == GRID_UNSTRUCTURED) )
cdoAbort("Bilinear/bicubic interpolation doesn't support unstructured source grids!"); cdoAbort("Bilinear/bicubic interpolation doesn't support unstructured source grids!");
/* initialize grid information for both grids */ // Initialize grid information for both grids
if ( cdoTimer ) timer_start(timer_remap_init); if ( cdoTimer ) timer_start(timer_remap_init);
remap_grids_init(map_type, remap_extrapolate, gridID1, &remaps[r].src_grid, gridID2, &remaps[r].tgt_grid); remap_grids_init(map_type, remap_extrapolate, gridID1, &remaps[r].src_grid, gridID2, &remaps[r].tgt_grid);
if ( cdoTimer ) timer_stop(timer_remap_init); if ( cdoTimer ) timer_stop(timer_remap_init);
...@@ -1130,7 +1125,7 @@ void *Remap(void *argument) ...@@ -1130,7 +1125,7 @@ void *Remap(void *argument)
if ( gridInqType(gridID1) == GRID_GME ) if ( gridInqType(gridID1) == GRID_GME )
{ {
j = 0; size_t j = 0;
for ( size_t i = 0; i < gridsize; i++ ) for ( size_t i = 0; i < gridsize; i++ )
if ( remaps[r].src_grid.vgpm[i] ) imask[j++] = imask[i]; if ( remaps[r].src_grid.vgpm[i] ) imask[j++] = imask[i];
} }
...@@ -1174,7 +1169,7 @@ void *Remap(void *argument) ...@@ -1174,7 +1169,7 @@ void *Remap(void *argument)
if ( gridInqType(gridID1) == GRID_GME ) if ( gridInqType(gridID1) == GRID_GME )
{ {
j = 0; size_t j = 0;
for ( size_t i = 0; i < gridsize; i++ ) for ( size_t i = 0; i < gridsize; i++ )
if ( remaps[r].src_grid.vgpm[i] ) array1[j++] = array1[i]; if ( remaps[r].src_grid.vgpm[i] ) array1[j++] = array1[i];
} }
...@@ -1252,8 +1247,8 @@ void *Remap(void *argument) ...@@ -1252,8 +1247,8 @@ void *Remap(void *argument)
{ {
int nd, ni, ni2, ni3; int nd, ni, ni2, ni3;
gridInqParamGME(gridID2, &nd, &ni, &ni2, &ni3); gridInqParamGME(gridID2, &nd, &ni, &ni2, &ni3);
j = remaps[r].tgt_grid.size;
size_t j = remaps[r].tgt_grid.size;
for ( size_t i = gridsize2; i > 0 ; i-- ) for ( size_t i = gridsize2; i > 0 ; i-- )
if ( remaps[r].tgt_grid.vgpm[i-1] ) array2[i-1] = array2[--j]; if ( remaps[r].tgt_grid.vgpm[i-1] ) array2[i-1] = array2[--j];
......
...@@ -562,8 +562,7 @@ void remap_define_reg2d(int gridID, remapgrid_t *grid) ...@@ -562,8 +562,7 @@ void remap_define_reg2d(int gridID, remapgrid_t *grid)
gridInqXvals(gridID, grid->reg2d_center_lon); gridInqXvals(gridID, grid->reg2d_center_lon);
gridInqYvals(gridID, grid->reg2d_center_lat); gridInqYvals(gridID, grid->reg2d_center_lat);
/* Convert lat/lon units if required */ // Convert lat/lon units if required
char yunits[CDI_MAX_NAME]; yunits[0] = 0; char yunits[CDI_MAX_NAME]; yunits[0] = 0;
cdiGridInqKeyStr(gridID, CDI_KEY_YUNITS, CDI_MAX_NAME, yunits); cdiGridInqKeyStr(gridID, CDI_KEY_YUNITS, CDI_MAX_NAME, yunits);
...@@ -584,9 +583,8 @@ void remap_define_reg2d(int gridID, remapgrid_t *grid) ...@@ -584,9 +583,8 @@ void remap_define_reg2d(int gridID, remapgrid_t *grid)
grid_gen_corners(ny, grid->reg2d_center_lat, grid->reg2d_corner_lat); grid_gen_corners(ny, grid->reg2d_center_lat, grid->reg2d_corner_lat);
grid_check_lat_borders_rad(ny+1, grid->reg2d_corner_lat); grid_check_lat_borders_rad(ny+1, grid->reg2d_corner_lat);
//for ( size_t i = 0; i < nxp1; ++i ) printf("lon %zu %g\n", i, grid->reg2d_corner_lon[i]); // for ( size_t i = 0; i < nxp1; ++i ) printf("lon %zu %g\n", i, grid->reg2d_corner_lon[i]);
//for ( size_t i = 0; i < nyp1; ++i ) printf("lat %zu %g\n", i, grid->reg2d_corner_lat[i]); // for ( size_t i = 0; i < nyp1; ++i ) printf("lat %zu %g\n", i, grid->reg2d_corner_lat[i]);
} }
static static
...@@ -884,9 +882,7 @@ void remap_grids_init(int map_type, bool lextrapolate, int gridID1, remapgrid_t ...@@ -884,9 +882,7 @@ void remap_grids_init(int map_type, bool lextrapolate, int gridID1, remapgrid_t
{ {
cell_bounding_boxes(src_grid, REMAP_GRID_BASIS_SRC); cell_bounding_boxes(src_grid, REMAP_GRID_BASIS_SRC);
cell_bounding_boxes(tgt_grid, REMAP_GRID_BASIS_TGT); cell_bounding_boxes(tgt_grid, REMAP_GRID_BASIS_TGT);
/* // Set up and assign address ranges to search bins in order to further restrict later searches
Set up and assign address ranges to search bins in order to further restrict later searches
*/
calc_lat_bins(src_grid, tgt_grid, map_type); calc_lat_bins(src_grid, tgt_grid, map_type);
} }
...@@ -1515,10 +1511,6 @@ void remap_stat(int remap_order, remapgrid_t src_grid, remapgrid_t tgt_grid, rem ...@@ -1515,10 +1511,6 @@ void remap_stat(int remap_order, remapgrid_t src_grid, remapgrid_t tgt_grid, rem
void remap_gradients(remapgrid_t grid, const double *restrict array, double *restrict grad_lat, void remap_gradients(remapgrid_t grid, const double *restrict array, double *restrict grad_lat,
double *restrict grad_lon, double *restrict grad_latlon) double *restrict grad_lon, double *restrict grad_latlon)
{ {
size_t i, j, ip1, im1, jp1, jm1, in, is, ie, iw, ine, inw, ise, isw;
double delew, delns;
double grad_lat_zero, grad_lon_zero;
if ( grid.rank != 2 ) if ( grid.rank != 2 )
cdoAbort("Internal problem (remap_gradients), grid rank = %d!", grid.rank); cdoAbort("Internal problem (remap_gradients), grid rank = %d!", grid.rank);
...@@ -1528,8 +1520,7 @@ void remap_gradients(remapgrid_t grid, const double *restrict array, double *res ...@@ -1528,8 +1520,7 @@ void remap_gradients(remapgrid_t grid, const double *restrict array, double *res
#if defined(_OPENMP) #if defined(_OPENMP)
#pragma omp parallel for default(none) \ #pragma omp parallel for default(none) \
shared(grid_size, grad_lat, grad_lon, grad_latlon, grid, nx, ny, array) \ shared(grid_size, grad_lat, grad_lon, grad_latlon, grid, nx, ny, array)
private(i, j, ip1, im1, jp1, jm1, in, is, ie, iw, ine, inw, ise, isw, delew, delns, grad_lat_zero, grad_lon_zero)
#endif #endif
for ( size_t n = 0; n < grid_size; ++n ) for ( size_t n = 0; n < grid_size; ++n )
{ {
...@@ -1539,77 +1530,47 @@ void remap_gradients(remapgrid_t grid, const double *restrict array, double *res ...@@ -1539,77 +1530,47 @@ void remap_gradients(remapgrid_t grid, const double *restrict array, double *res
if ( grid.mask[n] ) if ( grid.mask[n] )
{ {
delew = HALF; double delew = HALF;
delns = HALF; double delns = HALF;
j = n/nx + 1; size_t j = n/nx + 1;
i = n - (j-1)*nx + 1; size_t i = n - (j-1)*nx + 1;
ip1 = i + 1; size_t ip1 = i + 1;
im1 = i - 1; size_t im1 = i - 1;
jp1 = j + 1; size_t jp1 = j + 1;
jm1 = j - 1; size_t jm1 = j - 1;
if ( ip1 > nx ) ip1 = ip1 - nx; if ( ip1 > nx ) ip1 = ip1 - nx;
if ( im1 < 1 ) im1 = nx; if ( im1 < 1 ) im1 = nx;
if ( jp1 > ny ) if ( jp1 > ny ) { jp1 = j; delns = ONE; }
{ if ( jm1 < 1 ) { jm1 = j; delns = ONE; }
jp1 = j;
delns = ONE; size_t in = (jp1-1)*nx + i - 1;
} size_t is = (jm1-1)*nx + i - 1;
if ( jm1 < 1 ) size_t ie = (j -1)*nx + ip1 - 1;
{ size_t iw = (j -1)*nx + im1 - 1;
jm1 = j;
delns = ONE; size_t ine = (jp1-1)*nx + ip1 - 1;
} size_t inw = (jp1-1)*nx + im1 - 1;
size_t ise = (jm1-1)*nx + ip1 - 1;
in = (jp1-1)*nx + i - 1; size_t isw = (jm1-1)*nx + im1 - 1;
is = (jm1-1)*nx + i - 1;
ie = (j -1)*nx + ip1 - 1; // Compute i-gradient
iw = (j -1)*nx + im1 - 1; if ( ! grid.mask[ie] ) { ie = n; delew = ONE; }
if ( ! grid.mask[iw] ) { iw = n; delew = ONE; }
ine = (jp1-1)*nx + ip1 - 1;
inw = (jp1-1)*nx + im1 - 1;
ise = (jm1-1)*nx + ip1 - 1;
isw = (jm1-1)*nx + im1 - 1;
/* Compute i-gradient */
if ( ! grid.mask[ie] )
{
ie = n;
delew = ONE;
}
if ( ! grid.mask[iw] )
{
iw = n;
delew = ONE;
}
grad_lat[n] = delew*(array[ie] - array[iw]); grad_lat[n] = delew*(array[ie] - array[iw]);
/* Compute j-gradient */ // Compute j-gradient
if ( ! grid.mask[in] ) { in = n; delns = ONE; }
if ( ! grid.mask[in] ) if ( ! grid.mask[is] ) { is = n; delns = ONE; }
{
in = n;
delns = ONE;
}
if ( ! grid.mask[is] )
{
is = n;
delns = ONE;
}
grad_lon[n] = delns*(array[in] - array[is]); grad_lon[n] = delns*(array[in] - array[is]);
/* Compute ij-gradient */ // Compute ij-gradient
delew = HALF; delew = HALF;
if ( jp1 == j || jm1 == j ) delns = (jp1 == j || jm1 == j) ? ONE : HALF;
delns = ONE;
else
delns = HALF;
if ( ! grid.mask[ine] ) if ( ! grid.mask[ine] )
{ {
...@@ -1657,7 +1618,7 @@ void remap_gradients(remapgrid_t grid, const double *restrict array, double *res ...@@ -1657,7 +1618,7 @@ void remap_gradients(remapgrid_t grid, const double *restrict array, double *res
} }
} }
grad_lat_zero = delew*(array[ine] - array[inw]); double grad_lat_zero = delew*(array[ine] - array[inw]);
if ( ! grid.mask[ise] ) if ( ! grid.mask[ise] )
{ {
...@@ -1705,7 +1666,7 @@ void remap_gradients(remapgrid_t grid, const double *restrict array, double *res ...@@ -1705,7 +1666,7 @@ void remap_gradients(remapgrid_t grid, const double *restrict array, double *res
} }
} }
grad_lon_zero = delew*(array[ise] - array[isw]); double grad_lon_zero = delew*(array[ise] - array[isw]);
grad_latlon[n] = delns*(grad_lat_zero - grad_lon_zero); grad_latlon[n] = delns*(grad_lat_zero - grad_lon_zero);
} }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment