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

Contents 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 - (show 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 % Plot cube fields in six panel layout
2 function [] = t18_to_global_inputs()
3
4 % $Header: $
5 % $Name: $
6
7 % "facet" sizes
8 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 % Fifteent level llc_cs90x90x360
17 nz=50;
18
19 dojpeg = 'y'
20 doeps = 'n'
21
22 disp('This is version for correct, non-transposed input!')
23 nz
24
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 subplot(6,7,ploc(nbt));imagesc(fliplr(flipud(-phi')));
49 shading flat; axis equal;axis off;axis tight;caxis([0 6000])
50 end
51 end
52 if doeps =='y'; print -depsc bathy_tile.eps; end
53
54 % Read in potential temperature
55 figure(2);clf
56 bdir='.';
57 bpref=['temp_',int2str(nz),'_lev'];
58 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 subplot(6,7,ploc(npt));imagesc(fliplr(flipud(phi')));
79 shading flat; axis equal;axis off;axis tight;caxis([0 30])
80 end
81 end
82 if doeps =='y'; print -depsc ptemp_tile.eps; end
83
84 % Read in sailinity
85 figure(3);clf
86 bdir='.';
87 bpref=['salt_',int2str(nz),'_lev'];
88 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 subplot(6,7,ploc(nsalt));imagesc(fliplr(flipud(phi')));
109 shading flat; axis equal;axis off;axis tight;caxis([31 37])
110 end
111 end
112 if doeps =='y'; print -depsc salt_tile.eps; end
113
114 %% pause('wait')
115
116 if mapIO >= -1,
117 % Make a global file from the 18 (or 17) tile files.
118 lgx=2*nb+2*ng+2*nr;lgy=nr;
119 % 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 end
137
138 % Bathymetry
139 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
147 %With transposed input
148 %for c=1:size(phi,1)
149 %With correct input
150 %ph(
151 for c=1:size(phi,2)
152 %ph)
153
154 ry=oy(f)+c-1;
155 fsloc=ry*lgx+ox(f);
156
157 % With transposed input
158 % foo(fsloc+1:fsloc+90)=phi(c,:);
159 % With correct input
160 foo(fsloc+1:fsloc+90)=phi(:,c);
161
162 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 end
182 ccB=[-6000 0]; plot_faces(4,bF,0,ccB);
183 end
184 title('bathy')
185 if dojpeg =='y'; print -djpeg100 bathy.jpeg; end
186 if doeps =='y'; print -depsc bathy.eps; end
187 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 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
205 % With transposed input
206 % for c=1:size(phi,1)
207 % With correct input
208 for c=1:size(phi,2)
209
210 ry=oy(f)+c-1;
211 fsloc=(k-1)*lgx*lgy+ry*lgx+ox(f);
212 % With transposed input
213 % foo(fsloc+1:fsloc+90)=phi(c,:,k);
214 % With correct input
215 foo(fsloc+1:fsloc+90)=phi(:,c,k);
216
217 end
218 end
219 end
220 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 end
240 title('ptemp')
241 if dojpeg =='y'; print -djpeg100 ptemp.jpeg; end
242 if doeps =='y'; print -depsc ptemp.eps; end
243 fid=fopen(['llc90x90x360_',int2str(nz),'lev_ptemp.bin'],'w','ieee-be');
244 fwrite(fid,foo,'float64');
245 fclose(fid);
246
247 % Salinity (adds face 6 too)
248 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
255 for k=1:nz
256 %ph( do not mask, keep zonal mean values everywhere
257 %ph phi(:,:,k)=phi(:,:,k).*phiMsk;
258 phi(:,:,k)=phi(:,:,k);
259 %ph)
260
261 % With transposed input
262 % for c=1:size(phi,1)
263 % With correct input
264 for c=1:size(phi,2)
265
266 ry=oy(f)+c-1;
267 fsloc=(k-1)*lgx*lgy+ry*lgx+ox(f);
268 % With transposed input
269 % foo(fsloc+1:fsloc+90)=phi(c,:,k);
270 % With correct input
271 foo(fsloc+1:fsloc+90)=phi(:,c,k);
272
273 end
274 end
275 end
276 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 end
296 title('salt')
297 if dojpeg =='y'; print -djpeg100 salt.jpeg; end
298 if doeps =='y'; print -depsc salt.eps; end
299 fid=fopen(['llc90x90x360_',int2str(nz),'lev_salt.bin'],'w','ieee-be');
300 fwrite(fid,foo,'float64');
301 fclose(fid);
302
303 end

  ViewVC Help
Powered by ViewVC 1.1.22