ecco-godae/climatology/t18_to_global_inputs.m

% Plot cube fields in six panel layout
function [] = t18_to_global_inputs()

% "facet" sizes
% Fifteent level llc_cs90x90x360
nr=90;
ng=90;
nb=360;
nz=50;

dojpeg = 'y'
doeps  = 'n'

disp('This is version for correct, non-transposed input!')
nz

% Read in the bathymetry
figure(1);clf
bdir='.';
bpref='bathy';
nbx=90;
nby=90;
nbt=0;
bdata=[];
ploc=[36 29 22 15 37 30 23 16 9 10 11 12 13 3 4 5 6 7];
fl=dir(bdir);
for f=1:length(fl)
 n=fl(f).name;
 rpat=sprintf('^%s\\..*',bpref);
 a=regexp(n,rpat,'match');
 if length(a) == 1 & nbt <= 17
  nbt=nbt+1;
  bfile(nbt)=a;
  fn=sprintf('%s/%s',bdir,a{:});
  fid=fopen(fn,'r','ieee-be');
  bdata(nbt,:,:)=fread(fid,[nbx nby],'float64');
  fclose(fid);
  phi=flipud(squeeze(bdata(nbt,:,:)));
  phi(find(phi==0))=NaN;
  subplot(6,7,ploc(nbt));imagesc(fliplr(flipud(-phi')));
  shading flat; axis equal;axis off;axis tight;caxis([0 6000])
 end
end
if doeps =='y'; print -depsc bathy_tile.eps; end

% Read in potential temperature
figure(2);clf
bdir='.';
bpref=['temp_',int2str(nz),'_lev'];
nbx=90;
nby=90;
npt=0;
ptdata=[];
ploc=[36 29 22 15 37 30 23 16 9 10 11 12 13 3 4 5 6 7];
fl=dir(bdir);
for f=1:length(fl)
 n=fl(f).name;
 rpat=sprintf('^%s\\..*',bpref);
 a=regexp(n,rpat,'match');
 if length(a) == 1 & npt <= 17
  npt=npt+1;
  ptfile(npt)=a;
  fn=sprintf('%s/%s',bdir,a{:});
  fid=fopen(fn,'r','ieee-be');
  ptdata(npt,:,:,:)=reshape(fread(fid,nbx*nby*nz,'float64'),[nbx nby nz]);
  fclose(fid);
  phi=flipud(squeeze(ptdata(npt,:,:,1)));
  phiMsk=flipud(squeeze(bdata(npt,:,:)));
  phi(find(phiMsk==0))=NaN;
  subplot(6,7,ploc(npt));imagesc(fliplr(flipud(phi')));
  shading flat; axis equal;axis off;axis tight;caxis([0 30])
 end
end
if doeps =='y'; print -depsc ptemp_tile.eps; end

% Read in sailinity
figure(3);clf
bdir='.';
bpref=['salt_',int2str(nz),'_lev'];
nbx=90;
nby=90;
nsalt=0;
sdata=[];
ploc=[36 29 22 15 37 30 23 16 9 10 11 12 13 3 4 5 6 7];
fl=dir(bdir);
for f=1:length(fl)
 n=fl(f).name;
 rpat=sprintf('^%s\\..*',bpref);
 a=regexp(n,rpat,'match');
 if length(a) == 1 & nsalt <= 17
  nsalt=nsalt+1;
  sfile(nsalt)=a;
  fn=sprintf('%s/%s',bdir,a{:});
  fid=fopen(fn,'r','ieee-be');
  sdata(nsalt,:,:,:)=reshape(fread(fid,nbx*nby*nz,'float64'),[nbx nby nz]);
  fclose(fid);
  phi=flipud(squeeze(sdata(nsalt,:,:,1)));
  phiMsk=flipud(squeeze(bdata(nsalt,:,:)));
  phi(find(phiMsk==0))=NaN;
  subplot(6,7,ploc(nsalt));imagesc(fliplr(flipud(phi')));
  shading flat; axis equal;axis off;axis tight;caxis([31 37])
 end
end
if doeps =='y'; print -depsc salt_tile.eps; end

%%  pause('wait')

% Make a global file from the 18 (or 17) tile files.
lgx=2*nr+2*ng+2*nb;lgy=nb;
% Fseek offsets for the 18 90x90 tiles (numbers are specific to this setup)
tx=90;ty=90;
%   Facet 1 ...
%   Facet 2 ...
%   etc...
ox=[0 0 0 0                       ...
    tx   tx   tx   tx             ...
    2*tx                          ...
    3*tx   4*tx   5*tx    6*tx    ...
    7*tx   8*tx   9*tx   10*tx    ...
    11*tx  ];
oy=[0 ty   2*ty   3*ty            ...
    0 ty   2*ty   3*ty            ...
    0                             ...
    0 0    0      0               ...
    0 0    0      0               ...
    0];

% Bathymetry
foo=zeros(lgx,lgy);
for f=1:size(bdata,1)
 phi=squeeze(bdata(f,:,:));
 if f == 18
  phi=0;
 end

%With transposed input
%for c=1:size(phi,1)
%With correct input
%ph(
for c=1:size(phi,2)
%ph)

  ry=oy(f)+c-1;
  fsloc=ry*lgx+ox(f);

% With transposed input
%  foo(fsloc+1:fsloc+90)=phi(c,:);
% With correct input
  foo(fsloc+1:fsloc+90)=phi(:,c);

 end
end
figure(4);clf
pcolor(-foo');axis equal; shading flat;
axis([0 990 0 360])
caxis([0 6000]); colorbar('South');
title('bathy')
if dojpeg =='y'; print -djpeg100 bathy.jpeg; end
if doeps =='y'; print -depsc bathy.eps; end
fid=fopen('llc90x90x360_bathy.bin','w','ieee-be');
fwrite(fid,foo,'float64');
fclose(fid);

% Potential temperature (adds face 6 too)
foo=zeros(lgx,lgy,nz);
for f=1:size(ptdata,1)
 phi=squeeze(ptdata(f,:,:,:));
 phiMsk=squeeze(bdata(f,:,:));
 phiMsk(find(phiMsk~=0))=1.;

 for k=1:nz
%ph( do not mask, keep zonal mean values everywhere
%ph  phi(:,:,k)=phi(:,:,k).*phiMsk;
  phi(:,:,k)=phi(:,:,k);
%ph)

% With transposed input
%  for c=1:size(phi,1)
% With correct input
  for c=1:size(phi,2)

   ry=oy(f)+c-1;
   fsloc=(k-1)*lgx*lgy+ry*lgx+ox(f);
%  With transposed input
%   foo(fsloc+1:fsloc+90)=phi(c,:,k);
%  With correct input
  foo(fsloc+1:fsloc+90)=phi(:,c,k);

  end
 end
end
figure(5);clf
pcolor(squeeze(foo(:,:,1))');axis equal;  shading flat;
axis([0 990 0 360])
caxis([0. 30.]); colorbar('South');
title('ptemp')
if dojpeg =='y'; print -djpeg100 ptemp.jpeg; end
if doeps =='y'; print -depsc ptemp.eps; end
fid=fopen(['llc90x90x360_',int2str(nz),'lev_ptemp.bin'],'w','ieee-be');
fwrite(fid,foo,'float64');
fclose(fid);

% Salinity (adds face 6 too)
foo=zeros(lgx,lgy,nz);
for f=1:size(sdata,1)
 phi=squeeze(sdata(f,:,:,:));
 phiMsk=squeeze(bdata(f,:,:));
 phiMsk(find(phiMsk~=0))=1.;

 for k=1:nz
%ph( do not mask, keep zonal mean values everywhere
%ph  phi(:,:,k)=phi(:,:,k).*phiMsk;
  phi(:,:,k)=phi(:,:,k);
%ph)

% With transposed input
%  for c=1:size(phi,1)
% With correct input
   for c=1:size(phi,2)

   ry=oy(f)+c-1;
   fsloc=(k-1)*lgx*lgy+ry*lgx+ox(f);
%  With transposed input
%   foo(fsloc+1:fsloc+90)=phi(c,:,k);
%  With correct input
   foo(fsloc+1:fsloc+90)=phi(:,c,k);

  end
 end
end
figure(6);clf
pcolor(squeeze(foo(:,:,1))');axis equal;  shading flat;
axis([0 990 0 360])
caxis([31. 38.]); colorbar('South');
title('salt')
if dojpeg =='y'; print -djpeg100 salt.jpeg; end
if doeps =='y'; print -depsc salt.eps; end
fid=fopen(['llc90x90x360_',int2str(nz),'lev_salt.bin'],'w','ieee-be');
fwrite(fid,foo,'float64');
fclose(fid);

end

function [] = mkcubeplot(ptit,temp,sf)
ploc=[9 10 6 7 3 4];
phil=[];
for g=1:length(temp)
 phi=temp{g}*sf;
 phil=[ phil phi(:)'];
end
for g=1:length(temp)
 phi=temp{g}*sf;
 phi=flipud(phi');
 subplot(3,4,ploc(g));
 imagesc(phi);caxis([min(phil) max(phil)]);
end
subplot(3,4,1);
axis off
caxis([min(phil) max(phil)]);colorbar;
title(ptit);
drawnow

pnam=sprintf('%s.jpeg',ptit);
print('-djpeg100',pnam);

end
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	heimbach	1.3	doeps = 'n'
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,nbxnbynz,'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,nbxnbynz,'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=2nr+2ng+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			3tx 4tx 5tx 6tx ...
119			7tx 8tx 9tx 10tx ...
120			11*tx ];
121			oy=[0 ty 2ty 3ty ...
122			0 ty 2ty 3ty ...
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	heimbach	1.3	%ph( do not mask, keep zonal mean values everywhere
173			%ph phi(:,:,k)=phi(:,:,k).*phiMsk;
174			phi(:,:,k)=phi(:,:,k);
175			%ph)
176	heimbach	1.1
177			% With transposed input
178			% for c=1:size(phi,1)
179			% With correct input
180			for c=1:size(phi,2)
181
182			ry=oy(f)+c-1;
183			fsloc=(k-1)lgxlgy+ry*lgx+ox(f);
184			% With transposed input
185			% foo(fsloc+1:fsloc+90)=phi(c,:,k);
186			% With correct input
187			foo(fsloc+1:fsloc+90)=phi(:,c,k);
188
189			end
190			end
191			end
192			figure(5);clf
193	heimbach	1.2	pcolor(squeeze(foo(:,:,1))');axis equal; shading flat;
194			axis([0 990 0 360])
195			caxis([0. 30.]); colorbar('South');
196			title('ptemp')
197			if dojpeg =='y'; print -djpeg100 ptemp.jpeg; end
198			if doeps =='y'; print -depsc ptemp.eps; end
199	heimbach	1.3	fid=fopen(['llc90x90x360_',int2str(nz),'lev_ptemp.bin'],'w','ieee-be');
200	heimbach	1.1	fwrite(fid,foo,'float64');
201			fclose(fid);
202
203			% Salinity (adds face 6 too)
204			foo=zeros(lgx,lgy,nz);
205			for f=1:size(sdata,1)
206			phi=squeeze(sdata(f,:,:,:));
207			phiMsk=squeeze(bdata(f,:,:));
208			phiMsk(find(phiMsk~=0))=1.;
209
210			for k=1:nz
211	heimbach	1.3	%ph( do not mask, keep zonal mean values everywhere
212			%ph phi(:,:,k)=phi(:,:,k).*phiMsk;
213			phi(:,:,k)=phi(:,:,k);
214			%ph)
215	heimbach	1.1
216			% With transposed input
217			% for c=1:size(phi,1)
218			% With correct input
219			for c=1:size(phi,2)
220
221			ry=oy(f)+c-1;
222			fsloc=(k-1)lgxlgy+ry*lgx+ox(f);
223			% With transposed input
224			% foo(fsloc+1:fsloc+90)=phi(c,:,k);
225			% With correct input
226			foo(fsloc+1:fsloc+90)=phi(:,c,k);
227
228			end
229			end
230			end
231			figure(6);clf
232	heimbach	1.2	pcolor(squeeze(foo(:,:,1))');axis equal; shading flat;
233			axis([0 990 0 360])
234			caxis([31. 38.]); colorbar('South');
235			title('salt')
236			if dojpeg =='y'; print -djpeg100 salt.jpeg; end
237			if doeps =='y'; print -depsc salt.eps; end
238	heimbach	1.3	fid=fopen(['llc90x90x360_',int2str(nz),'lev_salt.bin'],'w','ieee-be');
239	heimbach	1.1	fwrite(fid,foo,'float64');
240			fclose(fid);
241
242			end
243
244			function [] = mkcubeplot(ptit,temp,sf)
245			ploc=[9 10 6 7 3 4];
246			phil=[];
247			for g=1:length(temp)
248			phi=temp{g}*sf;
249			phil=[ phil phi(:)'];
250			end
251			for g=1:length(temp)
252			phi=temp{g}*sf;
253			phi=flipud(phi');
254			subplot(3,4,ploc(g));
255			imagesc(phi);caxis([min(phil) max(phil)]);
256			end
257			subplot(3,4,1);
258			axis off
259			caxis([min(phil) max(phil)]);colorbar;
260			title(ptit);
261			drawnow
262
263			pnam=sprintf('%s.jpeg',ptit);
264			print('-djpeg100',pnam);
265
266			end