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