function [xplot_mf,yplot_mf,zplot_mf]=diag_mfCOORD(cur_Vflux,cur_RefField,cur_RefValues);

%input:	
%	cur_Vflux 	-> Vdxdz 3D field
%	cur_RefField 	-> Depth, Temperature, or Density 3D field
%	cur_RefValues 	-> 1D vector; values defining the classes center (e.g. density values)
%
%by hypothesis:
%	both cur_Vflux and cur_RefField are provided on V points, with NaN-masked land points
%	the cur_RefValues vector increases
%
%to plot the result: figure; contourf(xplot_mf,-yplot_mf,zplot_mf,[-48:4:48]);

jpi=size(cur_Vflux,1); %number of zonal grid points
jpj=size(cur_Vflux,2); %number of meridional grid points
jpk=size(cur_Vflux,3); %number of vertical grid points

jpk2=length(cur_RefValues);

xplot_mf=[1:jpj]'*ones(1,jpk2);
yplot_mf=ones(jpj,1)*cur_RefValues;
zplot_mf=NaN*ones(jpj,jpk2);

for kcur=1:jpk2
tmp1=cur_Vflux;
if kcur==1
	tmp1(find(cur_RefField>(cur_RefValues(1)+cur_RefValues(2))/2))=NaN; tmp1(find(isnan(tmp1)))=0;
elseif kcur==jpk2
        tmp1(find(cur_RefField<=(cur_RefValues(end)+cur_RefValues(end-1))/2))=NaN; tmp1(find(isnan(tmp1)))=0;
else
        tmp1(find(cur_RefField<=(cur_RefValues(kcur)+cur_RefValues(kcur-1))/2))=NaN;
        tmp1(find(cur_RefField>(cur_RefValues(kcur)+cur_RefValues(kcur+1))/2))=NaN; 
	tmp1(find(isnan(tmp1)))=0;
end
zplot_mf(:,kcur)=squeeze(sum(sum(tmp1,3),1));
end

mask_mf=zplot_mf; mask_mf(find(mask_mf==0))=NaN; mask_mf(find(~isnan(mask_mf)))=1;
%zplot_mf=cumsum(zplot_mf,2)*1e-6;
zplot_mf=flipdim(cumsum(flipdim(zplot_mf,2),2)*1e-6,2);
zplot_mf=zplot_mf.*mask_mf;