source: trunk/tools/density_binning/binning_neutral_and_co/bn2.pro

;+
;
; Compute the local Brunt-Vaisala frequency 
;       
; Method :
; -------
;       neos = 0  : UNESCO sea water properties
;          The brunt-vaisala frequency is computed using the polynomial
;       polynomial expression of McDougall (1987):
;               N^2 = g * beta * ( alpha/beta*dk[ tn ] - dk[ sn ] )/e3w
;
;       neos = 1  : linear equation of state (temperature only)
;               N^2 = g * ralpha * dk[ tn ]/e3w
;
;       neos = 2  : linear equation of state (temperature & salinity)
;               N^2 = g * (ralpha * dk[ tn ] - rbeta * dk[ sn ] ) / e3w
;
;       The use of potential density to compute N^2 introduces e r r o r
;      in the sign of N^2 at great depths. We recommand the use of 
;      neos = 0, except for academical studies.
;
;      N.B. N^2 is set to zero at the first level (JK=1) in inidtr
;      and is never used at this level.
;
; References :
; -----------
; McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
;
; @param TN {in} 
; potential temperature
;
; @param SN {in} 
; salinity sn
;
; @returns
; brunt-vaisala frequency
;
; @uses
; <pro>common</pro>
;
; @history
; Original  : 94-07 (G. Madec, M. Imbard)
; Additions : 97-07 (G. Madec) introduction of statement functions
; Idl version : 98-12 (M. Levy)
;
; @version
; $Id$
;
;-
FUNCTION bn2, tn, sn, neos
;
  compile_opt idl2, strictarrsubs
;
@common
;
   case n_params() of
      0:begin
         print, 'rentrer nom des champs de temperature et/ou salinite'
         goto, sortie
      end
      1:begin
         print, 'pas de salinite : valeur de neos imposee a 1'
         neos = 1
      end
      2:begin
         print, 'calcul Unesco par default'
         neos = 0
      END
      3:
   endcase
domdef 
;  coordonnees z : on s'assure qu'on a autant de points W que de points 
;  T sur le domaine defini par domdef.
;  on va se servir dans les calculs de e3w (e3) et gdepw (gdep)
   IF (nzt NE nzw ) THEN BEGIN
      print, 'le domaine doit etre choisi avec autant de points' 
      print, 'sur  niveaux T que sur les niveaux W pour calculer N2'
      z = -1
      GOTO, sortie
   endif
   IF (nzt LE 1 ) THEN BEGIN
      print, 'il faut un minimum de deux points sur la verticale' 
      print, ' pour calculer N2'
      z = -1
      GOTO, sortie
   endif
   cdgrid = 'W'
   g = 9.81
   grille, mask, glam, gphi, gdep, nx, ny,nz, premierx, premiery, premierz, dernierx, derniery, dernierz ,e1,e2,e3, TRI=tri
   zgde3w = reform( (fltarr(nx*ny)+1)#(g/e3), nx, ny, nz)
;
varname='Brunt-Vaisala frequency'
varunit='s-1'
;vargrid='W'
;
; mask tn and sn
        tn=tn*tmask
        sn=sn*tmask
;
;  Interior points ( 2=< jk =< jpkm1 )
;  ----------------
;
   case 1 of
      (neos EQ 0): begin 
; ... UNESCO seawater equation of state
;   ... temperature, salinity anomalie 
         zzt = 0.5*( tn + shift( tn, 0, 0, 1 ))
         zzs = 0.5*( sn + shift( sn, 0, 0, 1) ) - 35.0
         zzp = reform( (fltarr(nx*ny)+1)#gdep, nx , ny, nz)
;   ... ratio alpha/beta and beta
         zalbet = fsalbt( zzt, zzs, zzp )
         zbeta = fsbeta( zzt, zzs, zzp )
;   ... N^2 = g/e3w * beta * ( alpha/beta * dk[tn] - dk[sn])
         z = zgde3w * zbeta * mask $ 
          * ( zalbet * ( shift(tn, 0, 0, 1) - tn ) $ 
              - ( shift(sn, 0, 0, 1) - sn ) )
      end
      ( neos EQ 1 ): begin
; ... First Linear density formulation (function of temperature only)
;
        ralpha = 1.176e-4
        z = zgde3w * ralpha * mask* (shift(tn, 0, 0, 1) - tn)
         
      end
      ( neos EQ 2 ): begin
; ... Second linear density formulation (function of temp. and salinity)
;
         ralpha = 2.e-4
         rbeta = 0.001      
         z = zgde3w * mask $ 
          * (  ralpha * (shift(tn, 0, 0, 1) - tn) $ 
               - rbeta  * (shift(sn, 0, 0, 1) - sn)  )
      end   
      (neos NE 0) AND (neos NE 1) AND (neos NE 2): BEGIN
         print, 'neos doit etre egal a 0, 1 ou 2'
         GOTO, sortie
      end
   endcase
;
;
;  first and last levels
; ------------------------
; bn^2=0. at first vertical point and jk=jpk 
   z[*, *, 0] = 0
   IF (gdep[nz-1] eq gdepw[jpk-1]) THEN z[*, *, nz-1] = 0
sortie:
   return, z
END