



REstar = [2400 1800 1200 600 0 -600 -1200 -1800];
nR = length(REstar);

Bo = [ 200 1000 4000];
nB = length(Bo);

prod_mon = 11;
prod_day = [9 ];
nday = length(prod_day);
ny =8;

nt = 539;
Umean = nan*ones(nB,nR,ny);
Umean2 = nan*ones(nB,nR,ny);
Ustd = nan*ones(nB,nR,ny);
Umeanx = nan*ones(nB,nR,nt,ny);
UU = cell(nB,nR);
mU = nan*ones(nB,nR,560);
for k=1:nB,

    for i=1:nR,
        disp(i)
        % note U is U(nx, ny, nt)
        [X,Y,time,U] = read_wind_along_file(prod_day,REstar(i),Bo(k));
        disp(sprintf('Bo=%d, REstar=%d,  min time = %.2f,  max time = %.2f ',Bo(k),REstar(i), min(time),max(time)))
        itime = find(time>1002 & time<1021);
        UU{k,i}=U;
        %mU(k,i,:) = squeeze(mean(U(:,1,itime)));
      U2 = U(:,:,itime);
      U2 = permute(U2,[3 1 2]);
      for l=1:8,
          Utt = squeeze(U2(:,:,l));
          Utt = Utt(:);
          Umean2(k,i,l) = mean(Utt);
          Ustd(k,i,l) = std(Utt);
      end
      Umtmp = squeeze(mean(U2,[1 2]));  % this is the x and time mean of U
      Umean(k,i,:) = Umtmp;
      Um = squeeze(mean(U2,2));   % this is just the x-mean of U
                                  %      Umvar = var(Um,1);
                                  %      Ustd(k,i,:) = sqrt(Umvar);
      Umeanx(k,i,:,:) = Um;

      % now Umean is (day=9, REstar, y)
  end;

end;


C = 1.2549;
UC = Umean/C;

%U18 = squeeze(Umean(2,:,:))/C;

a0 = 0.6;
ndy = [1:8]/a0;

save ../mat/wind_1109.mat  C  Umean Ustd Umeanx UU  time itime REstar

figure(10),clf,orient tall
for iy=1:length(Y),
    subplot(4,2,iy)
    plot(REstar,squeeze(Umean(:,:,iy)),'.','markersize',12),grid
    title(sprintf('y=%.0f',Y(iy)))
    if (iy==1),
        legend('Bo=200','Bo=1000','Bo=4000','location','northwest')
    end;
    if (mod(iy,2)==1),
          ylabel('<Uwind>_{tx}')
    end;
    if (i>6),
          xlabel('RE^*')
    end
end;                       
print('-dpdf','Figures/Uwind_REstar_9.pdf')

% Make some plots vs. time at different Bo,REstar
for k=1:nB,
    figure(k),clf,orient tall
    for i=1:nR,
        subplot(4,2,i)
      plot(time(itime),squeeze(Umeanx(k,i,:,:))),grid
      title(sprintf('Bo=%d, REstar=%d, y=[1:8]',Bo(k),REstar(i)))
      tmp = squeeze(Umeanx(k,i,:,:));
      maxU = max(tmp(:));  minU = min(tmp(:));
      if (i<5)
          axis([1001 1022 -0.2 1.1*maxU])
      elseif (i==5)
          axis([1001 1022 -0.05 0.05])
      else
          axis([1001 1022 1.1*minU 0.2])
      end;
      if (i>6),
          xlabel('time')
      end
      if (mod(i,2)==1),
          ylabel('Uwind')
      end;
    end;
    print('-dpdf',sprintf('Figures/Wind_Bo%d_REstar_time.pdf',Bo(k)))
end;


% Make some plots vs. time at different Bo,REstar

for iy=1:5,
    figure(iy),clf,orient tall
    for i=1:nR,
        subplot(4,2,i)
      plot(time(itime),squeeze(Umeanx(:,i,:,iy))),grid
      title(sprintf('y=%.1f REstar=%d',Y(iy),REstar(i)))
%       tmp = squeeze(Umeanx(k,i,:,:));
%       maxU = max(tmp(:));  minU = min(tmp(:));
%       if (i<5)
%           axis([1001 1022 -0.2 1.1*maxU])
%       elseif (i==5)
%           axis([1001 1022 -0.05 0.05])
%       else
%           axis([1001 1022 1.1*minU 0.2])
%       end;
      if (i==1)
          legend('200','1000','4000','location','northwest')
      end
      if (i>6),
          xlabel('time')
      end
      if (mod(i,2)==1),
          ylabel('Uwind')
      end;
    end;
    print('-dpdf',sprintf('Figures/Wind_y%.0f_REstar_time.pdf',Y(iy)))
end;
g



return;

figure(1),clf,orient tall

for i=1:6,
    subplot(3,2,i)
    plot(U9(:,i),U18(:,i),'*'),grid
    ylabel('11-18 run <U>/C')
    xlabel('11-09 run <U>/C')
    if (i==1),
        title(sprintf('RE*=[-18 -12 6 0 6 12 18 24]*100, y/a0 = %.2f',ndy(i)))
    elseif (i==2)
        title(sprintf('Bo=1000, y/a0 = %.2f',ndy(i)))
    else
        title(sprintf('y/a0 = %.2f',ndy(i)))
    end;
    axis([-4 4 -4 4])
end;

%toptitle('REstar = [2400 1800 1200 600 0 -600 -1200 -1800], Bo=1000')
