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

Annotation of /MITgcm_contrib/timour_matlab/mscripts/DIAGW1.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 t0=input('Enter temperature : ')
2    
3     % Interpolate Dw on T=T0 surface
4    
5     h=zeros(Nx,Ny);
6     hh=zeros(Nx,Ny);
7     hhh=zeros(Nx,Ny);
8     Dh=zeros(Nx,Ny);
9    
10     Nz1=4;
11     he=dz*(Nz1-1);
12    
13     meantheta(:,:,Nz)=0.;
14     for i=1:Nx
15     for j=1:Ny
16     kk=find(meantheta(i,j,:)<t0);
17    
18     if kk(1)>1
19     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)));
20    
21     Dh(i,j)=Dw(i,j,kk(1)-1) ...
22     +(Dw(i,j,kk(1))-Dw(i,j,kk(1)-1)) ...
23     *(meantheta(i,j,kk(1)-1)-t0)/(meantheta(i,j,kk(1)-1)-meantheta(i,j,kk(1)));
24     else
25     h(i,j)=0;
26     Dh(i,j)=0;
27     end
28    
29     end
30     end
31    
32     hmax=max(max(h))
33    
34    
35     figure
36    
37     %pcolor(Dh.');shading flat; colorbar; axis square
38     % title('mean Wstar');
39    
40     %figure
41    
42     %pcolor(hhh.');shading flat; colorbar; axis square
43     %pcolor(flipud(h.'));shading flat; colorbar; axis square
44     %title(['mean H from ' num2str(eval(itstart)) ' to 'num2str(eval(itend)) ' level of t=' num2str(t0) ]);
45     %figure
46    
47    
48     we=0;
49     ws=0;
50     WE=zeros(Nx,Ny);
51     WS=zeros(Nx,Ny);
52     hh=zeros(Nx,Ny);
53    
54     %for i=Nx/2:Nx
55     for i=1:Nx
56     for j=1:Ny
57    
58     r=sqrt((i-3*Nx/4)*(i-3*Nx/4)+(j-Ny/2)*(j-Ny/2))*dx;
59     WE(i,j)=meanw(i,j,Nz1);
60     WS(i,j)=Dh(i,j);
61    
62     % if r<0.40
63     if h(i,j)>he
64     we=meanw(i,j,Nz1)*dx*dx+we;
65     ws=Dh(i,j)*dx*dx+ws;
66     % end
67     end
68     end
69     end
70    
71     hh=h;
72    
73     for k=1:3
74     WE1(:,:,k)=WE(:,:);
75     WS1(:,:,k)=WS(:,:);
76     hh1(:,:,k)=hh(:,:);
77     end
78    
79    
80     for k=1:5
81     WE1=smooth3(WE1);
82     WS1=smooth3(WS1);
83     hh1=smooth3(hh1);
84     end
85    
86     hhh(2:Nx-1,2:Ny-1)=10000*(hh1(1:Nx-2,2:Ny-1,1)+hh1(3:Nx,2:Ny-1,1)+...
87     hh1(2:Nx-1,1:Ny-2,1)+hh1(2:Nx-1,3:Ny,1)-4*hh1(2:Nx-1,2:Ny-1,1));
88    
89     %hhh(2:Nx-1,2:Ny-1)=10000*(hh1(1:Nx-2,2:Ny-1,1)-hh1(3:Nx,2:Ny-1,1)).*...
90     %(hh1(1:Nx-2,2:Ny-1,1)-hh1(3:Nx,2:Ny-1,1))+...
91     %10000*(hh1(2:Nx-1,1:Ny-2,1)-hh1(2:Nx-1,3:Ny,1)).*...
92     %(hh1(2:Nx-1,1:Ny-2,1)-hh1(2:Nx-1,3:Ny,1));
93    
94     i=find(h<0.015);
95     h(i)=0;
96     % i=find(abs(hhh)>0.001);
97     % hhh(i)=0;
98    
99     for i=1:Nx
100     for j=1:Ny
101     %r=sqrt((i-3*Nx/4)*(i-3*Nx/4)+(j-Ny/2)*(j-Ny/2))*dx;
102    
103     %if r>0.44
104     %hhh(i,j)=NaN;
105     %hh(i,j)=NaN;
106     %WS1(i,j,1)=NaN;
107     %end
108    
109     if h(i,j)<he
110     WE1(i,j,1)=0;
111     WS1(i,j,1)=NaN;
112     hh(i,j)=NaN;
113     hhh(i,j)=NaN;
114     end
115    
116     end
117     end
118    
119    
120     %pcolor(hh(1:Nx,:).');shading flat; colorbar; axis square
121     % set(gca,'DataAspectRatio',[2,2,2])
122     % title('h','FontSize',16);
123     %figure
124    
125     x=1:3:Nx-2;
126     y=1:3:Ny;
127     hh(:,:)=max(he,hh(:,:));
128     subplot('position',[0.30 0.7 0.4 0.3])
129     surf(x,y,200*(-hh(1:3:Nx-2,1:3:Ny)'+he));
130     axis([1 Nx-4 0 Ny-2 -10.0 0])
131     a=0:1:63;
132     a=a/63;
133     cmap=[a' a' a'];
134     colormap(cmap);
135     set(gca,'DataAspectRatio',[20,20,7])
136     set(gca,'XtickLabel','|')
137     set(gca,'YtickLabel','|')
138     set(gca,'ZtickLabel','|')
139     xlabel('x','FontSize',20)
140     ylabel('y','FontSize',20)
141     zlabel('z','FontSize',20,'Rotation',[0])
142     %title('h','FontSize',19);
143     view(-20,10);
144    
145     hhh(1:Nx,:)=-10*hhh(1:Nx,:);
146     % hhh(Nx/2:Nx,:)=100*hhh(Nx/2:Nx,:);
147     a=max(max(hhh(1:Nx,:)));
148    
149     %pcolor(hhh(1:Nx,:).');caxis([0 a]);shading flat; colorbar; axis square
150     %set(gca,'DataAspectRatio',[2,2,2])
151     %title('-\nabla^2 h','FontSize',12);
152     %xlabel('x','Fontsize',13);
153     %ylabel('y','Fontsize',13);
154     %text(50,-10,'x');
155     %text(-10,60,'y');
156    
157     %pcolor(squeeze(WE1(Nx/2:Nx,:,1)).');shading flat; colorbar; axis square
158     % set(gca,'DataAspectRatio',[2,2,2])
159     %title('Ek, smooth');
160     %figure
161    
162     %pcolor(meanw(Nx/2:Nx,:,3).');shading flat; colorbar; axis square
163     %set(gca,'DataAspectRatio',[2,2,2])
164     %title('Ek');
165     figure
166    
167     subplot('position',[0.15 0.05 0.7 0.5])
168     WS1(1:Nx,:,1)=-WS1(1:Nx,:,1)*16667;
169     a=max(max(WS1(1:Nx,:,1)));
170     pcolor(squeeze(WS1(1:Nx,:,1)).');caxis([0 a]);shading flat;colormap('default');
171     set(gca,'DataAspectRatio',[2,2,2])
172     set(gca,'Xtick',[1 100],'XtickLabel','|')
173     set(gca,'Ytick',[1 120],'YtickLabel','|')
174     title('w^*','FontSize',20);
175     xlabel('x','Fontsize',20);
176     ylabel('y','Fontsize',20,'Rotation',[0]);
177    
178    
179     we
180     ws
181    

  ViewVC Help
Powered by ViewVC 1.1.22