MITgcm_contrib/gmaze_pv/intbet2outcrops.m

%
% I  = intbet2outcrops(TRACER,LIMITS,LAT,LONG)
%
% This function computes the horizontal surface integral between two 
% outcrops of the TRACER field, given fixed limits eastward, westward 
% and southward.
%
% TRACER = TRACER(LAT,LONG) : surface tracer variable in 2D
% LIMITS = [OUTCROP1 OUTCROP2 MAX_LAT1 MAX_LAT2 MAX_LONG1 MAX_LONG2]
%          : limit's values (MAX_LAT2 is used only if
%            the outcrop's surfaces reach them).
% LAT                             : latitude axis (1D), degrees northward
% LONG                            : longitude axis (1D), degrees east
% I                               : single surface integral value
%
% 06/15/2006
% gmaze@mit.edu
% 


function varargout = intbet2outcrops(TRACER,LIMITS,LAT,LONG)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRE-PROCESS and ERROR CHECK %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pv_checkpath

% Check number of input:
if nargin ~= 4
  help intbet2outcrops.m
  error('intbet2outcrops.m : Wrong number of parameters')
  return
end %if
  
% Check dimensions:
n = size(TRACER);
if length(n)==2
  [ny nx] = size(TRACER);
  if ny~=length(LAT) | nx~=length(LONG)
     help intbet2outcrops.m
     error('intbet2outcrops.m : Axis must have same dimensions than TRACER field');
     return
  end %if
else
  help intbet2outcrops.m
  error('intbet2outcrops.m : TRACER must be a 2D field')
  return
end %if

% Ensure that axis are of dim: (1,N) and well sorted (increasing values):
a=size(LAT);
if a(1) ~= 1,  LAT=LAT'; end 
S = sort(LAT);
if S ~= LAT
  help intbet2outcrops.m
  error('intbet2outcrops.m : LAT must be increasing values')
  return
end %if
a=size(LONG);
if a(1) ~= 1,  LONG=LONG'; end 
S = sort(LONG);
if S ~= LONG
  help intbet2outcrops.m
  error('intbet2outcrops.m : LONG must be increasing values')
  return
end %if

% LIMITS definition:
if length(LIMITS) ~= 6
  help intbet2outcrops.m
  error('intbet2outcrops.m : LIMITS must contains 6 values')
  return
end %if
OUTCROPS = sort( LIMITS(1:2) );
LAT_MAX  = sort( LIMITS(3:4) );
LONG_MAX = sort( LIMITS(5:6) ); 


%%%%%%%%%%%%%%%%%%%%
% COMPUTE INTEGRAL %
%%%%%%%%%%%%%%%%%%%%
% We first determine the element surface matrix and points to integrate:
[I1 I1mat dI1] = subfct_getsurf(TRACER,LAT,LONG,[OUTCROPS(1) LAT_MAX LONG_MAX]);
[I2 I2mat dI2] = subfct_getsurf(TRACER,LAT,LONG,[OUTCROPS(2) LAT_MAX LONG_MAX]);

% Then we determine the outcrop surface limits:
I1mat = abs(I1mat - 1);
Imat  = (I1mat + I2mat)./2;
Imat(find(Imat<1)) = 0;
Imat = logical(Imat);

% And the integral of the TRACER on it:
I = sum(TRACER(Imat).*dI1(Imat));


%%%%%%%%%%%
% OUTPUTS %
%%%%%%%%%%%
switch nargout
 case {0,1}
  varargout(1) = {I};
 case 2
  varargout(1) = {I};
  varargout(2) = {Imat};
 case 3
  varargout(1) = {I};
  varargout(2) = {Imat};
  varargout(3) = {dI1};
end %switch nargout


1	%
2	% I = intbet2outcrops(TRACER,LIMITS,LAT,LONG)
3	%
4	% This function computes the horizontal surface integral between two
5	% outcrops of the TRACER field, given fixed limits eastward, westward
6	% and southward.
7	%
8	% TRACER = TRACER(LAT,LONG) : surface tracer variable in 2D
9	% LIMITS = [OUTCROP1 OUTCROP2 MAX_LAT1 MAX_LAT2 MAX_LONG1 MAX_LONG2]
10	% : limit's values (MAX_LAT2 is used only if
11	% the outcrop's surfaces reach them).
12	% LAT : latitude axis (1D), degrees northward
13	% LONG : longitude axis (1D), degrees east
14	% I : single surface integral value
15	%
16	% 06/15/2006
17	% gmaze@mit.edu
18	%
19
20
21	function varargout = intbet2outcrops(TRACER,LIMITS,LAT,LONG)
22
23
24	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
25	% PRE-PROCESS and ERROR CHECK %
26	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27	pv_checkpath
28
29	% Check number of input:
30	if nargin ~= 4
31	help intbet2outcrops.m
32	error('intbet2outcrops.m : Wrong number of parameters')
33	return
34	end %if
35
36	% Check dimensions:
37	n = size(TRACER);
38	if length(n)==2
39	[ny nx] = size(TRACER);
40	if ny~=length(LAT) \| nx~=length(LONG)
41	help intbet2outcrops.m
42	error('intbet2outcrops.m : Axis must have same dimensions than TRACER field');
43	return
44	end %if
45	else
46	help intbet2outcrops.m
47	error('intbet2outcrops.m : TRACER must be a 2D field')
48	return
49	end %if
50
51	% Ensure that axis are of dim: (1,N) and well sorted (increasing values):
52	a=size(LAT);
53	if a(1) ~= 1, LAT=LAT'; end
54	S = sort(LAT);
55	if S ~= LAT
56	help intbet2outcrops.m
57	error('intbet2outcrops.m : LAT must be increasing values')
58	return
59	end %if
60	a=size(LONG);
61	if a(1) ~= 1, LONG=LONG'; end
62	S = sort(LONG);
63	if S ~= LONG
64	help intbet2outcrops.m
65	error('intbet2outcrops.m : LONG must be increasing values')
66	return
67	end %if
68
69	% LIMITS definition:
70	if length(LIMITS) ~= 6
71	help intbet2outcrops.m
72	error('intbet2outcrops.m : LIMITS must contains 6 values')
73	return
74	end %if
75	OUTCROPS = sort( LIMITS(1:2) );
76	LAT_MAX = sort( LIMITS(3:4) );
77	LONG_MAX = sort( LIMITS(5:6) );
78
79
80
81	%%%%%%%%%%%%%%%%%%%%
82	% COMPUTE INTEGRAL %
83	%%%%%%%%%%%%%%%%%%%%
84	% We first determine the element surface matrix and points to integrate:
85	[I1 I1mat dI1] = subfct_getsurf(TRACER,LAT,LONG,[OUTCROPS(1) LAT_MAX LONG_MAX]);
86	[I2 I2mat dI2] = subfct_getsurf(TRACER,LAT,LONG,[OUTCROPS(2) LAT_MAX LONG_MAX]);
87
88	% Then we determine the outcrop surface limits:
89	I1mat = abs(I1mat - 1);
90	Imat = (I1mat + I2mat)./2;
91	Imat(find(Imat<1)) = 0;
92	Imat = logical(Imat);
93
94	% And the integral of the TRACER on it:
95	I = sum(TRACER(Imat).*dI1(Imat));
96
97
98
99
100	%%%%%%%%%%%
101	% OUTPUTS %
102	%%%%%%%%%%%
103	switch nargout
104	case {0,1}
105	varargout(1) = {I};
106	case 2
107	varargout(1) = {I};
108	varargout(2) = {Imat};
109	case 3
110	varargout(1) = {I};
111	varargout(2) = {Imat};
112	varargout(3) = {dI1};
113	end %switch nargout
114
115
116