Documentation of plot_climatology

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

  Plot the climatology for the backward rotating earth run.
  This set of files is for the T42 version

Cross-Reference Information

Listing of script plot_climatology

clean

filin1 = 'earth3.nc';
filin2 = 'htrae3.nc';

%  Set up data
cd /home/disk/hayes2/dvimont/ccm/ccm3.6/run/sun/HT/data
nc1 = netcdf(filin1, 'nowrite');
nc2 = netcdf(filin2, 'nowrite');

lat = nc1{'lat',1}(:);
lon = nc1{'lon',1}(:);
hyam = nc1{'hyam',1}(:);
hybm = nc1{'hybm',1}(:);
P0 = nc1{'P0',1}(:);

ps1 = nc1{'PS',1}(:,:,:);
ps2 = nc2{'PS',1}(:,:,:);

%  Get variables to plot
newlev = [1000:-50:50];

z1 = nc1{'Z3',1}(:,:,:,:);
z2 = nc2{'Z3',1}(:,:,:,:);
var1 = atlev(z1, newlev, ps1, hyam, hybm, P0);
var2 = atlev(z2, newlev, ps2, hyam, hybm, P0);
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
save z_atlev_T42_2.mat var1 var2 lat lon newlev

z1 = nc1{'U',1}(:,:,:,:);
z2 = nc2{'U',1}(:,:,:,:);
var1 = atlev(z1, newlev, ps1, hyam, hybm, P0);
var2 = atlev(z2, newlev, ps2, hyam, hybm, P0);
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
save u_atlev_T42_2.mat var1 var2 lat lon newlev

z1 = nc1{'V',1}(:,:,:,:);
z2 = nc2{'V',1}(:,:,:,:);
var1 = atlev(z1, newlev, ps1, hyam, hybm, P0);
var2 = atlev(z2, newlev, ps2, hyam, hybm, P0);
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
save v_atlev_T42_2.mat var1 var2 lat lon newlev

z1 = nc1{'OMEGA',1}(:,:,:,:);
z2 = nc2{'OMEGA',1}(:,:,:,:);
var1 = atlev(z1, newlev, ps1, hyam, hybm, P0);
var2 = atlev(z2, newlev, ps2, hyam, hybm, P0);
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
save omega_atlev_T42.mat var1 var2 lat lon newlev

z1 = nc1{'AIRT',1}(:,:,:,:);
z2 = nc2{'AIRT',1}(:,:,:,:);
var1 = atlev(z1, newlev, ps1, hyam, hybm, P0);
var2 = atlev(z2, newlev, ps2, hyam, hybm, P0);
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
save airt_atlev_T42.mat var1 var2 lat lon newlev

z1 = nc1{'VT',1}(:,:,:,:);
z2 = nc2{'VT',1}(:,:,:,:);
var1 = atlev(z1, newlev, ps1, hyam, hybm, P0);
var2 = atlev(z2, newlev, ps2, hyam, hybm, P0);
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
save vt_atlev_T42_2.mat var1 var2 lat lon newlev

nc1 = close(nc1);
nc2 = close(nc2);



%  Plot data
clean

color_trans; colormap(cmap);

default_global;  FRAME = [0 360 -90 90];
figure(1); fo(1); clf;
sz = get(gcf, 'PaperSize'); 
hw = 1/sz(1); vw = 1/sz(2); hsz = hw*(6.5); vsz = vw*(3.5);
bmarg = (sz(2) - 8.5)*0.5; tmarg = (sz(2) - 8)*0.5;
midv = 0.5; midh = 0.5;

cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
load z_atlev_T42.mat; var1a = var1; var2a = var2;
load z_atlev_T42_2.mat;
var1 = 0.5*(var1+var1a); var2 = 0.5*(var2+var2a);
[nmo, nlev, nlat, nlon] = size(var1);

lind = find(newlev == 500); cint = 60; tit = '500 mb Height';
cint2 = 5850; elim = 25; slim = [-150 150];

lind = find(newlev == 250); cint = 120; tit = '250 mb Height';
cint2 = 10860; elim = 60; slim = [-200 200];

win = [12 1 2]; mon = 'DJF';
win = [3 4 5]; mon = 'MAM';
win = [6 7 8]; mon = 'JJA';
win = [9 10 11]; mon = 'SON';

get_global; dg(lat, lon);
FRAME = [0 360 -90 90];
figure(1); fo(1); clf;

tem = squeeze(mean(squeeze(var1(win,lind,:,:))));
tem2 = (tem' - ones(nlon,1)*mean(tem'))';
tem2(abs(tem2)<elim) = NaN; tem2 = [tem2 tem2(:,1)];
lat2 = YAX*ones(1, nlon); lat2 = [lat2 lat2(:,1)];
lon2 = ones(nlat, 1)*XAX'; lon2 = [lon2 lon2(:,1)+360];

FRAME = [0 360 -90 90];
subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg)-vsz hsz vsz]); cla
%  [c,h] = mcont3(tem, cint, 'stereo', [90 0]);
  [c,h] = mcont2(tem, cint, 'giso');
  hold on;  
    surfacem(lat2, lon2, tem2, -0.5*ones(nlat, nlon+1));
    [c2, h2] = contorm(lat2, lon2, [tem tem(:,1)], ...
		       [0 cint2 12000], '--k');
  hold off;
  caxis([slim]);
  dcmfill(-1, 0.8);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 30);
  cl = clabelm(c, h, [0:4*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  cl2 = clabelm(c2, h2, [cint2]);
  set(cl2, 'fontsize', 8);
  set(gca, 'visible', 'on', 'linewidth', 2);
  t1 = title('EARTH', 'visible', 'on', 'fontsize', 10);
  xx = get(gca, 'XLim'); yy = get(gca, 'YLim');
  t2 = text(mean(xx), yy(2)+.5, ['\bf ' mon ':  ' tit], ...
	    'horizontalalignment', 'center', 'fontsize', 14);
  
tem = fliplr(squeeze(mean(squeeze(var2(win,lind,:,:)))));
tem2 = (tem' - ones(nlon,1)*mean(tem'))';
tem2(abs(tem2)<elim) = NaN; tem2 = [tem2 tem2(:,1)];
lat2 = YAX*ones(1, nlon); lat2 = [lat2 lat2(:,1)];
lon2 = ones(nlat, 1)*XAX'; lon2 = [lon2 lon2(:,1)];

FRAME = [0 360 -90 90];
subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg-midv)-2*vsz hsz vsz]); cla
%  [c,h] = mcont3(tem, cint, 'stereo', [90 0]);
  [c,h] = mcont2(tem, cint, 'giso');
  hold on;  
    surfacem(lat2, lon2, tem2, -0.5*ones(nlat, nlon+1));
    [c2, h2] = contorm(lat2, lon2, [tem tem(:,1)], ...
		       [0 cint2 12000], '--k');
  hold off;
  caxis([slim]);
  dcmfill(-1, 0.8);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 30);
  cl = clabelm(c, h, [0:4*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  cl2 = clabelm(c2, h2, [cint2]);
  set(cl2, 'fontsize', 8);
  set(gca, 'visible', 'on', 'linewidth', 2);
  t1 = title(['HTRAE'], 'visible', 'on', 'fontsize', 10);

subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg-2.5*midv)-2*vsz hsz 0.5*midv*vw]); cla
  cb = colorbar(gca);
  set(gca, 'XTick', 1:2:21, 'XTickLabel', ...
	   slim(1):(diff(slim)/10):slim(2), 'fontsize', 8);
  
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Figs2
eval(['print -dpsc2 ' tit(1:3) 'mb_height_and_eddy_' mon '.ps']);
%print -dpsc2 250mb_height_and_eddy_djf.ps
%print -dpsc2 500mb_height_and_eddy_djf.ps






%  Plot data
clean

color_trans; colormap(cmap);

default_global;  FRAME = [0 360 -90 90];
figure(1); fo(1); clf;
sz = get(gcf, 'PaperSize'); 
hw = 1/sz(1); vw = 1/sz(2); hsz = hw*(6.5); vsz = vw*(3.5);
bmarg = (sz(2) - 8.5)*0.5; tmarg = (sz(2) - 8)*0.5;
midv = 0.5; midh = 0.5;

cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
load u_atlev_T42.mat; var1a = var1; var2a = var2;
load u_atlev_T42_2.mat;
var1 = 0.5*(var1+var1a); var2 = 0.5*(var2+var2a);
[nmo, nlev, nlat, nlon] = size(var1);

lind = find(newlev == 500); cint = 5; tit = '500 mb Zonal Wind';
cint2 = 5850; elim = 10; slim = [-30 30];

lind = find(newlev == 250); cint = 120; tit = '250 mb Height';
cint2 = 10860; elim = 60; slim = [-200 200];

win = [12 1 2]; mon = 'DJF';
win = [3 4 5]; mon = 'MAM';
win = [6 7 8]; mon = 'JJA';
win = [9 10 11]; mon = 'SON';

get_global; dg(lat, lon);
FRAME = [0 360 -90 90];
figure(1); fl(1); clf;

tem = squeeze(mean(squeeze(var1(win,lind,:,:))));
tem2 = tem;
tem2(abs(tem2)<elim) = NaN; tem2 = [tem2 tem2(:,1)];
lat2 = YAX*ones(1, nlon); lat2 = [lat2 lat2(:,1)];
lon2 = ones(nlat, 1)*XAX'; lon2 = [lon2 lon2(:,1)+360];

FRAME = [0 360 -90 90];
subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg)-vsz hsz vsz]); cla
%  [c,h] = mcont3(tem, cint, 'stereo', [90 0]);
  [c,h] = mcont2(tem, cint, 'giso');
  hold on;  
    surfacem(lat2, lon2, tem2, -0.5*ones(nlat, nlon+1));
    [c2, h2] = contorm(lat2, lon2, [tem tem(:,1)], ...
		       [0 cint2 12000], '--k');
  hold off;
  caxis([slim]);
  dcmfill(-1, 0.8);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 30);
  cl = clabelm(c, h, [0:4*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  cl2 = clabelm(c2, h2, [cint2]);
  set(cl2, 'fontsize', 8);
  set(gca, 'visible', 'on', 'linewidth', 2);
  t1 = title('EARTH', 'visible', 'on', 'fontsize', 10);
  xx = get(gca, 'XLim'); yy = get(gca, 'YLim');
  t2 = text(mean(xx), yy(2)+.5, ['\bf ' mon ':  ' tit], ...
	    'horizontalalignment', 'center', 'fontsize', 14);
  
tem = -fliplr(squeeze(mean(squeeze(var2(win,lind,:,:)))));
tem2 = tem;
tem2(abs(tem2)<elim) = NaN; tem2 = [tem2 tem2(:,1)];
lat2 = YAX*ones(1, nlon); lat2 = [lat2 lat2(:,1)];
lon2 = ones(nlat, 1)*XAX'; lon2 = [lon2 lon2(:,1)];

FRAME = [0 360 -90 90];
subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg-midv)-2*vsz hsz vsz]); cla
%  [c,h] = mcont3(tem, cint, 'stereo', [90 0]);
  [c,h] = mcont2(tem, cint, 'giso');
  hold on;  
    surfacem(lat2, lon2, tem2, -0.5*ones(nlat, nlon+1));
    [c2, h2] = contorm(lat2, lon2, [tem tem(:,1)], ...
		       [0 cint2 12000], '--k');
  hold off;
  caxis([slim]);
  dcmfill(-1, 0.8);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 30);
  cl = clabelm(c, h, [0:4*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  cl2 = clabelm(c2, h2, [cint2]);
  set(cl2, 'fontsize', 8);
  set(gca, 'visible', 'on', 'linewidth', 2);
  t1 = title(['HTRAE'], 'visible', 'on', 'fontsize', 10);

subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg-2.5*midv)-2*vsz hsz 0.5*midv*vw]); cla
  cb = colorbar(gca);
  set(gca, 'XTick', 1:2:21, 'XTickLabel', ...
	   slim(1):(diff(slim)/10):slim(2), 'fontsize', 8);
  
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Figs2
eval(['print -dpsc2 ' tit(1:3) 'mb_Zonal_Wind_' mon '.ps']);
%print -dpsc2 250mb_height_and_eddy_djf.ps
%print -dpsc2 500mb_height_and_eddy_djf.ps

  
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Figs2
eval(['print -dpsc2 ' tit(1:3) 'mb_height_and_eddy_' mon '.ps']);
%print -dpsc2 250mb_height_and_eddy_djf.ps
%print -dpsc2 500mb_height_and_eddy_djf.ps






%  Eddy heights:  
cint = 20;
nlat = length(YAX); nlon = length(XAX);
figure(1); fo(1); clf;
subplot(2,1,1); cla;
  tem = squeeze(mean(squeeze(var1(win,lind,:,:))));
  tem = tem - mean(tem')'*ones(1,nlon);
  [c,h] = mcont2(tem, cint, 'giso');
  dcmfill(-1, 0.7);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 45);
  cl = clabelm(c, h, [0:4*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  set(gca, 'visible', 'off');
  t1 = title([mon ':  Eddy ' tit ' - EARTH'], 'visible', 'on');
  
subplot(2,1,2); cla;
  tem = fliplr(squeeze(mean(squeeze(var2(win,lind,:,:)))));
  tem = tem - mean(tem')'*ones(1,nlon);
  [c,h] = mcont2(tem, cint, 'giso');
  dcmfill(-1, 0.7);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 45);
  cl = clabelm(c, h, [0:4*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  set(gca, 'visible', 'off');
  t1 = title([mon ':  Eddy ' tit ' - HTRAE'], 'visible', 'on');

  

%  Plot data:  Wind
win = [12 1 2];
get_global; dg(lat, lon);
FRAME = [0 360 -90 90];

cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
load u_atlev_T42.mat; var1u = var1; var2u = var2;
load u_atlev_T42_2.mat;
var1u = 0.5*(var1u+var1); var2u = 0.5*(var2u+var2);
load v_atlev_T42.mat; var1a = var1; var2a = var2;
load v_atlev_T42_2.mat;
var1 = 0.5*(var1+var1a); var2 = 0.5*(var2+var2a);

lind = find(newlev == 950); cint = 2.5; tit = '500 mb Height';
%lind = find(newlev == 500); cint = 5; tit = '500 mb Height';
%lind = find(newlev == 250); cint = 100; tit = '250 mb Height';

figure(1); fo(1); clf;
lat1 = lat*ones(1, nlon);
lon1 = ones(nlat,1)*lon';
subplot(2,1,1); cla;
  tem = squeeze(mean(squeeze(var1u(win,lind,:,:))));
  tema = squeeze(mean(squeeze(var1(win,lind,:,:))));
  [c,h] = mcont2(tem, cint, 'giso');
  hold on;
    [hh] = quiverm(thin(lat1, 3), thin(lon1, 3), ...
		       thin(tema, 3), thin(tem, 3), 1.5);
  hold off;
  set(hh, 'color', 'k');
  set(h, 'color', 0.5*[1 1 1], 'linewidth', 2);
  dcmfill(-1, 0.8);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 45);
  cl = clabelm(c, h, [0:4*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  set(gca, 'visible', 'off');
  t1 = title(['DJF:  ' tit ' - EARTH'], 'visible', 'on');
  
subplot(2,1,2); cla;
  tem = -fliplr(squeeze(mean(squeeze(var2u(win,lind,:,:)))));
  tema = fliplr(squeeze(mean(squeeze(var2(win,lind,:,:)))));
  [c,h] = mcont2(-tem, cint, 'giso');
  hold on;
    [hh] = quiverm(thin(lat1, 3), thin(lon1, 3), ...
		       thin(tema, 3), thin(tem, 3), 1.5);
  hold off;
  set(hh, 'color', 'k');
  set(h, 'color', 0.5*[1 1 1], 'linewidth', 2);
  dcmfill(-1, 0.8);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 45);
  cl = clabelm(c, h, [0:4*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  set(gca, 'visible', 'off');
  t1 = title(['DJF:  ' tit ' - HTRAE'], 'visible', 'on');


%  Precip
var1 = nc1{'PRECL',1}(:,:,:)*24*3600*1000;
var1a = nc1{'PRECC',1}(:,:,:)*24*3600*1000;
var2 = nc2{'PRECL',1}(:,:,:)*24*3600*1000;
var2a = nc2{'PRECC',1}(:,:,:)*24*3600*1000;
win = [6 7 8];
tit = 'Precipitation (convective contoured, large scale shaded)'
cint = 2;

figure(2); fo(1); clf;
dg(lat, lon);
nlat = length(YAX); nlon = length(XAX);
lat1 = (lat-mean(diff(YAX))/2)*ones(1, nlon);
lon1 = ones(nlat,1)*(lon-mean(diff(XAX))/2)';
subplot(2,1,1); cla;
  tem = squeeze(mean(var1(win,:,:)));
  tema = squeeze(mean(var1a(win,:,:)));
  [c,h] = mcont2(tema, cint, 'giso');
  hold on;
    tem(tem < 1) = NaN;
    [hh, cc] = surfacem(lat1, lon1, tem);
  hold off;
  zd = get(hh, 'ZData'); set(hh, 'ZData', zd-0.1);
  caxis([-5 5]);
  dcmfill(-1, 0.7);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 45);
  cl = clabelm(c, h, [0:cint:max(max(tema))]);
  set(cl, 'fontsize', 8);
  set(gca, 'visible', 'off');
  t1 = title(['DJF:  ' tit ' - EARTH'], 'visible', 'on');
  
subplot(2,1,2); cla;
  tem = fliplr(squeeze(mean(var2(win,:,:))));
  tema = fliplr(squeeze(mean(var2a(win,:,:))));
  [c,h] = mcont2(tema, cint, 'giso');
  hold on;
    tem(tem < 1) = NaN;
    [hh, cc] = surfacem(lat1, lon1, tem);
  hold off;
  zd = get(hh, 'ZData'); set(hh, 'ZData', zd-0.1);
  caxis([-5 5]);
  dcmfill(-1, 0.7);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 45);
  cl = clabelm(c, h, [0:cint:max(max(tema))]);
  set(cl, 'fontsize', 8);
  set(gca, 'visible', 'off');
  t1 = title(['DJF:  ' tit ' - HTRAE'], 'visible', 'on');
  

%  Precip
var1 = nc1{'PSL',1}(:,:,:)/100;
var2 = nc2{'PSL',1}(:,:,:)/100;
win = [12 1 2];
tit = 'SLP'
cint = 5;

figure(2); fo(1); clf;
dg(lat, lon);
nlat = length(YAX); nlon = length(XAX);
lat1 = (lat-mean(diff(YAX))/2)*ones(1, nlon);
lon1 = ones(nlat,1)*(lon-mean(diff(XAX))/2)';
subplot(2,1,1); cla;
  tem = squeeze(mean(var1(win,:,:)));
  [c,h] = mcont2(tem, cint, 'giso');
  dcmfill(-1, 0.7);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 45);
  cl = clabelm(c, h, [0:2*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  set(gca, 'visible', 'off');
  t1 = title(['DJF:  ' tit ' - EARTH'], 'visible', 'on');
  
subplot(2,1,2); cla;
  tem = fliplr(squeeze(mean(var2(win,:,:))));
  [c,h] = mcont2(tem, cint, 'giso');
  dcmfill(-1, 0.7);
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 45);
  cl = clabelm(c, h, [0:2*cint:max(max(tem))]);
  set(cl, 'fontsize', 8);
  set(gca, 'visible', 'off');
  t1 = title(['DJF:  ' tit ' - HTRAE'], 'visible', 'on');
  

  
%  Playing around
clean

cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
load u_atlev_T42.mat; var1u = var1; var2u = var2;
load v_atlev_T42.mat;
lind = find(newlev == 500); cint = 5; tit = '500 mb Height';
%lind = find(newlev == 250); cint = 100; tit = '250 mb Height';

win = [12 1 2];
get_global; dg(lat, lon);
FRAME = [0 360 -90 90];

spd = shiftdim(squeeze(mean(var1u(win,:,:,:))), 1);

[x, y, z] = meshgrid(lon, lat, -newlev+1000);
global GRDX_SPACING; GRDX_SPACING = 45;
global GRDY_SPACING; GRDY_SPACING = 30;


figure(1); fl(1); clf;
subplot(1,1,1);
  pp = patch(isosurface(x,y,z,spd,35));
  isonormals(x,y,z,spd, pp);
  set(pp, 'FaceColor', 'red', 'EdgeColor', 'none');
  pp2 = patch(isocaps(x,y,z,spd,35))
  set(pp2, 'FaceColor', 'interp', 'EdgeColor', 'none');
  
  daspect([1.5 1 8])
  view([-27.5 20])
  axis([0 360 -90 90 0 1000]);
  camlight(0, -10); lighting phong
  box on
  camproj perspective;
  hh = dc;
  set(gca, 'ZTick', [0:250:1000])
  set(gca, 'ZTickLabel', 1000:-250:0)
  
load z_atlev_T42.mat;
lind = find(newlev == 500);
  hold on;
  [cc2, hh2] = contour(lon, lat, squeeze(mean(var1(win,lind,:,:))), ...
		       4000:150:6000);
  set(hh2, 'edgecolor', 'k', 'linewidth', 2);
  
  
  
clean

cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
load u_atlev_T42.mat; var1u = var1; var2u = var2;
load z_atlev_T42.mat; var1z = var1; var2z = var2;
load v_atlev_T42.mat;
lind = find(newlev == 500); cint = 5; tit = '500 mb Height';
%lind = find(newlev == 250); cint = 100; tit = '250 mb Height';

win = [12 1 2];
get_global; dg(lat, lon);
FRAME = [0 360 -90 90];

spd = shiftdim(squeeze(mean(var1u(win,:,:,:))), 1);

global GRDX_SPACING; GRDX_SPACING = 45;
global GRDY_SPACING; GRDY_SPACING = 30;

tem1 = shiftdim(squeeze(mean(var1u(win,:,:,:))), 1);
tem2 = shiftdim(squeeze(mean(var1(win,:,:,:))), 1);
tem3 = shiftdim(squeeze(mean(var1z(win,:,:,:))), 1);

[nlat, nlon, nlev] = size(tem1);
latthin = 4; lonthin = 8; levthin = 5;
tem1 = tem1(3:latthin:nlat, 3:lonthin:nlon, 1:levthin:nlev);
tem2 = tem2(3:latthin:nlat, 3:lonthin:nlon, 1:levthin:nlev);
tem3 = tem3(3:latthin:nlat, 3:lonthin:nlon, 1:levthin:nlev);
lat2 = lat(1:latthin:nlat);
lon2 = lon(1:lonthin:nlon);
lev2 = newlev(1:levthin:nlev);
[x, y, z] = meshgrid(lon2, lat2, 1000-lev2);
spd = sqrt(tem1.^2+tem2.^2);
tem = ones(size(spd));
tem(spd < 3) = NaN;

figure(1); fl(1); clf;
  daspect([1.5, 1, 8]);
  hh = coneplot(x, y, z, ...
		tem.*tem1, tem.*tem2, tem.*zeros(size(tem2)).*tem2, ...
		x, y, z, tem.*spd);
  set(hh, 'EdgeColor', 'none');
  view([0 80])
  axis([0 360 -90 90 0 1000]);
  camlight(0, -10); lighting phong
  box on
  camproj perspective;
  hh2 = dc;

  
  
  
%  Surface winds
%  Plot data
clean

colormap default

default_global;  FRAME = [0 360 -90 90];
figure(1); fo(1); clf;
sz = get(gcf, 'PaperSize'); 
hw = 1/sz(1); vw = 1/sz(2); hsz = hw*(6.5); vsz = vw*(3.5);
bmarg = (sz(2) - 8.5)*0.5; tmarg = (sz(2) - 8)*0.5;
midv = 0.5; midh = 0.5;

cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
load u_atlev_T42.mat; var1a = var1; var2a = var2;
load u_atlev_T42_2.mat;
var1u = 0.5*(var1+var1a); var2u = 0.5*(var2+var2a);
[nmo, nlev, nlat, nlon] = size(var1);

cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Data2
load v_atlev_T42.mat; var1a = var1; var2a = var2;
load v_atlev_T42_2.mat;
var1v = 0.5*(var1+var1a); var2v = 0.5*(var2+var2a);
clear var1a var2a

%  Get precip
cd /home/disk/hayes2/dvimont/ccm/ccm3.6/run/sun/HT/data
filin1 = 'earth2.nc'; filin2 = 'earth3.nc';
nc1 = netcdf(filin1, 'nowrite'); nc2 = netcdf(filin2, 'nowrite');
var1 = 0.5*(nc1{'PRECL',1}(:,:,:)+nc2{'PRECL',1}(:,:,:)+ ...
	nc1{'PRECC',1}(:,:,:)+nc2{'PRECC',1}(:,:,:))*24*3600*1000;
nc1 = close(nc1); nc2 = close(nc2);
filin1 = 'htrae2.nc'; filin2 = 'htrae3.nc';
nc1 = netcdf(filin1, 'nowrite'); nc2 = netcdf(filin2, 'nowrite');
var2 = 0.5*(nc1{'PRECL',1}(:,:,:)+nc2{'PRECL',1}(:,:,:)+ ...
	nc1{'PRECC',1}(:,:,:)+nc2{'PRECC',1}(:,:,:))*24*3600*1000;
nc1 = close(nc1); nc2 = close(nc2);

lind = find(newlev == 900); tit = '900 mb Winds and Precip';

win = [12 1 2]; mon = 'DJF';
win = [3 4 5]; mon = 'MAM';
win = [6 7 8]; mon = 'JJA';
win = [9 10 11]; mon = 'SON';

get_global;
FRAME = [0 360 -40 40];
figure(1); fo(1); clf;

nthin = 2; dg(lat, lon);
temx = thin(XAX2, nthin); temy = thin(YAX2, nthin);
temu = thin(squeeze(mean(squeeze(var1u(win,lind,:,:)))), nthin);
temv = thin(squeeze(mean(squeeze(var1v(win,lind,:,:)))), nthin);
temp = squeeze(mean(var1(win,:,:))); temp = [temp temp(:,1)];
XAX2 = [XAX2 XAX2(:,1)]; YAX2 = [YAX2 YAX2(:,1)];
tem2 = temp; tem2(abs(tem2)<4)=NaN;

subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg)-vsz hsz vsz]); cla
  maxes('giso', [0 30]);
  h = quiverm(temy, temx, temv, temu);
  hold on;  
    surfacem(YAX2, XAX2, tem2, -0.5*ones(nlat, nlon+1));
  hold off;
  caxis([0 15]);
  dcmfill(-1, 0.8);
  set(h, 'color', 'k');
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 30);
  set(gca, 'visible', 'on', 'linewidth', 2);
  t1 = title('EARTH', 'visible', 'on', 'fontsize', 10);
  xx = get(gca, 'XLim'); yy = get(gca, 'YLim');
  t2 = text(mean(xx), yy(2)+.5, ['\bf ' mon ':  ' tit], ...
	    'horizontalalignment', 'center', 'fontsize', 14);
  
dg(lat, lon);
temx = thin(XAX2, nthin); temy = thin(YAX2, nthin);
temu = -fliplr(thin(squeeze(mean(squeeze(var2u(win,lind,:,:)))), nthin));
temv = fliplr(thin(squeeze(mean(squeeze(var2v(win,lind,:,:)))), nthin));
temp = fliplr(squeeze(mean(var2(win,:,:)))); temp = [temp temp(:,1)];
XAX2 = [XAX2 XAX2(:,1)]; YAX2 = [YAX2 YAX2(:,1)];
tem2 = temp; tem2(abs(tem2)<4)=NaN;

subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg-midv)-2*vsz hsz vsz]); cla
  maxes('giso', [0 30]);
  h = quiverm(temy, temx, temv, temu);
  hold on;  
    surfacem(YAX2, XAX2, tem2, -0.5*ones(nlat, nlon+1));
  hold off;
  caxis([0 15]);
  dcmfill(-1, 0.8);
  set(h, 'color', 'k');
  gridm on;
  axis_limits(100);
  setm(gca, 'plinelocation', 20, 'mlinelocation', 30);
  set(gca, 'visible', 'on', 'linewidth', 2);
  t1 = title('HTRAE', 'visible', 'on', 'fontsize', 10);

subplot('position', [0.5-hsz/2 ...
        vw*(sz(2)-tmarg-2.5*midv)-2*vsz hsz 0.5*midv*vw]); cla
  cb = colorbar(gca);
  set(gca, 'XTick', 1:(64/15):65, 'XTickLabel', ...
	   0:15, 'fontsize', 8);
  
cd /home/disk/tao/dvimont/matlab/CCM/Htrea/Figs2
eval(['print -dpsc2 ' tit(1:3) 'mb_Zonal_Wind_' mon '.ps']);
%print -dpsc2 250mb_height_and_eddy_djf.ps
%print -dpsc2 500mb_height_and_eddy_djf.ps