/*

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "globals.h"
#include "mosaic_util.h"
#include "bilinear_interp.h"
#include "mpp_io.h"
#include "mpp.h"

#define min(a,b) (a<b ? a:b)
#define max(a,b) (a>b ? a:b)
#define sign(a,b)(b<0 ? -fabs(a):fabs(a))


int max_weight_index( double *var, int nvar);
double normalize_great_circle_distance(const double *v1, const double *v2);
/*double spherical_angle(double *v1, double *v2, double *v3);*/
double dist2side(const double *v1, const double *v2, const double *point);
void redu2x(const double *varfin, const double *yfin, int nxfin, int nyfin, double *varcrs,
	    int nxcrs, int nycrs, int nz, int has_missing, double missvalue);
void do_latlon_coarsening(const double *var_latlon, const double *ylat, int nlon, int nlat, int nz,    
			  double *var_latlon_crs, int finer_steps, int has_missing, double missvalue);
void do_c2l_interp(const Interp_config *interp, int nx_in, int ny_in, int nz, const Field_config *field_in,
		   int nx_out, int ny_out, double *data_out, int has_missing, double missing, int fill_missing );
void sort_index(int ntiles, int *index, double *shortest);
int get_index(const Grid_config *grid_in, const Grid_config *grid_out, int *index,
	       int i_in, int j_in, int l_in, int i_out, int j_out);
int get_closest_index(const Grid_config *grid_in, const Grid_config *grid_out, int *index,
		      int i_in, int j_in, int l_in, int i_out, int j_out);

/*******************************************************************************
  void setup_bilinear_interp( )
    !------------------------------------------------------------------!
    ! calculate weights for bilinear interpolation                     !
    ! from cubed sphere to latlon grid                                 !
    !                                                                  !
    ! input:                                                           !
    ! sph_corner      cubed sphere corner location in spherical coor   !
    ! npx, npy        number of corners per tile                       !
    ! ntiles          number of tiles                                  !
    ! xlon, ylat      latlon grid coor                                 !
    ! nlon, nlat      latlon grid dimension                            !
    !                                                                  !
    ! output:                                                          !
    ! c2l_index       cubed sphere index for latlon interpolation      !
    ! c2l_weight      weights for cubsph_to_latlon interpolation       !
    ! elon_cubsph     lon unit vector for cubed sphere center          !
    ! elat_cubsph     lat unit vector for cubed sphere center          !
    ! elon_latlon     lon unit vector for latlon grid                  !
    ! elat_latlon     lat unit vector for latlon grid                  !
    !------------------------------------------------------------------!
*******************************************************************************/
void setup_bilinear_interp(int ntiles_in, const Grid_config *grid_in, int ntiles_out, const Grid_config *grid_out, 
                      Interp_config *interp, unsigned int opcode)
{
  const int max_iter = 10;
  double abs_center, dcub, dlon, dlat, coslat, distance;
  double dist1, dist2, dist3, dist4, sum;
  double *shortest;
  int    i, j, n, l, ic, jc, lc, icc, jcc, i_min, i_max, j_min, j_max, iter;
  int    n0, n1, n2, n3, n4, m0, m1;
  int    nx_in, ny_in, nx_out, ny_out, nxd, nyd;
  double v0[3], v1[3], v2[3], v3[3], v4[3];
  int    all_done;
  int    *found, *index;
  
  /* ntiles_in must be six and ntiles_out must be one */
  if(ntiles_in != 6) mpp_error("Error from bilinear_interp: source mosaic should be cubic mosaic "
			       "and have six tiles when using bilinear option");
  if(ntiles_out != 1) mpp_error("Error from bilinear_interp: destination mosaic should be "
				"one tile lat-lon grid when using bilinear option"); 
  /*-----------------------------------------------------------------!
    ! cubed sphere: cartesian coordinates of cell corners,             !
    !               cell lenghts between corners,                      !
    !               cartesian and spherical coordinates of cell centers! 
    !               calculate latlon unit vector                       !
    !-----------------------------------------------------------------*/

  /* calculation is done on the fine grid */
  nx_out   = grid_out->nx_fine;
  ny_out   = grid_out->ny_fine;
  nx_in    = grid_in->nx;       /* the cubic grid has same resolution on each face */
  ny_in    = grid_in->ny;
  nxd      = nx_in + 2;
  nyd      = ny_in + 2;

  interp->index  = (int    *)malloc(3*nx_out*ny_out*sizeof(int   ));
  interp->weight = (double *)malloc(4*nx_out*ny_out*sizeof(double));

  if( (opcode & READ) && interp->file_exist ) { /* reading from file */
    int nx2, ny2, fid, vid;

    /* check the size of the grid matching the size in remapping file */
    printf("NOTE: reading index and weight for bilinear interpolation from file.\n");
    fid = mpp_open(interp->remap_file, MPP_READ);
    nx2 = mpp_get_dimlen(fid, "nlon");
    ny2 = mpp_get_dimlen(fid, "nlat");
    printf("grid size is nx=%d, ny=%d, remap file size is nx=%d, ny=%d.\n", nx_out, ny_out, nx2, ny2);
    if(nx2 != nx_out || ny2 != ny_out ) mpp_error("bilinear_interp: size mismatch between grid size and remap file size");
    vid = mpp_get_varid(fid, "index");
    mpp_get_var_value(fid, vid, interp->index);
    vid = mpp_get_varid(fid, "weight");
    mpp_get_var_value(fid, vid, interp->weight);
    mpp_close(fid);
    return;
  }
  
  /*------------------------------------------------------------------
    find lower left corner on cubed sphere for given latlon location 
    ------------------------------------------------------------------*/
  found    = (int *)malloc(nx_out*ny_out*sizeof(int));
  index    = (int *)malloc(ntiles_in*3*sizeof(int));
  shortest = (double *)malloc(ntiles_in*sizeof(double));
  for(i=0; i<nx_out*ny_out; i++) found[i] = 0;
    
  dlon=(M_PI+M_PI)/nx_out;
  dlat=M_PI/(ny_out-1);

  for(iter=1; iter<=max_iter; iter++) {
    for(l=0; l<ntiles_in; l++) {
      for(jc=1; jc<=ny_in; jc++) for(ic=1; ic<=nx_in; ic++) {
	/*------------------------------------------------------
	  guess latlon indexes for given cubed sphere cell     
	  ------------------------------------------------------*/
	n1 = jc*nxd+ic;
	n2 = (jc+1)*nxd+ic+1;
	v1[0] = grid_in[l].xt[n1];
	v1[1] = grid_in[l].yt[n1];
	v1[2] = grid_in[l].zt[n1];
	v2[0] = grid_in[l].xt[n2];
	v2[1] = grid_in[l].yt[n2];
	v2[2] = grid_in[l].zt[n2];
	dcub=iter*normalize_great_circle_distance(v1, v2);
	j_min=max(   1,  floor((grid_in[l].latt[n1]-dcub+0.5*M_PI)/dlat)-iter+1);
	j_max=min(ny_out,ceil((grid_in[l].latt[n1]+dcub+0.5*M_PI)/dlat)+iter-1);
	if (j_min==1 || j_max==ny_out) {
	  i_min=1;
	  i_max=nx_out;
	}
        else {
	  i_min=max(   1,  floor((grid_in[l].lont[n1]-dcub)/dlon-iter+1));
	  i_max=min(nx_out,ceil((grid_in[l].lont[n1]+dcub)/dlon+iter-1));
	}
	for(j=j_min-1; j<j_max; j++) for(i=i_min-1; i<i_max; i++) {
	  n0 = j*nx_out + i;
	  /*--------------------------------------------------------------------
            for latlon cell find nearest cubed sphere cell 
            ------------------------------------------------------------------*/
	  if (!found[n0]) {
	    shortest[l]=M_PI+M_PI;
	    for(jcc=jc; jcc<=min(ny_in, jc+1); jcc++) for(icc=ic; icc<=min(nx_in, ic+1); icc++) {
	      n1 = jcc*nxd + icc;
	      n2 = j*nx_out + i;
	      v1[0] = grid_in[l].xt[n1];
	      v1[1] = grid_in[l].yt[n1];
	      v1[2] = grid_in[l].zt[n1];
	      v2[0] = grid_out->xt[n2];
	      v2[1] = grid_out->yt[n2];
	      v2[2] = grid_out->zt[n2];	      
	      distance=normalize_great_circle_distance(v1, v2);
	      if (distance < shortest[l]) {
		shortest[l]=distance;
		index[3*l]  =icc;
		index[3*l+1]=jcc;
		index[3*l+2]=l;
	      }
	    }
	    /*------------------------------------------------
              determine lower left corner                    
              ------------------------------------------------*/
	    found[n0] = get_closest_index(&(grid_in[l]), grid_out, &(interp->index[3*(j*nx_out+i)]), index[3*l],
			      index[3*l+1], index[3*l+2], i, j);
	  }
	}
      }
    }
    if (iter>1) {
      all_done = 1;
      for(i=0; i<nx_out*ny_out; i++) {
	if( !found[i]) {
	  all_done = 0;
	  break;
	}
      };
      if (all_done) break;
    }
  }
  
  /*------------------------------------------------------------------
    double check if lower left corner was found                      
    calculate weights for interpolation                              
    ------------------------------------------------------------------*/
  for(j=0; j<ny_out; j++) for(i=0; i<nx_out; i++) {
    n0 = j*nx_out + i;
    m0 = 3*n0;
    m1 = 4*n0;
    if (!found[n0]) {
      printf("**************************************************************\n");
      printf("WARNING: didn't find lower left corner for (ilon,jlat) = (%d,%d)\n", i,j);
      printf("will perform expensive global sweep\n");
      printf("**************************************************************\n");
      /*---------------------------------------------------------
	for latlon cell find nearby cubed sphere cell           
	---------------------------------------------------------*/
      for(l=0; l<3*ntiles_in; l++) index[l] = 0;	
      for(l=0; l<ntiles_in; l++) shortest[l] = M_PI + M_PI;
      for(l=0; l<ntiles_in; l++) {
	for(jc=0; jc<ny_in; j++) for(ic=0; ic<nx_in; ic++) {
	  n1 = jc*nxd+ic;
	  v1[0] = grid_in[l].xt[n1];
	  v1[1] = grid_in[l].yt[n1];
	  v1[2] = grid_in[l].zt[n1];
	  v0[0] = grid_out[l].xt[n0];
	  v0[1] = grid_out[l].yt[n0];
	  v0[2] = grid_out[l].zt[n0];	    
	  distance=normalize_great_circle_distance(v1, v2);
	  if (distance<shortest[l]) {
	    shortest[l]=distance;
	    index[3*l  ]=ic;
	    index[3*l+1]=jc;
	    index[3*l+2]=l;
	  }
	}
      }

      /*---------------------------------------------------------
	determine lower left corner                             
	---------------------------------------------------------*/
      sort_index(ntiles_in, index, shortest);
      found[n0]=0;
      for(l=0; l<ntiles_in; l++) {
	if (!found[n0]) {
	  found[n0]=get_index(&(grid_in[l]), grid_out, &(interp->index[3*(j*nx_out+i)]), index[3*l],
			      index[3*l+1], index[3*l+2], i, j);
	  if (found[n0]) break;
	}
      }
      if (! found[n0] ) mpp_error("error from bilinear_interp: couldn't find lower left corner");
    }
    /*------------------------------------------------------------
      calculate shortest distance to each side of rectangle      
      formed by cubed sphere cell centers                        
      special corner treatment                                   
      ------------------------------------------------------------*/
    ic=interp->index[m0];
    jc=interp->index[m0+1];
    l =interp->index[m0+2];
    if (ic==nx_in && jc==ny_in) {
      /*------------------------------------------------------------
	calculate weights for bilinear interpolation near corner   
	------------------------------------------------------------*/
      n1 = jc*nxd+ic;
      n2 = jc*nxd+ic+1;
      n3 = (jc+1)*nxd+ic;
      v1[0] = grid_in[l].xt[n1];
      v1[1] = grid_in[l].yt[n1];
      v1[2] = grid_in[l].zt[n1];
      v2[0] = grid_in[l].xt[n2];
      v2[1] = grid_in[l].yt[n2];
      v2[2] = grid_in[l].zt[n2];
      v3[0] = grid_in[l].xt[n3];
      v3[1] = grid_in[l].yt[n3];
      v3[2] = grid_in[l].zt[n3];
      v0[0] = grid_out->xt[n0];
      v0[1] = grid_out->yt[n0];
      v0[2] = grid_out->zt[n0];	  
      dist1=dist2side(v2, v3, v0);
      dist2=dist2side(v2, v1, v0);
      dist3=dist2side(v1, v3, v0);
      interp->weight[m1]  =dist1;      /* ic,   jc    weight */
      interp->weight[m1+1]=dist2;      /* ic,   jc+1  weight */
      interp->weight[m1+2]=0.;         /* ic+1, jc+1  weight */
      interp->weight[m1+3]=dist3;      /* ic+1, jc    weight */
             
      sum=interp->weight[m1]+interp->weight[m1+1]+interp->weight[m1+3];
      interp->weight[m1]  /=sum;
      interp->weight[m1+1]/=sum;
      interp->weight[m1+3]/=sum;
    }
    else if (ic==0 && jc==ny_in) {
      /*------------------------------------------------------------
	calculate weights for bilinear interpolation near corner   
	------------------------------------------------------------*/
	
      n1 = jc*nxd+ic;
      n2 = jc*nxd+ic+1;
      n3 = (jc+1)*nxd+ic+1;
      v1[0] = grid_in[l].xt[n1];
      v1[1] = grid_in[l].yt[n1];
      v1[2] = grid_in[l].zt[n1];
      v2[0] = grid_in[l].xt[n2];
      v2[1] = grid_in[l].yt[n2];
      v2[2] = grid_in[l].zt[n2];
      v3[0] = grid_in[l].xt[n3];
      v3[1] = grid_in[l].yt[n3];
      v3[2] = grid_in[l].zt[n3];
      v0[0] = grid_out->xt[n0];
      v0[1] = grid_out->yt[n0];
      v0[2] = grid_out->zt[n0];	  
      dist1=dist2side(v3, v2, v0);
      dist2=dist2side(v2, v1, v0);
      dist3=dist2side(v3, v1, v0);             
      interp->weight[m1]  =dist1;      /* ic,   jc    weight */
      interp->weight[m1+1]=0.;         /* ic,   jc+1  weight */
      interp->weight[m1+2]=dist2;      /* ic+1, jc+1  weight */
      interp->weight[m1+3]=dist3;      /* ic+1, jc    weight */
             
      sum=interp->weight[m1]+interp->weight[m1+2]+interp->weight[m1+3];
      interp->weight[m1]  /=sum;
      interp->weight[m1+2]/=sum;
      interp->weight[m1+3]/=sum;
    }
    else if (jc==0 && ic==nx_in) {
      /*------------------------------------------------------------
	calculate weights for bilinear interpolation near corner   
	------------------------------------------------------------*/
      n1 = jc*nxd+ic;
      n2 = (jc+1)*nxd+ic;
      n3 = (jc+1)*nxd+ic+1;
      v1[0] = grid_in[l].xt[n1];
      v1[1] = grid_in[l].yt[n1];
      v1[2] = grid_in[l].zt[n1];
      v2[0] = grid_in[l].xt[n2];
      v2[1] = grid_in[l].yt[n2];
      v2[2] = grid_in[l].zt[n2];
      v3[0] = grid_in[l].xt[n3];
      v3[1] = grid_in[l].yt[n3];
      v3[2] = grid_in[l].zt[n3];
      v0[0] = grid_out->xt[n0];
      v0[1] = grid_out->yt[n0];
      v0[2] = grid_out->zt[n0];	  
      dist1=dist2side(v2, v3, v0);
      dist2=dist2side(v1, v3, v0);
      dist3=dist2side(v1, v2, v0);             	     
             
      interp->weight[m1]  =dist1;      /* ic,   jc    weight */
      interp->weight[m1+1]=dist2;      /* ic,   jc+1  weight */
      interp->weight[m1+2]=dist3;      /* ic+1, jc+1  weight */
      interp->weight[m1+3]=0.;         /* ic+1, jc    weight */
             
      sum=interp->weight[m1]+interp->weight[m1+1]+interp->weight[m1+2];
      interp->weight[m1]  /=sum;
      interp->weight[m1+1]/=sum;
      interp->weight[m1+2]/=sum;
    }
    else {
      /*------------------------------------------------------------
	calculate weights for bilinear interpolation if no corner  
	------------------------------------------------------------*/
      n1 = jc*nxd+ic;
      n2 = jc*nxd+ic+1;
      n3 = (jc+1)*nxd+ic;
      n4 = (jc+1)*nxd+ic+1;
      v1[0] = grid_in[l].xt[n1];
      v1[1] = grid_in[l].yt[n1];
      v1[2] = grid_in[l].zt[n1];
      v2[0] = grid_in[l].xt[n2];
      v2[1] = grid_in[l].yt[n2];
      v2[2] = grid_in[l].zt[n2];
      v3[0] = grid_in[l].xt[n3];
      v3[1] = grid_in[l].yt[n3];
      v3[2] = grid_in[l].zt[n3];
      v4[0] = grid_in[l].xt[n4];
      v4[1] = grid_in[l].yt[n4];
      v4[2] = grid_in[l].zt[n4];	
      v0[0] = grid_out->xt[n0];
      v0[1] = grid_out->yt[n0];
      v0[2] = grid_out->zt[n0];	  
      dist1=dist2side(v1, v3, v0);
      dist2=dist2side(v3, v4, v0);
      dist3=dist2side(v4, v2, v0);
      dist4=dist2side(v2, v1, v0);
	
      interp->weight[m1]  =dist2*dist3;      /* ic,   jc    weight */
      interp->weight[m1+1]=dist3*dist4;      /* ic,   jc+1  weight */
      interp->weight[m1+2]=dist4*dist1;      /* ic+1, jc+1  weight */
      interp->weight[m1+3]=dist1*dist2;      /* ic+1, jc    weight */
             
      sum=interp->weight[m1]+interp->weight[m1+1]+interp->weight[m1+2]+interp->weight[m1+3];
      interp->weight[m1]  /=sum;
      interp->weight[m1+1]/=sum;
      interp->weight[m1+2]/=sum;
      interp->weight[m1+3]/=sum;
    }
  }

  /* write out weight information if needed */
  if( opcode & WRITE ) {
    int fid, dim_three, dim_four, dim_nlon, dim_nlat, dims[3];
    int fld_index, fld_weight;
    
    fid = mpp_open( interp->remap_file, MPP_WRITE);
    dim_nlon = mpp_def_dim(fid, "nlon", nx_out);
    dim_nlat = mpp_def_dim(fid, "nlat", ny_out);
    dim_three = mpp_def_dim(fid, "three", 3);
    dim_four  = mpp_def_dim(fid, "four", 4);
    
    dims[0] = dim_three; dims[1] = dim_nlat; dims[2] = dim_nlon;
    fld_index = mpp_def_var(fid, "index", NC_INT, 3, dims, 0);
    dims[0] = dim_four; dims[1] = dim_nlat; dims[2] = dim_nlon;
    fld_weight = mpp_def_var(fid, "weight", NC_DOUBLE, 3, dims, 0);
    mpp_end_def(fid);
    mpp_put_var_value(fid, fld_index, interp->index);
    mpp_put_var_value(fid, fld_weight, interp->weight);
    mpp_close(fid);
  }

  /* release the memory */
  free(found);
  free(shortest);
  free(index);
  
  printf("\n done calculating interp_index and interp_weight\n");
}; /* setup_bilinear_interp */

/*----------------------------------------------------------------------------
   void do_scalar_bilinear_interp(Mosaic_config *input, Mosaic_config *output, int varid )
   interpolate scalar data to latlon,                               !
   --------------------------------------------------------------------------*/
void do_scalar_bilinear_interp(const Interp_config *interp, int vid, int ntiles_in, const Grid_config *grid_in, const Grid_config *grid_out,
			  const Field_config *field_in, Field_config *field_out, int finer_step, int fill_missing)
{
  int    nx_in, ny_in, nx_out, ny_out, nz;
  int    n, ts, tn, tw, te;
  int    has_missing;
  double missing;
  double *data_fine;
  
  /*------------------------------------------------------------------
    determine target grid resolution                                 
    ------------------------------------------------------------------*/
  nx_out     = grid_out->nx_fine;
  ny_out     = grid_out->ny_fine;
  nx_in       = grid_in->nx;
  ny_in       = grid_in->ny;
  /* currently we are regridding one vertical level for each call to reduce the memory usage */
  nz = 1;  
  missing     = field_in[0].var[vid].missing;
  has_missing = field_in[0].var[vid].has_missing;

  data_fine = (double *)malloc(nx_out*ny_out*nz*sizeof(double));
  
  do_c2l_interp(interp, nx_in, ny_in, nz, field_in, nx_out, ny_out, data_fine, has_missing, missing, fill_missing);
  do_latlon_coarsening(data_fine, grid_out->latt1D_fine, nx_out, ny_out, nz, field_out->data,
		       finer_step, has_missing, missing);
  free(data_fine);
  
}; /* do_c2l_scalar_interp */


   
/*----------------------------------------------------------------------------
   void do_vector_bilinear_interp()
   interpolate vector data to latlon,                               !
   --------------------------------------------------------------------------*/
void do_vector_bilinear_interp(Interp_config *interp, int vid, int ntiles_in, const Grid_config *grid_in, int ntiles_out, 
			  const Grid_config *grid_out, const Field_config *u_in,  const Field_config *v_in,
			  Field_config *u_out, Field_config *v_out, int finer_step, int fill_missing)
{
  Field_config *var_cubsph;
  int          nx_in, ny_in, nx_out, ny_out, nxd, nyd, nz, has_missing;
  int          i, j, k, n, n1, n2, ts, tn, tw, te;
  double       missing;
  double       *x_latlon, *y_latlon, *z_latlon, *var_latlon;
  
  nx_out      = grid_out->nx_fine;
  ny_out      = grid_out->ny_fine;
  nx_in       = grid_in->nx;
  ny_in       = grid_in->ny;
  nxd         = nx_in + 2;
  nyd         = ny_in + 2;
  /* currently we are regridding one vertical level for each call to reduce the memory usage */
  nz = 1;  
  missing     = u_in[0].var[vid].missing;
  has_missing = u_in[0].var[vid].has_missing;

  
  x_latlon   = (double *)malloc(nx_out*ny_out*nz*sizeof(double));
  y_latlon   = (double *)malloc(nx_out*ny_out*nz*sizeof(double));
  z_latlon   = (double *)malloc(nx_out*ny_out*nz*sizeof(double));
  var_latlon = (double *)malloc(nx_out*ny_out*nz*sizeof(double));
  var_cubsph = (Field_config *)malloc(ntiles_in*sizeof(Field_config));
  for(n=0; n<ntiles_in; n++) var_cubsph[n].data = (double *)malloc((nx_in+2)*(ny_in+2)*nz*sizeof(double));

  for(n=0; n<ntiles_in; n++) {
    for(k=0; k<nz; k++) for(j=0; j<nyd; j++) for(i=0; i<nxd; i++) {
      n1 = k*nxd*nyd + j*nxd + i;
      n2 = j*nxd + i;
      var_cubsph[n].data[n1] = u_in[n].data[n1]*grid_in[n].vlon_t[3*n2]+v_in[n].data[n1]*grid_in[n].vlat_t[3*n2];
    }
  }

  do_c2l_interp(interp, nx_in, ny_in, nz, var_cubsph, nx_out, ny_out, x_latlon, has_missing, missing, fill_missing);  

  for(n=0; n<ntiles_in; n++) {
    for(k=0; k<nz; k++) for(j=0; j<nyd; j++) for(i=0; i<nxd; i++) {
      n1 = k*nxd*nyd + j*nxd + i;
      n2 = j*nxd + i;
      var_cubsph[n].data[n1] = u_in[n].data[n1]*grid_in[n].vlon_t[3*n2+1]+v_in[n].data[n1]*grid_in[n].vlat_t[3*n2+1];
    }
  }
  do_c2l_interp(interp, nx_in, ny_in, nz, var_cubsph, nx_out, ny_out, y_latlon, has_missing, missing, fill_missing);  
  
  for(n=0; n<ntiles_in; n++) {
    for(k=0; k<nz; k++) for(j=0; j<nyd; j++) for(i=0; i<nxd; i++) {
      n1 = k*nxd*nyd + j*nxd + i;
      n2 = j*nxd + i;
      var_cubsph[n].data[n1] = u_in[n].data[n1]*grid_in[n].vlon_t[3*n2+2]+v_in[n].data[n1]*grid_in[n].vlat_t[3*n2+2];
    }
  }

  do_c2l_interp(interp, nx_in, ny_in, nz, var_cubsph, nx_out, ny_out, z_latlon, has_missing, missing, fill_missing);  

  for(n=0; n<ntiles_in; n++) free(var_cubsph[n].data);
  free(var_cubsph);

  
  for(k=0; k<nz; k++) for(j=0; j<ny_out; j++) for(i=0; i<nx_out; i++) {
    n1 = k*nx_out*ny_out + j*nx_out + i;
    n2 = j*nx_out + i;
    var_latlon[n1] = x_latlon[n1]*grid_out->vlon_t[3*n2] +  y_latlon[n1]*grid_out->vlon_t[3*n2+1] +  z_latlon[n1]*grid_out->vlon_t[3*n2+2];
  }
  do_latlon_coarsening(var_latlon, grid_out->latt1D_fine, nx_out, ny_out, nz, u_out->data,
		       finer_step, has_missing, missing);

  for(k=0; k<nz; k++) for(j=0; j<ny_out; j++) for(i=0; i<nx_out; i++) {
    n1 = k*nx_out*ny_out + j*nx_out + i;
    n2 = j*nx_out + i;
    var_latlon[n1] = x_latlon[n1]*grid_out->vlat_t[3*n2] +  y_latlon[n1]*grid_out->vlat_t[3*n2+1] +  z_latlon[n1]*grid_out->vlat_t[3*n2+2];
  }

  do_latlon_coarsening(var_latlon, grid_out->latt1D_fine, nx_out, ny_out, nz, v_out->data,
		       finer_step, has_missing, missing);

  free(x_latlon);
  free(y_latlon);
  free(z_latlon);
    
}; /* do_vector_bilinear_interp */


void do_c2l_interp(const Interp_config *interp, int nx_in, int ny_in, int nz, const Field_config *field_in,
		   int nx_out, int ny_out, double *data_out, int has_missing, double missing, int fill_missing )
{
  int i, j, k, nxd, nyd, ic, jc, ind, n1, tile;
  double d_in[4];
  
  nxd = nx_in + 2;
  nyd = ny_in + 2;
  
  if (has_missing) {
    for(k=0; k<nz; k++) for(j=0; j<ny_out; j++) for(i=0; i<nx_out; i++) {
      n1      = j*nx_out+i;
      ic      = interp->index[3*n1];
      jc      = interp->index[3*n1+1];
      tile    = interp->index[3*n1+2];
      d_in[0] = field_in[tile].data[k*nxd*nyd+jc    *nxd+ic];
      d_in[1] = field_in[tile].data[k*nxd*nyd+(jc+1)*nxd+ic];
      d_in[2] = field_in[tile].data[k*nxd*nyd+(jc+1)*nxd+ic+1];
      d_in[3] = field_in[tile].data[k*nxd*nyd+jc    *nxd+ic+1];
      if (d_in[0] == missing || d_in[1] == missing || d_in[2] == missing || d_in[3] == missing ) {
	if (fill_missing) {
	  ind = max_weight_index( &(interp->weight[4*n1]), 4);
	  data_out[k*nx_out*ny_out+n1] = d_in[ind];
	}
	else {
	  data_out[k*nx_out*ny_out+n1] = missing;
	}
      }
      else {
	data_out[k*nx_out*ny_out+n1] = d_in[0]*interp->weight[4*n1] + d_in[1]*interp->weight[4*n1+1]
	  + d_in[2]*interp->weight[4*n1+2] + d_in[3]*interp->weight[4*n1+3];
      }
    }
  }
  else {
    for(k=0; k<nz; k++) for(j=0; j<ny_out; j++) for(i=0; i<nx_out; i++) {
      n1      = j*nx_out+i;
      ic      = interp->index[3*n1];
      jc      = interp->index[3*n1+1];
      tile    = interp->index[3*n1+2];
      d_in[0] = field_in[tile].data[k*nxd*nyd+jc    *nxd+ic];
      d_in[1] = field_in[tile].data[k*nxd*nyd+(jc+1)*nxd+ic];
      d_in[2] = field_in[tile].data[k*nxd*nyd+(jc+1)*nxd+ic+1];
      d_in[3] = field_in[tile].data[k*nxd*nyd+jc    *nxd+ic+1];
      data_out[k*nx_out*ny_out+n1] = d_in[0]*interp->weight[4*n1] + d_in[1]*interp->weight[4*n1+1]
	+ d_in[2]*interp->weight[4*n1+2] + d_in[3]*interp->weight[4*n1+3];
    }
  }
  
}; /* do_c2l_interp */

   
/*------------------------------------------------------------------
  void sort_index()
  sort index by shortest                                         
  ----------------------------------------------------------------*/
void sort_index(int ntiles, int *index, double *shortest)
{
  int l, ll, lll, i;
  double *shortest_sort;
  int    *index_sort;

  shortest_sort = (double *)malloc(3*ntiles*sizeof(double));
  index_sort    = (int    *)malloc(  ntiles*sizeof(int   ));
  
  for(l=0; l<3*ntiles; l++)index_sort[l] = 0;
  for(l=0; l<ntiles; l++)shortest_sort[l] = M_PI+M_PI;
  for(l=0; l<ntiles; l++) {
    for(ll=0; ll<ntiles; ll++) {
      if (shortest[l]<shortest_sort[ll]) {
	for(lll=ntiles-2; lll>=ll; lll--) {
	  index_sort[3*lll+1]=index_sort[3*lll];
	  shortest_sort[lll+1]=shortest_sort[lll];
	}
	for(i=0; i<3; i++) index_sort[3*ll+i]=index[3*l+i];
	shortest_sort[ll]=shortest[l];
	break;
      }
    }
  }

  for(l=0; l<3*ntiles; l++) index[l]    = index_sort[l];
  for(l=0; l<  ntiles; l++) shortest[l] = shortest_sort[l];

  free(shortest_sort);
  free(index_sort);
  
}; /* sort_index */


/*------------------------------------------------------------------
  void get_index(ig, jg, lg)
  determine lower left corner                                    
  ----------------------------------------------------------------*/
int get_index(const Grid_config *grid_in, const Grid_config *grid_out, int *index,
	       int i_in, int j_in, int l_in, int i_out, int j_out)
{
  int    ok, n0, n1, n2, n3, n4, n5;
  int    nx_in, ny_in, nx_out, ny_out;
  double v0[3], v1[3], v2[3], v3[3], v4[3], v5[3];
  double angle_1, angle_1a, angle_1b, angle_2, angle_2a, angle_2b;
  double angle_3, angle_3a, angle_3b, angle_4, angle_4a, angle_4b;
  
  ok=1;
  nx_in  = grid_in->nx_fine;
  ny_in  = grid_in->nx_fine;
  nx_out = grid_out->nx;
  ny_out = grid_out->nx;
  n0 = j_out*nx_out + i_out;
  n1 = j_in*nx_in + i_in;
  n2 = j_in*nx_in + i_in+1;
  n3 = (j_in+1)*nx_in + i_in;
  v0[0] = grid_out->xt[n1];
  v0[1] = grid_out->yt[n1];
  v0[2] = grid_out->zt[n1];
  v1[0] = grid_in->xt[n1];
  v1[1] = grid_in->yt[n1];
  v1[2] = grid_in->zt[n1];
  v2[0] = grid_in->xt[n2];
  v2[1] = grid_in->yt[n2];
  v2[2] = grid_in->zt[n2];
  v3[0] = grid_in->xt[n3];
  v3[1] = grid_in->yt[n3];
  v3[2] = grid_in->zt[n3];
  angle_1 = spherical_angle(v1, v2, v3);
  angle_1a= spherical_angle(v1, v2, v0); 
  angle_1b= spherical_angle(v1, v3, v0);
  
  if (max(angle_1a,angle_1b)<angle_1) {
    index[0]=i_in;
    index[1]=j_in;
    index[2]=l_in;
  }
  else {
    n4 = j_in*nx_in + i_in-1;
    v4[0] = grid_in->xt[n4];
    v4[1] = grid_in->yt[n4];
    v4[2] = grid_in->zt[n4];
    angle_2 =spherical_angle(v1, v3, v4);
    angle_2a=angle_1b;
    angle_2b=spherical_angle(v1, v4, v0);
    if (max(angle_2a,angle_2b)<angle_2) {
      index[0]=i_in-1;
      index[1]=j_in;
      index[2]=l_in;
    }
    else {
      n5 = (j_in-1)*nx_in + i_in;
      v5[0] = grid_in->xt[n5];
      v5[1] = grid_in->yt[n5];
      v5[2] = grid_in->zt[n5];      
      angle_3 =spherical_angle(v1, v4, v5);
      angle_3a=angle_2b;
      angle_3b=spherical_angle(v1, v5, v0);
      if (max(angle_3a,angle_3b)<angle_3 && i_in>1 && j_in>1) {
	index[0]=i_in-1;
	index[1]=j_in-1;
	index[2]=l_in;
      }
      else {
	angle_4 =spherical_angle(v1, v5, v2);
	angle_4a=angle_3b;
	angle_4b=spherical_angle(v1, v2, v0);
	if (max(angle_4a,angle_4b)<angle_4) {
	  index[0]=i_in;
	  index[1]=j_in-1;
	  index[2]=l_in;
	}
        else {
	  ok=0;
	}
      }
    }
  }
  return ok;
  
}; /* get_index */


/*-------------------------------------------------------------
  determine lower left corner
  void get_closest_index(int ig, int jg, int lg, int ok)
  --------------------------------------------------------------*/
int get_closest_index(const Grid_config *grid_in, const Grid_config *grid_out, int *index,
		      int i_in, int j_in, int l_in, int i_out, int j_out)
{
  int found;
  double angle_11, angle_11a, angle_11b;
  double angle_22, angle_22a, angle_22b;
  double angle_33, angle_33a, angle_33b;
  double angle_44, angle_44a, angle_44b;
  double angle_1,  angle_1a,  angle_1b;
  double angle_2,  angle_2a,  angle_2b;
  double angle_3,  angle_3a,  angle_3b;
  double angle_4,  angle_4a,  angle_4b;  
  int    n0, n1, n2, n3, n4, n5, n6, n7, n8;
  int    nx_in, ny_in, nx_out, ny_out, nxd;
  double v0[3], v1[3], v2[3], v3[3], v4[3], v5[3], v6[3], v7[3], v8[3];

  found = 0;
  nx_in  = grid_in->nx;
  ny_in  = grid_in->ny;
  nxd    = nx_in + 2;
  nx_out = grid_out->nx_fine;
  ny_out = grid_out->ny_fine;
  n0     = j_out*nx_out+i_out;
  n1     = j_in*nxd+i_in;
  n2     = j_in*nxd+i_in+1;
  n3     = (j_in+1)*nxd+i_in;
  v1[0]  = grid_in->xt[n1];
  v1[1]  = grid_in->yt[n1];
  v1[2]  = grid_in->zt[n1];
  v2[0]  = grid_in->xt[n2];
  v2[1]  = grid_in->yt[n2];
  v2[2]  = grid_in->zt[n2];
  v3[0]  = grid_in->xt[n3];
  v3[1]  = grid_in->yt[n3];
  v3[2]  = grid_in->zt[n3];
  v0[0]  = grid_out->xt[n0];
  v0[1]  = grid_out->yt[n0];
  v0[2]  = grid_out->zt[n0];
  angle_1 =spherical_angle(v1, v2, v3);
  angle_1a=spherical_angle(v1, v2, v0);
  angle_1b=spherical_angle(v1, v3, v0);
  if (max(angle_1a,angle_1b) <= angle_1) {
    if (i_in==nx_in && j_in==ny_in) {
      angle_11 =spherical_angle(v2, v3, v1);
      angle_11a=spherical_angle(v2, v1, v0);
      angle_11b=spherical_angle(v2, v3, v0);
    }
    else {
      n4 = (j_in+1)*nxd+i_in+1;
      v4[0]  = grid_in->xt[n4];
      v4[1]  = grid_in->yt[n4];
      v4[2]  = grid_in->zt[n4];
      angle_11 =spherical_angle(v4, v3, v2);
      angle_11a=spherical_angle(v4, v2, v0);
      angle_11b=spherical_angle(v4, v3, v0);
    }
    if (max(angle_11a,angle_11b)<=angle_11) {
      found = 1;
      index[0]=i_in;
      index[1]=j_in;
      index[2]=l_in;
    }
  }
  else {
    n4 = j_in*nxd+i_in-1;
    v4[0]  = grid_in->xt[n4];
    v4[1]  = grid_in->yt[n4];
    v4[2]  = grid_in->zt[n4]; 
    angle_2 =spherical_angle(v1,v3,v4);
    angle_2a=angle_1b;
    angle_2b=spherical_angle(v1,v4,v0);
    if (max(angle_2a,angle_2b)<=angle_2) {
      if (i_in==1 && j_in==ny_in) {
	angle_22 =spherical_angle(v3, v1, v4);
	angle_22a=spherical_angle(v3, v4, v0);
	angle_22b=spherical_angle(v3, v1, v0);
      }
      else {
	n5 = (j_in+1)*nxd+i_in-1;
	n6 = j_in    *nxd+i_in-1;
	v5[0]  = grid_in->xt[n5];
	v5[1]  = grid_in->yt[n5];
	v5[2]  = grid_in->zt[n5];
	v6[0]  = grid_in->xt[n6];
	v6[1]  = grid_in->yt[n6];
	v6[2]  = grid_in->zt[n6];
	angle_22 =spherical_angle(v5, v3, v6);
	angle_22a=spherical_angle(v5, v6, v0);
	angle_22b=spherical_angle(v5, v3, v0);
      }
      if (max(angle_22a,angle_22b)<=angle_22) {
	found=1;
	index[0]=i_in-1;
	index[1]=j_in;
	index[2]=l_in;
      }
    }
    else {
      n5 = j_in*nxd+i_in-1;
      n6 = (j_in-1)*nxd+i_in;
      v5[0]  = grid_in->xt[n5];
      v5[1]  = grid_in->yt[n5];
      v5[2]  = grid_in->zt[n5];
      v6[0]  = grid_in->xt[n6];
      v6[1]  = grid_in->yt[n6];
      v6[2]  = grid_in->zt[n6];	  
      angle_3 =spherical_angle(v1, v5, v6);
      angle_3a=angle_2b;
      angle_3b=spherical_angle(v1, v6, v0);
      if (max(angle_3a,angle_3b)<=angle_3 && i_in>1 && j_in>1) {
	n7 = (j_in-1)*nxd+i_in-1;
	v7[0]  = grid_in->xt[n7];
	v7[1]  = grid_in->yt[n7];
	v7[2]  = grid_in->zt[n7];    
	angle_33 =spherical_angle(v7, v6, v5);
	angle_33a=spherical_angle(v7, v5, v0);
	angle_33b=spherical_angle(v7, v6, v0);
	if (max(angle_33a,angle_33b)<=angle_33) {
	  found=1;
	  index[0]=i_in-1;
	  index[1]=j_in-1;
	  index[2]=l_in;
	}
      }
      else {
	angle_4 =spherical_angle(v1, v6, v2);
	angle_4a=angle_3b;
	angle_4b=spherical_angle(v1, v2, v0);
	if (max(angle_4a,angle_4b)<=angle_4) {
	  if (i_in==nx_in && j_in==1) {
	    angle_44 =spherical_angle(v2, v1, v6);
	    angle_44a=spherical_angle(v2, v6, v0);
	    angle_44b=spherical_angle(v2, v1, v0);
	  }
	  else {
	    n8 = (j_in-1)*nxd+i_in+1;
	    v8[0]  = grid_in->xt[n8];
	    v8[1]  = grid_in->yt[n8];
	    v8[2]  = grid_in->zt[n8];   
	    angle_44 =spherical_angle(v8, v2, v6);
	    angle_44a=spherical_angle(v8, v6, v0);
	    angle_44b=spherical_angle(v8, v2, v0);
	  }
	  if (max(angle_44a,angle_44b)<=angle_44) {
	    found=1;
	    index[0]=i_in;
	    index[1]=j_in-1;
	    index[2]=l_in;
	  }
	}
      }
    }
  }
  return found;
  
}; /* get_closest_index */
	      


/*--------------------------------------------------------------------------

calculate normalized great circle distance between v1 and v2 
double normalize_great_circle_distance(v1, v2)
---------------------------------------------------------------------------*/
double normalize_great_circle_distance(const double *v1, const double *v2)
{
  double dist;

  dist=(v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
    /sqrt((v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2])                 
	  *(v2[0]*v2[0]+v2[1]*v2[1]+v2[2]*v2[2]));
  dist = sign(min(1.,fabs(dist)),dist);
  dist = acos(dist);
  return dist;
  
}; /* normalize_great_circle_distance */

/*------------------------------------------------------------------
  double spherical_angle(v1, v2, v3)

  calculate spherical angle of a triangle formed by v1, v2 and v3 at v1                                                            
  ------------------------------------------------------------------*/
/* double spherical_angle(double *v1, double *v2, double *v3) */
/* { */
/*   double angle; */
/*   double px, py, pz, qx, qy, qz, abs_p, abs_q; */

/*   /* vector product between v1 and v2 */
/*   px = v1[1]*v2[2] - v1[2]*v2[1]; */
/*   py = v1[2]*v2[0] - v1[0]*v2[2]; */
/*   pz = v1[0]*v2[1] - v1[1]*v2[0]; */
/*   /* vector product between v1 and v3  */
/*   qx = v1[1]*v3[2] - v1[2]*v3[1]; */
/*   qy = v1[2]*v3[0] - v1[0]*v3[2]; */
/*   qz = v1[0]*v3[1] - v1[1]*v3[0]; */
    
/*   /* angle between p and q */ 
/*   abs_p=px*px+py*py+pz*pz; */
/*   abs_q=qx*qx+qy*qy+qz*qz; */
/*   if (abs_p*abs_q==0.) */
/*     angle=0.; */
/*   else { */
/*     angle = (px*qx+py*qy+pz*qz)/sqrt(abs_p*abs_q); */
/*     angle = sign(min(1.,fabs(angle)),angle); */
/*     angle = acos(angle); */
/*   } */

/*   return angle; */
/* }; /* spherical_angle */

/*---------------------------------------------------------------------
  double dist2side(v1, v2, point)
  calculate shortest normalized distance on sphere                 
  from point to straight line defined by v1 and v2                 
  ------------------------------------------------------------------*/
double dist2side(const double *v1, const double *v2, const double *point)
{
  double angle, side;

  angle = spherical_angle(v1, v2, point);
  side  = normalize_great_circle_distance(v1, point);

  return (asin(sin(side)*sin(angle)));

};/* dist2side */

  
int max_weight_index( double *var, int nvar)
{

  int ind, i;

  ind = 0;

  for(i=1; i<nvar; i++) {
    if(var[i]>var[ind]) ind = i;
  }

  return ind;
}

/*------------------------------------------------------------------------------
  void do_latlon_coarsening(var_latlon, ylat, nlon, nlat, nz,     
                            var_latlon_crs, nlon_crs, nlat_crs,   
                            finer_steps, misval, varmisval)

  calculate variable on coarser latlon grid                        
  by doubling spatial resolution and preserving volume means       
  ---------------------------------------------------------------------------*/
void do_latlon_coarsening(const double *var_latlon, const double *ylat, int nlon, int nlat, int nz,    
			  double *var_latlon_crs, int finer_steps, int has_missing, double missvalue)
{

  double  *var_latlon_old, *ylat_old, *var_latlon_new;
  double  dlat;
  int     nlon_old, nlat_old, nlon_new, nlat_new, steps, i, j;
  int     nlon_crs, nlat_crs;
  
  nlon_crs=nlon/pow(2,finer_steps);
  nlat_crs=(nlat-1)/pow(2,finer_steps)+1;
  switch (finer_steps) {
  case 0:
    if (nlon_crs !=nlon || nlat_crs != nlat) mpp_error("bilinear_interp(do_latlon_coarsening): grid dimensions don't match");
    for(i=0; i<nlon*nlat*nz; i++) var_latlon_crs[i] = var_latlon[i];
    break;
  case 1:
    redu2x(var_latlon, ylat, nlon, nlat, var_latlon_crs, nlon_crs, nlat_crs, nz, has_missing, missvalue);
    break;
  default:
    nlon_new=nlon;
    nlat_new=nlat;
    for(steps=1; steps<=finer_steps; steps++) {
      nlon_old=nlon_new;
      nlat_old=nlat_new;
      var_latlon_old = (double *)malloc(nlon_old*nlat_old*nz*sizeof(double));
      ylat_old       = (double *)malloc(nlat_old*sizeof(double));
      if (steps==1) {
	for(i=0; i<nlat; i++) ylat_old[i] = ylat[i]; 
	for(i=0; i<nlon_old*nlat_old*nz; i++) var_latlon_old[i] = var_latlon[i];
      }
      else {
	dlat=M_PI/(nlat_old-1);
	ylat_old[0]=-0.5*M_PI;
	ylat_old[nlat_old-1]= 0.5*M_PI;
	for(j=1; j<nlat_old-1; j++) ylat_old[j] = ylat_old[0] + j*dlat;
	for(i=0; i<nlon_old*nlat_old*nz; i++) var_latlon_old[i] = var_latlon_new[i];
	free(var_latlon_new);
      }
          
      nlon_new=nlon_new/2;
      nlat_new=(nlat_new-1)/2+1;
      var_latlon_new = (double *)malloc(nlon_new*nlat_new*nz*sizeof(double));
      redu2x(var_latlon_old, ylat_old, nlon_old, nlat_old, var_latlon_new, nlon_new, nlat_new, nz, has_missing, missvalue);
      free(var_latlon_old);
      free(ylat_old);
    }
    for(i=0; i<nlon_new*nlat_new*nz; i++) var_latlon_crs[i] = var_latlon_new[i];
    free(var_latlon_new);
  }
}; /* do_latlon_coarsening */

/*------------------------------------------------------------------------------
  void redu2x(varfin, yfin, nxfin, nyfin, varcrs, nxcrs, nycrs)
  this routine is for reducing fvccm data by a factor of 2       
  volume averaging for all data except at the poles              
  original developer: S.-J. Lin                                  
  ----------------------------------------------------------------------------*/
void redu2x(const double *varfin, const double *yfin, int nxfin, int nyfin, double *varcrs,
	    int nxcrs, int nycrs, int nz, int has_missing, double missvalue)
{
  double  *cosp, *acosp, *vartmp;
  int     i, j, k, i2, j2, n1, n2;

  /*------------------------------------------------------------------
    calculate cosine of latitude                                   
    trick in cosp needed to maintain a constant field              
    ----------------------------------------------------------------*/
  cosp  = (double *)malloc(nyfin*sizeof(double));
  acosp = (double *)malloc(nyfin*sizeof(double));
  vartmp = (double *)malloc(nxcrs*nyfin*nz*sizeof(double));
  cosp[0]=0.;
  cosp[nyfin-1]=0.;
  for(j=1; j<nyfin-1; j++) cosp[j]  = cos(yfin[j]);
  for(j=1; j<nyfin-1; j++) acosp[j] = 1./(cosp[j]+0.5*(cosp[j-1]+cosp[j+1]));

  /*----------------------------------------------------------------
    x-sweep                                                        
    ----------------------------------------------------------------*/
  if(has_missing) {
    for(k=0; k<nz; k++) for(j=1; j<nyfin-1; j++) {
      n1 = k*nxfin*nyfin+j*nxfin;
      n2 = k*nxcrs*nyfin+j*nxcrs;
      if (varfin[n1+nxfin-1] == missvalue || varfin[n1] == missvalue || varfin[n1+1] == missvalue)
	vartmp[n2] = missvalue;
      else
	vartmp[n2] = 0.25*(varfin[n1+nxfin-1]+2.*varfin[n1]+varfin[n1+1]);
      for(i=2; i<nxfin-1; i+=2) {
	i2 = i/2;
	if (varfin[n1+i-1] == missvalue || varfin[n1+i] == missvalue || varfin[n1+i+1] == missvalue) 
	  vartmp[n2+i2] = missvalue;
        else
	  vartmp[n2+i2] = 0.25*(varfin[n1+i-1]+2.*varfin[n1+i]+varfin[n1+i+1]);
      }
    }
  }
  else {
    for(k=0; k<nz; k++) for(j=1; j<nyfin-1; j++) {
      n1 = k*nxfin*nyfin+j*nxfin;
      n2 = k*nxcrs*nyfin+j*nxcrs;
      vartmp[n2] = 0.25*(varfin[n1+nxfin-1]+2.*varfin[n1]+varfin[n1+1]);
      for(i=2; i<nxfin-1; i+=2) {
	i2 = i/2;
	vartmp[n2+i2] = 0.25*(varfin[n1+i-1]+2.*varfin[n1+i]+varfin[n1+i+1]);
      }
    }
  }
  /*---------------------------------------------------------------------
    poles:                                                         
    this code segment works for both the scalar and vector fields. 
    Winds at poles are wave-1; the follwoing is quick & dirty yet the correct way
    The skipping method. A more rigorous way is to                 
    recompute the wave-1 components for the coarser grid.          
    --------------------------------------------------------------------*/
  for(k=0; k<nz; k++) for(i=0; i<nxcrs; i++) {
    i2 = i*2;
    n1 = k*nxcrs*nycrs;
    n2 = k*nxfin*nyfin; 
    varcrs[n1+i] = varfin[n2+i2];
    varcrs[n1+(nycrs-1)*nxcrs+i] = varfin[n1+(nyfin-1)*nxfin+i2];
  }
  /*----------------------------------------------------------------
    y-sweep                                                        
    ----------------------------------------------------------------*/
  if (has_missing) {
    for(k=0; k<nz; k++) for(j=1; j<nyfin-1; j++) for(i=0; i<nxcrs; i++) {
      n1 = k*nxcrs*nyfin+j*nxcrs+i;
      if (vartmp[n1] /= missvalue) vartmp[n1] *= cosp[j];
    }
    for(k=0; k<nz; k++) for(j=2; j<nyfin-2; j+=2) {
      j2 = j/2;
      for(i=0; i<nxcrs; i++) {
	n1 = k*nxcrs*nyfin + i;
	n2 = k*nxcrs*nycrs + i;
	if (vartmp[n1+j*nxcrs] == missvalue || vartmp[n1+(j-1)*nxcrs] == missvalue
	    || vartmp[n1+(j+1)*nxcrs] == missvalue )
	  varcrs[n2+j2*nxcrs] = missvalue;
	else
	  varcrs[n2+j2*nxcrs] = acosp[j]*(vartmp[n1+j*nxcrs] + 0.5*(vartmp[n1+(j-1)*nxcrs]+ vartmp[n1+(j+1)*nxcrs]));
      }
    }
  }
  else {
    for(k=0; k<nz; k++) for(j=1; j<nyfin-1; j++) for(i=0; i<nxcrs; i++) {
      n1 = k*nxcrs*nyfin+j*nxcrs+i;
      vartmp[n1] *= cosp[j];
    }
    for(k=0; k<nz; k++) for(j=2; j<nyfin-2; j+=2) {
      j2 = j/2;
      for(i=0; i<nxcrs; i++) {
	n1 = k*nxcrs*nyfin + i;
	n2 = k*nxcrs*nycrs + i;
	varcrs[n2+j2*nxcrs] = acosp[j]*(vartmp[n1+j*nxcrs] + 0.5*(vartmp[n1+(j-1)*nxcrs]+ vartmp[n1+(j+1)*nxcrs]));
      }
    }
  }

  free(cosp);
  free(acosp);
  free(vartmp);
  
}; /*redu2x*/