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

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

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


Revision 1.1 - (show 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 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