/[MITgcm]/MITgcm/verification/advect_xz/input/gendata.m
ViewVC logotype

Diff of /MITgcm/verification/advect_xz/input/gendata.m

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

revision 1.1 by jmc, Tue Jan 3 19:04:51 2006 UTC revision 1.2 by jmc, Mon Dec 5 22:03:52 2011 UTC
# Line 12  ny=1; Line 12  ny=1;
12  nz=20;  nz=20;
13    
14  %- Horizontal & vertical resolution (m):  %- Horizontal & vertical resolution (m):
15  dx=1.e4; dy=dx; dz=1.e4;  dx=1.e4; dy=dx; dz=1.e2;
16  % full size of the domain:  % full size of the domain:
17  Lx=nx*dx;  Lx=nx*dx;
18  H=nz*dz ;  H=nz*dz ;
# Line 46  rhFc=reshape(rhFc,[nx nz]); Line 46  rhFc=reshape(rhFc,[nx nz]);
46    
47  %- horizontal flow field is defined from stream-function psi :  %- horizontal flow field is defined from stream-function psi :
48  psi=zeros(nx+1,nz+1);  psi=zeros(nx+1,nz+1);
49  psi1=Lx*sin(pi*xf'/Lx)*ones(1,nz+1);  psi1=sin(pi*xf'/Lx)*ones(1,nz+1);
50  zu=zeros(nx+1,1); zu(2:nx)=max(zb(2:nx),zb(1:nx-1));  zu=zeros(nx+1,1); zu(2:nx)=max(zb(2:nx),zb(1:nx-1));
51  rHu=zeros(nx+1,1); rHu(2:nx)=-1./zu(2:nx);  rHu=zeros(nx+1,1); rHu(2:nx)=-1./zu(2:nx);
52  psi2=max( ones(nx+1,1)*zf, zu*ones(1,nz+1) );  psi2=max( ones(nx+1,1)*zf, zu*ones(1,nz+1) );
53  psi2=psi2.*(rHu*ones(1,nz+1));  psi2=psi2.*(rHu*ones(1,nz+1));
54  psi2=sin(pi*psi2);  psi2=sin(pi*psi2);
55  psi=psi1.*psi2;  psi=H*psi1.*psi2;
56    
57  uTrans=psi(:,[1:nz])-psi(:,[2:nz+1]);  uTrans=psi(:,[1:nz])-psi(:,[2:nz+1]);
58  u0=uTrans(1:nx,:).*rhFc/dz;  u0=uTrans(1:nx,:).*rhFc/dz;
59    
60    %- add small, vertically uniform, divergent component:
61    ud=psi1(:,1).*rHu;
62    ud=ud*H/160.;
63    u1=ud(1:nx,1)*ones(1,nz);
64    u1(find(hFac==0.))=0.;
65    
66  if w2file,  if w2file,
67   uname='Uvel.bin';   uname='Uvel.bin';
68   fid=fopen(uname,'w',ieee); fwrite(fid,u0,prec); fclose(fid);   fid=fopen(uname,'w',ieee); fwrite(fid,u0,prec); fclose(fid);
69   fprintf([wf,uname,'\n']);   fprintf([wf,uname,'\n']);
70     uname='Udiv.bin';
71     fid=fopen(uname,'w',ieee); fwrite(fid,u0+u1,prec); fclose(fid);
72     fprintf([wf,uname,'\n']);
73  end  end
74    
75  %- Initial Temperature : Gaussian bump :  %- Initial Temperature : Gaussian bump :
# Line 88  end Line 97  end
97  %- make plots to check:  %- make plots to check:
98  figure(1);clf;  figure(1);clf;
99   subplot(211)   subplot(211)
100    %imagesc(xf,zf,psi1'); set(gca,'YDir','normal'); colorbar;
101   imagesc(xf,zf,psi'); set(gca,'YDir','normal'); colorbar;   imagesc(xf,zf,psi'); set(gca,'YDir','normal'); colorbar;
102   hold on ; plot(xc,zb,'b-'); hold off ; grid   hold on ; plot(xc,zb,'b-'); hold off ; grid
103   title('Stream-Function [m^2/s]');   title('Stream-Function [m^2/s]');
104   subplot(212)   subplot(212)
105    %imagesc(xf(1:nx),zc,u1'); set(gca,'YDir','normal'); colorbar;
106   imagesc(xf(1:nx),zc,u0'); set(gca,'YDir','normal'); colorbar;   imagesc(xf(1:nx),zc,u0'); set(gca,'YDir','normal'); colorbar;
107   hold on ; plot(xc,zb,'r-'); hold off ; grid   hold on ; plot(xc,zb,'r-'); hold off ; grid
108   title('U velocity [m/s]');   title('U velocity [m/s]');
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22