

addpath /work/uo0122/u241272/seawater_ver3_3.1/

dep_atl=ncread('/mnt/lustre01/work/uo0122/u241089/ATL/ATL24/grid_new.nc','depth_t',1,40,1);

ILD=zeros(360,173,108);
MLD=zeros(360,173,108);
BLT=zeros(360,173,108);



%Compute ILD

tempcEN_atldep=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/interp_EN4_2_1_TS.nc','temp',[1 1 1 1],[360 173 40 108],[1 1 1 1]);

for k=1:40
aux(:,:,k,:)=tempcEN_atldep(:,:,1,:) - 0.2;
end

tem2=tempcEN_atldep-aux;

for l=1:108
tic
l

 for i=1:360
  for j=1:173
    %ind=find(squeeze(tem2(i,j,:,l))>-0.2 & squeeze(tem2(i,j,:,l))<0.2);
    ind=find(squeeze(tem2(i,j,:,l))<=0);
    if isempty(ind)
       if isnan(tempcEN_atldep(i,j,1,l))
       ILD(i,j,l)=NaN;
       else
       tm=find(isnan(tempcEN_atldep(i,j,:,l)));
           if isempty(tm)
           ILD(i,j,l)=dep_atl(40);
           else
           ILD(i,j,l)=dep_atl(tm(1)-1);
           end
       end
    else
    ILD(i,j,l)=dep_atl(ind(1));
    end
    clear ind
  end %j
 end %i
 
toc
end %l


%Compute MLD

saltEN_atldep=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/interp_EN4_2_1_TS.nc','salt',[1 1 1 1],[360 173 40 108],[1 1 1 1]);


for l=1:108
for k=1:40

pden(:,:,k,l)=sw_dens(squeeze(saltEN_atldep(:,:,k,l)),squeeze(tempcEN_atldep(:,:,k,l)),0);

end %k
end %l


for l=1:108

corden(:,:,1,l)=sw_dens(squeeze(saltEN_atldep(:,:,1,l)),squeeze(tempcEN_atldep(:,:,1,l)-0.2),0);

end %l

clear saltEN_atldep

for k=1:40
auxm(:,:,k,:)=corden(:,:,1,:);
end

sigma=pden-auxm;

for l=1:108
tic
l

 for i=1:360
  for j=1:173
    %ind=find(squeeze(tem2(i,j,:,l))>-0.2 & squeeze(tem2(i,j,:,l))<0.2);
    ind=find(squeeze(sigma(i,j,:,l))>=0);
    if isempty(ind)
        if isnan(tempcEN_atldep(i,j,1,l))
       MLD(i,j,l)=NaN;
       else
       dm=find(isnan(pden(i,j,:,l)));
           if isempty(dm)
           MLD(i,j,l)=dep_atl(40);
           else
           MLD(i,j,l)=dep_atl(dm(1)-1);
           end
       end
    %MLD(i,j)=NaN;
    else
    MLD(i,j,l)=dep_atl(ind(1));
    end
    clear ind
  end %j
 end %i
 
toc
end %l


%Compute BLT

BLT = ILD - MLD; 


BLT_grid=BLT;
BLT_grid(361:360+50,:,:)=BLT(1:50,:,:);
BLT_grid(1:50,:,:)=[];

ILD_grid=ILD;
ILD_grid(361:360+50,:,:)=ILD(1:50,:,:);
ILD_grid(1:50,:,:)=[];

MLD_grid=MLD;
MLD_grid(361:360+50,:,:)=MLD(1:50,:,:);
MLD_grid(1:50,:,:)=[];


x_new=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/interp_EN4_2_1_TS.nc','longitude');
y=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/interp_EN4_2_1_TS.nc','latitude');

%saving


nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','longitude','dimensions',{'longitude' 360 'latitude' 173});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','latitude','dimensions',{'longitude' 360 'latitude' 173});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','BLT','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','BLT_grid','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','MLD_grid','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','MLD','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','ILD','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','ILD_grid','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','longitude',x_new);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','latitude',y);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','BLT',BLT);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','BLT_grid',BLT_grid);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','MLD',MLD);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','ILD',ILD);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','ILD_grid',ILD_grid);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','MLD_grid',MLD_grid);






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%{
%load newfor_interp_dep_TS_en4.mat dep_ENnew
load grid_EN4.mat lat_EN lon_EN
load /work/uo0122/u241272/Compare_EN4_GECCO2/interp_EN4_T_atldep.mat %tempcEN_atldep
load /work/uo0122/u241272/Compare_EN4_GECCO2/interp_EN4_S_atldep.mat %saltEN_atldep
addpath /work/uo0122/u241272/seawater_ver3_3.1/

dep_atl=ncread('/mnt/lustre01/work/uo0122/u241089/ATL/ATL24/grid_new.nc','depth_t',1,40,1);

ILD=zeros(360,173,108);
MLD=zeros(360,173,108);
BLT=zeros(360,173,108);

%Compute ILD

for k=1:40
aux(:,:,k,:)=tempcEN_atldep(:,:,1,:) - 0.2;
end

tem2=tempcEN_atldep(:,:,1:40,:)-aux;

for l=1:108
tic
l

 for i=1:360
  for j=1:173
    %ind=find(squeeze(tem2(i,j,:,l))>-0.2 & squeeze(tem2(i,j,:,l))<0.2);
    ind=find(squeeze(tem2(i,j,:,l))<=0);
    if isempty(ind)
       if isnan(tempcEN_atldep(i,j,1,l))
       ILD(i,j,l)=NaN;
       else
       tm=find(isnan(tempcEN_atldep(i,j,1:40,l)));
           if isempty(tm)
           ILD(i,j,l)=dep_atl(40);
           else
           ILD(i,j,l)=dep_atl(tm(1)-1);
           end
       end
    else
    ILD(i,j,l)=dep_atl(ind(1));
    end
    clear ind
  end %j
 end %i
 
toc
end %l


%Compute MLD

for l=1:108
for k=1:40

pden(:,:,k,l)=sw_dens(squeeze(saltEN_atldep(:,:,k,l)),squeeze(tempcEN_atldep(:,:,k,l)),0);

end %k
end %l


for l=1:108

corden(:,:,1,l)=sw_dens(squeeze(saltEN_atldep(:,:,1,l)),squeeze(tempcEN_atldep(:,:,1,l)-0.2),0);

end %l

clear saltEN_atldep

for k=1:40
auxm(:,:,k,:)=corden(:,:,1,:);
end

sigma=pden-auxm;

for l=1:108
tic
l

 for i=1:360
  for j=1:173
    %ind=find(squeeze(tem2(i,j,:,l))>-0.2 & squeeze(tem2(i,j,:,l))<0.2);
    ind=find(squeeze(sigma(i,j,:,l))>=0);
    if isempty(ind)
        if isnan(tempcEN_atldep(i,j,1,l))
       MLD(i,j,l)=NaN;
       else
       dm=find(isnan(pden(i,j,:,l)));
           if isempty(dm)
           MLD(i,j,l)=dep_atl(40);
           else
           MLD(i,j,l)=dep_atl(dm(1)-1);
           end
       end
    %MLD(i,j)=NaN;
    else
    MLD(i,j,l)=dep_atl(ind(1));
    end
    clear ind
  end %j
 end %i
 
toc
end %l


%Compute BLT

BLT = ILD - MLD; 

new_lon_EN=lon_EN;
new_lon_EN(new_lon_EN>180)=new_lon_EN(new_lon_EN>180)-360;
[y,x]=meshgrid(lat_EN,new_lon_EN);

BLT_grid=BLT;
BLT_grid(361:360+50,:,:)=BLT(1:50,:,:);
BLT_grid(1:50,:,:)=[];

ILD_grid=ILD;
ILD_grid(361:360+50,:,:)=ILD(1:50,:,:);
ILD_grid(1:50,:,:)=[];

MLD_grid=MLD;
MLD_grid(361:360+50,:,:)=MLD(1:50,:,:);
MLD_grid(1:50,:,:)=[];

x_new=x;
x_new(361:360+50,:)=x(1:50,:);
x_new(1:50,:)=[];


%saving


nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','longitude','dimensions',{'longitude' 360 'latitude' 173});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','latitude','dimensions',{'longitude' 360 'latitude' 173});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','BLT','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','BLT_grid','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','MLD_grid','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','MLD','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','ILD','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
nccreate('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','ILD_grid','dimensions',{'longitude' 360 'latitude' 173 'time' 108});
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','longitude',x_new);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','latitude',y);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','BLT',BLT);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','BLT_grid',BLT_grid);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','MLD',MLD);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','ILD',ILD);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','ILD_grid',ILD_grid);
ncwrite('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4toATLdep.nc','MLD_grid',MLD_grid);


save('blt_ild02_mld_interpEN4toATLdep.mat','BLT_grid','-v7.3'); 
save('blt_ild02_mld_interpEN4toATLdep.mat','BLT','-append');
save('blt_ild02_mld_interpEN4toATLdep.mat','ILD_grid','-append');
save('blt_ild02_mld_interpEN4toATLdep.mat','ILD','-append');
save('blt_ild02_mld_interpEN4toATLdep.mat','MLD_grid','-append');
save('blt_ild02_mld_interpEN4toATLdep.mat','MLD','-append');
save('blt_ild02_mld_interpEN4toATLdep.mat','x_new','-append');
save('blt_ild02_mld_interpEN4toATLdep.mat','y','-append');
%}

