--- MITgcm_contrib/bbl/lookat_output.m	2010/11/18 04:00:04	1.1
+++ MITgcm_contrib/bbl/lookat_output.m	2010/11/22 22:29:55	1.2
@@ -1,47 +1,69 @@
 % 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:21600
- s1=35+readbin(['run/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
- s2=35+readbin(['run.down_slope/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
- s3=35+readbin(['run.bbl/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
+ 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(1:15,25:-1:1,squeeze(s1(10,:,:))');
- caxis([30 30.1]), thincolorbar
+ 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(1:15,25:-1:1,squeeze(s2(10,:,:))');
- caxis([30 30.1]), thincolorbar
+ 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(1:15,25:-1:1,squeeze(s3(10,:,:))');
- caxis([30 30.1]), thincolorbar
+ mypcolor(lat,-dpt/1e3,squeeze(s3(10,:,:))'); caxis(cx), thincolorbar
  title(['bbl salinity section on day ' int2str(ts*1200/60/60/24)])
- pause(.1)
+ hold on, plot(lat,-Depth(10,:)/1e3,'r.-'), ylabel('depth (km)')
+ xlabel('latitude (deg)'), pause(.1)
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % look at some diagnostics from bbl integration
-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
 for ts=216:216:21600
  s1=35+readbin(['run.bbl/SALTanom.' myint2str(ts,10) '.data'],[25 15 25]);
- clf, subplot(321), caxis([30 30.1]), mypcolor(1:15,25:-1:1,squeeze(s1(10,:,:))');
- caxis([30 30.1]), thincolorbar
+ 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([30 30.1]), thincolorbar
- title('salinity at bottom grid')
+ 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([30 30.1]), thincolorbar
- title('bbl salinity')
+ subplot(323), mypcolor(s3'); caxis(cx), thincolorbar
+ title('salinity of bbl')
  subplot(324), mypcolor(s3'-s2'); thincolorbar
  title('bbl salinity minus bottom salinity')
- ts=readbin(['run.bbl/BBLtendS.' myint2str(ts,10) '.data'],[25 15]);
- subplot(325), mypcolor(ts'); thincolorbar
- title('salinity tendency'), pause(.1)
+ 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:21600
+ 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