FMS  2025.04
Flexible Modeling System
gradient_c2l.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 <math.h>
19 #include <stdlib.h>
20 #include "constant.h"
21 #include "grid_utils.h"
22 #include "gradient_c2l.h"
23 #include <stdio.h>
24 
25 /** \file
26  * \ingroup mosaic
27  * \ Grid utility functions for use in @ref mosaic_mod
28  */
29 
30 /*------------------------------------------------------------------------------
31  Routine to compute gradient terms for SCRIP:
32  SJL: Oct 5, 2007
33  NOTe: pin has halo size = 1.
34  the size of pin will be (nx+2,ny+2), T-cell center, with halo = 1
35  the size of dx will be (nx, ny+1), N-cell center
36  the size of dy will be (nx+1, ny), E-cell center
37  the size of area will be (nx, ny), T-cell center.
38  The size of edge_w will be (ny+1), C-cell center
39  The size of edge_e will be (ny+1), C-cell center
40  The size of edge_s will be (nx+1), C-cell center
41  The size of edge_n will be (nx+1), C-cell center
42  The size of en_n will be (nx, ny+1,3),N-cell center
43  The size of en_e will be (nx+1,ny,3), E-cell center
44  The size of vlon will be (nx, ny, 3) T-cell center
45  The size of vlat will be (nx, ny, 3), T-cell center
46  ----------------------------------------------------------------------------*/
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)
52 {
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);
55 }
56 
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)
62 {
63 
64  double *pb, *pdx, *pdy, *grad3;
65  int nx, ny, nxp, nyp, i, j, m0, m1, n;
66 
67  nx = *nlon;
68  ny = *nlat;
69  nxp = nx+1;
70  nyp = ny+1;
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);
76 
77  for(j=0; j<nyp; j++) for(i=0; i<nx; i++) {
78  m0 = j*nx+i;
79  m1 = j*nxp+i;
80  for(n=0; n<3; n++) {
81  pdx[3*m0+n] = 0.5*(pb[m1]+pb[m1+1])*dx[m0]*en_n[3*m0+n];
82  }
83  }
84 
85  for(j=0; j<ny; j++) for(i=0; i<nxp; i++) {
86  m0 = j*nxp+i;
87  for(n=0; n<3; n++) {
88  pdy[3*m0+n] = 0.5*(pb[m0]+pb[m0+nxp])*dy[m0]*en_e[3*m0+n];
89  }
90  }
91 
92  /* Compute 3D grad of the input scalar field by Green's theorem */
93  for(j=0; j<ny; j++) for(i=0; i<nx; i++) {
94  m0 = 3*(j*nx+i);
95  for(n=0; n<3; n++) {
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];
97  }
98  }
99 
100  /* Compute inner product: V3 * grad (pe) */
101  for(j=0; j<ny; j++) for(i=0; i<nx; i++) {
102  m0 = j*nx+i;
103  m1 = 3*m0;
104  /* dq / d(Lamda)*/
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;
107  /* dq / d(theta) */
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;
110  }
111 
112  free(pb);
113  free(pdx);
114  free(pdy);
115  free(grad3);
116 
117 } /* grad_c2l */
118 
119 /*------------------------------------------------------------------------------
120  qin: A-grid field, size (nx+2, ny+2)
121  qout: B-grid field, size (nx+1, ny+1)
122  ----------------------------------------------------------------------------*/
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)
126 {
127  int nxp, nyp, i, j;
128  int istart, iend, jstart, jend;
129  double *q1, *q2;
130  const double r3 = 1./3.;
131 
132  nxp = nx+1;
133  nyp = ny+1;
134  q1 = (double *)malloc((nx+2)*sizeof(double));
135  q2 = (double *)malloc((ny+2)*sizeof(double));
136 
137 
138  if(on_west_edge)
139  istart = 1;
140  else
141  istart = 0;
142  if(on_east_edge)
143  iend = nx;
144  else
145  iend = nxp;
146  if(on_south_edge)
147  jstart = 1;
148  else
149  jstart = 0;
150  if(on_north_edge)
151  jend = ny;
152  else
153  jend = nyp;
154 
155  /* internal region ( 1: nx-1, 1:ny-1) */
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] );
159  }
160 
161  /* Fix the 4 Corners */
162  if(on_west_edge && on_south_edge)qout[ 0] = r3*(qin[1* (nx+2)+1 ]+qin[1* (nx+2) ]+qin[ 1]); /* sw_corner */
163  if(on_east_edge && on_south_edge)qout[ nx] = r3*(qin[1* (nx+2)+nx]+qin[ nx ]+qin[1* (nx+2)+nxp]); /* se_corner */
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 ]); /* ne_corner */
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 ]); /* nw_corner */
166 
167  /* West Edges */
168  if(on_west_edge) {
169  for(j=jstart; j<=jend; j++){
170  q2[j] = 0.5*(qin[j*(nx+2)] + qin[j*(nx+2)+1]);
171  }
172  for(j=jstart; j<jend; j++){
173  qout[j*nxp] = edge_w[j]*q2[j] + (1-edge_w[j])*q2[j+1];
174  }
175  }
176 
177  /* East Edges */
178  if(on_east_edge) {
179  for(j=jstart; j<=jend; j++){
180  q2[j] = 0.5*(qin[j*(nx+2)+nx] + qin[j*(nx+2)+nxp]);
181  }
182  for(j=jstart; j<jend; j++){
183  qout[j*nxp+nx] = edge_e[j]*q2[j] + (1-edge_e[j])*q2[j+1];
184  }
185  }
186 
187  /* south edge */
188  if(on_south_edge) {
189  for(i=istart; i<=iend; i++){
190  q1[i] = 0.5*(qin[i] + qin[(nx+2)+i]);
191  }
192  for(i=istart; i<iend; i++){
193  qout[i] = edge_s[i]*q1[i] + (1 - edge_s[i])*q1[i+1];
194  }
195  }
196 
197  /* north edge */
198  if(on_north_edge) {
199  for(i=istart; i<=iend; i++){
200  q1[i] = 0.5*(qin[ny*(nx+2)+i] + qin[nyp*(nx+2)+i]);
201  }
202  for(i=istart; i<iend; i++){
203  qout[ny*nxp+i] = edge_n[i]*q1[i] + (1 - edge_n[i])*q1[i+1];
204  }
205  }
206 
207  free(q1);
208  free(q2);
209 
210 } /* a2b_ord2 */
211 
212 
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)
216 {
217  int i, j, nxp, nyp;
218  int istart, iend, jstart, jend;
219  double p1[2], p2[2];
220  double *py, *px;
221  double d1, d2;
222 
223  nxp = nx + 1;
224  nyp = ny + 1;
225 
226  for(i=0; i<nxp; i++) {
227  edge_s[i] = 0.5; /* dummy value */
228  edge_n[i] = 0.5; /* dummy value */
229  }
230  for(j=0; j<nyp; j++) {
231  edge_w[j] = 0.5; /* dummy value */
232  edge_e[j] = 0.5; /* dummy value */
233  }
234 
235  px = (double *)malloc(2*(nx+2)*sizeof(double));
236  py = (double *)malloc(2*(ny+2)*sizeof(double));
237 
238  if(on_west_edge)
239  istart = 1;
240  else
241  istart = 0;
242  if(on_east_edge)
243  iend = nx;
244  else
245  iend = nxp;
246  if(on_south_edge)
247  jstart = 1;
248  else
249  jstart = 0;
250  if(on_north_edge)
251  jend = ny;
252  else
253  jend = nyp;
254  /* west edge */
255 
256  if(on_west_edge) {
257  i=0;
258  for(j=jstart; j<=jend; j++) {
259  /* get mid point sphere */
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]));
263  }
264 
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);
271  }
272  }
273  /* east edge */
274  if(on_east_edge) {
275  i=nx;
276  for(j=jstart; j<=jend; j++) {
277  /* get mid point sphere */
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]));
281  }
282 
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);
289  }
290  }
291 
292  /* south edge */
293  if(on_south_edge) {
294  j=0;
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]));
299  }
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);
306  }
307  }
308  /* north edge */
309  if(on_north_edge) {
310  j=ny;
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]));
315  }
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);
322  }
323  }
324 
325  free(px);
326  free(py);
327 
328 } /* get_edge */
329 
330 void mid_pt_sphere(const double *p1, const double *p2, double *pm)
331 {
332  double e1[3], e2[3], e3[3];
333 
334  latlon2xyz(1, p1, p1+1, e1, e1+1, e1+2);
335  latlon2xyz(1, p2, p2+1, e2, e2+1, e2+2);
336  mid_pt3_cart(e1, e2, e3);
337  xyz2latlon(1, e3, e3+1, e3+2, pm, pm+1);
338 
339 }
340 
341 void mid_pt3_cart(const double *p1, const double *p2, double *e)
342 {
343  double dd;
344 
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] );
349  e[0] /= dd;
350  e[1] /= dd;
351  e[2] /= dd;
352 }
353 
354 /**********************************************************************************************
355  This routine is used to calculate grid information for second order conservative interpolation
356  from cubic grid to other grid
357  the size of xt will be (nx+2,ny+2), T-cell center, with halo = 1
358  the size of yt will be (nx+2,ny+2), T-cell center, with halo = 1
359  the size of xc will be (nx+1,ny+1), C-cell center
360  the size of yc will be (nx+1,ny+1), C-cell center
361  the size of dx will be (nx, ny+1), N-cell center
362  the size of dy will be (nx+1, ny), E-cell center
363  the size of area will be (nx, ny), T-cell center.
364  The size of edge_w will be (ny-1), C-cell center, without two end point
365  The size of edge_e will be (ny-1), C-cell center, without two end point
366  The size of edge_s will be (nx-1), C-cell center, without two end point
367  The size of edge_n will be (nx-1), C-cell center, without two end point
368  The size of en_n will be (nx, ny+1,3),N-cell center
369  The size of en_e will be (nx+1,ny,3), E-cell center
370  The size of vlon will be (nx, ny) T-cell center
371  The size of vlat will be (nx, ny), T-cell center
372 **********************************************************************************************/
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)
377 {
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);
380 
381 }
382 
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)
387 {
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];
391 
392 
393  nx = *nx_pt;
394  ny = *ny_pt;
395  nxp = nx+1;
396  nyp = ny+1;
397 
398  for(j=0; j<nyp; j++) for(i=0; i<nx; i++) {
399  p1[0] = xc[j*nxp+i];
400  p1[1] = yc[j*nxp+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);
404  }
405 
406  for(j=0; j<ny; j++) for(i=0; i<nxp; i++) {
407  p1[0] = xc[j*nxp+i];
408  p1[1] = yc[j*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);
412  }
413 
414  for(j=0; j<ny; j++) for(i=0; i<nx; i++) {
415  p1[0] = xc[j*nxp+i]; /* ll lon */
416  p1[1] = yc[j*nxp+i]; /* ll lat */
417  p2[0] = xc[(j+1)*nxp+i]; /* ul lon */
418  p2[1] = yc[(j+1)*nxp+i]; /* ul lat */
419  p3[0] = xc[j*nxp+i+1]; /* lr lon */
420  p3[1] = yc[j*nxp+i+1]; /* lr lat */
421  p4[0] = xc[(j+1)*nxp+i+1]; /* ur lon */
422  p4[1] = yc[(j+1)*nxp+i+1]; /* ur lat */
423  area[j*nx+i] = spherical_excess_area(p1, p2, p3, p4, RADIUS);
424  }
425 
426  x = (double *)malloc(nxp*nyp*sizeof(double));
427  y = (double *)malloc(nxp*nyp*sizeof(double));
428  z = (double *)malloc(nxp*nyp*sizeof(double));
429 
430  latlon2xyz(nxp*nyp, xc, yc, x, y, z);
431  for(j=0; j<nyp; j++) for(i=0; i<nx; i++) {
432  p1[0] = x[j*nxp+i];
433  p1[1] = y[j*nxp+i];
434  p1[2] = z[j*nxp+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));
440  }
441 
442  for(j=0; j<ny; j++) for(i=0; i<nxp; i++) {
443  p2[0] = x[j*nxp+i];
444  p2[1] = y[j*nxp+i];
445  p2[2] = z[j*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));
451  }
452 
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];
458  }
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);
462 
463  free(x);
464  free(y);
465  free(z);
466  free(xt_tmp);
467  free(yt_tmp);
468 
469 }
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
real(r8_kind), dimension(:,:), allocatable area
area of each grid box
integer nlat
No description.
integer nlon
No description.