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	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