/[MITgcm]/MITgcm_contrib/eh3/llc/ecco-godae/climatology/t18_to_global_inputs.m
ViewVC logotype

Annotation of /MITgcm_contrib/eh3/llc/ecco-godae/climatology/t18_to_global_inputs.m

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


Revision 1.2 - (hide annotations) (download)
Fri Nov 17 00:28:05 2006 UTC (18 years, 8 months ago) by heimbach
Branch: MAIN
Changes since 1.1: +38 -14 lines
Generalizing to nz param., changed plots.

1 heimbach 1.1 % Plot cube fields in six panel layout
2     function [] = t18_to_global_inputs()
3    
4     % "facet" sizes
5     % Fifteent level llc_cs90x90x360
6     nr=90;
7     ng=90;
8     nb=360;
9 heimbach 1.2 nz=50;
10    
11     dojpeg = 'y'
12     doeps = 'y'
13 heimbach 1.1
14     disp('This is version for correct, non-transposed input!')
15 heimbach 1.2 nz
16 heimbach 1.1
17     % Read in the bathymetry
18     figure(1);clf
19     bdir='.';
20     bpref='bathy';
21     nbx=90;
22     nby=90;
23     nbt=0;
24     bdata=[];
25     ploc=[36 29 22 15 37 30 23 16 9 10 11 12 13 3 4 5 6 7];
26     fl=dir(bdir);
27     for f=1:length(fl)
28     n=fl(f).name;
29     rpat=sprintf('^%s\\..*',bpref);
30     a=regexp(n,rpat,'match');
31     if length(a) == 1 & nbt <= 17
32     nbt=nbt+1;
33     bfile(nbt)=a;
34     fn=sprintf('%s/%s',bdir,a{:});
35     fid=fopen(fn,'r','ieee-be');
36     bdata(nbt,:,:)=fread(fid,[nbx nby],'float64');
37     fclose(fid);
38     phi=flipud(squeeze(bdata(nbt,:,:)));
39     phi(find(phi==0))=NaN;
40 heimbach 1.2 subplot(6,7,ploc(nbt));imagesc(fliplr(flipud(-phi')));
41     shading flat; axis equal;axis off;axis tight;caxis([0 6000])
42 heimbach 1.1 end
43     end
44 heimbach 1.2 if doeps =='y'; print -depsc bathy_tile.eps; end
45 heimbach 1.1
46     % Read in potential temperature
47     figure(2);clf
48     bdir='.';
49 heimbach 1.2 bpref=['temp_',int2str(nz),'_lev'];
50 heimbach 1.1 nbx=90;
51     nby=90;
52     npt=0;
53     ptdata=[];
54     ploc=[36 29 22 15 37 30 23 16 9 10 11 12 13 3 4 5 6 7];
55     fl=dir(bdir);
56     for f=1:length(fl)
57     n=fl(f).name;
58     rpat=sprintf('^%s\\..*',bpref);
59     a=regexp(n,rpat,'match');
60     if length(a) == 1 & npt <= 17
61     npt=npt+1;
62     ptfile(npt)=a;
63     fn=sprintf('%s/%s',bdir,a{:});
64     fid=fopen(fn,'r','ieee-be');
65     ptdata(npt,:,:,:)=reshape(fread(fid,nbx*nby*nz,'float64'),[nbx nby nz]);
66     fclose(fid);
67     phi=flipud(squeeze(ptdata(npt,:,:,1)));
68     phiMsk=flipud(squeeze(bdata(npt,:,:)));
69     phi(find(phiMsk==0))=NaN;
70 heimbach 1.2 subplot(6,7,ploc(npt));imagesc(fliplr(flipud(phi')));
71     shading flat; axis equal;axis off;axis tight;caxis([0 30])
72 heimbach 1.1 end
73     end
74 heimbach 1.2 if doeps =='y'; print -depsc ptemp_tile.eps; end
75 heimbach 1.1
76     % Read in sailinity
77     figure(3);clf
78     bdir='.';
79 heimbach 1.2 bpref=['salt_',int2str(nz),'_lev'];
80 heimbach 1.1 nbx=90;
81     nby=90;
82     nsalt=0;
83     sdata=[];
84     ploc=[36 29 22 15 37 30 23 16 9 10 11 12 13 3 4 5 6 7];
85     fl=dir(bdir);
86     for f=1:length(fl)
87     n=fl(f).name;
88     rpat=sprintf('^%s\\..*',bpref);
89     a=regexp(n,rpat,'match');
90     if length(a) == 1 & nsalt <= 17
91     nsalt=nsalt+1;
92     sfile(nsalt)=a;
93     fn=sprintf('%s/%s',bdir,a{:});
94     fid=fopen(fn,'r','ieee-be');
95     sdata(nsalt,:,:,:)=reshape(fread(fid,nbx*nby*nz,'float64'),[nbx nby nz]);
96     fclose(fid);
97     phi=flipud(squeeze(sdata(nsalt,:,:,1)));
98     phiMsk=flipud(squeeze(bdata(nsalt,:,:)));
99     phi(find(phiMsk==0))=NaN;
100 heimbach 1.2 subplot(6,7,ploc(nsalt));imagesc(fliplr(flipud(phi')));
101     shading flat; axis equal;axis off;axis tight;caxis([31 37])
102 heimbach 1.1 end
103     end
104 heimbach 1.2 if doeps =='y'; print -depsc salt_tile.eps; end
105    
106     %% pause('wait')
107 heimbach 1.1
108     % Make a global file from the 18 (or 17) tile files.
109     lgx=2*nr+2*ng+2*nb;lgy=nb;
110     % Fseek offsets for the 18 90x90 tiles (numbers are specific to this setup)
111     tx=90;ty=90;
112     % Facet 1 ...
113     % Facet 2 ...
114     % etc...
115     ox=[0 0 0 0 ...
116     tx tx tx tx ...
117     2*tx ...
118     3*tx 4*tx 5*tx 6*tx ...
119     7*tx 8*tx 9*tx 10*tx ...
120     11*tx ];
121     oy=[0 ty 2*ty 3*ty ...
122     0 ty 2*ty 3*ty ...
123     0 ...
124     0 0 0 0 ...
125     0 0 0 0 ...
126     0];
127    
128     % Bathymetry
129     foo=zeros(lgx,lgy);
130     for f=1:size(bdata,1)
131     phi=squeeze(bdata(f,:,:));
132     if f == 18
133     phi=0;
134     end
135    
136     %With transposed input
137     %for c=1:size(phi,1)
138     %With correct input
139     %ph(
140     for c=1:size(phi,2)
141     %ph)
142    
143     ry=oy(f)+c-1;
144     fsloc=ry*lgx+ox(f);
145    
146     % With transposed input
147     % foo(fsloc+1:fsloc+90)=phi(c,:);
148     % With correct input
149     foo(fsloc+1:fsloc+90)=phi(:,c);
150    
151     end
152     end
153     figure(4);clf
154 heimbach 1.2 pcolor(-foo');axis equal; shading flat;
155     axis([0 990 0 360])
156     caxis([0 6000]); colorbar('South');
157     title('bathy')
158     if dojpeg =='y'; print -djpeg100 bathy.jpeg; end
159     if doeps =='y'; print -depsc bathy.eps; end
160 heimbach 1.1 fid=fopen('llc90x90x360_bathy.bin','w','ieee-be');
161     fwrite(fid,foo,'float64');
162     fclose(fid);
163    
164     % Potential temperature (adds face 6 too)
165     foo=zeros(lgx,lgy,nz);
166     for f=1:size(ptdata,1)
167     phi=squeeze(ptdata(f,:,:,:));
168     phiMsk=squeeze(bdata(f,:,:));
169     phiMsk(find(phiMsk~=0))=1.;
170    
171     for k=1:nz
172     phi(:,:,k)=phi(:,:,k).*phiMsk;
173    
174     % With transposed input
175     % for c=1:size(phi,1)
176     % With correct input
177     for c=1:size(phi,2)
178    
179     ry=oy(f)+c-1;
180     fsloc=(k-1)*lgx*lgy+ry*lgx+ox(f);
181     % With transposed input
182     % foo(fsloc+1:fsloc+90)=phi(c,:,k);
183     % With correct input
184     foo(fsloc+1:fsloc+90)=phi(:,c,k);
185    
186     end
187     end
188     end
189     figure(5);clf
190 heimbach 1.2 pcolor(squeeze(foo(:,:,1))');axis equal; shading flat;
191     axis([0 990 0 360])
192     caxis([0. 30.]); colorbar('South');
193     title('ptemp')
194     if dojpeg =='y'; print -djpeg100 ptemp.jpeg; end
195     if doeps =='y'; print -depsc ptemp.eps; end
196     fid=fopen(['llc90x90x360_',int2str(nz),'l1v_ptemp.bin'],'w','ieee-be');
197 heimbach 1.1 fwrite(fid,foo,'float64');
198     fclose(fid);
199    
200     % Salinity (adds face 6 too)
201     foo=zeros(lgx,lgy,nz);
202     for f=1:size(sdata,1)
203     phi=squeeze(sdata(f,:,:,:));
204     phiMsk=squeeze(bdata(f,:,:));
205     phiMsk(find(phiMsk~=0))=1.;
206    
207     for k=1:nz
208     phi(:,:,k)=phi(:,:,k).*phiMsk;
209    
210     % With transposed input
211     % for c=1:size(phi,1)
212     % With correct input
213     for c=1:size(phi,2)
214    
215     ry=oy(f)+c-1;
216     fsloc=(k-1)*lgx*lgy+ry*lgx+ox(f);
217     % With transposed input
218     % foo(fsloc+1:fsloc+90)=phi(c,:,k);
219     % With correct input
220     foo(fsloc+1:fsloc+90)=phi(:,c,k);
221    
222     end
223     end
224     end
225     figure(6);clf
226 heimbach 1.2 pcolor(squeeze(foo(:,:,1))');axis equal; shading flat;
227     axis([0 990 0 360])
228     caxis([31. 38.]); colorbar('South');
229     title('salt')
230     if dojpeg =='y'; print -djpeg100 salt.jpeg; end
231     if doeps =='y'; print -depsc salt.eps; end
232     fid=fopen(['llc90x90x360_',int2str(nz),'l1v_salt.bin'],'w','ieee-be');
233 heimbach 1.1 fwrite(fid,foo,'float64');
234     fclose(fid);
235    
236     end
237    
238     function [] = mkcubeplot(ptit,temp,sf)
239     ploc=[9 10 6 7 3 4];
240     phil=[];
241     for g=1:length(temp)
242     phi=temp{g}*sf;
243     phil=[ phil phi(:)'];
244     end
245     for g=1:length(temp)
246     phi=temp{g}*sf;
247     phi=flipud(phi');
248     subplot(3,4,ploc(g));
249     imagesc(phi);caxis([min(phil) max(phil)]);
250     end
251     subplot(3,4,1);
252     axis off
253     caxis([min(phil) max(phil)]);colorbar;
254     title(ptit);
255     drawnow
256    
257     pnam=sprintf('%s.jpeg',ptit);
258     print('-djpeg100',pnam);
259    
260     end

  ViewVC Help
Powered by ViewVC 1.1.22