/[MITgcm]/MITgcm_contrib/gmaze_pv/subfct/diagHatisoC.m
ViewVC logotype

Contents of /MITgcm_contrib/gmaze_pv/subfct/diagHatisoC.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Wed Sep 19 15:24:41 2007 UTC (17 years, 10 months ago) by gmaze
Branch: MAIN
CVS Tags: HEAD
General Update

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

  ViewVC Help
Powered by ViewVC 1.1.22