clear path global Nx Ny Nz global lat long dz dm mdep global delt_su su_its t_su delt global descriptor this_path global f deltaf Q beta r_expt r_heat H global time rots it global g Cp rho_bar alpha global u v t w global iterations param_file_name = ... input(' Please enter the name of the m-file with the parameters for this run : ','s') ; feval(param_file_name) ; % iterations itstart = input(' Please enter start iteration : ','s') itend = input(' Please enter end iteration : ','s') sizeit=size(iterations); for i=1:sizeit(1) iter(i)=eval(iterations(i,1:10)); end nitstart=find(iter==eval(itstart)) nitend=find(iter==eval(itend)) path = this_path cmdstr=['cd ' path ]; eval(cmdstr); path=pwd sumtheta=zeros(Nx,Ny,Nz); sumu=zeros(Nx,Ny,Nz); sumv=zeros(Nx,Ny,Nz); sumw=zeros(Nx,Ny,Nz); counter=0; for i=nitstart:nitend i,0 tfilename=(['T.' iterations((i),1:10) ]) ; t=rdmds(tfilename,'b'); ufilename=(['U.' iterations((i),1:10) ]) ; u=rdmds(ufilename,'b'); u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:)); vfilename=(['V.' iterations((i),1:10) ]) ; v=rdmds(vfilename,'b'); v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:)); wfilename=(['W.' iterations((i),1:10) ]) ; w=rdmds(wfilename,'b'); sumtheta=sumtheta+t; sumu=sumu+u; sumv=sumv+v; sumw=sumw+w; counter=counter+1; end meantheta=sumtheta/counter; meanu=sumu/counter; meanv=sumv/counter; meanw=sumw/counter; sumut=zeros(Nx,Ny,Nz); sumvt=zeros(Nx,Ny,Nz); sumwt=zeros(Nx,Ny,Nz); meanu=smooth3(meanu); meanv=smooth3(meanv); meanw=smooth3(meanw); meantheta=smooth3(meantheta); for i=nitstart:nitend i tfilename=(['T.' iterations((i),1:10) ]) ; t=rdmds(tfilename,'b'); ufilename=(['U.' iterations((i),1:10) ]) ; u=rdmds(ufilename,'b'); u(1:Nx-1,:,:)=0.5*(u(1:Nx-1,:,:)+u(2:Nx,:,:)); vfilename=(['V.' iterations((i),1:10) ]) ; v=rdmds(vfilename,'b'); v(:,1:Ny-1,:)=0.5*(v(:,1:Ny-1,:)+v(:,2:Ny,:)); wfilename=(['W.' iterations((i),1:10) ]) ; w=rdmds(wfilename,'b'); t=t-meantheta; u=u-meanu; v=v-meanv; w=w-meanw; sumut=sumut+u.*t; sumvt=sumvt+v.*t; sumwt=sumwt+w.*t; end sumut=sumut/counter; sumvt=sumvt/counter; sumwt=sumwt/counter; D=zeros(Nx,Ny,Nz); dx=dm; dy=dm; dz=-dz; D(2:Nx-1,2:Ny-1,2:Nz-1)=-(sumut(3:Nx,2:Ny-1,2:Nz-1)-sumut(1:Nx-2,2:Ny-1,2:Nz-1))/(2*dx) ... -(sumvt(2:Nx-1,3:Ny,2:Nz-1)-sumvt(2:Nx-1,1:Ny-2,2:Nz-1))/(2*dy) ... -(sumwt(2:Nx-1,2:Ny-1,3:Nz)-sumwt(2:Nx-1,2:Ny-1,1:Nz-2))/(2*dz); Dw=zeros(Nx,Ny,Nz); ff=zeros(Nx,Ny,Nz); for j=1:Ny ff(:,j,:)=f+(j-1)*dy*beta; end Dw(:,:,2:Nz-1)=D(:,:,2:Nz-1)*2*dz./(meantheta(:,:,3:Nz)-meantheta(:,:,1:Nz-2)); i=find(meantheta(:,:,3:Nz)==meantheta(:,:,1:Nz-2)); Dw(i)=0; t0=input('Enter temperature : ') % Interpolate Dw on T=T0 surface h=zeros(Nx,Ny); Dh=zeros(Nx,Ny); dz=-dz meantheta(:,:,Nz)=0.; for i=1:Nx for j=1:Ny kk=find(meantheta(i,j,:)1 h(i,j)=(kk(1)-1)*dz+dz*(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1))); Dh(i,j)=Dw(i,j,kk(1)-1) ... +(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ... *(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1))); else h(i,j)=0; Dh(i,j)=0; end end end hmax=max(max(h)) figure pcolor(Dh.');shading flat; colorbar; axis square title('mean Wstar'); % title(['mean Wstar from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]); figure pcolor(h.');shading flat; colorbar; axis square %pcolor(flipud(h.'));shading flat; colorbar; axis square %title(['mean H from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]); we=0; ws=0; for i=1:Nx for j=1:Ny r=sqrt((i-3*Nx/4)*(i-3*Nx/4)+(j-Ny/2)*(j-Ny/2))*dx; %if r<0.44 we=meanw(i,j,3)*dx*dx+we; ws=Dh(i,j)*dx*dx+ws; %end end end we ws