FMS 2025.01-dev
Flexible Modeling System
Loading...
Searching...
No Matches
gradient_c2l.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 <math.h>
20#include <stdlib.h>
21#include "constant.h"
22#include "grid_utils.h"
23#include "gradient_c2l.h"
24#include <stdio.h>
25
26/** \file
27 * \ingroup mosaic
28 * \ Grid utility functions for use in @ref mosaic_mod
29 */
30
31/*------------------------------------------------------------------------------
32 Routine to compute gradient terms for SCRIP:
33 SJL: Oct 5, 2007
34 NOTe: pin has halo size = 1.
35 the size of pin will be (nx+2,ny+2), T-cell center, with halo = 1
36 the size of dx will be (nx, ny+1), N-cell center
37 the size of dy will be (nx+1, ny), E-cell center
38 the size of area will be (nx, ny), T-cell center.
39 The size of edge_w will be (ny+1), C-cell center
40 The size of edge_e will be (ny+1), C-cell center
41 The size of edge_s will be (nx+1), C-cell center
42 The size of edge_n will be (nx+1), C-cell center
43 The size of en_n will be (nx, ny+1,3),N-cell center
44 The size of en_e will be (nx+1,ny,3), E-cell center
45 The size of vlon will be (nx, ny, 3) T-cell center
46 The size of vlat will be (nx, ny, 3), T-cell center
47 ----------------------------------------------------------------------------*/
48void grad_c2l_(const int *nlon, const int *nlat, const double *pin, const double *dx, const double *dy, const double *area,
49 const double *edge_w, const double *edge_e, const double *edge_s, const double *edge_n,
50 const double *en_n, const double *en_e, const double *vlon, const double *vlat,
51 double *grad_x, double *grad_y, const int *on_west_edge, const int *on_east_edge,
52 const int *on_south_edge, const int *on_north_edge)
53{
54 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,
55 on_west_edge, on_east_edge, on_south_edge, on_north_edge);
56}
57
58void grad_c2l(const int *nlon, const int *nlat, const double *pin, const double *dx, const double *dy, const double *area,
59 const double *edge_w, const double *edge_e, const double *edge_s, const double *edge_n,
60 const double *en_n, const double *en_e, const double *vlon, const double *vlat,
61 double *grad_x, double *grad_y, const int *on_west_edge, const int *on_east_edge,
62 const int *on_south_edge, const int *on_north_edge)
63{
64
65 double *pb, *pdx, *pdy, *grad3;
66 int nx, ny, nxp, nyp, i, j, m0, m1, n;
67
68 nx = *nlon;
69 ny = *nlat;
70 nxp = nx+1;
71 nyp = ny+1;
72 pb = (double *)malloc(nxp*nyp*sizeof(double));
73 pdx = (double *)malloc(3*nx*(ny+1)*sizeof(double));
74 pdy = (double *)malloc(3*(nx+1)*ny*sizeof(double));
75 grad3 = (double *)malloc(3*nx*ny*sizeof(double));
76 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
78 for(j=0; j<nyp; j++) for(i=0; i<nx; i++) {
79 m0 = j*nx+i;
80 m1 = j*nxp+i;
81 for(n=0; n<3; n++) {
82 pdx[3*m0+n] = 0.5*(pb[m1]+pb[m1+1])*dx[m0]*en_n[3*m0+n];
83 }
84 }
85
86 for(j=0; j<ny; j++) for(i=0; i<nxp; i++) {
87 m0 = j*nxp+i;
88 for(n=0; n<3; n++) {
89 pdy[3*m0+n] = 0.5*(pb[m0]+pb[m0+nxp])*dy[m0]*en_e[3*m0+n];
90 }
91 }
92
93 /* Compute 3D grad of the input scalar field by Green's theorem */
94 for(j=0; j<ny; j++) for(i=0; i<nx; i++) {
95 m0 = 3*(j*nx+i);
96 for(n=0; n<3; n++) {
97 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];
98 }
99 }
100
101 /* Compute inner product: V3 * grad (pe) */
102 for(j=0; j<ny; j++) for(i=0; i<nx; i++) {
103 m0 = j*nx+i;
104 m1 = 3*m0;
105 /* dq / d(Lamda)*/
106 grad_x[m0] = (vlon[m1]*grad3[m1] + vlon[m1+1]*grad3[m1+1] + vlon[m1+2]*grad3[m1+2])/area[m0];
107 grad_x[m0] *= RADIUS;
108 /* dq / d(theta) */
109 grad_y[m0] = (vlat[m1]*grad3[m1] + vlat[m1+1]*grad3[m1+1] + vlat[m1+2]*grad3[m1+2] )/area[m0];
110 grad_y[m0] *= RADIUS;
111 }
112
113 free(pb);
114 free(pdx);
115 free(pdy);
116 free(grad3);
117
118} /* grad_c2l */
119
120/*------------------------------------------------------------------------------
121 qin: A-grid field, size (nx+2, ny+2)
122 qout: B-grid field, size (nx+1, ny+1)
123 ----------------------------------------------------------------------------*/
124void a2b_ord2(int nx, int ny, const double *qin, const double *edge_w, const double *edge_e,
125 const double *edge_s, const double *edge_n, double *qout,
126 int on_west_edge, int on_east_edge, int on_south_edge, int on_north_edge)
127{
128 int nxp, nyp, i, j;
129 int istart, iend, jstart, jend;
130 double *q1, *q2;
131 const double r3 = 1./3.;
132
133 nxp = nx+1;
134 nyp = ny+1;
135 q1 = (double *)malloc((nx+2)*sizeof(double));
136 q2 = (double *)malloc((ny+2)*sizeof(double));
137
138
139 if(on_west_edge)
140 istart = 1;
141 else
142 istart = 0;
143 if(on_east_edge)
144 iend = nx;
145 else
146 iend = nxp;
147 if(on_south_edge)
148 jstart = 1;
149 else
150 jstart = 0;
151 if(on_north_edge)
152 jend = ny;
153 else
154 jend = nyp;
155
156 /* internal region ( 1: nx-1, 1:ny-1) */
157 for(j=jstart; j<jend; j++) for(i=istart; i<iend; i++) {
158 qout[j*nxp+i] = 0.25*(qin[j*(nx+2)+i] + qin[j*(nx+2)+i+1] +
159 qin[(j+1)*(nx+2)+i] + qin[(j+1)*(nx+2)+i+1] );
160 }
161
162 /* Fix the 4 Corners */
163 if(on_west_edge && on_south_edge)qout[ 0] = r3*(qin[1* (nx+2)+1 ]+qin[1* (nx+2) ]+qin[ 1]); /* sw_corner */
164 if(on_east_edge && on_south_edge)qout[ nx] = r3*(qin[1* (nx+2)+nx]+qin[ nx ]+qin[1* (nx+2)+nxp]); /* se_corner */
165 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 ]); /* ne_corner */
166 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 ]); /* nw_corner */
167
168 /* West Edges */
169 if(on_west_edge) {
170 for(j=jstart; j<=jend; j++){
171 q2[j] = 0.5*(qin[j*(nx+2)] + qin[j*(nx+2)+1]);
172 }
173 for(j=jstart; j<jend; j++){
174 qout[j*nxp] = edge_w[j]*q2[j] + (1-edge_w[j])*q2[j+1];
175 }
176 }
177
178 /* East Edges */
179 if(on_east_edge) {
180 for(j=jstart; j<=jend; j++){
181 q2[j] = 0.5*(qin[j*(nx+2)+nx] + qin[j*(nx+2)+nxp]);
182 }
183 for(j=jstart; j<jend; j++){
184 qout[j*nxp+nx] = edge_e[j]*q2[j] + (1-edge_e[j])*q2[j+1];
185 }
186 }
187
188 /* south edge */
189 if(on_south_edge) {
190 for(i=istart; i<=iend; i++){
191 q1[i] = 0.5*(qin[i] + qin[(nx+2)+i]);
192 }
193 for(i=istart; i<iend; i++){
194 qout[i] = edge_s[i]*q1[i] + (1 - edge_s[i])*q1[i+1];
195 }
196 }
197
198 /* north edge */
199 if(on_north_edge) {
200 for(i=istart; i<=iend; i++){
201 q1[i] = 0.5*(qin[ny*(nx+2)+i] + qin[nyp*(nx+2)+i]);
202 }
203 for(i=istart; i<iend; i++){
204 qout[ny*nxp+i] = edge_n[i]*q1[i] + (1 - edge_n[i])*q1[i+1];
205 }
206 }
207
208 free(q1);
209 free(q2);
210
211} /* a2b_ord2 */
212
213
214void get_edge(int nx, int ny, const double *lont, const double *latt,
215 const double *lonc, const double *latc, double *edge_w, double *edge_e, double *edge_s, double *edge_n,
216 int on_west_edge, int on_east_edge, int on_south_edge, int on_north_edge)
217{
218 int i, j, nxp, nyp;
219 int istart, iend, jstart, jend;
220 double p1[2], p2[2];
221 double *py, *px;
222 double d1, d2;
223
224 nxp = nx + 1;
225 nyp = ny + 1;
226
227 for(i=0; i<nxp; i++) {
228 edge_s[i] = 0.5; /* dummy value */
229 edge_n[i] = 0.5; /* dummy value */
230 }
231 for(j=0; j<nyp; j++) {
232 edge_w[j] = 0.5; /* dummy value */
233 edge_e[j] = 0.5; /* dummy value */
234 }
235
236 px = (double *)malloc(2*(nx+2)*sizeof(double));
237 py = (double *)malloc(2*(ny+2)*sizeof(double));
238
239 if(on_west_edge)
240 istart = 1;
241 else
242 istart = 0;
243 if(on_east_edge)
244 iend = nx;
245 else
246 iend = nxp;
247 if(on_south_edge)
248 jstart = 1;
249 else
250 jstart = 0;
251 if(on_north_edge)
252 jend = ny;
253 else
254 jend = nyp;
255 /* west edge */
256
257 if(on_west_edge) {
258 i=0;
259 for(j=jstart; j<=jend; j++) {
260 /* get mid point sphere */
261 p1[0] = lont[j*(nx+2)+i ]; p1[1] = latt[j*(nx+2)+i ];
262 p2[0] = lont[j*(nx+2)+i+1]; p2[1] = latt[j*(nx+2)+i+1];
263 mid_pt_sphere(p1, p2, &(py[2*j]));
264 }
265
266 for(j=jstart; j<jend; j++) {
267 p1[0] = lonc[j*nxp+i];
268 p1[1] = latc[j*nxp+i];
269 d1 = great_circle_distance(py+2*j, p1);
270 d2 = great_circle_distance(py+2*(j+1), p1);
271 edge_w[j] = d2/(d1+d2);
272 }
273 }
274 /* east edge */
275 if(on_east_edge) {
276 i=nx;
277 for(j=jstart; j<=jend; j++) {
278 /* get mid point sphere */
279 p1[0] = lont[j*(nx+2)+i ]; p1[1] = latt[j*(nx+2)+i ];
280 p2[0] = lont[j*(nx+2)+i+1]; p2[1] = latt[j*(nx+2)+i+1];
281 mid_pt_sphere(p1, p2, &(py[2*j]));
282 }
283
284 for(j=jstart; j<jend; j++) {
285 p1[0] = lonc[j*nxp+i];
286 p1[1] = latc[j*nxp+i];
287 d1 = great_circle_distance(&(py[2*j]), p1);
288 d2 = great_circle_distance(&(py[2*(j+1)]), p1);
289 edge_e[j] = d2/(d1+d2);
290 }
291 }
292
293 /* south edge */
294 if(on_south_edge) {
295 j=0;
296 for(i=istart; i<=iend; i++) {
297 p1[0] = lont[j *(nx+2)+i]; p1[1] = latt[j *(nx+2)+i];
298 p2[0] = lont[(j+1)*(nx+2)+i]; p2[1] = latt[(j+1)*(nx+2)+i];
299 mid_pt_sphere(p1, p2, &(px[2*i]));
300 }
301 for(i=istart; i<iend; i++) {
302 p1[0] = lonc[j*nxp+i];
303 p1[1] = latc[j*nxp+i];
304 d1 = great_circle_distance(&(px[2*i]), p1);
305 d2 = great_circle_distance(&(px[2*(i+1)]), p1);
306 edge_s[i] = d2/(d1+d2);
307 }
308 }
309 /* north edge */
310 if(on_north_edge) {
311 j=ny;
312 for(i=istart; i<=iend; i++) {
313 p1[0] = lont[j *(nx+2)+i]; p1[1] = latt[j *(nx+2)+i];
314 p2[0] = lont[(j+1)*(nx+2)+i]; p2[1] = latt[(j+1)*(nx+2)+i];
315 mid_pt_sphere(p1, p2, &(px[2*i]));
316 }
317 for(i=istart; i<iend; i++) {
318 p1[0] = lonc[j*nxp+i];
319 p1[1] = latc[j*nxp+i];
320 d1 = great_circle_distance(&(px[2*i]), p1);
321 d2 = great_circle_distance(&(px[2*(i+1)]), p1);
322 edge_n[i] = d2/(d1+d2);
323 }
324 }
325
326 free(px);
327 free(py);
328
329} /* get_edge */
330
331void mid_pt_sphere(const double *p1, const double *p2, double *pm)
332{
333 double e1[3], e2[3], e3[3];
334
335 latlon2xyz(1, p1, p1+1, e1, e1+1, e1+2);
336 latlon2xyz(1, p2, p2+1, e2, e2+1, e2+2);
337 mid_pt3_cart(e1, e2, e3);
338 xyz2latlon(1, e3, e3+1, e3+2, pm, pm+1);
339
340}
341
342void mid_pt3_cart(const double *p1, const double *p2, double *e)
343{
344 double dd;
345
346 e[0] = p1[0] + p2[0];
347 e[1] = p1[1] + p2[1];
348 e[2] = p1[2] + p2[2];
349 dd = sqrt( e[0]*e[0] + e[1]*e[1] + e[2]*e[2] );
350 e[0] /= dd;
351 e[1] /= dd;
352 e[2] /= dd;
353}
354
355/**********************************************************************************************
356 This routine is used to calculate grid information for second order conservative interpolation
357 from cubic grid to other grid
358 the size of xt will be (nx+2,ny+2), T-cell center, with halo = 1
359 the size of yt will be (nx+2,ny+2), T-cell center, with halo = 1
360 the size of xc will be (nx+1,ny+1), C-cell center
361 the size of yc will be (nx+1,ny+1), C-cell center
362 the size of dx will be (nx, ny+1), N-cell center
363 the size of dy will be (nx+1, ny), E-cell center
364 the size of area will be (nx, ny), T-cell center.
365 The size of edge_w will be (ny-1), C-cell center, without two end point
366 The size of edge_e will be (ny-1), C-cell center, without two end point
367 The size of edge_s will be (nx-1), C-cell center, without two end point
368 The size of edge_n will be (nx-1), C-cell center, without two end point
369 The size of en_n will be (nx, ny+1,3),N-cell center
370 The size of en_e will be (nx+1,ny,3), E-cell center
371 The size of vlon will be (nx, ny) T-cell center
372 The size of vlat will be (nx, ny), T-cell center
373**********************************************************************************************/
374void calc_c2l_grid_info_(int *nx_pt, int *ny_pt, const double *xt, const double *yt, const double *xc, const double *yc,
375 double *dx, double *dy, double *area, double *edge_w, double *edge_e, double *edge_s,
376 double *edge_n, double *en_n, double *en_e, double *vlon, double *vlat,
377 int *on_west_edge, int *on_east_edge, int *on_south_edge, int *on_north_edge)
378{
379 calc_c2l_grid_info(nx_pt, ny_pt, xt, yt, xc, yc, dx, dy, area, edge_w, edge_e, edge_s, edge_n,
380 en_n, en_e, vlon, vlat, on_west_edge, on_east_edge, on_south_edge, on_north_edge);
381
382}
383
384void calc_c2l_grid_info(int *nx_pt, int *ny_pt, const double *xt, const double *yt, const double *xc, const double *yc,
385 double *dx, double *dy, double *area, double *edge_w, double *edge_e, double *edge_s,
386 double *edge_n, double *en_n, double *en_e, double *vlon, double *vlat,
387 int *on_west_edge, int *on_east_edge, int *on_south_edge, int *on_north_edge)
388{
389 double *x, *y, *z, *xt_tmp, *yt_tmp;
390 int nx, ny, nxp, nyp, i, j;
391 double p1[3], p2[3], p3[3], p4[3];
392
393
394 nx = *nx_pt;
395 ny = *ny_pt;
396 nxp = nx+1;
397 nyp = ny+1;
398
399 for(j=0; j<nyp; j++) for(i=0; i<nx; i++) {
400 p1[0] = xc[j*nxp+i];
401 p1[1] = yc[j*nxp+i];
402 p2[0] = xc[j*nxp+i+1];
403 p2[1] = yc[j*nxp+i+1];
404 dx[j*nx+i] = great_circle_distance(p1, p2);
405 }
406
407 for(j=0; j<ny; j++) for(i=0; i<nxp; i++) {
408 p1[0] = xc[j*nxp+i];
409 p1[1] = yc[j*nxp+i];
410 p2[0] = xc[(j+1)*nxp+i];
411 p2[1] = yc[(j+1)*nxp+i];
412 dy[j*nxp+i] = great_circle_distance(p1, p2);
413 }
414
415 for(j=0; j<ny; j++) for(i=0; i<nx; i++) {
416 p1[0] = xc[j*nxp+i]; /* ll lon */
417 p1[1] = yc[j*nxp+i]; /* ll lat */
418 p2[0] = xc[(j+1)*nxp+i]; /* ul lon */
419 p2[1] = yc[(j+1)*nxp+i]; /* ul lat */
420 p3[0] = xc[j*nxp+i+1]; /* lr lon */
421 p3[1] = yc[j*nxp+i+1]; /* lr lat */
422 p4[0] = xc[(j+1)*nxp+i+1]; /* ur lon */
423 p4[1] = yc[(j+1)*nxp+i+1]; /* ur lat */
424 area[j*nx+i] = spherical_excess_area(p1, p2, p3, p4, RADIUS);
425 }
426
427 x = (double *)malloc(nxp*nyp*sizeof(double));
428 y = (double *)malloc(nxp*nyp*sizeof(double));
429 z = (double *)malloc(nxp*nyp*sizeof(double));
430
431 latlon2xyz(nxp*nyp, xc, yc, x, y, z);
432 for(j=0; j<nyp; j++) for(i=0; i<nx; i++) {
433 p1[0] = x[j*nxp+i];
434 p1[1] = y[j*nxp+i];
435 p1[2] = z[j*nxp+i];
436 p2[0] = x[j*nxp+i+1];
437 p2[1] = y[j*nxp+i+1];
438 p2[2] = z[j*nxp+i+1];
439 vect_cross(p1, p2, en_n+3*(j*nx+i) );
440 normalize_vect(en_n+3*(j*nx+i));
441 }
442
443 for(j=0; j<ny; j++) for(i=0; i<nxp; i++) {
444 p2[0] = x[j*nxp+i];
445 p2[1] = y[j*nxp+i];
446 p2[2] = z[j*nxp+i];
447 p1[0] = x[(j+1)*nxp+i];
448 p1[1] = y[(j+1)*nxp+i];
449 p1[2] = z[(j+1)*nxp+i];
450 vect_cross(p1, p2, en_e+3*(j*nxp+i) );
451 normalize_vect(en_e+3*(j*nxp+i));
452 }
453
454 xt_tmp = (double *)malloc(nx*ny*sizeof(double));
455 yt_tmp = (double *)malloc(nx*ny*sizeof(double));
456 for(j=0; j<ny; j++)for(i=0; i<nx; i++) {
457 xt_tmp[j*nx+i] = xt[(j+1)*(nx+2)+i+1];
458 yt_tmp[j*nx+i] = yt[(j+1)*(nx+2)+i+1];
459 }
460 unit_vect_latlon(nx*ny, xt_tmp, yt_tmp, vlon, vlat);
461 get_edge(nx, ny, xt, yt, xc, yc, edge_w, edge_e, edge_s, edge_n, *on_west_edge, *on_east_edge,
462 *on_south_edge, *on_north_edge);
463
464 free(x);
465 free(y);
466 free(z);
467 free(xt_tmp);
468 free(yt_tmp);
469
470}
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.