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

% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/eh3/llc/ecco-godae/climatology/t18_to_global_inputs.m,v 1.4 2007/03/28 18:13:38 jmc Exp $
% $Name:  $

% "facet" sizes
nr=360;nb=90;ng=90;
%- set all 6 faces dimensions
nf=ones(6,2);
nf(1,:)=[nb nr]; nf(2,:)=[ng nr]; nf(3,:)=[ng nb];
nf(4,:)=[nr nb]; nf(5,:)=[nr ng]; nf(6,:)=[nb ng];
fdim=prod(nf,2); fd2= cumsum(fdim); fd1=fd2-fdim+1;
mapIO=1; Nfaces=5;

% Fifteent level llc_cs90x90x360
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')

if mapIO >= -1,
% Make a global file from the 18 (or 17) tile files.
lgx=2*nb+2*ng+2*nr;lgy=nr;
% 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];
end

% Bathymetry
if mapIO == -1,
 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');
%-----------------------------------
else
 phi=permute(bdata,[2 3 1]);
 foo=zeros(fd2(end),1);
 ntb=0;
 for n=1:Nfaces
  nt=fdim(n)/nbx/nby;
  vv=phi(:,:,ntb+1:ntb+nt);
  if nf(n,1) > nbx, vv=permute(vv,[1 3 2]); end
 %- save face in structure "bF":
  vv=reshape(vv,nf(n,:)); nvar=sprintf('bF.f%3.3i=vv;',n); eval(nvar)
 %- save compact format in "foo":
  foo(fd1(n):fd2(n),1)=reshape(vv,[fdim(n) 1]);
  ntb=ntb+nt;
 end
 ccB=[-6000 0]; plot_faces(4,bF,0,ccB);
end
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)
if mapIO == -1,
 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');
else
 phi=permute(ptdata,[2 3 1 4]);
 foo=zeros(fd2(end),nz);
 ntb=0;
 for n=1:Nfaces
  nt=fdim(n)/nbx/nby;
  vv=phi(:,:,ntb+1:ntb+nt,:);
  if nf(n,1) > nbx, vv=permute(vv,[1 3 2 4]); end
 %- save face in structure "tF":
  vv=reshape(vv,[nf(n,:) nz]); nvar=sprintf('tF.f%3.3i=vv;',n); eval(nvar)
 %- save compact format in "foo":
  foo(fd1(n):fd2(n),:)=reshape(vv,[fdim(n) nz]);
  ntb=ntb+nt;
 end
 ccB=[0 30]; plot_faces(5,tF,1,ccB);
end
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)
if mapIO == -1,
 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');
else
 phi=permute(sdata,[2 3 1 4]);
 foo=zeros(fd2(end),nz);
 ntb=0;
 for n=1:Nfaces
  nt=fdim(n)/nbx/nby;
  vv=phi(:,:,ntb+1:ntb+nt,:);
  if nf(n,1) > nbx, vv=permute(vv,[1 3 2 4]); end
 %- save face in structure "sF":
  vv=reshape(vv,[nf(n,:) nz]); nvar=sprintf('sF.f%3.3i=vv;',n); eval(nvar)
 %- save compact format in "foo":
  foo(fd1(n):fd2(n),:)=reshape(vv,[fdim(n) nz]);
  ntb=ntb+nt;
 end
 ccB=[31 38]; plot_faces(6,sF,1,ccB);
end
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