%computing MLD-10m up to ILD , balterms integration.
%nohup matlab -nodisplay -nosplash -nodesktop -r "run('find_integration_ild40.m')" > nohup_ild40_t1.txt 2>&1 &

lati=ncread('/mnt/lustre01/work/uo0122/u241089/ATL/ATL24/grid_new.nc','yc'); %2D
long=ncread('/mnt/lustre01/work/uo0122/u241089/ATL/ATL24/grid_new.nc','xc'); %2D
depth=ncread('/mnt/lustre01/work/uo0122/u241089/ATL/ATL24/grid_new.nc','depth_t',[2],[38],[1]);  %1D
dzo=ncread('/mnt/lustre01/work/uo0122/u241089/ATL/ATL24/grid_new.nc','drf',[2],[38],[1]); %2D

dzor=repmat(dzo,[1 755 2698]);dzor=permute(dzor,[3 2 1]); %depth spacing


% 1 - Term 1 (Hor. Adv.) | lon 2:2699 | lat 662:1416
% 2 - Term 2             | lon 1:2700 | lat 661:1417
% 3 - Term 3 (Tilting)   | lon 2:2699 | lat 662:1416
% 4 - Term 4 (Stretching)| lon 1:2700 | lat 661:1417
% 5 - Term 5 (turbmix)      | lon 1:2700 | lat 661:1417 
% 6 - Term LHS           | lon 1:2700 | lat 661:1417
% 7 - ILD / MLD / BLT    | lon 1:2700 | lat 517:1416
lat_start = [662 661 662 661 661 661 662 662 517];
lat_end   = [1416 1417 1416 1417 1417 1417 1416 1416 1416];
lon_start = [2 1 2 1 1 1 2 2 1];
lon_end   = [2699 2700 2699 2700 2700 2700 2699 2699 2700];
lat_range = [662 1416];
lon_range = [2 2699];
%lat_range = [latlim1 latlim2];
%lon_range = [lonlim1 lonlim2];
Term_Labels = ["Term 1 (Horizontal Advection)", "Term 2", "Term 3 (Tilting)", "Term 4 (Stretching)", "Term 5 (Turbulent Mixing)","LHS - sum term1to4","term1 + term2", "LHS"];

% Path = Folder/Filename_%4d<timestep>_DD-MM-YYYY.nc
Term_Filenames=["RHS_ATL24_horadv_term1", "RHS_ATL24_term2", "RHS_ATL24_tilt_term3", "RHS_ATL24_stretch_term4", "RHS_ATL24_term5_turbmix", "LHS_ATL24", "residual5_sum1to4", "RHS_term1n2sum"];
Term_Folders=["HorAdv", "term2", "Tilting", "Stretching", "term5_turbmix", "LHS", "Residual_term5", "term1n2sum"];
Term_Variables_T=["horadvT", "wTzz", "tiltingT", "stretchingT", "turbmixT", "dTzt", "residualT_term5", "sum12T"];
Term_Variables_S=["horadvS", "wSzz", "tiltingS", "stretchingS", "turbmixS", "dSzt", "residualS_term5", "sum12S"];


% Read lat from Term:
%  -> lat_j - lat_start + 1 =
% Read lon from Term:
%  -> lon_i - lon_start + 1 =
% Plot for range:
%  lower bound = lat_range(1) - lat_start(t) + 1 = 
%  upper bound = lower_bound + lat_range(2) - lat_range(1)
%  upper bound -> lat_range(2) - lat_range(1) + lat_range(1) - lat_start(t) + 1 = lat_range(2) + 1 - lat_start(t)
% -> lat_range(1) - lat_start(t) + 1 : lat_range(2) + 1 - lat_start(t)
%  lower bound -> lon_range(1) - lon_start(t) + 1 = 
%  upper bound = lower_bound + lon_range(2) - lon_range(1)
%  upper bound -> lon_range(2) - lon_range(1) + lon_range(1) - lon_start(t) + 1 = lon_range(2) + 1 - lon_start(t)
% -> lon_range(1) - lon_start(t) + 1 : lon_range(2) + 1 - lon_start(t)


%Dates for output (plot) filenames
DatesOut=strings(3287,1);
idate=1;
MonthDays=[31 28 31 30 31 30 31 31 30 31 30 31];
MonthDaysLeap=[31 29 31 30 31 30 31 31 30 31 30 31];
for iyear=2003:2011
  for imonth=1:12
    if iyear==2004 || iyear==2008
      for iday=1:MonthDaysLeap(imonth)
        DatesOut(idate)=sprintf('%02d-%02d-%d',iday, imonth, iyear);
        idate=idate+1;
      end
    else
      for iday=1:MonthDays(imonth)
        DatesOut(idate)=sprintf('%02d-%02d-%d',iday, imonth, iyear);
        idate=idate+1;
      end
    end
  end
end


    
for t=6 %1:8 %ERMX %1:6   %terms
%tic
t
  

for l=1892 %:1656 %VALUEX %1:3286
tic
l

filenum=sprintf('%04d', l);
FileName=char(strcat('/work/uo0122/u241272/ATL24/Sal_temp_balance_terms/',Term_Folders(t),'/',Term_Filenames(t),'_',filenum,'_',DatesOut(l), '.nc'));

T=squeeze(ncread(FileName,char(Term_Variables_T(t)),[(lon_range(1)-lon_start(t)+1) (lat_range(1)-lat_start(t)+1) 1 1],[(lon_range(2)-lon_range(1)+1) (lat_range(2)-lat_range(1)+1) 38 1],[1 1 1 1]));
S=squeeze(ncread(FileName,char(Term_Variables_S(t)),[(lon_range(1)-lon_start(t)+1) (lat_range(1)-lat_start(t)+1) 1 1],[(lon_range(2)-lon_range(1)+1) (lat_range(2)-lat_range(1)+1) 38 1],[1 1 1 1]));

ild=ncread('/mnt/lustre01/work/uo0122/u241272/ATL24/ATL24_lt_leap.nc','ILD',[lon_range(1)-lon_start(9)+1 lat_range(1)-lat_start(9)+1 l],[lon_range(2)-lon_range(1)+1 lat_range(2)-lat_range(1)+1 1],[1 1 1]);
mld=ncread('/mnt/lustre01/work/uo0122/u241272/ATL24/ATL24_lt_leap.nc','MLD',[lon_range(1)-lon_start(9)+1 lat_range(1)-lat_start(9)+1 l],[lon_range(2)-lon_range(1)+1 lat_range(2)-lat_range(1)+1 1],[1 1 1]);


T_10mld_ild=zeros(size(ild,1),size(ild,2));
S_10mld_ild=zeros(size(ild,1),size(ild,2));

Td=T.*dzor;
Sd=S.*dzor;


for i=1:size(ild,1)
  for j=1:size(ild,2)

    indm=find(depth==mld(i,j));
    ind=find(depth==ild(i,j));
    
    if isempty(ind) || isempty(indm)
    T_10mld_ild(i,j)=NaN;
    S_10mld_ild(i,j)=NaN;
    elseif indm<=2
    T_10mld_ild(i,j)=nansum(Td(i,j,1:ind),3);
    S_10mld_ild(i,j)=nansum(Sd(i,j,1:ind),3);
    else
    T_10mld_ild(i,j)=nansum(Td(i,j,indm-2:ind),3);
    S_10mld_ild(i,j)=nansum(Sd(i,j,indm-2:ind),3);
    end  

    clear ind indm %ind2
  end %j
end %i

      
mask=S(:,:,1);
mask(~isnan(mask))=1;
T_10mld_ild=T_10mld_ild.*mask;
S_10mld_ild=S_10mld_ild.*mask;


%T_10mld_ild=single(T_10mld_ild);
%S_10mld_ild=single(S_10mld_ild);


outfile=char(strcat('/work/uo0122/u241272/ATL24/Sal_temp_balance_terms/',Term_Folders(t),'/Integrated_10mld_ild/',Term_Filenames(t),'_10mld_ild_',filenum,'_',DatesOut(l),'.nc'));

nccreate(outfile,'T_10mld_ild','dimensions',{'longitude' lon_end(1)-lon_start(1)+1 'latitude' lat_end(1)-lat_start(1)+1 'time' 1},'Datatype','single','ChunkSize',[270 76 1],'DeflateLevel',2);
ncwriteatt(outfile,'T_10mld_ild','units','psu/sec');
ncwriteatt(outfile,'T_10mld_ild','description','integrated mld-10m up to ild');

nccreate(outfile,'S_10mld_ild','dimensions',{'longitude' lon_end(1)-lon_start(1)+1 'latitude' lat_end(1)-lat_start(1)+1 'time' 1},'Datatype','single','ChunkSize',[270 76 1],'DeflateLevel',2);
ncwriteatt(outfile,'S_10mld_ild','units','psu/sec');
ncwriteatt(outfile,'S_10mld_ild','description','integrated mld-10m up to ild');

nccreate(outfile,'time','dimensions',{'time' 1},'Datatype','single');
ncwriteatt(outfile,'time','units','days since 2002-12-31 00:00');
ncwriteatt(outfile,'time','calendar','gregorian');

ncwrite(outfile,'T_10mld_ild',T_10mld_ild);
ncwrite(outfile,'S_10mld_ild',S_10mld_ild);
ncwrite(outfile,'time',l);


clear T_10mld_ild S_10mld_ild T S Td Sd ild 

toc
end %l


%toc
end %t

exit
