/[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.4 - (hide annotations) (download)
Wed Mar 28 18:13:38 2007 UTC (18 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +113 -76 lines
- use standard convention for nr,ng,nb (same as in utils/exch2);
- remove function mkcubeplot (not used; could put it back in a separated .m file)
- add option to write global input files in compact format (=the default);
  (set mapIO=-1 to return to old sparse format).

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

  ViewVC Help
Powered by ViewVC 1.1.22