FMS 2025.01.02-dev
Flexible Modeling System
Loading...
Searching...
No Matches
horiz_interp_conserve_xgrid.c
Go to the documentation of this file.
1/***********************************************************************
2 * GNU Lesser General Public License
3 *
4 * This file is part of the GFDL Flexible Modeling System (FMS).
5 *
6 * FMS 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 * FMS 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 FMS. If not, see <http://www.gnu.org/licenses/>.
18 **********************************************************************/
19#include <stdlib.h>
20#include <stdio.h>
21#include <math.h>
22#include "grid_utils.h"
23#include "tree_utils.h"
24#include "horiz_interp_conserve_xgrid.h"
25#include "constant.h"
26
27#if defined(_OPENMP)
28#include <omp.h>
29#endif
30
31/** \file
32 * \ingroup mosaic
33 * \brief Grid creation and calculation functions for use in @ref mosaic_mod
34 * /
35
36/*******************************************************************************
37 void create_xgrid_1dx2d_order1
38 This routine generate exchange grids between two grids for the first order
39 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
40 and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds.
41*******************************************************************************/
42int create_xgrid_1dx2d_order1_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
43 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
44 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area)
45{
46 int nxgrid;
47
48 nxgrid = create_xgrid_1dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in,
49 i_in, j_in, i_out, j_out, xgrid_area);
50 return nxgrid;
51
52}
53
54int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in,
55 const double *lat_in, const double *lon_out, const double *lat_out,
56 const double *mask_in, int *i_in, int *j_in, int *i_out,
57 int *j_out, double *xgrid_area)
58{
59
60 int nx1, ny1, nx2, ny2, nx1p, nx2p;
61 int i1, j1, i2, j2, nxgrid;
62 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
63 double *area_in, *area_out, min_area;
64 double *tmpx, *tmpy;
65
66 nx1 = *nlon_in;
67 ny1 = *nlat_in;
68 nx2 = *nlon_out;
69 ny2 = *nlat_out;
70
71 nxgrid = 0;
72 nx1p = nx1 + 1;
73 nx2p = nx2 + 1;
74
75 area_in = (double *)malloc(nx1*ny1*sizeof(double));
76 area_out = (double *)malloc(nx2*ny2*sizeof(double));
77 tmpx = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
78 tmpy = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
79 for(j1=0; j1<=ny1; j1++) for(i1=0; i1<=nx1; i1++) {
80 tmpx[j1*nx1p+i1] = lon_in[i1];
81 tmpy[j1*nx1p+i1] = lat_in[j1];
82 }
83 /* This is just a temporary fix to solve the issue that there is one point in zonal direction */
84 if(nx1 > 1)
85 get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in);
86 else
87 get_grid_area_no_adjust(nlon_in, nlat_in, tmpx, tmpy, area_in);
88
89 get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out);
90 free(tmpx);
91 free(tmpy);
92
93 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
94
95 ll_lon = lon_in[i1]; ll_lat = lat_in[j1];
96 ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1];
97 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
98 int n_in, n_out;
99 double Xarea;
100
101 y_in[0] = lat_out[j2*nx2p+i2];
102 y_in[1] = lat_out[j2*nx2p+i2+1];
103 y_in[2] = lat_out[(j2+1)*nx2p+i2+1];
104 y_in[3] = lat_out[(j2+1)*nx2p+i2];
105 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
106 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
107 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
108 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
109
110 x_in[0] = lon_out[j2*nx2p+i2];
111 x_in[1] = lon_out[j2*nx2p+i2+1];
112 x_in[2] = lon_out[(j2+1)*nx2p+i2+1];
113 x_in[3] = lon_out[(j2+1)*nx2p+i2];
114 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
115
116 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
117 Xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
118 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
119 if( Xarea/min_area > AREA_RATIO_THRESH ) {
120 xgrid_area[nxgrid] = Xarea;
121 i_in[nxgrid] = i1;
122 j_in[nxgrid] = j1;
123 i_out[nxgrid] = i2;
124 j_out[nxgrid] = j2;
125 ++nxgrid;
126 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
127 }
128 }
129 }
130 }
131
132 free(area_in);
133 free(area_out);
134
135 return nxgrid;
136
137} /* create_xgrid_1dx2d_order1 */
138
139
140/*******************************************************************************
141 void create_xgrid_1dx2d_order1_ug
142 This routine generate exchange grids between two grids for the first order
143 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
144 and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds.
145*******************************************************************************/
146int create_xgrid_1dx2d_order1_ug_(const int *nlon_in, const int *nlat_in, const int *npts_out,
147 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
148 const double *mask_in, int *i_in, int *j_in, int *l_out, double *xgrid_area)
149{
150 int nxgrid;
151
152 nxgrid = create_xgrid_1dx2d_order1_ug(nlon_in, nlat_in, npts_out, lon_in, lat_in, lon_out, lat_out, mask_in,
153 i_in, j_in, l_out, xgrid_area);
154 return nxgrid;
155
156}
157
158int create_xgrid_1dx2d_order1_ug(const int *nlon_in, const int *nlat_in, const int *npts_out, const double *lon_in,
159 const double *lat_in, const double *lon_out, const double *lat_out,
160 const double *mask_in, int *i_in, int *j_in, int *l_out, double *xgrid_area)
161{
162
163 int nx1, ny1, nx1p, nv, npts2;
164 int i1, j1, l2, nxgrid;
165 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
166 double *area_in, *area_out, min_area;
167 double *tmpx, *tmpy;
168
169 nx1 = *nlon_in;
170 ny1 = *nlat_in;
171 nv = 4;
172 npts2 = *npts_out;
173
174 nxgrid = 0;
175 nx1p = nx1 + 1;
176
177 area_in = (double *)malloc(nx1*ny1*sizeof(double));
178 area_out = (double *)malloc(npts2*sizeof(double));
179 tmpx = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
180 tmpy = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
181 for(j1=0; j1<=ny1; j1++) for(i1=0; i1<=nx1; i1++) {
182 tmpx[j1*nx1p+i1] = lon_in[i1];
183 tmpy[j1*nx1p+i1] = lat_in[j1];
184 }
185 /* This is just a temporary fix to solve the issue that there is one point in zonal direction */
186 if(nx1 > 1)
187 get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in);
188 else
189 get_grid_area_no_adjust(nlon_in, nlat_in, tmpx, tmpy, area_in);
190
191 get_grid_area_ug(npts_out, lon_out, lat_out, area_out);
192 free(tmpx);
193 free(tmpy);
194
195 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
196
197 ll_lon = lon_in[i1]; ll_lat = lat_in[j1];
198 ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1];
199 for(l2=0; l2<npts2; l2++) {
200 int n_in, n_out;
201 double Xarea;
202
203 y_in[0] = lat_out[l2*nv];
204 y_in[1] = lat_out[l2*nv+1];
205 y_in[2] = lat_out[l2*nv+2];
206 y_in[3] = lat_out[l2*nv+3];
207 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
208 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
209 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
210 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
211
212 x_in[0] = lon_out[l2*nv];
213 x_in[1] = lon_out[l2*nv+1];
214 x_in[2] = lon_out[l2*nv+2];
215 x_in[3] = lon_out[l2*nv+3];
216 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
217
218 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
219 Xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
220 min_area = min(area_in[j1*nx1+i1], area_out[l2]);
221 if( Xarea/min_area > AREA_RATIO_THRESH ) {
222 xgrid_area[nxgrid] = Xarea;
223 i_in[nxgrid] = i1;
224 j_in[nxgrid] = j1;
225 l_out[nxgrid] = l2;
226 ++nxgrid;
227 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
228 }
229 }
230 }
231 }
232
233 free(area_in);
234 free(area_out);
235
236 return nxgrid;
237
238} /* create_xgrid_1dx2d_order1_ug */
239
240/********************************************************************************
241 void create_xgrid_1dx2d_order2
242 This routine generate exchange grids between two grids for the second order
243 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
244 and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds.
245********************************************************************************/
246int create_xgrid_1dx2d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
247 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
248 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
249 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
250{
251 int nxgrid;
252 nxgrid = create_xgrid_1dx2d_order2(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in,
253 j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat);
254 return nxgrid;
255
256}
257int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
258 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
259 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
260 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
261{
262
263 int nx1, ny1, nx2, ny2, nx1p, nx2p;
264 int i1, j1, i2, j2, nxgrid;
265 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
266 double *area_in, *area_out, min_area;
267 double *tmpx, *tmpy;
268
269 nx1 = *nlon_in;
270 ny1 = *nlat_in;
271 nx2 = *nlon_out;
272 ny2 = *nlat_out;
273
274 nxgrid = 0;
275 nx1p = nx1 + 1;
276 nx2p = nx2 + 1;
277
278 area_in = (double *)malloc(nx1*ny1*sizeof(double));
279 area_out = (double *)malloc(nx2*ny2*sizeof(double));
280 tmpx = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
281 tmpy = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
282 for(j1=0; j1<=ny1; j1++) for(i1=0; i1<=nx1; i1++) {
283 tmpx[j1*nx1p+i1] = lon_in[i1];
284 tmpy[j1*nx1p+i1] = lat_in[j1];
285 }
286 get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in);
287 get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out);
288 free(tmpx);
289 free(tmpy);
290
291 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
292
293 ll_lon = lon_in[i1]; ll_lat = lat_in[j1];
294 ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1];
295 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
296 int n_in, n_out;
297 double xarea, lon_in_avg;
298
299 y_in[0] = lat_out[j2*nx2p+i2];
300 y_in[1] = lat_out[j2*nx2p+i2+1];
301 y_in[2] = lat_out[(j2+1)*nx2p+i2+1];
302 y_in[3] = lat_out[(j2+1)*nx2p+i2];
303 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
304 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
305 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
306 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
307
308 x_in[0] = lon_out[j2*nx2p+i2];
309 x_in[1] = lon_out[j2*nx2p+i2+1];
310 x_in[2] = lon_out[(j2+1)*nx2p+i2+1];
311 x_in[3] = lon_out[(j2+1)*nx2p+i2];
312 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
313 lon_in_avg = avgval_double(n_in, x_in);
314
315 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
316 xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
317 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
318 if(xarea/min_area > AREA_RATIO_THRESH ) {
319 xgrid_area[nxgrid] = xarea;
320 xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg);
321 xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out );
322 i_in[nxgrid] = i1;
323 j_in[nxgrid] = j1;
324 i_out[nxgrid] = i2;
325 j_out[nxgrid] = j2;
326 ++nxgrid;
327 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
328 }
329 }
330 }
331 }
332 free(area_in);
333 free(area_out);
334
335 return nxgrid;
336
337} /* create_xgrid_1dx2d_order2 */
338
339/*******************************************************************************
340 void create_xgrid_2dx1d_order1
341 This routine generate exchange grids between two grids for the first order
342 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
343 and lon_out,lat_out are 1-D grid bounds, lon_in,lat_in are geographic grid location of grid cell bounds.
344 mask is on grid lon_in/lat_in.
345*******************************************************************************/
346int create_xgrid_2dx1d_order1_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
347 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
348 const double *mask_in, int *i_in, int *j_in, int *i_out,
349 int *j_out, double *xgrid_area)
350{
351 int nxgrid;
352
353 nxgrid = create_xgrid_2dx1d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in,
354 i_in, j_in, i_out, j_out, xgrid_area);
355 return nxgrid;
356
357}
358int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in,
359 const double *lat_in, const double *lon_out, const double *lat_out,
360 const double *mask_in, int *i_in, int *j_in, int *i_out,
361 int *j_out, double *xgrid_area)
362{
363
364 int nx1, ny1, nx2, ny2, nx1p, nx2p;
365 int i1, j1, i2, j2, nxgrid;
366 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
367 double *area_in, *area_out, min_area;
368 double *tmpx, *tmpy;
369 int n_in, n_out;
370 double Xarea;
371
372
373 nx1 = *nlon_in;
374 ny1 = *nlat_in;
375 nx2 = *nlon_out;
376 ny2 = *nlat_out;
377
378 nxgrid = 0;
379 nx1p = nx1 + 1;
380 nx2p = nx2 + 1;
381 area_in = (double *)malloc(nx1*ny1*sizeof(double));
382 area_out = (double *)malloc(nx2*ny2*sizeof(double));
383 tmpx = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double));
384 tmpy = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double));
385 for(j2=0; j2<=ny2; j2++) for(i2=0; i2<=nx2; i2++) {
386 tmpx[j2*nx2p+i2] = lon_out[i2];
387 tmpy[j2*nx2p+i2] = lat_out[j2];
388 }
389 get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in);
390 get_grid_area(nlon_out, nlat_out, tmpx, tmpy, area_out);
391
392 free(tmpx);
393 free(tmpy);
394
395 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
396
397 ll_lon = lon_out[i2]; ll_lat = lat_out[j2];
398 ur_lon = lon_out[i2+1]; ur_lat = lat_out[j2+1];
399 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
400
401 y_in[0] = lat_in[j1*nx1p+i1];
402 y_in[1] = lat_in[j1*nx1p+i1+1];
403 y_in[2] = lat_in[(j1+1)*nx1p+i1+1];
404 y_in[3] = lat_in[(j1+1)*nx1p+i1];
405 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
406 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
407 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
408 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
409
410 x_in[0] = lon_in[j1*nx1p+i1];
411 x_in[1] = lon_in[j1*nx1p+i1+1];
412 x_in[2] = lon_in[(j1+1)*nx1p+i1+1];
413 x_in[3] = lon_in[(j1+1)*nx1p+i1];
414
415 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
416
417 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
418 Xarea = poly_area ( x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
419 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
420 if( Xarea/min_area > AREA_RATIO_THRESH ) {
421 xgrid_area[nxgrid] = Xarea;
422 i_in[nxgrid] = i1;
423 j_in[nxgrid] = j1;
424 i_out[nxgrid] = i2;
425 j_out[nxgrid] = j2;
426 ++nxgrid;
427 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
428 }
429 }
430 }
431 }
432
433 free(area_in);
434 free(area_out);
435
436 return nxgrid;
437
438} /* create_xgrid_2dx1d_order1 */
439
440
441/********************************************************************************
442 void create_xgrid_2dx1d_order2
443 This routine generate exchange grids between two grids for the second order
444 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
445 and lon_out,lat_out are 1-D grid bounds, lon_in,lat_in are geographic grid location of grid cell bounds.
446 mask is on grid lon_in/lat_in.
447********************************************************************************/
448int create_xgrid_2dx1d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
449 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
450 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
451 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
452{
453 int nxgrid;
454 nxgrid = create_xgrid_2dx1d_order2(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in,
455 j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat);
456 return nxgrid;
457
458}
459
460int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
461 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
462 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
463 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
464{
465
466 int nx1, ny1, nx2, ny2, nx1p, nx2p;
467 int i1, j1, i2, j2, nxgrid;
468 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
469 double *tmpx, *tmpy;
470 double *area_in, *area_out, min_area;
471 double lon_in_avg;
472 int n_in, n_out;
473 double xarea;
474
475
476 nx1 = *nlon_in;
477 ny1 = *nlat_in;
478 nx2 = *nlon_out;
479 ny2 = *nlat_out;
480
481 nxgrid = 0;
482 nx1p = nx1 + 1;
483 nx2p = nx2 + 1;
484
485 area_in = (double *)malloc(nx1*ny1*sizeof(double));
486 area_out = (double *)malloc(nx2*ny2*sizeof(double));
487 tmpx = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double));
488 tmpy = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double));
489 for(j2=0; j2<=ny2; j2++) for(i2=0; i2<=nx2; i2++) {
490 tmpx[j2*nx2p+i2] = lon_out[i2];
491 tmpy[j2*nx2p+i2] = lat_out[j2];
492 }
493 get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in);
494 get_grid_area(nlon_out, nlat_out, tmpx, tmpy, area_out);
495
496 free(tmpx);
497 free(tmpy);
498
499 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
500
501 ll_lon = lon_out[i2]; ll_lat = lat_out[j2];
502 ur_lon = lon_out[i2+1]; ur_lat = lat_out[j2+1];
503 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
504
505 y_in[0] = lat_in[j1*nx1p+i1];
506 y_in[1] = lat_in[j1*nx1p+i1+1];
507 y_in[2] = lat_in[(j1+1)*nx1p+i1+1];
508 y_in[3] = lat_in[(j1+1)*nx1p+i1];
509 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
510 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
511 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
512 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
513
514 x_in[0] = lon_in[j1*nx1p+i1];
515 x_in[1] = lon_in[j1*nx1p+i1+1];
516 x_in[2] = lon_in[(j1+1)*nx1p+i1+1];
517 x_in[3] = lon_in[(j1+1)*nx1p+i1];
518
519 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
520 lon_in_avg = avgval_double(n_in, x_in);
521
522 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
523 xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
524 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
525 if(xarea/min_area > AREA_RATIO_THRESH ) {
526 xgrid_area[nxgrid] = xarea;
527 xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg);
528 xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out );
529 i_in[nxgrid] = i1;
530 j_in[nxgrid] = j1;
531 i_out[nxgrid] = i2;
532 j_out[nxgrid] = j2;
533 ++nxgrid;
534 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
535 }
536 }
537 }
538 }
539
540 free(area_in);
541 free(area_out);
542
543 return nxgrid;
544
545} /* create_xgrid_2dx1d_order2 */
546
547/*******************************************************************************
548 void create_xgrid_2DX2D_order1
549 This routine generate exchange grids between two grids for the first order
550 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
551 and lon_in,lat_in, lon_out,lat_out are geographic grid location of grid cell bounds.
552 mask is on grid lon_in/lat_in.
553*******************************************************************************/
554int create_xgrid_2dx2d_order1_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
555 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
556 const double *mask_in, int *i_in, int *j_in, int *i_out,
557 int *j_out, double *xgrid_area)
558{
559 int nxgrid;
560
561 nxgrid = create_xgrid_2dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in,
562 i_in, j_in, i_out, j_out, xgrid_area);
563 return nxgrid;
564
565}
566int create_xgrid_2dx2d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
567 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
568 const double *mask_in, int *i_in, int *j_in, int *i_out,
569 int *j_out, double *xgrid_area)
570{
571
572 int nx1, nx2, ny1, ny2, nx1p, nx2p, nxgrid;
573 double *area_in, *area_out;
574 int nblocks =1;
575 int *istart2=NULL, *iend2=NULL;
576 int npts_left, nblks_left, pos, m, npts_my, ij;
577 double *lon_out_min_list,*lon_out_max_list,*lon_out_avg,*lat_out_min_list,*lat_out_max_list;
578 double *lon_out_list, *lat_out_list;
579 int *pnxgrid=NULL, *pstart;
580 int *pi_in=NULL, *pj_in=NULL, *pi_out=NULL, *pj_out=NULL;
581 double *pxgrid_area=NULL;
582 int *n2_list;
583 int nthreads, nxgrid_block_max;
584
585 nx1 = *nlon_in;
586 ny1 = *nlat_in;
587 nx2 = *nlon_out;
588 ny2 = *nlat_out;
589 nx1p = nx1 + 1;
590 nx2p = nx2 + 1;
591
592 area_in = (double *)malloc(nx1*ny1*sizeof(double));
593 area_out = (double *)malloc(nx2*ny2*sizeof(double));
594 get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in);
595 get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out);
596
597 nthreads = 1;
598#if defined(_OPENMP)
599#pragma omp parallel
600 nthreads = omp_get_num_threads();
601#endif
602
603 nblocks = nthreads;
604
605 istart2 = (int *)malloc(nblocks*sizeof(int));
606 iend2 = (int *)malloc(nblocks*sizeof(int));
607
608 pstart = (int *)malloc(nblocks*sizeof(int));
609 pnxgrid = (int *)malloc(nblocks*sizeof(int));
610
611 nxgrid_block_max = MAXXGRID/nblocks;
612
613 for(m=0; m<nblocks; m++) {
614 pnxgrid[m] = 0;
615 pstart[m] = m*nxgrid_block_max;
616 }
617
618 if(nblocks == 1) {
619 pi_in = i_in;
620 pj_in = j_in;
621 pi_out = i_out;
622 pj_out = j_out;
623 pxgrid_area = xgrid_area;
624 }
625 else {
626 pi_in = (int *)malloc(MAXXGRID*sizeof(int));
627 pj_in = (int *)malloc(MAXXGRID*sizeof(int));
628 pi_out = (int *)malloc(MAXXGRID*sizeof(int));
629 pj_out = (int *)malloc(MAXXGRID*sizeof(int));
630 pxgrid_area = (double *)malloc(MAXXGRID*sizeof(double));
631 }
632
633 npts_left = nx2*ny2;
634 nblks_left = nblocks;
635 pos = 0;
636 for(m=0; m<nblocks; m++) {
637 istart2[m] = pos;
638 npts_my = npts_left/nblks_left;
639 iend2[m] = istart2[m] + npts_my - 1;
640 pos = iend2[m] + 1;
641 npts_left -= npts_my;
642 nblks_left--;
643 }
644
645 lon_out_min_list = (double *)malloc(nx2*ny2*sizeof(double));
646 lon_out_max_list = (double *)malloc(nx2*ny2*sizeof(double));
647 lat_out_min_list = (double *)malloc(nx2*ny2*sizeof(double));
648 lat_out_max_list = (double *)malloc(nx2*ny2*sizeof(double));
649 lon_out_avg = (double *)malloc(nx2*ny2*sizeof(double));
650 n2_list = (int *)malloc(nx2*ny2*sizeof(int));
651 lon_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double));
652 lat_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double));
653#if defined(_OPENMP)
654#pragma omp parallel for default(none) shared(nx2,ny2,nx2p,lon_out,lat_out,lat_out_min_list, \
655 lat_out_max_list,lon_out_min_list,lon_out_max_list, \
656 lon_out_avg,n2_list,lon_out_list,lat_out_list)
657#endif
658 for(ij=0; ij<nx2*ny2; ij++){
659 int i2, j2, n, n0, n1, n2, n3, n2_in, l;
660 double x2_in[MV], y2_in[MV];
661 i2 = ij%nx2;
662 j2 = ij/nx2;
663 n = j2*nx2+i2;
664 n0 = j2*nx2p+i2; n1 = j2*nx2p+i2+1;
665 n2 = (j2+1)*nx2p+i2+1; n3 = (j2+1)*nx2p+i2;
666 x2_in[0] = lon_out[n0]; y2_in[0] = lat_out[n0];
667 x2_in[1] = lon_out[n1]; y2_in[1] = lat_out[n1];
668 x2_in[2] = lon_out[n2]; y2_in[2] = lat_out[n2];
669 x2_in[3] = lon_out[n3]; y2_in[3] = lat_out[n3];
670
671 lat_out_min_list[n] = minval_double(4, y2_in);
672 lat_out_max_list[n] = maxval_double(4, y2_in);
673 n2_in = fix_lon(x2_in, y2_in, 4, M_PI);
674 if(n2_in > MAX_V) error_handler("create_xgrid.c: n2_in is greater than MAX_V");
675 lon_out_min_list[n] = minval_double(n2_in, x2_in);
676 lon_out_max_list[n] = maxval_double(n2_in, x2_in);
677 lon_out_avg[n] = avgval_double(n2_in, x2_in);
678 n2_list[n] = n2_in;
679 for(l=0; l<n2_in; l++) {
680 lon_out_list[n*MAX_V+l] = x2_in[l];
681 lat_out_list[n*MAX_V+l] = y2_in[l];
682 }
683 }
684
685 nxgrid = 0;
686
687#if defined(_OPENMP)
688#pragma omp parallel for default(none) shared(nblocks,nx1,ny1,nx1p,mask_in,lon_in,lat_in, \
689 istart2,iend2,nx2,lat_out_min_list,lat_out_max_list, \
690 n2_list,lon_out_list,lat_out_list,lon_out_min_list, \
691 lon_out_max_list,lon_out_avg,area_in,area_out, \
692 pxgrid_area,pnxgrid,pi_in,pj_in,pi_out,pj_out,pstart,nthreads)
693#endif
694 for(m=0; m<nblocks; m++) {
695 int i1, j1, ij;
696 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
697 int n0, n1, n2, n3, l,n1_in;
698 double lat_in_min,lat_in_max,lon_in_min,lon_in_max,lon_in_avg;
699 double x1_in[MV], y1_in[MV], x_out[MV], y_out[MV];
700
701 n0 = j1*nx1p+i1; n1 = j1*nx1p+i1+1;
702 n2 = (j1+1)*nx1p+i1+1; n3 = (j1+1)*nx1p+i1;
703 x1_in[0] = lon_in[n0]; y1_in[0] = lat_in[n0];
704 x1_in[1] = lon_in[n1]; y1_in[1] = lat_in[n1];
705 x1_in[2] = lon_in[n2]; y1_in[2] = lat_in[n2];
706 x1_in[3] = lon_in[n3]; y1_in[3] = lat_in[n3];
707 lat_in_min = minval_double(4, y1_in);
708 lat_in_max = maxval_double(4, y1_in);
709 n1_in = fix_lon(x1_in, y1_in, 4, M_PI);
710 lon_in_min = minval_double(n1_in, x1_in);
711 lon_in_max = maxval_double(n1_in, x1_in);
712 lon_in_avg = avgval_double(n1_in, x1_in);
713 for(ij=istart2[m]; ij<=iend2[m]; ij++) {
714 int n_out, i2, j2, n2_in;
715 double xarea, dx, lon_out_min, lon_out_max;
716 double x2_in[MAX_V], y2_in[MAX_V];
717
718 i2 = ij%nx2;
719 j2 = ij/nx2;
720
721 if(lat_out_min_list[ij] >= lat_in_max || lat_out_max_list[ij] <= lat_in_min ) continue;
722 /* adjust x2_in according to lon_in_avg*/
723 n2_in = n2_list[ij];
724 for(l=0; l<n2_in; l++) {
725 x2_in[l] = lon_out_list[ij*MAX_V+l];
726 y2_in[l] = lat_out_list[ij*MAX_V+l];
727 }
728 lon_out_min = lon_out_min_list[ij];
729 lon_out_max = lon_out_max_list[ij];
730 dx = lon_out_avg[ij] - lon_in_avg;
731 if(dx < -M_PI ) {
732 lon_out_min += TPI;
733 lon_out_max += TPI;
734 for (l=0; l<n2_in; l++) x2_in[l] += TPI;
735 }
736 else if (dx > M_PI) {
737 lon_out_min -= TPI;
738 lon_out_max -= TPI;
739 for (l=0; l<n2_in; l++) x2_in[l] -= TPI;
740 }
741
742 /* x2_in should in the same range as x1_in after lon_fix, so no need to
743 consider cyclic condition
744 */
745 if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue;
746 if ( (n_out = clip_2dx2d( x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out )) > 0) {
747 double min_area;
748 int nn;
749 xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
750 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
751 if( xarea/min_area > AREA_RATIO_THRESH ) {
752 pnxgrid[m]++;
753 if(pnxgrid[m]>= MAXXGRID/nthreads)
754 error_handler("nxgrid is greater than MAXXGRID/nthreads, increase MAXXGRID, decrease nthreads, or increase number of MPI ranks");
755 nn = pstart[m] + pnxgrid[m]-1;
756
757 pxgrid_area[nn] = xarea;
758 pi_in[nn] = i1;
759 pj_in[nn] = j1;
760 pi_out[nn] = i2;
761 pj_out[nn] = j2;
762 }
763
764 }
765
766 }
767 }
768 }
769
770 /*copy data if nblocks > 1 */
771 if(nblocks == 1) {
772 nxgrid = pnxgrid[0];
773 pi_in = NULL;
774 pj_in = NULL;
775 pi_out = NULL;
776 pj_out = NULL;
777 pxgrid_area = NULL;
778 }
779 else {
780 int nn, i;
781 nxgrid = 0;
782 for(m=0; m<nblocks; m++) {
783 for(i=0; i<pnxgrid[m]; i++) {
784 nn = pstart[m] + i;
785 i_in[nxgrid] = pi_in[nn];
786 j_in[nxgrid] = pj_in[nn];
787 i_out[nxgrid] = pi_out[nn];
788 j_out[nxgrid] = pj_out[nn];
789 xgrid_area[nxgrid] = pxgrid_area[nn];
790 nxgrid++;
791 }
792 }
793 free(pi_in);
794 free(pj_in);
795 free(pi_out);
796 free(pj_out);
797 free(pxgrid_area);
798 }
799
800 free(area_in);
801 free(area_out);
802 free(lon_out_min_list);
803 free(lon_out_max_list);
804 free(lat_out_min_list);
805 free(lat_out_max_list);
806 free(lon_out_avg);
807 free(n2_list);
808 free(lon_out_list);
809 free(lat_out_list);
810
811 return nxgrid;
812
813}/* get_xgrid_2Dx2D_order1 */
814
815/********************************************************************************
816 void create_xgrid_2dx1d_order2
817 This routine generate exchange grids between two grids for the second order
818 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
819 and lon_in,lat_in, lon_out,lat_out are geographic grid location of grid cell bounds.
820 mask is on grid lon_in/lat_in.
821********************************************************************************/
822int create_xgrid_2dx2d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
823 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
824 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
825 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
826{
827 int nxgrid;
828 nxgrid = create_xgrid_2dx2d_order2(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in,
829 j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat);
830 return nxgrid;
831
832}
833int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
834 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
835 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
836 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
837{
838
839 int nx1, nx2, ny1, ny2, nx1p, nx2p, nxgrid;
840 double *area_in, *area_out;
841 int nblocks =1;
842 int *istart2=NULL, *iend2=NULL;
843 int npts_left, nblks_left, pos, m, npts_my, ij;
844 double *lon_out_min_list,*lon_out_max_list,*lon_out_avg,*lat_out_min_list,*lat_out_max_list;
845 double *lon_out_list, *lat_out_list;
846 int *pnxgrid=NULL, *pstart;
847 int *pi_in=NULL, *pj_in=NULL, *pi_out=NULL, *pj_out=NULL;
848 double *pxgrid_area=NULL, *pxgrid_clon=NULL, *pxgrid_clat=NULL;
849 int *n2_list;
850 int nthreads, nxgrid_block_max;
851
852 nx1 = *nlon_in;
853 ny1 = *nlat_in;
854 nx2 = *nlon_out;
855 ny2 = *nlat_out;
856 nx1p = nx1 + 1;
857 nx2p = nx2 + 1;
858
859 area_in = (double *)malloc(nx1*ny1*sizeof(double));
860 area_out = (double *)malloc(nx2*ny2*sizeof(double));
861 get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in);
862 get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out);
863
864 nthreads = 1;
865#if defined(_OPENMP)
866#pragma omp parallel
867 nthreads = omp_get_num_threads();
868#endif
869
870 nblocks = nthreads;
871
872 istart2 = (int *)malloc(nblocks*sizeof(int));
873 iend2 = (int *)malloc(nblocks*sizeof(int));
874
875 pstart = (int *)malloc(nblocks*sizeof(int));
876 pnxgrid = (int *)malloc(nblocks*sizeof(int));
877
878 nxgrid_block_max = MAXXGRID/nblocks;
879
880 for(m=0; m<nblocks; m++) {
881 pnxgrid[m] = 0;
882 pstart[m] = m*nxgrid_block_max;
883 }
884
885 if(nblocks == 1) {
886 pi_in = i_in;
887 pj_in = j_in;
888 pi_out = i_out;
889 pj_out = j_out;
890 pxgrid_area = xgrid_area;
891 pxgrid_clon = xgrid_clon;
892 pxgrid_clat = xgrid_clat;
893 }
894 else {
895 pi_in = (int *)malloc(MAXXGRID*sizeof(int));
896 pj_in = (int *)malloc(MAXXGRID*sizeof(int));
897 pi_out = (int *)malloc(MAXXGRID*sizeof(int));
898 pj_out = (int *)malloc(MAXXGRID*sizeof(int));
899 pxgrid_area = (double *)malloc(MAXXGRID*sizeof(double));
900 pxgrid_clon = (double *)malloc(MAXXGRID*sizeof(double));
901 pxgrid_clat = (double *)malloc(MAXXGRID*sizeof(double));
902 }
903
904 npts_left = nx2*ny2;
905 nblks_left = nblocks;
906 pos = 0;
907 for(m=0; m<nblocks; m++) {
908 istart2[m] = pos;
909 npts_my = npts_left/nblks_left;
910 iend2[m] = istart2[m] + npts_my - 1;
911 pos = iend2[m] + 1;
912 npts_left -= npts_my;
913 nblks_left--;
914 }
915
916 lon_out_min_list = (double *)malloc(nx2*ny2*sizeof(double));
917 lon_out_max_list = (double *)malloc(nx2*ny2*sizeof(double));
918 lat_out_min_list = (double *)malloc(nx2*ny2*sizeof(double));
919 lat_out_max_list = (double *)malloc(nx2*ny2*sizeof(double));
920 lon_out_avg = (double *)malloc(nx2*ny2*sizeof(double));
921 n2_list = (int *)malloc(nx2*ny2*sizeof(int));
922 lon_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double));
923 lat_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double));
924#if defined(_OPENMP)
925#pragma omp parallel for default(none) shared(nx2,ny2,nx2p,lon_out,lat_out,lat_out_min_list, \
926 lat_out_max_list,lon_out_min_list,lon_out_max_list, \
927 lon_out_avg,n2_list,lon_out_list,lat_out_list)
928#endif
929 for(ij=0; ij<nx2*ny2; ij++){
930 int i2, j2, n, n0, n1, n2, n3, n2_in, l;
931 double x2_in[MV], y2_in[MV];
932 i2 = ij%nx2;
933 j2 = ij/nx2;
934 n = j2*nx2+i2;
935 n0 = j2*nx2p+i2; n1 = j2*nx2p+i2+1;
936 n2 = (j2+1)*nx2p+i2+1; n3 = (j2+1)*nx2p+i2;
937 x2_in[0] = lon_out[n0]; y2_in[0] = lat_out[n0];
938 x2_in[1] = lon_out[n1]; y2_in[1] = lat_out[n1];
939 x2_in[2] = lon_out[n2]; y2_in[2] = lat_out[n2];
940 x2_in[3] = lon_out[n3]; y2_in[3] = lat_out[n3];
941
942 lat_out_min_list[n] = minval_double(4, y2_in);
943 lat_out_max_list[n] = maxval_double(4, y2_in);
944 n2_in = fix_lon(x2_in, y2_in, 4, M_PI);
945 if(n2_in > MAX_V) error_handler("create_xgrid.c: n2_in is greater than MAX_V");
946 lon_out_min_list[n] = minval_double(n2_in, x2_in);
947 lon_out_max_list[n] = maxval_double(n2_in, x2_in);
948 lon_out_avg[n] = avgval_double(n2_in, x2_in);
949 n2_list[n] = n2_in;
950 for(l=0; l<n2_in; l++) {
951 lon_out_list[n*MAX_V+l] = x2_in[l];
952 lat_out_list[n*MAX_V+l] = y2_in[l];
953 }
954 }
955
956 nxgrid = 0;
957
958#if defined(_OPENMP)
959#pragma omp parallel for default(none) shared(nblocks,nx1,ny1,nx1p,mask_in,lon_in,lat_in, \
960 istart2,iend2,nx2,lat_out_min_list,lat_out_max_list, \
961 n2_list,lon_out_list,lat_out_list,lon_out_min_list, \
962 lon_out_max_list,lon_out_avg,area_in,area_out, \
963 pxgrid_area,pnxgrid,pxgrid_clon,pxgrid_clat,pi_in, \
964 pj_in,pi_out,pj_out,pstart,nthreads)
965#endif
966 for(m=0; m<nblocks; m++) {
967 int i1, j1, ij;
968 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
969 int n0, n1, n2, n3, l,n1_in;
970 double lat_in_min,lat_in_max,lon_in_min,lon_in_max,lon_in_avg;
971 double x1_in[MV], y1_in[MV], x_out[MV], y_out[MV];
972
973 n0 = j1*nx1p+i1; n1 = j1*nx1p+i1+1;
974 n2 = (j1+1)*nx1p+i1+1; n3 = (j1+1)*nx1p+i1;
975 x1_in[0] = lon_in[n0]; y1_in[0] = lat_in[n0];
976 x1_in[1] = lon_in[n1]; y1_in[1] = lat_in[n1];
977 x1_in[2] = lon_in[n2]; y1_in[2] = lat_in[n2];
978 x1_in[3] = lon_in[n3]; y1_in[3] = lat_in[n3];
979 lat_in_min = minval_double(4, y1_in);
980 lat_in_max = maxval_double(4, y1_in);
981 n1_in = fix_lon(x1_in, y1_in, 4, M_PI);
982 lon_in_min = minval_double(n1_in, x1_in);
983 lon_in_max = maxval_double(n1_in, x1_in);
984 lon_in_avg = avgval_double(n1_in, x1_in);
985 for(ij=istart2[m]; ij<=iend2[m]; ij++) {
986 int n_out, i2, j2, n2_in;
987 double xarea, dx, lon_out_min, lon_out_max;
988 double x2_in[MAX_V], y2_in[MAX_V];
989
990 i2 = ij%nx2;
991 j2 = ij/nx2;
992
993 if(lat_out_min_list[ij] >= lat_in_max || lat_out_max_list[ij] <= lat_in_min ) continue;
994 /* adjust x2_in according to lon_in_avg*/
995 n2_in = n2_list[ij];
996 for(l=0; l<n2_in; l++) {
997 x2_in[l] = lon_out_list[ij*MAX_V+l];
998 y2_in[l] = lat_out_list[ij*MAX_V+l];
999 }
1000 lon_out_min = lon_out_min_list[ij];
1001 lon_out_max = lon_out_max_list[ij];
1002 dx = lon_out_avg[ij] - lon_in_avg;
1003 if(dx < -M_PI ) {
1004 lon_out_min += TPI;
1005 lon_out_max += TPI;
1006 for (l=0; l<n2_in; l++) x2_in[l] += TPI;
1007 }
1008 else if (dx > M_PI) {
1009 lon_out_min -= TPI;
1010 lon_out_max -= TPI;
1011 for (l=0; l<n2_in; l++) x2_in[l] -= TPI;
1012 }
1013
1014 /* x2_in should in the same range as x1_in after lon_fix, so no need to
1015 consider cyclic condition
1016 */
1017 if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue;
1018 if ( (n_out = clip_2dx2d( x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out )) > 0) {
1019 double min_area;
1020 int nn;
1021 xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
1022 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
1023 if( xarea/min_area > AREA_RATIO_THRESH ) {
1024 pnxgrid[m]++;
1025 if(pnxgrid[m]>= MAXXGRID/nthreads)
1026 error_handler("nxgrid is greater than MAXXGRID/nthreads, increase MAXXGRID, decrease nthreads, or increase number of MPI ranks");
1027 nn = pstart[m] + pnxgrid[m]-1;
1028 pxgrid_area[nn] = xarea;
1029 pxgrid_clon[nn] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg);
1030 pxgrid_clat[nn] = poly_ctrlat (x_out, y_out, n_out );
1031 pi_in[nn] = i1;
1032 pj_in[nn] = j1;
1033 pi_out[nn] = i2;
1034 pj_out[nn] = j2;
1035 }
1036 }
1037 }
1038 }
1039 }
1040
1041 /*copy data if nblocks > 1 */
1042 if(nblocks == 1) {
1043 nxgrid = pnxgrid[0];
1044 pi_in = NULL;
1045 pj_in = NULL;
1046 pi_out = NULL;
1047 pj_out = NULL;
1048 pxgrid_area = NULL;
1049 pxgrid_clon = NULL;
1050 pxgrid_clat = NULL;
1051 }
1052 else {
1053 int nn, i;
1054 nxgrid = 0;
1055 for(m=0; m<nblocks; m++) {
1056 for(i=0; i<pnxgrid[m]; i++) {
1057 nn = pstart[m] + i;
1058 i_in[nxgrid] = pi_in[nn];
1059 j_in[nxgrid] = pj_in[nn];
1060 i_out[nxgrid] = pi_out[nn];
1061 j_out[nxgrid] = pj_out[nn];
1062 xgrid_area[nxgrid] = pxgrid_area[nn];
1063 xgrid_clon[nxgrid] = pxgrid_clon[nn];
1064 xgrid_clat[nxgrid] = pxgrid_clat[nn];
1065 nxgrid++;
1066 }
1067 }
1068 free(pi_in);
1069 free(pj_in);
1070 free(pi_out);
1071 free(pj_out);
1072 free(pxgrid_area);
1073 free(pxgrid_clon);
1074 free(pxgrid_clat);
1075 }
1076
1077 free(area_in);
1078 free(area_out);
1079 free(lon_out_min_list);
1080 free(lon_out_max_list);
1081 free(lat_out_min_list);
1082 free(lat_out_max_list);
1083 free(lon_out_avg);
1084 free(n2_list);
1085 free(lon_out_list);
1086 free(lat_out_list);
1087
1088 return nxgrid;
1089
1090}/* get_xgrid_2Dx2D_order2 */
1091
1092int create_xgrid_great_circle_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
1093 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1094 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
1095 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1096{
1097 int nxgrid;
1098 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out,
1099 mask_in, i_in, j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat);
1100
1101 return nxgrid;
1102}
1103
1104int create_xgrid_great_circle(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
1105 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1106 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
1107 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1108{
1109
1110 int nx1, nx2, ny1, ny2, nx1p, nx2p, ny1p, ny2p, nxgrid, n1_in, n2_in;
1111 int n0, n1, n2, n3, i1, j1, i2, j2;
1112 double x1_in[MV], y1_in[MV], z1_in[MV];
1113 double x2_in[MV], y2_in[MV], z2_in[MV];
1114 double x_out[MV], y_out[MV], z_out[MV];
1115 double *x1=NULL, *y1=NULL, *z1=NULL;
1116 double *x2=NULL, *y2=NULL, *z2=NULL;
1117
1118 double *area1, *area2, min_area;
1119
1120 nx1 = *nlon_in;
1121 ny1 = *nlat_in;
1122 nx2 = *nlon_out;
1123 ny2 = *nlat_out;
1124 nxgrid = 0;
1125 nx1p = nx1 + 1;
1126 nx2p = nx2 + 1;
1127 ny1p = ny1 + 1;
1128 ny2p = ny2 + 1;
1129
1130 /* first convert lon-lat to cartesian coordinates */
1131 x1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1132 y1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1133 z1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1134 x2 = (double *)malloc(nx2p*ny2p*sizeof(double));
1135 y2 = (double *)malloc(nx2p*ny2p*sizeof(double));
1136 z2 = (double *)malloc(nx2p*ny2p*sizeof(double));
1137
1138 latlon2xyz(nx1p*ny1p, lon_in, lat_in, x1, y1, z1);
1139 latlon2xyz(nx2p*ny2p, lon_out, lat_out, x2, y2, z2);
1140
1141 area1 = (double *)malloc(nx1*ny1*sizeof(double));
1142 area2 = (double *)malloc(nx2*ny2*sizeof(double));
1143 get_grid_great_circle_area(nlon_in, nlat_in, lon_in, lat_in, area1);
1144 get_grid_great_circle_area(nlon_out, nlat_out, lon_out, lat_out, area2);
1145 n1_in = 4;
1146 n2_in = 4;
1147
1148 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
1149 /* clockwise */
1150 n0 = j1*nx1p+i1; n1 = (j1+1)*nx1p+i1;
1151 n2 = (j1+1)*nx1p+i1+1; n3 = j1*nx1p+i1+1;
1152 x1_in[0] = x1[n0]; y1_in[0] = y1[n0]; z1_in[0] = z1[n0];
1153 x1_in[1] = x1[n1]; y1_in[1] = y1[n1]; z1_in[1] = z1[n1];
1154 x1_in[2] = x1[n2]; y1_in[2] = y1[n2]; z1_in[2] = z1[n2];
1155 x1_in[3] = x1[n3]; y1_in[3] = y1[n3]; z1_in[3] = z1[n3];
1156
1157 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
1158 int n_out;
1159 double xarea;
1160
1161 n0 = j2*nx2p+i2; n1 = (j2+1)*nx2p+i2;
1162 n2 = (j2+1)*nx2p+i2+1; n3 = j2*nx2p+i2+1;
1163 x2_in[0] = x2[n0]; y2_in[0] = y2[n0]; z2_in[0] = z2[n0];
1164 x2_in[1] = x2[n1]; y2_in[1] = y2[n1]; z2_in[1] = z2[n1];
1165 x2_in[2] = x2[n2]; y2_in[2] = y2[n2]; z2_in[2] = z2[n2];
1166 x2_in[3] = x2[n3]; y2_in[3] = y2[n3]; z2_in[3] = z2[n3];
1167
1168 if ( (n_out = clip_2dx2d_great_circle( x1_in, y1_in, z1_in, n1_in, x2_in, y2_in, z2_in, n2_in,
1169 x_out, y_out, z_out)) > 0) {
1170 xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1];
1171 min_area = min(area1[j1*nx1+i1], area2[j2*nx2+i2]);
1172 if( xarea/min_area > AREA_RATIO_THRESH ) {
1173 xgrid_area[nxgrid] = xarea;
1174 xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */
1175 xgrid_clat[nxgrid] = 0;
1176 i_in[nxgrid] = i1;
1177 j_in[nxgrid] = j1;
1178 i_out[nxgrid] = i2;
1179 j_out[nxgrid] = j2;
1180 ++nxgrid;
1181 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
1182 }
1183 }
1184 }
1185 }
1186
1187
1188 free(area1);
1189 free(area2);
1190
1191 free(x1);
1192 free(y1);
1193 free(z1);
1194 free(x2);
1195 free(y2);
1196 free(z2);
1197
1198 return nxgrid;
1199
1200}/* create_xgrid_great_circle */
1201
1202int create_xgrid_great_circle_ug_(const int *nlon_in, const int *nlat_in, const int *npts_out,
1203 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1204 const double *mask_in, int *i_in, int *j_in, int *l_out,
1205 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1206{
1207 int nxgrid;
1208 nxgrid = create_xgrid_great_circle_ug(nlon_in, nlat_in, npts_out, lon_in, lat_in, lon_out, lat_out,
1209 mask_in, i_in, j_in, l_out, xgrid_area, xgrid_clon, xgrid_clat);
1210
1211 return nxgrid;
1212}
1213
1214int create_xgrid_great_circle_ug(const int *nlon_in, const int *nlat_in, const int *npts_out,
1215 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1216 const double *mask_in, int *i_in, int *j_in, int *l_out,
1217 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1218{
1219
1220 int nx1, ny1, npts2, nx1p, ny1p, nxgrid, n1_in, n2_in, nv;
1221 int n0, n1, n2, n3, i1, j1, l2;
1222 double x1_in[MV], y1_in[MV], z1_in[MV];
1223 double x2_in[MV], y2_in[MV], z2_in[MV];
1224 double x_out[MV], y_out[MV], z_out[MV];
1225 double *x1=NULL, *y1=NULL, *z1=NULL;
1226 double *x2=NULL, *y2=NULL, *z2=NULL;
1227
1228 double *area1, *area2, min_area;
1229
1230 nx1 = *nlon_in;
1231 ny1 = *nlat_in;
1232 nv = 4;
1233 npts2 = *npts_out;
1234 nxgrid = 0;
1235 nx1p = nx1 + 1;
1236 ny1p = ny1 + 1;
1237
1238 /* first convert lon-lat to cartesian coordinates */
1239 x1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1240 y1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1241 z1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1242 x2 = (double *)malloc(npts2*nv*sizeof(double));
1243 y2 = (double *)malloc(npts2*nv*sizeof(double));
1244 z2 = (double *)malloc(npts2*nv*sizeof(double));
1245
1246 latlon2xyz(nx1p*ny1p, lon_in, lat_in, x1, y1, z1);
1247 latlon2xyz(npts2*nv, lon_out, lat_out, x2, y2, z2);
1248
1249 area1 = (double *)malloc(nx1*ny1*sizeof(double));
1250 area2 = (double *)malloc(npts2*sizeof(double));
1251 get_grid_great_circle_area(nlon_in, nlat_in, lon_in, lat_in, area1);
1252 get_grid_great_circle_area_ug(npts_out, lon_out, lat_out, area2);
1253 n1_in = 4;
1254 n2_in = 4;
1255
1256 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
1257 /* clockwise */
1258 n0 = j1*nx1p+i1; n1 = (j1+1)*nx1p+i1;
1259 n2 = (j1+1)*nx1p+i1+1; n3 = j1*nx1p+i1+1;
1260 x1_in[0] = x1[n0]; y1_in[0] = y1[n0]; z1_in[0] = z1[n0];
1261 x1_in[1] = x1[n1]; y1_in[1] = y1[n1]; z1_in[1] = z1[n1];
1262 x1_in[2] = x1[n2]; y1_in[2] = y1[n2]; z1_in[2] = z1[n2];
1263 x1_in[3] = x1[n3]; y1_in[3] = y1[n3]; z1_in[3] = z1[n3];
1264
1265 for(l2=0; l2<npts2; l2++) {
1266 int n_out;
1267 double xarea;
1268
1269 n0 = l2*nv; n1 = l2*nv+1;
1270 n2 = l2*nv+2; n3 = l2*nv+3;
1271 x2_in[0] = x2[n0]; y2_in[0] = y2[n0]; z2_in[0] = z2[n0];
1272 x2_in[1] = x2[n1]; y2_in[1] = y2[n1]; z2_in[1] = z2[n1];
1273 x2_in[2] = x2[n2]; y2_in[2] = y2[n2]; z2_in[2] = z2[n2];
1274 x2_in[3] = x2[n3]; y2_in[3] = y2[n3]; z2_in[3] = z2[n3];
1275
1276 if ( (n_out = clip_2dx2d_great_circle( x1_in, y1_in, z1_in, n1_in, x2_in, y2_in, z2_in, n2_in,
1277 x_out, y_out, z_out)) > 0) {
1278 xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1];
1279 min_area = min(area1[j1*nx1+i1], area2[l2]);
1280 if( xarea/min_area > AREA_RATIO_THRESH ) {
1281 xgrid_area[nxgrid] = xarea;
1282 xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */
1283 xgrid_clat[nxgrid] = 0;
1284 i_in[nxgrid] = i1;
1285 j_in[nxgrid] = j1;
1286 l_out[nxgrid] = l2;
1287 ++nxgrid;
1288 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
1289 }
1290 }
1291 }
1292 }
1293
1294
1295 free(area1);
1296 free(area2);
1297
1298 free(x1);
1299 free(y1);
1300 free(z1);
1301 free(x2);
1302 free(y2);
1303 free(z2);
1304
1305 return nxgrid;
1306
1307}/* create_xgrid_great_circle_ug */
1308
1309/*******************************************************************************
1310 int get_maxxgrid
1311 return constants MAXXGRID.
1312*******************************************************************************/
1313int get_maxxgrid(void)
1314{
1315 return MAXXGRID;
1316}
1317
1318int get_maxxgrid_(void)
1319{
1320 return get_maxxgrid();
1321}
diag_grid_mod::latlon2xyz
pure elemental type(point) function latlon2xyz(lat, lon)
Return the (x,y,z) position of a given (lat,lon) point.
Definition diag_grid.F90:1019