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

% [psi,mskG,ylat] = calcEulerPsiCube(d,g,flu,rstar,blkFile,[mask]);
%
% Input arguments:
%       d  [structure]    Velocity field (Mass-weighted if rstar=1):
%                         UVELMASS,VVELMASS for rstar=1
%                         UVEL,VVEL         for rstar=0
%       g  [structure]    drF,dxG,dyG,HFacW,HFacS
%       flu [string]      'O' or 'A' for ocean or atmosphere
%       rstar [integer]   1 or 0 if you are using r* coordinates or not
%       blkFile [string]  Broken line file ('isoLat_cs32_59.mat')
%
% Optional inputs:
%       mask [structure]  Optional: Mask field for computation per basin, it assumes that
%                         maskW and maskS are provided in a structure
%
% Output fields:
%       psi     Overturning
%       mskG    Land mask
%       ylat    Latitude coordinate of psi
%
% Description:
%   Calculates overturning streamfunction (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 rstar
%   is on (=1), UVELMASS and VVELMASS must be used. If off (rstar=0),
%   the hfacw*.UVEL and hfacs*.VVEL 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;
nBas=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.
load(blkFile);
ydim = length(bkl_Ylat);
ylat = [-90,bkl_Ylat,90];

% 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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   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