gmaze_pv/subfct/diagHatisoC.m

% [H,h,[dH,dh]] = diagHatisoC(C,Z,isoC,[dC])
%
% Get depth of C(depth,lat,lon) = isoC
% Z < 0
%
% OUTPUTS:
% H(lat,lon) is the depth determine with the input resolution
% h(lat,lon) is a more accurate depth (determined with interpolation)
% dH(lat,lon) is the thickness of the layer: isoC-dC < C < isoC+dC from H
% dh(lat,lon) is the thickness of the layer: isoC-dC < C < isoC+dC from h
%
% G. Maze, MIT, June 2007
%

function varargout = diagHatisoC(C,Z,isoC,varargin)


% 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPROC
[nz,ny,nx] = size(C);
H = zeros(ny,nx).*NaN;
if nargout >= 2
  h = zeros(ny,nx).*NaN;
  z = [0:-1:Z(end)]; % Vertical axis of the interpolated depth
  if nargin == 4
    dh = zeros(ny,nx).*NaN;
  end
end
if nargin == 4
  dC = varargin{1};
  dH = zeros(ny,nx).*NaN;
end

% 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMPUTING
warning off
 for ix = 1 : nx
   for iy = 1 : ny
     c = squeeze(C(:,iy,ix))';   
     if isnan(c(1)) ~= 1
     if length(find(c>=isoC))>0 & length(find(c<=isoC))>0
        
        % Raw value:
        [cm icm] = min(abs(abs(c)-abs(isoC)));
        H(iy,ix) = Z(icm);
        
        if nargout >= 2
          % Interp guess:
          cc = feval(@interp1,Z,c,z,'linear');
          [cm icm] = min(abs(abs(cc)-abs(isoC)));
          h(iy,ix) = z(icm);
        end % if 2 outputs
        
        if nargin == 4
           [cm icm1] = min(abs(abs(c)-abs(isoC+dC)));
           [cm icm2] = min(abs(abs(c)-abs(isoC-dC)));
           dH(iy,ix) = max(Z([icm1 icm2])) - min(Z([icm1 icm2]));
           
           if nargout >= 2
              [cm icm1] = min(abs(abs(cc)-abs(isoC+dC)));
              [cm icm2] = min(abs(abs(cc)-abs(isoC-dC)));
              dh(iy,ix) = max(z([icm1 icm2])) - min(z([icm1 icm2]));
           end % if 2 outputs
        end % if thickness
        
     end % if found value in the profile
     end % if point n ocean
   end
 end
warning on 
 
% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS
switch nargout
 case 1
  varargout(1) = {H};
 case 2
  varargout(1) = {H};
  varargout(2) = {h};
 case 3
  varargout(1) = {H};
  varargout(2) = {h};
  varargout(3) = {dH};
 case 4
  varargout(1) = {H};
  varargout(2) = {h};
  varargout(3) = {dH};
  varargout(4) = {dh};
end
1	% [H,h,[dH,dh]] = diagHatisoC(C,Z,isoC,[dC])
2	%
3	% Get depth of C(depth,lat,lon) = isoC
4	% Z < 0
5	%
6	% OUTPUTS:
7	% H(lat,lon) is the depth determine with the input resolution
8	% h(lat,lon) is a more accurate depth (determined with interpolation)
9	% dH(lat,lon) is the thickness of the layer: isoC-dC < C < isoC+dC from H
10	% dh(lat,lon) is the thickness of the layer: isoC-dC < C < isoC+dC from h
11	%
12	% G. Maze, MIT, June 2007
13	%
14
15	function varargout = diagHatisoC(C,Z,isoC,varargin)
16
17
18	% 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPROC
19	[nz,ny,nx] = size(C);
20	H = zeros(ny,nx).*NaN;
21	if nargout >= 2
22	h = zeros(ny,nx).*NaN;
23	z = [0:-1:Z(end)]; % Vertical axis of the interpolated depth
24	if nargin == 4
25	dh = zeros(ny,nx).*NaN;
26	end
27	end
28	if nargin == 4
29	dC = varargin{1};
30	dH = zeros(ny,nx).*NaN;
31	end
32
33	% 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMPUTING
34	warning off
35	for ix = 1 : nx
36	for iy = 1 : ny
37	c = squeeze(C(:,iy,ix))';
38	if isnan(c(1)) ~= 1
39	if length(find(c>=isoC))>0 & length(find(c<=isoC))>0
40
41	% Raw value:
42	[cm icm] = min(abs(abs(c)-abs(isoC)));
43	H(iy,ix) = Z(icm);
44
45	if nargout >= 2
46	% Interp guess:
47	cc = feval(@interp1,Z,c,z,'linear');
48	[cm icm] = min(abs(abs(cc)-abs(isoC)));
49	h(iy,ix) = z(icm);
50	end % if 2 outputs
51
52	if nargin == 4
53	[cm icm1] = min(abs(abs(c)-abs(isoC+dC)));
54	[cm icm2] = min(abs(abs(c)-abs(isoC-dC)));
55	dH(iy,ix) = max(Z([icm1 icm2])) - min(Z([icm1 icm2]));
56
57	if nargout >= 2
58	[cm icm1] = min(abs(abs(cc)-abs(isoC+dC)));
59	[cm icm2] = min(abs(abs(cc)-abs(isoC-dC)));
60	dh(iy,ix) = max(z([icm1 icm2])) - min(z([icm1 icm2]));
61	end % if 2 outputs
62	end % if thickness
63
64	end % if found value in the profile
65	end % if point n ocean
66	end
67	end
68	warning on
69
70	% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS
71	switch nargout
72	case 1
73	varargout(1) = {H};
74	case 2
75	varargout(1) = {H};
76	varargout(2) = {h};
77	case 3
78	varargout(1) = {H};
79	varargout(2) = {h};
80	varargout(3) = {dH};
81	case 4
82	varargout(1) = {H};
83	varargout(2) = {h};
84	varargout(3) = {dH};
85	varargout(4) = {dh};
86	end