## etendency1.m -- Compute the bulk energy tendency of the earth
## for variable albedo as a function of temp characteristic of ice ages.

function etendency2

  ## set constants
  fs = 1370.;
  sigma = 5.67e-8;
  tcc = 5418.;     # tcc = L/R_v
  tr = 288.;
  greenhouse = 33;
  calbedo = 0.6;
  walbedo = 0.3;
  malbedo = 0.5*(calbedo + walbedo);
  dalbedo = (calbedo - walbedo)/3.14159;
  salbedo = 5;

  ## plot prelims
  closeplot;
  hold off;
  grid;
  xlabel("ocean temperature (K)");
  ylabel("net flux imbalance (W/m^2)");
  title("Cold albedo = 0.6; warm albedo = 0.3; Fixed greenhouse gain");

  ## compute tendency
  max = 59;
  oceantemp0 = 240;
  doceantemp = 1;
  oceantemp = oceantemp0 + doceantemp*[1:max];
  albedo = malbedo - dalbedo*atan((oceantemp - 273)/salbedo);
  tendin = fs*(1 - albedo)/4;
  trad = oceantemp - greenhouse;
  for i = 1:max
    if (trad(i) < 0)
      trad(i) = 0;
    endif
  endfor
  tendout = sigma*trad.^4;
  tendency = tendin - tendout;

  ## plot results
  gset nokey;
  plot(oceantemp,tendency);
  hold on;
  plot([240,300],[0,0]);

  ## do hard copy
  gset terminal postscript eps;
  gset output "etendency2.eps";
  replot;

endfunction
