| 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 |
Parent Directory
| 
Revision Log
| 
Revision Graph