MITgcm_contrib/gmaze_pv/subfct_getvol.m

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