MITgcm_contrib/gmaze_pv/subfct_getsurf.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Master function:

function varargout = subfct_getsurf(CHP,Y,X,LIMITS)

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

% Compute the surface:
[S Smat dS] = getsurf(Y,X,O,MY,MX,CHP);

% Outputs:
switch nargout
 case 1
  varargout(1) = {S};
 case 2
  varargout(1) = {S};
  varargout(2) = {Smat};
 case 3
  varargout(1) = {S};
  varargout(2) = {Smat};
  varargout(3) = {dS};
end %switch nargout


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the surface 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) 
function varargout = getsurf(Y,X,O,MY,MX,CHP)

%% Dim:
ny = length(Y);
nx = length(X);
%disp(num2str([ny nx]));
%% Indices:
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([iymin iymax ixmin ixmax]));

%% 1- determine the 2D matrix of surface elements defined by
% the grid:
dS = getdS(Y,X);

%% 2- compute the 2D surface matrix where 1 means dS must be 
% counted and 0 must not:

S = ones(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 iso-O limit
% a: determine test type:
N = iymax - iymin + 1; % Number of Y points in the domain
CHPsouth = nanmean(nanmean(squeeze(CHP(iymin:iymin+fix(N/2),ixmin:ixmax))));
CHPnorth = nanmean(nanmean(squeeze(CHP(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
S = double(feval(testype,CHP,O));

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

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

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

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


%% 3- Then compute the surface by summing dS elements 
%     for non nul S points
Skeep = S.*dS;
Skeep = sum(sum(sum(Skeep)));

%% 4- Outputs:
switch nargout
 case 1
  varargout(1) = {Skeep}; % Surface single value
 case 2
  varargout(1) = {Skeep};
  varargout(2) = {S};     % Logical S matrix with included/excluded points
 case 3
  varargout(1) = {Skeep};
  varargout(2) = {S};
  varargout(3) = {dS};    % Surface elements matrix
end %switch nargout


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the 2D dS surface elements.
function ds = getdS(Y,X);

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;

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