MITgcm_contrib/bbl/lookat_output.m

% cd ~/mitgcm/bbl/MITgcm
lon=0:.48:11.6; lat=-70:.16:-67.7; dpt=72:144:3600;
Depth=readbin('run/Depth.data',[25 15]);
colormap(cmap), cx=[29.99 30.1];

% find bottom
kbot=ones(25,15);
s1=35+readbin(['run.bbl/SALTanom.' myint2str(216,10) '.data'],[25 15 25]);
for k=2:25
 kbot(find(s1(:,:,k)>0))=k;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compare baseline, downslope, and bbl
colormap(cmap), cx=[29.99 30.1];
s1=zeros(25,15,25); s2=s1; s3=s1;
for ts=216:216:7200
 s1(:,:,1:25)=35+readbin(['run/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
 s2(:,:,1:25)=35+readbin(['run.down_slope/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
 s3(:,:,1:25)=35+readbin(['run.bbl/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
 clf, subplot(311)
 mypcolor(lat,-dpt/1e3,squeeze(s1(10,:,:))'); caxis(cx), thincolorbar
 title(['baseline salinity section on day ' int2str(ts*1200/60/60/24)])
 hold on, plot(lat,-Depth(10,:)/1e3,'r.-'), ylabel('depth (km)')
 subplot(312)
 mypcolor(lat,-dpt/1e3,squeeze(s2(10,:,:))'); caxis(cx), thincolorbar
 title(['downslope salinity section on day ' int2str(ts*1200/60/60/24)])
 hold on, plot(lat,-Depth(10,:)/1e3,'r.-'), ylabel('depth (km)')
 subplot(313)
 mypcolor(lat,-dpt/1e3,squeeze(s3(10,:,:))'); caxis(cx), thincolorbar
 title(['bbl salinity section on day ' int2str(ts*1200/60/60/24)])
 hold on, plot(lat,-Depth(10,:)/1e3,'r.-'), ylabel('depth (km)')
 xlabel('latitude (deg)'), pause(.5)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% look at some diagnostics from bbl integration
for ts=216:216:7200
 s1=35+readbin(['run.bbl/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
 clf, subplot(321), mypcolor(lat,-dpt/1e3,squeeze(s1(10,:,:))'); caxis(cx), thincolorbar
 title(['salinity section on day ' int2str(ts*1200/60/60/24)])
 s2=kbot; for i=1:25, for j=1:15, s2(i,j)=s1(i,j,kbot(i,j)); end, end
 subplot(322), mypcolor(s2'); caxis(cx), thincolorbar
 title('salinity of bottom grid box')
 s3=readbin(['run.bbl/BBLsalt.' myint2str(ts,10) '.data'],[25 15]);
 subplot(323), mypcolor(s3'); caxis(cx), thincolorbar
 title('salinity of bbl')
 subplot(324), mypcolor(s3'-s2'); thincolorbar
 title('bbl salinity minus bottom salinity')
 tendS=readbin(['run.bbl/BBLtendS.' myint2str(ts,10) '.data'],[25 15]);
 subplot(325), mypcolor(tendS'); caxis([-1 1]*3e-9), thincolorbar
 title('bottom salinity tendency due to bbl'), pause(.1)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check that tendency integrates to zero
RAC=readbin('run.bbl/RAC.data',[25 15]);
hFacC=readbin('run.bbl/hFacC.data',[25 15 25]);
for i=1:25
 for j=1:15
  RAC(i,j)=RAC(i,j)*hFacC(i,j,kbot(i,j));
 end
end
for ts=216:216:7200
 eta=readbin(['run.bbl/BBLeta.' myint2str(ts,10) '.data'],[25 15]);
 tendS=readbin(['run.bbl/BBLtendS.' myint2str(ts,10) '.data'],[25 15]);
 tmp=RAC.*tendS;
 disp([ts sum(tmp(:))]), pause(1)
end
1	dimitri	1.1	% cd ~/mitgcm/bbl/MITgcm
2	dimitri	1.2	lon=0:.48:11.6; lat=-70:.16:-67.7; dpt=72:144:3600;
3			Depth=readbin('run/Depth.data',[25 15]);
4			colormap(cmap), cx=[29.99 30.1];
5
6			% find bottom
7			kbot=ones(25,15);
8			s1=35+readbin(['run.bbl/SALTanom.' myint2str(216,10) '.data'],[25 15 25]);
9			for k=2:25
10			kbot(find(s1(:,:,k)>0))=k;
11			end
12	dimitri	1.1
13			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14			% compare baseline, downslope, and bbl
15	dimitri	1.2	colormap(cmap), cx=[29.99 30.1];
16			s1=zeros(25,15,25); s2=s1; s3=s1;
17	dimitri	1.3	for ts=216:216:7200
18	dimitri	1.2	s1(:,:,1:25)=35+readbin(['run/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
19			s2(:,:,1:25)=35+readbin(['run.down_slope/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
20			s3(:,:,1:25)=35+readbin(['run.bbl/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
21	dimitri	1.1	clf, subplot(311)
22	dimitri	1.2	mypcolor(lat,-dpt/1e3,squeeze(s1(10,:,:))'); caxis(cx), thincolorbar
23	dimitri	1.1	title(['baseline salinity section on day ' int2str(ts*1200/60/60/24)])
24	dimitri	1.2	hold on, plot(lat,-Depth(10,:)/1e3,'r.-'), ylabel('depth (km)')
25	dimitri	1.1	subplot(312)
26	dimitri	1.2	mypcolor(lat,-dpt/1e3,squeeze(s2(10,:,:))'); caxis(cx), thincolorbar
27	dimitri	1.1	title(['downslope salinity section on day ' int2str(ts*1200/60/60/24)])
28	dimitri	1.2	hold on, plot(lat,-Depth(10,:)/1e3,'r.-'), ylabel('depth (km)')
29	dimitri	1.1	subplot(313)
30	dimitri	1.2	mypcolor(lat,-dpt/1e3,squeeze(s3(10,:,:))'); caxis(cx), thincolorbar
31	dimitri	1.1	title(['bbl salinity section on day ' int2str(ts*1200/60/60/24)])
32	dimitri	1.2	hold on, plot(lat,-Depth(10,:)/1e3,'r.-'), ylabel('depth (km)')
33	dimitri	1.3	xlabel('latitude (deg)'), pause(.5)
34	dimitri	1.1	end
35
36			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37			% look at some diagnostics from bbl integration
38	dimitri	1.3	for ts=216:216:7200
39	dimitri	1.1	s1=35+readbin(['run.bbl/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
40	dimitri	1.2	clf, subplot(321), mypcolor(lat,-dpt/1e3,squeeze(s1(10,:,:))'); caxis(cx), thincolorbar
41	dimitri	1.1	title(['salinity section on day ' int2str(ts*1200/60/60/24)])
42			s2=kbot; for i=1:25, for j=1:15, s2(i,j)=s1(i,j,kbot(i,j)); end, end
43	dimitri	1.2	subplot(322), mypcolor(s2'); caxis(cx), thincolorbar
44			title('salinity of bottom grid box')
45	dimitri	1.1	s3=readbin(['run.bbl/BBLsalt.' myint2str(ts,10) '.data'],[25 15]);
46	dimitri	1.2	subplot(323), mypcolor(s3'); caxis(cx), thincolorbar
47			title('salinity of bbl')
48	dimitri	1.1	subplot(324), mypcolor(s3'-s2'); thincolorbar
49			title('bbl salinity minus bottom salinity')
50	dimitri	1.2	tendS=readbin(['run.bbl/BBLtendS.' myint2str(ts,10) '.data'],[25 15]);
51			subplot(325), mypcolor(tendS'); caxis([-1 1]*3e-9), thincolorbar
52			title('bottom salinity tendency due to bbl'), pause(.1)
53			end
54
55			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56			% check that tendency integrates to zero
57			RAC=readbin('run.bbl/RAC.data',[25 15]);
58			hFacC=readbin('run.bbl/hFacC.data',[25 15 25]);
59			for i=1:25
60			for j=1:15
61			RAC(i,j)=RAC(i,j)*hFacC(i,j,kbot(i,j));
62			end
63			end
64	dimitri	1.3	for ts=216:216:7200
65	dimitri	1.2	eta=readbin(['run.bbl/BBLeta.' myint2str(ts,10) '.data'],[25 15]);
66			tendS=readbin(['run.bbl/BBLtendS.' myint2str(ts,10) '.data'],[25 15]);
67			tmp=RAC.*tendS;
68			disp([ts sum(tmp(:))]), pause(1)
69	dimitri	1.1	end