|
| 1 | +/*********************************************************************** |
| 2 | + * GNU Lesser General Public License |
| 3 | + * |
| 4 | + * This file is part of the GFDL FRE NetCDF tools package (FRE-NCTools). |
| 5 | + * |
| 6 | + * FRE-NCtools is free software: you can redistribute it and/or modify it under |
| 7 | + * the terms of the GNU Lesser General Public License as published by |
| 8 | + * the Free Software Foundation, either version 3 of the License, or (at |
| 9 | + * your option) any later version. |
| 10 | + * |
| 11 | + * FRE-NCtools is distributed in the hope that it will be useful, but WITHOUT |
| 12 | + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
| 13 | + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| 14 | + * for more details. |
| 15 | + * |
| 16 | + * You should have received a copy of the GNU Lesser General Public |
| 17 | + * License along with FRE-NCTools. If not, see |
| 18 | + * <http://www.gnu.org/licenses/>. |
| 19 | + **********************************************************************/ |
| 20 | +/* |
| 21 | + Copyright (C) 2011 NOAA Geophysical Fluid Dynamics Lab, Princeton, NJ |
| 22 | +*/ |
| 23 | +#include <stdio.h> |
| 24 | +#include <stdlib.h> |
| 25 | +#include <math.h> |
| 26 | +#include <openacc.h> |
| 27 | +#include "mosaic_util.h" |
| 28 | +#include "interp_gpu.h" |
| 29 | +#include "interp.h" |
| 30 | +#include "create_xgrid.h" |
| 31 | +#include "create_xgrid_gpu.h" |
| 32 | +#include "create_xgrid_utils_gpu.h" |
| 33 | + |
| 34 | +/*------------------------------------------------------------------------------ |
| 35 | + void conserve_interp() |
| 36 | + conservative interpolation through exchange grid. |
| 37 | + Currently only first order interpolation are implemented here. |
| 38 | + ----------------------------------------------------------------------------*/ |
| 39 | +void conserve_interp_gpu(int nx_src, int ny_src, int nx_dst, int ny_dst, double *x_src, |
| 40 | + double *y_src, double *x_dst, double *y_dst, |
| 41 | + double *mask_src, double *data_src, double *data_dst ) |
| 42 | +{ |
| 43 | + Grid_cells_struct_config output_grid_cells; |
| 44 | + Interp_per_input_tile interp_gpu; |
| 45 | + |
| 46 | + int ncells_src = nx_src * ny_src; |
| 47 | + int ncells_dst = nx_dst * ny_dst; |
| 48 | + int ngridpts_src = (nx_src+1) * (ny_src+1); |
| 49 | + int ngridpts_dst = (nx_dst+1) * (ny_dst+1); |
| 50 | + int jstart = 0; |
| 51 | + int jend = ny_src-1; |
| 52 | + |
| 53 | + int *approx_nxcells; approx_nxcells = (int *)malloc(ncells_src*sizeof(int)); |
| 54 | + int *ij2_start; ij2_start = (int *)malloc(ncells_src*sizeof(int)); |
| 55 | + int *ij2_end; ij2_end = (int *)malloc(ncells_src*sizeof(int)); |
| 56 | + |
| 57 | +#pragma acc enter data copyin(x_src[:ngridpts_src], y_src[:ngridpts_src], \ |
| 58 | + x_dst[:ngridpts_dst], y_dst[:ngridpts_dst],\ |
| 59 | + mask_src[:ncells_src]) |
| 60 | + |
| 61 | + get_grid_cell_struct_gpu(nx_dst, ny_dst, x_dst, y_dst, &output_grid_cells); |
| 62 | + |
| 63 | +#pragma acc enter data create(approx_nxcells[:ncells_src], ij2_start[:ncells_src], ij2_end[:ncells_src]) |
| 64 | + |
| 65 | + int upbound_nxcells = get_upbound_nxcells_2dx2d_gpu(nx_src, ny_src, nx_dst, ny_dst, jstart, jend, |
| 66 | + x_src, y_src, x_dst, y_dst, mask_src, &output_grid_cells, |
| 67 | + approx_nxcells, ij2_start, ij2_end); |
| 68 | + |
| 69 | + int nxgrid = create_xgrid_2dx2d_order1_gpu(nx_src, ny_src, nx_dst, ny_dst, jstart, jend, x_src, y_src, |
| 70 | + x_dst, y_dst, upbound_nxcells, mask_src, &output_grid_cells, |
| 71 | + approx_nxcells, ij2_start, ij2_end, &interp_gpu); |
| 72 | + |
| 73 | +#pragma acc exit data copyout(interp_gpu.input_parent_cell_index[:nxgrid], \ |
| 74 | + interp_gpu.output_parent_cell_index[:nxgrid], \ |
| 75 | + interp_gpu.xcell_area[:nxgrid]) |
| 76 | + |
| 77 | + int *xgrid_ij1 = interp_gpu.input_parent_cell_index; |
| 78 | + int *xgrid_ij2 = interp_gpu.output_parent_cell_index; |
| 79 | + double *xgrid_area = interp_gpu.xcell_area; |
| 80 | + double *dst_area; dst_area = (double *)malloc(nx_dst*ny_dst*sizeof(double)); |
| 81 | + |
| 82 | + for(int n=0; n<nx_dst*ny_dst; n++) { |
| 83 | + dst_area[n] = 0.0; |
| 84 | + data_dst[n] = 0.0; |
| 85 | + } |
| 86 | + |
| 87 | + /* The source grid may not cover the destination grid |
| 88 | + so need to sum of exchange grid area to get dst_area |
| 89 | + get_grid_area(&nx_dst, &ny_dst, x_dst, y_dst, dst_area); |
| 90 | + */ |
| 91 | + |
| 92 | + for(int n=0; n<nxgrid; n++) dst_area[xgrid_ij2[n]] += xgrid_area[n]; |
| 93 | + for(int n=0; n<nxgrid; n++) { |
| 94 | + int ij2 = xgrid_ij2[n]; |
| 95 | + int ij1 = xgrid_ij1[n]; |
| 96 | + double area_frac = xgrid_area[n]/dst_area[ij2]; |
| 97 | + data_dst[ij2] += data_src[ij1]*area_frac; |
| 98 | + } |
| 99 | + |
| 100 | + free(xgrid_ij1); |
| 101 | + free(xgrid_ij2); |
| 102 | + free(xgrid_area); |
| 103 | + free(approx_nxcells); |
| 104 | + free(ij2_start); |
| 105 | + free(ij2_end); |
| 106 | + |
| 107 | +}; /* conserve_interp */ |
0 commit comments