MITgcm_contrib/gag_subduction_scripts/mldepth.m

function [mld,rho] = mldepth(T,S,depth,epsilon)
%function [mld,rho] = mldepth(T,S,depth,epsilon)
%
% Solve for mixed layer depth on a 1-meter resolution grid.
% 
% Handles input temperature and salinity of any dimension, i.e. 2-D, 3-D,
% 4-D, with time and space in any order.
% 
% Returns mixed layer depth in same dimension as T,S, except without 
% the vertical dimension.
%
% depth = depths on which theta and S are defined.
% epsilon = threshold for density difference, surface to mixld.
%
% Method: Solve for potential density with the surface reference pressure.
%         Interpolate density onto a 1-meter resolution grid.
%         Search for the depth where surface density differs by some
%         threshold. This depth is the mixed layer depth.
%
% G. Gebbie, MIT-WHOI, August 22, 2001. on board the R/V Oceanus.
%
% Vectorized for Matlab, November 2003. GG. MIT-WHOI.
%
% required: seawater toolbox. WARNING: SEAWATER TOOLBOX SYNTAX
%  MAY HAVE CHANGED.

 mldlimit = 500 ;%  a priori maximum limit of mixed layer depth

 S( S == 0) = NaN;
 T( T == 0) = NaN;

% mldlimit is the limit of mixed layer depth here.
 grrid = (2*depth(1)):1:mldlimit;
 
% Set reference pressure to zero. Should not make a difference if mixld < 500.
 pr =0;

%%  The vertical direction is special. Its dimension is specified by "depth". 
 nz = length(depth);

 nn = size(T);

%% Find vertical dimension.
 zindex = find(nn==nz);

 oindex = 1:ndims(T);
 oindex(oindex==zindex)=[];
 nx = prod(nn(oindex));

%% Put the vertical direction at the end. Squeeze the rest.
 temp = permute(T,[oindex zindex]);
 temp = reshape(temp,nx,nz);
% temp (temp==0) = nan;

 salt = permute(S,[1 2 4 3]);
 salt = reshape(salt,nx,nz);
% salt (salt==0) = nan;
  
 pden = sw_pden(salt,temp,depth,pr);

 if nargout ==2
   rho = reshape(pden,[nn(oindex) nz]) ;  
   rho = permute(rho,[1 2 4 3]);
 end
  
 temphi = interp1( depth', pden', grrid);
 differ = cumsum(diff(temphi));

 %% preallocate memory.
 mld = zeros(nx,1);
 
 % how would one vectorize this section?
 for i = 1:nx
    index =find(differ(:,i)> epsilon);
   if( isempty ( index) ==1)
     tmpmld  = NaN;
   else
     tmpmld = grrid( index(1));
   end 
   mld(i) = tmpmld;
 end

 % Make the user happy. Return mixed layer depth in the same form as the
 % input T,S, except vertical dimension is gone.
 
 mld = reshape(mld,[nn(oindex) 1]);
 mld = squeeze(mld);

 mld(isnan(mld)) = 0;
 
 return


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