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