MITgcm_contrib/gmaze_pv/subfct_getvol.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%        SUB-FUNCTIONS          %%
%% USED BY: VOLBET2ISO           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% USE DIRECTLY AT YOUR OWN RISK ! %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Master sub-function:
function varargout = subfct_getvol(CHP,Z,Y,X,LIMITS)

% Limits:
%disp(strcat('Limits: ',num2str(LIMITS)));
O  = LIMITS(1);
MZ = LIMITS(2);
MY = sort( LIMITS(3:4) );
MX = sort( LIMITS(5:6) );

% Compute the volume:
[V Vmat dV] = getvol(Z,Y,X,O,MZ,MY,MX,CHP);

% Outputs:
switch nargout
 case 1
  varargout(1) = {V};
 case 2
  varargout(1) = {V};
  varargout(2) = {Vmat};
 case 3
  varargout(1) = {V};
  varargout(2) = {Vmat};
  varargout(3) = {dV};
end %switch nargout


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the volume limited southward by
% MY(1), northward by iso-O (or MY(2) if iso-O reaches it),
% eastward by MX(1), westward by MX(2) and downward by iso-O
% (or MZ if iso-O reaches it).
%
% Updates: 
% 20060615: add a default meridional gradient (negative)
%           test on x,y,z limits
%
function varargout = getvol(Z,Y,X,O,MZ,MY,MX,CHP)

%% Dim:
nz = length(Z);
ny = length(Y);
nx = length(X);
%disp(num2str([nz ny nx]));

%% Indices:
izmax = min( find( Z>MZ    ) ); if isempty(izmax),izmax=nz;end;
iymin = max( find( Y<MY(1) ) ); if isempty(iymin),iymin=1; end;
iymax = min( find( Y>MY(2) ) ); if isempty(iymax),iymax=ny;end;
ixmin = max( find( X<MX(1) ) ); if isempty(ixmin),ixmin=1; end;
ixmax = min( find( X>MX(2) ) ); if isempty(ixmax),ixmax=nx;end;
%disp(num2str([1 izmax iymin iymax ixmin ixmax]));

%% 1- determine the 3D matrix of volume elements defined by
% the grid:
dV = getdV(Z,Y,X);

%% 2- compute the 3D volume matrix where 1 means dV must be 
% counted and 0 must not:

V = ones(nz,ny,nx); % initialy keep all points

% Exclude northward iso-O limits:
% NB: here the test depends on the meridional gradient of CHP
%     if CHP increase (resp. decreases) with LAT then we must
%     keep lower (resp. higher) values than O limit
% a: determine test type:
N = iymax - iymin + 1; % Number of Y points in the domain
CHPsouth = nanmean(nanmean(squeeze(CHP(1,iymin:iymin+fix(N/2),ixmin:ixmax))));
CHPnorth = nanmean(nanmean(squeeze(CHP(1,iymin+fix(N/2):iymax,ixmin:ixmax))));
SNgrad = (CHPnorth - CHPsouth)./abs(CHPnorth - CHPsouth);
if isnan(SNgrad), SNgrad=-1; end % Assume negative gradient by default
%disp(strcat('Northward gradient sign is:',num2str(SNgrad)));
switch SNgrad
 case  1, testype = 'le'; % Less than  or equal 
 case -1, testype = 'ge'; % Greater than or equal
end %switch
% b: exclude points
for iz=1:izmax
  chpZ = squeeze(CHP(iz,:,:));
  V(iz,:,:)=double(feval(testype,chpZ,O));
end %for iz

% Exclude southward limit:
V(:,1:iymin,:) = zeros(nz,iymin,nx);

% Exclude northward limit:
V(:,iymax:ny,:) = zeros(nz,(ny-iymax)+1,nx);

% Exclude westward limit:
V(:,:,1:ixmin) = zeros(nz,ny,ixmin);

% Exclude eastward limit:
V(:,:,ixmax:nx) = zeros(nz,ny,(nx-ixmax)+1);

% Exclude downward limit:
V(izmax:nz,:,:) = zeros((nz-izmax)+1,ny,nx);


%% 3- Then compute the volume by summing dV elements 
%     for non 0 V points
Vkeep = V.*dV;
Vkeep = sum(sum(sum(Vkeep)));

%% 4- Outputs:
switch nargout
 case 1
  varargout(1) = {Vkeep}; % Volume single value
 case 2
  varargout(1) = {Vkeep};
  varargout(2) = {V};     % Logical V matrix with included/excluded points
 case 3
  varargout(1) = {Vkeep};
  varargout(2) = {V};
  varargout(3) = {dV};    % Volume elements matrix
end %switch nargout


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the 3D dV volume elements.
function dv = getdV(Z,Y,X);

nz = length(Z);
ny = length(Y);
nx = length(X);

%%% Compute the DY:
% Assuming Y is independant of ix:
d  = m_lldist([1 1]*X(1),Y);
dy = [d(1)/2  (d(2:length(d))+d(1:length(d)-1))/2 d(length(d))/2];
dy = meshgrid(dy,X)';

%%% Compute the DX:
clear d
for iy = 1 : ny
   d(:,iy) = m_lldist(X,Y([iy iy]));
end
dx = [d(1,:)/2 ;  ( d(2:size(d,1),:) + d(1:size(d,1)-1,:) )./2 ; d(size(d,1),:)/2];
dx = dx';

%% Compute the horizontal DS surface element:
ds = dx.*dy;

%% Compute the DZ:
d = diff(Z);
dz = [d(1)/2 ( d(2:length(d)) + d(1:length(d)-1) )./2 d(length(d))/2];

%% Then compute the DV volume elements:
dv = ones(nz,ny,nx)*NaN;
for iz=1:nz
  dv(iz,:,:) = dz(iz).*ds;
end

1	gmaze	1.2	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2			%% SUB-FUNCTIONS %%
3			%% USED BY: VOLBET2ISO %%
4			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5			% USE DIRECTLY AT YOUR OWN RISK ! %
6			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7
8
9	gmaze	1.1	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10	gmaze	1.2	% Master sub-function:
11	gmaze	1.1	function varargout = subfct_getvol(CHP,Z,Y,X,LIMITS)
12
13			% Limits:
14			%disp(strcat('Limits: ',num2str(LIMITS)));
15			O = LIMITS(1);
16			MZ = LIMITS(2);
17			MY = sort( LIMITS(3:4) );
18			MX = sort( LIMITS(5:6) );
19
20			% Compute the volume:
21			[V Vmat dV] = getvol(Z,Y,X,O,MZ,MY,MX,CHP);
22
23			% Outputs:
24			switch nargout
25			case 1
26			varargout(1) = {V};
27			case 2
28			varargout(1) = {V};
29			varargout(2) = {Vmat};
30			case 3
31			varargout(1) = {V};
32			varargout(2) = {Vmat};
33			varargout(3) = {dV};
34			end %switch nargout
35
36
37			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38			% This function computes the volume limited southward by
39			% MY(1), northward by iso-O (or MY(2) if iso-O reaches it),
40			% eastward by MX(1), westward by MX(2) and downward by iso-O
41			% (or MZ if iso-O reaches it).
42	gmaze	1.2	%
43			% Updates:
44			% 20060615: add a default meridional gradient (negative)
45			% test on x,y,z limits
46			%
47	gmaze	1.1	function varargout = getvol(Z,Y,X,O,MZ,MY,MX,CHP)
48
49			%% Dim:
50			nz = length(Z);
51			ny = length(Y);
52			nx = length(X);
53			%disp(num2str([nz ny nx]));
54	gmaze	1.2
55	gmaze	1.1	%% Indices:
56	gmaze	1.2	izmax = min( find( Z>MZ ) ); if isempty(izmax),izmax=nz;end;
57			iymin = max( find( Y<MY(1) ) ); if isempty(iymin),iymin=1; end;
58			iymax = min( find( Y>MY(2) ) ); if isempty(iymax),iymax=ny;end;
59			ixmin = max( find( X<MX(1) ) ); if isempty(ixmin),ixmin=1; end;
60			ixmax = min( find( X>MX(2) ) ); if isempty(ixmax),ixmax=nx;end;
61	gmaze	1.1	%disp(num2str([1 izmax iymin iymax ixmin ixmax]));
62
63			%% 1- determine the 3D matrix of volume elements defined by
64			% the grid:
65			dV = getdV(Z,Y,X);
66
67			%% 2- compute the 3D volume matrix where 1 means dV must be
68			% counted and 0 must not:
69
70			V = ones(nz,ny,nx); % initialy keep all points
71
72			% Exclude northward iso-O limits:
73			% NB: here the test depends on the meridional gradient of CHP
74			% if CHP increase (resp. decreases) with LAT then we must
75			% keep lower (resp. higher) values than O limit
76			% a: determine test type:
77			N = iymax - iymin + 1; % Number of Y points in the domain
78			CHPsouth = nanmean(nanmean(squeeze(CHP(1,iymin:iymin+fix(N/2),ixmin:ixmax))));
79			CHPnorth = nanmean(nanmean(squeeze(CHP(1,iymin+fix(N/2):iymax,ixmin:ixmax))));
80			SNgrad = (CHPnorth - CHPsouth)./abs(CHPnorth - CHPsouth);
81	gmaze	1.2	if isnan(SNgrad), SNgrad=-1; end % Assume negative gradient by default
82	gmaze	1.1	%disp(strcat('Northward gradient sign is:',num2str(SNgrad)));
83			switch SNgrad
84			case 1, testype = 'le'; % Less than or equal
85			case -1, testype = 'ge'; % Greater than or equal
86			end %switch
87			% b: exclude points
88			for iz=1:izmax
89			chpZ = squeeze(CHP(iz,:,:));
90			V(iz,:,:)=double(feval(testype,chpZ,O));
91			end %for iz
92
93			% Exclude southward limit:
94			V(:,1:iymin,:) = zeros(nz,iymin,nx);
95
96			% Exclude northward limit:
97			V(:,iymax:ny,:) = zeros(nz,(ny-iymax)+1,nx);
98
99			% Exclude westward limit:
100			V(:,:,1:ixmin) = zeros(nz,ny,ixmin);
101
102			% Exclude eastward limit:
103			V(:,:,ixmax:nx) = zeros(nz,ny,(nx-ixmax)+1);
104
105			% Exclude downward limit:
106			V(izmax:nz,:,:) = zeros((nz-izmax)+1,ny,nx);
107
108
109			%% 3- Then compute the volume by summing dV elements
110			% for non 0 V points
111			Vkeep = V.*dV;
112			Vkeep = sum(sum(sum(Vkeep)));
113
114			%% 4- Outputs:
115			switch nargout
116			case 1
117			varargout(1) = {Vkeep}; % Volume single value
118			case 2
119			varargout(1) = {Vkeep};
120			varargout(2) = {V}; % Logical V matrix with included/excluded points
121			case 3
122			varargout(1) = {Vkeep};
123			varargout(2) = {V};
124			varargout(3) = {dV}; % Volume elements matrix
125			end %switch nargout
126
127
128
129			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
130			% This function computes the 3D dV volume elements.
131			function dv = getdV(Z,Y,X);
132
133			nz = length(Z);
134			ny = length(Y);
135			nx = length(X);
136
137			%%% Compute the DY:
138			% Assuming Y is independant of ix:
139			d = m_lldist([1 1]*X(1),Y);
140			dy = [d(1)/2 (d(2:length(d))+d(1:length(d)-1))/2 d(length(d))/2];
141			dy = meshgrid(dy,X)';
142
143			%%% Compute the DX:
144			clear d
145			for iy = 1 : ny
146			d(:,iy) = m_lldist(X,Y([iy iy]));
147			end
148			dx = [d(1,:)/2 ; ( d(2:size(d,1),:) + d(1:size(d,1)-1,:) )./2 ; d(size(d,1),:)/2];
149			dx = dx';
150
151			%% Compute the horizontal DS surface element:
152			ds = dx.*dy;
153
154			%% Compute the DZ:
155			d = diff(Z);
156			dz = [d(1)/2 ( d(2:length(d)) + d(1:length(d)-1) )./2 d(length(d))/2];
157
158			%% Then compute the DV volume elements:
159			dv = ones(nz,ny,nx)*NaN;
160			for iz=1:nz
161			dv(iz,:,:) = dz(iz).*ds;
162			end
163