MITgcm_contrib/gmaze_pv/subfct_getsurf.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Master sub-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) 
%
% Updates: 
% 20060615: add a default meridional gradient (negative)
%           test on x,y limits
%
function varargout = getsurf(Y,X,O,MY,MX,CHP)

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

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