/****************************************************************
                     87-30-87

FUNCTION:  INTGRID   (derived from REGRID)

This function is needed for calculating the transport in isothermal
layers, for example.

Integrate an array f(zo) over each depth segment specified by
zn(i) to zn(i+1).
This routine approximates the function
as a set of line segments connecting the input data points.


MODIFIED: Tue  09-27-1988
   Removed a bug causing extrapolation when the function is
   called with no overlap between the new and old grids.
   Changed the location of initialization of imin, imax,
   for readability.


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

#include "common.h"
#include "vector.h"

/************      intgrid     *****************************/

#if PROTOTYPE_ALLOWED
void intgrid(float zo[], float fo[], float zn[], float fn[], int m, int n,
  int *n0_ptr, int *nout_ptr)
#else
void intgrid(zo, fo, zn, fn, m, n, n0_ptr, nout_ptr)
float     zo[],          /* input array of m depths: "z_old"  */
          fo[],          /* input array of m; f(z_old)        */
          zn[],          /* input array of n depths: "z_new". */
          fn[];          /* output array of n-1; f(z_new)     */
int       m,             /* number of points in input arrays  */
          n;             /*    "   "    "     " output  "     */
                         /* additional output:                */
                         /* calling routine supplies the space*/
int       *n0_ptr,       /* index of first good point in fn   */
          *nout_ptr;     /* number of good points in fn       */
                         /*  (should be n-1)                  */
#endif
{
     int  imin,          /* first new gridpoint index for fn  */
          imax;          /* last new gridpoint index for fn   */
     int  i, j;          /* index counters; i:new; j:old      */
     int  jlow, jhigh;
     int  k, bad_flag;
     double    zlow,     /* zn(i)                             */
               zhigh,    /* zn(i+1)                           */
               flow,     /* f(zlow)  (interpolated)           */
               fhigh;    /* f(zhigh) (interpolated)           */
     double    integral;

/****************************************************************
                 zo(jlow)              zo(jhigh)
         ----------|--------------------|-------------       zo
        -------------|-----------------|---------------   zn
                    zn(i)            zn(i+1)
******************************************************************/

/****Find limits of gridpoint index for new grid. *****************

           |-----------------------------------------|        zn (case a)
        -----------------------------------------------       zo
     ---|---------------------------------------------|------ zn (case b)
       imin                                         imax

     There can be no extrapolation in this version of regrid;
     data on the old grid must match or overhang the part of the
     new grid for which interpolation will be done.
*********************************************************************/

     imin = 0;
     imax = n-1;

     while( (zn[imin] < zo[0]) && (imin < n-1) ) imin++;
     while( (zn[imax] > zo[m-1]) && (imax > 0) )  imax--;

     /*** If there is no overlap between the old and
          new grids, there can be no interpolation.
          Set number of points regridded to zero, and
          quit.  If this occurs, then either imin = n-1
          or imax = 0; the converse is NOT true.
          Note also that since we are integrating from one
          new grid point to the next, we need at least 2 valid
          points.  Hence we exit if imax = imin.
     ***/
     if ( (zn[imin] < zo[0]) || (zn[imax] > zo[m-1]) || (imax == imin) )
     {
          *nout_ptr = 0;
          return;
     }

     /*** Otherwise, the first interpolated point will be
          at index imin (to imin+1), and the last will be
          at imax-1 (to imax).
     ***/
     *n0_ptr = imin;
     *nout_ptr = imax - imin;


     /*** Initialize ?high to start the loop.
          In the loop ?high is immediately replaced by ?low,
          and a new ?high is calculated.
     ***/
     zhigh = zn[imin];
     jhigh = 1;
     while( (zo[jhigh] < zhigh) && (jhigh < m-1) ) jhigh++;
     /* This is one past initial jlow.  */
     fhigh = interpolate( zo+jhigh-1, fo+jhigh-1, zhigh);

     for ( i=imin; i<imax; i++)
     {          /* I think this could be rearranged for compactness */
          zlow = zhigh;
          zhigh = zn[i+1];
          if (zhigh<=zlow)
          {
               printf("\nError in INTGRID: i= %d  zhigh= %f  zlow= %f\n",
                                            i, zhigh, zlow);
               exit(-1);
          }

          jlow = jhigh-1;      /* new jlow is 1 to left of old jhigh */
          while ( (zo[jhigh] < zhigh) && (jhigh < m-1) ) jhigh++;

          for (bad_flag=0, k=jlow; k<=jhigh; k++)
               if (fo[k]>=ADJ_BADFLOAT) bad_flag = 1;

          if (bad_flag)  
          {
               fn[i] = BADFLOAT;
               printf("\nBADFLOAT found in input to INTGRID, i= %d\n", i);
          }
          else
          {
               flow = fhigh;
               fhigh = interpolate(zo+jhigh-1, fo+jhigh-1, zhigh);
               if ( jhigh == (jlow+1) )
                    fn[i] = average(flow, fhigh) * (zhigh-zlow);
               else                               /* all other cases       */
               {
                    /** end pieces are needed first **/
                    integral = average(flow,fo[jlow+1]) * (zo[jlow+1]-zlow)
                             + average(fhigh,fo[jhigh-1]) * (zhigh-zo[jhigh-1]);
                    /** now the middle sections  **/
                    for (j=jlow+1; j<jhigh-1; j++)
                         integral += average(fo[j], fo[j+1]) * (zo[j+1]-zo[j]);
                    fn[i] = integral;
               }
          }
     }
     return;
}
/**************** end of intgrid()  ************************************/