;+
;
; computes rho  (in situ volumic mass) 
;
; REFERENCE:
;	Compute the potential volumic mass (Kg/m3) from known potential 
;       temperature and salinity fields using Jackett and McDougall 
;       (1994) equation of state.
;       It is computed directly as a function of
;       potential temperature relative to the surface (the opa tn
;       variable), salt and pressure (assuming no pressure variation
;       along geopotential surfaces, i.e. the pressure p in decibars
;       is approximated by the depth in meters.
;              rhop(t,s) = rho(t,s,0)
;         with pressure                      p        decibars
;              potential temperature         t        deg celsius
;              salinity                      s        psu
;              reference volumic mass        rau0     kg/m**3
;              in situ volumic mass          rho      kg/m**3
;
;         Check value: rho = 1059.8204 kg/m**3 for p=10000 dbar,
;          t = 40 deg celcius, s=40 psu
;
;       Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
;
; @categories
;       calculation
;
; @examples 
;       tableau=eos(t,s,sigma)
;
; @param T {in}{required}
; temperature
;
; @param S  {in}{required}
; salinity
;
; @param SIGMA {in}{required}
; reference level (real number belonging to [0,5])
; for instance, sigma=3. means that rho will be referenced to 3000m  
;
; @uses
; <pro>common</pro>
;
; @history
;       Maurice Imbard
;       Eric Guilyardi - adaptation to post_it
;       Gurvan Madec (04/13/2001) 
; @version
; $Id$
;
;-
FUNCTION eos, t, s, sigma
;
  compile_opt idl2, strictarrsubs
;
@common
;
IF (sigma GT 10) THEN BEGIN
print, 'Caution ! Sigma must be smaller than 5'
ENDIF ELSE BEGIN
	; mask t and s fields
	  zt=t*tmask
	  zs=s*tmask 
	; potential volumic mass
	;
	;   ... depth
	            zh = 1000. * sigma
	            
	;   ... square root salinity
     	
	            zsr= sqrt( abs( zs ) )
	;   ... compute volumic mass pure water at atm pressure
	            zr1= ( ( ( ( 6.536332E-9*zt-1.120083E-6 )*zt+1.001685E-4)*zt  $
	                        -9.095290E-3 )*zt+6.793952E-2 )*zt+999.842594
	;   ... seawater volumic mass atm pressure
	            zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt       $
	                      -4.0899e-3 ) *zt+0.824493
	            zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3
	            zr4= 4.8314e-4
	;
	;   ... masked potential volumic mass (referenced to the surface)
	            zrhop= (( zr4*zs + zr3*zsr + zr2 ) *zs + zr1) * tmask
	
     IF (sigma EQ 0) THEN BEGIN 
     rho=zrhop
     ENDIF ELSE BEGIN
	
	
	;   ... add the compression terms
	            ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6
	            zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4
	            zb = zbw + ze * zs
	;
	            zd = -2.042967e-2
	            zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896
	            zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt		$
	                -4.721788
	            za = ( zd*zsr + zc ) *zs + zaw
	;
	            zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545
	            za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077
	            zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt	$
	                      +2098.925 ) *zt+190925.6
	            zk0= ( zb1*zsr + za1 )*zs + zkw
	;
	;   ... masked in situ volumic mass referenced to sigma level
	            rho = zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) * tmask
	ENDELSE            
ENDELSE
return, rho 
end