1 /***********************************************************************
2 * GNU Lesser General Public License
4 * This file is part of the GFDL Flexible Modeling System (FMS).
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.
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
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 **********************************************************************/
22 #include "mosaic_util.h"
23 #include "create_xgrid.h"
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)
39 * \brief Grid creation and calculation functions for use in @ref mosaic_mod
42 /*******************************************************************************
44 return constants MAXXGRID.
45 *******************************************************************************/
46 int get_maxxgrid(void)
51 int get_maxxgrid_(void)
53 return get_maxxgrid();
57 /*******************************************************************************
58 void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, const double *area)
60 *******************************************************************************/
61 void get_grid_area_(const int *nlon
, const int *nlat
, const double *lon
, const double *lat
, double *area
)
63 get_grid_area(nlon
, nlat
, lon
, lat
, area
);
66 void get_grid_area(const int *nlon
, const int *nlat
, const double *lon
, const double *lat
, double *area
)
68 int nx
, ny
, nxp
, i
, j
, n_in
;
69 double x_in
[20], y_in
[20];
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
);
91 /*******************************************************************************
92 void get_grid_area_ug(const int *npts, const double *lon, const double *lat, const double *area)
94 *******************************************************************************/
95 void get_grid_area_ug_(const int *npts
, const double *lon
, const double *lat
, double *area
)
97 get_grid_area_ug(npts
, lon
, lat
, area
);
100 void get_grid_area_ug(const int *npts
, const double *lon
, const double *lat
, double *area
)
103 double x_in
[20], y_in
[20];
108 for(l
=0; l
<nl
; l
++) {
110 x_in
[1] = lon
[l
*nv
+1];
111 x_in
[2] = lon
[l
*nv
+2];
112 x_in
[3] = lon
[l
*nv
+3];
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
);
121 } /* get_grid_area_ug */
124 void get_grid_great_circle_area_(const int *nlon
, const int *nlat
, const double *lon
, const double *lat
, double *area
)
126 get_grid_great_circle_area(nlon
, nlat
, lon
, lat
, area
);
130 void get_grid_great_circle_area(const int *nlon
, const int *nlat
, const double *lon
, const double *lat
, double *area
)
132 int nx
, ny
, nxp
, nyp
, i
, j
;
134 struct Node
*grid
=NULL
;
135 double *x
=NULL
, *y
=NULL
, *z
=NULL
;
143 x
= (double *)malloc(nxp
*nyp
*sizeof(double));
144 y
= (double *)malloc(nxp
*nyp
*sizeof(double));
145 z
= (double *)malloc(nxp
*nyp
*sizeof(double));
147 latlon2xyz(nxp
*nyp
, lon
, lat
, x
, y
, z
);
149 for(j
=0; j
<ny
; j
++) for(i
=0; i
< nx
; i
++) {
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
);
168 } /* get_grid_great_circle_area */
170 void get_grid_great_circle_area_ug_(const int *npts
, const double *lon
, const double *lat
, double *area
)
172 get_grid_great_circle_area_ug(npts
, lon
, lat
, area
);
176 void get_grid_great_circle_area_ug(const int *npts
, const double *lon
, const double *lat
, double *area
)
180 struct Node
*grid
=NULL
;
181 double *x
=NULL
, *y
=NULL
, *z
=NULL
;
186 x
= (double *)malloc(nl
*nv
*sizeof(double));
187 y
= (double *)malloc(nl
*nv
*sizeof(double));
188 z
= (double *)malloc(nl
*nv
*sizeof(double));
190 latlon2xyz(nl
*nv
, lon
, lat
, x
, y
, z
);
192 for(l
=0; l
<nv
; l
++) {
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
);
211 } /* get_grid_great_circle_area_ug */
213 void get_grid_area_dimensionless(const int *nlon
, const int *nlat
, const double *lon
, const double *lat
, double *area
)
215 int nx
, ny
, nxp
, i
, j
, n_in
;
216 double x_in
[20], y_in
[20];
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
);
235 } /* get_grid_area */
239 void get_grid_area_no_adjust(const int *nlon
, const int *nlat
, const double *lon
, const double *lat
, double *area
)
241 int nx
, ny
, nxp
, i
, j
, n_in
;
242 double x_in
[20], y_in
[20];
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
];
258 area
[j
*nx
+i
] = poly_area_no_adjust(x_in
, y_in
, n_in
);
261 } /* get_grid_area_no_adjust */
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
)
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
);
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
)
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
;
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
];
310 /* This is just a temporary fix to solve the issue that there is one point in zonal direction */
312 get_grid_area(nlon_in
, nlat_in
, tmpx
, tmpy
, area_in
);
314 get_grid_area_no_adjust(nlon_in
, nlat_in
, tmpx
, tmpy
, area_in
);
316 get_grid_area(nlon_out
, nlat_out
, lon_out
, lat_out
, area_out
);
320 for(j1
=0; j1
<ny1
; j1
++) for(i1
=0; i1
<nx1
; i1
++) if( mask_in
[j1
*nx1
+i1
] > MASK_THRESH
) {
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
++) {
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;
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);
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
;
353 if(nxgrid
> MAXXGRID
) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
364 } /* create_xgrid_1dx2d_order1 */
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
)
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
);
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
)
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
;
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
];
412 /* This is just a temporary fix to solve the issue that there is one point in zonal direction */
414 get_grid_area(nlon_in
, nlat_in
, tmpx
, tmpy
, area_in
);
416 get_grid_area_no_adjust(nlon_in
, nlat_in
, tmpx
, tmpy
, area_in
);
418 get_grid_area_ug(npts_out
, lon_out
, lat_out
, area_out
);
422 for(j1
=0; j1
<ny1
; j1
++) for(i1
=0; i1
<nx1
; i1
++) if( mask_in
[j1
*nx1
+i1
] > MASK_THRESH
) {
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
++) {
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;
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);
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
;
454 if(nxgrid
> MAXXGRID
) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
465 } /* create_xgrid_1dx2d_order1_ug */
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
)
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
);
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
)
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
;
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
];
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
);
518 for(j1
=0; j1
<ny1
; j1
++) for(i1
=0; i1
<nx1
; i1
++) if( mask_in
[j1
*nx1
+i1
] > MASK_THRESH
) {
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
++) {
524 double xarea
, lon_in_avg
;
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;
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
);
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
);
554 if(nxgrid
> MAXXGRID
) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
564 } /* create_xgrid_1dx2d_order2 */
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
)
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
);
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
)
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
;
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
];
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
);
622 for(j2
=0; j2
<ny2
; j2
++) for(i2
=0; i2
<nx2
; i2
++) {
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
) {
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;
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
];
642 n_in
= fix_lon(x_in
, y_in
, 4, (ll_lon
+ur_lon
)/2);
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
;
654 if(nxgrid
> MAXXGRID
) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
665 } /* create_xgrid_2dx1d_order1 */
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
)
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
);
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
)
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
];
697 double *area_in
, *area_out
, min_area
;
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
];
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
);
726 for(j2
=0; j2
<ny2
; j2
++) for(i2
=0; i2
<nx2
; i2
++) {
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
) {
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;
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
];
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
);
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
);
761 if(nxgrid
> MAXXGRID
) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
772 } /* create_xgrid_2dx1d_order2 */
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
)
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
);
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
)
800 int nx1
, nx2
, ny1
, ny2
, nx1p
, nx2p
, nxgrid
;
801 double *area_in
, *area_out
;
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
;
811 int nthreads
, nxgrid_block_max
;
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
);
828 nthreads
= omp_get_num_threads();
833 istart2
= (int *)malloc(nblocks
*sizeof(int));
834 iend2
= (int *)malloc(nblocks
*sizeof(int));
836 pstart
= (int *)malloc(nblocks
*sizeof(int));
837 pnxgrid
= (int *)malloc(nblocks
*sizeof(int));
839 nxgrid_block_max
= MAXXGRID
/nblocks
;
841 for(m
=0; m
<nblocks
; m
++) {
843 pstart
[m
] = m
*nxgrid_block_max
;
851 pxgrid_area
= xgrid_area
;
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));
862 nblks_left
= nblocks
;
864 for(m
=0; m
<nblocks
; m
++) {
866 npts_my
= npts_left
/nblks_left
;
867 iend2
[m
] = istart2
[m
] + npts_my
- 1;
869 npts_left
-= npts_my
;
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));
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)
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
];
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
];
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
);
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
];
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)
922 for(m
=0; m
<nblocks
; m
++) {
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
];
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
];
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*/
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
];
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
;
962 for (l
=0; l
<n2_in
; l
++) x2_in
[l
] += TPI
;
964 else if (dx
> M_PI
) {
967 for (l
=0; l
<n2_in
; l
++) x2_in
[l
] -= TPI
;
970 /* x2_in should in the same range as x1_in after lon_fix, so no need to
971 consider cyclic condition
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) {
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
) {
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;
985 pxgrid_area
[nn
] = xarea
;
998 /*copy data if nblocks > 1 */
1000 nxgrid
= pnxgrid
[0];
1010 for(m
=0; m
<nblocks
; m
++) {
1011 for(i
=0; i
<pnxgrid
[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
];
1030 free(lon_out_min_list
);
1031 free(lon_out_max_list
);
1032 free(lat_out_min_list
);
1033 free(lat_out_max_list
);
1041 }/* get_xgrid_2Dx2D_order1 */
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
)
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
);
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
)
1068 int nx1
, nx2
, ny1
, ny2
, nx1p
, nx2p
, nxgrid
;
1069 double *area_in
, *area_out
;
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
;
1079 int nthreads
, nxgrid_block_max
;
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
);
1094 #if defined(_OPENMP)
1095 #pragma omp parallel
1096 nthreads
= omp_get_num_threads();
1101 istart2
= (int *)malloc(nblocks
*sizeof(int));
1102 iend2
= (int *)malloc(nblocks
*sizeof(int));
1104 pstart
= (int *)malloc(nblocks
*sizeof(int));
1105 pnxgrid
= (int *)malloc(nblocks
*sizeof(int));
1107 nxgrid_block_max
= MAXXGRID
/nblocks
;
1109 for(m
=0; m
<nblocks
; m
++) {
1111 pstart
[m
] = m
*nxgrid_block_max
;
1119 pxgrid_area
= xgrid_area
;
1120 pxgrid_clon
= xgrid_clon
;
1121 pxgrid_clat
= xgrid_clat
;
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));
1133 npts_left
= nx2
*ny2
;
1134 nblks_left
= nblocks
;
1136 for(m
=0; m
<nblocks
; m
++) {
1138 npts_my
= npts_left
/nblks_left
;
1139 iend2
[m
] = istart2
[m
] + npts_my
- 1;
1141 npts_left
-= npts_my
;
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)
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
];
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
];
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
);
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
];
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)
1195 for(m
=0; m
<nblocks
; m
++) {
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
];
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
];
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
];
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
;
1235 for (l
=0; l
<n2_in
; l
++) x2_in
[l
] += TPI
;
1237 else if (dx
> M_PI
) {
1240 for (l
=0; l
<n2_in
; l
++) x2_in
[l
] -= TPI
;
1243 /* x2_in should in the same range as x1_in after lon_fix, so no need to
1244 consider cyclic condition
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) {
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
) {
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
);
1270 /*copy data if nblocks > 1 */
1272 nxgrid
= pnxgrid
[0];
1284 for(m
=0; m
<nblocks
; m
++) {
1285 for(i
=0; i
<pnxgrid
[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
];
1308 free(lon_out_min_list
);
1309 free(lon_out_max_list
);
1310 free(lat_out_min_list
);
1311 free(lat_out_max_list
);
1319 }/* get_xgrid_2Dx2D_order2 */
1322 /*******************************************************************************
1323 Sutherland-Hodgeman algorithm sequentially clips parts outside 4 boundaries
1324 *******************************************************************************/
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
[])
1329 double x_tmp
[MV
], y_tmp
[MV
], x_last
, y_last
;
1330 int i_in
, i_out
, n_out
, inside_last
, inside
;
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
++) {
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
);
1344 /* if "to" point is right of LEFT boundary, output it */
1346 x_tmp
[i_out
] = lon_in
[i_in
];
1347 y_tmp
[i_out
++] = lat_in
[i_in
];
1349 x_last
= lon_in
[i_in
];
1350 y_last
= lat_in
[i_in
];
1351 inside_last
= inside
;
1353 if (!(n_out
=i_out
)) return(0);
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
++) {
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
);
1368 /* if "to" point is left of RIGHT boundary, output it */
1370 lon_out
[i_out
] = x_tmp
[i_in
];
1371 lat_out
[i_out
++] = y_tmp
[i_in
];
1374 x_last
= x_tmp
[i_in
];
1375 y_last
= y_tmp
[i_in
];
1376 inside_last
= inside
;
1378 if (!(n_out
=i_out
)) return(0);
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
++) {
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
);
1392 /* if "to" point is above BOTTOM boundary, output it */
1394 x_tmp
[i_out
] = lon_out
[i_in
];
1395 y_tmp
[i_out
++] = lat_out
[i_in
];
1397 x_last
= lon_out
[i_in
];
1398 y_last
= lat_out
[i_in
];
1399 inside_last
= inside
;
1401 if (!(n_out
=i_out
)) return(0);
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
++) {
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
);
1415 /* if "to" point is below TOP boundary, output it */
1417 lon_out
[i_out
] = x_tmp
[i_in
];
1418 lat_out
[i_out
++] = y_tmp
[i_in
];
1420 x_last
= x_tmp
[i_in
];
1421 y_last
= y_tmp
[i_in
];
1422 inside_last
= inside
;
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 *******************************************************************************/
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
[])
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
;
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 */
1445 for(i1
=0; i1
<n1_in
; i1
++) {
1446 lon_tmp
[i1
] = lon1_in
[i1
];
1447 lat_tmp
[i1
] = lat1_in
[i1
];
1449 x2_0
= lon2_in
[n2_in
-1];
1450 y2_0
= lat2_in
[n2_in
-1];
1451 for(i2
=0; i2
<n2_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
++) {
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 */
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>");
1475 lon_out
[i_out
] = (dx2
*ds1
- dx1
*ds2
)/determ
;
1476 lat_out
[i_out
++] = (dy2
*ds1
- dy1
*ds2
)/determ
;
1481 lon_out
[i_out
] = x1_1
;
1482 lat_out
[i_out
++] = y1_1
;
1486 inside_last
= inside
;
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
];
1493 /* shift the starting point */
1500 /*#define debug_test_create_xgrid*/
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
)
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
);
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
)
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
;
1528 double *area1
, *area2
, min_area
;
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));
1548 latlon2xyz(nx1p
*ny1p
, lon_in
, lat_in
, x1
, y1
, z1
);
1549 latlon2xyz(nx2p
*ny2p
, lon_out
, lat_out
, x2
, y2
, z2
);
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
);
1558 for(j1
=0; j1
<ny1
; j1
++) for(i1
=0; i1
<nx1
; i1
++) if( mask_in
[j1
*nx1
+i1
] > MASK_THRESH
) {
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
];
1567 for(j2
=0; j2
<ny2
; j2
++) for(i2
=0; i2
<nx2
; i2
++) {
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
];
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
);
1586 xgrid_area
[nxgrid
] = xarea
;
1587 xgrid_clon
[nxgrid
] = 0; /*z1l: will be developed very soon */
1588 xgrid_clat
[nxgrid
] = 0;
1594 if(nxgrid
> MAXXGRID
) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
1613 }/* create_xgrid_great_circle */
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
)
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
);
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
)
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
;
1641 double *area1
, *area2
, min_area
;
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));
1659 latlon2xyz(nx1p
*ny1p
, lon_in
, lat_in
, x1
, y1
, z1
);
1660 latlon2xyz(npts2
*nv
, lon_out
, lat_out
, x2
, y2
, z2
);
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
);
1669 for(j1
=0; j1
<ny1
; j1
++) for(i1
=0; i1
<nx1
; i1
++) if( mask_in
[j1
*nx1
+i1
] > MASK_THRESH
) {
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
];
1678 for(l2
=0; l2
<npts2
; l2
++) {
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
];
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
);
1697 xgrid_area
[nxgrid
] = xarea
;
1698 xgrid_clon
[nxgrid
] = 0; /*z1l: will be developed very soon */
1699 xgrid_clat
[nxgrid
] = 0;
1704 if(nxgrid
> MAXXGRID
) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID");
1723 }/* create_xgrid_great_circle_ug */
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 *******************************************************************************/
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
[])
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
;
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];
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
;
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;
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;
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;
1787 grid1List
= getNext();
1788 grid2List
= getNext();
1789 intersectList
= getNext();
1790 polyList
= getNext();
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
);
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");
1803 /* first check number of points in grid1 is inside grid2 */
1807 if(insidePolygon(temp
, grid2List
))
1811 temp
= getNextNode(temp
);
1814 #ifdef debug_test_create_xgrid
1815 printf("\nNOTE from clip_2dx2d_great_circle: begin to set inside value of grid2List\n");
1817 /* check if grid2List is inside grid1List */
1821 if(insidePolygon(temp
, grid1List
))
1825 temp
= getNextNode(temp
);
1828 /* make sure the grid box is clockwise */
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");
1836 #ifdef debug_test_create_xgrid
1837 printNode(grid1List
, "grid1List");
1838 printNode(grid2List
, "grid2List");
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
1846 for(i1
=0; i1
<npts1
; i1
++) {
1847 getCoordinates(temp
, pt1
[i1
]);
1851 for(i2
=0; i2
<npts2
; i2
++) {
1852 getCoordinates(temp
, pt2
[i2
]);
1856 firstIntersect
=getNext();
1857 curIntersect
= getNext();
1859 #ifdef debug_test_create_xgrid
1860 printf("\n\n************************ Start line_intersect_2D_3D ******************************\n");
1862 /* first find all the intersection points */
1864 for(i1
=0; i1
<npts1
; i1
++) {
1868 for(i2
=0; i2
<npts2
; i2
++) {
1870 i2p2
= (i2
+2)%npts2
;
1874 #ifdef debug_test_create_xgrid
1875 printf("\n******************************************************************************\n");
1876 printf(" i1 = %d, i2 = %d \n", i1
, i2
);
1877 printf("********************************************************************************\n");
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 */
1882 /* add the intersection into intersetList, The intersection might already be in
1883 intersectList and will be taken care addIntersect
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 */
1889 insertIntersect(grid1List
, intersect
[0], intersect
[1], intersect
[2], 0.0, u2
, inbound
, p1_1
[0], p1_1
[1], p1_1
[2]);
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 */
1895 p1_1
[0] = intersect
[0];
1896 p1_1
[1] = intersect
[1];
1897 p1_1
[2] = intersect
[2];
1900 p1_0
[0] = intersect
[0];
1901 p1_0
[1] = intersect
[1];
1902 p1_0
[2] = intersect
[2];
1904 /* add the intersection into the grid2List */
1906 insertIntersect(grid2List
, intersect
[0], intersect
[1], intersect
[2], 0.0, u1
, 0, p2_1
[0], p2_1
[1], p2_1
[2]);
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 */
1911 p2_1
[0] = intersect
[0];
1912 p2_1
[1] = intersect
[1];
1913 p2_1
[2] = intersect
[2];
1916 p2_0
[0] = intersect
[0];
1917 p2_0
[1] = intersect
[1];
1918 p2_0
[2] = intersect
[2];
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
1929 /* get the first point in intersectList has inbound = 2, if not, set inbound value */
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
) {
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;
1948 /* if has_inbound = 1, find the overlapping */
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");
1961 temp1
= getNode(grid1List
, *firstIntersect
);
1962 if( temp1
== NULL
) {
1963 double lon
[10], lat
[10];
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
);
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
);
1972 error_handler("firstIntersect is not in the grid1List");
1974 addNode(polyList
, *firstIntersect
);
1976 #ifdef debug_test_create_xgrid
1977 printNode(polyList
, "polyList at stage 1");
1980 /* Loop over the grid1List and grid2List to find again the firstIntersect */
1981 curList
= grid1List
;
1984 /* Loop through curList to find the next intersection, the loop will end
1985 when come back to firstIntersect
1987 copyNode(curIntersect
, *firstIntersect
);
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");
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
;
2001 maxiter2
= length(curList
);
2004 /* Loop until find the next intersection */
2005 while( iter2
< maxiter2
) {
2006 int temp2IsIntersect
;
2008 temp2IsIntersect
= 0;
2009 if( isIntersect( *temp2
) ) { /* copy the point and switch to the grid2List */
2012 /* first check if temp2 is the firstIntersect */
2013 if( sameNode( *temp2
, *firstIntersect
) ) {
2018 temp3
= temp2
->Next
;
2019 if( temp3
== NULL
) temp3
= curList
;
2020 if( temp3
== NULL
) error_handler("creat_xgrid.c: temp3 can not be NULL");
2022 /* if next node is inside or an intersection,
2023 need to keep on curList
2025 temp2IsIntersect
= 1;
2026 if( isIntersect(*temp3
) || (temp3
->isInside
== 1) ) found2
= 0;
2029 copyNode(curIntersect
, *temp2
);
2033 addNode(polyList
, *temp2
);
2034 #ifdef debug_test_create_xgrid
2035 printNode(polyList
, "polyList at stage 2");
2037 if(temp2IsIntersect
) {
2041 temp2
= temp2
->Next
;
2042 if( temp2
== NULL
) temp2
= curList
;
2047 if( !found2
) error_handler(" not found the next intersection ");
2049 /* if find the first intersection, the poly found */
2050 if( sameNode( *curIntersect
, *firstIntersect
) ) {
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");
2063 /* switch curList */
2064 if( curListNum
== 0) {
2065 curList
= grid2List
;
2069 curList
= grid1List
;
2074 if(!found1
) error_handler("not return back to the first intersection");
2076 /* currently we are only clipping convex polygon to convex polygon */
2077 if( nintersect
> 0) error_handler("After clipping, nintersect should be 0");
2079 /* copy the polygon to x_out, y_out, z_out */
2081 while (temp1
!= NULL
) {
2082 getCoordinate(*temp1
, x_out
+n_out
, y_out
+n_out
, z_out
+n_out
);
2083 temp1
= temp1
->Next
;
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");
2094 /* check if grid1 is inside grid2 */
2096 /* first check number of points in grid1 is inside grid2 */
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");
2105 if(temp
->intersect
!= 1) {
2106 #ifdef debug_test_create_xgrid
2107 printf("grid1->isInside = %d\n", temp
->isInside
);
2109 if( temp
->isInside
== 1) n1in2
++;
2111 temp
= getNextNode(temp
);
2113 if(npts1
==n1in2
) { /* grid1 is inside grid2 */
2118 getCoordinate(*temp
, &x_out
[n
], &y_out
[n
], &z_out
[n
]);
2120 temp
= getNextNode(temp
);
2123 if(n_out
>0) return n_out
;
2126 /* check if grid2List is inside grid1List */
2129 #ifdef debug_test_create_xgrid
2130 printf("\nNOTE from clip_2dx2d_great_circle: check if grid2 is inside grid1\n");
2136 if(temp
->intersect
!= 1) {
2137 #ifdef debug_test_create_xgrid
2138 printf("grid2->isInside = %d\n", temp
->isInside
);
2140 if( temp
->isInside
== 1) n2in1
++;
2142 temp
= getNextNode(temp
);
2145 if(npts2
==n2in1
) { /* grid2 is inside grid1 */
2150 getCoordinate(*temp
, &x_out
[n
], &y_out
[n
], &z_out
[n
]);
2152 temp
= getNextNode(temp
);
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)
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
){
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.
2181 double p1
[3], v1
[3], v2
[3];
2182 double c1
[3], c2
[3], c3
[3];
2183 double coincident
, sense
, norm
;
2185 int is_inter1
, is_inter2
;
2189 /* first check if any vertices are the same */
2190 if(samePoint(a1
[0], a1
[1], a1
[2], q1
[0], q1
[1], q1
[2])) {
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
);
2201 else if (samePoint(a1
[0], a1
[1], a1
[2], q2
[0], q2
[1], q2
[2])) {
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
);
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
);
2218 intersect
[0] = a2
[0];
2219 intersect
[1] = a2
[1];
2220 intersect
[2] = a2
[2];
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
);
2229 intersect
[0] = a2
[0];
2230 intersect
[1] = a2
[1];
2231 intersect
[2] = a2
[2];
2236 /* Load points defining plane into variable (these are supposed to be in counterclockwise order) */
2247 /* Intersect the segment with the plane */
2248 is_inter1
= intersect_tri_with_line(plane
, a1
, a2
, plane_p
, u_a
);
2253 if(fabs(*u_a
) < EPSLN8
) *u_a
= 0;
2254 if(fabs(*u_a
-1) < EPSLN8
) *u_a
= 1;
2257 #ifdef debug_test_create_xgrid
2258 printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f\n", *u_a
);
2262 if( (*u_a
< 0) || (*u_a
> 1) ) return 0;
2264 /* Load points defining plane into variable (these are supposed to be in counterclockwise order) */
2275 /* Intersect the segment with the plane */
2276 is_inter2
= intersect_tri_with_line(plane
, q1
, q2
, plane_p
, u_q
);
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
);
2288 if( (*u_q
< 0) || (*u_q
> 1) ) return 0;
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
);
2298 if(fabs(coincident
) < EPSLN30
) return 0;
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]);
2305 norm
= metric( intersect
);
2306 for(i
= 0; i
< 3; i
++) intersect
[i
] /= norm
;
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){
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];
2321 vect_cross(v1
, v2
, c1
);
2322 vect_cross(v1
, p1
, c2
);
2324 sense
= dot(c1
, c2
);
2326 if(sense
> 0) *inbound
= 2; /* v1 going into v2 in CCW sense */
2328 #ifdef debug_test_create_xgrid
2329 printf("\nNOTE from line_intersect_2D_3D: inbound=%d\n", *inbound
);
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 ---------------------------------------------------------------------------*/
2341 double poly_ctrlat(const double x
[], const double y
[], int n
)
2343 double ctrlat
= 0.0;
2348 double dx
= (x
[ip
]-x
[i
]);
2349 double dy
, avg_y
, hdy
;
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
;
2360 if ( fabs(hdy
)< SMALL_VALUE
) /* cheap area calculation along latitude */
2361 ctrlat
-= dx
*(2*cos(avg_y
) + lat2
*sin(avg_y
) - cos(lat1
) );
2363 ctrlat
-= dx
*( (sin(hdy
)/hdy
)*(2*cos(avg_y
) + lat2
*sin(avg_y
)) - cos(lat1
) );
2365 return (ctrlat
*RADIUS
*RADIUS
);
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
)
2374 double ctrlon
= 0.0;
2379 double phi1
, phi2
, dphi
, lat1
, lat2
, dphi1
, dphi2
;
2380 double f1
, f2
, fac
, fint
;
2387 if (dphi
==0.0) continue;
2389 f1
= 0.5*(cos(lat1
)*sin(lat1
)+lat1
);
2390 f2
= 0.5*(cos(lat2
)*sin(lat2
)+lat2
);
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
;
2400 if( dphi2
> M_PI
) dphi2
-= 2.0*M_PI
;
2401 if( dphi2
<-M_PI
) dphi2
+= 2.0*M_PI
;
2403 if(fabs(dphi2
-dphi1
) < M_PI
) {
2404 ctrlon
-= dphi
* (dphi1
*f1
+dphi2
*f2
)/2.0;
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
;
2417 return (ctrlon
*RADIUS
*RADIUS
);
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
)
2426 double dphi
= ur_lon
-ll_lon
;
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
);
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
)
2441 double phi1
, phi2
, dphi
, lat1
, lat2
, dphi1
, dphi2
;
2442 double f1
, f2
, fac
, fint
;
2443 double ctrlon
= 0.0;
2445 for( i
=0; i
<2; i
++) {
2449 lat1
= lat2
= ll_lat
;
2454 lat1
= lat2
= ur_lat
;
2457 f1
= 0.5*(cos(lat1
)*sin(lat1
)+lat1
);
2458 f2
= 0.5*(cos(lat2
)*sin(lat2
)+lat2
);
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
;
2467 if( dphi2
> M_PI
) dphi2
-= 2.0*M_PI
;
2468 if( dphi2
<-M_PI
) dphi2
+= 2.0*M_PI
;
2470 if(fabs(dphi2
-dphi1
) < M_PI
) {
2471 ctrlon
-= dphi
* (dphi1
*f1
+dphi2
*f2
)/2.0;
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
;
2483 return (ctrlon
*RADIUS
*RADIUS
);
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
)
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.));
2503 radius
= sqrt(radius
);
2507 } /* grid_box_radius */
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
)
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.));
2531 } /* dist_between_boxes */
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
)
2544 const double SMALL
= 1.e
-12;
2547 product
= ( x
-x0
)*(y1
-y0
) + (x0
-x1
)*(y
-y0
);
2548 return (product
<=SMALL
) ? 1:0;
2553 /* The following is a test program to test subroutines in create_xgrid.c */
2555 #ifdef test_create_xgrid
2557 #include "create_xgrid.h"
2560 #define D2R (M_PI/180)
2561 #define R2D (180/M_PI)
2562 #define MAXPOINT 1000
2564 int main(int argc
, char* argv
[])
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;
2579 for(n
=11; n
<=ntest
; n
++) {
2583 /****************************************************************
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)
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;
2604 /****************************************************************
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)
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;
2617 for(i
=0; i
<n2_in
; i
++) {
2618 lon2_in
[i
] = lon1_in
[i
];
2619 lat2_in
[i
] = lat1_in
[i
];
2624 /****************************************************************
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)
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;
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;
2644 for(i=0; i<4; i++) {
2645 lon2_in[i] = lon1_in[i];
2646 lat2_in[i] = lat1_in[i];
2652 /****************************************************************
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)
2659 ****************************************************************/
2660 n1_in
= 4; n2_in
= 4;
2661 /* first a simple lat-lon grid box to clip another lat-lon grid box */
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;
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;
2675 /****************************************************************
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)
2682 ****************************************************************/
2683 n1_in
= 4; n2_in
= 4;
2684 /* first a simple lat-lon grid box to clip another lat-lon grid box */
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;
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;
2698 /****************************************************************
2700 test clip_2dx2d_great_cirle case 6: One cubic sphere grid arounc the pole with one tripolar grid box
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)
2707 ****************************************************************/
2708 n1_in
= 4; n2_in
= 4;
2709 /* first a simple lat-lon grid box to clip another lat-lon grid box */
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;
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;
2723 /****************************************************************
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)
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;
2744 /****************************************************************
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)
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;
2767 /****************************************************************
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)
2775 ****************************************************************/
2776 n1_in
= 4; n2_in
= 4;
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;
2789 /****************************************************************
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)
2797 ****************************************************************/
2798 n1_in
= 4; n2_in
= 4;
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;
2811 /****************************************************************
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)
2819 ****************************************************************/
2824 n1_in
= (nlon1
+1)*(nlat1
+1);
2825 n2_in
= (nlon2
+1)*(nlat2
+1);
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;
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; */
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;
2851 /* n1_in = (nlon1+1)*(nlat1+1); */
2852 /* n2_in = (nlon2+1)*(nlat2+1); */
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; */
2875 /* n1_in = (nlon1+1)*(nlat1+1); */
2876 /* n2_in = (nlon2+1)*(nlat2+1); */
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; */
2891 n1_in
= (nlon1
+1)*(nlat1
+1);
2892 n2_in
= (nlon2
+1)*(nlat2
+1);
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;
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;
2912 /****************************************************************
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)
2919 ****************************************************************/
2924 n1_in
= (nlon1
+1)*(nlat1
+1);
2925 n2_in
= (nlon2
+1)*(nlat2
+1);
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;
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;
2945 n1_in
= (nlon1
+1)*(nlat1
+1);
2946 n2_in
= (nlon2
+1)*(nlat2
+1);
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] = ; */
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; */
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; */
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; */
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;
2988 error_handler("test_create_xgrid: incorrect case number");
2991 /* convert to radian */
2993 for(i
=0; i
<n1_in
; i
++) {
2994 lon1_in
[i
] *= D2R
; lat1_in
[i
] *=D2R
;
2996 for(i
=0; i
<n2_in
; i
++) {
2997 lon2_in
[i
] *= D2R
; lat2_in
[i
] *=D2R
;
3001 printf("\n*********************************************************\n");
3002 printf("\n Case %d \n", n
);
3003 printf("\n*********************************************************\n");
3008 int *i1
, *j1
, *i2
, *j2
;
3009 double *xarea
, *xclon
, *xclat
, *mask1
;
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));
3020 for(i
=0; i
<nlon1
*nlat1
; i
++) mask1
[i
] = 1.0;
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
);
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
);
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
]);
3038 /* comparing the area sum of exchange grid and grid1 area */
3040 double *x1
, *y1
, *z1
, *area1
;
3045 for(i
=0; i
<nxgrid
; i
++) {
3046 area_sum
+= xarea
[i
];
3049 area1
= (double *)malloc((nlon1
)*(nlat1
)*sizeof(double));
3050 get_grid_great_circle_area_(&nlon1
, &nlat1
, lon1_in
, lat1_in
, area1
);
3052 printf("xgrid area sum is %g, grid 1 area is %g\n", area_sum
, area1
[0]);
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
);
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
);
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
);
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
);
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
);