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