FUNCTION alpha,  t, s, p
;
;   compute alpha=-1/rho (d rho / d T)
;
   dt = 0.05
   siga = eos(t, s)-1000.
   sigb = eos(t+dt, s)-1000.
   
   al = -0.001*(sigb-siga)/dt/(1.+1.e-3*siga)

   return,  al

END