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