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

Contents of /MITgcm_contrib/timour_matlab/mscripts/DIAG1.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
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