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

   VGRID.C module for ADCPSECT

   This module maps the ADCP velocity field onto a specified
   vertical grid by one of several means.

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

#include "geninc.h"
#include "ioserv.h"
#include "vector.h"
#include "use_db.h"
#include "vstat.h"
#include "pos_to_m.h"
#include "adcp.h"
#include "adcpsect.h"

#define MAX_SURFACE_EXTENSION 3  /* This could be made into a parameter;
                                    it is the maximum number of bins
                                    at the top of the profile which will
                                    be filled if the integrate_top option
                                    is used.
                                  */
extern int verbose;

/*--> regridding variables and constants: */
   float sigma_sva_grid[NDMAX_NG];
   float adj_depth[NDMAX_NG];
   float adj_u[NDMAX_NG];
   float adj_v[NDMAX_NG];
   int adj_ndepth;  /* new grid; = rg_number if regridding is done */
   int adj_nu;      /* adj_ndepth or one less; npts in array of centers */
   VERTICAL_GRID_TYPE vgrid;
   int regrid_var = RV_NONE;
   int regrid_method;

static int rg_number;
static double rg_origin, rg_increment;
static int top = 0;


static VERTICAL_GRID_TYPE vgrid_array[] =
{
   /* array     cnt_array   ca_alloc incr. down  label precision   #defined index */
   {adj_depth     ,  NULL , 0, TRUE   , "depth m"    , 0 }, /* RV_NONE        */
   {adj_depth     ,  NULL , 0, TRUE   , "depth m"    , 0 }, /* RV_DEPTH       */
   {sigma_sva_grid,  NULL , 0, FALSE  , "svan cl/t"  , 0 }, /* RV_SVAN        */
   {sigma_sva_grid,  NULL , 0, FALSE  , "p svan cl/t", 0 }, /* RV_POT_SVAN    */
   {sigma_sva_grid,  NULL , 0, TRUE   , "- sigma-th" , 2 }, /* RV_SIGMA_THETA */
   {sigma_sva_grid,  NULL , 0, TRUE   , "- sigma-t"  , 2 }  /* RV_SIGMA_T     */
};


NAME_LIST_ENTRY_TYPE regrid_methods[] =
{
   {"interpolate"       , RM_INTERPOLATE},
   {"centered_average"  , RM_CENTERED_AV},
   {"integrate"         , RM_INTEGRATE},
   {"integrate_top"     , RM_INTEGRATE_TOP},
   {"average"           , RM_AVERAGE},
   {"per_grid_interval" , RM_TRANSDERIV},
   {NULL                , 0}
};
NAME_LIST_ENTRY_TYPE regrid_variables[] =
{
   {"none"              , RV_NONE},
   {"svan"              , RV_SVAN},
   {"pot_svan"          , RV_POT_SVAN},
   {"depth"             , RV_DEPTH},
   {"sigma_theta"       , RV_SIGMA_THETA},
   {"sigma_t"           , RV_SIGMA_T},
   {NULL                , 0}
};
static NAME_LIST_TYPE regrid_name_lists[] =
{
   regrid_methods,
   regrid_variables
};

static PARAMETER_LIST_ENTRY_TYPE regrid_params[] =
{
   {""   , " %s"     , &regrid_method  , TYPE_TRANS},
   {""   , " %s"     , &regrid_var     , TYPE_TRANS + 256},
   {NULL , NULL      , NULL            , 0}
};

#if PROTOTYPE_ALLOWED
int init_regrid(FILE *fp_cnt)
#else
int init_regrid(fp_cnt)
FILE *fp_cnt;
#endif
{
   char sub_option[20];
   double first, last, other;
   float  *grid_array;
   int i;


   check_error(get_parameters(fp_cnt, regrid_params, regrid_name_lists),
                  "regrid parameters");

   check_error( (fscanf(fp_cnt, " %10s",
             sub_option) != 1),
            "reading regrid options");


   grid_array = ((regrid_var == RV_NONE || regrid_var == RV_DEPTH) ?
                  adj_depth : sigma_sva_grid);

   if (regrid_method == RM_INTEGRATE_TOP)
   {
      grid_array[0] = 0.0;
      top = 1;
   }
   if (strcmp(sub_option, "grid") == 0)
   {
      check_error( (fscanf(fp_cnt,
               " number= %d origin= %lf increment= %lf",
               &rg_number, &rg_origin, &rg_increment) != 3),
               "reading grid parameters");
      rg_number = min_val(rg_number, NDMAX_NG-1);
      /* We want the grid to go from large svan (shallow depths)
         to small svan (greater depths), so that it will be in
         the same sense as the svan (or depth) array from the
         ctd profile.
      */
      first = (rg_origin + rg_number * rg_increment);
      last  = rg_origin;
      if (((regrid_var == RV_SVAN)        && (first < last)) ||
          ((regrid_var == RV_POT_SVAN)    && (first < last)) ||
          ((regrid_var == RV_DEPTH)       && (first > last)) ||
          ((regrid_var == RV_SIGMA_THETA) && (first > last)) )
                                /* Wrong order; swap them. */
      {
         other = first;
         first = last;
         last = other;
      }
      make_uniform_grid(first, last, rg_number, grid_array+top);
   }

   else if (strcmp(sub_option, "grid_list") == 0)
   {
      check_error( (fscanf(fp_cnt, " number= %d boundaries:",
         &rg_number) != 1),
               "reading grid parameters");
      check_error( rg_number > NDMAX_NG-top, "too many grid boundaries");

      /* Read in the grid boundaries. */
      for (i=top; i<rg_number+top; i++)
         check_error( (fscanf(fp_cnt, " %f", grid_array+i) != 1),
                     "reading grid array");

      /* Make sure that depth increases, or svan decreases. */
      if (regrid_var == RV_SVAN || regrid_var == RV_POT_SVAN)
         for (i=1+top; i<rg_number+top; i++)
            check_error( grid_array[i] >= grid_array[i-1],
               "boundaries equal or out of order");
      if (regrid_var == RV_DEPTH || regrid_var == RV_SIGMA_THETA)
         for (i=1+top; i<rg_number+top; i++)
            check_error( grid_array[i] <= grid_array[i-1],
               "boundaries equal or out of order");
   }
   else check_error( 1, "regrid: ... grid or grid_list");
   if (regrid_var == RV_DEPTH)
      for (i=top; i<rg_number+top; i++)
         check_error((grid_array[i] < 0.0),
         "Negative depth in grid; depths must be positive.");


   /* If necessary, shift the depth grid to eliminate a
      duplicate value of zero. This could occur if the depth
      grid started at zero and the integrate_top option were
      selected.
   */
   if (top && (regrid_var == RV_DEPTH) &&
               (grid_array[0] == grid_array[1]))
   {
      for (i=1+top; i<rg_number+top; i++)
      {
         grid_array[i-1] = grid_array[i];
      }
      rg_number--;
   }
   return(0);
}                       /* init_regrid() */

#if PROTOTYPE_ALLOWED
void write_vgrid_params(FILE *fp)
#else
void write_vgrid_params(fp)
FILE *fp;
#endif
{
    fprintf(fp, "\n");
    check_error(print_parameters(fp, regrid_params, regrid_name_lists, "   "),
                  "regrid parameters");
}


/*
   Since init_vgrid is not called unless there is to be regridding,
   the vertical grid setup must be done by the main program after
   reading in all the parameters.
*/
void set_vgrid(ndepth)
int ndepth;
{
   int i;
   extern VERTICAL_GRID_TYPE vgrid_array[];
   extern int adj_ndepth, adj_nu, regrid_method, regrid_var;

   adj_ndepth = regrid_var ? rg_number+top : ndepth+1;  /*ef*/

   vgrid = vgrid_array[regrid_var];  /* structure assignment */
   if (regrid_method == RM_INTERPOLATE ||
       regrid_method == RM_CENTERED_AV ||
       regrid_var    == RV_NONE )
       /* grid point is middle of increment */
   {
      vgrid.center_array = vgrid.array;
      vgrid.ca_allocated = 0; /* no calloc, so don't try to free it! */
      adj_nu = adj_ndepth;
   }
   else    /* any other regridding averages over the increment */
   {
      adj_nu = adj_ndepth - 1;
      vgrid.center_array = (float *)calloc(adj_nu, sizeof(float));
      vgrid.ca_allocated = 1;  /* flag so we know to free it */
      check_error(vgrid.center_array == NULL,
                  "init_regrid; memory allocation for center_array");
      for (i=0; i<adj_nu; i++)
      {
         vgrid.center_array[i] = (0.5) *
                        (vgrid.array[i] + vgrid.array[i+1]);
      }
      if (regrid_var != RV_NONE &&
          regrid_var != RV_DEPTH && top)
         vgrid.center_array[0] = 0.0;  /* for integrate_top option */
   }
}                       /* set_vgrid() */


void regrid_velocity(adcp, ctd)
ADCP_PROFILE_TYPE *adcp;
CTD_PROFILE_TYPE *ctd;
{
   int i, igood;
   int max_extend = MAX_SURFACE_EXTENSION;
   double length;
   static YMDHMS_TIME_TYPE last_ctd_time;
   float *grid_ptr;
   extern int top, adj_nu, adj_ndepth;

   if (top == 1)
   {
      /* The checking for igood starts at zero, just in case
         some other part of the system has used another method
         of assigning the first ([0]) velocity.
      */
      for (igood = 0; (igood < adcp->nu)
                   && (adcp->u[igood] >= ADJ_BADFLOAT); igood++) ;
      if (igood <= max_extend && igood < adcp->nu)
      {
         for (i=0; i<igood; i++)
         {
            adcp->u[i] = adcp->u[igood];
            adcp->v[i] = adcp->v[igood];
         }
      }
   }
   if (regrid_var != RV_NONE)
   {
      /* The CTD profile must be used to calculate adj_depth
         unless depth is the grid variable.
      */
      if (regrid_var != RV_DEPTH)
      {
         /* Calculate a new adj_depth grid from the CTD
            profile only if it has changed since the
            previous call to this routine.
         */
         if (TIMDIF(&last_ctd_time, &ctd->time))
         {
            if (verbose)
            {
               printf("\n ctd->sigma_sva: \n");
               print_array_float(stdout, ctd->sigma_sva, ctd->n, 8, " %7.3g ");
               printf("\n ctd->depth: \n");
               print_array_float(stdout, ctd->depth, ctd->n, 8, " %7.3g ");
               printf("\n sigma_sva_grid: \n");
               print_array_float(stdout, sigma_sva_grid, adj_ndepth, 8, " %7.3g ");
            }
            s_regridli(ctd->sigma_sva, ctd->depth,
                           sigma_sva_grid+top, adj_depth+top,
                              ctd->n-top, adj_ndepth-top);
            last_ctd_time = ctd->time;
         }
         /* If top = 1 (integrate_top option), ensure that
            adj_depth consists of a single sequence of good depth
            values, starting from zero.  Previous values should be
            flagged with BADFLOAT.  It is assumed that if adj_depth
            needs to be changed, it will have a first value of zero
            followed by a single sequence of BADFLOAT and then the
            good depth values.  Therefore it is necessary only to
            find the transition from bad to good, set the last
            bad to zero (if the first good is not already zero),
            and change the first value from zero to bad.
         */
         if (top)
         {
            for (i=1; i<adj_ndepth-1; i++)
            {
               if (adj_depth[i] >= ADJ_BADFLOAT &&
                   adj_depth[i+1] <= ADJ_BADFLOAT)
               {
                  if (adj_depth[i+1] > 0.0)
                     adj_depth[i] = 0.0;
                  adj_depth[0] = BADFLOAT;    /* (It was 0.) */
                  break;
               }
            } /* end loop; exit here with break */
         }

      }  /* adj_depth has now been calculated. */

      if (verbose)
      {
         printf("\n depth: \n");
         print_array_float(stdout, adcp->depth, adcp->nu, 8, " %7.3g ");
         printf("\n u: \n");
         print_array_float(stdout, adcp->u, adcp->nu, 8, " %7.3g ");
         printf("\n adj_depth: \n");
         print_array_float(stdout, adj_depth, adj_ndepth, 8, " %7.3g ");
      }

      /* The following regrid methods require integration
         from one grid point to the next.
      */
      if (     regrid_method == RM_AVERAGE
            || regrid_method == RM_TRANSDERIV
            || regrid_method == RM_INTEGRATE
            || regrid_method == RM_INTEGRATE_TOP )
      {
         s_intgrid(adcp->depth, adcp->u, adj_depth, adj_u, adcp->nu, adj_ndepth);
         s_intgrid(adcp->depth, adcp->v, adj_depth, adj_v, adcp->nu, adj_ndepth);
         if ((regrid_method == RM_AVERAGE) ||
             (regrid_method == RM_TRANSDERIV))
         {
            /* Divide the depth integral by a suitable length.
               For averaging, this length is always the depth
               interval over which the average was taken.
               For RM_TRANSDERIV, the length is the grid variable
               interval.  We want the result to be the derivative of
               the transport per unit width with respect to the
               grid variable.  If the variable is depth, then this
               is the same as the average.  If it is sigma_sva, then
               integrating the "transport_derivative" wrt sigma_sva
               gives the transport per unit width.
            */
            if ((regrid_method == RM_TRANSDERIV) &&
                (regrid_var != RV_DEPTH))
               grid_ptr = sigma_sva_grid;
            else
               grid_ptr = adj_depth;

            for (i=0; i<adj_nu; i++)
            {
               if (good_float(adj_u[i]))
               {
                  length = grid_ptr[i+1] - grid_ptr[i];
                  adj_u[i] /= length;
                  adj_v[i] /= length;
               }
            }
         } /* end of RM_AVERAGE or RM_TRANSDERIV block */
      }
      /* A "centered average" is the average value in a
         region surrounding the grid point.
      */
      else  if (regrid_method == RM_CENTERED_AV)
      {
         s_regrid(adcp->depth, adcp->u, adj_depth, adj_u, adcp->nu, adj_ndepth);
         s_regrid(adcp->depth, adcp->v, adj_depth, adj_v, adcp->nu, adj_ndepth);
      }
      /* Interpolation is simple linear interpolation from
         the old grid to the new one.
      */
      else    /* RM_INTERPOLATE */
      {
         s_regridli(adcp->depth, adcp->u, adj_depth, adj_u,
                     adcp->nu, adj_ndepth);
         s_regridli(adcp->depth, adcp->v, adj_depth, adj_v,
                     adcp->nu, adj_ndepth);
      }
      /* End of the regridding operation. */

      /* we are relying on the data being flagged bad, instead
         of trying to keep track of the range of good data.
      */

      if (verbose)
      {
         printf("\n adj_u: \n");
         print_array_float(stdout, adj_u, adj_nu, 4, " %7.3g ");
      }

   }
   else    /* no regridding */
   {
      /* If there is no regridding, just copy without adjusting. */
      memcpy(adj_depth, adcp->depth, adj_ndepth*sizeof(float));
      memcpy(adj_u, adcp->u, adcp->nu*sizeof(float));
      memcpy(adj_v, adcp->v, adcp->nu*sizeof(float));
      /* s_regrid() and related functions fill gaps with BADFLOAT;
         this must be done manually here where we are not using
         one of those functions.  Note that adcp->nu has been
         changed by the editing operations and the contents
         of the access variable structure, while adj_nu is set
         once at the start of the program and left alone.
         Therefore we do not want to just set adj_nu equal to
         adcp->nu.
      */
      for (i=adcp->nu; i<adj_nu; i++)
      {
         adj_u[i] = BADFLOAT;
         adj_v[i] = BADFLOAT;
      }
   }
}                       /* regrid_velocity() */


void close_vgrid()
{
   if (vgrid.ca_allocated == 1) free(vgrid.center_array);
}