mosaic/interp.c

   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 <stdio.h>
  20 #include <stdlib.h>
  21 #include <math.h>
  22 #include "mosaic_util.h"
  23 #include "interp.h"
  24 #include "create_xgrid.h"
  25
  26 /** \file
  27  *  \ingroup mosaic
  28  *  \brief Grid interpolation functions for use in @ref mosaic_mod
  29  */
  30
  31 /*********************************************************************
  32    void cublic_spline_sp(size1, size2, grid1, grid2, data1, data2)
  33
  34    Calculate a shape preserving cubic spline. Monotonicity is ensured over each subinterval
  35    unlike classic cubic spline interpolation.
  36    It will be used to interpolation data in 1-D space.
  37
  38    INPUT Arguments:
  39     grid1:    grid for input data grid.
  40     grid2:    grid for output data grid.
  41     size1:    size of input grid.
  42     size2:    size of output grid.
  43     data1:    input data associated with grid1.
  44
  45    OUTPUT ARGUMENTS:
  46     data2:    output data associated with grid2. (OUTPUT)
  47
  48 *********************************************************************/
  49
  50 void cubic_spline_sp(int size1, int size2, const double *grid1, const double *grid2, const double *data1,
  51                      double *data2 )
  52 {
  53   double *delta=NULL, *d=NULL, *dh=NULL, *b=NULL, *c = NULL;
  54   double s, w1, w2, p;
  55   int i, k, n, klo, khi, kmax;
  56
  57   for(i=1; i<size1; i++) {
  58     if( grid1[i] <= grid1[i-1] ) error_handler("cubic_spline_sp: grid1 is not monotonic increasing");
  59   }
  60
  61   for(i=0; i<size2; i++) {
  62     if( grid2[i] < grid1[0] || grid2[i] > grid1[size1-1]) error_handler("cubic_spline_sp: grid2 lies outside grid1");
  63   }
  64
  65   if(size1 < 2) error_handler("cubic_spline_sp: the size of input grid should be at least 2");
  66   if(size1 == 2) {  /* when size1 is 2, it just reduced to a linear interpolation */
  67     p = (data1[1]-data1[0])/(grid1[1]-grid1[0]);
  68     for(i=0; i< size2; i++) data2[i] = p*(grid2[i] - grid1[0]) + data1[0];
  69     return;
  70   }
  71   delta = (double *)malloc((size1-1)*sizeof(double));
  72   dh = (double *)malloc((size1-1)*sizeof(double));
  73   d = (double *)malloc(size1*sizeof(double));
  74   for(k=0;k<size1-1;k++) {
  75     dh[k] = grid1[k+1]-grid1[k];
  76     delta[k] = (data1[k+1]-data1[k])/(dh[k]);
  77   }
  78   /*
  79     Interior slopes
  80   */
  81   for(k=1;k<size1-1;k++) {
  82     if( delta[k]*delta[k-1] > 0.0 ) {
  83       w1 = 2.0*dh[k] + dh[k-1];
  84       w2 = dh[k] + 2.0*dh[k-1];
  85       d[k] = (w1+w2)/(w1/delta[k-1]+w2/delta[k]);
  86     }
  87     else {
  88       d[k] = 0.0;
  89     }
  90   }
  91   /*
  92     End slopes
  93   */
  94   kmax = size1-1;
  95   d[0] = ((2.0*dh[0] + dh[1])*delta[0] - dh[0]*delta[1])/(dh[0]+dh[1]);
  96
  97   if ( d[0]*delta[0] < 0.0 ) {
  98     d[0] = 0.0;
  99   }
 100   else {
 101     if ( delta[0]*delta[1] < 0.0 && fabs(d[0]) > fabs(3.0*delta[0])) {
 102       d[0]=3.0*delta[0];
 103     }
 104   }
 105
 106   d[kmax] = ((2.0*dh[kmax-1] + dh[kmax-2])*delta[kmax-1] - dh[kmax-1]*delta[kmax-2])/(dh[kmax-1]+dh[kmax-2]);
 107   if ( d[kmax]*delta[kmax-1] < 0.0 ) {
 108     d[kmax] = 0.0;
 109   }
 110   else {
 111     if ( delta[kmax-1]*delta[kmax-2] < 0.0 && fabs(d[kmax]) > fabs(3.0*delta[kmax-1])) {
 112       d[kmax]=3.0*delta[kmax-1];
 113     }
 114   }
 115
 116   /* Precalculate coefficients */
 117   b = (double *)malloc((size1-1)*sizeof(double));
 118   c = (double *)malloc((size1-1)*sizeof(double));
 119   for (k=0; k<size1-1; k++) {
 120     c[k]   = (3.0*delta[k]-2.0*d[k]-d[k+1])/dh[k];
 121     b[k]   = (d[k]-2.0*delta[k]+d[k+1])/(dh[k]*dh[k]);
 122   }
 123   /* interpolate data onto grid2 */
 124   for(k=0; k<size2; k++) {
 125     n = nearest_index(grid2[k],grid1, size1);
 126     if (grid1[n] < grid2[k]) {
 127       klo = n;
 128     }
 129     else {
 130       if(n==0) {
 131         klo = n;
 132       }
 133       else {
 134         klo = n -1;
 135       }
 136     }
 137     khi = klo+1;
 138     s   = grid2[k] - grid1[klo];
 139     data2[k] = data1[klo]+s*(d[klo]+s*(c[klo]+s*b[klo]));
 140   }
 141
 142   free(c);
 143   free(b);
 144   free(d);
 145   free(dh);
 146   free(delta);
 147
 148 }/* cubic spline sp */
 149
 150
 151 /*********************************************************************
 152    void cublic_spline(size1, size2, grid1, grid2, data1, data2, yp1, ypn  )
 153
 154    This alborithm is to get an interpolation formula that is smooth in the first
 155    derivative, and continuous in the second derivative, both within an interval
 156    and its boundaries. It will be used to interpolation data in 1-D space.
 157
 158    INPUT Arguments:
 159     grid1:    grid for input data grid.
 160     grid2:    grid for output data grid.
 161     size1:    size of input grid.
 162     size2:    size of output grid.
 163     data1:    input data associated with grid1.
 164     yp1:      first derivative of starting point.
 165               Set to 0 to be "natural"           (INPUT)
 166     ypn:      first derivative of ending point.
 167               Set to 0 to be "natural"           (INPUT)
 168
 169    OUTPUT ARGUMENTS:
 170     data2:    output data associated with grid2. (OUTPUT)
 171
 172 *********************************************************************/
 173
 174
 175 void cubic_spline(int size1, int size2, const double *grid1, const double *grid2, const double *data1,
 176                   double *data2, double yp1, double ypn  )
 177 {
 178   double *y2=NULL, *u=NULL;
 179   double p, sig, qn, un, h, a, b;
 180   int i, k, n, klo, khi;
 181
 182   for(i=1; i<size1; i++) {
 183     if( grid1[i] <= grid1[i-1] ) error_handler("cubic_spline: grid1 is not monotonic increasing");
 184   }
 185
 186   for(i=0; i<size2; i++) {
 187     if( grid2[i] < grid1[0] || grid2[i] > grid1[size1-1]) error_handler("cubic_spline: grid2 lies outside grid1");
 188   }
 189
 190   if(size1 < 2) error_handler("cubic_spline: the size of input grid should be at least 2");
 191   if(size1 == 2) {  /* when size1 is 2, it just reduced to a linear interpolation */
 192     p = (data1[1]-data1[0])/(grid1[1]-grid1[0]);
 193     for(i=0; i< size2; i++) data2[i] = p*(grid2[i] - grid1[0]) + data1[0];
 194     return;
 195   }
 196   y2 = (double *)malloc(size1*sizeof(double));
 197   u = (double *)malloc(size1*sizeof(double));
 198   if (yp1 >.99e30) {
 199     y2[0]=0.;
 200     u[0]=0.;
 201   }
 202   else {
 203     y2[0]=-0.5;
 204     u[0]=(3./(grid1[1]-grid1[0]))*((data1[1]-data1[0])/(grid1[1]-grid1[0])-yp1);
 205   }
 206
 207   for(i=1; i<size1-1; i++) {
 208     sig=(grid1[i]-grid1[i-1])/(grid1[i+1]-grid1[i-1]);
 209     p=sig*y2[i-1]+2.;
 210     y2[i]=(sig-1.)/p;
 211     u[i]=(6.*((data1[i+1]-data1[i])/(grid1[i+1]-grid1[i])-(data1[i]-data1[i-1])
 212               /(grid1[i]-grid1[i-1]))/(grid1[i+1]-grid1[i-1])-sig*u[i-1])/p;
 213
 214   }
 215
 216   if (ypn > .99e30) {
 217     qn=0.;
 218     un=0.;
 219   }
 220   else {
 221     qn=0.5;
 222     un=(3./(grid1[size1-1]-grid1[size1-2]))*(ypn-(data1[size1-1]-data1[size1-2])/(grid1[size1-1]-grid1[size1-2]));
 223   }
 224
 225   y2[size1-1]=(un-qn*u[size1-2])/(qn*y2[size1-2]+1.);
 226
 227   for(k=size1-2; k>=0; k--) y2[k] = y2[k]*y2[k+1]+u[k];
 228
 229   /* interpolate data onto grid2 */
 230   for(k=0; k<size2; k++) {
 231     n = nearest_index(grid2[k],grid1, size1);
 232     if (grid1[n] < grid2[k]) {
 233       klo = n;
 234     }
 235     else {
 236       if(n==0) {
 237         klo = n;
 238       }
 239       else {
 240         klo = n -1;
 241       }
 242     }
 243     khi = klo+1;
 244     h   = grid1[khi]-grid1[klo];
 245     a   = (grid1[khi] - grid2[k])/h;
 246     b   = (grid2[k] - grid1[klo])/h;
 247     data2[k] = a*data1[klo] + b*data1[khi]+ ((pow(a,3.0)-a)*y2[klo] + (pow(b,3.0)-b)*y2[khi])*(pow(h,2.0))/6;
 248   }
 249
 250   free(y2);
 251   free(u);
 252
 253 }/* cubic spline */
 254
 255
 256 /*------------------------------------------------------------------------------
 257   void conserve_interp()
 258   conservative interpolation through exchange grid.
 259   Currently only first order interpolation are implemented here.
 260   ----------------------------------------------------------------------------*/
 261 void conserve_interp(int nx_src, int ny_src, int nx_dst, int ny_dst, const double *x_src,
 262                      const double *y_src, const double *x_dst, const double *y_dst,
 263                      const double *mask_src, const double *data_src, double *data_dst )
 264 {
 265   int n, nxgrid;
 266   int *xgrid_i1, *xgrid_j1, *xgrid_i2, *xgrid_j2;
 267   double *xgrid_area, *dst_area, *area_frac;
 268
 269   /* get the exchange grid between source and destination grid. */
 270   xgrid_i1   = (int    *)malloc(MAXXGRID*sizeof(int));
 271   xgrid_j1   = (int    *)malloc(MAXXGRID*sizeof(int));
 272   xgrid_i2   = (int    *)malloc(MAXXGRID*sizeof(int));
 273   xgrid_j2   = (int    *)malloc(MAXXGRID*sizeof(int));
 274   xgrid_area = (double *)malloc(MAXXGRID*sizeof(double));
 275   dst_area   = (double *)malloc(nx_dst*ny_dst*sizeof(double));
 276   nxgrid = create_xgrid_2dx2d_order1(&nx_src, &ny_src, &nx_dst, &ny_dst, x_src, y_src, x_dst, y_dst, mask_src,
 277                                      xgrid_i1, xgrid_j1, xgrid_i2, xgrid_j2, xgrid_area );
 278   /* The source grid may not cover the destination grid
 279      so need to sum of exchange grid area to get dst_area
 280      get_grid_area(&nx_dst, &ny_dst, x_dst, y_dst, dst_area);
 281   */
 282   for(n=0; n<nx_dst*ny_dst; n++) dst_area[n] = 0;
 283   for(n=0; n<nxgrid; n++) dst_area[xgrid_j2[n]*nx_dst+xgrid_i2[n]] += xgrid_area[n];
 284
 285   area_frac = (double *)malloc(nxgrid*sizeof(double));
 286   for(n=0; n<nxgrid; n++) area_frac[n] = xgrid_area[n]/dst_area[xgrid_j2[n]*nx_dst+xgrid_i2[n]];
 287
 288   for(n=0; n<nx_dst*ny_dst; n++) {
 289     data_dst[n] = 0;
 290   }
 291   for(n=0; n<nxgrid; n++) {
 292     data_dst[xgrid_j2[n]*nx_dst+xgrid_i2[n]] += data_src[xgrid_j1[n]*nx_src+xgrid_i1[n]]*area_frac[n];
 293   }
 294
 295   free(xgrid_i1);
 296   free(xgrid_j1);
 297   free(xgrid_i2);
 298   free(xgrid_j2);
 299   free(xgrid_area);
 300   free(dst_area);
 301   free(area_frac);
 302
 303 } /* conserve_interp */
 304
 305 /*------------------------------------------------------------------------------
 306   void conserve_interp_great_circle()
 307   conservative interpolation through exchange grid.
 308   Currently only first order interpolation are implemented here.
 309   great_circle algorithm is used for clipping and interpolation.
 310   ----------------------------------------------------------------------------*/
 311 void conserve_interp_great_circle(int nx_src, int ny_src, int nx_dst, int ny_dst, const double *x_src,
 312                                   const double *y_src, const double *x_dst, const double *y_dst,
 313                                   const double *mask_src, const double *data_src, double *data_dst )
 314 {
 315   int n, nxgrid;
 316   int *xgrid_i1, *xgrid_j1, *xgrid_i2, *xgrid_j2;
 317   double *xgrid_area, *dst_area, *area_frac, *xgrid_di, *xgrid_dj;
 318
 319   /* get the exchange grid between source and destination grid. */
 320   xgrid_i1   = (int    *)malloc(MAXXGRID*sizeof(int));
 321   xgrid_j1   = (int    *)malloc(MAXXGRID*sizeof(int));
 322   xgrid_i2   = (int    *)malloc(MAXXGRID*sizeof(int));
 323   xgrid_j2   = (int    *)malloc(MAXXGRID*sizeof(int));
 324   xgrid_area = (double *)malloc(MAXXGRID*sizeof(double));
 325   xgrid_di   = (double *)malloc(MAXXGRID*sizeof(double));
 326   xgrid_dj   = (double *)malloc(MAXXGRID*sizeof(double));
 327   dst_area   = (double *)malloc(nx_dst*ny_dst*sizeof(double));
 328   nxgrid = create_xgrid_great_circle(&nx_src, &ny_src, &nx_dst, &ny_dst, x_src, y_src, x_dst, y_dst, mask_src,
 329                                      xgrid_i1, xgrid_j1, xgrid_i2, xgrid_j2, xgrid_area, xgrid_di, xgrid_dj );
 330   /* The source grid may not cover the destination grid
 331      so need to sum of exchange grid area to get dst_area
 332      get_grid_area(&nx_dst, &ny_dst, x_dst, y_dst, dst_area);
 333   */
 334   for(n=0; n<nx_dst*ny_dst; n++) dst_area[n] = 0;
 335   for(n=0; n<nxgrid; n++) dst_area[xgrid_j2[n]*nx_dst+xgrid_i2[n]] += xgrid_area[n];
 336
 337   area_frac = (double *)malloc(nxgrid*sizeof(double));
 338   for(n=0; n<nxgrid; n++) area_frac[n] = xgrid_area[n]/dst_area[xgrid_j2[n]*nx_dst+xgrid_i2[n]];
 339
 340   for(n=0; n<nx_dst*ny_dst; n++) {
 341     data_dst[n] = 0;
 342   }
 343   for(n=0; n<nxgrid; n++) {
 344     data_dst[xgrid_j2[n]*nx_dst+xgrid_i2[n]] += data_src[xgrid_j1[n]*nx_src+xgrid_i1[n]]*area_frac[n];
 345   }
 346
 347   free(xgrid_i1);
 348   free(xgrid_j1);
 349   free(xgrid_i2);
 350   free(xgrid_j2);
 351   free(xgrid_area);
 352   free(dst_area);
 353   free(area_frac);
 354
 355 } /* conserve_interp_great_circle */
 356
 357
 358
 359 void linear_vertical_interp(int nx, int ny, int nk1, int nk2, const double *grid1, const double *grid2,
 360                             double *data1, double *data2)
 361 {
 362   int n, k, l;
 363   double w;
 364
 365   for(k=1; k<nk1; k++) {
 366     if(grid1[k] <= grid1[k-1]) error_handler("interp.c: grid1 not monotonic");
 367   }
 368   for(k=1; k<nk2; k++) {
 369     if(grid2[k] <= grid2[k-1]) error_handler("interp.c: grid2 not monotonic");
 370   }
 371
 372   if (grid1[0] > grid2[0] ) error_handler("interp.c: grid2 lies outside grid1");
 373   if (grid1[nk1-1] < grid2[nk2-1] ) error_handler("interp.c: grid2 lies outside grid1");
 374
 375   for(k=0; k<nk2; k++) {
 376     n = nearest_index(grid2[k],grid1,nk1);
 377     if (grid1[n] < grid2[k]) {
 378       w = (grid2[k]-grid1[n])/(grid1[n+1]-grid1[n]);
 379       for(l=0; l<nx*ny; l++) {
 380         data2[k*nx*ny+l] = (1.-w)*data1[n*nx*ny+l] + w*data1[(n+1)*nx*ny+l];
 381       }
 382     }
 383     else {
 384       if(n==0)
 385         for(l=0;l<nx*ny;l++) data2[k*nx*ny+l] = data1[n*nx*ny+l];
 386       else {
 387         w = (grid2[k]-grid1[n-1])/(grid1[n]-grid1[n-1]);
 388         for(l=0; l<nx*ny; l++) {
 389           data2[k*nx*ny+l] = (1.-w)*data1[(n-1)*nx*ny+l] + w*data1[n*nx*ny+l];
 390         }
 391       }
 392     }
 393   }
 394 }