function Eddy=eddyDetection_output(lon,lat,sla,upperFolderName,whichYear,whichMonth,whichDay,whichHour,varargin)

if nargin<8

error('Not Enough Inputs. 8 inputs are required.');

if nargin>10

error('Too Many Inputs. No more than 10 inputs please.');

p = inputParser;

p.addRequired('lon',@(x)validateattributes(x,{'numeric'},...

{'nonempty'},'eddyDetection_output','lon',1));

p.addRequired('lat',@(x)validateattributes(x,{'numeric'},...

{'nonempty'},'eddyDetection_output','lat',2));

p.addRequired('sla',@(x)validateattributes(x,{'numeric'},...

{'2d'},'eddyDetection_output','sla',3));

p.addRequired('upperFolderName',@(x)validateattributes(x,{'char'},...

{'nonempty'},'eddyDetection_output','upperFolderName',4));

p.addRequired('whichYear',@(x)validateattributes(x,{'numeric'},...

{'nonempty','positive','integer'},'eddyDetection_output','whichYear',5));

p.addRequired('whichMonth',@(x)validateattributes(x,{'numeric'},...

{'nonempty','positive','integer','>=',1,'<=',12},'eddyDetection_output','whichMonth',6));

p.addRequired('whichDay',@(x)validateattributes(x,{'numeric'},...

{'positive','integer'},'eddyDetection_output','whichDay',7));

p.addRequired('whichHour',@(x)validateattributes(x,{'numeric'},...

{'integer'},'eddyDetection_output','whichHour',8));

defaultAmplitudeThreshold=3; % unit: cm

p.addOptional('amplitudeThreshold',defaultAmplitudeThreshold,...

@(x)validateattributes(x,{'numeric'},...

{'nonzero','integer','positive'},'eddyDetection_output','amplitudeThreshold',9));

defaultRadiusThreshold=45; % unit: Km

p.addOptional('radiusThreshold',defaultRadiusThreshold,...

@(x)validateattributes(x,{'numeric'},...

{'nonzero','positive'},'eddyDetection_output','radiusThreshold',10));

p.parse(lon,lat,sla,upperFolderName,whichYear,whichMonth,whichDay,whichHour,varargin{:});

if size(p.Results.lon,2)==1&size(p.Results.lat,2)==1

[X,Y]=meshgrid(lon,lat);

elseif size(p.Results.lon,1)==1&size(p.Results.lat,1)==1

[X,Y]=meshgrid(lon,lat);

X=lon;Y=lat;

if isequal(size(p.Results.sla),size(X))==0

sla=p.Results.sla;

sla=sla';

sla=p.Results.sla;

if isequal(size(sla),size(X))==0

error('The zonal or meridional grid points of "sla" do not match those of "lon" or "lat".');

m_proj('miller','lon',[min(p.Results.lon(:)) max(p.Results.lon(:))],...

'lat',[min(p.Results.lat(:)) max(p.Results.lat(:))]);

omg=2*pi/24/3600;

f=2*omg*sin(Y/180*pi); % coriolis parameter

f(f<=1.2676e-005&f>=0)=1.2676e-005;

f(f<0&f>-1.2676e-005)=-1.2676e-005;

disy2(1:size(X,1),1:size(X,2))=m_lldist([X(1,1) X(3,1)],[Y(1,1) Y(3,1)])*1000;

disy2(1,:)=m_lldist([X(1,1) X(2,1)],[Y(1,1) Y(2,1)])*1000;

disy2(size(X,1),:)=disy2(1,:);

for cc=1:length(X(:,1))

disx2(cc,1:size(X,2))=m_lldist([X(cc,1) X(cc,3)],[Y(cc,1) Y(cc,3)]).*1000;

i=1; % index for the scenario that 'sla' is a 3-D matrix. Here as it is 2-D matrix, we set 'i' to 1;

clear hhh1

hhh1=sla(:,:,i)*100; % sea level anomaly unit from meter to centimeter

ampThreshold=p.Results.amplitudeThreshold;

radThreshold=p.Results.radiusThreshold;

upperFolderPath=p.Results.upperFolderName;

n=0; % stands for the numbers of detected eddies that meet all the identification thresholds. At first, set it to zero.

for j=100:-1:-100 % for cyclonic eddy detection, scan from positive 100 cm sea level anomaly to negative 100 cm.

levels=[j 1000];cf=figure('visible','off');

[c1,ch1]=contour(X,Y,hhh1,levels); % draw the contours of sea level anomalies that equal to j.

% The elements of each contour are stored in 'c1', including how many location points there are in a certain contour and associated 'lon' and 'lat'.

wz=find(c1(1,:)==j); % Each contour's information starts with c1(1,?)==j. Thus, 'wz' is an index for each contour.

for k=1:length(wz) % length(wz) tells how many contours have been generated.

clear in ky kx bx by in hmeanw EKEmean EKEkg area XX YY U V UV fuqi depth_thermo rr WOeke1 WOgpe1 WOeke2 WOgpe2 tempe jd wd EKEunitmass AGPEunitmass alpha con Dis kmianx kmiany

if c1(2,wz(k))==floor(c1(2,wz(k))) % to insure 'c1(2,wz(k))' is an integer.

if c1(1,wz(k)+1)==c1(1,wz(k)+c1(2,wz(k))) & c1(2,wz(k)+1)==c1(2,wz(k)+c1(2,wz(k)))&c1(2,wz(k))~=2 % to ensure the contour is closed and not only made of two [lon,lat] points.

% to ensure the zonal and meridinal range of this closed contour not less than 0.5 degree. Note that this is an artificial threshold.

if (max(c1(1,wz(k)+1:wz(k)+c1(2,wz(k))-1))-min(c1(1,wz(k)+1:wz(k)+c1(2,wz(k))-1)))>=0.5&(max(c1(2,wz(k)+1:wz(k)+c1(2,wz(k))-1))-min(c1(2,wz(k)+1:wz(k)+c1(2,wz(k))-1)))>=0.5

% Then record the lon/lat coordinates of the contour meet the criterion above.

kmianx=(c1(1,wz(k)+1:wz(k)+c1(2,wz(k))));

kmiany=(c1(2,wz(k)+1:wz(k)+c1(2,wz(k))));

% Calculate the distance between arbitrary pair of points on the contour.

for p=0:length(kmianx)-1;

for q=1:length(kmianx)

Dis(p*length(kmianx)+q)=m_lldist([kmianx(p+1) kmianx(q)],[kmiany(p+1) kmiany(q)]);

end

end

if max(Dis)<=400 % the distance between arbitrary pair of points on the contour should be not longer than 400 km. Note that this is an artificial threshold.

% "b" here means "boundary". Generate the vector contains the location coordinates of closed contour.

bx=[kmianx';kmianx(1)]';

by=[kmiany';kmiany(1)]';

in=inpolygon(X,Y,bx,by); % if value at certain grid point in matrix "in" equals 1, it means that grid point falls within the range of area enclosed by the closed contour.

n1=length(X(in)); % how many grid points have been enclosed by the closed contour.

% sumupx=0;

% sumdownx=0;

% sumupy=0;

% sumdowny=0;

% for k1=wz(k)+1:(wz(k)+c1(2,wz(k))-1)

% sumupx=sumupx+c1(1,k1);

% sumdownx=sumdownx+1;

% sumupy=sumupy+c1(2,k1);

% sumdowny=sumdowny+1;

% end

% xxjd=sumupx/sumdownx;% longitude of eddy centroid

% xxwd=sumupy/sumdowny;% latitude of eddy centroid

% As no eddy exsits in the band of equator due to the requirement of geostrophic balance, the lon/lat of eddy centroid is simply the average of coordinates of the closed countour.

xxjd=mean(kmianx);% longitude of eddy centroid

xxwd=mean(kmiany);% latitude of eddy centroid

if inpolygon(xxjd,xxwd,bx,by)~=0 % to ensure the eddy centroid is enclosed by the closed contour.

XX=X(in); % "XX" contains the longitudes of grid points that enclosed by the closed contour.

YY=Y(in); % "YY" contains the latitudes of grid points that enclosed by the closed contour.

for i1=1:length(XX)

jd(i1)=find(X(1,:)==XX(i1)); % index (of lon/lat grid) for longitudes of grid points that enclosed by the closed contour.

wd(i1)=find(Y(:,1)==YY(i1)); % index (of lon/lat grid) for latitudes of grid points that enclosed by the closed contour.

end

fuqi=hhh1(in); % sea surface anomaly vector enclosed by the closed contour. unit: cm

if min(fuqi)<=j % for cyclonic eddy, its extreme value of sea surface height anomaly should be smaller than the SSHA value of the enclosed contour.

tempe=find(fuqi==min(fuqi)); % find the index for the extreme value of SSHA in this cyclonic eddy.

if length(tempe)>1 % in case 2 or more extremes with same value are found.

tempe=floor(mean(tempe));

end

if j-fuqi(tempe)>=ampThreshold % to ensure the amplitude of this eddy is not less than the input threshold. Note that this is an artificial threshold with default value of 3 cm.

area=disx2(in).*disy2(in)/4;% calculat the area of ocean surface enclosed by the closed contour.

if sqrt(sum(sum(area/1000/1000))/pi)>radThreshold % to ensure the equivalent radius of the eddy is larger than 45 km. This threshold is established due to the spatial resolving capability of altimetry data.

n=n+1; % If all these thresholds are met, finally the procedure successfully found a mesoscale cyclonic eddy.

disp(sprintf('%d%s',n,' eddy(eddies) found.')); % procedure can't wait to tell you the good news.

hhh1(in)=NaN; % wipe out the SSHA information that covered by this new detected eddy for next-round scanning.

Eddy(n).type=-1; % -1 for cyclonic eddy

Eddy(n).center(1)=mean(kmianx); % longitude of eddy centroid

Eddy(n).center(2)=mean(kmiany); % latitude of eddy centroid

Eddy(n).amplitude(1)=XX(tempe); % longitude coordinate of extreme value of eddy's SSHA.

Eddy(n).amplitude(2)=YY(tempe); % latitude coordinate of extreme value of eddy's SSHA.

Eddy(n).amplitude(3)=j-fuqi(tempe); % the amplitude of eddy. unit: cm

Eddy(n).amplitude(4)=ampThreshold; % the amplitude threshold value. unit: cm

Eddy(n).radius(1)=sqrt(sum(sum(area/1000/1000))/pi); % the radius (unit: Km) of a circle which has the same area as the enclosed region by the eddy boundary.

Eddy(n).radius(2)=radThreshold; % the radius threshold value (unit: Km)

Eddy(n).edge(1,1:length(bx))=bx; % the longitude coordinates of points on the eddy boundary.

Eddy(n).edge(2,1:length(by))=by; % the latitude coordinates of points on the eddy boundary.

end

end

end

end

end

end

end

end

end

close(cf);close all; % because the procedure produced an invisible picture of contours.

close all hidden

clear hhh1

hhh1=sla(:,:,i)*100; % sea level anomaly unit from meter to centimeter

clear c1 ch1 wz

for j=-100:100 % this time we scan from -100 cm sea level anomaly to 100 cm.

levels=[j 1000];cf=figure('visible','off');

[c1,ch1]=contour(X,Y,hhh1,levels);

wz=find(c1(1,:)==j);

for k=1:length(wz)

clear in ky kx bx by in hmeanw EKEmean EKEkg area XX YY U V UV fuqi depth_thermo rr WOeke1 WOgpe1 WOeke2 WOgpe2 tempe jd wd EKEunitmass AGPEunitmass alpha con Dis

if c1(2,wz(k))==floor(c1(2,wz(k)))

if c1(1,wz(k)+1)==c1(1,wz(k)+c1(2,wz(k)))&c1(2,wz(k)+1)==c1(2,wz(k)+c1(2,wz(k)))&c1(2,wz(k))~=2

if (max(c1(1,wz(k)+1:wz(k)+c1(2,wz(k))-1))-min(c1(1,wz(k)+1:wz(k)+c1(2,wz(k))-1)))>=0.5&(max(c1(2,wz(k)+1:wz(k)+c1(2,wz(k))-1))-min(c1(2,wz(k)+1:wz(k)+c1(2,wz(k))-1)))>=0.5%2*2

kmianx=(c1(1,wz(k)+1:wz(k)+c1(2,wz(k))-1));

kmiany=(c1(2,wz(k)+1:wz(k)+c1(2,wz(k))-1));

for p=0:length(kmianx)-1;

for q=1:length(kmianx)

Dis(p*length(kmianx)+q)=m_lldist([kmianx(p+1) kmianx(q)],[kmiany(p+1) kmiany(q)]);

end

end

if max(Dis)<=400

bx=[kmianx';kmianx(1)]';

by=[kmiany';kmiany(1)]';

in=inpolygon(X,Y,bx,by);

n1=length(X(in));

% sumupx=0;

% sumdownx=0;

% sumupy=0;

% sumdowny=0;

% for k1=wz(k)+1:(wz(k)+c1(2,wz(k))-1)

% sumupx=sumupx+c1(1,k1);

% sumdownx=sumdownx+1;

% sumupy=sumupy+c1(2,k1);

% sumdowny=sumdowny+1;

% end

% xxjd=sumupx/sumdownx;

% xxwd=sumupy/sumdowny;

xxjd=mean(kmianx);

xxwd=mean(kmiany);

if inpolygon(xxjd,xxwd,bx,by)~=0

XX=X(in);

YY=Y(in);

for i1=1:length(XX)

jd(i1)=find(X(1,:)==XX(i1));

wd(i1)=find(Y(:,1)==YY(i1));

end

fuqi=hhh1(in);

if min(fuqi)>=j

tempe=find(fuqi==max(fuqi)); % for anticyclonic eddy, its extreme value of sea surface height anomaly should be larger than the SSHA value of the enclosed contour.

if length(tempe)>1

tempe=floor(mean(tempe));

end

if fuqi(tempe)-j>=ampThreshold

area=disx2(in).*disy2(in)/4;

if sqrt(sum(sum(area/1000/1000))/pi)>radThreshold

n=n+1;

disp(sprintf('%d%s',n,' eddy(eddies) found.'));

hhh1(in)=NaN;

Eddy(n).type=1; % 1 for eddy polarity: anticyclonic

Eddy(n).center(1)=mean(kmianx);

Eddy(n).center(2)=mean(kmiany);

Eddy(n).amplitude(1)=XX(tempe);

Eddy(n).amplitude(2)=YY(tempe);

Eddy(n).amplitude(3)=fuqi(tempe)-j;

Eddy(n).amplitude(4)=ampThreshold; % the amplitude threshold value. unit: cm

Eddy(n).radius(1)=sqrt(sum(sum(area/1000/1000))/pi);

Eddy(n).radius(2)=radThreshold;

Eddy(n).edge(1,1:length(bx))=bx;

Eddy(n).edge(2,1:length(by))=by;

end

end

end

end

end

end

end

end

end

close(cf);close all;

close all hidden

if isempty(whichDay)==0 & isempty(whichHour)==0

if whichMonth<10 & whichDay<10 & whichHour<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

'0',num2str(whichMonth),'0',num2str(whichDay),'0',num2str(whichHour)]);

elseif whichMonth<10 & whichDay<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

'0',num2str(whichMonth),'0',num2str(whichDay),num2str(whichHour)]);

elseif whichMonth<10 & whichHour<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

'0',num2str(whichMonth),num2str(whichDay),'0',num2str(whichHour)]);

elseif whichDay<10 & whichHour<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

num2str(whichMonth),'0',num2str(whichDay),'0',num2str(whichHour)]);

elseif whichMonth<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

'0',num2str(whichMonth),num2str(whichDay),num2str(whichHour)]);

elseif whichDay<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

num2str(whichMonth),'0',num2str(whichDay),num2str(whichHour)]);

elseif whichHour<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

num2str(whichMonth),num2str(whichDay),'0',num2str(whichHour)]);

end

elseif isempty(whichDay)==0 & isempty(whichHour)==1

if whichMonth<10 & whichDay<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

'0',num2str(whichMonth),'0',num2str(whichDay)]);

elseif whichMonth<10 & whichDay<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

'0',num2str(whichMonth),'0',num2str(whichDay)]);

elseif whichMonth<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

'0',num2str(whichMonth),num2str(whichDay)]);

elseif whichDay<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

num2str(whichMonth),'0',num2str(whichDay)]);

end

elseif isempty(whichDay)==1 & isempty(whichHour)==1

if whichMonth<10

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

'0',num2str(whichMonth)]);

else

outputFolderPath=fullfile(upperFolderPath,[num2str(whichYear),...

num2str(whichMonth)]);

end

mkdir(outputFolderPath);

for i=1:length(Eddy)

nccreate(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'type','Dimensions',{'type',1});

ncwrite(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'type',Eddy(i).type);

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'type', 'standard_name', 'Eddy polarity');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'type', 'long_name', 'Eddy polarity, cyclonic or anticyclonic');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'type', 'description', '-1 for cyclonic eddy while 1 for anticyclonic eddy');

nccreate(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'center','Dimensions',{'center',2});

ncwrite(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'center',Eddy(i).center);

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'center', 'standard_name', 'Eddy centroid');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'center', 'long_name', 'Longitude and latitude coordinate of eddy centroid');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'center', 'description', '[longitude latitude]');

nccreate(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'amplitude','Dimensions',{'amplitude',4});

ncwrite(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'amplitude',[Eddy(i).amplitude]);

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'amplitude', 'standard_name', 'Eddy amplitude');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'amplitude', 'long_name', 'The sea suface height (SSH) difference (positive with unit of centimeters) between extreme value of SSH falls within the range of eddy and the SSH on eddy edge. Company with the lon/lat coordinate of that extreme value of SSH and value of the minimum amplitude threshold.');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'amplitude', 'description', '[longitude latitude amplitude amplitude_threshold]');

nccreate(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'radius','Dimensions',{'radius',2});

ncwrite(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'radius',Eddy(i).radius);

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'radius', 'standard_name', 'Eddy radius');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'radius', 'long_name', 'Radius of an circle that has the same area as the eddy area enclosed by the eddy edge as well as value of the minimum radius threshold. Unit: Km');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'radius', 'description', '[radius radius_threshold]');

nccreate(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'edge','Dimensions',{'edge_row',2,'edge_col',length(Eddy(i).edge)});

ncwrite(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']),'edge',Eddy(i).edge);

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'edge', 'standard_name', 'Eddy edge');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'edge', 'long_name', 'The longitudes and latitudes of eddy edge defined by the largest closed sla contour.');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), 'edge', 'description', 'Longitudes are in the 1st row while the 2nd row is filled with latitudes');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), '/', 'dataset_name', 'Preliminary Mesoscale Eddy Detection Results based on Sea Surface Height');

ncwriteatt(fullfile(outputFolderPath,['eddy',num2str(i),'.nc']), '/', 'dataset_description', 'Outputs are eddy polarity, coordinates of eddy centroid, eddy amplitude and according coordinates, eddy radius and coordinates of eddy edge.');