close all
clear all

nu=1/1000;
% node coordinates (cell center)
load x.dat

% node coordinates (cell center)
load z.dat

% node coordinates (cell center)
load y_half.dat
y=y_half;

% face coordinates
load yc_half.dat
yc=yc_half;

ni=length(x);
nk=length(z);
nj=length(y_half);

% read u1 & transform u1 to a 3D array
load u3d_re1000.mat % the name of the array is u3d
load v3d_re1000.mat % the name of the array is v3d
load w3d_re1000.mat % the name of the array is w3d
load p3d_re1000.mat % the name of the array is p3d
% x coordinate direction = index i, first index
% y coordinate direction = index j, second index
% z coordinate direction = index k, third index

% compute average u by averaging in x (i) and z (k) directions
umean(1:nj)=0;
for k=1:nk
for j=1:nj
for i=1:ni
   umean(j)=umean(j)+u3d(i,j,k)/ni/nk;
end
end
end

for j=2:nj-1
% dy = cell width
   dy(j)=yc(j)-yc(j-1);
end

figure(1)
% plot u
plot(y,umean)
axis([0 1 0 22])
xlabel('y','fontsize',[20])
ylabel('U+','fontsize',[20])
handle=gca;
set(handle,'fontsize',[20])
print u.ps  -depsc
