/[MITgcm]/MITgcm_contrib/osse/utils/ensavg.m
ViewVC logotype

Diff of /MITgcm_contrib/osse/utils/ensavg.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.5 by afe, Wed Jun 16 15:29:50 2004 UTC revision 1.6 by afe, Wed Jun 16 21:03:04 2004 UTC
# Line 2  nx=120;ny=31;nz=29;n=nx*ny*nz;nens=30; Line 2  nx=120;ny=31;nz=29;n=nx*ny*nz;nens=30;
2  %iter=68;  %iter=68;
3  z=15;  z=15;
4  disp=0;  disp=0;
5  disprms=0;  disprms=1;
6    jpg=0;
7    vars=['U' 'V' 'W' 'P'];
8    U=find(vars=='U');
9    V=find(vars=='V');
10    W=find(vars=='W');
11    P=find(vars=='P');
12    varnum= size(vars,2);
13    
14    
15  %filename=sprintf('%s','../da/00/assimilate/pickup.in')  %filename=sprintf('%s','../da/00/assimilate/pickup.in')
16  filename=sprintf('%s%03i%s','../da/00/assimilate/pickup.',iter,'.out')  filename=sprintf('%s%03i%s','../da/00/assimilate/pickup.',iter,'.out');
17  fid=fopen(filename,'r','ieee-be');  %filename=sprintf('%s','../da/i90/00/assimilate/pickup.0000016820.001.001.data')
18     fid=fopen(filename,'r','ieee-be');
19  foo=fread(fid,13*n+ny*nx,'float64');  foo=fread(fid,13*n+ny*nx,'float64');
20  %s=0; mn(k+1)=mean(foo(s*n+1:s*n+n));  i=0;tru(:,:,:,U)=reshape(foo(i*n+1:(i+1)*n),[nx ny nz]);
21  i=0;truu=reshape(foo(i*n+1:(i+1)*n),[nx ny nz]);  i=3;tru(:,:,:,V)=reshape(foo(i*n+1:(i+1)*n),[nx ny nz]);
22  i=3;truv=reshape(foo(i*n+1:(i+1)*n),[nx ny nz]);  i=6;tru(:,:,:,W)=reshape(foo(i*n+1:(i+1)*n),[nx ny nz]);
 i=6;truw=reshape(foo(i*n+1:(i+1)*n),[nx ny nz]);  
23    
24    
25    
26  for k=1:nens  for k=1:nens
27  %k=0;  %k=0;
 filename=sprintf('%s%02i%s%03i%s','../da/',k,'/assimilate/pickup.',iter,'.out')  
28  %filename=sprintf('%s%02i%s','../da/',k,'/assimilate/pickup.in')  %filename=sprintf('%s%02i%s','../da/',k,'/assimilate/pickup.in')
29    filename=sprintf('%s%02i%s%03i%s','../da/',k,'/assimilate/pickup.',iter,'.out')
30    %filename=sprintf('%s%02i%s','../da/i90/',k,'/assimilate/pickup.0000016820.001.001.data')
31  fid=fopen(filename,'r','ieee-be');  fid=fopen(filename,'r','ieee-be');
32  foo=fread(fid,13*n+ny*nx,'float64');  foo=fread(fid,13*n+ny*nx,'float64');
33  fclose(fid);  fclose(fid);
34  %filename=sprintf('%s%02i%s%03i%s','../da/',k,'/assimilate/pickup_nh.',iter,'.out')  filename=sprintf('%s%02i%s%03i%s','../da/',k,'/assimilate/pickup_nh.',iter,'.out')
35  %s=0; mn(k+1)=mean(foo(s*n+1:s*n+n));  fid=fopen(filename,'r','ieee-be');
36  %fid=fopen(filename,'r','ieee-be');  foonh=fread(fid,2*n,'float64');
37  %foonh=fread(fid,2*n,'float64');  fclose(fid);
38  %fclose(fid);  i=0;mem(:,U,k)=foo(i*n+1:(i+1)*n);
39  i=0;memu(:,k)=foo(i*n+1:(i+1)*n);  i=3;mem(:,V,k)=foo(i*n+1:(i+1)*n);
40  i=3;memv(:,k)=foo(i*n+1:(i+1)*n);  i=6;mem(:,W,k)=foo(i*n+1:(i+1)*n);
41  i=6;memw(:,k)=foo(i*n+1:(i+1)*n);  i=0;mem(:,P,k)=foonh(i*n+1:(i+1)*n);
 i=0;memp(:,k)=foonh(i*n+1:(i+1)*n);  
42  end  end
43    
44  ensu(:,:,:)=reshape(mean(memu,2),[nx ny nz]);  ens=reshape(mean(mem,3),[nx ny nz varnum]);
 ensv(:,:,:)=reshape(mean(memv,2),[nx ny nz]);  
 ensw(:,:,:)=reshape(mean(memw,2),[nx ny nz]);  
 ensp(:,:,:)=reshape(mean(memp,2),[nx ny nz]);  
45    
46    rms=reshape(std(mem,0,3),[nx ny nz varnum]);
47    
 rmsu=reshape(std(memu,0,2),[nx ny nz]);  
 rmsv=reshape(std(memv,0,2),[nx ny nz]);  
 rmsw=reshape(std(memw,0,2),[nx ny nz]);  
 rmsp=reshape(std(memp,0,2),[nx ny nz]);  
   
 if(disprms)  
48  iterstr=sprintf('%i',iter);  iterstr=sprintf('%i',iter);
49  figure(1);imagesc(rmsu(:,:,z));colorbar;  if(disprms)
50  title(strcat('U ensemble rms, iter ',iterstr));     for i=1:varnum
51  figure(2);imagesc(rmsv(:,:,z));colorbar;        figure(i);imagesc(rms(:,:,z,i));colorbar;
52  title(strcat('V ensemble rms, iter ',iterstr));        title(strcat(vars(i),' ensemble rms, iter ',iterstr));
53  figure(3);imagesc(rmsw(:,:,z));colorbar;     end
 title(strcat('W ensemble rms, iter ',iterstr));  
 figure(4);imagesc(rmsp(:,:,z));colorbar;  
 title(strcat('Phi ensemble rms, iter ',iterstr));  
54  end  end
55    
56  if (jgp)  
57    for i=1:3  if (jpg)
58      for i=1:varnum
59      set(i,'PaperPositionMode', 'manual');      set(i,'PaperPositionMode', 'manual');
60      set(i,'PaperPosition',[0.25 0.25 8 10.5]);      set(i,'PaperPosition',[0.25 0.25 8 10.5]);
61      t=get(get(get(i,'CurrentAxes'),'Title'),'String');      t=get(get(get(i,'CurrentAxes'),'Title'),'String');
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22