/*  destag.c
   
    This plug-in takes a staggered grid in raw wrfout data
    (e.g., u) and interpolates it to the mass grid, 
    effectively "destaggering" it. Properly destaggered
    arrays will have sizes equal to their mass-grid
    counterparts for given ranges of x,y,z.
    In principle, this routine could work
    on staggered grids in general.

   *  Arguments expected:
   *    1) Field Q to destagger [expression, grid]
   *    2) Dimension along which to destagger [char]
   * 
   *  Returns:
   *    1) Field Q interpolated to the destaggered grid

    Example use in a GrADS script:
        'reinit'

        'xdfopen raw_m.ctl'   ;* read wrfout data on the mass grid
        'xdfopen raw_u.ctl'   ;* read wrfout data on the U grid
        'xdfopen raw_v.ctl'   ;* read wrfout data on the V grid

        ;* destagger the u-component of the wind
        'define ucen = destag(u.2,x)'
        
        ;* destagger the v-component of the wind
        'define vcen = destag(v.3,y)'
        
        ;* close raw_v and raw_u files
        'close 3'; 'close 2'
        [...]

************************************************************************/

#include "grads.h"    

gaint destag (struct gafunc *pfc, struct gastat *pst, struct gaudpinfo *pinfo) {
  gaint i,j,rc;
  struct gagrid *pgr;
  gadouble *val;
  char *valu;
  char stagdim;

  /* Check if this plug-in is still compatible with the GrADS interface */
  if (pinfo->version != UDPVERS) printf("Error: version mismatch. Recompile this plug-in\n");

  /* Check for the proper number of arguments */
  if (pfc->argnum!=2) {
    gaprnt (0,"Error from DESTAG: Two arguments expected. \n");
    return (1);
  }

  /* Invoke gaexpr to evaluate the first argument, 
   * which should be a GrADS expression.
     The result gets stored in the gastat structure 'pst' */
  rc = gaexpr(pfc->argpnt[0],pst);
  if (rc) return (rc);

  /* If the expression is a grid ... */
  if (pst->type==1) {
    /* Set up pointers to the gridded data structure, the grid of data, and the undef mask  */
    pgr = pst->result.pgr;
    val = pgr->grid;
    valu = pgr->umask;
   
    /* two arguments are expected: the name of the staggered field (e.g., u)
     * and the staggered dimension (x,y,z). */
    stagdim = pfc->argpnt[1][0];
    switch(stagdim) {
      case 'x': 
        gaprnt (0,"Note from DESTAG: destaggering along x. \n");
        for (j=0; j < pgr->jsiz ; j++) {
          for (i=0; i < pgr->isiz ; i++) {
            if (i == pgr->isiz-1) {
              /* we're at the end of the staggered grid's row.
               * the unstaggered array is undefined here. */
              *valu=0;
            } else {
              if (*valu!=0 && *(valu+1)!=0) *val = (*val + *(val+1)) / 2.0;
            }
            val++; valu++;
          } 
        }
        break;
      case 'y': 
        gaprnt (0,"Note from DESTAG: destaggering along y. \n");
        for (j=0; j < pgr->jsiz ; j++) {
          for (i=0; i < pgr->isiz ; i++) {
            if (j == pgr->jsiz-1) {
              /* we're at the end of the staggered grid's column.
               * the unstaggered array is undefined here. */
              *valu=0;
            } else {
              if (*valu!=0 && *(valu+pgr->isiz)!=0) *val = (*val + *(val+pgr->isiz)) / 2.0;
            }
            val++; valu++;
          } 
        }
        break;
      case 'z': 
        gaprnt (0,"Note from DESTAG: destaggering along z NOT IMPLEMENTED. \n");
        break;
      default: 
        gaprnt (0,"Error from DESTAG: unrecognized stagger dimension. \n");
        gaprnt (0,"Options include x, y, and z. \n");
        return (1);
        break;
    }

  } 
  /* ...the expression is station data */
  else {              
    gaprnt (0,"Error from DESTAG: bad data type (non-grid). \n");
    return (1);
  }

  /* finish without errors */
  return (0);
}