% to find areal average of a region EN4 timemean
%nohup matlab -nodisplay -nosplash -nodesktop -r "run('make_area_avg_EN4_analysis.m')" > nohup_en4.txt 2>&1 &


% read lat lon and spacings

%lat=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','latitude');
%lon=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/blt_ild02_mld_interpEN4_2_1_toATLdep.nc','longitude');

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


% 173 lat y
%360 lon x

%{
%Area definition
lone=-46; 
lonw=-66;
latn=24;
lats=6;
%}

%{
lone=-46; 
lonw=-58;
latn=20;
lats=14;

  lone=-51; 
  lonw=-59;
  latn=10;
  lats=6;
%}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%other regions%%%%%%%%%east atlantic%%%%%%%%%%%%%%%%%%%%%%%%%%




%{
lone=15; 
lonw=-1;
latn=9;
lats=-1;
%}





%{
lone=15; 
lonw=7;
latn=-1;
lats=-11;


lone=15; 
lonw=-1;
latn=9;
lats=-12;

lone=15; 
lonw=-1;
latn=9;
lats=-11;

  Area_lone=[15   15   15   15   15    -1   -1];
  Area_lonw=[-1    7    7   -1   -1   -35  -35];
  Area_latn=[ 9   -1   -1    9    9    13   13];
  Area_lats=[-1  -12  -10  -12  -10    -9  -10];


  Area_lone=[15   15   15   15   15   15    -1   -1];
  Area_lonw=[ 7   -1    7    7   -1   -1   -35  -35];
  Area_latn=[-1    9   -1   -1    9    9    13   13];
  Area_lats=[-8   -1  -12  -10  -12  -10    -9  -10];
%}
  
  Area_lone=[-6     -6   15   15   15   15   15   15   15    -1   -1];
  Area_lonw=[-34    -30  -1    7   -1    7    7   -1   -1   -35  -35];
  Area_latn=[13     13    9   -1    9   -1   -1    9    9    13   13];
  Area_lats=[-4     -4   -8   -8   -1  -12  -10  -12  -10    -9  -10];


blt_abv0=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/EN4_BLT_interpEN4_2_1_toATLdep_monthly_climatology_200301-201112.nc','BLT');
blt_abv5=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/EN4_BLT_interpEN4_2_1_toATLdep_monthly_climatology_200301-201112.nc','BLT');
 %blt_abv0=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/EN4_BLT_reflev2_interpEN4_2_1_toATLdep_monthly_climatology_200301-201112.nc','BLT');

blt_abv0(blt_abv0<=0)=NaN;

blt_abv5(blt_abv5<=5)=NaN;



for areai=3 %:7 %:2
tic
areai

  lone=Area_lone(areai); 
  lonw=Area_lonw(areai);
  latn=Area_latn(areai);
  lats=Area_lats(areai);

  %Find grid values closest to latn,lats,lone,lonw
  %minlon=lone;
  %maxlon=lonw; 
  %minlat=latn;
  %maxlat=lats;

%{
%Find grid values closest to latn,lats,lone,lonw
minlat=latn;
maxlat=lats;
minlon=lone;
maxlon=lonw; 

%i,j of the grid values closest to latn,lats,lone,lonw
minlati=[0 0];
minloni=[0 0];
maxlati=[0 0];
maxloni=[0 0];
%}


%Count all the grid points within the area defined by latn,lats,lone,lonw
NumPoints=0;
for i=1:360
for j=1:173
if lat(i,j)>=lats && lat(i,j)<=latn && lon(i,j)>=lonw && lon(i,j)<=lone
NumPoints=NumPoints+1;
%i
%j
%lat(i,j)
%lon(i,j)
end
end
end

%Find i and j of all points in that area, find the values of the edges.
Ti=zeros(NumPoints,2); %Indices of all points in rectangle: i -> Ti(x,1); j-> Ti(x,2)
NumPoints=0;
for i=1:360
for j=1:173
if lat(i,j)>=lats && lat(i,j)<=latn && lon(i,j)>=lonw && lon(i,j)<=lone
%if lat(i,j)>=lats & lat(i,j)<=latn & lon(i,j)>=lonw & lon(i,j)<=lone
NumPoints=NumPoints+1;
Ti(NumPoints,1)=i;
Ti(NumPoints,2)=j;
%{
if lat(i,j)<minlat
minlat=lat(i,j);
minlati=[i j];
end
if lat(i,j)>maxlat
maxlat=lat(i,j);
maxlati=[i j];
end
if lon(i,j)<minlon
minlon=lon(i,j);
minloni=[i j];
end
if lon(i,j)>maxlon
maxlon=lon(i,j);
maxloni=[i j];
end
%}
end
end
end
%{
maxlon
maxloni
minlon
minloni
maxlat
maxlati
minlat
minlati
%}


% regional mean of blt above 10m.
%blt=ncread('ATL24_lt_leap_mnth.nc','BLT_mnth');
%blt_abv10=blt;
%blt_abv10(blt_abv10<=10)=NaN;
%mask=blt_abv10;
%mask=mask*0+1;

%maxblt_abv10m=zeros(108,1);
%Tavg_abv10m=zeros(108,1);
%Tsum=zeros(NumPoints,1);
%dsum=zeros(NumPoints,1);
%bltnan=zeros(NumPoints,1);

%for l=1:108
%  for ij=1:NumPoints
%    bltnan(ij)=blt_abv10(Ti(ij,1), Ti(ij,2), l);
%    if blt_abv10(Ti(ij,1), Ti(ij,2), l)>maxblt_abv10m(l)
%      maxblt_abv10m(l)=blt_abv10(Ti(ij,1), Ti(ij,2), l);
%    end
%    Tsum(ij)=blt_abv10(Ti(ij,1), Ti(ij,2), l)*dx(Ti(ij,1), Ti(ij,2))*dy(Ti(ij,1), Ti(ij,2));
%    %dsum(ij)=mask(Ti(ij,1), Ti(ij,2), l)*dx(Ti(ij,1), Ti(ij,2))*dy(Ti(ij,1), Ti(ij,2));
%    if ~isnan(Tsum(ij))
%      dsum(ij)=dx(Ti(ij,1), Ti(ij,2))*dy(Ti(ij,1), Ti(ij,2));
%    else
%      dsum(ij)=NaN;
%    end    
%  end
%  min(bltnan)
%  Tavg_abv10m(l)=nansum(Tsum)./nansum(dsum);
%end
%save('/work/uo0122/u241272/ATL24/area_avg_blt_8Nto13N_50Wto60W.mat','-v7.3','Tavg_abv10m');
%sum(isnan(Tsum))
%sum(isnan(dsum))
%clear Tsum dsum %mask


%regional mean of blt above 0m.
%blt_abv0=blt;
%blt_abv0(blt_abv0<=0)=NaN;   % try also without putting NaNs. 
%mask=blt_abv0;
%mask=mask*0+1;

%load atl24_clim_area.mat newblt2
%blt_abv0=newblt2;



%maxblt_abv0=zeros(12,1);
clim_en4_blt_avg=zeros(12,1);
clim_en4_blt5_avg=zeros(12,1);


for l=1:12
  Tsum=zeros(NumPoints,1);
  dsum=zeros(NumPoints,1);
  %bltnan_abv0=zeros(NumPoints,1);
  for ij=1:NumPoints
    %bltnan_abv0(ij)=blt_abv0(Ti(ij,1), Ti(ij,2), l);
    %if blt_abv0(Ti(ij,1), Ti(ij,2), l)>maxblt_abv0(l)
    %  maxblt_abv0(l)=blt_abv0(Ti(ij,1), Ti(ij,2), l);
    %end
    Tsum(ij)=blt_abv0(Ti(ij,1), Ti(ij,2), l)*111000*111000;  % 1deg = 111km = 111000m
    %dsum(ij)=mask(Ti(ij,1), Ti(ij,2), l)*111000*111000;
    if ~isnan(Tsum(ij))
      dsum(ij)=111000*111000;
    else
      dsum(ij)=NaN;
    end
  end
  %min(bltnan_abv0);
  clim_en4_blt_avg(l)=nansum(Tsum)./nansum(dsum);
  clear Tsum dsum 


  Tsum=zeros(NumPoints,1);
  dsum=zeros(NumPoints,1);
  %bltnan_abv0=zeros(NumPoints,1);
  for ij=1:NumPoints
    %bltnan_abv0(ij)=blt_abv0(Ti(ij,1), Ti(ij,2), l);
    %if blt_abv0(Ti(ij,1), Ti(ij,2), l)>maxblt_abv0(l)
    %  maxblt_abv0(l)=blt_abv0(Ti(ij,1), Ti(ij,2), l);
    %end
    Tsum(ij)=blt_abv5(Ti(ij,1), Ti(ij,2), l)*111000*111000;  % 1deg = 111km = 111000m
    %dsum(ij)=mask(Ti(ij,1), Ti(ij,2), l)*111000*111000;
    if ~isnan(Tsum(ij))
      dsum(ij)=111000*111000;
    else
      dsum(ij)=NaN;
    end
  end
  %min(bltnan_abv0);
  clim_en4_blt5_avg(l)=nansum(Tsum)./nansum(dsum);
  clear Tsum dsum 

end
%sum(isnan(Tsum));
%sum(isnan(dsum));

%save('atl24_clim_area.mat','clim_atl_6_10N_51_59W','-append');

%save('atl24_clim2003to2011_abv5_areavg.mat','clim_atl_6_10N_51_59W','-append');

%save('atl24_clim2003to2011_abv10_areavg.mat','clim_atl_6_10N_51_59W','-append');

%save('EN4_2_1_clim2003to2011_abv0_area_avg.mat','clim_en4_6_10N_51_59W','-append');

%save('EN4_2_1_reflev2_clim2003to2011_abv0_area_avg.mat','clim_en4_6_10N_51_59W','-append');

%save('/work/uo0122/u241272/ATL24/area_avg_blt_8Nto13N_50Wto60W.mat','Tavg_abv0m','-append');

outfile2=char(strcat('/work/uo0122/u241272/ATL24/Sal_temp_balance_terms/','en4_area_avg_',int2str(abs(lonw)),'_',int2str(abs(lone)),'W_',int2str(abs(lats)),'_',int2str(abs(latn)),'N.mat'));
save(outfile2,'clim_en4_blt_avg','clim_en4_blt5_avg');


outfile=char(strcat('/work/uo0122/u241272/ATL24/Sal_temp_balance_terms/','en4_area_avg_',int2str(abs(lonw)),'_',int2str(abs(lone)),'W_',int2str(abs(lats)),'_',int2str(abs(latn)),'N.nc'));



nccreate(outfile,'clim_en4_blt_avg','dimensions',{'time' 3287 'na' 1});     
ncwriteatt(outfile,'clim_en4_blt_avg','units','m'); 
ncwrite(outfile,'clim_en4_blt_avg',clim_en4_blt_avg);


nccreate(outfile,'clim_en4_blt5_avg','dimensions',{'time' 3287 'na' 1});      
ncwriteatt(outfile,'clim_en4_blt5_avg','units','m'); 
ncwrite(outfile,'clim_en4_blt5_avg',clim_en4_blt5_avg);



clear lone lonw lats latn clim_db_blt_avg clim_db_blt5_avg
%clear maxlon minlon maxlat minlat



toc
end %areai

exit


%mask(isnan(blt))=NaN;
%for l=1:108
%Tsum=0;
%dsum=0;
%for ij=1:NumPoints
%Tsum=Tsum+blt(Ti(ij,1), Ti(ij,2), l)*dx(Ti(ij,1), Ti(ij,2))*dy(Ti(ij,1), Ti(ij,2));
%dsum=dsum+dx(Ti(ij,1), Ti(ij,2))*dy(Ti(ij,1), Ti(ij,2));
%end
%Tavg(l)=Tsum/dsum;
%end
%regmean_blt = squeeze(nanmean(nanmean(blt_abv10(374:554,241:400,:),2),1)); %65 to 50W, 10 to 22N :Box1
