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

Contents of /MITgcm_contrib/timour_matlab/mscripts/working.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 clear
2 clear path
3
4 global Nx Ny Nz
5 global lat long dz dm mdep
6 global delt_su su_its t_su delt
7 global descriptor this_path
8 global f deltaf Q beta r_expt r_heat H
9 global time rots it
10 global g Cp rho_bar alpha
11 global u v t w
12 global iterations
13
14
15 param_file_name = ...
16 input(' Please enter the name of the m-file with the parameters for this run : ','s') ;
17 feval(param_file_name) ;
18 path = this_path
19 cmdstr=['cd ' path ];
20 eval(cmdstr);
21 path=pwd
22
23 % create circular sampling rings
24
25
26 Nn=Nx/2;
27
28 rad=[1:1:(Nn-1)];
29 ring=ones(Nx,Ny,Nz);
30 weight=zeros(1,Nn-1);
31 radialmean=zeros(1,Nn-1);
32 XCenter = (Nx+1)/3 ;
33 YCenter = (Ny+1)/2 ;
34
35
36 for k=1:Nn-1
37 for i=1:Nx
38 for j=1:Ny
39 radius(i,j) = sqrt((i-XCenter)^2+(j-YCenter)^2);
40 end
41 end
42 radius(find( round(radius) > (k+1) ))=radius(find (round(radius) > (k+1) ))*0;
43 radius(find( round(radius) < (k) ))=radius(find (round(radius) < (k) ))*0;
44 radius(find(radius>0))=1.0;
45 ring(:,:,k)=radius;
46 weight(:,k)=length(find(ring(:,:,k) > 0));
47 end
48
49 sizeit=size(iterations);
50 theta20=zeros(sizeit(1),1);
51 figure;hold on
52 for i=3:sizeit(1)-1
53 tfilename=(['T.' iterations((i),1:10) ]);
54 t=rdmds(tfilename,'b');
55
56 tring=t(:,:,2).*ring(:,:,41);
57 meantee=sum(sum(tring))/weight(:,41);
58
59 tval(i,:)=meantee-20;
60 centralt=tval;
61 plot(eval(iterations((i),1:10)),centralt(i,:),'.');drawnow
62 tincol(i)=sum((centralt(i,:)).*dz);
63
64 %time(i)=(eval(iterations(i,1:10))-su_its)*delt;
65 time(i)=eval(iterations(i,1:10));
66
67
68 end
69
70
71
72 tcol=tval;
73 title(descriptor)
74 xlabel('timestep')
75 ylabel('central temperature')
76 grid
77
78 computertime=time;
79
80 %save centralt centralt
81 %save iterations iterations
82 %save tcol tcol
83 %save tincol tincol
84 %save mdep mdep
85 save working tcol tincol mdep computertime

  ViewVC Help
Powered by ViewVC 1.1.22