1 |
function [mld,rho] = mldepth(T,S,depth,epsilon) |
2 |
%function [mld,rho] = mldepth(T,S,depth,epsilon) |
3 |
% |
4 |
% Solve for mixed layer depth on a 1-meter resolution grid. |
5 |
% |
6 |
% Handles input temperature and salinity of any dimension, i.e. 2-D, 3-D, |
7 |
% 4-D, with time and space in any order. |
8 |
% |
9 |
% Returns mixed layer depth in same dimension as T,S, except without |
10 |
% the vertical dimension. |
11 |
% |
12 |
% depth = depths on which theta and S are defined. |
13 |
% epsilon = threshold for density difference, surface to mixld. |
14 |
% |
15 |
% Method: Solve for potential density with the surface reference pressure. |
16 |
% Interpolate density onto a 1-meter resolution grid. |
17 |
% Search for the depth where surface density differs by some |
18 |
% threshold. This depth is the mixed layer depth. |
19 |
% |
20 |
% G. Gebbie, MIT-WHOI, August 22, 2001. on board the R/V Oceanus. |
21 |
% |
22 |
% Vectorized for Matlab, November 2003. GG. MIT-WHOI. |
23 |
% |
24 |
% required: seawater toolbox. WARNING: SEAWATER TOOLBOX SYNTAX |
25 |
% MAY HAVE CHANGED. |
26 |
|
27 |
mldlimit = 500 ;% a priori maximum limit of mixed layer depth |
28 |
|
29 |
S( S == 0) = NaN; |
30 |
T( T == 0) = NaN; |
31 |
|
32 |
% mldlimit is the limit of mixed layer depth here. |
33 |
grrid = (2*depth(1)):1:mldlimit; |
34 |
|
35 |
% Set reference pressure to zero. Should not make a difference if mixld < 500. |
36 |
pr =0; |
37 |
|
38 |
%% The vertical direction is special. Its dimension is specified by "depth". |
39 |
nz = length(depth); |
40 |
|
41 |
nn = size(T); |
42 |
|
43 |
%% Find vertical dimension. |
44 |
zindex = find(nn==nz); |
45 |
|
46 |
oindex = 1:ndims(T); |
47 |
oindex(oindex==zindex)=[]; |
48 |
nx = prod(nn(oindex)); |
49 |
|
50 |
%% Put the vertical direction at the end. Squeeze the rest. |
51 |
temp = permute(T,[oindex zindex]); |
52 |
temp = reshape(temp,nx,nz); |
53 |
% temp (temp==0) = nan; |
54 |
|
55 |
salt = permute(S,[1 2 4 3]); |
56 |
salt = reshape(salt,nx,nz); |
57 |
% salt (salt==0) = nan; |
58 |
|
59 |
pden = sw_pden(salt,temp,depth,pr); |
60 |
|
61 |
if nargout ==2 |
62 |
rho = reshape(pden,[nn(oindex) nz]) ; |
63 |
rho = permute(rho,[1 2 4 3]); |
64 |
end |
65 |
|
66 |
temphi = interp1( depth', pden', grrid); |
67 |
differ = cumsum(diff(temphi)); |
68 |
|
69 |
%% preallocate memory. |
70 |
mld = zeros(nx,1); |
71 |
|
72 |
% how would one vectorize this section? |
73 |
for i = 1:nx |
74 |
index =find(differ(:,i)> epsilon); |
75 |
if( isempty ( index) ==1) |
76 |
tmpmld = NaN; |
77 |
else |
78 |
tmpmld = grrid( index(1)); |
79 |
end |
80 |
mld(i) = tmpmld; |
81 |
end |
82 |
|
83 |
% Make the user happy. Return mixed layer depth in the same form as the |
84 |
% input T,S, except vertical dimension is gone. |
85 |
|
86 |
mld = reshape(mld,[nn(oindex) 1]); |
87 |
mld = squeeze(mld); |
88 |
|
89 |
mld(isnan(mld)) = 0; |
90 |
|
91 |
return |
92 |
|
93 |
|
94 |
|
95 |
|
96 |
|
97 |
|