gmaze_pv/subfct/getFLUXbudgetV.m

% [D1,D2] = getFLUXbudgetV(z,y,x,Fx,Fy,Fz,box,mask)
%
% Compute the two terms:
% D1 as the volume integral of the flux divergence
% D2 as the surface integral of the normal flux across the volume's boundary
%
% Given a 3D flux vector ie:
%  Fx(z,y,x)
%  Fy(z,y,x)
%  Fz(z,y,x)
%
% Defined on the C-grid at U,V,W locations (bounding the tracer point)
% given by:
% z ( = W detph )
% y ( = V latitude)
% x ( = U longitude)
%
% box is a 0/1 3D matrix defined on the tracer grid
% ie, of dimension: z-1 , y-1 , x-1
% 
% All fluxes are supposed to be scaled by the surface of the cell tile they
% account for.
%
% Each D is decomposed as: 
%  D(1) = Total integral (Vertical+Horizontal)
%  D(2) = Vertical contribution
%  D(3) = Horizontal contribution
%
% Rq:
% The divergence theorem is thus a conservation law which states that 
% the volume total of all sinks and sources, the volume integral of 
% the divergence, is equal to the net flow across the volume's boundary.
%
% gmaze@mit.edu 2007/07/19
%
%

function varargout = getFLUXbudgetV(varargin)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PRELIM

dptw = varargin{1}; ndptw = length(dptw);
latg = varargin{2}; nlatg = length(latg);
long = varargin{3}; nlong = length(long);

ndpt = ndptw - 1;
nlon = nlong - 1;
nlat = nlatg - 1;

Fx = varargin{4};
Fy = varargin{5};
Fz = varargin{6};

if size(Fx,1) ~= ndpt
  disp('Error, Fx(1) wrong dim');
  return
end
if size(Fx,2) ~= nlatg-1
  disp('Error, Fx(2) wrong dim');
  whos Fx 
  return
end
if size(Fx,3) ~= nlong
  disp('Error, Fx(3) wrong dim');
  return
end

pii  = varargin{7};
mask = varargin{8};

% Ensure we're not gonna missed points cause is messy around:
if 1
  Fx(isnan(Fx)) = 0;
  Fy(isnan(Fy)) = 0;
  Fz(isnan(Fz)) = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute the volume integral of flux divergence:
% (gonna be on the tracer grid)
dFdx = ( Fx(:,:,1:nlong-1) - Fx(:,:,2:nlong) );
dFdy = ( Fy(:,1:nlatg-1,:) - Fy(:,2:nlatg,:) );
dFdz = ( Fz(2:ndptw,:,:)   - Fz(1:ndptw-1,:,:) );
%whos dFdx dFdy dFdz

dFdx = dFdx.*mask;
dFdy = dFdy.*mask;
dFdz = dFdz.*mask;

% And sum it over the box:
D1(1) = nansum(nansum(nansum( dFdx.*pii + dFdy.*pii + dFdz.*pii )));
D1(2) = nansum(nansum(nansum( dFdz.*pii )));
D1(3) = nansum(nansum(nansum( dFdy.*pii + dFdx.*pii )));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute the surface integral of the flux:
if nargout > 1
if exist('getVOLbounds')
  method = 3;
else
  method = 2;
end

switch method 
%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
 case 2
  bounds_W = zeros(ndpt,nlat,nlon);
  bounds_E = zeros(ndpt,nlat,nlon);
  bounds_S = zeros(ndpt,nlat,nlon);
  bounds_N = zeros(ndpt,nlat,nlon);
  bounds_T = zeros(ndpt,nlat,nlon);
  bounds_B = zeros(ndpt,nlat,nlon);
  Zflux = 0;
  Mflux = 0;
  Vflux = 0;

  for iz = 1 : ndpt
    for iy = 1 : nlat
      for ix = 1 : nlon
        if pii(iz,iy,ix) == 1
          
          % Is it a western boundary ?
          if ix-1 <= 0 % Reach the domain limit
            bounds_W(iz,iy,ix) = 1;
            Zflux = Zflux + Fx(iz,iy,ix);
          elseif pii(iz,iy,ix-1) == 0 
            bounds_W(iz,iy,ix) = 1;
            Zflux = Zflux + Fx(iz,iy,ix);
          end
          % Is it a eastern boundary ?
          if ix+1 >= nlon % Reach the domain limit
            bounds_E(iz,iy,ix) = 1;
            Zflux = Zflux - Fx(iz,iy,ix+1);
          elseif pii(iz,iy,ix+1) == 0
            bounds_E(iz,iy,ix) = 1;
            Zflux = Zflux - Fx(iz,iy,ix+1);
          end
          
          % Is it a southern boundary ?
          if iy-1 <= 0 % Reach the domain limit
            bounds_S(iz,iy,ix) = 1;
            Mflux = Mflux + Fy(iz,iy,ix);
          elseif pii(iz,iy-1,ix) == 0
            bounds_S(iz,iy,ix) = 1;
            Mflux = Mflux + Fy(iz,iy,ix);
          end
          % Is it a northern boundary ?
          if iy+1 >= nlat % Reach the domain limit
            bounds_N(iz,iy,ix) = 1;
            Mflux = Mflux - Fy(iz,iy+1,ix);
          elseif pii(iz,iy+1,ix) == 0
            bounds_N(iz,iy,ix) = 1;
            Mflux = Mflux - Fy(iz,iy+1,ix);
          end
          
          % Is it a top boundary ?
          if iz-1 <= 0 % Reach the domain limit
            bounds_T(iz,iy,ix) = 1;
            Vflux = Vflux - Fz(iz,iy,ix);
          elseif pii(iz-1,iy,ix) == 0
            bounds_T(iz,iy,ix) = 1;
            Vflux = Vflux - Fz(iz,iy,ix);
          end
          % Is it a bottom boundary ?
          if iz+1 >= ndpt % Reach the domain limit
            bounds_B(iz,iy,ix) = 1;
            Vflux = Vflux + Fz(iz+1,iy,ix);
          elseif pii(iz+1,iy,ix) == 0
            bounds_B(iz,iy,ix) = 1;
            Vflux = Vflux + Fz(iz+1,iy,ix);
          end
          
        end %for iy
      end %for ix       
       
    end
  end  
  
D2(1) = Vflux+Mflux+Zflux;
D2(2) = Vflux;
D2(3) = Mflux+Zflux;


%%%%%%%%%%%%%%%%%%%%
  case 3
  [bounds_N bounds_S bounds_W bounds_E bounds_T bounds_B] = getVOLbounds(pii.*mask);
  Mflux = nansum(nansum(nansum(...
      bounds_S.*squeeze(Fy(:,1:nlat,:)) - bounds_N.*squeeze(Fy(:,2:nlat+1,:)) )));
  Zflux = nansum(nansum(nansum(...
      bounds_W.*squeeze(Fx(:,:,1:nlon)) - bounds_E.*squeeze(Fx(:,:,2:nlon+1)) )));
  Vflux = nansum(nansum(nansum(...
      bounds_B.*squeeze(Fz(2:ndpt+1,:,:))-bounds_T.*squeeze(Fz(1:ndpt,:,:)) )));

  D2(1) = Vflux+Mflux+Zflux;
  D2(2) = Vflux;
  D2(3) = Mflux+Zflux;
  
end %switch method surface flux
end %if we realy need to compute this ?


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS


switch nargout
  case 1
varargout(1) = {D1};
  case 2
varargout(1) = {D1};
varargout(2) = {D2};
  case 3
varargout(1) = {D1};
varargout(2) = {D2};
varargout(3) = {bounds_N};
  case 4
varargout(1) = {D1};
varargout(2) = {D2};
varargout(3) = {bounds_N};
varargout(4) = {bounds_S};
  case 5
varargout(1) = {D1};
varargout(2) = {D2};
varargout(3) = {bounds_N};
varargout(4) = {bounds_S};
varargout(5) = {bounds_W};
  case 6
varargout(1) = {D1};
varargout(2) = {D2};
varargout(3) = {bounds_N};
varargout(4) = {bounds_S};
varargout(5) = {bounds_W};
varargout(6) = {bounds_E};
  case 7
varargout(1) = {D1};
varargout(2) = {D2};
varargout(3) = {bounds_N};
varargout(4) = {bounds_S};
varargout(5) = {bounds_W};
varargout(6) = {bounds_E};
varargout(7) = {bounds_T};
  case 8
varargout(1) = {D1};
varargout(2) = {D2};
varargout(3) = {bounds_N};
varargout(4) = {bounds_S};
varargout(5) = {bounds_W};
varargout(6) = {bounds_E};
varargout(7) = {bounds_T};
varargout(8) = {bounds_B};

end %switch
1	gmaze	1.1	% [D1,D2] = getFLUXbudgetV(z,y,x,Fx,Fy,Fz,box,mask)
2			%
3			% Compute the two terms:
4			% D1 as the volume integral of the flux divergence
5			% D2 as the surface integral of the normal flux across the volume's boundary
6			%
7			% Given a 3D flux vector ie:
8			% Fx(z,y,x)
9			% Fy(z,y,x)
10			% Fz(z,y,x)
11			%
12			% Defined on the C-grid at U,V,W locations (bounding the tracer point)
13			% given by:
14			% z ( = W detph )
15			% y ( = V latitude)
16			% x ( = U longitude)
17			%
18			% box is a 0/1 3D matrix defined on the tracer grid
19			% ie, of dimension: z-1 , y-1 , x-1
20			%
21			% All fluxes are supposed to be scaled by the surface of the cell tile they
22			% account for.
23			%
24			% Each D is decomposed as:
25			% D(1) = Total integral (Vertical+Horizontal)
26			% D(2) = Vertical contribution
27			% D(3) = Horizontal contribution
28			%
29			% Rq:
30			% The divergence theorem is thus a conservation law which states that
31			% the volume total of all sinks and sources, the volume integral of
32			% the divergence, is equal to the net flow across the volume's boundary.
33			%
34			% gmaze@mit.edu 2007/07/19
35			%
36			%
37
38			function varargout = getFLUXbudgetV(varargin)
39
40
41			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PRELIM
42
43			dptw = varargin{1}; ndptw = length(dptw);
44			latg = varargin{2}; nlatg = length(latg);
45			long = varargin{3}; nlong = length(long);
46
47			ndpt = ndptw - 1;
48			nlon = nlong - 1;
49			nlat = nlatg - 1;
50
51			Fx = varargin{4};
52			Fy = varargin{5};
53			Fz = varargin{6};
54
55			if size(Fx,1) ~= ndpt
56			disp('Error, Fx(1) wrong dim');
57			return
58			end
59			if size(Fx,2) ~= nlatg-1
60			disp('Error, Fx(2) wrong dim');
61			whos Fx
62			return
63			end
64			if size(Fx,3) ~= nlong
65			disp('Error, Fx(3) wrong dim');
66			return
67			end
68
69			pii = varargin{7};
70			mask = varargin{8};
71
72			% Ensure we're not gonna missed points cause is messy around:
73			if 1
74			Fx(isnan(Fx)) = 0;
75			Fy(isnan(Fy)) = 0;
76			Fz(isnan(Fz)) = 0;
77			end
78
79			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute the volume integral of flux divergence:
80			% (gonna be on the tracer grid)
81			dFdx = ( Fx(:,:,1:nlong-1) - Fx(:,:,2:nlong) );
82			dFdy = ( Fy(:,1:nlatg-1,:) - Fy(:,2:nlatg,:) );
83			dFdz = ( Fz(2:ndptw,:,:) - Fz(1:ndptw-1,:,:) );
84			%whos dFdx dFdy dFdz
85
86			dFdx = dFdx.*mask;
87			dFdy = dFdy.*mask;
88			dFdz = dFdz.*mask;
89
90			% And sum it over the box:
91			D1(1) = nansum(nansum(nansum( dFdx.pii + dFdy.pii + dFdz.*pii )));
92			D1(2) = nansum(nansum(nansum( dFdz.*pii )));
93			D1(3) = nansum(nansum(nansum( dFdy.pii + dFdx.pii )));
94
95
96			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Compute the surface integral of the flux:
97			if nargout > 1
98			if exist('getVOLbounds')
99			method = 3;
100			else
101			method = 2;
102			end
103
104			switch method
105			%%%%%%%%%%%%%%%%%%%%
106
107			%%%%%%%%%%%%%%%%%%%%
108			case 2
109			bounds_W = zeros(ndpt,nlat,nlon);
110			bounds_E = zeros(ndpt,nlat,nlon);
111			bounds_S = zeros(ndpt,nlat,nlon);
112			bounds_N = zeros(ndpt,nlat,nlon);
113			bounds_T = zeros(ndpt,nlat,nlon);
114			bounds_B = zeros(ndpt,nlat,nlon);
115			Zflux = 0;
116			Mflux = 0;
117			Vflux = 0;
118
119			for iz = 1 : ndpt
120			for iy = 1 : nlat
121			for ix = 1 : nlon
122			if pii(iz,iy,ix) == 1
123
124			% Is it a western boundary ?
125			if ix-1 <= 0 % Reach the domain limit
126			bounds_W(iz,iy,ix) = 1;
127			Zflux = Zflux + Fx(iz,iy,ix);
128			elseif pii(iz,iy,ix-1) == 0
129			bounds_W(iz,iy,ix) = 1;
130			Zflux = Zflux + Fx(iz,iy,ix);
131			end
132			% Is it a eastern boundary ?
133			if ix+1 >= nlon % Reach the domain limit
134			bounds_E(iz,iy,ix) = 1;
135			Zflux = Zflux - Fx(iz,iy,ix+1);
136			elseif pii(iz,iy,ix+1) == 0
137			bounds_E(iz,iy,ix) = 1;
138			Zflux = Zflux - Fx(iz,iy,ix+1);
139			end
140
141			% Is it a southern boundary ?
142			if iy-1 <= 0 % Reach the domain limit
143			bounds_S(iz,iy,ix) = 1;
144			Mflux = Mflux + Fy(iz,iy,ix);
145			elseif pii(iz,iy-1,ix) == 0
146			bounds_S(iz,iy,ix) = 1;
147			Mflux = Mflux + Fy(iz,iy,ix);
148			end
149			% Is it a northern boundary ?
150			if iy+1 >= nlat % Reach the domain limit
151			bounds_N(iz,iy,ix) = 1;
152			Mflux = Mflux - Fy(iz,iy+1,ix);
153			elseif pii(iz,iy+1,ix) == 0
154			bounds_N(iz,iy,ix) = 1;
155			Mflux = Mflux - Fy(iz,iy+1,ix);
156			end
157
158			% Is it a top boundary ?
159			if iz-1 <= 0 % Reach the domain limit
160			bounds_T(iz,iy,ix) = 1;
161			Vflux = Vflux - Fz(iz,iy,ix);
162			elseif pii(iz-1,iy,ix) == 0
163			bounds_T(iz,iy,ix) = 1;
164			Vflux = Vflux - Fz(iz,iy,ix);
165			end
166			% Is it a bottom boundary ?
167			if iz+1 >= ndpt % Reach the domain limit
168			bounds_B(iz,iy,ix) = 1;
169			Vflux = Vflux + Fz(iz+1,iy,ix);
170			elseif pii(iz+1,iy,ix) == 0
171			bounds_B(iz,iy,ix) = 1;
172			Vflux = Vflux + Fz(iz+1,iy,ix);
173			end
174
175			end %for iy
176			end %for ix
177
178			end
179			end
180
181			D2(1) = Vflux+Mflux+Zflux;
182			D2(2) = Vflux;
183			D2(3) = Mflux+Zflux;
184
185
186			%%%%%%%%%%%%%%%%%%%%
187			case 3
188			[bounds_N bounds_S bounds_W bounds_E bounds_T bounds_B] = getVOLbounds(pii.*mask);
189			Mflux = nansum(nansum(nansum(...
190			bounds_S.squeeze(Fy(:,1:nlat,:)) - bounds_N.squeeze(Fy(:,2:nlat+1,:)) )));
191			Zflux = nansum(nansum(nansum(...
192			bounds_W.squeeze(Fx(:,:,1:nlon)) - bounds_E.squeeze(Fx(:,:,2:nlon+1)) )));
193			Vflux = nansum(nansum(nansum(...
194			bounds_B.squeeze(Fz(2:ndpt+1,:,:))-bounds_T.squeeze(Fz(1:ndpt,:,:)) )));
195
196			D2(1) = Vflux+Mflux+Zflux;
197			D2(2) = Vflux;
198			D2(3) = Mflux+Zflux;
199
200			end %switch method surface flux
201			end %if we realy need to compute this ?
202
203
204
205
206			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS
207
208
209			switch nargout
210			case 1
211			varargout(1) = {D1};
212			case 2
213			varargout(1) = {D1};
214			varargout(2) = {D2};
215			case 3
216			varargout(1) = {D1};
217			varargout(2) = {D2};
218			varargout(3) = {bounds_N};
219			case 4
220			varargout(1) = {D1};
221			varargout(2) = {D2};
222			varargout(3) = {bounds_N};
223			varargout(4) = {bounds_S};
224			case 5
225			varargout(1) = {D1};
226			varargout(2) = {D2};
227			varargout(3) = {bounds_N};
228			varargout(4) = {bounds_S};
229			varargout(5) = {bounds_W};
230			case 6
231			varargout(1) = {D1};
232			varargout(2) = {D2};
233			varargout(3) = {bounds_N};
234			varargout(4) = {bounds_S};
235			varargout(5) = {bounds_W};
236			varargout(6) = {bounds_E};
237			case 7
238			varargout(1) = {D1};
239			varargout(2) = {D2};
240			varargout(3) = {bounds_N};
241			varargout(4) = {bounds_S};
242			varargout(5) = {bounds_W};
243			varargout(6) = {bounds_E};
244			varargout(7) = {bounds_T};
245			case 8
246			varargout(1) = {D1};
247			varargout(2) = {D2};
248			varargout(3) = {bounds_N};
249			varargout(4) = {bounds_S};
250			varargout(5) = {bounds_W};
251			varargout(6) = {bounds_E};
252			varargout(7) = {bounds_T};
253			varargout(8) = {bounds_B};
254
255			end %switch