% make for all TC locations the S raw, budget terms, for eyery year 2010-2021 

%Dates for input filenames
DatesIn2=strings(365,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=2005 %:2016
  for imonth=1:12
    if iyear==2004 || iyear==2008 || iyear==2012 || iyear==2016 || iyear==2020
      for iday=1:MonthDaysLeap(imonth)
        DatesIn2(idate)=sprintf('%02d%02d', imonth, iday);%sprintf('%02d-%02d-%d',iday, imonth, iyear);
        idate=idate+1;
      end
    else
      for iday=1:MonthDays(imonth)
        DatesIn2(idate)=sprintf('%02d%02d', imonth, iday);
        idate=idate+1;
      end
    end
  end
end

DatesIn3=strings(8401,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=2000:2022
  for imonth=1:12
    if iyear==1996 || iyear==2000 || iyear==2004 || iyear==2008 || iyear==2012 || iyear==2016 || iyear==2020
      for iday=1:MonthDaysLeap(imonth)
        DatesIn3(idate)=sprintf('%04d%02d%02d%d', iyear, imonth, iday);%sprintf('%02d-%02d-%d',iday, imonth, iyear);
        idate=idate+1;
      end
    else
      for iday=1:MonthDays(imonth)
        DatesIn3(idate)=sprintf('%04d%02d%02d%d', iyear, imonth, iday);
        idate=idate+1;
      end
    end
  end
end

%for OSCAR uv currents- horizontal advection term
DatesIn4=strings(11323,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=1992:2022
  for imonth=1:12
    if iyear==1992 || iyear==1996 || iyear==2000 || iyear==2004 || iyear==2008 || iyear==2012 || iyear==2016 || iyear==2020
      for iday=1:MonthDaysLeap(imonth)
        DatesIn4(idate)=sprintf('%04d%02d%02d%d', iyear, imonth, iday);%sprintf('%02d-%02d-%d',iday, imonth, iyear);
        idate=idate+1;
      end
    else
      for iday=1:MonthDays(imonth)
        DatesIn4(idate)=sprintf('%04d%02d%02d%d', iyear, imonth, iday);
        idate=idate+1;
      end
    end
  end
end
DatesIn5=DatesIn4(279:end);

%loading grids for interpolation to oscar

%sat SST 
Tlon1=[-179.956:0.08787:179.956];
Tlat1=[-89.956:0.08787:89.956];
%x_dim=[-179.956:0.08787:179.956];
%y_dim=[-89.956:0.08787:89.956];
%Tlon=x_dim(900:1860); clear x_dim
%Tlat=y_dim(1000:1550); clear y_dim
%{
%%%%%%Oscar uv spatial% oscar
Olat1=ncread('/phoddat/share/asaha/data/oscar_data/oscar_vel2010.nc','latitude');
Olon=ncread('/phoddat/share/asaha/data/oscar_data/oscar_vel2010.nc','longitude');
new_lon=Olon; new_lon(new_lon>180)=new_lon(new_lon>180)-360; clear Olon
Olon1=new_lon; clear new_lon
Olat=Olat1(127:226); clear Olat1
Olon=Olon1(718:976); clear Olon1
%}
fid=fopen('/Users/aurpita.saha/Documents/asaha/data/sst_mw_ir_3dy.dat','r');
xdim=4096;
ydim=2048;

%loading TC locs
%load('/Users/aurpita.saha/Documents/asaha/data/mm_loc_TC33_2000_2022_atlantic.mat') %,'wind_TC','transpeed_TC');%lat_TC(lat_TC==0)=NaN;lon_TC(lon_TC==0)=NaN;
%load('/Users/aurpita.saha/Documents/asaha/data/mm_loc_TC33_2000_2022_Nindian.mat') %,'wind_TC','transpeed_TC');%lat_TC(lat_TC==0)=NaN;lon_TC(lon_TC==0)=NaN;
load('/Users/aurpita.saha/Documents/asaha/data/mm_loc_TC33_2000_2022_atlantic.mat') %,'wind_TC','transpeed_TC');%lat_TC(lat_TC==0)=NaN;lon_TC(lon_TC==0)=NaN;

a=nanmean(lat_TC);
aj=find(a~=0);
clear a


bigyear=2010 %:2020

ii=121 %121 %start year and month = %2010-  :252 = 2020

for yrr=2010:2020 %133=2011; 145=2012; 157=2013; 169=2014; 181=2015; 193=2016; 205=2017; 217=2018; 229=2019; 241=2020
yrr

for m=1:12 
%m
%year=sprintf('%04d', yr);
%month=sprintf('%02d', m);

for jj=1:length(aj)

outfile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/TCall/balterms/Natlantic_T/TCloc_Tbudget_year',int2str(bigyear),'_i',int2str(ii),'_j',int2str(jj),'.mat'));
if isfile(outfile)
continue
end %if

  if new_profilenoi(ii,jj)==0
     continue
  end %if
  jj

datet=sprintf('%04d%02d%02d%d', bigyear, m, day(ii,jj)); %yrr
indmaha=find(DatesIn3==datet);


%define TC box
latN=lat_TC(ii,jj)+3; latS=lat_TC(ii,jj)-3;	
	if lon_TC(ii,jj)==177
		lonW=177-3; lonE=179.8;
	elseif lon_TC(ii,jj)==178
		lonW=178-3; lonE=-179;
	elseif lon_TC(ii,jj)==179
		lonW=179-3; lonE=-178;
	elseif lon_TC(ii,jj)==180
		lonW=180-3; lonE=-177;
	elseif lon_TC(ii,jj)==-179
		lonW=178; lonE=-179+3;
	elseif lon_TC(ii,jj)==-178
		lonW=179; lonE=-178+3;
	elseif lon_TC(ii,jj)==-177
		lonW=-179.8; lonE=-177+3;
	else 
		lonW=lon_TC(ii,jj)-3; lonE=lon_TC(ii,jj)+3; 
	end %if

r=find(Tlat1>=latS); indlaS=r(1); clear r %indlaS(indlaS==1)=2; 
r=find(Tlat1<=latN); indlaN=r(end); clear r %indlaN(indlaN==1)=2; 
r=find(Tlon1>=lonW); indloW=r(1); clear r %indloW(indloW==1)=2; 
r=find(Tlon1<=lonE); indloE=r(end); clear r %indloE(indloE==1)=2; 
%{
[Ty,Tx]=meshgrid(Tlat1(indlaN:indlaS),Tlon1(indloW:indloE));
r=find(Olat1>=latS); indlaSo=r(end); indlaSo(indlaSo==1)=2; clear r
r=find(Olat1<=latN); indlaNo=r(1); indlaNo(indlaNo==1)=2; clear r
r=find(Olon1>=lonW); indloWo=r(1); indloWo(indloWo==1)=2; clear r
r=find(Olon1<=lonE); indloEo=r(end); indloEo(indloEo==1)=2; clear r
[Oy,Ox]=meshgrid(Olat1(indlaNo:indlaSo),Olon1(indloWo:indloEo));
clear indlaNo indlaSo indloWo indloEo latS latN lonW lonE
good_valuesT=find(~isnan(Tx));
%}
%6 days episode for S- 2 days before TC
p=0;
for j=-122:6:-8 %126 %366 %for 1/2 a year span
a=char(DatesIn3(indmaha+j)); 
ind=find(DatesIn3==a); 
for i=1:6
a=char(DatesIn3(ind+i));
daynew=strcat(a(7),a(8));
mnew=strcat(a(5),a(6));
yrnew=strcat(a(1),a(2),a(3),a(4)); clear a
%read SST
daynew=str2num(daynew);
mnew=str2num(mnew);
yrnew=str2num(yrnew);
t_get=datenum(yrnew,mnew,daynew)-datenum(2002,6,1)+1;
fseek(fid,4*xdim*ydim*(t_get-1),-1);
a=fread(fid,xdim*ydim,'float');
sst=reshape(a,[xdim,ydim]);
%satT=sst(900:1860,1000:1550); clear sst a t_get;
satT=sst(indloW:indloE,indlaS:indlaN); clear sst a t_get;
satT(satT>32)=NaN; satT(satT<=0)=NaN;
satT(satT==-9999)=NaN;
beforeSST(:,:,i)=satT;
clear satT 
end %i
clear ind 
p=p+1
SST_sat(:,:,p)=nanmean(beforeSST(1:68,1:68,:),3); clear beforeSST
end %j

%6 days episode for S after TC
indf=indmaha-1; %find(DatesIn3=='20100904'); p=12;
for j=1:6:366 %126 %366 %for 1/2 a year span
a=char(DatesIn3(indf+j)); 
ind=find(DatesIn3==a); 
for i=1:6
a=char(DatesIn3(ind+i));
daynew=strcat(a(7),a(8));
mnew=strcat(a(5),a(6));
yrnew=strcat(a(1),a(2),a(3),a(4)); clear a
%read SST
daynew=str2num(daynew);
mnew=str2num(mnew);
yrnew=str2num(yrnew);
t_get=datenum(yrnew,mnew,daynew)-datenum(2002,6,1)+1;
fseek(fid,4*xdim*ydim*(t_get-1),-1);
a=fread(fid,xdim*ydim,'float');
sst=reshape(a,[xdim,ydim]);
%satT=sst(900:1860,1000:1550); clear sst a t_get;
satT=sst(indloW:indloE,indlaS:indlaN); clear sst a t_get;
satT(satT>32)=NaN; satT(satT<=0)=NaN;
satT(satT==-9999)=NaN;
after(:,:,i)=satT;
clear satT
end %i
clear ind %ind2
p=p+1
SST_sat(:,:,p)=nanmean(after(1:68,1:68,:),3); clear after
end %j
%outfile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/TCall/balterms/atlantic_T/TCloc_SSTsatellite_oscar_year',int2str(bigyear),'_i',int2str(ii),'_j',int2str(jj),'.mat'));
%save(outfile,'-v7.3','SST_sat');
%clear SST_sat

%{
%%%%%%%%%%%%%%%%%%%%%%%% interpolate to oscar grid
%load(outfile,'SST_sat');
for n=1:61+20 %21 p-1
n
%%%%%% sat T to oscar
orig_fieldT=squeeze(SST_sat(:,:,n));
F=scatteredInterpolant(Tx(good_valuesT),Ty(good_valuesT),orig_fieldT(good_valuesT),'linear');
dest_field=F(Ox,Oy);
satT_in_oscar(:,:,n)=dest_field;
clear F dest_field orig_fieldT 
end %n
clear n 
%save interpolated fields in mat file save in matlab format 
save(outfile,'satT_in_oscar','-append'); 
%}

%%% make the balance terms for that lat-lon box
%SSS tendency and no tendency, and wake
%LHS= satT_in_oscar(:,:,2:end)-satT_in_oscar(:,:,1:end-1);  
%tLHS= SST_sat(:,:,2:end)-SST_sat(:,:,1:end-1);

%for l=1:22 %p %16
%S_wake(:,:,l)= satT_in_oscar(indloW:indloE,indlaN:indlaS,l)-satT_in_oscar(indloW:indloE,indlaN:indlaS,1); 
%end %l

%T_raw=squeeze(nanmean(nanmean(satT_in_oscar)));
T_raw=squeeze(nanmean(nanmean(SST_sat)));

%save
%outfile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/TCall/balterms/Natlantic_T/TCloc_Tbudget_year',int2str(bigyear),'_i',int2str(ii),'_j',int2str(jj),'.mat'));
save(outfile,'-v7.3','T_raw');

%save(outfile,'-v7.3','LHS','tLHS','T_raw','tT_raw');

clear outfile SST_sat T_raw

%clear satT_in_oscar T_raw outfile LHS SST_sat tLHS tT_raw

end %jj

ii=ii+1
end %m
end %yr


















