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