gael/matlab/diag_mfCOORD.m

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;


1	gforget	1.1	function [xplot_mf,yplot_mf,zplot_mf]=diag_mfCOORD(cur_Vflux,cur_RefField,cur_RefValues);
2
3			%input:
4			% cur_Vflux -> Vdxdz 3D field
5			% cur_RefField -> Depth, Temperature, or Density 3D field
6			% cur_RefValues -> 1D vector; values defining the classes center (e.g. density values)
7			%
8			%by hypothesis:
9			% both cur_Vflux and cur_RefField are provided on V points, with NaN-masked land points
10			% the cur_RefValues vector increases
11			%
12			%to plot the result: figure; contourf(xplot_mf,-yplot_mf,zplot_mf,[-48:4:48]);
13
14			jpi=size(cur_Vflux,1); %number of zonal grid points
15			jpj=size(cur_Vflux,2); %number of meridional grid points
16			jpk=size(cur_Vflux,3); %number of vertical grid points
17
18			jpk2=length(cur_RefValues);
19
20			xplot_mf=[1:jpj]'*ones(1,jpk2);
21			yplot_mf=ones(jpj,1)*cur_RefValues;
22			zplot_mf=NaN*ones(jpj,jpk2);
23
24			for kcur=1:jpk2
25			tmp1=cur_Vflux;
26			if kcur==1
27			tmp1(find(cur_RefField>(cur_RefValues(1)+cur_RefValues(2))/2))=NaN; tmp1(find(isnan(tmp1)))=0;
28			elseif kcur==jpk2
29			tmp1(find(cur_RefField<=(cur_RefValues(end)+cur_RefValues(end-1))/2))=NaN; tmp1(find(isnan(tmp1)))=0;
30			else
31			tmp1(find(cur_RefField<=(cur_RefValues(kcur)+cur_RefValues(kcur-1))/2))=NaN;
32			tmp1(find(cur_RefField>(cur_RefValues(kcur)+cur_RefValues(kcur+1))/2))=NaN;
33			tmp1(find(isnan(tmp1)))=0;
34			end
35			zplot_mf(:,kcur)=squeeze(sum(sum(tmp1,3),1));
36			end
37
38			mask_mf=zplot_mf; mask_mf(find(mask_mf==0))=NaN; mask_mf(find(~isnan(mask_mf)))=1;
39			%zplot_mf=cumsum(zplot_mf,2)*1e-6;
40			zplot_mf=flipdim(cumsum(flipdim(zplot_mf,2),2)*1e-6,2);
41			zplot_mf=zplot_mf.*mask_mf;
42
43