clear all
close all
clc

new=pwd;
old='/phoddat/share/pita/EN4/';
lat=ncread('/phoddat/share/pita/EN4/2005/EN.4.2.2.f.analysis.g10.200501.nc','lat');
lon=ncread('/phoddat/share/pita/EN4/2005/EN.4.2.2.f.analysis.g10.200501.nc','lon');
depth=ncread('/phoddat/share/pita/EN4/2005/EN.4.2.2.f.analysis.g10.200501.nc','depth');
depth_bnds=ncread('/phoddat/share/pita/EN4/2005/EN.4.2.2.f.analysis.g10.200501.nc','depth_bnds');

%% load EN4 data

year=2000:2023;
T=nan(length(lon),length(lat(lat<=0)),length(depth),length(year).*12);
S=T;
cont=0;

for yy=year(1):year(end)
    
    for mm=1:12
        display([num2str(mm) '/' num2str(yy)])
        if mm<10
            file=[old num2str(yy) '/EN.4.2.2.f.analysis.g10.' num2str(yy) '0' num2str(mm) '.nc' ];
        else
            file=[old num2str(yy) '/EN.4.2.2.f.analysis.g10.' num2str(yy) num2str(mm) '.nc' ];
        end
        cont=cont+1;
        T_file=ncread(file,'temperature');
T(:,:,:,cont)=T_file(:,lat<=0,:);
        S_file=ncread(file,'salinity');    
     S(:,:,:,cont)=S_file(:,lat<=0,:);
    time(cont)=julian(yy,mm,15,0);

    end
end
lat=lat(lat<=0);

%organize longitude data
lon(lon>=180)=lon(lon>=180)-360;

T=cat(1,T(lon<=0,:,:,:),T(lon>0,:,:,:));
S=cat(1,S(lon<=0,:,:,:),S(lon>0,:,:,:));
lon=cat(1,lon(lon<=0),lon(lon>0));
ind_lon=lon>=-55 & lon<=20;
ind_lat=lat<=0& lat>=-40;

T=T(ind_lon,ind_lat,:,:);
S=S(ind_lon,ind_lat,:,:);
lon=lon(ind_lon);
lat=lat(ind_lat);

save('/phoddat/share/pita/EN4/EN4_T_S_SA.mat','time','lat','lon','depth','depth_bnds','T','S');
