function [psi,mskG,ylat] = calcEulerPsiCube(varargin);

% [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[optional]);
%
% Input arguements:  
%       d  [Field data]  Velocity field (Mass-weighted if rstar=1)
%       g  [Grid data ]  drF,dxG,dyG,HFacW,HFacS
%       flu      (str)   'O' or 'A' for ocean or atmosphere
%       rstar    (int)   0 or 1 if you are using r* coordinates or not
%       blkFile  (str)   Broken line file ('isoLat_cs32_59.mat')
%       The incoming field data (d) and grid data (g) must be in a
%       structured array format
%     Optional parameters:
%       mask   (struct)  mask (structured array including maskW and maskS)
%
% Output fields:
%       psi     Overturning (eg [61,6,nt])
%       mskG    Land mask (eg [60,5])
%       ylat    Latitude coordinate of psi (eg [61,1])
%
% Description:
%   Caculates overturning stream function (psi).  For the atmosphere, data
%   is must be in p-coordinates and the output is the mass transport psi
%   [10^9 kg/s].  For the ocean, data should be in z-coordinates and the
%   output is the volume transport psi [10^6 m^3/s = Sv].  If the rstar
%   parameters is on, hu and hv are used, if off, the hfacw*.u and hfacs*.v
%   are used (the multiplication being done inside the function).
%
%   'psi' is tabulated on broken lines at the interface between cells in
%   the vertical.  'mskG' is for the area between broken lines and between
%   the cell interfaces in the vertical.
%

% Defaults that can be overriden.
grav = 9.81;
masking=0;

% Read input parameters.
d = varargin{1};
g = varargin{2};
flu = varargin{3};
rstar = varargin{4};
blkFile = varargin{5};
if length(varargin) == 6
   mask = varargin{6};
   masking = 1;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Prepare / reform  incoming data                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nc = size(g.XC,2);
nr = length(g.drF);

delM = g.drF;
dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]);
dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]);
if rstar
    nt = size(d.UVELMASS,4);
    hu = reshape(d.UVELMASS(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
    hv = reshape(d.VVELMASS(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
else
    nt = size(d.UVEL,4);
    hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
    hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
    hu = reshape(d.UVEL(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
    hv = reshape(d.VVEL(1:6*nc,1:nc,1:nr,1:nt),[6*nc*nc,nr,nt]);
    for it = 1:nt
        hu(:,:,it) = hw.*hu(:,:,it);
        hv(:,:,it) = hs.*hv(:,:,it);
    end
end

mskWloc = ones(6*nc*nc,1);
mskSloc = ones(6*nc*nc,1);

if masking == 1
   mskWloc=reshape(mask.maskW(:,:,1),6*nc*nc,1);
   mskSloc=reshape(mask.maskS(:,:,1),6*nc*nc,1);
   %hu = repmat(reshape(mask.maskW,6*nc*nc,1),[1 nr nt]) .* hu;
   %hv = repmat(reshape(mask.maskS,6*nc*nc,1),[1 nr nt]) .* hv;
end

% Load broken information.
% I looked at calc_psiH_CS.m, and did not find it very clear.
% May be you can try to see what is in 
% MITgcm/utils/matlab/cs_grid/bk_line/use_psiLine.m
% it s shorter, and slightly better.
load(blkFile);
ydim = length(bkl_Ylat);
ylat = [-90,bkl_Ylat,90];

% kMsep=1;
% if (nargin < 6), kfac=0;
% else kfac=1; end;
nBas=0;

% Prepare arrays.
psi = zeros(ydim+2,nr+1,1+nBas,nt); 
mskZ = zeros(ydim+2,nr+1,1+nBas);  % Mask of psi
mskV = zeros(ydim+2,nr,1+nBas);    % Mask of the merid. transport
mskG = zeros(ydim+1,nr,1+nBas);    % Mask of the ground

% The variable "bkl_Flg" is -1/1 if edge (on a given broken) has a u point
% and -2/2 if it has a v point.  Positive/negative values contribute
% positively/negatively to northward heat transport (this depends on the
% oreientation of the cell).  A zero value indicates an end of edges that
% contribute to a broken line.  The u and v information is parced into two
% seperate fields, ufac and vfac (-2/2 are reduced to -1/1 for vfac).
ufac = zeros([size(bkl_Flg),1+nBas]);
vfac = zeros([size(bkl_Flg),1+nBas]);
ufac(:,:,1) = rem(bkl_Flg,2);
vfac(:,:,1) = fix(bkl_Flg/2);
% for jl=1:ydim,
%     ie=bkl_Npts(jl);
%     for b=1:nBas, 
%         ufac(1:ie,jl,1+b)=mskBw(bkl_IJuv(1:ie,jl),b).*ufac(1:ie,jl,1);
%         vfac(1:ie,jl,1+b)=mskBs(bkl_IJuv(1:ie,jl),b).*vfac(1:ie,jl,1);
%     end;
% end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Compute mass/volume stream function                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Compute volume transport through broken lines a hence psi.  ut/vt is the
%  velocity times the edge length it is passing through.  The sum of this
%  quantity along a broken line (vz) times the cell height is the volume
%  transport through broken line at one layer (delM(k)*vz).  psi is then
%  the value of the volume transport through the level above subtracted
%  from the value of psi above.
for it = 1:nt
    for k = nr:-1:1
        ut = dyg.*hu(:,k,it).*mskWloc;
        vt = dxg.*hv(:,k,it).*mskSloc; 
        for jl = 1:ydim
            ie = bkl_Npts(jl);
            for b = 1:1+nBas
                vz = sum(   ufac(1:ie,jl,b).*ut(bkl_IJuv(1:ie,jl)) ...
                          + vfac(1:ie,jl,b).*vt(bkl_IJuv(1:ie,jl)) );
                psi(jl+1,k,b,it) = psi(jl+1,k+1,b,it) - delM(k)*vz;
            end
        end
    end
end

psi = squeeze(psi);

%% For Ocean, result in Sv (10^6 m3/s)
%% For Atmos, results in 10^9 kg/s
if isequal(flu,'O'), psi = 1e-6*squeeze(psi); end
if isequal(flu,'A'), psi =-1e-9/grav*squeeze(psi); end



% % Compute the mask :
% if kfac == 1, 
%     ufac=abs(ufac) ; vfac=abs(vfac);
%     for jl=1:ydim, 
%         ie=bkl_Npts(jl); 
%         hw=zeros(ie,nr); hs=zeros(ie,nr);
%         hw=hw(bkl_IJuv(1:ie,jl),:); % Would need correction!
%         hs=hs(bkl_IJuv(1:ie,jl),:);
%         for b=1:1+nBas,
%             for k=1:nr,
%                 %   for ii=1:bkl_Npts(jl);
%                 %    ij=bkl_IJuv(ii,jl);
%                 %    mskV(jl+1,k,b)=mskV(jl+1,k,b)+ufac(ii,jl,b)*hw(ij,k)+vfac(ii,jl,b)*hs(ij,k);
%                 %   end ;
%                  tmpv=ufac(1:ie,jl,b).*hw(:,k)+vfac(1:ie,jl,b).*hs(:,k);
%                  mskV(jl+1,k,b)=mskV(jl+1,k,b)+max(tmpv);
%             end
%         end
%     end
%     mskV=ceil(mskV); mskV=min(1,mskV);
%     %- build the real mask (=mskG, ground) used to draw the continent with "surf":
%     %   position=centered , dim= ydim+1 x nr
%     mskG=mskV(1:ydim+1,:,:)+mskV(2:ydim+2,:,:); mskG=min(1,mskG);
% 
%     if kMsep & nBas > 0,
%         mskW=1+min(1,ceil(hw));
%         mskS=1+min(1,ceil(hs));
%         for b=1:nBas,
%             bs=b; b1=1+bs; b2=2+rem(bs,nBas);
%             if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end
%             for j=1:ydim+1,
%                 for i=1:np_Sep(bs,j),
%                     ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i));
%                     if typ == 1, 
%                         mskG(j,:,b1)=mskG(j,:,b1).*mskW(ij,:);
%                         mskG(j,:,b2)=mskG(j,:,b2).*mskW(ij,:);
%                     elseif typ == 2,
%                         mskG(j,:,b1)=mskG(j,:,b1).*mskS(ij,:);
%                         mskG(j,:,b2)=mskG(j,:,b2).*mskS(ij,:);
%                     end
%                 end
%             end
%         end
%         mskG=min(2,mskG);
%     end
% 
%     %- to keep psi=0 on top & bottom
%     mskZ(:,[2:nr+1],:)=mskV; 
%     mskZ(:,[1:nr],:)=mskZ(:,[1:nr],:)+mskV;
%     %- to keep psi=0 on lateral boundaries :
%     mskZ([1:ydim],:,:)=mskZ([1:ydim],:,:)+mskZ([2:ydim+1],:,:);
%     mskZ([2:ydim+1],:,:)=mskZ([2:ydim+1],:,:)+mskZ([3:ydim+2],:,:);
%     mskZ=ceil(mskZ); mskZ=min(1,mskZ);
%     if kMsep & nBas > 0,
%         mskM=zeros(ydim+2,nr,1+nBas); mskM(2:ydim+2,:,:)=min(2-mskG,1); 
%         mskM(1:ydim+1,:,:)=mskM(1:ydim+1,:,:)+mskM(2:ydim+2,:,:);
%         mskZ(:,1:nr,:)=min(mskZ(:,1:nr,:),mskM); 
%     end
%     %- apply the mask (and remove dim = 1) :
%     if nt == 1,
%         psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ);
%         psi( find(mskZ==0) )=NaN ;
%     else
%         for nt=1:nt,
%             psi1=psi(:,:,:,nt); psi1( find(mskZ==0) )=NaN ; psi(:,:,:,nt)=psi1;
%         end
%         if nBas < 1, psi=squeeze(psi); mskV=squeeze(mskV); mskZ=squeeze(mskZ); end
%     end
% else
%     if nBas < 1 | nt == 1, psi=squeeze(psi); end
% end