% This is a script that executes a series of commands for
% extracting Levitus data, interpolating, and gap-filling
% according to the mask-file generated by "topo.m".
%
% Uses data files: bathymetry.bin pmask.bin
% Creates data files:  lev_clim_*.bin lev_monthly_*.bin
% Creates arrays: n/a
% Uses m-scripts: setgrid grid_sphere finegrid gen_hxhy xyrecur_*
%                 hfac
%
% Created  08/15/99 by adcroft@mit.edu
% Modified 11/09/99 by adcroft@mit.edu
% Maintained by adcroft@mit.edu, abiastoch@ucsd.edu

skip_interp=0

if ~skip_interp
%clear all
format compact
more off

load FMT
load HGRID.mat lon_lo lon_hi lat_lo lat_hi xc yc GRID
dd=0;
load VGRID.mat zc
load MASKS
GRID.mskc=msk;
clear msk

disp('Reading Levitus climatological T');
[TC,xl,yl,zlc]=...
  extract_levitus_clim('temp',lon_lo-dd,lon_hi+dd,lat_lo-dd,lat_hi+dd);
if size(xc,3)==1
  tc=interp3ddata(xl,yl,zlc,TC,xc,yc,zc,GRID.mskc);
  tc2=interp3ddata_int2(xl,yl,zlc,TC,xc,yc,zc,GRID.mskc);
else
 for t=1:size(xc,3);
  tc(:,:,t,:)=interp3ddata(xl,yl,zlc,TC,xc(:,:,t),yc(:,:,t),zc,squeeze(GRID.mskc(:,:,t,:)));
  tc2(:,:,t,:)=interp3ddata_int2(xl,yl,zlc,TC,xc(:,:,t),yc(:,:,t),zc,squeeze(GRID.mskc(:,:,t,:)));
 end
end
check3dmask(tc);
check3dmask(tc2);
save tc tc tc2
%wrda('temp.bin',tc2,1,fmt,Ieee);

disp('Reading Levitus climatological S');
[SC,xl,yl,zlc]=...
  extract_levitus_clim('salt',lon_lo-dd,lon_hi+dd,lat_lo-dd,lat_hi+dd);
if size(xc,3)==1
  sc=interp3ddata(xl,yl,zlc,SC,xc,yc,zc,GRID.mskc);
  sc2=interp3ddata_int2(xl,yl,zlc,SC,xc,yc,zc,GRID.mskc);
else
 for t=1:size(xc,3);
  sc(:,:,t,:)=interp3ddata(xl,yl,zlc,SC,xc(:,:,t),yc(:,:,t),zc,squeeze(GRID.mskc(:,:,t,:)));
  sc2(:,:,t,:)=interp3ddata_int2(xl,yl,zlc,SC,xc(:,:,t),yc(:,:,t),zc,squeeze(GRID.mskc(:,:,t,:)));
 end
end
check3dmask(sc);
check3dmask(sc2);
save sc sc sc2
%wrda('salt.bin',sc2,1,fmt,Ieee);

stats( mask(tc(:))-mask(sc(:)) )
%clear tc sc tc2 sc2


% Levitus climatology grid
[LCGRID.xc,LCGRID.yc,LCGRID.xg,LCGRID.yg,LCGRID.rac]=grid_sphere( ...
  (xl(2)-xl(1))*ones(size(xl)),(yl(2)-yl(1))*ones(size(yl)), ...
  1.5*xl(1)-0.5*xl(2),1.5*yl(1)-0.5*yl(2));
LCGRID.mskc=mask(TC);
[TCm,TCsd] = meanprofile(LCGRID,TC);
[SCm,SCsd] = meanprofile(LCGRID,SC);
stats( mask(TC(:))-mask(SC(:)) )
clear TC SC

disp('Reading Levitus monthly T');
[TM,xl,yl,zlm]=...
  extract_levitus_monthly('temp',1:12,lon_lo-dd,lon_hi+dd,lat_lo-dd,lat_hi+dd);
kk=find(zc>=min(zlm));
for m=1:12;
 if size(xc,3)==1
  tm(:,:,:,m)=interp3ddata(xl,yl,zlm,TM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk));
  tm2(:,:,:,m)=interp3ddata_int2(xl,yl,zlm,TM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk));
 else
  for t=1:size(xc,3);
   tm(:,:,t,:,m)=interp3ddata(xl,yl,zlm,TM(:,:,:,m),xc(:,:,t),yc(:,:,t),zc(kk),squeeze(GRID.mskc(:,:,kk,t)));
   tm2(:,:,t,:,m)=interp3ddata_int2(xl,yl,zlm,TM(:,:,:,m),xc(:,:,t),yc(:,:,t),zc(kk),squeeze(GRID.mskc(:,:,kk,t)));
%tm(:,:,:,m)=interp3ddata(xl,yl,zlm,TM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk));
%tm2(:,:,:,m)=interp3ddata_int2(xl,yl,zlm,TM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk));
  end
 end
end
%wrda('sst.bin',tm2(:,:,1,:),1,fmt,Ieee);
save tm tm tm2

disp('Reading Levitus monthly S');
[SM,xl,yl,zlm]=...
  extract_levitus_monthly('salt',1:12,lon_lo-dd,lon_hi+dd,lat_lo-dd,lat_hi+dd);
for m=1:12;
 if size(xc,3)==1
  sm(:,:,:,m)=interp3ddata(xl,yl,zlm,SM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk));
  sm2(:,:,:,m)=interp3ddata_int2(xl,yl,zlm,SM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk));
 else
  for t=1:size(xc,3);
   sm(:,:,t,:,m)=interp3ddata(xl,yl,zlm,SM(:,:,:,m),xc(:,:,t),yc(:,:,t),zc(kk),squeeze(GRID.mskc(:,:,kk,t)));
   sm2(:,:,t,:,m)=interp3ddata_int2(xl,yl,zlm,SM(:,:,:,m),xc(:,:,t),yc(:,:,t),zc(kk),squeeze(GRID.mskc(:,:,kk,t)));
%sm(:,:,:,m)=interp3ddata(xl,yl,zlm,SM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk));
%sm2(:,:,:,m)=interp3ddata_int2(xl,yl,zlm,SM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk));
  end
 end
end
%wrda('sss.bin',sm2(:,:,1,:),1,fmt,Ieee);
save sm sm sm2

% Levitus monthly grid
LMGRID=LCGRID;
LMGRID.mskc=mask(TM(:,:,:,1));
MM=mean(TM,4);
[TMm,TMsd] = meanprofile(LMGRID,MM);
clear TM MM
MGRID=GRID;
MGRID.mskc=MGRID.mskc(:,:,1:size(tm,3));
[tmm,tmsd] = meanprofile(MGRID,mean(tm,4));
[tm2m,tm2sd] = meanprofile(MGRID,mean(tm2,4));

MM=mean(SM,4);
[SMm,SMsd] = meanprofile(LMGRID,MM);
clear SM MM
[smm,smsd] = meanprofile(MGRID,mean(sm,4));
[sm2m,sm2sd] = meanprofile(MGRID,mean(sm2,4));

km=size(tmm,2);
end

load tc
[tcm,tcsd] = meanprofile(GRID,tc);
[tc2m,tc2sd] = meanprofile(GRID,tc2);
load sc
[scm,scsd] = meanprofile(GRID,sc);
[sc2m,sc2sd] = meanprofile(GRID,sc2);

% $$$ figure(1)
% $$$ subplot(121)
% $$$ plot(TCm,zlc,'k',...
% $$$      TMm,zlm,'r',...
% $$$      tcm,zc,'g',...
% $$$      tc2m,zc,'b',...
% $$$      tmm,zc(1:km),'c',...
% $$$      tm2m,zc(1:km),'m',...
% $$$      [TCm'-TCsd' TCm'+TCsd'],zlc,'k--',...
% $$$      [TMm'-TMsd' TMm'+TMsd'],zlm,'r--',...
% $$$      [tcm'-tcsd' tcm'+tcsd'],zc,'g--',...
% $$$      [tc2m'-tc2sd' tc2m'+tc2sd'],zc,'b--')
% $$$ xlabel('\theta (^oC)');ylabel('Depth (m)');
% $$$ legend('Levitus Clim.','Levitus Monthly','Interpolated','Integrated',4)
% $$$ 
% $$$ subplot(122)
% $$$ plot(SCm,zlc,'k',...
% $$$      SMm,zlm,'r',...
% $$$      scm,zc,'g',...
% $$$      sc2m,zc,'b',...
% $$$      smm,zc(1:km),'c',...
% $$$      sm2m,zc(1:km),'m',...
% $$$      [SCm'-SCsd' SCm'+SCsd'],zlc,'k--',...
% $$$      [SMm'-SMsd' SMm'+SMsd'],zlm,'r--',...
% $$$      [scm'-scsd' scm'+scsd'],zc,'g--',...
% $$$      [sc2m'-sc2sd' sc2m'+sc2sd'],zc,'b--')
% $$$ xlabel('S (psu)');ylabel('Depth (m)');
% $$$ legend('Levitus Clim.','Levitus Monthly','Interpolated','Integrated',4)
% $$$ 
% $$$ supertext(.5,1,...
% $$$  'Levitus climatology, monthly and interpolated datasets',...
% $$$    'FontSize',14,...
% $$$    'VerticalAlignment','Bottom',...
% $$$    'HorizontalAlignment','Center')
% $$$ print -depsc2 comparison_profiles.eps
% $$$ print -djpeg100 -r90 comparison_profiles.jpg
% $$$ 
% $$$ figure(2)
% $$$ subplot(121)
% $$$ plot(TMm-TCm(1:19),zlm,'k',tc2m-tcm,zc,'r')
% $$$ xlabel('\Delta \theta (^oC)');ylabel('Depth (m)');
% $$$ legend('Clim. - Monthly','Integr. - Intep.',4);
% $$$ subplot(122)
% $$$ plot(SMm-SCm(1:19),zlm,'k',sc2m-scm,zc,'r')
% $$$ xlabel('\Delta S (psu)');ylabel('Depth (m)');
% $$$ legend('Clim. - Monthly','Integr. - Intep.',4);
% $$$ 
% $$$ supertext(.5,1,...
% $$$  'Difference between datasets',...
% $$$    'FontSize',14,...
% $$$    'VerticalAlignment','Bottom',...
% $$$    'HorizontalAlignment','Center')
% $$$ print -depsc2 levitus_problems.eps
% $$$ print -djpeg100 -r90 levitus_problems.jpg
% $$$ 
% $$$ 
% $$$ % Test EOS
% $$$ P3=ini_poly3;
% $$$ load tc tc2;T=tc2;clear tc2
% $$$ %T=rdda('temp.bin',[90 40 15],1,fmt,Ieee);
% $$$ load sc sc2;S=sc2;clear sc2
% $$$ %S=rdda('salt.bin',[90 40 15],1,fmt,Ieee);
% $$$ D=dens_poly3(P3,T,S);
% $$$ 
% $$$ z=ones(90*40,1)*zc;
% $$$ y=ones(90,1)*yc;y=y(:);
% $$$ p=sw_pres(-z',y')';
% $$$ p=reshape(p,[90 40 15]);
% $$$ clear z y
% $$$ d=(sw_dens(S,T,p)-1000).*mask(T);
% $$$ Dd=D-d;
% $$$ 
% $$$ [Dm,Dsd]=meanprofile(GRID,D);
% $$$ [dm,dsd]=meanprofile(GRID,d);
% $$$ [Ddm,Ddsd]=meanprofile(GRID,D,d);
% $$$ 
% $$$ figure(3)
% $$$ subplot(121)
% $$$ plot(dm,zc,'k',...
% $$$      Dm,zc,'g',...
% $$$      [dm'-dsd' dm'+dsd'],zc,'k--',...
% $$$      [Dm'-Dsd' Dm'+Dsd'],zc,'g--')
% $$$ axis([22 51 -5000 0])
% $$$ xlabel('\rho-1000 (kg/m^3)');ylabel('Depth (m)');
% $$$ legend('Seawater EOS','POLY3',3)
% $$$ 
% $$$ subplot(122)
% $$$ plot(Ddm,zc,'k',...
% $$$      [Ddm'-dsd' Ddm'+dsd'],zc,'g',...
% $$$      [Ddm'-Ddsd' Ddm'+Ddsd'],zc,'k--')
% $$$ axis([[-1 1]*.5 -5000 0])
% $$$ xlabel('\Delta \rho (kg/m^3)');ylabel('Depth (m)');
% $$$ title('POLY3-EOS80')
% $$$ legend('Mean difference','Observed SD of \rho','SD of P3-EOS80',4)
% $$$ print -depsc2 eos_problems.eps
% $$$ print -djpeg100 -r90 eos_problems.jpg
% $$$ 
% $$$ figure(4)
% $$$ subplot(221)
% $$$ pcol(xc,yc, sq(D(:,:,8))' );
% $$$ ylabel('Latitude (^oN)');
% $$$ colorbar('h');
% $$$ title('\rho @ z=-1250m');
% $$$ 
% $$$ subplot(222)
% $$$ pcol(xc,yc, sq(Dd(:,:,8))' );
% $$$ ylabel('Latitude (^oN)');
% $$$ %cax=caxis;caxis([-1 1]*max(abs(cax)));
% $$$ colorbar('h');
% $$$ title('\Delta \rho (POLY3-EOS80) @ z=-1250m');
% $$$ 
% $$$ subplot(223)
% $$$ pcol(xc,yc, sq(D(:,:,12))' );
% $$$ ylabel('Latitude (^oN)');
% $$$ colorbar('h');
% $$$ title('\rho @ z=-3010m');
% $$$ 
% $$$ subplot(224)
% $$$ pcol(xc,yc, sq(Dd(:,:,12))' );
% $$$ ylabel('Latitude (^oN)');
% $$$ %cax=caxis;caxis([-1 1]*max(abs(cax)));
% $$$ colorbar('h');
% $$$ title('\Delta \rho (POLY3-ESO80) @ z=-3010m');
% $$$ print -depsc2 eos_problems2.eps
% $$$ print -djpeg100 -r90 eos_problems2.jpg
% $$$ 
% $$$ figure(5)
% $$$ subplot(221)
% $$$ pcol(xc,yc, sq(mean(tm2(:,:,1,:)-tm(:,:,1,:),4))')
% $$$ ylabel('Latitude (^oN)');
% $$$ caxis([-1 1]*max(abs(caxis)));colorbar('h')
% $$$ title('\Delta\theta  (Integrated-Interpolated) k=1')
% $$$ 
% $$$ subplot(222)
% $$$ pcol(xc,yc, sq(mean(sm2(:,:,1,:)-sm(:,:,1,:),4))')
% $$$ ylabel('Latitude (^oN)');
% $$$ caxis([-1 1]*0.25)
% $$$ caxis([-1 1]*max(abs(caxis)));colorbar('h')
% $$$ title('{\Delta}S  (Integrated-Interpolated) k=1')
% $$$ 
% $$$ subplot(223)
% $$$ pcol(xc,yc, sq(mean(tc2(:,:,10,:)-tc(:,:,10,:),4))')
% $$$ ylabel('Latitude (^oN)');
% $$$ caxis([-1 1]*max(abs(caxis)));colorbar('h')
% $$$ title('\Delta\theta  (Integrated-Interpolated) k=10')
% $$$ 
% $$$ subplot(224)
% $$$ pcol(xc,yc, sq(mean(sc2(:,:,10,:)-sc(:,:,10,:),4))')
% $$$ ylabel('Latitude (^oN)');
% $$$ caxis([-1 1]*max(abs(caxis)));colorbar('h')
% $$$ title('{\Delta}S  (Integrated-Interpolated) k=10')
% $$$ 
% $$$ print -depsc2 intgr-interp.eps
% $$$ print -djpeg100 -r90 intgr-interp.jpg