%time average of all terms 9 yrs.
%nohup matlab -nodisplay -nosplash -nodesktop -r "run('monthly_mean_balterms.m')" > nohup_mm_t3.txt 2>&1 &

% 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 (Flux)      | 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 517];
lat_end   = [1416 1417 1416 1417 1417 1417 1416];
lon_start = [2 1 2 1 1 1 1];
lon_end   = [2699 2700 2699 2700 2700 2700 2700];
lat_range = [662 1416];
lon_range = [2 2699];
Term_Labels = ["Term 1 (Horizontal Advection)", "Term 2", "Term 3 (Tilting)", "Term 4 (Stretching)", "Term 5 (Turbulent Mixing)", "LHS"];

% 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)

% 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_turbmix_term5", "LHS_ATL24", "RHS_ATL24_wSp_top_epr"];
Term_Folders=["HorAdv", "term2", "Tilting", "Stretching", "term5_turbmix", "LHS", "term5_timemean_top_epr"];
Term_Variables_T=["horadvT", "wTzz", "tiltingT", "stretchingT", "turbmixT", "dTzt","wT"];
Term_Variables_S=["horadvS", "wSzz", "tiltingS", "stretchingS", "turbmixS", "dSzt","wS"];

%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

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

Months_364 = [31 28 31 30 31 30 31 31 30 31 30 30]';
Months_365 = [31 28 31 30 31 30 31 31 30 31 30 31]';
Months_366 = [31 29 31 30 31 30 31 31 30 31 30 31]';
Months = [Months_365, Months_366, Months_365, Months_365, Months_365, Months_366, Months_365, Months_365, Months_364];

for t=TERMX
 
  Term_Filenames(t)

  for year=YEARX %2003:2011
    YearString=sprintf('%04d', year);

    for i=MONTHX %  1:12      
 
      MonthString=sprintf('%02d', i);
      outfile=char(strcat('/work/uo0122/u241272/ATL24/Sal_temp_balance_terms/',Term_Folders(t),'/',Term_Filenames(t),'_',MonthString,'-',YearString,'_mm', '.nc'));  
      if isfile(outfile)
  	outfile
 	continue
      end

      tic
      i

      mmS=zeros(lon_range(2)-lon_range(1)+1,lat_range(2)-lat_range(1)+1,46,1); %lon_range(2)-lon_range(1)+1,lat_range(2)-lat_range(1)+1
      mmT=zeros(lon_range(2)-lon_range(1)+1,lat_range(2)-lat_range(1)+1,46,1); %lon_range(2)-lon_range(1)+1,lat_range(2)-lat_range(1)+1
      
      for j=1:Months(i, year-2002)
        DateIndex=strcat(sprintf('%02d', j),'-',sprintf('%02d', i),'-',sprintf('%04d', year))
	index=find(DatesOut==DateIndex);
        filenum=sprintf('%04d', index)

        FileName=char(strcat('/work/uo0122/u241272/ATL24/Sal_temp_balance_terms/',Term_Folders(t),'/',Term_Filenames(t),'_',filenum,'_',DatesOut(index), '.nc'));            

        DataT=ncread(FileName,char(Term_Variables_T(t)),[(lon_range(1)-lon_start(2)+1) (lat_range(1)-lat_start(2)+1) 1 1],[(lon_range(2)-lon_range(1)+1) (lat_range(2)-lat_range(1)+1) 46 1],[1 1 1 1]);
        DataS=ncread(FileName,char(Term_Variables_S(t)),[(lon_range(1)-lon_start(2)+1) (lat_range(1)-lat_start(2)+1) 1 1],[(lon_range(2)-lon_range(1)+1) (lat_range(2)-lat_range(1)+1) 46 1],[1 1 1 1]);
            
        mmS(:,:,:,1)=mmS(:,:,:,1)+DataS; 
        mmT(:,:,:,1)=mmT(:,:,:,1)+DataT;

        clear DataS DataT

      end %j day

     
      mmS(:,:,:,1)=mmS(:,:,:,1)./Months(i, year-2002);
      mmT(:,:,:,1)=mmT(:,:,:,1)./Months(i, year-2002);

      nccreate(outfile,'mmS','dimensions',{'longitude' 2698 'latitude' 755 'depth' 46 'time' 1},'Datatype','single','ChunkSize',[270 76 4 1],'DeflateLevel',2);
      ncwriteatt(outfile,'mmS','units','psu/m.sec');
      ncwriteatt(outfile,'mmS','description','monthly mean of each term of RHS z derivative of Sal balance');

      nccreate(outfile,'mmT','dimensions',{'longitude' 2698 'latitude' 755 'depth' 46 'time' 1},'Datatype','single','ChunkSize',[270 76 4 1],'DeflateLevel',2);
      ncwriteatt(outfile,'mmT','units','deg.C/m.sec');
      ncwriteatt(outfile,'mmT','description','monthly mean of each term of RHS z derivative of Temp balance');

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

      % writing the output to .nc file
      ncwrite(outfile,'mmS',mmS);
      ncwrite(outfile,'mmT',mmT);
      ncwrite(outfile,'time',index-Months(i, year-2002)/2);

      clear mmS mmT
    toc
    end %i month
  end %year

end %t  

exit
