/*

      SUBROUTINE BVFREQ(Z,S,T,P,G,N,DENS,SV,SNDVEL,BVF)
C===========================================================
C
C    Subroutine BVFREQ
C
C       It generates a Brunt-Vaisala frequency profile from salinity,
C       temperature, pressure and gravity profiles.
C
C   INPUT:      Z --- depth array in m. (depths should be negative)
C               S --- real array of N salinity values in ppt
C               T ---   "    "    " " temperature values in deg C
C               P ---   "    "    " " pressure values in db
C               L --- double, latitude in degrees (for gravity calc.)
C               N --- number of points in profiles
C
C   OUTPUT:
C             BVF --- real array with corresponding N values of the bouyancy
C                     frequency in radians/s
C
C============================================================================

      Modified 88-05-30 by removing sound speed, density, and gravity
      from the argument list.  This permits sound speed and density
      to be allocated dynamically as double precision variables,
      and makes it easier to call the routine.  The rationale is
      that one usually does not really want the other variables as
      output when one is calculating BVF.  Calculating the acceleration
      of gravity internally also saves trouble.

      Fri  02-10-1989  revised, debugged, and checked; alternative
      method of calculation added.

      Fri  08-31-1990  removed error-checking bug in BVFREQ.

      WARNING: THIS VERSION OF BVFREQ HAS NOT BEEN
      ADEQUATELY CHECKED.  THE TWO ALGORITHMS DO NOT EXACTLY
      AGREE, AND I AM NOT SURE THAT EITHER OF THEM IS A GOOD
      ONE TO USE.

*/
#include "common.h"
#include "ocean.h"
#include "dbext.h"
#include "ioserv.h"

/* next function was not originally in this set of routines */
/* K is a ceiling above which the gradient is not allowed to rise */
/* Note that with Z negative down, the density gradient is negative.
   We will assume that increasing index means increasing depth,
   regardless of whether Z is taken as positive or negative down.
   GRAD here uses a simple centered difference scheme, which should
     be OK for typical situations.  If Z is not a uniform grid then
     perhaps a fancier  estimate of the derivative should be used.
*/

void GRAD(double *DENS, float *Z, double *DENSGR, int N)
{
     int i;

     for (i = 1; i < N-1; i++)
          DENSGR[i] = -(DENS[i+1] - DENS[i-1])/fabs(Z[i+1] - Z[i-1]);
     DENSGR[0] = -(DENS[1] - DENS[0])/fabs(Z[1] - Z[0]);
     DENSGR[N-1] = -(DENS[N-1] - DENS[N-2])/fabs(Z[N-1] - Z[N-2]);
     return;
}

/*  input:  Z, S, T, P, L, N */
void BVFREQ(float Z[], float S[], float T[], float P[], double L, int N, float BVF[])
{
     double *DENSGR;
     double BVFI;
     double *SIGMA;
     double sndvel, gi;  /* temporary variables used in a loop */
     int i;

     check_error( ((DENSGR = (double *)(calloc(N, sizeof(double)))) == 0),
                   "BVFREQ; out of memory");
     check_error( ((SIGMA = (double *)(calloc(N, sizeof(double)))) == 0),
                   "BVFREQ; out of memory");

/* ====> OBTAIN DENSITIES AND SOUND SPEEDS */

     for (i = 0; i < N; i++)
     {
        SVAN(S[i],T[i],P[i], SIGMA+i);
     }

/*====> OBTAIN DENSITY GRADIENT */

     GRAD(SIGMA,Z,DENSGR,N);

/*===>  CALCULATE BV FREQUENCIES */

     for (i=0; i<N; i++)
     {
          sndvel = SVEL(S[i],T[i],P[i]);
          gi = GRAV((double)P[i], L);
          BVFI = -gi*( DENSGR[i]/(1000.0 + SIGMA[i]) + gi/(sndvel * sndvel) );
          if (BVFI > 0) BVF[i] = sqrt(BVFI);
          else  BVF[i] = 0.0;
     }
     free(DENSGR);
     free(SIGMA);
     return;
}

double bvf2(double Z0, double T0, double S0, double Z1, double T1, double S1, double lat)
{
   double p0, p1, pm;
   double sigma0, sigma1;
   double bvfsq;

   Z0 = fabs(Z0);
   Z1 = fabs(Z1);
   p0 = PRESS(Z0);
   p1 = PRESS(Z1);
   pm = 0.5 * (p0 + p1);
   SVAN(S0, THETA(S0, T0, p0, pm), pm, &sigma0);
   SVAN(S1, THETA(S1, T1, p1, pm), pm, &sigma1);

   bvfsq = (GRAV( pm, lat) / (1000.0 + sigma0))   /* g/rho */
               * (sigma1 - sigma0) / (Z1 - Z0);         /* d rho / dz */
   if (bvfsq < 0)
   {
      printf("\nwarning: unstable stratification\n");
      return(0.0);
   }
   return(sqrt(bvfsq));
}

#ifdef EXE

void main(int argc, char *argv[])
{
   float T[2], S[2], Z[2];      /* input */
   double L;
   float P[2];                  /* intermediate */
   float B[2];                  /* output */

   if (argc != 8)
   {
      printf("\n USAGE: bvfreq <Z1> <T1> <S1> <Z2> <T2> <S2> <Lat>");
      printf("\n NOTE: Z2 must be deeper than Z1\n\n");
      exit(0);
   }
   else
   {
      Z[0] = atof(argv[1]);
      T[0] = atof(argv[2]);
      S[0] = atof(argv[3]);
      Z[1] = atof(argv[4]);
      T[1] = atof(argv[5]);
      S[1] = atof(argv[6]);
      L    = atof(argv[7]);
   }
   P[0] = PRESS(fabs(Z[0]));
   P[1] = PRESS(fabs(Z[1]));
   BVFREQ(Z, S, T, P, L, 2, B);
   printf("\n BVF  = %e /s",  B[0]);
   printf("\n BVF2 = %le /s",  bvf2(Z[0], T[0], S[0], Z[1], T[1], S[1], L));
   printf("\n for S1 = %g ppt, T1 = %g C, Z1 = %g m (%g dbar),", 
      S[0], T[0], Z[0], P[0]);
   printf("\n     S2 = %g ppt, T2 = %g C, Z2 = %g m (%g dbar),", S[1], T[1], Z[1], P[1]);
   printf("\n at %g degrees latitude\n\n", L);
}
#endif