krd=2; kgr=1; kwr=0; Tprt=0;

% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/high_res_cube/matlab/use_psiLine.m,v 1.1 2007/03/28 18:07:27 dimitri Exp $
% $Name:  $

if krd > 0,
 Rac='grid_files/';
 drF=squeeze(rdmds([Rac,'DRF']))'; nr=length(drF);

%- to make a plot with grph_CS:
 fprintf(' load xc,yc,xg,yg file ...');
 xcs=rdmds([Rac,'XC']);
 ycs=rdmds([Rac,'YC']);
 xcg=rdmds([Rac,'XG']);
 ycg=rdmds([Rac,'YG']);
 fprintf(' done\n');
 nc=size(xcs,2); nPg=nc*nc*6;
%- add the 2 missing corners:
 xcg=reshape(xcg,[nPg 1]); ycg=reshape(ycg,[nPg 1]);
 xcg(nPg+1)=xcg(1); ycg(nPg+1)=ycg(1+2*nc);
 xcg(nPg+2)=xcg(1+3*nc); ycg(nPg+2)=ycg(1);

 psiLineF=[Rac,'psiLine_N2S_cs',int2str(nc)];
%psiLineF=['psiLine_N2S_cs',int2str(nc)];

 fprintf([' load psiLine from file=',psiLineF,' ...']);
%- load : 'psi_N','psi_C','psi_P','psiUV','psi_T' :
 load(psiLineF);
 fprintf(' done\n');
%- load the grid dx,dy , convert to 10^6 m :
%dxg=rdmds([Rac,'DXG']); dxg=dxg*1.e-6;
%dyg=rdmds([Rac,'DYG']); dyg=dyg*1.e-6;

end

if krd > 1,
%-- load Vertically Integrated transport :
 rac='output_files/';
%nit=72020;
%uu=rdmds([rac,'hUtave'],nit);
%vv=rdmds([rac,'hVtave'],nit);
 namf=[rac,'UintTrsYr.bin'];
   fprintf([' read file: ',namf,' ...']);
 ut=rdda(namf,[nc*6*nc 1],1,'real*8','b');
   fprintf(' done\n');
 namf=[rac,'VintTrsYr.bin'];
   fprintf([' read file: ',namf,' ...']);
 vt=rdda(namf,[nc*6*nc 1],1,'real*8','b');
   fprintf(' done\n');
end

%- compute detph integrated flow :
%delR=[50 70 100 140 190 240 290 340 390 440 490 540 590 640 690]; % = delR in m
%ut=reshape(uu,nPg,nr); vt=reshape(vv,nPg,nr);
%dxg=reshape(dxg,nPg,1); dyg=reshape(dyg,nPg,1);
%var=dyg*delR; ut=var.*ut; ut=sum(ut,2);
%var=dxg*delR; vt=var.*vt; vt=sum(vt,2);

%- compute Barotropic Stream Function :
psiNx=size(psi_C,1);psiNy=size(psi_C,2); nPc2=nPg+2;
psiBa=zeros(nPc2,1);
ufac=rem(psi_T,2) ; vfac=fix(psi_T/2) ;
fprintf(' start computing stream-function ...');
if Tprt, TimeT0=clock; end
for j=2:psiNy,
 for i=1:psi_N(j),
   i1=psi_C(i,j); i0=psi_P(i,j); i2=psiUV(i,j);
   psiBa(i1)=psiBa(i0)+ufac(i,j)*ut(i2)+vfac(i,j)*vt(i2);
 end
end
if Tprt, TimeT1=clock; end
fprintf(' done\n');
if Tprt, fprintf(' Total time: %9.6f \n',etime(TimeT1,TimeT0) ); end

if kwr==1,
 namfil='Psi_Horiz.bin';
 fprintf([' write to file: ',namfil,' ...']);
 fid=fopen(namfil,'w','b'); fwrite(fid,psiBa,'real*8'); fclose(fid);
 fprintf(' done\n');
end

if kgr > 0,
 shift=0; ccB=[0 0]; cbV=3; AxBx=[-180 180 -90 90]; kEnv=1;
 ccB=[-240 240];
%AxBx=[-210 210 30 90];
 grph_CSz(psiBa,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv)
end