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

Annotation of /MITgcm_contrib/timour_matlab/mscripts/DIAG1.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
2     colormap gray
3     cmap=colormap;
4     cmap=flipud(cmap);
5     %cmap=sqrt(cmap);
6     cmap=0.5*(1-cos(3.1415*cmap));
7     colormap(cmap)
8    
9     sumM=zeros(Nx,Ny);
10     sumD=zeros(Nx,Ny);
11     sumW=zeros(Nx,Ny);
12    
13     dz=0.005;
14    
15     % DEFINE BOX OF INTEGRATION
16     Nx1=4;
17     Nx2=Nx-3;
18     Ny1=4;
19     Ny2=Ny-3;
20     Nz1=5;
21     Nz2=Nz-4;
22    
23     pv1=pv;
24     Mp1=Mp;
25     Dp1=Dp;
26    
27     % pv1=smooth3(pv);
28     % Mp1=smooth3(Mp);
29     % Dp1=smooth3(Dp);
30    
31     % for kk=1:4
32     % pv1=smooth3(pv1);
33     % Mp1=smooth3(Mp1);
34     % Dp1=smooth3(Dp1);
35     % end
36    
37     for k=Nz1:Nz2
38     if k==Nz1
39     a=0.5;
40     elseif k==Nz2
41     a=0.5;
42     else
43     a=1;
44     end
45    
46     sumM(3:Nx-2,3:Ny-2)=sumM(3:Nx-2,3:Ny-2)+a*Mp1(3:Nx-2,3:Ny-2,k).*pv1(3:Nx-2,3:Ny-2,k)*dz;
47     sumD(3:Nx-2,3:Ny-2)=sumD(3:Nx-2,3:Ny-2)+a*Dp1(3:Nx-2,3:Ny-2,k).*pv1(3:Nx-2,3:Ny-2,k)*dz;
48     end
49     sumW(3:Nx-2,3:Ny-2)=0.5*meanw(3:Nx-2,3:Ny-2,Nz1).*pv1(3:Nx-2,3:Ny-2,Nz1).*pv1(3:Nx-2,3:Ny-2,Nz1);
50    
51     for k=1:3
52     sumM1(:,:,k)=sumM(:,:);
53     sumD1(:,:,k)=sumD(:,:);
54     sumW1(:,:,k)=sumW(:,:);
55     end
56    
57     for k=1:25
58    
59     sumM1=smooth3(sumM1);
60     sumD1=smooth3(sumD1);
61     sumW1=smooth3(sumW1);
62     end
63    
64     V=[0.25*min(min(sumW(Nx1+2:Nx2,Ny1:Ny2-4))) 0];
65    
66     % title='Vorticity input';
67     % imagesc(lat,long,sumW(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
68     % set(gca,'ydir','norm')
69     % text(0,110,title);
70     %figure
71     v=zeros(10,1);
72     v1=zeros(10,1);
73     for i=1:10
74     %v(i)=-0.1*(i-0.5)*7;
75     %v1(i)=0.1*(i-0.5)*7;
76     v(i)=-0.1*(i)*10;
77     v1(i)=0.1*(i)*10;
78     end
79    
80     %figure
81    
82     % contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v)
83     % hold on
84     % contour(squeeze(sumW1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
85     % hold off
86     % text(20,0,'countour interval is 0.7 in the non-dimensional units');
87     % text(20,95,'Vorticity input','Fontsize',16);
88    
89     %figure
90    
91     %V=[-0.5*max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4)))) 0.5*max(max(abs(sumD(Nx1:Nx2,Ny1:Ny2-4))))];
92     % title='Buoyancy diffusion integral';
93     % imagesc(lat,long,sumD(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
94     % set(gca,'ydir','norm')
95     % text(0,110,title);
96    
97    
98     figure
99     x=Nx1:Nx2;
100     y=Ny1:Ny2;
101    
102     subplot(2,1,1)
103     contourf(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',10)
104     % contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v)
105     % hold on
106     % contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',10);colorbar;
107     % contour(x,y,squeeze(sumD1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
108     text(10,Ny2+5,'Buoyancy eddy-transfer integral','Fontsize',13);
109     hold off
110     xlabel('X ')
111     ylabel('Y ')
112     set(gca,'DataAspectRatio',[2,2,2])
113     set(gca,'XtickLabel','|||||||')
114    
115     %figure
116     %V=[-0.5*max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4)))) 0.5*max(max(abs(sumM(Nx1:Nx2,Ny1:Ny2-4))))];
117     % title='Momentum eddy-transfer integral';
118     % imagesc(lat,long,sumM(Nx1:Nx2,Ny1:Ny2)');shading flat;caxis(V);axis image;colorbar('vertical');
119     % set(gca,'ydir','norm')
120     % text(0,110,title);
121    
122     subplot(2,1,2)
123     xx=Nx1:Nx2;
124     yy=Ny1:Ny2;
125     % contour(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v)
126     contourf(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',10)
127    
128     text(10,Ny2+5,'Momentum eddy-transfer integral','Fontsize',13);
129     %hold on
130     % contour(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',[0,0],' ')
131     % contour(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',v1,'--')
132     % contour(xx,yy,squeeze(sumM1(Nx1:Nx2,Ny1:Ny2,1))',10);colorbar;
133     %hold off
134     xlabel('X (gridpoints)')
135     ylabel('Y (gridpoints)')
136     set(gca,'DataAspectRatio',[2,2,2])
137     % text(40,-25,'countour interval - 0.7 (non-dimensional units)','FontSize',7);
138     % text(40,-30,'negative values - solid line','FontSize',7);
139     % text(40,-35,'positive values - dashed line','FontSize',7);
140    
141    
142    
143     %-----TEST-------------------------------------------------
144    
145     sum=0;
146     sumT=0;
147     sumY=0;
148     for i=Nx1:Nx2
149     for j=Ny1:Ny2
150    
151    
152     if j==Ny1
153     a=0.5;
154     elseif j==Ny2
155     a=0.5;
156     else
157     a=1;
158     end
159    
160     if i==Nx1
161     b=0.5;
162     elseif i==Nx2
163     b=0.5;
164     else
165     b=1;
166     end
167    
168    
169    
170     sum=sum+a*b*0.5*meanw(i,j,Nz1)*pv(i,j,Nz1)*pv(i,j,Nz1)*dx*dy;
171     sumT=sum-0.5*meanw(i,j,Nz2)*pv(i,j,Nz2)*pv(i,j,Nz2)*dx*dy;
172     end
173     end
174    
175     for j=Ny1:Ny2
176     for k=Nz1:Nz2
177     sumT=sumT+0.5*meanu(Nx2,j,k)*pv(Nx2,j,k)*pv(Nx2,j,k)*dz*dy;
178     sumT=sumT-0.5*meanu(Nx1,j,k)*pv(Nx1,j,k)*pv(Nx1,j,k)*dz*dy;
179     end
180     end
181    
182     for i=Nx1:Nx2
183     for k=Nz1:Nz2
184     sumT=sumT+0.5*meanv(i,Ny2,k)*pv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
185     sumY=sumY+0.5*meanv(i,Ny2,k)*pv(i,Ny2,k)*pv(i,Ny2,k)*dz*dx;
186     sumT=sumT-0.5*meanv(i,Ny1,k)*pv(i,Ny1,k)*pv(i,Ny1,k)*dz*dx;
187     end
188     end
189     'VORTICITY GENERATION'
190     sum
191     sumY
192     sumT
193    
194    
195     sum=0;
196     for i=Nx1:Nx2
197     for j=Ny1:Ny2
198    
199     if j==Ny1
200     a=0.5;
201     elseif j==Ny2
202     a=0.5;
203     else
204     a=1;
205     end
206    
207     if i==Nx1
208     b=0.5;
209     elseif i==Nx2
210     b=0.5;
211     else
212     b=1;
213     end
214    
215     sum=sum+a*b*sumD(i,j)*dx*dy;
216     end
217     end
218     'DISSIPATION BY EDDY-DIFFUSIVITY'
219     sum
220    
221     sum=0;
222     for i=Nx1:Nx2
223     for j=Ny1:Ny2
224    
225     if j==Ny1
226     a=0.5;
227     elseif j==Ny2
228     a=0.5;
229     else
230     a=1;
231     end
232    
233     if i==Nx1
234     b=0.5;
235     elseif i==Nx2
236     b=0.5;
237     else
238     b=1;
239     end
240    
241     sum=sum+a*b*sumM(i,j)*dx*dy;
242     end
243     end
244     'DISSIPATION BY EDDY-VISCOSITY'
245     sum
  ViewVC Help
Powered by ViewVC 1.1.22