function eos_nemo, t_in, psal
;------------------------------------------------------------------------------
;
; Purpose:
;       To calculate the density using the equation of state used in NEMO.
;       This is based on the Jackett and McDougall (1994) equation of state
;       for calculating the in situ density based on potential temperature
;       and salinity.
;
; Inputs:
;       temperature    =>  1d array potential temperature (deg C)
;       salinity       =>  1d array salinity (PSU)
;
; Outputs:
;       rhop            =>  1d array pot density (kg/m**3)
;
; could add use a difference reference pressure
; observations are in-situ temperature?
;
; Author:
;       D. J. Lea. Dec 2006.
;------------------------------------------------------------------------------

          zws = SQRT( ABS( psal ) )

	sz=size(t_in)
	ndim=sz(0)
	NO_LEVS=sz(ndim)
	
	prhop=psal*0.

;                 
          for jk = 0L, NO_LEVS-1 do begin  

                   zt = T_IN(jk)
                   zs = psal(jk)
; * depth
;                   zh = O_DEP_LEVS(jk)		;used in calculating insitu density only
                   zsr= zws(jk)
; * 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

; * potential volumic mass (reference to the surface)
                   zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1

; * save potential volumic mass
                   prhop(jk) = zrhop

; * 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

;                   ; in situ density anomaly
;                   prd(jk) = zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )
;
          endfor  



return, prhop

END