dfer/matlab_stuff/calcBolusPsiCube.m

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

% [PsiB,ylat]=calcBolusPsiCube(d,g,GMform,blkFile);
%
% Compute bolus streamfunction from GM 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       [Field data]  Kwx,Kwy
%       g       [Grid data ]  drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS
%       GMform  [string]      GM form 'Skew' or 'Advc'
%       blkFile [file name]   Broken line file  
% 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);

%if masking == 1
%   psiX_all = repmat(reshape(mask.maskW,6*nc*nc,1),[1 nr nt]) .* psiX_all;
%   psiY_all = repmat(reshape(mask.maskS,6*nc*nc,1),[1 nr nt]) .* psiY_all;
%end

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