## baroclinic.m -- computes and plots the baroclinic instability
## dispersion relation in isentropic coordinates

function baroclinic

  kappa = [1:100]/40;
  ef = exp(-2*kappa);
  a = (2 - kappa).*(2 - kappa);
  b = (2 + kappa).*(2 + kappa);
  omsq = (a - b.*ef)./(1 - ef);

  id = ncoct("create", "baro.nc");
  dim = ncoct("putdim", id, "kappa", kappa);
  ncoct("putvar", id, "omegasq" ,-1 , dim, omsq);
  ncoct("close", id);

endfunction
