function [depth, Dd] = kara_mld_1(D,z,DT,D0)
% [depth, Dd] = kara_mld_1(D,z,DT,D0)
% 27 March 2004		David S. Ciochetto
%	
% A Matlab version of the main loop in Kara et al. (2000)
% program (interph.f90) from http://www7320.nrlssc.navy.mil/nmld/nmld.html
%	It has been modified here to take advantage of matlab functions and places more emphasis on the
% calling program to send it just the data to evaluate, not some entire column that may be empty.
%
% The program that calls this m-file should select the section of data to search for the ILD or MLD.
% Thus cut any data above and below the profile that is not desired out of the data that is sent (see
% Note 1 below). Along with the data, the T difference criteria is also sent. Then this program sets
% up an initial estimate of the uniform property depth based on pairs of points in the data profile
% that is sent (D) and a Tolerance = +/-0.1*DT. This is the first point whos difference is greater than 
% DT/10 in absolute value. Then the program trys to estimate the MLD with the method of Kara et al. 2000.
% It compares successive points working downward in the profile intil the lower point exceeds the reference
% value +/- DT. If so it checks if the difference between the current and next point falls under DT/10
% in absolute value. If it does, it defines this next point as the new refernce value and repeats the search
% until the next value exceeds the last reference value +/- DT and abs(diff(D(k:k+1))) > DT/10. The last
% depth that meets these criteria. If a layer is not found, it goes through the profile again and checks
% if it should set the MLD to the entire depth and usually will for uniform profiles.
%
% Basically this is a sophisticated method to estimate a layer depth based on the depth for which the
% absolute difference between the reference value and a local value with DT recommended to be 0.8 deg C.
% For density, the condition is the difference in sigma_T that corresponds to DT = 0.8. For CTD profile
% data to date, this is no different than a simple absolute difference criteria.
%
% INPUTS:
%		D			= the data to be checked for a uniform layer (either T [deg C] or sigma_T [kg/m^3])
%		z			= the depth coordinates associated with the data [m]
%		DT		= the temperature difference criteria to select the depth of the uniform layer
%		D0		= the initial value for the variable in D (for T: T(1)-DT and for dens: dens(1) + Dst)
%
% OUTPUTS:
%		depth = the interpolated estimate of the MLD or ILD if T or sigma_T is sent as D
%		Dd		= the value of the input data D at depth
%
% NOTES:
%		1. D is not the entire water column. Kara et al. (2000) recommend starting the layer search
%				at 10 m depth. Thomson and Vine (2003) suggest using 2.5 m for profile data. This differs
%				from the original Kara et al. (2000) implimentation where they search for the indices of the
%				profile max and min values. The data sent should extend to the max depth of interest.
%
%		2. This is made to work only with vectors of data in columns, not arrays of data.
%	
% REFERENCES:
% 	Kara, A. B., Rochford, P. A. and Hurlburt, H. E. 2000 "An optimal definition
%			for ocean mixed layer depth," J. Geophys. Res., 105, 16803-16821.
%	
%		Thomson, R.E. and I. Fine, 2003 "Estimating Mixed Layer Depth from Oceanic Profile Data,"
%			Journal of Atmospheric and Oceanic Technology, Vol. 20, No. 2, pp. 319-329. 
%	
% USES:
%	
% SEE ALSO: mld.m, MLD_snm.m
%	

%% Things that may be desired to adjust %%
% The tolerance to which the data is checked. Kara et al. (2000) recommend 1/10 of Delta T
TOL = 0.1*DT;
CHK = 'n';			% ['n'|'y'|'yf'] make a check plot / animation, yf is a faster annimation.
verbose = 'n';	% turn debugging comments on/off y/n
depthinit = -5;	% ic for depth searching (later we set the depths to (+) for the analysis)
%%

% Make sure the data is in columns
D=D(:); z=z(:);

% Make sure that the depth is positive (Kara assumes positive depths)
pos = 'y';
idx = z<0;
if sum(idx)>0.1*length(z)								% If over 10% of the values are (+) then invert it.
	pos = 'n';
	z = -z;
	if verbose == 'y'
		disp('Negative numbers were sent, set to positive for analysis')
	end
end

% Check if there is only one point in the profile, if so, that is the desired depth.
%	(after Kara - in practice here, this should be screened for in the calling program.)
if length(D) == 1
	depth = z(end);
	Dd = D(end);
	if verbose == 'n'
		disp(' '); disp('Only one point was found in the profile'); disp(' ');
	end
	return
end

% Kara sends this to the program, in case no value is found.
% Be careful, it differes if T or sigmaT is sent
Dd = D0;

% Initialize vars
imax=0; imin=0; imld=0;

% First Loop of Kara et al. (2000)
% The max and min indicies and values are the 1st and last elements in the arrays of data sent here.
% Find the index of the max and min value for the profile.
imax = find(D == max(D)); imin = find(D == min(D));
% Find the index for the bottom of the first uniform property layer
imld = abs(diff(D)) <= abs(TOL);
if sum(imld)==length(D)-1			% If the values all fall within the range
	imld = find(imld==1);				% Take the bottom value
	imld = imld(end)+1;
else													% If there are some that fall outside of +/- TOL
	imld = find(abs(imld-1));		% Take the first value that falls outside of +/- TOL
	imld = imld(1);
end
if CHK == 'y' % CHECK
	plot(D,z,'.k-'); hold on; grid on; plot(D(imld),z(imld),'ro');
	line([D(1)-TOL D(1)+TOL],[z(1) z(1)],'color','g','LineWidth',3);
	line([D(1)-TOL D(1)-TOL],get(gca,'YLim'),'color','g','LineWidth',2);
	line([D(1)+TOL D(1)+TOL],get(gca,'YLim'),'color','g','LineWidth',2); zoom on; hold off;
end

if imin == 0	% Kara had this, but it does not seem to make sense, imax and imin are used later so it is retained
	imin = imax;
elseif imax == 0
	imax = imin;
end

% Second Loop of Kara - Search for the MLD allow for 2 passes to set it to the entire layer if nothing is found.
% NEED to optimize this (Coded just as Kara, slightly modified)
Lfound = 0; ifound = 1; Lnewsearchval = 1; depth = depthinit;
valref=D(1); zref=z(1);
while ifound<3 & ~Lfound
	if verbose == 'y'
		disp(['ifound = ' num2str(ifound)]) %% A feature to help debug the program.
	end
	% 1st pass - search for MLD taking into account the near surface uniform layer (will move the ref value)
	% From the initial reference value (first point in the data passed here) look for the first value that
	%	falls outside of the bounds ref +/- TOL, then set the reference value to the lower of the two points if
	%	the next value is less than TOL from the point outside, repeat until the value is equal to or below imld.
	for k=2:length(D)
		% Find the min and max values and slope sign for a pair of points k-1,k
		seg=D(k-1:k);
		fldmin=min(seg);
		fldmax=max(seg);
		zseg=z(k-1:k);
		zmin=zseg(find(seg==min(seg)));
		zmax=zseg(find(seg==max(seg)));
		if diff(seg) == 0
			s = 1;
		else
			s = diff(seg)/abs(diff(seg));
		end
		
		% Set valsearch based on the last valref +/- DT
		if fldmin <= valref + s*DT
			valsearch = valref + s*DT;
		elseif fldmax >= valref - s*DT
			valsearch = valref - s*DT;
    elseif Lnewsearchval
			valsearch = valref + s*DT;
			Lnewsearchval = 0;
		end
		
		% CHECK PLOTTING
		if CHK == 'yf' % Can't maintian a zoom level with this one but faster
			plot(D,z,'.k-',D(imld),z(imld),'ro'); grid on; 
			line([D(1)-TOL D(1)+TOL],[zref zref],'color','g','LineWidth',3);
			line([D(1)-TOL D(1)-TOL],get(gca,'YLim'),'color','g','LineWidth',1);
			line([D(1)+TOL D(1)+TOL],get(gca,'YLim'),'color','g','LineWidth',1);
		end
		if CHK == 'y' | CHK == 'yf'	
			hold on;
			plot(fldmin,zmin,'m<',fldmax,zmax,'m>'); plot(valref,zref,'ms');
			line([valsearch valsearch],get(gca,'YLim'),'color','r'); zoom on;
			pause(0.0001)
			if CHK == 'y'
				plot(fldmin,zmin,'<',fldmax,zmax,'>',valref,z(find(D==valref)),'s','color',[1 1 1]);
				line([valsearch valsearch],get(gca,'YLim'),'color',[1 1 1]);
			end
			hold off
		end
		
		% Determine if current pair of data contain valsearch
		%if valsearch < fldmax & valsearch >= fldmin
		% Had to modify this for cases where valsearch was not found before the end of the data vector!!!
		% I hope that it will find the bottom point... yes it does.
		if valsearch < fldmax & valsearch >= fldmin & k ~= length(z)
			% Check if field is uniform below this depth (should be able to let matlab do this better!)
			if abs(diff(D(k:k+1))) <= abs(TOL) & k <= imld
				% Move the reference value to the current (k) value
				valref = D(k); zref=z(k);
				Lnewsearchval = 1;
				if CHK == 'y' % CHECK
					hold on;
					line([D(k)-TOL D(k)+TOL],[zref zref],'color','g');
					line([D(k)-TOL D(k)-TOL],get(gca,'YLim'),'color','g');
					line([D(k)+TOL D(k)+TOL],get(gca,'YLim'),'color','g');
					hold off;
				end
			else		% Determine the mixed layer depth by (linear) interpolation
				depth = (((zmax - zmin)/(fldmax - fldmin))*(valsearch - fldmin)) + zmin;
				Lfound = 1;
				if verbose == 'y'
					disp(['break from loop 1: depth = ' num2str(depth)]); disp(' ');
				end
				break
			end
		end
	end % 1st pass
	
	% 2nd pass - search for fixed changes from initial reference value (NOT TESTED YET)
	if depth == depthinit		% Repeat search using value of D(1)
		valref = D(1);
		if ifound == 2
			depth = z(end);
			if verbose == 'y'
				disp(['Uniform layer not found. Depth set to max depth = ' num2str(depth) ])
			end
		end
	else
		% I do not understand what they intend to do here. I thought they were looking for values that
		% fell outside of the depth range. But they recommend DT=0.8 without discussion of its sign. If
		% I send -0.8, then it does not find a ILD, it selects the whole depth. When I send 0.8, the following
		% if statement fails. It does not hurt the estimate as the first part selects the same depth again.
		if depth > z(imax) & DT > 0		% It seems that they make DT sign specific (but see above)
			valref = D(1);
			Lfound = 0;
			if verbose == 'y'
				disp('Fail on depth > z(imax) & DT > 0');
				disp(['for depth = ' num2str(depth) ', z(imax) = ' num2str(z(imax)) ' : and : DT = ' num2str(DT)])
			end
		elseif depth < z(imin) & DT < 0
			valref = D(1);
			Lfound = 0;
			if verbose == 'y'
				disp('Fail on depth < z(imin) & DT < 0');
				disp(['for depth = ' num2str(depth) ', z(imin) = ' num2str(z(imin)) ' : and : DT = ' num2str(DT)])
			end
		end
	end		% 2nd pass
	ifound = ifound + 1;
	
end

Dd = valsearch;
if verbose == 'y'
	disp(['the value of the data = ' num2str(Dd)])
end

% Convert the depth back to a (+) number if (+) was fed to the routine.
if pos == 'n'
	depth = -depth;
end

if verbose == 'y'
	disp(['The final depth = ' num2str(depth)])
end

if CHK == 'y' | CHK == 'yf'  % CHECK
	hold on; plot(Dd,depth,'bo'); hold off;
end
