/[MITgcm]/MITgcm_contrib/timour_matlab/mscripts/pvort.m
ViewVC logotype

Annotation of /MITgcm_contrib/timour_matlab/mscripts/pvort.m

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


Revision 1.1 - (hide annotations) (download)
Wed Sep 3 21:22:22 2003 UTC (21 years, 10 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
initial checkin of Timour's MatLAB scripts

1 edhill 1.1 function plotpv( iter)
2     clear
3     clear path
4    
5     global Nx Ny Nz
6     global lat long dz dm mdep
7     global delt_su su_its t_su delt
8     global descriptor this_path
9     global f deltaf Q beta r_expt r_heat H
10     global time rots it
11     global g Cp rho_bar alpha
12     global u v t w
13     global iterations
14    
15    
16     param_file_name = ...
17     input(' Please enter the name of the m-file with the parameters for this run : ','s') ;
18     feval(param_file_name) ;
19    
20     iterations
21    
22     it = input(' Please enter iteration : ','s')
23    
24     path = this_path
25     cmdstr=['cd ' path ];
26     eval(cmdstr);
27     path=pwd
28    
29     ufilename=(['U.' it ]);
30     vfilename=(['V.' it ]);
31     tfilename=(['T.' it ]);
32     u=rdmds(ufilename,'b');
33     v=rdmds(vfilename,'b');
34     t=rdmds(tfilename,'b');
35    
36     pv=zeros(Nx,Ny,Nz);
37     ff=zeros(Nx,Ny,Nz);
38     dx=dm;
39     dy=dm;
40    
41     for j=1:Ny
42     ff(:,j,:)=f+(j-1)*dy*beta;
43     end
44    
45     pv(1:Nx-1,1:Ny-1,1:Nz-1)=(v(2:Nx,1:Ny-1,1:Nz-1)/dx -v(1:Nx-1,1:Ny-1,1:Nz-1)/dx ...
46     -u(1:Nx-1,2:Ny,1:Nz-1)/dy +u(1:Nx-1,1:Ny-1,1:Nz-1)/dy ...
47     +ff(1:Nx-1,1:Ny-1,1:Nz-1)) ...
48     .*(t(1:Nx-1,1:Ny-1,1:Nz-1)-t(1:Nx-1,1:Ny-1,2:Nz))/dz;
49    
50     t0=input('Enter temperature : ')
51    
52     % Interpolate PV on T=T0 surface
53    
54     h=zeros(Nz,Ny);
55     pvh=zeros(Nz,Ny);
56    
57     for i=1:Nx
58     for j=1:Ny
59     kk=find(t(i,j,:)<t0);
60    
61     if kk(1)>1
62     h(i,j)=(kk(1)-1)*dz+dz*(t(i,j,kk(1)-1)-t0)/(t(i,j,kk(1)-1)-t(i,j,kk(1)));
63     pvh(i,j)=pv(i,j,kk(1)-1) ...
64     +(pv(i,j,kk(1))-pv(i,j,kk(1)-1))*(t(i,j,kk(1)-1)-t0)/(t(i,j,kk(1)-1)-t(i,j,kk(1)));
65     else
66     h(i,j)=0;
67     pvh(i,j)=0;
68     end
69    
70     end
71     end
72    
73     hmax=max(max(h))
74    
75    
76     %cmin=0;
77     %cmax=0.1;
78     %V=[cmin,cmax];
79     %contourf(h,50);caxis(V);colorbar;grid
80    
81     % title='depth of a density interface';
82     % imagesc(lat,long,h);shading flat;axis image;caxis(V);colorbar('horizontal');
83     % text(0,-30,descriptor);
84     % text(0,140,title);
85    
86     figure
87     pcolor(pvh.');shading flat; colorbar; axis square
88     title(['pv timestep ' num2str(eval(it)) ' level of t=' num2str(t0) ]);
89     drawnow
90    
91     figure
92    
93     pcolor(h.');shading flat; colorbar; axis square
94     title(['h timestep' num2str(eval(it)) ' level of t=' num2str(t0) ]);
95     drawnow
96    
97     return

  ViewVC Help
Powered by ViewVC 1.1.22