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	gmaze	1.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