subroutine vfreq(grav,prsdb,prvprs,nxtprs,
     1                 prvrho,nxtrho,c,vfreq2)
c
c     title:
c     *****
c
c       vfreq  -- calculate brunt-vaisala frequency squared
c
c     system:
c     ******
c
c       pacodf hydrographic data library
c
c     purpose:
c     *******
c
c       to calculate brunt-vaisala frequency squared
c
c     method:
c     ******
c
c       reference:
c
c          eckart, c., 1962. the sea. interscience publishers.
c          m. n. hill general editor. chapter 2, page 37, equation 35.
c
c       v**2 = g**2*(drho/(dp*10)-1/c**2)*1.0e6
c
c       where:    g   = acceleration of gravity m/s**2
c                 drho= insitu density gradient tons/m**3
c                 dp  = pressure gradient in db
c                 c   = sound velocity m/s
c                                                      6
c                 v   = brunt-vaisala frequency hz x 10
c
c       brunt-vaisala frequency is calculated over an interval
c       centered around the pressure corresponding to the sound
c       velocity; i.e., the density gradient and pressure gradient
c       are calculated from a previous pressure to a next pressure.
c
c     usage:
c     *****
c
c       the first and last intervals calculated will have no brunt-vaisala
c       frequency, as there are respectively no previous or next density
c       or pressure values defined.
c
c     external references:
c     ******** **********
c
c       sin    -- intrinsic fortran trig function
c
c     parameters:
c     **********
c
c       grav    -> local gravity @ insitu pressure (m/sec/sec)
c       prsdb   -> pressure in decibars
c       prvprs  -> previous pressure in decibars (for pressure gradient)
c       nxtprs  -> next pressure in decibars (for pressure gradient)
c       prvrho  -> previous insitu density kg/m**3 @prvprs
c       nxtrho  -> next insitu density kg/m**3 @nxtprs
c       c       -> sound velocity (m/s) @ prsdb
c       vfreq2  <- brunt-vaisala frequency**2
c                  (hz*1.oe6)
c
        real grav,prsdb,prvprs,nxtprs,prvrho,nxtrho,c,vfreq2
c
c     variables:
c     *********
c
        real*8 terma,drho
        real*8 dgrav
c       *** note: this computation must be done in double-precision ***
c
c     constants:
c     *********
c
        integer isnan
c         /* null data value */
c
c     code:
c     ****
c
c     /* calculate brunt-vaisala frequency */
c
c     select
c       (no previous pressure or no next pressure):
          if(isnan(prvprs).ne.0.or.isnan(nxtprs).ne.0) then
          call setnan(vfreq2)
          go to 999
          endif
c       (otherwise):
   10     continue
          dgrav = dble(grav)
          drho = (dble(nxtrho)/1000.0d0-dble(prvrho)/1000.0d0)
c           /* density gradient in tons/m**3 */
          terma  = dgrav*dgrav*(drho/(dble(nxtprs-prvprs)*10.0d0)
     1                        -1.0d0/(dble(c)*dble(c)))*1.0d6
          vfreq2 = terma
c           /* brunt-vaisala frequency squared */
c     end select
c
c     /* return */
c
  999 continue
      return
c
c     end vfreq
c
      end