?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fix: out-of-bounds memory access in axis_utils2 (#1157)
[FMS.git] / mosaic / create_xgrid.c
blob7698303b92d9b5da79273441e9d326ed0f73c537
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 "mosaic_util.h"
23 #include "create_xgrid.h"
24 #include "constant.h"
25 #if defined(_OPENMP)
26 #include <omp.h>
27 #endif
28
29 #define AREA_RATIO_THRESH (1.e-6)
30 #define MASK_THRESH (0.5)
31 #define EPSLN8 (1.e-8)
32 #define EPSLN30 (1.0e-30)
33 #define EPSLN10 (1.0e-10)
34 #define R2D (180/M_PI)
35 #define TPI (2.0*M_PI)
36
37 /** \file
38 * \ingroup mosaic
39 * \brief Grid creation and calculation functions for use in @ref mosaic_mod
40 * /
41
42 /*******************************************************************************
43 int get_maxxgrid
44 return constants MAXXGRID.
45 *******************************************************************************/
46 int get_maxxgrid(void)
47 {
48 return MAXXGRID;
49 }
50
51 int get_maxxgrid_(void)
52 {
53 return get_maxxgrid();
54 }
55
56
57 /*******************************************************************************
58 void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, const double *area)
59 return the grid area.
60 *******************************************************************************/
61 void get_grid_area_(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area)
62 {
63 get_grid_area(nlon, nlat, lon, lat, area);
64 }
65
66 void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area)
67 {
68 int nx, ny, nxp, i, j, n_in;
69 double x_in[20], y_in[20];
70
71 nx = *nlon;
72 ny = *nlat;
73 nxp = nx + 1;
74
75 for(j=0; j<ny; j++) for(i=0; i < nx; i++) {
76 x_in[0] = lon[j*nxp+i];
77 x_in[1] = lon[j*nxp+i+1];
78 x_in[2] = lon[(j+1)*nxp+i+1];
79 x_in[3] = lon[(j+1)*nxp+i];
80 y_in[0] = lat[j*nxp+i];
81 y_in[1] = lat[j*nxp+i+1];
82 y_in[2] = lat[(j+1)*nxp+i+1];
83 y_in[3] = lat[(j+1)*nxp+i];
84 n_in = fix_lon(x_in, y_in, 4, M_PI);
85 area[j*nx+i] = poly_area(x_in, y_in, n_in);
86 }
87
88 } /* get_grid_area */
89
90
91 /*******************************************************************************
92 void get_grid_area_ug(const int *npts, const double *lon, const double *lat, const double *area)
93 return the grid area.
94 *******************************************************************************/
95 void get_grid_area_ug_(const int *npts, const double *lon, const double *lat, double *area)
96 {
97 get_grid_area_ug(npts, lon, lat, area);
98 }
99
100 void get_grid_area_ug(const int *npts, const double *lon, const double *lat, double *area)
101 {
102 int nl, l, n_in, nv;
103 double x_in[20], y_in[20];
104
105 nl = *npts;
106 nv = 4;
107
108 for(l=0; l<nl; l++) {
109 x_in[0] = lon[l*nv];
110 x_in[1] = lon[l*nv+1];
111 x_in[2] = lon[l*nv+2];
112 x_in[3] = lon[l*nv+3];
113 y_in[0] = lat[l*nv];
114 y_in[1] = lat[l*nv+1];
115 y_in[2] = lat[l*nv+2];
116 y_in[3] = lat[l*nv+3];
117 n_in = fix_lon(x_in, y_in, nv, M_PI);
118 area[l] = poly_area(x_in, y_in, n_in);
119 }
120
121 } /* get_grid_area_ug */
122
123
124 void get_grid_great_circle_area_(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area)
125 {
126 get_grid_great_circle_area(nlon, nlat, lon, lat, area);
127
128 }
129
130 void get_grid_great_circle_area(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area)
131 {
132 int nx, ny, nxp, nyp, i, j;
133 int n0, n1, n2, n3;
134 struct Node *grid=NULL;
135 double *x=NULL, *y=NULL, *z=NULL;
136
137
138 nx = *nlon;
139 ny = *nlat;
140 nxp = nx + 1;
141 nyp = ny + 1;
142
143 x = (double *)malloc(nxp*nyp*sizeof(double));
144 y = (double *)malloc(nxp*nyp*sizeof(double));
145 z = (double *)malloc(nxp*nyp*sizeof(double));
146
147 latlon2xyz(nxp*nyp, lon, lat, x, y, z);
148
149 for(j=0; j<ny; j++) for(i=0; i < nx; i++) {
150 /* clockwise */
151 n0 = j*nxp+i;
152 n1 = (j+1)*nxp+i;
153 n2 = (j+1)*nxp+i+1;
154 n3 = j*nxp+i+1;
155 rewindList();
156 grid = getNext();
157 addEnd(grid, x[n0], y[n0], z[n0], 0, 0, 0, -1);
158 addEnd(grid, x[n1], y[n1], z[n1], 0, 0, 0, -1);
159 addEnd(grid, x[n2], y[n2], z[n2], 0, 0, 0, -1);
160 addEnd(grid, x[n3], y[n3], z[n3], 0, 0, 0, -1);
161 area[j*nx+i] = gridArea(grid);
162 }
163
164 free(x);
165 free(y);
166 free(z);
167
168 } /* get_grid_great_circle_area */
169
170 void get_grid_great_circle_area_ug_(const int *npts, const double *lon, const double *lat, double *area)
171 {
172 get_grid_great_circle_area_ug(npts, lon, lat, area);
173
174 }
175
176 void get_grid_great_circle_area_ug(const int *npts, const double *lon, const double *lat, double *area)
177 {
178 int l, nl, nv;
179 int n0, n1, n2, n3;
180 struct Node *grid=NULL;
181 double *x=NULL, *y=NULL, *z=NULL;
182
183 nl = *npts;
184 nv = 4;
185
186 x = (double *)malloc(nl*nv*sizeof(double));
187 y = (double *)malloc(nl*nv*sizeof(double));
188 z = (double *)malloc(nl*nv*sizeof(double));
189
190 latlon2xyz(nl*nv, lon, lat, x, y, z);
191
192 for(l=0; l<nv; l++) {
193 /* clockwise */
194 n0 = l*nv;
195 n1 = l*nv+1;
196 n2 = l*nv+2;
197 n3 = l*nv+3;
198 rewindList();
199 grid = getNext();
200 addEnd(grid, x[n0], y[n0], z[n0], 0, 0, 0, -1);
201 addEnd(grid, x[n1], y[n1], z[n1], 0, 0, 0, -1);
202 addEnd(grid, x[n2], y[n2], z[n2], 0, 0, 0, -1);
203 addEnd(grid, x[n3], y[n3], z[n3], 0, 0, 0, -1);
204 area[l] = gridArea(grid);
205 }
206
207 free(x);
208 free(y);
209 free(z);
210
211 } /* get_grid_great_circle_area_ug */
212
213 void get_grid_area_dimensionless(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area)
214 {
215 int nx, ny, nxp, i, j, n_in;
216 double x_in[20], y_in[20];
217
218 nx = *nlon;
219 ny = *nlat;
220 nxp = nx + 1;
221
222 for(j=0; j<ny; j++) for(i=0; i < nx; i++) {
223 x_in[0] = lon[j*nxp+i];
224 x_in[1] = lon[j*nxp+i+1];
225 x_in[2] = lon[(j+1)*nxp+i+1];
226 x_in[3] = lon[(j+1)*nxp+i];
227 y_in[0] = lat[j*nxp+i];
228 y_in[1] = lat[j*nxp+i+1];
229 y_in[2] = lat[(j+1)*nxp+i+1];
230 y_in[3] = lat[(j+1)*nxp+i];
231 n_in = fix_lon(x_in, y_in, 4, M_PI);
232 area[j*nx+i] = poly_area_dimensionless(x_in, y_in, n_in);
233 }
234
235 } /* get_grid_area */
236
237
238
239 void get_grid_area_no_adjust(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area)
240 {
241 int nx, ny, nxp, i, j, n_in;
242 double x_in[20], y_in[20];
243
244 nx = *nlon;
245 ny = *nlat;
246 nxp = nx + 1;
247
248 for(j=0; j<ny; j++) for(i=0; i < nx; i++) {
249 x_in[0] = lon[j*nxp+i];
250 x_in[1] = lon[j*nxp+i+1];
251 x_in[2] = lon[(j+1)*nxp+i+1];
252 x_in[3] = lon[(j+1)*nxp+i];
253 y_in[0] = lat[j*nxp+i];
254 y_in[1] = lat[j*nxp+i+1];
255 y_in[2] = lat[(j+1)*nxp+i+1];
256 y_in[3] = lat[(j+1)*nxp+i];
257 n_in = 4;
258 area[j*nx+i] = poly_area_no_adjust(x_in, y_in, n_in);
259 }
260
261 } /* get_grid_area_no_adjust */
262
263 /*******************************************************************************
264 void create_xgrid_1dx2d_order1
265 This routine generate exchange grids between two grids for the first order
266 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
267 and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds.
268 *******************************************************************************/
269 int create_xgrid_1dx2d_order1_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
270 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
271 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area)
272 {
273 int nxgrid;
274
275 nxgrid = create_xgrid_1dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in,
276 i_in, j_in, i_out, j_out, xgrid_area);
277 return nxgrid;
278
279 }
280
281 int create_xgrid_1dx2d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in,
282 const double *lat_in, const double *lon_out, const double *lat_out,
283 const double *mask_in, int *i_in, int *j_in, int *i_out,
284 int *j_out, double *xgrid_area)
285 {
286
287 int nx1, ny1, nx2, ny2, nx1p, nx2p;
288 int i1, j1, i2, j2, nxgrid;
289 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
290 double *area_in, *area_out, min_area;
291 double *tmpx, *tmpy;
292
293 nx1 = *nlon_in;
294 ny1 = *nlat_in;
295 nx2 = *nlon_out;
296 ny2 = *nlat_out;
297
298 nxgrid = 0;
299 nx1p = nx1 + 1;
300 nx2p = nx2 + 1;
301
302 area_in = (double *)malloc(nx1*ny1*sizeof(double));
303 area_out = (double *)malloc(nx2*ny2*sizeof(double));
304 tmpx = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
305 tmpy = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
306 for(j1=0; j1<=ny1; j1++) for(i1=0; i1<=nx1; i1++) {
307 tmpx[j1*nx1p+i1] = lon_in[i1];
308 tmpy[j1*nx1p+i1] = lat_in[j1];
309 }
310 /* This is just a temporary fix to solve the issue that there is one point in zonal direction */
311 if(nx1 > 1)
312 get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in);
313 else
314 get_grid_area_no_adjust(nlon_in, nlat_in, tmpx, tmpy, area_in);
315
316 get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out);
317 free(tmpx);
318 free(tmpy);
319
320 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
321
322 ll_lon = lon_in[i1]; ll_lat = lat_in[j1];
323 ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1];
324 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
325 int n_in, n_out;
326 double Xarea;
327
328 y_in[0] = lat_out[j2*nx2p+i2];
329 y_in[1] = lat_out[j2*nx2p+i2+1];
330 y_in[2] = lat_out[(j2+1)*nx2p+i2+1];
331 y_in[3] = lat_out[(j2+1)*nx2p+i2];
332 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
333 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
334 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
335 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
336
337 x_in[0] = lon_out[j2*nx2p+i2];
338 x_in[1] = lon_out[j2*nx2p+i2+1];
339 x_in[2] = lon_out[(j2+1)*nx2p+i2+1];
340 x_in[3] = lon_out[(j2+1)*nx2p+i2];
341 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
342
343 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
344 Xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
345 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
346 if( Xarea/min_area > AREA_RATIO_THRESH ) {
347 xgrid_area[nxgrid] = Xarea;
348 i_in[nxgrid] = i1;
349 j_in[nxgrid] = j1;
350 i_out[nxgrid] = i2;
351 j_out[nxgrid] = j2;
352 ++nxgrid;
353 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
354 }
355 }
356 }
357 }
358
359 free(area_in);
360 free(area_out);
361
362 return nxgrid;
363
364 } /* create_xgrid_1dx2d_order1 */
365
366
367 /*******************************************************************************
368 void create_xgrid_1dx2d_order1_ug
369 This routine generate exchange grids between two grids for the first order
370 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
371 and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds.
372 *******************************************************************************/
373 int create_xgrid_1dx2d_order1_ug_(const int *nlon_in, const int *nlat_in, const int *npts_out,
374 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
375 const double *mask_in, int *i_in, int *j_in, int *l_out, double *xgrid_area)
376 {
377 int nxgrid;
378
379 nxgrid = create_xgrid_1dx2d_order1_ug(nlon_in, nlat_in, npts_out, lon_in, lat_in, lon_out, lat_out, mask_in,
380 i_in, j_in, l_out, xgrid_area);
381 return nxgrid;
382
383 }
384
385 int create_xgrid_1dx2d_order1_ug(const int *nlon_in, const int *nlat_in, const int *npts_out, const double *lon_in,
386 const double *lat_in, const double *lon_out, const double *lat_out,
387 const double *mask_in, int *i_in, int *j_in, int *l_out, double *xgrid_area)
388 {
389
390 int nx1, ny1, nx1p, nv, npts2;
391 int i1, j1, l2, nxgrid;
392 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
393 double *area_in, *area_out, min_area;
394 double *tmpx, *tmpy;
395
396 nx1 = *nlon_in;
397 ny1 = *nlat_in;
398 nv = 4;
399 npts2 = *npts_out;
400
401 nxgrid = 0;
402 nx1p = nx1 + 1;
403
404 area_in = (double *)malloc(nx1*ny1*sizeof(double));
405 area_out = (double *)malloc(npts2*sizeof(double));
406 tmpx = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
407 tmpy = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
408 for(j1=0; j1<=ny1; j1++) for(i1=0; i1<=nx1; i1++) {
409 tmpx[j1*nx1p+i1] = lon_in[i1];
410 tmpy[j1*nx1p+i1] = lat_in[j1];
411 }
412 /* This is just a temporary fix to solve the issue that there is one point in zonal direction */
413 if(nx1 > 1)
414 get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in);
415 else
416 get_grid_area_no_adjust(nlon_in, nlat_in, tmpx, tmpy, area_in);
417
418 get_grid_area_ug(npts_out, lon_out, lat_out, area_out);
419 free(tmpx);
420 free(tmpy);
421
422 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
423
424 ll_lon = lon_in[i1]; ll_lat = lat_in[j1];
425 ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1];
426 for(l2=0; l2<npts2; l2++) {
427 int n_in, n_out;
428 double Xarea;
429
430 y_in[0] = lat_out[l2*nv];
431 y_in[1] = lat_out[l2*nv+1];
432 y_in[2] = lat_out[l2*nv+2];
433 y_in[3] = lat_out[l2*nv+3];
434 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
435 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
436 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
437 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
438
439 x_in[0] = lon_out[l2*nv];
440 x_in[1] = lon_out[l2*nv+1];
441 x_in[2] = lon_out[l2*nv+2];
442 x_in[3] = lon_out[l2*nv+3];
443 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
444
445 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
446 Xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
447 min_area = min(area_in[j1*nx1+i1], area_out[l2]);
448 if( Xarea/min_area > AREA_RATIO_THRESH ) {
449 xgrid_area[nxgrid] = Xarea;
450 i_in[nxgrid] = i1;
451 j_in[nxgrid] = j1;
452 l_out[nxgrid] = l2;
453 ++nxgrid;
454 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
455 }
456 }
457 }
458 }
459
460 free(area_in);
461 free(area_out);
462
463 return nxgrid;
464
465 } /* create_xgrid_1dx2d_order1_ug */
466
467 /********************************************************************************
468 void create_xgrid_1dx2d_order2
469 This routine generate exchange grids between two grids for the second order
470 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
471 and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds.
472 ********************************************************************************/
473 int create_xgrid_1dx2d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
474 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
475 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
476 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
477 {
478 int nxgrid;
479 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,
480 j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat);
481 return nxgrid;
482
483 }
484 int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
485 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
486 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
487 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
488 {
489
490 int nx1, ny1, nx2, ny2, nx1p, nx2p;
491 int i1, j1, i2, j2, nxgrid;
492 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
493 double *area_in, *area_out, min_area;
494 double *tmpx, *tmpy;
495
496 nx1 = *nlon_in;
497 ny1 = *nlat_in;
498 nx2 = *nlon_out;
499 ny2 = *nlat_out;
500
501 nxgrid = 0;
502 nx1p = nx1 + 1;
503 nx2p = nx2 + 1;
504
505 area_in = (double *)malloc(nx1*ny1*sizeof(double));
506 area_out = (double *)malloc(nx2*ny2*sizeof(double));
507 tmpx = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
508 tmpy = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double));
509 for(j1=0; j1<=ny1; j1++) for(i1=0; i1<=nx1; i1++) {
510 tmpx[j1*nx1p+i1] = lon_in[i1];
511 tmpy[j1*nx1p+i1] = lat_in[j1];
512 }
513 get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in);
514 get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out);
515 free(tmpx);
516 free(tmpy);
517
518 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
519
520 ll_lon = lon_in[i1]; ll_lat = lat_in[j1];
521 ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1];
522 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
523 int n_in, n_out;
524 double xarea, lon_in_avg;
525
526 y_in[0] = lat_out[j2*nx2p+i2];
527 y_in[1] = lat_out[j2*nx2p+i2+1];
528 y_in[2] = lat_out[(j2+1)*nx2p+i2+1];
529 y_in[3] = lat_out[(j2+1)*nx2p+i2];
530 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
531 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
532 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
533 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
534
535 x_in[0] = lon_out[j2*nx2p+i2];
536 x_in[1] = lon_out[j2*nx2p+i2+1];
537 x_in[2] = lon_out[(j2+1)*nx2p+i2+1];
538 x_in[3] = lon_out[(j2+1)*nx2p+i2];
539 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
540 lon_in_avg = avgval_double(n_in, x_in);
541
542 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
543 xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
544 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
545 if(xarea/min_area > AREA_RATIO_THRESH ) {
546 xgrid_area[nxgrid] = xarea;
547 xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg);
548 xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out );
549 i_in[nxgrid] = i1;
550 j_in[nxgrid] = j1;
551 i_out[nxgrid] = i2;
552 j_out[nxgrid] = j2;
553 ++nxgrid;
554 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
555 }
556 }
557 }
558 }
559 free(area_in);
560 free(area_out);
561
562 return nxgrid;
563
564 } /* create_xgrid_1dx2d_order2 */
565
566 /*******************************************************************************
567 void create_xgrid_2dx1d_order1
568 This routine generate exchange grids between two grids for the first order
569 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
570 and lon_out,lat_out are 1-D grid bounds, lon_in,lat_in are geographic grid location of grid cell bounds.
571 mask is on grid lon_in/lat_in.
572 *******************************************************************************/
573 int create_xgrid_2dx1d_order1_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
574 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
575 const double *mask_in, int *i_in, int *j_in, int *i_out,
576 int *j_out, double *xgrid_area)
577 {
578 int nxgrid;
579
580 nxgrid = create_xgrid_2dx1d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in,
581 i_in, j_in, i_out, j_out, xgrid_area);
582 return nxgrid;
583
584 }
585 int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in,
586 const double *lat_in, const double *lon_out, const double *lat_out,
587 const double *mask_in, int *i_in, int *j_in, int *i_out,
588 int *j_out, double *xgrid_area)
589 {
590
591 int nx1, ny1, nx2, ny2, nx1p, nx2p;
592 int i1, j1, i2, j2, nxgrid;
593 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
594 double *area_in, *area_out, min_area;
595 double *tmpx, *tmpy;
596 int n_in, n_out;
597 double Xarea;
598
599
600 nx1 = *nlon_in;
601 ny1 = *nlat_in;
602 nx2 = *nlon_out;
603 ny2 = *nlat_out;
604
605 nxgrid = 0;
606 nx1p = nx1 + 1;
607 nx2p = nx2 + 1;
608 area_in = (double *)malloc(nx1*ny1*sizeof(double));
609 area_out = (double *)malloc(nx2*ny2*sizeof(double));
610 tmpx = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double));
611 tmpy = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double));
612 for(j2=0; j2<=ny2; j2++) for(i2=0; i2<=nx2; i2++) {
613 tmpx[j2*nx2p+i2] = lon_out[i2];
614 tmpy[j2*nx2p+i2] = lat_out[j2];
615 }
616 get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in);
617 get_grid_area(nlon_out, nlat_out, tmpx, tmpy, area_out);
618
619 free(tmpx);
620 free(tmpy);
621
622 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
623
624 ll_lon = lon_out[i2]; ll_lat = lat_out[j2];
625 ur_lon = lon_out[i2+1]; ur_lat = lat_out[j2+1];
626 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
627
628 y_in[0] = lat_in[j1*nx1p+i1];
629 y_in[1] = lat_in[j1*nx1p+i1+1];
630 y_in[2] = lat_in[(j1+1)*nx1p+i1+1];
631 y_in[3] = lat_in[(j1+1)*nx1p+i1];
632 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
633 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
634 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
635 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
636
637 x_in[0] = lon_in[j1*nx1p+i1];
638 x_in[1] = lon_in[j1*nx1p+i1+1];
639 x_in[2] = lon_in[(j1+1)*nx1p+i1+1];
640 x_in[3] = lon_in[(j1+1)*nx1p+i1];
641
642 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
643
644 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
645 Xarea = poly_area ( x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
646 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
647 if( Xarea/min_area > AREA_RATIO_THRESH ) {
648 xgrid_area[nxgrid] = Xarea;
649 i_in[nxgrid] = i1;
650 j_in[nxgrid] = j1;
651 i_out[nxgrid] = i2;
652 j_out[nxgrid] = j2;
653 ++nxgrid;
654 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
655 }
656 }
657 }
658 }
659
660 free(area_in);
661 free(area_out);
662
663 return nxgrid;
664
665 } /* create_xgrid_2dx1d_order1 */
666
667
668 /********************************************************************************
669 void create_xgrid_2dx1d_order2
670 This routine generate exchange grids between two grids for the second order
671 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
672 and lon_out,lat_out are 1-D grid bounds, lon_in,lat_in are geographic grid location of grid cell bounds.
673 mask is on grid lon_in/lat_in.
674 ********************************************************************************/
675 int create_xgrid_2dx1d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
676 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
677 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
678 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
679 {
680 int nxgrid;
681 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,
682 j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat);
683 return nxgrid;
684
685 }
686
687 int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
688 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
689 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
690 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
691 {
692
693 int nx1, ny1, nx2, ny2, nx1p, nx2p;
694 int i1, j1, i2, j2, nxgrid;
695 double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV];
696 double *tmpx, *tmpy;
697 double *area_in, *area_out, min_area;
698 double lon_in_avg;
699 int n_in, n_out;
700 double xarea;
701
702
703 nx1 = *nlon_in;
704 ny1 = *nlat_in;
705 nx2 = *nlon_out;
706 ny2 = *nlat_out;
707
708 nxgrid = 0;
709 nx1p = nx1 + 1;
710 nx2p = nx2 + 1;
711
712 area_in = (double *)malloc(nx1*ny1*sizeof(double));
713 area_out = (double *)malloc(nx2*ny2*sizeof(double));
714 tmpx = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double));
715 tmpy = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double));
716 for(j2=0; j2<=ny2; j2++) for(i2=0; i2<=nx2; i2++) {
717 tmpx[j2*nx2p+i2] = lon_out[i2];
718 tmpy[j2*nx2p+i2] = lat_out[j2];
719 }
720 get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in);
721 get_grid_area(nlon_out, nlat_out, tmpx, tmpy, area_out);
722
723 free(tmpx);
724 free(tmpy);
725
726 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
727
728 ll_lon = lon_out[i2]; ll_lat = lat_out[j2];
729 ur_lon = lon_out[i2+1]; ur_lat = lat_out[j2+1];
730 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
731
732 y_in[0] = lat_in[j1*nx1p+i1];
733 y_in[1] = lat_in[j1*nx1p+i1+1];
734 y_in[2] = lat_in[(j1+1)*nx1p+i1+1];
735 y_in[3] = lat_in[(j1+1)*nx1p+i1];
736 if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat)
737 && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue;
738 if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat)
739 && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue;
740
741 x_in[0] = lon_in[j1*nx1p+i1];
742 x_in[1] = lon_in[j1*nx1p+i1+1];
743 x_in[2] = lon_in[(j1+1)*nx1p+i1+1];
744 x_in[3] = lon_in[(j1+1)*nx1p+i1];
745
746 n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2);
747 lon_in_avg = avgval_double(n_in, x_in);
748
749 if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) {
750 xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
751 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
752 if(xarea/min_area > AREA_RATIO_THRESH ) {
753 xgrid_area[nxgrid] = xarea;
754 xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg);
755 xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out );
756 i_in[nxgrid] = i1;
757 j_in[nxgrid] = j1;
758 i_out[nxgrid] = i2;
759 j_out[nxgrid] = j2;
760 ++nxgrid;
761 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
762 }
763 }
764 }
765 }
766
767 free(area_in);
768 free(area_out);
769
770 return nxgrid;
771
772 } /* create_xgrid_2dx1d_order2 */
773
774 /*******************************************************************************
775 void create_xgrid_2DX2D_order1
776 This routine generate exchange grids between two grids for the first order
777 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
778 and lon_in,lat_in, lon_out,lat_out are geographic grid location of grid cell bounds.
779 mask is on grid lon_in/lat_in.
780 *******************************************************************************/
781 int create_xgrid_2dx2d_order1_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
782 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
783 const double *mask_in, int *i_in, int *j_in, int *i_out,
784 int *j_out, double *xgrid_area)
785 {
786 int nxgrid;
787
788 nxgrid = create_xgrid_2dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in,
789 i_in, j_in, i_out, j_out, xgrid_area);
790 return nxgrid;
791
792 }
793 int create_xgrid_2dx2d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
794 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
795 const double *mask_in, int *i_in, int *j_in, int *i_out,
796 int *j_out, double *xgrid_area)
797 {
798
799 #define MAX_V 8
800 int nx1, nx2, ny1, ny2, nx1p, nx2p, nxgrid;
801 double *area_in, *area_out;
802 int nblocks =1;
803 int *istart2=NULL, *iend2=NULL;
804 int npts_left, nblks_left, pos, m, npts_my, ij;
805 double *lon_out_min_list,*lon_out_max_list,*lon_out_avg,*lat_out_min_list,*lat_out_max_list;
806 double *lon_out_list, *lat_out_list;
807 int *pnxgrid=NULL, *pstart;
808 int *pi_in=NULL, *pj_in=NULL, *pi_out=NULL, *pj_out=NULL;
809 double *pxgrid_area=NULL;
810 int *n2_list;
811 int nthreads, nxgrid_block_max;
812
813 nx1 = *nlon_in;
814 ny1 = *nlat_in;
815 nx2 = *nlon_out;
816 ny2 = *nlat_out;
817 nx1p = nx1 + 1;
818 nx2p = nx2 + 1;
819
820 area_in = (double *)malloc(nx1*ny1*sizeof(double));
821 area_out = (double *)malloc(nx2*ny2*sizeof(double));
822 get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in);
823 get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out);
824
825 nthreads = 1;
826 #if defined(_OPENMP)
827 #pragma omp parallel
828 nthreads = omp_get_num_threads();
829 #endif
830
831 nblocks = nthreads;
832
833 istart2 = (int *)malloc(nblocks*sizeof(int));
834 iend2 = (int *)malloc(nblocks*sizeof(int));
835
836 pstart = (int *)malloc(nblocks*sizeof(int));
837 pnxgrid = (int *)malloc(nblocks*sizeof(int));
838
839 nxgrid_block_max = MAXXGRID/nblocks;
840
841 for(m=0; m<nblocks; m++) {
842 pnxgrid[m] = 0;
843 pstart[m] = m*nxgrid_block_max;
844 }
845
846 if(nblocks == 1) {
847 pi_in = i_in;
848 pj_in = j_in;
849 pi_out = i_out;
850 pj_out = j_out;
851 pxgrid_area = xgrid_area;
852 }
853 else {
854 pi_in = (int *)malloc(MAXXGRID*sizeof(int));
855 pj_in = (int *)malloc(MAXXGRID*sizeof(int));
856 pi_out = (int *)malloc(MAXXGRID*sizeof(int));
857 pj_out = (int *)malloc(MAXXGRID*sizeof(int));
858 pxgrid_area = (double *)malloc(MAXXGRID*sizeof(double));
859 }
860
861 npts_left = nx2*ny2;
862 nblks_left = nblocks;
863 pos = 0;
864 for(m=0; m<nblocks; m++) {
865 istart2[m] = pos;
866 npts_my = npts_left/nblks_left;
867 iend2[m] = istart2[m] + npts_my - 1;
868 pos = iend2[m] + 1;
869 npts_left -= npts_my;
870 nblks_left--;
871 }
872
873 lon_out_min_list = (double *)malloc(nx2*ny2*sizeof(double));
874 lon_out_max_list = (double *)malloc(nx2*ny2*sizeof(double));
875 lat_out_min_list = (double *)malloc(nx2*ny2*sizeof(double));
876 lat_out_max_list = (double *)malloc(nx2*ny2*sizeof(double));
877 lon_out_avg = (double *)malloc(nx2*ny2*sizeof(double));
878 n2_list = (int *)malloc(nx2*ny2*sizeof(int));
879 lon_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double));
880 lat_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double));
881 #if defined(_OPENMP)
882 #pragma omp parallel for default(none) shared(nx2,ny2,nx2p,lon_out,lat_out,lat_out_min_list, \
883 lat_out_max_list,lon_out_min_list,lon_out_max_list, \
884 lon_out_avg,n2_list,lon_out_list,lat_out_list)
885 #endif
886 for(ij=0; ij<nx2*ny2; ij++){
887 int i2, j2, n, n0, n1, n2, n3, n2_in, l;
888 double x2_in[MV], y2_in[MV];
889 i2 = ij%nx2;
890 j2 = ij/nx2;
891 n = j2*nx2+i2;
892 n0 = j2*nx2p+i2; n1 = j2*nx2p+i2+1;
893 n2 = (j2+1)*nx2p+i2+1; n3 = (j2+1)*nx2p+i2;
894 x2_in[0] = lon_out[n0]; y2_in[0] = lat_out[n0];
895 x2_in[1] = lon_out[n1]; y2_in[1] = lat_out[n1];
896 x2_in[2] = lon_out[n2]; y2_in[2] = lat_out[n2];
897 x2_in[3] = lon_out[n3]; y2_in[3] = lat_out[n3];
898
899 lat_out_min_list[n] = minval_double(4, y2_in);
900 lat_out_max_list[n] = maxval_double(4, y2_in);
901 n2_in = fix_lon(x2_in, y2_in, 4, M_PI);
902 if(n2_in > MAX_V) error_handler("create_xgrid.c: n2_in is greater than MAX_V");
903 lon_out_min_list[n] = minval_double(n2_in, x2_in);
904 lon_out_max_list[n] = maxval_double(n2_in, x2_in);
905 lon_out_avg[n] = avgval_double(n2_in, x2_in);
906 n2_list[n] = n2_in;
907 for(l=0; l<n2_in; l++) {
908 lon_out_list[n*MAX_V+l] = x2_in[l];
909 lat_out_list[n*MAX_V+l] = y2_in[l];
910 }
911 }
912
913 nxgrid = 0;
914
915 #if defined(_OPENMP)
916 #pragma omp parallel for default(none) shared(nblocks,nx1,ny1,nx1p,mask_in,lon_in,lat_in, \
917 istart2,iend2,nx2,lat_out_min_list,lat_out_max_list, \
918 n2_list,lon_out_list,lat_out_list,lon_out_min_list, \
919 lon_out_max_list,lon_out_avg,area_in,area_out, \
920 pxgrid_area,pnxgrid,pi_in,pj_in,pi_out,pj_out,pstart,nthreads)
921 #endif
922 for(m=0; m<nblocks; m++) {
923 int i1, j1, ij;
924 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
925 int n0, n1, n2, n3, l,n1_in;
926 double lat_in_min,lat_in_max,lon_in_min,lon_in_max,lon_in_avg;
927 double x1_in[MV], y1_in[MV], x_out[MV], y_out[MV];
928
929 n0 = j1*nx1p+i1; n1 = j1*nx1p+i1+1;
930 n2 = (j1+1)*nx1p+i1+1; n3 = (j1+1)*nx1p+i1;
931 x1_in[0] = lon_in[n0]; y1_in[0] = lat_in[n0];
932 x1_in[1] = lon_in[n1]; y1_in[1] = lat_in[n1];
933 x1_in[2] = lon_in[n2]; y1_in[2] = lat_in[n2];
934 x1_in[3] = lon_in[n3]; y1_in[3] = lat_in[n3];
935 lat_in_min = minval_double(4, y1_in);
936 lat_in_max = maxval_double(4, y1_in);
937 n1_in = fix_lon(x1_in, y1_in, 4, M_PI);
938 lon_in_min = minval_double(n1_in, x1_in);
939 lon_in_max = maxval_double(n1_in, x1_in);
940 lon_in_avg = avgval_double(n1_in, x1_in);
941 for(ij=istart2[m]; ij<=iend2[m]; ij++) {
942 int n_out, i2, j2, n2_in;
943 double xarea, dx, lon_out_min, lon_out_max;
944 double x2_in[MAX_V], y2_in[MAX_V];
945
946 i2 = ij%nx2;
947 j2 = ij/nx2;
948
949 if(lat_out_min_list[ij] >= lat_in_max || lat_out_max_list[ij] <= lat_in_min ) continue;
950 /* adjust x2_in according to lon_in_avg*/
951 n2_in = n2_list[ij];
952 for(l=0; l<n2_in; l++) {
953 x2_in[l] = lon_out_list[ij*MAX_V+l];
954 y2_in[l] = lat_out_list[ij*MAX_V+l];
955 }
956 lon_out_min = lon_out_min_list[ij];
957 lon_out_max = lon_out_max_list[ij];
958 dx = lon_out_avg[ij] - lon_in_avg;
959 if(dx < -M_PI ) {
960 lon_out_min += TPI;
961 lon_out_max += TPI;
962 for (l=0; l<n2_in; l++) x2_in[l] += TPI;
963 }
964 else if (dx > M_PI) {
965 lon_out_min -= TPI;
966 lon_out_max -= TPI;
967 for (l=0; l<n2_in; l++) x2_in[l] -= TPI;
968 }
969
970 /* x2_in should in the same range as x1_in after lon_fix, so no need to
971 consider cyclic condition
972 */
973 if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue;
974 if ( (n_out = clip_2dx2d( x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out )) > 0) {
975 double min_area;
976 int nn;
977 xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
978 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
979 if( xarea/min_area > AREA_RATIO_THRESH ) {
980 pnxgrid[m]++;
981 if(pnxgrid[m]>= MAXXGRID/nthreads)
982 error_handler("nxgrid is greater than MAXXGRID/nthreads, increase MAXXGRID, decrease nthreads, or increase number of MPI ranks");
983 nn = pstart[m] + pnxgrid[m]-1;
984
985 pxgrid_area[nn] = xarea;
986 pi_in[nn] = i1;
987 pj_in[nn] = j1;
988 pi_out[nn] = i2;
989 pj_out[nn] = j2;
990 }
991
992 }
993
994 }
995 }
996 }
997
998 /*copy data if nblocks > 1 */
999 if(nblocks == 1) {
1000 nxgrid = pnxgrid[0];
1001 pi_in = NULL;
1002 pj_in = NULL;
1003 pi_out = NULL;
1004 pj_out = NULL;
1005 pxgrid_area = NULL;
1006 }
1007 else {
1008 int nn, i;
1009 nxgrid = 0;
1010 for(m=0; m<nblocks; m++) {
1011 for(i=0; i<pnxgrid[m]; i++) {
1012 nn = pstart[m] + i;
1013 i_in[nxgrid] = pi_in[nn];
1014 j_in[nxgrid] = pj_in[nn];
1015 i_out[nxgrid] = pi_out[nn];
1016 j_out[nxgrid] = pj_out[nn];
1017 xgrid_area[nxgrid] = pxgrid_area[nn];
1018 nxgrid++;
1019 }
1020 }
1021 free(pi_in);
1022 free(pj_in);
1023 free(pi_out);
1024 free(pj_out);
1025 free(pxgrid_area);
1026 }
1027
1028 free(area_in);
1029 free(area_out);
1030 free(lon_out_min_list);
1031 free(lon_out_max_list);
1032 free(lat_out_min_list);
1033 free(lat_out_max_list);
1034 free(lon_out_avg);
1035 free(n2_list);
1036 free(lon_out_list);
1037 free(lat_out_list);
1038
1039 return nxgrid;
1040
1041 }/* get_xgrid_2Dx2D_order1 */
1042
1043 /********************************************************************************
1044 void create_xgrid_2dx1d_order2
1045 This routine generate exchange grids between two grids for the second order
1046 conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell
1047 and lon_in,lat_in, lon_out,lat_out are geographic grid location of grid cell bounds.
1048 mask is on grid lon_in/lat_in.
1049 ********************************************************************************/
1050 int create_xgrid_2dx2d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
1051 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1052 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
1053 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1054 {
1055 int nxgrid;
1056 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,
1057 j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat);
1058 return nxgrid;
1059
1060 }
1061 int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
1062 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1063 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
1064 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1065 {
1066
1067 #define MAX_V 8
1068 int nx1, nx2, ny1, ny2, nx1p, nx2p, nxgrid;
1069 double *area_in, *area_out;
1070 int nblocks =1;
1071 int *istart2=NULL, *iend2=NULL;
1072 int npts_left, nblks_left, pos, m, npts_my, ij;
1073 double *lon_out_min_list,*lon_out_max_list,*lon_out_avg,*lat_out_min_list,*lat_out_max_list;
1074 double *lon_out_list, *lat_out_list;
1075 int *pnxgrid=NULL, *pstart;
1076 int *pi_in=NULL, *pj_in=NULL, *pi_out=NULL, *pj_out=NULL;
1077 double *pxgrid_area=NULL, *pxgrid_clon=NULL, *pxgrid_clat=NULL;
1078 int *n2_list;
1079 int nthreads, nxgrid_block_max;
1080
1081 nx1 = *nlon_in;
1082 ny1 = *nlat_in;
1083 nx2 = *nlon_out;
1084 ny2 = *nlat_out;
1085 nx1p = nx1 + 1;
1086 nx2p = nx2 + 1;
1087
1088 area_in = (double *)malloc(nx1*ny1*sizeof(double));
1089 area_out = (double *)malloc(nx2*ny2*sizeof(double));
1090 get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in);
1091 get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out);
1092
1093 nthreads = 1;
1094 #if defined(_OPENMP)
1095 #pragma omp parallel
1096 nthreads = omp_get_num_threads();
1097 #endif
1098
1099 nblocks = nthreads;
1100
1101 istart2 = (int *)malloc(nblocks*sizeof(int));
1102 iend2 = (int *)malloc(nblocks*sizeof(int));
1103
1104 pstart = (int *)malloc(nblocks*sizeof(int));
1105 pnxgrid = (int *)malloc(nblocks*sizeof(int));
1106
1107 nxgrid_block_max = MAXXGRID/nblocks;
1108
1109 for(m=0; m<nblocks; m++) {
1110 pnxgrid[m] = 0;
1111 pstart[m] = m*nxgrid_block_max;
1112 }
1113
1114 if(nblocks == 1) {
1115 pi_in = i_in;
1116 pj_in = j_in;
1117 pi_out = i_out;
1118 pj_out = j_out;
1119 pxgrid_area = xgrid_area;
1120 pxgrid_clon = xgrid_clon;
1121 pxgrid_clat = xgrid_clat;
1122 }
1123 else {
1124 pi_in = (int *)malloc(MAXXGRID*sizeof(int));
1125 pj_in = (int *)malloc(MAXXGRID*sizeof(int));
1126 pi_out = (int *)malloc(MAXXGRID*sizeof(int));
1127 pj_out = (int *)malloc(MAXXGRID*sizeof(int));
1128 pxgrid_area = (double *)malloc(MAXXGRID*sizeof(double));
1129 pxgrid_clon = (double *)malloc(MAXXGRID*sizeof(double));
1130 pxgrid_clat = (double *)malloc(MAXXGRID*sizeof(double));
1131 }
1132
1133 npts_left = nx2*ny2;
1134 nblks_left = nblocks;
1135 pos = 0;
1136 for(m=0; m<nblocks; m++) {
1137 istart2[m] = pos;
1138 npts_my = npts_left/nblks_left;
1139 iend2[m] = istart2[m] + npts_my - 1;
1140 pos = iend2[m] + 1;
1141 npts_left -= npts_my;
1142 nblks_left--;
1143 }
1144
1145 lon_out_min_list = (double *)malloc(nx2*ny2*sizeof(double));
1146 lon_out_max_list = (double *)malloc(nx2*ny2*sizeof(double));
1147 lat_out_min_list = (double *)malloc(nx2*ny2*sizeof(double));
1148 lat_out_max_list = (double *)malloc(nx2*ny2*sizeof(double));
1149 lon_out_avg = (double *)malloc(nx2*ny2*sizeof(double));
1150 n2_list = (int *)malloc(nx2*ny2*sizeof(int));
1151 lon_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double));
1152 lat_out_list = (double *)malloc(MAX_V*nx2*ny2*sizeof(double));
1153 #if defined(_OPENMP)
1154 #pragma omp parallel for default(none) shared(nx2,ny2,nx2p,lon_out,lat_out,lat_out_min_list, \
1155 lat_out_max_list,lon_out_min_list,lon_out_max_list, \
1156 lon_out_avg,n2_list,lon_out_list,lat_out_list)
1157 #endif
1158 for(ij=0; ij<nx2*ny2; ij++){
1159 int i2, j2, n, n0, n1, n2, n3, n2_in, l;
1160 double x2_in[MV], y2_in[MV];
1161 i2 = ij%nx2;
1162 j2 = ij/nx2;
1163 n = j2*nx2+i2;
1164 n0 = j2*nx2p+i2; n1 = j2*nx2p+i2+1;
1165 n2 = (j2+1)*nx2p+i2+1; n3 = (j2+1)*nx2p+i2;
1166 x2_in[0] = lon_out[n0]; y2_in[0] = lat_out[n0];
1167 x2_in[1] = lon_out[n1]; y2_in[1] = lat_out[n1];
1168 x2_in[2] = lon_out[n2]; y2_in[2] = lat_out[n2];
1169 x2_in[3] = lon_out[n3]; y2_in[3] = lat_out[n3];
1170
1171 lat_out_min_list[n] = minval_double(4, y2_in);
1172 lat_out_max_list[n] = maxval_double(4, y2_in);
1173 n2_in = fix_lon(x2_in, y2_in, 4, M_PI);
1174 if(n2_in > MAX_V) error_handler("create_xgrid.c: n2_in is greater than MAX_V");
1175 lon_out_min_list[n] = minval_double(n2_in, x2_in);
1176 lon_out_max_list[n] = maxval_double(n2_in, x2_in);
1177 lon_out_avg[n] = avgval_double(n2_in, x2_in);
1178 n2_list[n] = n2_in;
1179 for(l=0; l<n2_in; l++) {
1180 lon_out_list[n*MAX_V+l] = x2_in[l];
1181 lat_out_list[n*MAX_V+l] = y2_in[l];
1182 }
1183 }
1184
1185 nxgrid = 0;
1186
1187 #if defined(_OPENMP)
1188 #pragma omp parallel for default(none) shared(nblocks,nx1,ny1,nx1p,mask_in,lon_in,lat_in, \
1189 istart2,iend2,nx2,lat_out_min_list,lat_out_max_list, \
1190 n2_list,lon_out_list,lat_out_list,lon_out_min_list, \
1191 lon_out_max_list,lon_out_avg,area_in,area_out, \
1192 pxgrid_area,pnxgrid,pxgrid_clon,pxgrid_clat,pi_in, \
1193 pj_in,pi_out,pj_out,pstart,nthreads)
1194 #endif
1195 for(m=0; m<nblocks; m++) {
1196 int i1, j1, ij;
1197 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
1198 int n0, n1, n2, n3, l,n1_in;
1199 double lat_in_min,lat_in_max,lon_in_min,lon_in_max,lon_in_avg;
1200 double x1_in[MV], y1_in[MV], x_out[MV], y_out[MV];
1201
1202 n0 = j1*nx1p+i1; n1 = j1*nx1p+i1+1;
1203 n2 = (j1+1)*nx1p+i1+1; n3 = (j1+1)*nx1p+i1;
1204 x1_in[0] = lon_in[n0]; y1_in[0] = lat_in[n0];
1205 x1_in[1] = lon_in[n1]; y1_in[1] = lat_in[n1];
1206 x1_in[2] = lon_in[n2]; y1_in[2] = lat_in[n2];
1207 x1_in[3] = lon_in[n3]; y1_in[3] = lat_in[n3];
1208 lat_in_min = minval_double(4, y1_in);
1209 lat_in_max = maxval_double(4, y1_in);
1210 n1_in = fix_lon(x1_in, y1_in, 4, M_PI);
1211 lon_in_min = minval_double(n1_in, x1_in);
1212 lon_in_max = maxval_double(n1_in, x1_in);
1213 lon_in_avg = avgval_double(n1_in, x1_in);
1214 for(ij=istart2[m]; ij<=iend2[m]; ij++) {
1215 int n_out, i2, j2, n2_in;
1216 double xarea, dx, lon_out_min, lon_out_max;
1217 double x2_in[MAX_V], y2_in[MAX_V];
1218
1219 i2 = ij%nx2;
1220 j2 = ij/nx2;
1221
1222 if(lat_out_min_list[ij] >= lat_in_max || lat_out_max_list[ij] <= lat_in_min ) continue;
1223 /* adjust x2_in according to lon_in_avg*/
1224 n2_in = n2_list[ij];
1225 for(l=0; l<n2_in; l++) {
1226 x2_in[l] = lon_out_list[ij*MAX_V+l];
1227 y2_in[l] = lat_out_list[ij*MAX_V+l];
1228 }
1229 lon_out_min = lon_out_min_list[ij];
1230 lon_out_max = lon_out_max_list[ij];
1231 dx = lon_out_avg[ij] - lon_in_avg;
1232 if(dx < -M_PI ) {
1233 lon_out_min += TPI;
1234 lon_out_max += TPI;
1235 for (l=0; l<n2_in; l++) x2_in[l] += TPI;
1236 }
1237 else if (dx > M_PI) {
1238 lon_out_min -= TPI;
1239 lon_out_max -= TPI;
1240 for (l=0; l<n2_in; l++) x2_in[l] -= TPI;
1241 }
1242
1243 /* x2_in should in the same range as x1_in after lon_fix, so no need to
1244 consider cyclic condition
1245 */
1246 if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue;
1247 if ( (n_out = clip_2dx2d( x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out )) > 0) {
1248 double min_area;
1249 int nn;
1250 xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1];
1251 min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]);
1252 if( xarea/min_area > AREA_RATIO_THRESH ) {
1253 pnxgrid[m]++;
1254 if(pnxgrid[m]>= MAXXGRID/nthreads)
1255 error_handler("nxgrid is greater than MAXXGRID/nthreads, increase MAXXGRID, decrease nthreads, or increase number of MPI ranks");
1256 nn = pstart[m] + pnxgrid[m]-1;
1257 pxgrid_area[nn] = xarea;
1258 pxgrid_clon[nn] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg);
1259 pxgrid_clat[nn] = poly_ctrlat (x_out, y_out, n_out );
1260 pi_in[nn] = i1;
1261 pj_in[nn] = j1;
1262 pi_out[nn] = i2;
1263 pj_out[nn] = j2;
1264 }
1265 }
1266 }
1267 }
1268 }
1269
1270 /*copy data if nblocks > 1 */
1271 if(nblocks == 1) {
1272 nxgrid = pnxgrid[0];
1273 pi_in = NULL;
1274 pj_in = NULL;
1275 pi_out = NULL;
1276 pj_out = NULL;
1277 pxgrid_area = NULL;
1278 pxgrid_clon = NULL;
1279 pxgrid_clat = NULL;
1280 }
1281 else {
1282 int nn, i;
1283 nxgrid = 0;
1284 for(m=0; m<nblocks; m++) {
1285 for(i=0; i<pnxgrid[m]; i++) {
1286 nn = pstart[m] + i;
1287 i_in[nxgrid] = pi_in[nn];
1288 j_in[nxgrid] = pj_in[nn];
1289 i_out[nxgrid] = pi_out[nn];
1290 j_out[nxgrid] = pj_out[nn];
1291 xgrid_area[nxgrid] = pxgrid_area[nn];
1292 xgrid_clon[nxgrid] = pxgrid_clon[nn];
1293 xgrid_clat[nxgrid] = pxgrid_clat[nn];
1294 nxgrid++;
1295 }
1296 }
1297 free(pi_in);
1298 free(pj_in);
1299 free(pi_out);
1300 free(pj_out);
1301 free(pxgrid_area);
1302 free(pxgrid_clon);
1303 free(pxgrid_clat);
1304 }
1305
1306 free(area_in);
1307 free(area_out);
1308 free(lon_out_min_list);
1309 free(lon_out_max_list);
1310 free(lat_out_min_list);
1311 free(lat_out_max_list);
1312 free(lon_out_avg);
1313 free(n2_list);
1314 free(lon_out_list);
1315 free(lat_out_list);
1316
1317 return nxgrid;
1318
1319 }/* get_xgrid_2Dx2D_order2 */
1320
1321
1322 /*******************************************************************************
1323 Sutherland-Hodgeman algorithm sequentially clips parts outside 4 boundaries
1324 *******************************************************************************/
1325
1326 int clip(const double lon_in[], const double lat_in[], int n_in, double ll_lon, double ll_lat,
1327 double ur_lon, double ur_lat, double lon_out[], double lat_out[])
1328 {
1329 double x_tmp[MV], y_tmp[MV], x_last, y_last;
1330 int i_in, i_out, n_out, inside_last, inside;
1331
1332 /* clip polygon with LEFT boundary - clip V_IN to V_TMP */
1333 x_last = lon_in[n_in-1];
1334 y_last = lat_in[n_in-1];
1335 inside_last = (x_last >= ll_lon);
1336 for (i_in=0,i_out=0;i_in<n_in;i_in++) {
1337
1338 /* if crossing LEFT boundary - output intersection */
1339 if ((inside=(lon_in[i_in] >= ll_lon))!=inside_last) {
1340 x_tmp[i_out] = ll_lon;
1341 y_tmp[i_out++] = y_last + (ll_lon - x_last) * (lat_in[i_in] - y_last) / (lon_in[i_in] - x_last);
1342 }
1343
1344 /* if "to" point is right of LEFT boundary, output it */
1345 if (inside) {
1346 x_tmp[i_out] = lon_in[i_in];
1347 y_tmp[i_out++] = lat_in[i_in];
1348 }
1349 x_last = lon_in[i_in];
1350 y_last = lat_in[i_in];
1351 inside_last = inside;
1352 }
1353 if (!(n_out=i_out)) return(0);
1354
1355 /* clip polygon with RIGHT boundary - clip V_TMP to V_OUT */
1356 x_last = x_tmp[n_out-1];
1357 y_last = y_tmp[n_out-1];
1358 inside_last = (x_last <= ur_lon);
1359 for (i_in=0,i_out=0;i_in<n_out;i_in++) {
1360
1361 /* if crossing RIGHT boundary - output intersection */
1362 if ((inside=(x_tmp[i_in] <= ur_lon))!=inside_last) {
1363 lon_out[i_out] = ur_lon;
1364 lat_out[i_out++] = y_last + (ur_lon - x_last) * (y_tmp[i_in] - y_last)
1365 / (x_tmp[i_in] - x_last);
1366 }
1367
1368 /* if "to" point is left of RIGHT boundary, output it */
1369 if (inside) {
1370 lon_out[i_out] = x_tmp[i_in];
1371 lat_out[i_out++] = y_tmp[i_in];
1372 }
1373
1374 x_last = x_tmp[i_in];
1375 y_last = y_tmp[i_in];
1376 inside_last = inside;
1377 }
1378 if (!(n_out=i_out)) return(0);
1379
1380 /* clip polygon with BOTTOM boundary - clip V_OUT to V_TMP */
1381 x_last = lon_out[n_out-1];
1382 y_last = lat_out[n_out-1];
1383 inside_last = (y_last >= ll_lat);
1384 for (i_in=0,i_out=0;i_in<n_out;i_in++) {
1385
1386 /* if crossing BOTTOM boundary - output intersection */
1387 if ((inside=(lat_out[i_in] >= ll_lat))!=inside_last) {
1388 y_tmp[i_out] = ll_lat;
1389 x_tmp[i_out++] = x_last + (ll_lat - y_last) * (lon_out[i_in] - x_last) / (lat_out[i_in] - y_last);
1390 }
1391
1392 /* if "to" point is above BOTTOM boundary, output it */
1393 if (inside) {
1394 x_tmp[i_out] = lon_out[i_in];
1395 y_tmp[i_out++] = lat_out[i_in];
1396 }
1397 x_last = lon_out[i_in];
1398 y_last = lat_out[i_in];
1399 inside_last = inside;
1400 }
1401 if (!(n_out=i_out)) return(0);
1402
1403 /* clip polygon with TOP boundary - clip V_TMP to V_OUT */
1404 x_last = x_tmp[n_out-1];
1405 y_last = y_tmp[n_out-1];
1406 inside_last = (y_last <= ur_lat);
1407 for (i_in=0,i_out=0;i_in<n_out;i_in++) {
1408
1409 /* if crossing TOP boundary - output intersection */
1410 if ((inside=(y_tmp[i_in] <= ur_lat))!=inside_last) {
1411 lat_out[i_out] = ur_lat;
1412 lon_out[i_out++] = x_last + (ur_lat - y_last) * (x_tmp[i_in] - x_last) / (y_tmp[i_in] - y_last);
1413 }
1414
1415 /* if "to" point is below TOP boundary, output it */
1416 if (inside) {
1417 lon_out[i_out] = x_tmp[i_in];
1418 lat_out[i_out++] = y_tmp[i_in];
1419 }
1420 x_last = x_tmp[i_in];
1421 y_last = y_tmp[i_in];
1422 inside_last = inside;
1423 }
1424 return(i_out);
1425 } /* clip */
1426
1427
1428 /*******************************************************************************
1429 Revise Sutherland-Hodgeman algorithm to find the vertices of the overlapping
1430 between any two grid boxes. It return the number of vertices for the exchange grid.
1431 *******************************************************************************/
1432
1433 int clip_2dx2d(const double lon1_in[], const double lat1_in[], int n1_in,
1434 const double lon2_in[], const double lat2_in[], int n2_in,
1435 double lon_out[], double lat_out[])
1436 {
1437 double lon_tmp[MV], lat_tmp[MV];
1438 double x1_0, y1_0, x1_1, y1_1, x2_0, y2_0, x2_1, y2_1;
1439 double dx1, dy1, dx2, dy2, determ, ds1, ds2;
1440 int i_out, n_out, inside_last, inside, i1, i2;
1441
1442 /* clip polygon with each boundary of the polygon */
1443 /* We treat lon1_in/lat1_in as clip polygon and lon2_in/lat2_in as subject polygon */
1444 n_out = n1_in;
1445 for(i1=0; i1<n1_in; i1++) {
1446 lon_tmp[i1] = lon1_in[i1];
1447 lat_tmp[i1] = lat1_in[i1];
1448 }
1449 x2_0 = lon2_in[n2_in-1];
1450 y2_0 = lat2_in[n2_in-1];
1451 for(i2=0; i2<n2_in; i2++) {
1452 x2_1 = lon2_in[i2];
1453 y2_1 = lat2_in[i2];
1454 x1_0 = lon_tmp[n_out-1];
1455 y1_0 = lat_tmp[n_out-1];
1456 inside_last = inside_edge( x2_0, y2_0, x2_1, y2_1, x1_0, y1_0);
1457 for(i1=0, i_out=0; i1<n_out; i1++) {
1458 x1_1 = lon_tmp[i1];
1459 y1_1 = lat_tmp[i1];
1460 if((inside = inside_edge(x2_0, y2_0, x2_1, y2_1, x1_1, y1_1)) != inside_last ) {
1461 /* there is intersection, the line between <x1_0,y1_0> and <x1_1,y1_1>
1462 should not parallel to the line between <x2_0,y2_0> and <x2_1,y2_1>
1463 may need to consider truncation error */
1464 dy1 = y1_1-y1_0;
1465 dy2 = y2_1-y2_0;
1466 dx1 = x1_1-x1_0;
1467 dx2 = x2_1-x2_0;
1468 ds1 = y1_0*x1_1 - y1_1*x1_0;
1469 ds2 = y2_0*x2_1 - y2_1*x2_0;
1470 determ = dy2*dx1 - dy1*dx2;
1471 if(fabs(determ) < EPSLN30) {
1472 error_handler("the line between <x1_0,y1_0> and <x1_1,y1_1> should not parallel to "
1473 "the line between <x2_0,y2_0> and <x2_1,y2_1>");
1474 }
1475 lon_out[i_out] = (dx2*ds1 - dx1*ds2)/determ;
1476 lat_out[i_out++] = (dy2*ds1 - dy1*ds2)/determ;
1477
1478
1479 }
1480 if(inside) {
1481 lon_out[i_out] = x1_1;
1482 lat_out[i_out++] = y1_1;
1483 }
1484 x1_0 = x1_1;
1485 y1_0 = y1_1;
1486 inside_last = inside;
1487 }
1488 if(!(n_out=i_out)) return 0;
1489 for(i1=0; i1<n_out; i1++) {
1490 lon_tmp[i1] = lon_out[i1];
1491 lat_tmp[i1] = lat_out[i1];
1492 }
1493 /* shift the starting point */
1494 x2_0 = x2_1;
1495 y2_0 = y2_1;
1496 }
1497 return(n_out);
1498 } /* clip */
1499
1500 /*#define debug_test_create_xgrid*/
1501
1502 int create_xgrid_great_circle_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
1503 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1504 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
1505 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1506 {
1507 int nxgrid;
1508 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out,
1509 mask_in, i_in, j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat);
1510
1511 return nxgrid;
1512 }
1513
1514 int create_xgrid_great_circle(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
1515 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1516 const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
1517 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1518 {
1519
1520 int nx1, nx2, ny1, ny2, nx1p, nx2p, ny1p, ny2p, nxgrid, n1_in, n2_in;
1521 int n0, n1, n2, n3, i1, j1, i2, j2;
1522 double x1_in[MV], y1_in[MV], z1_in[MV];
1523 double x2_in[MV], y2_in[MV], z2_in[MV];
1524 double x_out[MV], y_out[MV], z_out[MV];
1525 double *x1=NULL, *y1=NULL, *z1=NULL;
1526 double *x2=NULL, *y2=NULL, *z2=NULL;
1527
1528 double *area1, *area2, min_area;
1529
1530 nx1 = *nlon_in;
1531 ny1 = *nlat_in;
1532 nx2 = *nlon_out;
1533 ny2 = *nlat_out;
1534 nxgrid = 0;
1535 nx1p = nx1 + 1;
1536 nx2p = nx2 + 1;
1537 ny1p = ny1 + 1;
1538 ny2p = ny2 + 1;
1539
1540 /* first convert lon-lat to cartesian coordinates */
1541 x1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1542 y1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1543 z1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1544 x2 = (double *)malloc(nx2p*ny2p*sizeof(double));
1545 y2 = (double *)malloc(nx2p*ny2p*sizeof(double));
1546 z2 = (double *)malloc(nx2p*ny2p*sizeof(double));
1547
1548 latlon2xyz(nx1p*ny1p, lon_in, lat_in, x1, y1, z1);
1549 latlon2xyz(nx2p*ny2p, lon_out, lat_out, x2, y2, z2);
1550
1551 area1 = (double *)malloc(nx1*ny1*sizeof(double));
1552 area2 = (double *)malloc(nx2*ny2*sizeof(double));
1553 get_grid_great_circle_area(nlon_in, nlat_in, lon_in, lat_in, area1);
1554 get_grid_great_circle_area(nlon_out, nlat_out, lon_out, lat_out, area2);
1555 n1_in = 4;
1556 n2_in = 4;
1557
1558 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
1559 /* clockwise */
1560 n0 = j1*nx1p+i1; n1 = (j1+1)*nx1p+i1;
1561 n2 = (j1+1)*nx1p+i1+1; n3 = j1*nx1p+i1+1;
1562 x1_in[0] = x1[n0]; y1_in[0] = y1[n0]; z1_in[0] = z1[n0];
1563 x1_in[1] = x1[n1]; y1_in[1] = y1[n1]; z1_in[1] = z1[n1];
1564 x1_in[2] = x1[n2]; y1_in[2] = y1[n2]; z1_in[2] = z1[n2];
1565 x1_in[3] = x1[n3]; y1_in[3] = y1[n3]; z1_in[3] = z1[n3];
1566
1567 for(j2=0; j2<ny2; j2++) for(i2=0; i2<nx2; i2++) {
1568 int n_out;
1569 double xarea;
1570
1571 n0 = j2*nx2p+i2; n1 = (j2+1)*nx2p+i2;
1572 n2 = (j2+1)*nx2p+i2+1; n3 = j2*nx2p+i2+1;
1573 x2_in[0] = x2[n0]; y2_in[0] = y2[n0]; z2_in[0] = z2[n0];
1574 x2_in[1] = x2[n1]; y2_in[1] = y2[n1]; z2_in[1] = z2[n1];
1575 x2_in[2] = x2[n2]; y2_in[2] = y2[n2]; z2_in[2] = z2[n2];
1576 x2_in[3] = x2[n3]; y2_in[3] = y2[n3]; z2_in[3] = z2[n3];
1577
1578 if ( (n_out = clip_2dx2d_great_circle( x1_in, y1_in, z1_in, n1_in, x2_in, y2_in, z2_in, n2_in,
1579 x_out, y_out, z_out)) > 0) {
1580 xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1];
1581 min_area = min(area1[j1*nx1+i1], area2[j2*nx2+i2]);
1582 if( xarea/min_area > AREA_RATIO_THRESH ) {
1583 #ifdef debug_test_create_xgrid
1584 printf("(i2,j2)=(%d,%d), (i1,j1)=(%d,%d), xarea=%g\n", i2, j2, i1, j1, xarea);
1585 #endif
1586 xgrid_area[nxgrid] = xarea;
1587 xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */
1588 xgrid_clat[nxgrid] = 0;
1589 i_in[nxgrid] = i1;
1590 j_in[nxgrid] = j1;
1591 i_out[nxgrid] = i2;
1592 j_out[nxgrid] = j2;
1593 ++nxgrid;
1594 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
1595 }
1596 }
1597 }
1598 }
1599
1600
1601 free(area1);
1602 free(area2);
1603
1604 free(x1);
1605 free(y1);
1606 free(z1);
1607 free(x2);
1608 free(y2);
1609 free(z2);
1610
1611 return nxgrid;
1612
1613 }/* create_xgrid_great_circle */
1614
1615 int create_xgrid_great_circle_ug_(const int *nlon_in, const int *nlat_in, const int *npts_out,
1616 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1617 const double *mask_in, int *i_in, int *j_in, int *l_out,
1618 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1619 {
1620 int nxgrid;
1621 nxgrid = create_xgrid_great_circle_ug(nlon_in, nlat_in, npts_out, lon_in, lat_in, lon_out, lat_out,
1622 mask_in, i_in, j_in, l_out, xgrid_area, xgrid_clon, xgrid_clat);
1623
1624 return nxgrid;
1625 }
1626
1627 int create_xgrid_great_circle_ug(const int *nlon_in, const int *nlat_in, const int *npts_out,
1628 const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
1629 const double *mask_in, int *i_in, int *j_in, int *l_out,
1630 double *xgrid_area, double *xgrid_clon, double *xgrid_clat)
1631 {
1632
1633 int nx1, ny1, npts2, nx1p, ny1p, nxgrid, n1_in, n2_in, nv;
1634 int n0, n1, n2, n3, i1, j1, l2;
1635 double x1_in[MV], y1_in[MV], z1_in[MV];
1636 double x2_in[MV], y2_in[MV], z2_in[MV];
1637 double x_out[MV], y_out[MV], z_out[MV];
1638 double *x1=NULL, *y1=NULL, *z1=NULL;
1639 double *x2=NULL, *y2=NULL, *z2=NULL;
1640
1641 double *area1, *area2, min_area;
1642
1643 nx1 = *nlon_in;
1644 ny1 = *nlat_in;
1645 nv = 4;
1646 npts2 = *npts_out;
1647 nxgrid = 0;
1648 nx1p = nx1 + 1;
1649 ny1p = ny1 + 1;
1650
1651 /* first convert lon-lat to cartesian coordinates */
1652 x1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1653 y1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1654 z1 = (double *)malloc(nx1p*ny1p*sizeof(double));
1655 x2 = (double *)malloc(npts2*nv*sizeof(double));
1656 y2 = (double *)malloc(npts2*nv*sizeof(double));
1657 z2 = (double *)malloc(npts2*nv*sizeof(double));
1658
1659 latlon2xyz(nx1p*ny1p, lon_in, lat_in, x1, y1, z1);
1660 latlon2xyz(npts2*nv, lon_out, lat_out, x2, y2, z2);
1661
1662 area1 = (double *)malloc(nx1*ny1*sizeof(double));
1663 area2 = (double *)malloc(npts2*sizeof(double));
1664 get_grid_great_circle_area(nlon_in, nlat_in, lon_in, lat_in, area1);
1665 get_grid_great_circle_area_ug(npts_out, lon_out, lat_out, area2);
1666 n1_in = 4;
1667 n2_in = 4;
1668
1669 for(j1=0; j1<ny1; j1++) for(i1=0; i1<nx1; i1++) if( mask_in[j1*nx1+i1] > MASK_THRESH ) {
1670 /* clockwise */
1671 n0 = j1*nx1p+i1; n1 = (j1+1)*nx1p+i1;
1672 n2 = (j1+1)*nx1p+i1+1; n3 = j1*nx1p+i1+1;
1673 x1_in[0] = x1[n0]; y1_in[0] = y1[n0]; z1_in[0] = z1[n0];
1674 x1_in[1] = x1[n1]; y1_in[1] = y1[n1]; z1_in[1] = z1[n1];
1675 x1_in[2] = x1[n2]; y1_in[2] = y1[n2]; z1_in[2] = z1[n2];
1676 x1_in[3] = x1[n3]; y1_in[3] = y1[n3]; z1_in[3] = z1[n3];
1677
1678 for(l2=0; l2<npts2; l2++) {
1679 int n_out;
1680 double xarea;
1681
1682 n0 = l2*nv; n1 = l2*nv+1;
1683 n2 = l2*nv+2; n3 = l2*nv+3;
1684 x2_in[0] = x2[n0]; y2_in[0] = y2[n0]; z2_in[0] = z2[n0];
1685 x2_in[1] = x2[n1]; y2_in[1] = y2[n1]; z2_in[1] = z2[n1];
1686 x2_in[2] = x2[n2]; y2_in[2] = y2[n2]; z2_in[2] = z2[n2];
1687 x2_in[3] = x2[n3]; y2_in[3] = y2[n3]; z2_in[3] = z2[n3];
1688
1689 if ( (n_out = clip_2dx2d_great_circle( x1_in, y1_in, z1_in, n1_in, x2_in, y2_in, z2_in, n2_in,
1690 x_out, y_out, z_out)) > 0) {
1691 xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1];
1692 min_area = min(area1[j1*nx1+i1], area2[l2]);
1693 if( xarea/min_area > AREA_RATIO_THRESH ) {
1694 #ifdef debug_test_create_xgrid
1695 printf("(l2)=(%d,%d), (i1,j1)=(%d,%d), xarea=%g\n", l2, i1, j1, xarea);
1696 #endif
1697 xgrid_area[nxgrid] = xarea;
1698 xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */
1699 xgrid_clat[nxgrid] = 0;
1700 i_in[nxgrid] = i1;
1701 j_in[nxgrid] = j1;
1702 l_out[nxgrid] = l2;
1703 ++nxgrid;
1704 if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
1705 }
1706 }
1707 }
1708 }
1709
1710
1711 free(area1);
1712 free(area2);
1713
1714 free(x1);
1715 free(y1);
1716 free(z1);
1717 free(x2);
1718 free(y2);
1719 free(z2);
1720
1721 return nxgrid;
1722
1723 }/* create_xgrid_great_circle_ug */
1724
1725
1726 /*******************************************************************************
1727 Revise Sutherland-Hodgeman algorithm to find the vertices of the overlapping
1728 between any two grid boxes. It return the number of vertices for the exchange grid.
1729 Each edge of grid box is a part of great circle. All the points are cartesian
1730 coordinates. Here we are assuming each polygon is convex.
1731 RANGE_CHECK_CRITERIA is used to determine if the two grid boxes are possible to be
1732 overlap. The size should be between 0 and 0.5. The larger the range_check_criteria,
1733 the more expensive of the computatioin. When the value is close to 0,
1734 some small exchange grid might be lost. Suggest to use value 0.05 for C48.
1735 *******************************************************************************/
1736
1737 int clip_2dx2d_great_circle(const double x1_in[], const double y1_in[], const double z1_in[], int n1_in,
1738 const double x2_in[], const double y2_in[], const double z2_in [], int n2_in,
1739 double x_out[], double y_out[], double z_out[])
1740 {
1741 struct Node *grid1List=NULL;
1742 struct Node *grid2List=NULL;
1743 struct Node *intersectList=NULL;
1744 struct Node *polyList=NULL;
1745 struct Node *curList=NULL;
1746 struct Node *firstIntersect=NULL, *curIntersect=NULL;
1747 struct Node *temp1=NULL, *temp2=NULL, *temp=NULL;
1748
1749 int i1, i2, i1p, i2p, i2p2, npts1, npts2;
1750 int nintersect, n_out;
1751 int maxiter1, maxiter2, iter1, iter2;
1752 int found1, found2, curListNum;
1753 int has_inbound, inbound;
1754 double pt1[MV][3], pt2[MV][3];
1755 double *p1_0=NULL, *p1_1=NULL;
1756 double *p2_0=NULL, *p2_1=NULL, *p2_2=NULL;
1757 double intersect[3];
1758 double u1, u2;
1759 double min_x1, max_x1, min_y1, max_y1, min_z1, max_z1;
1760 double min_x2, max_x2, min_y2, max_y2, min_z2, max_z2;
1761
1762
1763 /* first check the min and max of (x1_in, y1_in, z1_in) with (x2_in, y2_in, z2_in) */
1764 min_x1 = minval_double(n1_in, x1_in);
1765 max_x2 = maxval_double(n2_in, x2_in);
1766 if(min_x1 >= max_x2+RANGE_CHECK_CRITERIA) return 0;
1767 max_x1 = maxval_double(n1_in, x1_in);
1768 min_x2 = minval_double(n2_in, x2_in);
1769 if(min_x2 >= max_x1+RANGE_CHECK_CRITERIA) return 0;
1770
1771 min_y1 = minval_double(n1_in, y1_in);
1772 max_y2 = maxval_double(n2_in, y2_in);
1773 if(min_y1 >= max_y2+RANGE_CHECK_CRITERIA) return 0;
1774 max_y1 = maxval_double(n1_in, y1_in);
1775 min_y2 = minval_double(n2_in, y2_in);
1776 if(min_y2 >= max_y1+RANGE_CHECK_CRITERIA) return 0;
1777
1778 min_z1 = minval_double(n1_in, z1_in);
1779 max_z2 = maxval_double(n2_in, z2_in);
1780 if(min_z1 >= max_z2+RANGE_CHECK_CRITERIA) return 0;
1781 max_z1 = maxval_double(n1_in, z1_in);
1782 min_z2 = minval_double(n2_in, z2_in);
1783 if(min_z2 >= max_z1+RANGE_CHECK_CRITERIA) return 0;
1784
1785 rewindList();
1786
1787 grid1List = getNext();
1788 grid2List = getNext();
1789 intersectList = getNext();
1790 polyList = getNext();
1791
1792 /* insert points into SubjList and ClipList */
1793 for(i1=0; i1<n1_in; i1++) addEnd(grid1List, x1_in[i1], y1_in[i1], z1_in[i1], 0, 0, 0, -1);
1794 for(i2=0; i2<n2_in; i2++) addEnd(grid2List, x2_in[i2], y2_in[i2], z2_in[i2], 0, 0, 0, -1);
1795 npts1 = length(grid1List);
1796 npts2 = length(grid2List);
1797
1798 n_out = 0;
1799 /* set the inside value */
1800 #ifdef debug_test_create_xgrid
1801 printf("\nNOTE from clip_2dx2d_great_circle: begin to set inside value grid1List\n");
1802 #endif
1803 /* first check number of points in grid1 is inside grid2 */
1804
1805 temp = grid1List;
1806 while(temp) {
1807 if(insidePolygon(temp, grid2List))
1808 temp->isInside = 1;
1809 else
1810 temp->isInside = 0;
1811 temp = getNextNode(temp);
1812 }
1813
1814 #ifdef debug_test_create_xgrid
1815 printf("\nNOTE from clip_2dx2d_great_circle: begin to set inside value of grid2List\n");
1816 #endif
1817 /* check if grid2List is inside grid1List */
1818 temp = grid2List;
1819
1820 while(temp) {
1821 if(insidePolygon(temp, grid1List))
1822 temp->isInside = 1;
1823 else
1824 temp->isInside = 0;
1825 temp = getNextNode(temp);
1826 }
1827
1828 /* make sure the grid box is clockwise */
1829
1830 /*make sure each polygon is convex, which is equivalent that the great_circle_area is positive */
1831 if( gridArea(grid1List) <= 0 )
1832 error_handler("create_xgrid.c(clip_2dx2d_great_circle): grid box 1 is not convex");
1833 if( gridArea(grid2List) <= 0 )
1834 error_handler("create_xgrid.c(clip_2dx2d_great_circle): grid box 2 is not convex");
1835
1836 #ifdef debug_test_create_xgrid
1837 printNode(grid1List, "grid1List");
1838 printNode(grid2List, "grid2List");
1839 #endif
1840
1841 /* get the coordinates from grid1List and grid2List.
1842 Please not npts1 might not equal n1_in, npts2 might not equal n2_in because of pole
1843 */
1844
1845 temp = grid1List;
1846 for(i1=0; i1<npts1; i1++) {
1847 getCoordinates(temp, pt1[i1]);
1848 temp = temp->Next;
1849 }
1850 temp = grid2List;
1851 for(i2=0; i2<npts2; i2++) {
1852 getCoordinates(temp, pt2[i2]);
1853 temp = temp->Next;
1854 }
1855
1856 firstIntersect=getNext();
1857 curIntersect = getNext();
1858
1859 #ifdef debug_test_create_xgrid
1860 printf("\n\n************************ Start line_intersect_2D_3D ******************************\n");
1861 #endif
1862 /* first find all the intersection points */
1863 nintersect = 0;
1864 for(i1=0; i1<npts1; i1++) {
1865 i1p = (i1+1)%npts1;
1866 p1_0 = pt1[i1];
1867 p1_1 = pt1[i1p];
1868 for(i2=0; i2<npts2; i2++) {
1869 i2p = (i2+1)%npts2;
1870 i2p2 = (i2+2)%npts2;
1871 p2_0 = pt2[i2];
1872 p2_1 = pt2[i2p];
1873 p2_2 = pt2[i2p2];
1874 #ifdef debug_test_create_xgrid
1875 printf("\n******************************************************************************\n");
1876 printf(" i1 = %d, i2 = %d \n", i1, i2);
1877 printf("********************************************************************************\n");
1878 #endif
1879 if( line_intersect_2D_3D(p1_0, p1_1, p2_0, p2_1, p2_2, intersect, &u1, &u2, &inbound) ) {
1880 /* from the value of u1, u2 and inbound, we can partially decide if a point is inside or outside of polygon */
1881
1882 /* add the intersection into intersetList, The intersection might already be in
1883 intersectList and will be taken care addIntersect
1884 */
1885 if(addIntersect(intersectList, intersect[0], intersect[1], intersect[2], 1, u1, u2, inbound, i1, i1p, i2, i2p)) {
1886 /* add the intersection into the grid1List */
1887
1888 if(u1 == 1) {
1889 insertIntersect(grid1List, intersect[0], intersect[1], intersect[2], 0.0, u2, inbound, p1_1[0], p1_1[1], p1_1[2]);
1890 }
1891 else
1892 insertIntersect(grid1List, intersect[0], intersect[1], intersect[2], u1, u2, inbound, p1_0[0], p1_0[1], p1_0[2]);
1893 /* when u1 == 0 or 1, need to adjust the vertice to intersect value for roundoff error */
1894 if(u1==1) {
1895 p1_1[0] = intersect[0];
1896 p1_1[1] = intersect[1];
1897 p1_1[2] = intersect[2];
1898 }
1899 else if(u1 == 0) {
1900 p1_0[0] = intersect[0];
1901 p1_0[1] = intersect[1];
1902 p1_0[2] = intersect[2];
1903 }
1904 /* add the intersection into the grid2List */
1905 if(u2==1)
1906 insertIntersect(grid2List, intersect[0], intersect[1], intersect[2], 0.0, u1, 0, p2_1[0], p2_1[1], p2_1[2]);
1907 else
1908 insertIntersect(grid2List, intersect[0], intersect[1], intersect[2], u2, u1, 0, p2_0[0], p2_0[1], p2_0[2]);
1909 /* when u2 == 0 or 1, need to adjust the vertice to intersect value for roundoff error */
1910 if(u2==1) {
1911 p2_1[0] = intersect[0];
1912 p2_1[1] = intersect[1];
1913 p2_1[2] = intersect[2];
1914 }
1915 else if(u2 == 0) {
1916 p2_0[0] = intersect[0];
1917 p2_0[1] = intersect[1];
1918 p2_0[2] = intersect[2];
1919 }
1920 }
1921 }
1922 }
1923 }
1924
1925 /* set inbound value for the points in intersectList that has inbound == 0,
1926 this will also set some inbound value of the points in grid1List
1927 */
1928
1929 /* get the first point in intersectList has inbound = 2, if not, set inbound value */
1930 has_inbound = 0;
1931 /* loop through intersectList to see if there is any has inbound=1 or 2 */
1932 temp = intersectList;
1933 nintersect = length(intersectList);
1934 if(nintersect > 1) {
1935 getFirstInbound(intersectList, firstIntersect);
1936 if(firstIntersect->initialized) {
1937 has_inbound = 1;
1938 }
1939 }
1940
1941 /* when has_inbound == 0, get the grid1List and grid2List */
1942 if( !has_inbound && nintersect > 1) {
1943 setInbound(intersectList, grid1List);
1944 getFirstInbound(intersectList, firstIntersect);
1945 if(firstIntersect->initialized) has_inbound = 1;
1946 }
1947
1948 /* if has_inbound = 1, find the overlapping */
1949 n_out = 0;
1950
1951 if(has_inbound) {
1952 maxiter1 = nintersect;
1953 #ifdef debug_test_create_xgrid
1954 printf("\nNOTE from clip_2dx2d_great_circle: number of intersect is %d\n", nintersect);
1955 printf("\n size of grid2List is %d, size of grid1List is %d\n", length(grid2List), length(grid1List));
1956 printNode(intersectList, "beginning intersection list");
1957 printNode(grid2List, "beginning clip list");
1958 printNode(grid1List, "beginning subj list");
1959 printf("\n************************ End line_intersect_2D_3D **********************************\n\n");
1960 #endif
1961 temp1 = getNode(grid1List, *firstIntersect);
1962 if( temp1 == NULL) {
1963 double lon[10], lat[10];
1964 int i;
1965 xyz2latlon(n1_in, x1_in, y1_in, z1_in, lon, lat);
1966 for(i=0; i< n1_in; i++) printf("lon1 = %g, lat1 = %g\n", lon[i]*R2D, lat[i]*R2D);
1967 printf("\n");
1968 xyz2latlon(n2_in, x2_in, y2_in, z2_in, lon, lat);
1969 for(i=0; i< n2_in; i++) printf("lon2 = %g, lat2 = %g\n", lon[i]*R2D, lat[i]*R2D);
1970 printf("\n");
1971
1972 error_handler("firstIntersect is not in the grid1List");
1973 }
1974 addNode(polyList, *firstIntersect);
1975 nintersect--;
1976 #ifdef debug_test_create_xgrid
1977 printNode(polyList, "polyList at stage 1");
1978 #endif
1979
1980 /* Loop over the grid1List and grid2List to find again the firstIntersect */
1981 curList = grid1List;
1982 curListNum = 0;
1983
1984 /* Loop through curList to find the next intersection, the loop will end
1985 when come back to firstIntersect
1986 */
1987 copyNode(curIntersect, *firstIntersect);
1988 iter1 = 0;
1989 found1 = 0;
1990
1991 while( iter1 < maxiter1 ) {
1992 #ifdef debug_test_create_xgrid
1993 printf("\n----------- At iteration = %d\n\n", iter1+1 );
1994 printNode(curIntersect, "curIntersect at the begining of iter1");
1995 #endif
1996 /* find the curIntersect in curList and get the next intersection points */
1997 temp1 = getNode(curList, *curIntersect);
1998 temp2 = temp1->Next;
1999 if( temp2 == NULL ) temp2 = curList;
2000
2001 maxiter2 = length(curList);
2002 found2 = 0;
2003 iter2 = 0;
2004 /* Loop until find the next intersection */
2005 while( iter2 < maxiter2 ) {
2006 int temp2IsIntersect;
2007
2008 temp2IsIntersect = 0;
2009 if( isIntersect( *temp2 ) ) { /* copy the point and switch to the grid2List */
2010 struct Node *temp3;
2011
2012 /* first check if temp2 is the firstIntersect */
2013 if( sameNode( *temp2, *firstIntersect) ) {
2014 found1 = 1;
2015 break;
2016 }
2017
2018 temp3 = temp2->Next;
2019 if( temp3 == NULL) temp3 = curList;
2020 if( temp3 == NULL) error_handler("creat_xgrid.c: temp3 can not be NULL");
2021 found2 = 1;
2022 /* if next node is inside or an intersection,
2023 need to keep on curList
2024 */
2025 temp2IsIntersect = 1;
2026 if( isIntersect(*temp3) || (temp3->isInside == 1) ) found2 = 0;
2027 }
2028 if(found2) {
2029 copyNode(curIntersect, *temp2);
2030 break;
2031 }
2032 else {
2033 addNode(polyList, *temp2);
2034 #ifdef debug_test_create_xgrid
2035 printNode(polyList, "polyList at stage 2");
2036 #endif
2037 if(temp2IsIntersect) {
2038 nintersect--;
2039 }
2040 }
2041 temp2 = temp2->Next;
2042 if( temp2 == NULL ) temp2 = curList;
2043 iter2 ++;
2044 }
2045 if(found1) break;
2046
2047 if( !found2 ) error_handler(" not found the next intersection ");
2048
2049 /* if find the first intersection, the poly found */
2050 if( sameNode( *curIntersect, *firstIntersect) ) {
2051 found1 = 1;
2052 break;
2053 }
2054
2055 /* add curIntersect to polyList and remove it from intersectList and curList */
2056 addNode(polyList, *curIntersect);
2057 #ifdef debug_test_create_xgrid
2058 printNode(polyList, "polyList at stage 3");
2059 #endif
2060 nintersect--;
2061
2062
2063 /* switch curList */
2064 if( curListNum == 0) {
2065 curList = grid2List;
2066 curListNum = 1;
2067 }
2068 else {
2069 curList = grid1List;
2070 curListNum = 0;
2071 }
2072 iter1++;
2073 }
2074 if(!found1) error_handler("not return back to the first intersection");
2075
2076 /* currently we are only clipping convex polygon to convex polygon */
2077 if( nintersect > 0) error_handler("After clipping, nintersect should be 0");
2078
2079 /* copy the polygon to x_out, y_out, z_out */
2080 temp1 = polyList;
2081 while (temp1 != NULL) {
2082 getCoordinate(*temp1, x_out+n_out, y_out+n_out, z_out+n_out);
2083 temp1 = temp1->Next;
2084 n_out++;
2085 }
2086
2087 /* if(n_out < 3) error_handler(" The clipped region has < 3 vertices"); */
2088 if( n_out < 3) n_out = 0;
2089 #ifdef debug_test_create_xgrid
2090 printNode(polyList, "polyList after clipping");
2091 #endif
2092 }
2093
2094 /* check if grid1 is inside grid2 */
2095 if(n_out==0){
2096 /* first check number of points in grid1 is inside grid2 */
2097 int n, n1in2;
2098 /* One possible is that grid1List is inside grid2List */
2099 #ifdef debug_test_create_xgrid
2100 printf("\nNOTE from clip_2dx2d_great_circle: check if grid1 is inside grid2\n");
2101 #endif
2102 n1in2 = 0;
2103 temp = grid1List;
2104 while(temp) {
2105 if(temp->intersect != 1) {
2106 #ifdef debug_test_create_xgrid
2107 printf("grid1->isInside = %d\n", temp->isInside);
2108 #endif
2109 if( temp->isInside == 1) n1in2++;
2110 }
2111 temp = getNextNode(temp);
2112 }
2113 if(npts1==n1in2) { /* grid1 is inside grid2 */
2114 n_out = npts1;
2115 n = 0;
2116 temp = grid1List;
2117 while( temp ) {
2118 getCoordinate(*temp, &x_out[n], &y_out[n], &z_out[n]);
2119 n++;
2120 temp = getNextNode(temp);
2121 }
2122 }
2123 if(n_out>0) return n_out;
2124 }
2125
2126 /* check if grid2List is inside grid1List */
2127 if(n_out ==0){
2128 int n, n2in1;
2129 #ifdef debug_test_create_xgrid
2130 printf("\nNOTE from clip_2dx2d_great_circle: check if grid2 is inside grid1\n");
2131 #endif
2132
2133 temp = grid2List;
2134 n2in1 = 0;
2135 while(temp) {
2136 if(temp->intersect != 1) {
2137 #ifdef debug_test_create_xgrid
2138 printf("grid2->isInside = %d\n", temp->isInside);
2139 #endif
2140 if( temp->isInside == 1) n2in1++;
2141 }
2142 temp = getNextNode(temp);
2143 }
2144
2145 if(npts2==n2in1) { /* grid2 is inside grid1 */
2146 n_out = npts2;
2147 n = 0;
2148 temp = grid2List;
2149 while( temp ) {
2150 getCoordinate(*temp, &x_out[n], &y_out[n], &z_out[n]);
2151 n++;
2152 temp = getNextNode(temp);
2153 }
2154
2155 }
2156 }
2157
2158
2159 return n_out;
2160 }
2161
2162
2163 /* Intersects between the line a and the seqment s
2164 where both line and segment are great circle lines on the sphere represented by
2165 3D cartesian points.
2166 [sin sout] are the ends of a line segment
2167 returns true if the lines could be intersected, false otherwise.
2168 inbound means the direction of (a1,a2) go inside or outside of (q1,q2,q3)
2169 */
2170
2171 int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double *q3,
2172 double *intersect, double *u_a, double *u_q, int *inbound){
2173
2174 /* Do this intersection by reprsenting the line a1 to a2 as a plane through the
2175 two line points and the origin of the sphere (0,0,0). This is the
2176 definition of a great circle arc.
2177 */
2178 double plane[9];
2179 double plane_p[2];
2180 double u;
2181 double p1[3], v1[3], v2[3];
2182 double c1[3], c2[3], c3[3];
2183 double coincident, sense, norm;
2184 int i;
2185 int is_inter1, is_inter2;
2186
2187 *inbound = 0;
2188
2189 /* first check if any vertices are the same */
2190 if(samePoint(a1[0], a1[1], a1[2], q1[0], q1[1], q1[2])) {
2191 *u_a = 0;
2192 *u_q = 0;
2193 intersect[0] = a1[0];
2194 intersect[1] = a1[1];
2195 intersect[2] = a1[2];
2196 #ifdef debug_test_create_xgrid
2197 printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f, u_q=%19.15f, inbound=%d\n", *u_a, *u_q, *inbound);
2198 #endif
2199 return 1;
2200 }
2201 else if (samePoint(a1[0], a1[1], a1[2], q2[0], q2[1], q2[2])) {
2202 *u_a = 0;
2203 *u_q = 1;
2204 intersect[0] = a1[0];
2205 intersect[1] = a1[1];
2206 intersect[2] = a1[2];
2207 #ifdef debug_test_create_xgrid
2208 printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f, u_q=%19.15f, inbound=%d\n", *u_a, *u_q, *inbound);
2209 #endif
2210 return 1;
2211 }
2212 else if(samePoint(a2[0], a2[1], a2[2], q1[0], q1[1], q1[2])) {
2213 #ifdef debug_test_create_xgrid
2214 printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f, u_q=%19.15f, inbound=%d\n", *u_a, *u_q, *inbound);
2215 #endif
2216 *u_a = 1;
2217 *u_q = 0;
2218 intersect[0] = a2[0];
2219 intersect[1] = a2[1];
2220 intersect[2] = a2[2];
2221 return 1;
2222 }
2223 else if (samePoint(a2[0], a2[1], a2[2], q2[0], q2[1], q2[2])) {
2224 #ifdef debug_test_create_xgrid
2225 printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f, u_q=%19.15f, inbound=%d\n", *u_a, *u_q, *inbound);
2226 #endif
2227 *u_a = 1;
2228 *u_q = 1;
2229 intersect[0] = a2[0];
2230 intersect[1] = a2[1];
2231 intersect[2] = a2[2];
2232 return 1;
2233 }
2234
2235
2236 /* Load points defining plane into variable (these are supposed to be in counterclockwise order) */
2237 plane[0]=q1[0];
2238 plane[1]=q1[1];
2239 plane[2]=q1[2];
2240 plane[3]=q2[0];
2241 plane[4]=q2[1];
2242 plane[5]=q2[2];
2243 plane[6]=0.0;
2244 plane[7]=0.0;
2245 plane[8]=0.0;
2246
2247 /* Intersect the segment with the plane */
2248 is_inter1 = intersect_tri_with_line(plane, a1, a2, plane_p, u_a);
2249
2250 if(!is_inter1)
2251 return 0;
2252
2253 if(fabs(*u_a) < EPSLN8) *u_a = 0;
2254 if(fabs(*u_a-1) < EPSLN8) *u_a = 1;
2255
2256
2257 #ifdef debug_test_create_xgrid
2258 printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f\n", *u_a);
2259 #endif
2260
2261
2262 if( (*u_a < 0) || (*u_a > 1) ) return 0;
2263
2264 /* Load points defining plane into variable (these are supposed to be in counterclockwise order) */
2265 plane[0]=a1[0];
2266 plane[1]=a1[1];
2267 plane[2]=a1[2];
2268 plane[3]=a2[0];
2269 plane[4]=a2[1];
2270 plane[5]=a2[2];
2271 plane[6]=0.0;
2272 plane[7]=0.0;
2273 plane[8]=0.0;
2274
2275 /* Intersect the segment with the plane */
2276 is_inter2 = intersect_tri_with_line(plane, q1, q2, plane_p, u_q);
2277
2278 if(!is_inter2)
2279 return 0;
2280
2281 if(fabs(*u_q) < EPSLN8) *u_q = 0;
2282 if(fabs(*u_q-1) < EPSLN8) *u_q = 1;
2283 #ifdef debug_test_create_xgrid
2284 printf("\nNOTE from line_intersect_2D_3D: u_q = %19.15f\n", *u_q);
2285 #endif
2286
2287
2288 if( (*u_q < 0) || (*u_q > 1) ) return 0;
2289
2290 u =*u_a;
2291
2292 /* The two planes are coincidental */
2293 vect_cross(a1, a2, c1);
2294 vect_cross(q1, q2, c2);
2295 vect_cross(c1, c2, c3);
2296 coincident = metric(c3);
2297
2298 if(fabs(coincident) < EPSLN30) return 0;
2299
2300 /* Calculate point of intersection */
2301 intersect[0]=a1[0] + u*(a2[0]-a1[0]);
2302 intersect[1]=a1[1] + u*(a2[1]-a1[1]);
2303 intersect[2]=a1[2] + u*(a2[2]-a1[2]);
2304
2305 norm = metric( intersect );
2306 for(i = 0; i < 3; i ++) intersect[i] /= norm;
2307
2308 /* when u_q =0 or u_q =1, the following could not decide the inbound value */
2309 if(*u_q != 0 && *u_q != 1){
2310
2311 p1[0] = a2[0]-a1[0];
2312 p1[1] = a2[1]-a1[1];
2313 p1[2] = a2[2]-a1[2];
2314 v1[0] = q2[0]-q1[0];
2315 v1[1] = q2[1]-q1[1];
2316 v1[2] = q2[2]-q1[2];
2317 v2[0] = q3[0]-q2[0];
2318 v2[1] = q3[1]-q2[1];
2319 v2[2] = q3[2]-q2[2];
2320
2321 vect_cross(v1, v2, c1);
2322 vect_cross(v1, p1, c2);
2323
2324 sense = dot(c1, c2);
2325 *inbound = 1;
2326 if(sense > 0) *inbound = 2; /* v1 going into v2 in CCW sense */
2327 }
2328 #ifdef debug_test_create_xgrid
2329 printf("\nNOTE from line_intersect_2D_3D: inbound=%d\n", *inbound);
2330 #endif
2331
2332 return 1;
2333 }
2334
2335
2336 /*------------------------------------------------------------------------------
2337 double poly_ctrlat(const double x[], const double y[], int n)
2338 This routine is used to calculate the latitude of the centroid
2339 ---------------------------------------------------------------------------*/
2340
2341 double poly_ctrlat(const double x[], const double y[], int n)
2342 {
2343 double ctrlat = 0.0;
2344 int i;
2345
2346 for (i=0;i<n;i++) {
2347 int ip = (i+1) % n;
2348 double dx = (x[ip]-x[i]);
2349 double dy, avg_y, hdy;
2350 double lat1, lat2;
2351 lat1 = y[ip];
2352 lat2 = y[i];
2353 dy = lat2 - lat1;
2354 hdy = dy*0.5;
2355 avg_y = (lat1+lat2)*0.5;
2356 if (dx==0.0) continue;
2357 if(dx > M_PI) dx = dx - 2.0*M_PI;
2358 if(dx < -M_PI) dx = dx + 2.0*M_PI;
2359
2360 if ( fabs(hdy)< SMALL_VALUE ) /* cheap area calculation along latitude */
2361 ctrlat -= dx*(2*cos(avg_y) + lat2*sin(avg_y) - cos(lat1) );
2362 else
2363 ctrlat -= dx*( (sin(hdy)/hdy)*(2*cos(avg_y) + lat2*sin(avg_y)) - cos(lat1) );
2364 }
2365 return (ctrlat*RADIUS*RADIUS);
2366 } /* poly_ctrlat */
2367
2368 /*------------------------------------------------------------------------------
2369 double poly_ctrlon(const double x[], const double y[], int n, double clon)
2370 This routine is used to calculate the lontitude of the centroid.
2371 ---------------------------------------------------------------------------*/
2372 double poly_ctrlon(const double x[], const double y[], int n, double clon)
2373 {
2374 double ctrlon = 0.0;
2375 int i;
2376
2377 for (i=0;i<n;i++) {
2378 int ip = (i+1) % n;
2379 double phi1, phi2, dphi, lat1, lat2, dphi1, dphi2;
2380 double f1, f2, fac, fint;
2381 phi1 = x[ip];
2382 phi2 = x[i];
2383 lat1 = y[ip];
2384 lat2 = y[i];
2385 dphi = phi1 - phi2;
2386
2387 if (dphi==0.0) continue;
2388
2389 f1 = 0.5*(cos(lat1)*sin(lat1)+lat1);
2390 f2 = 0.5*(cos(lat2)*sin(lat2)+lat2);
2391
2392 /* this will make sure longitude of centroid is at
2393 the same interval as the center of any grid */
2394 if(dphi > M_PI) dphi = dphi - 2.0*M_PI;
2395 if(dphi < -M_PI) dphi = dphi + 2.0*M_PI;
2396 dphi1 = phi1 - clon;
2397 if( dphi1 > M_PI) dphi1 -= 2.0*M_PI;
2398 if( dphi1 <-M_PI) dphi1 += 2.0*M_PI;
2399 dphi2 = phi2 -clon;
2400 if( dphi2 > M_PI) dphi2 -= 2.0*M_PI;
2401 if( dphi2 <-M_PI) dphi2 += 2.0*M_PI;
2402
2403 if(fabs(dphi2 -dphi1) < M_PI) {
2404 ctrlon -= dphi * (dphi1*f1+dphi2*f2)/2.0;
2405 }
2406 else {
2407 if(dphi1 > 0.0)
2408 fac = M_PI;
2409 else
2410 fac = -M_PI;
2411 fint = f1 + (f2-f1)*(fac-dphi1)/fabs(dphi);
2412 ctrlon -= 0.5*dphi1*(dphi1-fac)*f1 - 0.5*dphi2*(dphi2+fac)*f2
2413 + 0.5*fac*(dphi1+dphi2)*fint;
2414 }
2415
2416 }
2417 return (ctrlon*RADIUS*RADIUS);
2418 } /* poly_ctrlon */
2419
2420 /* -----------------------------------------------------------------------------
2421 double box_ctrlat(double ll_lon, double ll_lat, double ur_lon, double ur_lat)
2422 This routine is used to calculate the latitude of the centroid.
2423 ---------------------------------------------------------------------------*/
2424 double box_ctrlat(double ll_lon, double ll_lat, double ur_lon, double ur_lat)
2425 {
2426 double dphi = ur_lon-ll_lon;
2427 double ctrlat;
2428
2429 if(dphi > M_PI) dphi = dphi - 2.0*M_PI;
2430 if(dphi < -M_PI) dphi = dphi + 2.0*M_PI;
2431 ctrlat = dphi*(cos(ur_lat) + ur_lat*sin(ur_lat)-(cos(ll_lat) + ll_lat*sin(ll_lat)));
2432 return (ctrlat*RADIUS*RADIUS);
2433 } /* box_ctrlat */
2434
2435 /*------------------------------------------------------------------------------
2436 double box_ctrlon(double ll_lon, double ll_lat, double ur_lon, double ur_lat, double clon)
2437 This routine is used to calculate the lontitude of the centroid
2438 ----------------------------------------------------------------------------*/
2439 double box_ctrlon(double ll_lon, double ll_lat, double ur_lon, double ur_lat, double clon)
2440 {
2441 double phi1, phi2, dphi, lat1, lat2, dphi1, dphi2;
2442 double f1, f2, fac, fint;
2443 double ctrlon = 0.0;
2444 int i;
2445 for( i =0; i<2; i++) {
2446 if(i == 0) {
2447 phi1 = ur_lon;
2448 phi2 = ll_lon;
2449 lat1 = lat2 = ll_lat;
2450 }
2451 else {
2452 phi1 = ll_lon;
2453 phi2 = ur_lon;
2454 lat1 = lat2 = ur_lat;
2455 }
2456 dphi = phi1 - phi2;
2457 f1 = 0.5*(cos(lat1)*sin(lat1)+lat1);
2458 f2 = 0.5*(cos(lat2)*sin(lat2)+lat2);
2459
2460 if(dphi > M_PI) dphi = dphi - 2.0*M_PI;
2461 if(dphi < -M_PI) dphi = dphi + 2.0*M_PI;
2462 /* make sure the center is in the same grid box. */
2463 dphi1 = phi1 - clon;
2464 if( dphi1 > M_PI) dphi1 -= 2.0*M_PI;
2465 if( dphi1 <-M_PI) dphi1 += 2.0*M_PI;
2466 dphi2 = phi2 -clon;
2467 if( dphi2 > M_PI) dphi2 -= 2.0*M_PI;
2468 if( dphi2 <-M_PI) dphi2 += 2.0*M_PI;
2469
2470 if(fabs(dphi2 -dphi1) < M_PI) {
2471 ctrlon -= dphi * (dphi1*f1+dphi2*f2)/2.0;
2472 }
2473 else {
2474 if(dphi1 > 0.0)
2475 fac = M_PI;
2476 else
2477 fac = -M_PI;
2478 fint = f1 + (f2-f1)*(fac-dphi1)/fabs(dphi);
2479 ctrlon -= 0.5*dphi1*(dphi1-fac)*f1 - 0.5*dphi2*(dphi2+fac)*f2
2480 + 0.5*fac*(dphi1+dphi2)*fint;
2481 }
2482 }
2483 return (ctrlon*RADIUS*RADIUS);
2484 } /* box_ctrlon */
2485
2486 /*******************************************************************************
2487 double grid_box_radius(double *x, double *y, double *z, int n);
2488 Find the radius of the grid box, the radius is defined the
2489 maximum distance between any two vertices
2490 *******************************************************************************/
2491 double grid_box_radius(const double *x, const double *y, const double *z, int n)
2492 {
2493 double radius;
2494 int i, j;
2495
2496 radius = 0;
2497 for(i=0; i<n-1; i++) {
2498 for(j=i+1; j<n; j++) {
2499 radius = max(radius, pow(x[i]-x[j],2.)+pow(y[i]-y[j],2.)+pow(z[i]-z[j],2.));
2500 }
2501 }
2502
2503 radius = sqrt(radius);
2504
2505 return (radius);
2506
2507 } /* grid_box_radius */
2508
2509 /*******************************************************************************
2510 double dist_between_boxes(const double *x1, const double *y1, const double *z1, int n1,
2511 const double *x2, const double *y2, const double *z2, int n2);
2512 Find the distance between any two grid boxes. The distance is defined by the maximum
2513 distance between any vertices of these two box
2514 *******************************************************************************/
2515 double dist_between_boxes(const double *x1, const double *y1, const double *z1, int n1,
2516 const double *x2, const double *y2, const double *z2, int n2)
2517 {
2518 double dist;
2519 int i, j;
2520
2521 dist = 0.0;
2522 for(i=0; i<n1; i++) {
2523 for(j=0; j<n2; j++) {
2524 dist = max(dist, pow(x1[i]-x2[j],2.)+pow(y1[i]-y2[j],2.)+pow(z1[i]-z2[j],2.));
2525 }
2526 }
2527
2528 dist = sqrt(dist);
2529 return (dist);
2530
2531 } /* dist_between_boxes */
2532
2533 /*******************************************************************************
2534 int inside_edge(double x0, double y0, double x1, double y1, double x, double y)
2535 determine a point(x,y) is inside or outside a given edge with vertex,
2536 (x0,y0) and (x1,y1). return 1 if inside and 0 if outside. <y1-y0, -(x1-x0)> is
2537 the outward edge normal from vertex <x0,y0> to <x1,y1>. <x-x0,y-y0> is the vector
2538 from <x0,y0> to <x,y>.
2539 if Inner produce <x-x0,y-y0>*<y1-y0, -(x1-x0)> > 0, outside, otherwise inside.
2540 inner product value = 0 also treate as inside.
2541 *******************************************************************************/
2542 int inside_edge(double x0, double y0, double x1, double y1, double x, double y)
2543 {
2544 const double SMALL = 1.e-12;
2545 double product;
2546
2547 product = ( x-x0 )*(y1-y0) + (x0-x1)*(y-y0);
2548 return (product<=SMALL) ? 1:0;
2549
2550 } /* inside_edge */
2551
2552
2553 /* The following is a test program to test subroutines in create_xgrid.c */
2554
2555 #ifdef test_create_xgrid
2556
2557 #include "create_xgrid.h"
2558 #include <math.h>
2559
2560 #define D2R (M_PI/180)
2561 #define R2D (180/M_PI)
2562 #define MAXPOINT 1000
2563
2564 int main(int argc, char* argv[])
2565 {
2566
2567 double lon1_in[MAXPOINT], lat1_in[MAXPOINT];
2568 double lon2_in[MAXPOINT], lat2_in[MAXPOINT];
2569 double x1_in[MAXPOINT], y1_in[MAXPOINT], z1_in[MAXPOINT];
2570 double x2_in[MAXPOINT], y2_in[MAXPOINT], z2_in[MAXPOINT];
2571 double lon_out[20], lat_out[20];
2572 double x_out[20], y_out[20], z_out[20];
2573 int n1_in, n2_in, n_out, i, j;
2574 int nlon1=0, nlat1=0, nlon2=0, nlat2=0;
2575 int n;
2576 int ntest = 11;
2577
2578
2579 for(n=11; n<=ntest; n++) {
2580
2581 switch (n) {
2582 case 1:
2583 /****************************************************************
2584
2585 test clip_2dx2d_great_cirle case 1:
2586 box 1: (20,10), (20,12), (22,12), (22,10)
2587 box 2: (21,11), (21,14), (24,14), (24,11)
2588 out : (21, 12.0018), (22, 12), (22, 11.0033), (21, 11)
2589
2590 ****************************************************************/
2591 n1_in = 4; n2_in = 4;
2592 /* first a simple lat-lon grid box to clip another lat-lon grid box */
2593 lon1_in[0] = 20; lat1_in[0] = 10;
2594 lon1_in[1] = 20; lat1_in[1] = 12;
2595 lon1_in[2] = 22; lat1_in[2] = 12;
2596 lon1_in[3] = 22; lat1_in[3] = 10;
2597 lon2_in[0] = 21; lat2_in[0] = 11;
2598 lon2_in[1] = 21; lat2_in[1] = 14;
2599 lon2_in[2] = 24; lat2_in[2] = 14;
2600 lon2_in[3] = 24; lat2_in[3] = 11;
2601 break;
2602
2603 case 2:
2604 /****************************************************************
2605
2606 test clip_2dx2d_great_cirle case 2: two identical box
2607 box 1: (20,10), (20,12), (22,12), (22,10)
2608 box 2: (20,10), (20,12), (22,12), (22,10)
2609 out : (20,10), (20,12), (22,12), (22,10)
2610
2611 ****************************************************************/
2612 lon1_in[0] = 20; lat1_in[0] = 10;
2613 lon1_in[1] = 20; lat1_in[1] = 12;
2614 lon1_in[2] = 22; lat1_in[2] = 12;
2615 lon1_in[3] = 22; lat1_in[3] = 10;
2616
2617 for(i=0; i<n2_in; i++) {
2618 lon2_in[i] = lon1_in[i];
2619 lat2_in[i] = lat1_in[i];
2620 }
2621 break;
2622
2623 case 3:
2624 /****************************************************************
2625
2626 test clip_2dx2d_great_cirle case 3: one cubic sphere grid close to the pole with lat-lon grid.
2627 box 1: (251.7, 88.98), (148.3, 88.98), (57.81, 88.72), (342.2, 88.72)
2628 box 2: (150, 88), (150, 90), (152.5, 90), (152.5, 88)
2629 out : (152.5, 89.0642), (150, 89.0165), (0, 90)
2630
2631 ****************************************************************/
2632 n1_in = 4; n2_in = 4;
2633 /* first a simple lat-lon grid box to clip another lat-lon grid box */
2634 lon1_in[0] = 251.7; lat1_in[0] = 88.98;
2635 lon1_in[1] = 148.3; lat1_in[1] = 88.98;
2636 lon1_in[2] = 57.81; lat1_in[2] = 88.72;
2637 lon1_in[3] = 342.2; lat1_in[3] = 88.72;
2638
2639 lon2_in[0] = 150; lat2_in[0] = 88;
2640 lon2_in[1] = 150; lat2_in[1] = 90;
2641 lon2_in[2] = 152.5; lat2_in[2] = 90;
2642 lon2_in[3] = 152.5; lat2_in[3] = 88;
2643 /*
2644 for(i=0; i<4; i++) {
2645 lon2_in[i] = lon1_in[i];
2646 lat2_in[i] = lat1_in[i];
2647 }
2648 */
2649 break;
2650
2651 case 4:
2652 /****************************************************************
2653
2654 test clip_2dx2d_great_cirle case 4: One box contains the pole
2655 box 1: (-160, 88.5354), (152.011, 87.8123) , (102.985, 88.4008), (20, 89.8047)
2656 box 2: (145,88), (145,90), (150,90), (150,88)
2657 out : (145.916, 88.0011), (145, 88.0249), (0, 90), (150, 88)
2658
2659 ****************************************************************/
2660 n1_in = 4; n2_in = 4;
2661 /* first a simple lat-lon grid box to clip another lat-lon grid box */
2662
2663 lon1_in[0] = -160; lat1_in[0] = 88.5354;
2664 lon1_in[1] = 152.011; lat1_in[1] = 87.8123;
2665 lon1_in[2] = 102.985; lat1_in[2] = 88.4008;
2666 lon1_in[3] = 20; lat1_in[3] = 89.8047;
2667
2668 lon2_in[0] = 145; lat2_in[0] = 88;
2669 lon2_in[1] = 145; lat2_in[1] = 90;
2670 lon2_in[2] = 150; lat2_in[2] = 90;
2671 lon2_in[3] = 150; lat2_in[3] = 88;
2672 break;
2673
2674 case 5:
2675 /****************************************************************
2676
2677 test clip_2dx2d_great_cirle case 5: One tripolar grid around the pole with lat-lon grid.
2678 box 1: (-202.6, 87.95), (-280, 89.56), (-100, 90), (-190, 88)
2679 box 2: (21,11), (21,14), (24,14), (24,11)
2680 out : (150, 88.7006), (145, 88.9507), (0, 90)
2681
2682 ****************************************************************/
2683 n1_in = 4; n2_in = 4;
2684 /* first a simple lat-lon grid box to clip another lat-lon grid box */
2685
2686 lon1_in[0] = -202.6; lat1_in[0] = 87.95;
2687 lon1_in[1] = -280.; lat1_in[1] = 89.56;
2688 lon1_in[2] = -100.0; lat1_in[2] = 90;
2689 lon1_in[3] = -190.; lat1_in[3] = 88;
2690
2691 lon2_in[0] = 145; lat2_in[0] = 88;
2692 lon2_in[1] = 145; lat2_in[1] = 90;
2693 lon2_in[2] = 150; lat2_in[2] = 90;
2694 lon2_in[3] = 150; lat2_in[3] = 88;
2695 break;
2696
2697 case 6:
2698 /****************************************************************
2699
2700 test clip_2dx2d_great_cirle case 6: One cubic sphere grid arounc the pole with one tripolar grid box
2701 around the pole.
2702 box 1: (-160, 88.5354), (152.011, 87.8123) , (102.985, 88.4008), (20, 89.8047)
2703 box 2: (-202.6, 87.95), (-280, 89.56), (-100, 90), (-190, 88)
2704 out : (170, 88.309), (157.082, 88.0005), (83.714, 89.559), (80, 89.6094), (0, 90), (200, 88.5354)
2705
2706
2707 ****************************************************************/
2708 n1_in = 4; n2_in = 4;
2709 /* first a simple lat-lon grid box to clip another lat-lon grid box */
2710
2711 lon1_in[0] = -160; lat1_in[0] = 88.5354;
2712 lon1_in[1] = 152.011; lat1_in[1] = 87.8123;
2713 lon1_in[2] = 102.985; lat1_in[2] = 88.4008;
2714 lon1_in[3] = 20; lat1_in[3] = 89.8047;
2715
2716 lon2_in[0] = -202.6; lat2_in[0] = 87.95;
2717 lon2_in[1] = -280.; lat2_in[1] = 89.56;
2718 lon2_in[2] = -100.0; lat2_in[2] = 90;
2719 lon2_in[3] = -190.; lat2_in[3] = 88;
2720 break;
2721
2722 case 7:
2723 /****************************************************************
2724
2725 test clip_2dx2d_great_cirle case 7: One small grid box inside a big grid box.
2726 box 1: (20,10), (20,12), (22,12), (22,10)
2727 box 2: (18,8), (18,14), (24,14), (24,8)
2728 out : (20,10), (20,12), (22,12), (22,10)
2729
2730 ****************************************************************/
2731 n1_in = 4; n2_in = 4;
2732 /* first a simple lat-lon grid box to clip another lat-lon grid box */
2733 lon1_in[0] = 20; lat1_in[0] = 10;
2734 lon1_in[1] = 20; lat1_in[1] = 12;
2735 lon1_in[2] = 22; lat1_in[2] = 12;
2736 lon1_in[3] = 22; lat1_in[3] = 10;
2737 lon2_in[0] = 18; lat2_in[0] = 8;
2738 lon2_in[1] = 18; lat2_in[1] = 14;
2739 lon2_in[2] = 24; lat2_in[2] = 14;
2740 lon2_in[3] = 24; lat2_in[3] = 8;
2741 break;
2742
2743 case 8:
2744 /****************************************************************
2745
2746 test clip_2dx2d_great_cirle case 8: Cubic sphere grid at tile = 1, point (i=25,j=1)
2747 with N45 at (i=141,j=23)
2748 box 1:
2749 box 2:
2750 out : None
2751
2752 ****************************************************************/
2753 n1_in = 4; n2_in = 4;
2754 /* first a simple lat-lo
2755 n grid box to clip another lat-lon grid box */
2756 lon1_in[0] = 350.0; lat1_in[0] = -45;
2757 lon1_in[1] = 350.0; lat1_in[1] = -43.43;
2758 lon1_in[2] = 352.1; lat1_in[2] = -43.41;
2759 lon1_in[3] = 352.1; lat1_in[3] = -44.98;
2760 lon2_in[0] = 350.0; lat2_in[0] = -46;
2761 lon2_in[1] = 350.0; lat2_in[1] = -44;
2762 lon2_in[2] = 352.5; lat2_in[2] = -44;
2763 lon2_in[3] = 352.5; lat2_in[3] = -46;
2764 break;
2765
2766 case 9:
2767 /****************************************************************
2768
2769 test clip_2dx2d_great_cirle case 9: Cubic sphere grid at tile = 1, point (i=1,j=1)
2770 with N45 at (i=51,j=61)
2771 box 1:
2772 box 2:
2773 out : None
2774
2775 ****************************************************************/
2776 n1_in = 4; n2_in = 4;
2777
2778 lon1_in[0] = 305.0; lat1_in[0] = -35.26;
2779 lon1_in[1] = 305.0; lat1_in[1] = -33.80;
2780 lon1_in[2] = 306.6; lat1_in[2] = -34.51;
2781 lon1_in[3] = 306.6; lat1_in[3] = -35.99;
2782 lon2_in[0] = 125; lat2_in[0] = 32;
2783 lon2_in[1] = 125; lat2_in[1] = 34;
2784 lon2_in[2] = 127.5; lat2_in[2] = 34;
2785 lon2_in[3] = 127.5; lat2_in[3] = 32;
2786 break;
2787
2788 case 10:
2789 /****************************************************************
2790
2791 test clip_2dx2d_great_cirle case 10: Cubic sphere grid at tile = 3, point (i=24,j=1)
2792 with N45 at (i=51,j=46)
2793 box 1:
2794 box 2:
2795 out : None
2796
2797 ****************************************************************/
2798 n1_in = 4; n2_in = 4;
2799
2800 lon1_in[0] = 125.0; lat1_in[0] = 1.46935;
2801 lon1_in[1] = 126.573; lat1_in[1] = 1.5091;
2802 lon1_in[2] = 126.573; lat1_in[2] = 0;
2803 lon1_in[3] = 125.0; lat1_in[3] = 0;
2804 lon2_in[0] = 125; lat2_in[0] = 0;
2805 lon2_in[1] = 125; lat2_in[1] = 2;
2806 lon2_in[2] = 127.5; lat2_in[2] = 2;
2807 lon2_in[3] = 127.5; lat2_in[3] = 0;
2808 break;
2809
2810 case 11:
2811 /****************************************************************
2812
2813 test clip_2dx2d_great_cirle case 10: Cubic sphere grid at tile = 3, point (i=24,j=1)
2814 with N45 at (i=51,j=46)
2815 box 1:
2816 box 2:
2817 out :
2818
2819 ****************************************************************/
2820 nlon1 = 1;
2821 nlat1 = 1;
2822 nlon2 = 1;
2823 nlat2 = 1;
2824 n1_in = (nlon1+1)*(nlat1+1);
2825 n2_in = (nlon2+1)*(nlat2+1);
2826
2827 lon1_in[0] = 350.0; lat1_in[0] = 90.00;
2828 lon1_in[1] = 170.0; lat1_in[1] = 87.92;
2829 lon1_in[2] = 260.0; lat1_in[2] = 87.92;
2830 lon1_in[3] = 215.0; lat1_in[3] = 87.06;
2831
2832 /* lon1_in[0] = 35.0; lat1_in[0] = 87.06; */
2833 /* lon1_in[1] = 80.0; lat1_in[1] = 87.92; */
2834 /* lon1_in[2] = 125.0; lat1_in[2] = 87.06; */
2835 /* lon1_in[3] = 350.0; lat1_in[3] = 87.92; */
2836 /* lon1_in[4] = 350.0; lat1_in[4] = 90.00; */
2837 /* lon1_in[5] = 170.0; lat1_in[5] = 87.92; */
2838 /* lon1_in[6] = 305.0; lat1_in[6] = 87.06; */
2839 /* lon1_in[7] = 260.0; lat1_in[7] = 87.92; */
2840 /* lon1_in[8] = 215.0; lat1_in[8] = 87.06; */
2841
2842 lon2_in[0] = 167.5; lat2_in[0] = 88;
2843 lon2_in[1] = 170; lat2_in[1] = 88;
2844 lon2_in[2] = 167.5; lat2_in[2] = 90;
2845 lon2_in[3] = 170; lat2_in[3] = 90;
2846
2847 /* nlon1 = 3; */
2848 /* nlat1 = 2; */
2849 /* nlon2 = 1; */
2850 /* nlat2 = 1; */
2851 /* n1_in = (nlon1+1)*(nlat1+1); */
2852 /* n2_in = (nlon2+1)*(nlat2+1); */
2853
2854 /* lon1_in[0] = 35.00; lat1_in[0] = -59.90; */
2855 /* lon1_in[1] = 37.64; lat1_in[1] = -58.69; */
2856 /* lon1_in[2] = 40.07; lat1_in[2] = -57.44; */
2857 /* lon1_in[3] = 42.32; lat1_in[3] = -56.15; */
2858 /* lon1_in[4] = 32.36; lat1_in[4] = -58.69; */
2859 /* lon1_in[5] = 35.00; lat1_in[5] = -57.56; */
2860 /* lon1_in[6] = 37.45; lat1_in[6] = -56.39; */
2861 /* lon1_in[7] = 39.74; lat1_in[7] = -55.18; */
2862 /* lon1_in[8] = 29.93; lat1_in[8] = -57.44; */
2863 /* lon1_in[9] = 32.55; lat1_in[9] = -56.39; */
2864 /* lon1_in[10] = 35.00; lat1_in[10] = -55.29; */
2865 /* lon1_in[11] = 37.30; lat1_in[11] = -54.16; */
2866 /* lon2_in[0] = 35; lat2_in[0] = -58; */
2867 /* lon2_in[1] = 37.5; lat2_in[1] = -58; */
2868 /* lon2_in[2] = 35; lat2_in[2] = -56; */
2869 /* lon2_in[3] = 37.5; lat2_in[3] = -56; */
2870
2871 /* nlon1 = 1; */
2872 /* nlat1 = 1; */
2873 /* nlon2 = 1; */
2874 /* nlat2 = 1; */
2875 /* n1_in = (nlon1+1)*(nlat1+1); */
2876 /* n2_in = (nlon2+1)*(nlat2+1); */
2877
2878 /* lon1_in[0] = 305; lat1_in[0] = -35.26; */
2879 /* lon1_in[1] = 306; lat1_in[1] = -35.99; */
2880 /* lon1_in[2] = 305; lat1_in[2] = -33.80; */
2881 /* lon1_in[3] = 306; lat1_in[3] = -34.51; */
2882 /* lon2_in[0] = 305; lat2_in[0] = -34; */
2883 /* lon2_in[1] = 307.5; lat2_in[1] = -34; */
2884 /* lon2_in[2] = 305; lat2_in[2] = -32; */
2885 /* lon2_in[3] = 307.5; lat2_in[3] = -32; */
2886
2887 nlon1 = 2;
2888 nlat1 = 2;
2889 nlon2 = 1;
2890 nlat2 = 1;
2891 n1_in = (nlon1+1)*(nlat1+1);
2892 n2_in = (nlon2+1)*(nlat2+1);
2893
2894 lon1_in[0] = 111.3; lat1_in[0] = 1.591;
2895 lon1_in[1] = 109.7; lat1_in[1] = 2.926;
2896 lon1_in[2] = 108.2; lat1_in[2] = 4.256;
2897 lon1_in[3] = 110.0; lat1_in[3] = 0.000;
2898 lon1_in[4] = 108.4; lat1_in[4] = 1.335;
2899 lon1_in[5] = 106.8; lat1_in[5] = 2.668;
2900 lon1_in[6] = 108.7; lat1_in[6] = -1.591;
2901 lon1_in[7] = 107.1; lat1_in[7] = -0.256;
2902 lon1_in[8] = 105.5; lat1_in[8] = 1.078;
2903
2904 lon2_in[0] = 107.5; lat2_in[0] = 0;
2905 lon2_in[1] = 110; lat2_in[1] = 0;
2906 lon2_in[2] = 107.5; lat2_in[2] = 2;
2907 lon2_in[3] = 110; lat2_in[3] = 2;
2908
2909 break;
2910
2911 case 12:
2912 /****************************************************************
2913
2914 test : create_xgrid_great_circle
2915 box 1: (20,10), (20,12), (22,12), (22,10)
2916 box 2: (21,11), (21,14), (24,14), (24,11)
2917 out : (21, 12.0018), (22, 12), (22, 11.0033), (21, 11)
2918
2919 ****************************************************************/
2920 nlon1 = 2;
2921 nlat1 = 2;
2922 nlon2 = 3;
2923 nlat2 = 3;
2924 n1_in = (nlon1+1)*(nlat1+1);
2925 n2_in = (nlon2+1)*(nlat2+1);
2926
2927 /* first a simple lat-lon grid box to clip another lat-lon grid box */
2928 for(j=0; j<=nlat1; j++) for(i=0; i<=nlon1; i++){
2929 lon1_in[j*(nlon1+1)+i] = 20.0 + (i-1)*2.0;
2930 lat1_in[j*(nlon1+1)+i] = 10.0 + (j-1)*2.0;
2931 }
2932 for(j=0; j<=nlat2; j++) for(i=0; i<=nlon2; i++){
2933 lon2_in[j*(nlon2+1)+i] = 19.0 + (i-1)*2.0;
2934 lat2_in[j*(nlon2+1)+i] = 9.0 + (j-1)*2.0;
2935 }
2936
2937 break;
2938
2939 case 13:
2940
2941 nlon1 = 1;
2942 nlat1 = 1;
2943 nlon2 = 1;
2944 nlat2 = 1;
2945 n1_in = (nlon1+1)*(nlat1+1);
2946 n2_in = (nlon2+1)*(nlat2+1);
2947
2948 /* lon1_in[0] = ; lat1_in[0] = ; */
2949 /* lon1_in[1] = ; lat1_in[1] = ; */
2950 /* lon1_in[2] = ; lat1_in[2] = ; */
2951 /* lon1_in[3] = ; lat1_in[3] = ; */
2952 /* lon2_in[0] = ; lat2_in[0] = ; */
2953 /* lon2_in[1] = ; lat2_in[1] = ; */
2954 /* lon2_in[2] = ; lat2_in[2] = ; */
2955 /* lon2_in[3] = ; lat2_in[3] = ; */
2956
2957 /* lon1_in[0] = 1.35536; lat1_in[0] = 1.16251; */
2958 /* lon1_in[1] = 1.36805; lat1_in[1] = 1.15369; */
2959 /* lon1_in[2] = 1.37843; lat1_in[2] = 1.16729; */
2960 /* lon1_in[3] = 1.39048; lat1_in[3] = 1.15826; */
2961 /* lon2_in[0] = 1.34611; lat2_in[0] = 1.16372; */
2962 /* lon2_in[1] = 1.35616; lat2_in[1] = 1.15802; */
2963 /* lon2_in[2] = 1.35143; lat2_in[2] = 1.16509; */
2964 /* lon2_in[3] = 1.36042; lat2_in[3] = 1.15913; */
2965
2966 /* lon1_in[0] = 12.508065121288551; lat1_in[0] = -87.445883646793547; */
2967 /* lon1_in[1] = 325.425637772; lat1_in[1] = -86.481216821859505; */
2968 /* lon1_in[2] = 97.5; lat1_in[2] = -89.802136057677174; */
2969 /* lon1_in[3] = 277.5; lat1_in[3] = -87.615232005344637; */
2970
2971 /* for(j=0; j<=nlat2; j++) for(i=0; i<=nlon2; i++) { */
2972 /* lon2_in[j*(nlon2+1)+i] = -280.0 + i*1.0; */
2973 /* lat2_in[j*(nlon2+1)+i] = -90.0 + j*8.0; */
2974 /* } */
2975 lon1_in[0] = 120.369397984526174; lat1_in[0] = 16.751543427495864;
2976 lon1_in[1] = 119.999999999999986; lat1_in[1] = 16.751871929590038;
2977 lon1_in[2] = 120.369397846883501; lat1_in[2] = 16.397797979598028;
2978 lon1_in[3] = 119.999999999999986; lat1_in[3] = 16.398120477217255;
2979 lon2_in[0] = 120.369415056522087; lat2_in[0] = 16.752176828509153;
2980 lon2_in[1] = 119.999999999999986; lat2_in[1] = 16.752505523196167;
2981 lon2_in[2] = 120.369415056522087; lat2_in[2] = 16.397797949548146;
2982 lon2_in[3] = 119.999999999999986; lat2_in[3] = 16.398120477217255;
2983
2984 break;
2985
2986
2987 default:
2988 error_handler("test_create_xgrid: incorrect case number");
2989 }
2990
2991 /* convert to radian */
2992
2993 for(i=0; i<n1_in; i++) {
2994 lon1_in[i] *= D2R; lat1_in[i] *=D2R;
2995 }
2996 for(i=0; i<n2_in; i++) {
2997 lon2_in[i] *= D2R; lat2_in[i] *=D2R;
2998 }
2999
3000
3001 printf("\n*********************************************************\n");
3002 printf("\n Case %d \n", n);
3003 printf("\n*********************************************************\n");
3004
3005
3006 if( n > 10 ) {
3007 int nxgrid;
3008 int *i1, *j1, *i2, *j2;
3009 double *xarea, *xclon, *xclat, *mask1;
3010
3011 mask1 = (double *)malloc(nlon1*nlat1*sizeof(double));
3012 i1 = (int *)malloc(MAXXGRID*sizeof(int));
3013 j1 = (int *)malloc(MAXXGRID*sizeof(int));
3014 i2 = (int *)malloc(MAXXGRID*sizeof(int));
3015 j2 = (int *)malloc(MAXXGRID*sizeof(int));
3016 xarea = (double *)malloc(MAXXGRID*sizeof(double));
3017 xclon = (double *)malloc(MAXXGRID*sizeof(double));
3018 xclat = (double *)malloc(MAXXGRID*sizeof(double));
3019
3020 for(i=0; i<nlon1*nlat1; i++) mask1[i] = 1.0;
3021
3022 nxgrid = create_xgrid_great_circle(&nlon1, &nlat1, &nlon2, &nlat2, lon1_in, lat1_in,
3023 lon2_in, lat2_in, mask1, i1, j1, i2, j2,
3024 xarea, xclon, xclat);
3025 printf("\n*********************************************************\n");
3026 printf("\n First input grid box longitude, latitude \n \n");
3027 for(i=0; i<n1_in; i++) printf(" %g, %g \n", lon1_in[i]*R2D, lat1_in[i]*R2D);
3028
3029 printf("\n Second input grid box longitude, latitude \n \n");
3030 for(i=0; i<n2_in; i++) printf(" %g, %g \n", lon2_in[i]*R2D, lat2_in[i]*R2D);
3031
3032 printf("\n Number of exchange grid is %d\n", nxgrid);
3033 for(i=0; i<nxgrid; i++) {
3034 printf("(i1,j1)=(%d,%d), (i2,j2)=(%d, %d), xgrid_area=%g, xgrid_clon=%g, xgrid_clat=%g\n",
3035 i1[i], j1[i], i2[i], j2[i], xarea[i], xclon[i], xclat[i]);
3036 }
3037
3038 /* comparing the area sum of exchange grid and grid1 area */
3039 {
3040 double *x1, *y1, *z1, *area1;
3041 double area_sum;
3042 int i;
3043 area_sum = 0.0;
3044
3045 for(i=0; i<nxgrid; i++) {
3046 area_sum+= xarea[i];
3047 }
3048
3049 area1 = (double *)malloc((nlon1)*(nlat1)*sizeof(double));
3050 get_grid_great_circle_area_(&nlon1, &nlat1, lon1_in, lat1_in, area1);
3051
3052 printf("xgrid area sum is %g, grid 1 area is %g\n", area_sum, area1[0]);
3053 }
3054
3055 printf("\n");
3056 free(i1);
3057 free(i2);
3058 free(j1);
3059 free(j2);
3060 free(xarea);
3061 free(xclon);
3062 free(xclat);
3063 free(mask1);
3064 }
3065 else {
3066 latlon2xyz(n1_in, lon1_in, lat1_in, x1_in, y1_in, z1_in);
3067 latlon2xyz(n2_in, lon2_in, lat2_in, x2_in, y2_in, z2_in);
3068
3069 n_out = clip_2dx2d_great_circle(x1_in, y1_in, z1_in, 4, x2_in, y2_in, z2_in, n2_in,
3070 x_out, y_out, z_out);
3071 xyz2latlon(n_out, x_out, y_out, z_out, lon_out, lat_out);
3072
3073 printf("\n*********************************************************\n");
3074 printf("\n First input grid box longitude, latitude \n \n");
3075 for(i=0; i<n1_in; i++) printf(" %g, %g \n", lon1_in[i]*R2D, lat1_in[i]*R2D);
3076
3077 printf("\n Second input grid box longitude, latitude \n \n");
3078 for(i=0; i<n2_in; i++) printf(" %g, %g \n", lon2_in[i]*R2D, lat2_in[i]*R2D);
3079
3080 printf("\n output clip grid box longitude, latitude for case 1 \n \n");
3081 for(i=0; i<n_out; i++) printf(" %g, %g \n", lon_out[i]*R2D, lat_out[i]*R2D);
3082 printf("\n");
3083 }
3084 }
3085 }
3086
3087
3088 #endif