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