dfer/matlab_stuff/calcBolusPsiCube.m

function [PsiB,ylat]=calcBolusPsiCube(varargin);

% [PsiB,ylat]=calcBolusPsiCube(d,g,GMform,blkFile);
%
% Compute eddy-induced streamfunction from Gent and McWilliams scheme
%
% Input arguments:  
%   The incoming field data (d) and grid data (g) must be in a structured
%   array format (which is the format that comes from rdmnc):
%       d       [structure]   Kwx, Kwy (Skew flux form) or GM_PsiX, GM_PsiY (advective form)
%       g       [structure]   drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS
%       GMform  [string]      GM form: 'Skew' or 'Advc'
%       blkFile [string]      Broken line file
%       mask    [structure]   Optional: Mask field for computation per basin, it assumes that
%                             maskW and maskS are provided in a structure
% Output arguments:
%       PsiB : bolus streamfunction at interface level (in Sv) 
%       ylat : meridional coordinate of PsiB
%
% Comments:
%   -For Skew-flux form:
%        PsiB computed from Kwx & Kwy divided by 2. 
%        first average Kwx and Kwy at u- and v-points:
%        psiX=(rAc*Kwx)_i / (dXc*dYg) ; psiY=(rAc*Kwy)_j / dYc ;
%        and then "zonally" average along broken lines
%   -For Advective form:
%        PsiB computed from PsiX and PsiY
%        just need to "zonally" average along broken lines
%
%---------------------------------------------------------------------

masking=0;

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Prepare grid stuff                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

switch GMform
case 'Skew'
    nt = size(d.GM_Kwx,4);
case 'Advc'
    nt = size(d.GM_PsiX,4);
end

%--- areas :
ra  = g.rA;
dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]);
dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]);
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]);

rAu=dxc.*dyg;
rAv=dyc.*dxg;

%--- masks :
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]);
mskw=ceil(hw); mskw=min(1,mskw);
msks=ceil(hs); msks=min(1,msks);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Read/Prepare GM fields                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
psiX_all = zeros(6*nc*nc,nr,nt);
psiY_all = zeros(6*nc*nc,nr,nt);

switch GMform
case 'Skew'

kwx_all = 0.5*d.GM_Kwx;
kwy_all = 0.5*d.GM_Kwy;

for it = 1:nt
    kwx = kwx_all(:,:,:,it);
    kwy = kwy_all(:,:,:,it);

    %-- K*ra + add 1 overlap :
    kwx = repmat(ra,[1 1 nr]).*kwx;
    kwy = repmat(ra,[1 1 nr]).*kwy;
    v6X = split_C_cub(kwx,1);
    v6Y = split_C_cub(kwy,1);
    k6x = v6X(:,[2:nc+1],:,:);
    k6y = v6Y([2:nc+1],:,:,:);

    %----------------- 
    clear v6X v6Y
    v6X = (k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:))/2;
    v6Y = (k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:))/2;

    psiX = zeros(6*nc,nc,nr);
    psiY = zeros(6*nc,nc,nr);

    for n = 1:6
        is = 1+nc*(n-1);ie=nc*n;
        psiX([is:ie],[1:nc],[1:nr]) = v6X(:,:,:,n);
        psiY([is:ie],[1:nc],[1:nr]) = v6Y(:,:,:,n);
    end

    psiX = reshape(psiX,6*nc*nc,nr);
    psiY = reshape(psiY,6*nc*nc,nr);

    psiX_all(:,:,it) = mskw .* psiX ./ repmat(rAu,[1,nr]);
    psiY_all(:,:,it) = msks .* psiY ./ repmat(rAv,[1,nr]);

end

case 'Advc'

psiX_all = reshape(d.GM_PsiX(1:6*nc,:,:,1:nt),6*nc*nc,nr,nt);
psiY_all = reshape(d.GM_PsiY(:,1:nc,:,1:nt)  ,6*nc*nc,nr,nt);

otherwise 
    disp(['C est Portnawak: GMform should be Skew or Advc: ',GMform])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     Zonally integrate along broken lines                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load(blkFile);
ydim = length(bkl_Ylat);
ylat = [-90,bkl_Ylat,90];
ufac= rem(bkl_Flg,2);
vfac= fix(bkl_Flg/2);

PsiB= zeros(ydim+2,nr+1,nt);

for it = 1:nt
    for k = 1:nr
        psixt=dyg.*psiX_all(:,k,it).*mskWloc;
        psiyt=dxg.*psiY_all(:,k,it).*mskSloc;
        for jl = 1:ydim
            ie = bkl_Npts(jl);
            PsiB(jl+1,k,it) = sum( ufac(1:ie,jl).*psixt(bkl_IJuv(1:ie,jl)) ...
                      + vfac(1:ie,jl).*psiyt(bkl_IJuv(1:ie,jl)) );
        end
    end
end
PsiB = 1e-6*PsiB;
1	dfer	1.1	function [PsiB,ylat]=calcBolusPsiCube(varargin);
2
3			% [PsiB,ylat]=calcBolusPsiCube(d,g,GMform,blkFile);
4			%
5	dfer	1.2	% Compute eddy-induced streamfunction from Gent and McWilliams scheme
6	dfer	1.1	%
7			% Input arguments:
8			% The incoming field data (d) and grid data (g) must be in a structured
9			% array format (which is the format that comes from rdmnc):
10	dfer	1.3	% d [structure] Kwx, Kwy (Skew flux form) or GM_PsiX, GM_PsiY (advective form)
11			% g [structure] drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS
12	dfer	1.2	% GMform [string] GM form: 'Skew' or 'Advc'
13	dfer	1.3	% blkFile [string] Broken line file
14	dfer	1.2	% mask [structure] Optional: Mask field for computation per basin, it assumes that
15	dfer	1.3	% maskW and maskS are provided in a structure
16	dfer	1.1	% Output arguments:
17			% PsiB : bolus streamfunction at interface level (in Sv)
18			% ylat : meridional coordinate of PsiB
19			%
20			% Comments:
21			% -For Skew-flux form:
22			% PsiB computed from Kwx & Kwy divided by 2.
23			% first average Kwx and Kwy at u- and v-points:
24			% psiX=(rAcKwx)_i / (dXcdYg) ; psiY=(rAc*Kwy)_j / dYc ;
25			% and then "zonally" average along broken lines
26			% -For Advective form:
27			% PsiB computed from PsiX and PsiY
28			% just need to "zonally" average along broken lines
29			%
30			%---------------------------------------------------------------------
31
32			masking=0;
33
34			% Read input parameters.
35			d = varargin{1};
36			g = varargin{2};
37			GMform = varargin{3};
38			blkFile = varargin{4};
39			if length(varargin) == 5
40			mask = varargin{5};
41			masking = 1;
42			end
43
44
45			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46			% Prepare grid stuff %
47			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48
49			nc = size(g.XC,2);
50			nr = length(g.drF);
51
52			switch GMform
53			case 'Skew'
54			nt = size(d.GM_Kwx,4);
55			case 'Advc'
56			nt = size(d.GM_PsiX,4);
57			end
58
59			%--- areas :
60			ra = g.rA;
61			dxc = reshape(g.dxC(1:6nc,1:nc),[6nc*nc,1]);
62			dyc = reshape(g.dyC(1:6nc,1:nc),[6nc*nc,1]);
63			dxg = reshape(g.dxG(1:6nc,1:nc),[6nc*nc,1]);
64			dyg = reshape(g.dyG(1:6nc,1:nc),[6nc*nc,1]);
65
66			rAu=dxc.*dyg;
67			rAv=dyc.*dxg;
68
69			%--- masks :
70			hw = reshape(g.HFacW(1:6nc,1:nc,1:nr),[6nc*nc,nr]);
71			hs = reshape(g.HFacS(1:6nc,1:nc,1:nr),[6nc*nc,nr]);
72			mskw=ceil(hw); mskw=min(1,mskw);
73			msks=ceil(hs); msks=min(1,msks);
74
75			mskWloc = ones(6ncnc,1);
76			mskSloc = ones(6ncnc,1);
77			if masking == 1
78			mskWloc=reshape(mask.maskW(:,:,1),6ncnc,1);
79			mskSloc=reshape(mask.maskS(:,:,1),6ncnc,1);
80			end
81
82			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
83			% Read/Prepare GM fields %
84			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85			psiX_all = zeros(6ncnc,nr,nt);
86			psiY_all = zeros(6ncnc,nr,nt);
87
88			switch GMform
89			case 'Skew'
90
91			kwx_all = 0.5*d.GM_Kwx;
92			kwy_all = 0.5*d.GM_Kwy;
93
94			for it = 1:nt
95			kwx = kwx_all(:,:,:,it);
96			kwy = kwy_all(:,:,:,it);
97
98			%-- K*ra + add 1 overlap :
99			kwx = repmat(ra,[1 1 nr]).*kwx;
100			kwy = repmat(ra,[1 1 nr]).*kwy;
101			v6X = split_C_cub(kwx,1);
102			v6Y = split_C_cub(kwy,1);
103			k6x = v6X(:,[2:nc+1],:,:);
104			k6y = v6Y([2:nc+1],:,:,:);
105
106			%-----------------
107			clear v6X v6Y
108			v6X = (k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:))/2;
109			v6Y = (k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:))/2;
110
111			psiX = zeros(6*nc,nc,nr);
112			psiY = zeros(6*nc,nc,nr);
113
114			for n = 1:6
115			is = 1+nc(n-1);ie=ncn;
116			psiX([is:ie],[1:nc],[1:nr]) = v6X(:,:,:,n);
117			psiY([is:ie],[1:nc],[1:nr]) = v6Y(:,:,:,n);
118			end
119
120			psiX = reshape(psiX,6ncnc,nr);
121			psiY = reshape(psiY,6ncnc,nr);
122
123			psiX_all(:,:,it) = mskw .* psiX ./ repmat(rAu,[1,nr]);
124			psiY_all(:,:,it) = msks .* psiY ./ repmat(rAv,[1,nr]);
125
126			end
127
128			case 'Advc'
129
130			psiX_all = reshape(d.GM_PsiX(1:6nc,:,:,1:nt),6nc*nc,nr,nt);
131			psiY_all = reshape(d.GM_PsiY(:,1:nc,:,1:nt) ,6ncnc,nr,nt);
132
133			otherwise
134			disp(['C est Portnawak: GMform should be Skew or Advc: ',GMform])
135			end
136
137			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
138			% Zonally integrate along broken lines %
139			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140
141			load(blkFile);
142			ydim = length(bkl_Ylat);
143			ylat = [-90,bkl_Ylat,90];
144			ufac= rem(bkl_Flg,2);
145			vfac= fix(bkl_Flg/2);
146
147			PsiB= zeros(ydim+2,nr+1,nt);
148
149			for it = 1:nt
150			for k = 1:nr
151			psixt=dyg.psiX_all(:,k,it).mskWloc;
152			psiyt=dxg.psiY_all(:,k,it).mskSloc;
153			for jl = 1:ydim
154			ie = bkl_Npts(jl);
155			PsiB(jl+1,k,it) = sum( ufac(1:ie,jl).*psixt(bkl_IJuv(1:ie,jl)) ...
156			+ vfac(1:ie,jl).*psiyt(bkl_IJuv(1:ie,jl)) );
157			end
158			end
159			end
160			PsiB = 1e-6*PsiB;