21 #include "grid_utils.h"
22 #include "gradient_c2l.h"
47 void grad_c2l_(
const int *
nlon,
const int *
nlat,
const double *pin,
const double *dx,
const double *dy,
const double *
area,
48 const double *edge_w,
const double *edge_e,
const double *edge_s,
const double *edge_n,
49 const double *en_n,
const double *en_e,
const double *vlon,
const double *vlat,
50 double *grad_x,
double *grad_y,
const int *on_west_edge,
const int *on_east_edge,
51 const int *on_south_edge,
const int *on_north_edge)
53 grad_c2l(
nlon,
nlat, pin, dx, dy,
area, edge_w, edge_e, edge_s, edge_n, en_n, en_e, vlon, vlat, grad_x, grad_y,
54 on_west_edge, on_east_edge, on_south_edge, on_north_edge);
57 void grad_c2l(
const int *
nlon,
const int *
nlat,
const double *pin,
const double *dx,
const double *dy,
const double *
area,
58 const double *edge_w,
const double *edge_e,
const double *edge_s,
const double *edge_n,
59 const double *en_n,
const double *en_e,
const double *vlon,
const double *vlat,
60 double *grad_x,
double *grad_y,
const int *on_west_edge,
const int *on_east_edge,
61 const int *on_south_edge,
const int *on_north_edge)
64 double *pb, *pdx, *pdy, *grad3;
65 int nx, ny, nxp, nyp, i, j, m0, m1, n;
71 pb = (
double *)malloc(nxp*nyp*
sizeof(
double));
72 pdx = (
double *)malloc(3*nx*(ny+1)*
sizeof(double));
73 pdy = (
double *)malloc(3*(nx+1)*ny*
sizeof(double));
74 grad3 = (
double *)malloc(3*nx*ny*
sizeof(
double));
75 a2b_ord2(nx, ny, pin, edge_w, edge_e, edge_s, edge_n, pb, *on_west_edge, *on_east_edge,*on_south_edge, *on_north_edge);
77 for(j=0; j<nyp; j++)
for(i=0; i<nx; i++) {
81 pdx[3*m0+n] = 0.5*(pb[m1]+pb[m1+1])*dx[m0]*en_n[3*m0+n];
85 for(j=0; j<ny; j++)
for(i=0; i<nxp; i++) {
88 pdy[3*m0+n] = 0.5*(pb[m0]+pb[m0+nxp])*dy[m0]*en_e[3*m0+n];
93 for(j=0; j<ny; j++)
for(i=0; i<nx; i++) {
96 grad3[m0+n] = pdx[3*((j+1)*nx+i)+n]-pdx[m0+n]-pdy[3*(j*nxp+i)+n]+pdy[3*(j*nxp+i+1)+n];
101 for(j=0; j<ny; j++)
for(i=0; i<nx; i++) {
105 grad_x[m0] = (vlon[m1]*grad3[m1] + vlon[m1+1]*grad3[m1+1] + vlon[m1+2]*grad3[m1+2])/
area[m0];
106 grad_x[m0] *= RADIUS;
108 grad_y[m0] = (vlat[m1]*grad3[m1] + vlat[m1+1]*grad3[m1+1] + vlat[m1+2]*grad3[m1+2] )/
area[m0];
109 grad_y[m0] *= RADIUS;
123 void a2b_ord2(
int nx,
int ny,
const double *qin,
const double *edge_w,
const double *edge_e,
124 const double *edge_s,
const double *edge_n,
double *qout,
125 int on_west_edge,
int on_east_edge,
int on_south_edge,
int on_north_edge)
128 int istart, iend, jstart, jend;
130 const double r3 = 1./3.;
134 q1 = (
double *)malloc((nx+2)*
sizeof(double));
135 q2 = (
double *)malloc((ny+2)*
sizeof(double));
156 for(j=jstart; j<jend; j++)
for(i=istart; i<iend; i++) {
157 qout[j*nxp+i] = 0.25*(qin[j*(nx+2)+i] + qin[j*(nx+2)+i+1] +
158 qin[(j+1)*(nx+2)+i] + qin[(j+1)*(nx+2)+i+1] );
162 if(on_west_edge && on_south_edge)qout[ 0] = r3*(qin[1* (nx+2)+1 ]+qin[1* (nx+2) ]+qin[ 1]);
163 if(on_east_edge && on_south_edge)qout[ nx] = r3*(qin[1* (nx+2)+nx]+qin[ nx ]+qin[1* (nx+2)+nxp]);
164 if(on_east_edge && on_north_edge)qout[ny*nxp+nx] = r3*(qin[ny*(nx+2)+nx]+qin[ny*(nx+2)+nxp]+qin[nyp*(nx+2)+nx ]);
165 if(on_west_edge && on_north_edge)qout[ny*nxp ] = r3*(qin[ny*(nx+2)+1 ]+qin[ny*(nx+2) ]+qin[nyp*(nx+2)+1 ]);
169 for(j=jstart; j<=jend; j++){
170 q2[j] = 0.5*(qin[j*(nx+2)] + qin[j*(nx+2)+1]);
172 for(j=jstart; j<jend; j++){
173 qout[j*nxp] = edge_w[j]*q2[j] + (1-edge_w[j])*q2[j+1];
179 for(j=jstart; j<=jend; j++){
180 q2[j] = 0.5*(qin[j*(nx+2)+nx] + qin[j*(nx+2)+nxp]);
182 for(j=jstart; j<jend; j++){
183 qout[j*nxp+nx] = edge_e[j]*q2[j] + (1-edge_e[j])*q2[j+1];
189 for(i=istart; i<=iend; i++){
190 q1[i] = 0.5*(qin[i] + qin[(nx+2)+i]);
192 for(i=istart; i<iend; i++){
193 qout[i] = edge_s[i]*q1[i] + (1 - edge_s[i])*q1[i+1];
199 for(i=istart; i<=iend; i++){
200 q1[i] = 0.5*(qin[ny*(nx+2)+i] + qin[nyp*(nx+2)+i]);
202 for(i=istart; i<iend; i++){
203 qout[ny*nxp+i] = edge_n[i]*q1[i] + (1 - edge_n[i])*q1[i+1];
213 void get_edge(
int nx,
int ny,
const double *lont,
const double *latt,
214 const double *lonc,
const double *latc,
double *edge_w,
double *edge_e,
double *edge_s,
double *edge_n,
215 int on_west_edge,
int on_east_edge,
int on_south_edge,
int on_north_edge)
218 int istart, iend, jstart, jend;
226 for(i=0; i<nxp; i++) {
230 for(j=0; j<nyp; j++) {
235 px = (
double *)malloc(2*(nx+2)*
sizeof(double));
236 py = (
double *)malloc(2*(ny+2)*
sizeof(double));
258 for(j=jstart; j<=jend; j++) {
260 p1[0] = lont[j*(nx+2)+i ]; p1[1] = latt[j*(nx+2)+i ];
261 p2[0] = lont[j*(nx+2)+i+1]; p2[1] = latt[j*(nx+2)+i+1];
262 mid_pt_sphere(p1, p2, &(py[2*j]));
265 for(j=jstart; j<jend; j++) {
266 p1[0] = lonc[j*nxp+i];
267 p1[1] = latc[j*nxp+i];
268 d1 = great_circle_distance(py+2*j, p1);
269 d2 = great_circle_distance(py+2*(j+1), p1);
270 edge_w[j] = d2/(d1+d2);
276 for(j=jstart; j<=jend; j++) {
278 p1[0] = lont[j*(nx+2)+i ]; p1[1] = latt[j*(nx+2)+i ];
279 p2[0] = lont[j*(nx+2)+i+1]; p2[1] = latt[j*(nx+2)+i+1];
280 mid_pt_sphere(p1, p2, &(py[2*j]));
283 for(j=jstart; j<jend; j++) {
284 p1[0] = lonc[j*nxp+i];
285 p1[1] = latc[j*nxp+i];
286 d1 = great_circle_distance(&(py[2*j]), p1);
287 d2 = great_circle_distance(&(py[2*(j+1)]), p1);
288 edge_e[j] = d2/(d1+d2);
295 for(i=istart; i<=iend; i++) {
296 p1[0] = lont[j *(nx+2)+i]; p1[1] = latt[j *(nx+2)+i];
297 p2[0] = lont[(j+1)*(nx+2)+i]; p2[1] = latt[(j+1)*(nx+2)+i];
298 mid_pt_sphere(p1, p2, &(px[2*i]));
300 for(i=istart; i<iend; i++) {
301 p1[0] = lonc[j*nxp+i];
302 p1[1] = latc[j*nxp+i];
303 d1 = great_circle_distance(&(px[2*i]), p1);
304 d2 = great_circle_distance(&(px[2*(i+1)]), p1);
305 edge_s[i] = d2/(d1+d2);
311 for(i=istart; i<=iend; i++) {
312 p1[0] = lont[j *(nx+2)+i]; p1[1] = latt[j *(nx+2)+i];
313 p2[0] = lont[(j+1)*(nx+2)+i]; p2[1] = latt[(j+1)*(nx+2)+i];
314 mid_pt_sphere(p1, p2, &(px[2*i]));
316 for(i=istart; i<iend; i++) {
317 p1[0] = lonc[j*nxp+i];
318 p1[1] = latc[j*nxp+i];
319 d1 = great_circle_distance(&(px[2*i]), p1);
320 d2 = great_circle_distance(&(px[2*(i+1)]), p1);
321 edge_n[i] = d2/(d1+d2);
330 void mid_pt_sphere(
const double *p1,
const double *p2,
double *pm)
332 double e1[3], e2[3], e3[3];
336 mid_pt3_cart(e1, e2, e3);
337 xyz2latlon(1, e3, e3+1, e3+2, pm, pm+1);
341 void mid_pt3_cart(
const double *p1,
const double *p2,
double *e)
345 e[0] = p1[0] + p2[0];
346 e[1] = p1[1] + p2[1];
347 e[2] = p1[2] + p2[2];
348 dd = sqrt( e[0]*e[0] + e[1]*e[1] + e[2]*e[2] );
373 void calc_c2l_grid_info_(
int *nx_pt,
int *ny_pt,
const double *xt,
const double *yt,
const double *xc,
const double *yc,
374 double *dx,
double *dy,
double *
area,
double *edge_w,
double *edge_e,
double *edge_s,
375 double *edge_n,
double *en_n,
double *en_e,
double *vlon,
double *vlat,
376 int *on_west_edge,
int *on_east_edge,
int *on_south_edge,
int *on_north_edge)
378 calc_c2l_grid_info(nx_pt, ny_pt, xt, yt, xc, yc, dx, dy,
area, edge_w, edge_e, edge_s, edge_n,
379 en_n, en_e, vlon, vlat, on_west_edge, on_east_edge, on_south_edge, on_north_edge);
383 void calc_c2l_grid_info(
int *nx_pt,
int *ny_pt,
const double *xt,
const double *yt,
const double *xc,
const double *yc,
384 double *dx,
double *dy,
double *
area,
double *edge_w,
double *edge_e,
double *edge_s,
385 double *edge_n,
double *en_n,
double *en_e,
double *vlon,
double *vlat,
386 int *on_west_edge,
int *on_east_edge,
int *on_south_edge,
int *on_north_edge)
388 double *x, *y, *z, *xt_tmp, *yt_tmp;
389 int nx, ny, nxp, nyp, i, j;
390 double p1[3], p2[3], p3[3], p4[3];
398 for(j=0; j<nyp; j++)
for(i=0; i<nx; i++) {
401 p2[0] = xc[j*nxp+i+1];
402 p2[1] = yc[j*nxp+i+1];
403 dx[j*nx+i] = great_circle_distance(p1, p2);
406 for(j=0; j<ny; j++)
for(i=0; i<nxp; i++) {
409 p2[0] = xc[(j+1)*nxp+i];
410 p2[1] = yc[(j+1)*nxp+i];
411 dy[j*nxp+i] = great_circle_distance(p1, p2);
414 for(j=0; j<ny; j++)
for(i=0; i<nx; i++) {
417 p2[0] = xc[(j+1)*nxp+i];
418 p2[1] = yc[(j+1)*nxp+i];
419 p3[0] = xc[j*nxp+i+1];
420 p3[1] = yc[j*nxp+i+1];
421 p4[0] = xc[(j+1)*nxp+i+1];
422 p4[1] = yc[(j+1)*nxp+i+1];
423 area[j*nx+i] = spherical_excess_area(p1, p2, p3, p4, RADIUS);
426 x = (
double *)malloc(nxp*nyp*
sizeof(
double));
427 y = (
double *)malloc(nxp*nyp*
sizeof(
double));
428 z = (
double *)malloc(nxp*nyp*
sizeof(
double));
431 for(j=0; j<nyp; j++)
for(i=0; i<nx; i++) {
435 p2[0] = x[j*nxp+i+1];
436 p2[1] = y[j*nxp+i+1];
437 p2[2] = z[j*nxp+i+1];
438 vect_cross(p1, p2, en_n+3*(j*nx+i) );
439 normalize_vect(en_n+3*(j*nx+i));
442 for(j=0; j<ny; j++)
for(i=0; i<nxp; i++) {
446 p1[0] = x[(j+1)*nxp+i];
447 p1[1] = y[(j+1)*nxp+i];
448 p1[2] = z[(j+1)*nxp+i];
449 vect_cross(p1, p2, en_e+3*(j*nxp+i) );
450 normalize_vect(en_e+3*(j*nxp+i));
453 xt_tmp = (
double *)malloc(nx*ny*
sizeof(
double));
454 yt_tmp = (
double *)malloc(nx*ny*
sizeof(
double));
455 for(j=0; j<ny; j++)
for(i=0; i<nx; i++) {
456 xt_tmp[j*nx+i] = xt[(j+1)*(nx+2)+i+1];
457 yt_tmp[j*nx+i] = yt[(j+1)*(nx+2)+i+1];
459 unit_vect_latlon(nx*ny, xt_tmp, yt_tmp, vlon, vlat);
460 get_edge(nx, ny, xt, yt, xc, yc, edge_w, edge_e, edge_s, edge_n, *on_west_edge, *on_east_edge,
461 *on_south_edge, *on_north_edge);
pure elemental type(point) function latlon2xyz(lat, lon)
Return the (x,y,z) position of a given (lat,lon) point.
real(r8_kind), dimension(:,:), allocatable area
area of each grid box
integer nlat
No description.
integer nlon
No description.