
%This script is created by Marlos
%Last change: 01/07/2013
 addpath('/data/mgoes/matlab/Seawater')
 addpath('/data/mgoes/AX8/Time_transform')
addpath ../.

clear
close all

dias = textread('bad.list','%s\n');
nz0 = 140; %3000;
dvec = 1:140; %[1:2:1000 1050:50:nz0];
nz=length(dvec);%3000;
T = nan*ones(nz,300,length(dias));
S=T;P=T;DH=T;vg=T;depth=T;
LAT = nan*ones(300,length(dias));
LON=LAT;time=LAT;year=LAT;

for i=1:length(dias)
    month = dias{i}; cruise = 'ax07';

load(['Data/data_',month,'.mat'],'data_full0','lat','lon','TIME','variables')
nx = length(lon)
P(1:nz,1:nx,i)     = data_full0(dvec,1,:);
T(1:nz,1:nx,i)     = data_full0(dvec,2,:);
S(1:nz,1:nx,i)     = data_full0(dvec,4,:);
DH(1:nz,1:nx,i)    = data_full0(dvec,5,:);
vg(1:nz,1:nx,i)    = data_full0(dvec,7,:);
depth(1:nz,1:nx,i) = data_full0(dvec,6,:);
LAT(1:nx,i) = lat;
LON(1:nx,i) = lon;
time(1:nx,i) = TIME;
year(1:nx,i) = month2yfraction(num2str(TIME'));
end
year(year==0)=nan;

%SORT(TIME)
[month,imon]=sort(nanmean(year));
year = year(:,imon);
P    = P(:,:,imon);
T    = T(:,:,imon);
S    = S(:,:,imon);
DH   = DH(:,:,imon);
vg   = vg(:,:,imon);
depth= depth(:,:,imon);
LAT  = LAT(:,imon);
LON  = LON(:,imon);
time = time(:,imon);
dia = yfraction2month(month);
Pn=depth(:,1,1);

save(['Lat_lon_time_',cruise,'.mat'],'LON','LAT','time')

%% PLOT FIELDS


% OHC = squeeze(nanmean(T(1:21,:,:)));OHC(OHC==0)=nan;
% SC  = squeeze(nanmean(S(1:21,:,:)));SC(SC<=0)=nan;
% DH0 = squeeze(DH(1,:,:));DH0(DH0<=0)=nan;

LONS = LON;LONS(LONS==0)=nan;
LATS = LAT;LATS(LATS==0)=nan;
% DH0S = DH0;DH0S(LAT==0)=nan;
% SCS = SC;SCS(LAT==0)=nan;
% OHCS=OHC;OHCS(LAT==0)=nan;

TIME2       = time;
TIME2(TIME2==0)=nan;
    dy = nanmedian(LONS(2:end,:)-LONS(1:end-1,:));
    ilat = find(dy>0);  %NOW IT IS >0, SO REVERT FROM SOUTH TO NORTH
    LATS(:,ilat) = LATS(end:-1:1,ilat);
    LONS(:,ilat) = LONS(end:-1:1,ilat);
%     DH0S(:,ilat) = DH0S(end:-1:1,ilat);
%     SCS(:,ilat)  =  SCS(end:-1:1,ilat);
%     OHCS(:,ilat) = OHCS(end:-1:1,ilat);
    S(:,:,ilat) = S(:,end:-1:1,ilat);
    T(:,:,ilat) = T(:,end:-1:1,ilat);

    nx2 = nnz(~isnan(nanmean(TIME2,2)));
    nx = nnz(~isnan(nanmean(LATS,2)));
    ns = size(LATS,2);
    
    lats=nan*ones(nx,ns);lons=lats;
 %   dh0=lats;sc=lats;ohc=lats;
    S0 = nan*ones(nz,nx,ns);T0=S0;
    
    ina = ~isnan(TIME2);
    nna = ~isnan(LATS);
    
    lats(ina) = LATS(nna);
    lons(ina) = LONS(nna);
%     sc(ina)   = SCS(nna);
%     ohc(ina)  = OHCS(nna);
%     dh0(ina)  = DH0S(nna);
    S0(1:nz,ina) = S(1:nz,nna);
    T0(1:nz,ina) = T(1:nz,nna);
    
    LON = lons(1:nx2,:);LAT=lats(1:nx2,:);time = TIME2(1:nx2,:);
    S2 = S0(:,1:nx2,:);T2=T0(:,1:nx2,:);
    clear T0 S0
   save(['Lat_lon_time_T_S_',cruise,'.mat'],'LON','LAT','time','S2','T2','Pn') 
%%
stop
    lon_int = 0:-1:-80;nlo=length(lon_int);
    year_int=nanmean(year);
      for jj=1:ns
      lonj = lons(:,jj);
      scj   = sc(:,jj);
      ohcj = ohc(:,jj);
      dh0j = dh0(:,jj);
      ilon=~isnan(lonj);
      lonj = lonj(ilon);
      
      scj   = scj(ilon);
      ohcj = ohcj(ilon);
      dh0j = dh0j(ilon);
    
      [~,ilon]=unique(lonj);
        sci(1:nlo,jj)  = interp1(lonj(ilon),scj(ilon),lon_int');
       ohci(1:nlo,jj) = interp1(lonj(ilon),ohcj(ilon),lon_int');
       dh0i(1:nlo,jj) = interp1(lonj(ilon),dh0j(ilon),lon_int');
      end
     
DH0C = make_clim_harmonico(dh0i',dia,[],'monthly');
DH0A = dh0i-DH0C';
OHCC = make_clim_harmonico(ohci',dia,[],'monthly');
OHCA = ohci-OHCC';
SCC = make_clim_harmonico(sci',dia,[],'monthly');
SA = sci-SCC';

load mapa
 figure,pcolor(year_int,lon_int,SA);shading flat   
 colormap(mapa)
 caxis([-1 1])
 ylabel('Longitude')
 set(gca,'fontsize',12)
  title('SAL')

 figure,pcolor(year_int,lon_int,DH0A);shading flat   
 colormap(mapa)
 caxis([-40 40])
 ylabel('Longitude')
 set(gca,'fontsize',12)
 title('DH0')

 figure,pcolor(year_int,lon_int,OHCA);shading flat   
 colormap(mapa)
 caxis([-5 5])
 ylabel('Longitude')
 set(gca,'fontsize',12)
 title('T100')

 
 
 % load ../mapa
% figure,pcolor(year,LAT,OHCA);shading flat
% colormap(mapa)
% caxis([-10 10])
% 
% figure,pcolor(year,LAT,DH0A);shading flat
% colormap(mapa)
% caxis([-20 20])
% 
% figure,pcolor(year,LAT,SA);shading flat
% colormap(mapa)
% caxis([-.5 .5])
