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

   CTD module for ADCPSECT

*/
#include "geninc.h"
#include "ioserv.h"
#include "vector.h"
#include "use_db.h"
#include "vstat.h"
#include "ocean.h"
#include "pos_to_m.h"
#include "adcp.h"
#include "adcpsect.h"       /* includes ctd.h */
#include "data_id.h"

#define TRUE  1
#define FALSE 0
#define sqr(x)  (x)*(x)
#define good_float(x)  ( (x) < ADJ_BADFLOAT )


         int   use_ctd = FALSE;            /* global flag; ctd data in use */
static   FILE_NAME_TYPE  ctd_dbname;          /* ctd option */
         double nominal_svel = NOMINAL_SVEL;

#if PROTOTYPE_ALLOWED
int init_ctd(FILE *fp_cnt)
#else
int init_ctd(fp_cnt)
FILE *fp_cnt;
#endif
{
   use_ctd = TRUE;
   check_error( (fscanf(fp_cnt," dbname: %79s", ctd_dbname) != 1),
                                 "Can't read ctd_dbname");
   check_dbopen(CTD_DB, ctd_dbname, READ_ONLY, DIR_IN_MEMORY);
   printf("\n\n DATABASE:  %s SUCCESSFULLY OPENED",ctd_dbname);

   check_dbset(ADCP_DB, "DBSET, after opening ctd db");
   return(0);
}

#if PROTOTYPE_ALLOWED
void write_ctd_params(FILE *fp)
#else
void write_ctd_params(fp)
FILE *fp;
#endif
{
   fprintf(fp, "\n CTD Database: %s", ctd_dbname);
   fprintf(fp, "\n    depth correction with nominal sound speed = %7.2f",
                                                         nominal_svel);
}

#if PROTOTYPE_ALLOWED
void set_ctd_grid(CTD_PROFILE_TYPE *ctd, int regrid_var)
#else
void set_ctd_grid(ctd, regrid_var)
CTD_PROFILE_TYPE *ctd;
int regrid_var;
#endif
{
   if (regrid_var == RV_SVAN)
      ctd->density_var = CTD_SVAN;
   else if (regrid_var == RV_POT_SVAN)
      ctd->density_var = CTD_POT_SVAN;
   else if (regrid_var == RV_SIGMA_THETA)
      ctd->density_var = CTD_SIGMA_THETA;
   else if (regrid_var == RV_SIGMA_T)
      ctd->density_var = CTD_SIGMA_T;
   else
      ctd->density_var = CTD_NONE;
}

#if PROTOTYPE_ALLOWED
void get_ctd_profile(CTD_PROFILE_TYPE *prof)
#else
void get_ctd_profile(prof)
CTD_PROFILE_TYPE *prof;
#endif
{
   double sigma;           /* dummy used in SVAN */
   double theta;
   int step = 1;
   int ngaps = 0;     /* max number of interpolations done */
   int ierr;
   int n, ng;
   int i;

   get_latlon(&(prof->position));
   check_dbget(TIME, (char *) &(prof->time), sizeof(prof->time),
                                       "DBGET, time");

   prof->n = check_dbget_f(P, prof->pressure, NPMAX, "DBGET_F, pressure");

   n = check_dbget_f(TEMPERATURE, prof->temperature, NPMAX,
                                       "DBGET_F, temperature");
   prof->n = min_val(prof->n, n);

   n = check_dbget_f(SALINITY, prof->salinity, NPMAX,
                                       "DBGET_F, salinity");
   prof->n = min_val(prof->n, n);

   /* It is assumed that the CTD data are nearly perfect, so this
      gap filling will only occasionally be called into play, and
      only for small gaps.
   */
   ng = fill_gaps(prof->salinity, prof->n);
   ngaps = ng;
   ng = fill_gaps(prof->temperature, prof->n);
   ngaps = max_val(ngaps, ng);
   ng = fill_gaps(prof->pressure, prof->n);
   ngaps = max_val(ngaps, ng);
   if (ngaps == prof->n)
      prof->n = 0;  /* don't try to use this profile */

   for (i=0; i < prof->n; i++)
   {
      prof->depth[i] = p_depth( prof->pressure[i], prof->position.lat);
      prof->svel[i]  = SVEL( prof->salinity[i], prof->temperature[i],
                              prof->pressure[i]);
      switch (prof->density_var)
      {
         case CTD_SVAN:
            prof->sigma_sva[i]  = SVAN( prof->salinity[i], prof->temperature[i],
                                    prof->pressure[i], &sigma);
            break;
         case CTD_SIGMA_THETA:
            theta = THETA(prof->salinity[i], prof->temperature[i],
                          0.0, prof->pressure[i]);
            SVAN( prof->salinity[i], theta, 0.0, &sigma);
            prof->sigma_sva[i] = sigma;
            break;
         case CTD_POT_SVAN:
            theta = THETA(prof->salinity[i], prof->temperature[i],
                          0.0, prof->pressure[i]);
            prof->sigma_sva[i] = SVAN( prof->salinity[i], theta, 0.0, &sigma);
                        break;
         case CTD_SIGMA_T:
            SVAN( prof->salinity[i], prof->temperature[i], 0.0, &sigma);
            prof->sigma_sva[i]  = sigma;
		        break;
         case CTD_NONE:
            prof->sigma_sva[i] = BADFLOAT;
                        break;
         default:
            printf("\n%d is an invalid density variable specification\n",
                    prof->density_var);
            exit(-1);
      }

   }
   if (prof->n == 0)   /* try the next one */
   {
      printf("\nBad CTD profile.");
      write_ymdhms_time(stdout, &(prof->time));
      printf("  position: %f  %f\n", prof->position.lat,
                                       prof->position.lon);
      DBMOVE(&step, &ierr);
      check_error(ierr, "DBMOVE in get_ctd_profile");
      get_ctd_profile(prof);      /* recursive call */
   }
   return;  /* only exit */
}                       /* get_ctd_profile */

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

   correct_depths()

   The ADCP depth correction factor, to be multiplied times
   the nominal depth, is the average sound speed from the
   transducer to the depth bin, divided by the nominal value
   used in the instrument.  Averaging from the center of the
   first bin, as done here, is close enough.  (After
   89/02/25 the first depth bin is at zero meters, rather
   than the actual position of the first recorded ADCP bin.
   This is generally closer to the transducer depth than the
   center of the first ADCP bin is.)

   Input: a pointer to a ctd profile with sound speed, the
   array of depths to be corrected, and a pointer to the
   number of elements in the depth array.

   Output: the depth array is corrected in place, and the
   number of elements may be reduced if the ctd profile did
   not extend deep enough.

   With each call, this routine checks to see if the ctd
   profile is the same as the last one (same time).  Only if
   it is different does the routine calculate the correction
   factor.

*/

#if PROTOTYPE_ALLOWED
void correct_depths(CTD_PROFILE_TYPE *ctd, float *depth, int *nd_ptr)
#else
void correct_depths(ctd, depth, nd_ptr)
CTD_PROFILE_TYPE *ctd;
float *depth;
int *nd_ptr;
#endif
{
   int n0, nout;
   int i;
   double svel0;
   static  YMDHMS_TIME_TYPE last_time;
   static  float depth_correction[NDMAX];
   static  int nd_ctd; /* 94/02/10 jr+ */

   if (TIMDIF(&last_time, &(ctd->time)))   /* true -> new ctd profile */
   {
      LONG long_0 = 0L;
      DIFTIM(&(ctd->time), &last_time, &long_0);  /* last equals present */
      /* Calculate soundspeed at the nominal ADCP depths. */
      regrid(ctd->depth, ctd->svel, depth, depth_correction, ctd->n, *nd_ptr,
                        &n0, &nout);
      /* If many points are lost, something is probably very wrong. */
      if( (n0 > 5) || ((n0 + nout) < (*nd_ptr-5)))
      {
         printf("\nWARNING: available depth range in correct_depths is");
         printf("\n   %f to %f,  at ctd time ",
                     depth[n0], depth[n0 + nout - 1] );
         write_ymdhms_time(stdout, &(ctd->time));
      }

      /* If there is a gap at the start, assume constant soundspeed */
      for (i=0; i<n0; i++)
         depth_correction[i] = depth_correction[n0];

      /* If there is a gap at the end, shorten the profile. */
      /* *nd_ptr = min_val(*nd_ptr, n0+nout); */
      /* Fri  98/07/17 No; extend the sound speed down. */
      for (i=n0+nout; i<*nd_ptr; i++)
         depth_correction[i] = depth_correction[n0+nout-1];


      svel0 = depth_correction[0];
      integrate_tr(depth_correction, depth, depth_correction, *nd_ptr);
      depth_correction[0] = svel0 / nominal_svel;
      for (i=1; i<*nd_ptr; i++)
      {
         depth_correction[i] /= (nominal_svel * (depth[i] - depth[0]) );
      }
      nd_ctd = *nd_ptr; /* 94/02/10 jr+: remember this for subsequent ensembles! */
   }
   *nd_ptr = min_val(*nd_ptr, nd_ctd); /* jr+: 94/02/10 */
   scale_elements(depth, depth_correction, depth, *nd_ptr); 
}                       /* correct_depths() */

/*  end of CTD module for ADCPSECT
****************************************************************/

