FUNCTION betar,  t, s, p
;
;   compute beta= 1/rho (d rho / d S)
;
   ds = 0.01
   siga = eos(t, s)-1000.
   sigb = eos(t, s+ds)-1000.
   
   be = 0.001*(sigb-siga)/ds/(1.+1.e-3*siga)

   return,  be

END