% 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
%%%%%Rainfall spatial
%{
lat1=double(ncread('/Users/aurpita.saha/Documents/asaha/data/rainfall_GES_DISC/3B-DAY.MS.MRG.3IMERG.20100122-S000000-E235959.V06.nc4','lat'));
lon1=double(ncread('/Users/aurpita.saha/Documents/asaha/data/rainfall_GES_DISC/3B-DAY.MS.MRG.3IMERG.20100122-S000000-E235959.V06.nc4','lon')); 
lat=lat1(950:1281); clear lat1
lon=lon1(790:1651); clear lon1
[Ry,Rx]=meshgrid(lat,lon);
%}
% satellite SSS %[y,x]=meshgrid(lat,lon);
Slon1=double(ncread('/Users/aurpita.saha/Documents/asaha/data/SMOS/SMOS_daily/SM_OPER_MIR_CSF2QD_20211231T000000_20211231T235959_331_001_7.DBL.nc','lon'));
Slat1=double(ncread('/Users/aurpita.saha/Documents/asaha/data/SMOS/SMOS_daily/SM_OPER_MIR_CSF2QD_20211231T000000_20211231T235959_331_001_7.DBL.nc','lat'));
%Slat=Slat1(317:473); clear Slat1
%Slon=Slon1(305:637); clear Slon1

%{
%%%%%%Oscar uv spatial% oscar
Olat1=ncread('/Users/aurpita.saha/Documents/asaha/data/oscar_data/oscar_vel2010.nc','latitude');
Olon=ncread('/Users/aurpita.saha/Documents/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
[Oy,Ox]=meshgrid(Olat,Olon);

%good_valuesR=find(~isnan(Rx));
good_valuesS=find(~isnan(Sx));
%}
%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 %:2013 %2021 %2022	%:2022


ii=121 %121 %193 	%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_S/TCloc_Sbudget_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(Slat1>=latS); indlaS=r(1); clear r %indlaS(indlaS==1)=2; clear r
r=find(Slat1<=latN); indlaN=r(end); clear r %indlaN(indlaN==1)=2; clear r
r=find(Slon1>=lonW); indloW=r(1); clear r %indloW(indloW==1)=2; clear r
r=find(Slon1<=lonE); indloE=r(end); clear r %indloE(indloE==1)=2; clear r

%[Sy,Sx]=meshgrid(Slat,Slon);
%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
infile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/satellite_Amazon/satelliteSSS_SmosSmap_',DatesIn3(ind+i),'.mat'));
if isfile(infile)
load(infile,'satSr'); 
%before(:,:,i)=satSr(305:637,317:473);
before(:,:,i)=satSr(indloW:indloE,indlaS:indlaN);
clear satSr 
else
before(:,:,i)=NaN(indloE-indloW+1,indlaN-indlaS+1);
end 
end %i
clear ind 
p=p+1
SSS_sat(:,:,p)=nanmean(before,3); clear before
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
infile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/satellite_Amazon/satelliteSSS_SmosSmap_',DatesIn3(ind+i),'.mat'));
if isfile(infile)
load(infile,'satSr'); 
%after(:,:,i)=satSr(305:637,317:473);
after(:,:,i)=satSr(indloW:indloE,indlaS:indlaN);
clear satSr 
else
%after(:,:,i)=NaN(637-305+1,473-317+1);
after(:,:,i)=NaN(indloE-indloW+1,indlaN-indlaS+1);
end
end %i
clear ind %ind2
p=p+1
SSS_sat(:,:,p)=nanmean(after,3); clear after
end %j
%outfile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/TCall/all/trial_onlyS_satellite_oscar_year',int2str(bigyear),'_i',int2str(ii),'_j',int2str(jj),'.mat'));
%save(outfile,'-v7.3','SSS_sat');
%clear SSS_sat

%{
%%%%%%%%%%%%%%%%%%%%%%%%%&&& centered &&&%%%%%%%%%%%%%%%%%%%%%%%%
%S before 
ind=indmaha-5; %find(DatesIn3=='20100831'); 
c=0;
ind2=find(DatesIn5==DatesIn3(ind));
for i=1:3
infile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/satellite_Amazon/satelliteSSS_SmosSmap_',DatesIn3(ind+i),'.mat'));
load(infile,'satSr');
before(:,:,i)=satSr(305:637,317:473);
clear satSr 
a=char(DatesIn5(ind2+i));
yr=strcat(a(1),a(2),a(3),a(4));
outfile1=char(strcat('/Users/aurpita.saha/Documents/asaha/data/oscar_data/oscar_vel',yr,'.nc'));
time=ncread(outfile1,'time');
days=DatesIn5(time+1);
indd=find(days==DatesIn5(ind2+i));	
if isempty(indd)
continue
end %if
c=c+1;
%DatesIn5(ind2+i)
uu=ncread(outfile1,'u',[718 127 1 indd],[976-718+1 226-127+1 1 1],[1 1 1 1]);
vv=ncread(outfile1,'v',[718 127 1 indd],[976-718+1 226-127+1 1 1],[1 1 1 1]);
beforeu(:,:,c)=uu;
beforev(:,:,c)=vv;
clear uu vv indd outfile1
end %i

ind=indmaha; %find(DatesIn3=='20100905');
ind2=find(DatesIn5==DatesIn3(ind));
for i=1:3
infile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/satellite_Amazon/satelliteSSS_SmosSmap_',DatesIn3(ind+i),'.mat'));
load(infile,'satSr');
before(:,:,i+3)=satSr(305:637,317:473);
clear satSr 
a=char(DatesIn5(ind2+i));
yr=strcat(a(1),a(2),a(3),a(4));
outfile1=char(strcat('/Users/aurpita.saha/Documents/asaha/data/oscar_data/oscar_vel',yr,'.nc'));
time=ncread(outfile1,'time');
days=DatesIn5(time+1);
indd=find(days==DatesIn5(ind2+i));	
if isempty(indd)
continue
end %if
c=c+1;
%DatesIn5(ind2+i)
uu=ncread(outfile1,'u',[718 127 1 indd],[976-718+1 226-127+1 1 1],[1 1 1 1]);
vv=ncread(outfile1,'v',[718 127 1 indd],[976-718+1 226-127+1 1 1],[1 1 1 1]);
beforeu(:,:,c)=uu;
beforev(:,:,c)=vv;
clear uu vv indd outfile1
end %i
%clear c indd ind ind2
SSSc_sat(:,:,1)=nanmean(before,3); clear before
if exist('beforeu','var')==0
u_oscar(:,:,1)=NaN(976-718+1,226-127+1); %mean(beforeu,3); clear beforeu
v_oscar(:,:,1)=NaN(976-718+1,226-127+1); %mean(beforev,3); clear beforev
else
u_oscar(:,:,1)=nanmean(beforeu,3); clear beforeu
v_oscar(:,:,1)=nanmean(beforev,3); clear beforev
end

indf=indmaha+2; p=1; %find(DatesIn3=='20100907'); 
for j=1:6:120 %360
a=char(DatesIn3(indf+j)); 
ind=find(DatesIn3==a); ind2=find(DatesIn5==DatesIn3(indf+j)); c=0;
for i=1:6
infile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/satellite_Amazon/satelliteSSS_SmosSmap_',DatesIn3(ind+i),'.mat'));
load(infile,'satSr');
a=char(DatesIn3(ind+i));
after(:,:,i)=satSr(305:637,317:473);
clear satSr 
yr=strcat(a(1),a(2),a(3),a(4));
outfile1=char(strcat('/Users/aurpita.saha/Documents/asaha/data/oscar_data/oscar_vel',yr,'.nc'));
time=ncread(outfile1,'time');
days=DatesIn5(time+1);
indd=find(days==DatesIn5(ind2+i));	
if isempty(indd)
continue
end %if
c=c+1;
%DatesIn5(ind2+i)
uu=ncread(outfile1,'u',[718 127 1 indd],[976-718+1 226-127+1 1 1],[1 1 1 1]);
vv=ncread(outfile1,'v',[718 127 1 indd],[976-718+1 226-127+1 1 1],[1 1 1 1]);
beforeu(:,:,c)=uu;
beforev(:,:,c)=vv;
clear uu vv indd outfile1
end %i
clear c indd ind ind2
p=p+1;
SSSc_sat(:,:,p)=nanmean(after,3); clear after
if exist('beforeu','var')==0
u_oscar(:,:,p)=NaN(976-718+1,226-127+1); %mean(beforeu,3); clear beforeu
v_oscar(:,:,p)=NaN(976-718+1,226-127+1); %mean(beforev,3); clear beforev
else
u_oscar(:,:,p)=nanmean(beforeu,3); clear beforeu
v_oscar(:,:,p)=nanmean(beforev,3); clear beforev
end
end %j

outfile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/TCall/SuvRain120_satellite_oscar_i',int2str(ii),'_j',int2str(jj),'.mat'));
save(outfile,'-v7.3','SSS_sat','SSSc_sat','u_oscar','v_oscar');
clear SSS_sat SSSc_sat u_oscar v_oscar

%%%%%%%%%%%rain
ind=indmaha-5; %ind=find(DatesIn3=='20100831'); %c=0;
for i=1:3
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));
%read rainfall
infile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/rainfall_GES_DISC/3B-DAY.MS.MRG.3IMERG.',yrnew,mnew,daynew,'-S000000-E235959.V06.nc4'));
if ~isfile(infile)
continue
end %if
rain=ncread(infile,'precipitationCal',[950 790 1],[1281-950+1 1651-790+1 1],[1 1 1]); 
beforeP(:,:,i)=rain'; clear rain
end %i
ind=indmaha; %ind=find(DatesIn3=='20100905'); %c=0;
for i=1:3
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));
%read rainfall
infile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/rainfall_GES_DISC/3B-DAY.MS.MRG.3IMERG.',yrnew,mnew,daynew,'-S000000-E235959.V06.nc4'));
if ~isfile(infile)
continue
end %if
rain=ncread(infile,'precipitationCal',[950 790 1],[1281-950+1 1651-790+1 1],[1 1 1]); 
beforeP(:,:,i+3)=rain'; clear rain
end %i
beforeP(beforeP==0)=NaN;
rain_sat(:,:,1)=nanmean(beforeP,3); clear beforeP

indf=indmaha+2; p=1; %find(DatesIn3=='20100907'); 
for j=1:6:120 %360
a=char(DatesIn3(indf+j));   %mday=strcat(a(5),a(6),a(7),a(8));
% 6-11th Sep after TC; - 6 days ---------1 map % every 6 days check after maps for 2 months
ind=find(DatesIn3==a); %c=0;
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));
%read rainfall
infile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/rainfall_GES_DISC/3B-DAY.MS.MRG.3IMERG.',yrnew,mnew,daynew,'-S000000-E235959.V06.nc4'));
if ~isfile(infile)
continue
end %if
rain=ncread(infile,'precipitationCal',[950 790 1],[1281-950+1 1651-790+1 1],[1 1 1]); %rain(rain==0)=NaN;
beforeP(:,:,i)=rain';
clear rain
end %i
%clear c
p=p+1;
beforeP(beforeP==0)=NaN;
rain_sat(:,:,p)=nanmean(beforeP,3); clear beforeP
end %j

outfile=char(strcat('/Users/aurpita.saha/Documents/asaha/data/TCall/SuvRain120_satellite_oscar_i',int2str(ii),'_j',int2str(jj),'.mat'));
save(outfile,'rain_sat','-append');
clear rain_sat
%}

%%%%%%%%%%%%%%%%%%%%%%%% interpolate to oscar grid
%load(outfile,'SSS_sat');
%for n=1:61+20 %21 p-1
%n
%{
%%%%%%%%%%% select time step  
orig_fieldR=squeeze(rain_sat(:,:,n));
%%%%%interpolate R to oscar grid % build interpolant % evaluate field in the destination grid % store in matrix % clear temporary variables
F=scatteredInterpolant(Rx(good_valuesR),Ry(good_valuesR),orig_fieldR(good_valuesR),'linear');
dest_field=F(Ox,Oy);
rain_in_oscar(:,:,n)=dest_field;
clear F dest_field orig_fieldR %rain_sat

orig_fieldS=double(squeeze(SSSc_sat(:,:,n)));
F=scatteredInterpolant(Sx(good_valuesS),Sy(good_valuesS),orig_fieldS(good_valuesS),'linear');
dest_field=F(Ox,Oy);
satSc_in_oscar(:,:,n)=dest_field;
clear F dest_field orig_fieldS %SSS_clim

%%%%%% sat S to oscar
orig_fieldS=squeeze(SSS_sat(:,:,n));
F=scatteredInterpolant(Sx(good_valuesS),Sy(good_valuesS),orig_fieldS(good_valuesS),'linear');
dest_field=F(Ox,Oy);
satS_in_oscar(:,:,n)=dest_field;
clear F dest_field orig_fieldS %SSS_sat
end %n

n=62; %22; %p
%%%%%% sat S to oscar
orig_fieldS=squeeze(SSS_sat(:,:,n));
F=scatteredInterpolant(Sx(good_valuesS),Sy(good_valuesS),orig_fieldS(good_valuesS),'linear');
dest_field=F(Ox,Oy);
satS_in_oscar(:,:,n)=dest_field;
clear F dest_field orig_fieldS %SSS_sat
%}
%clear SSS_sat n

%save interpolated fields in mat file save in matlab format %outfile=char('/Users/aurpita.saha/Documents/asaha/data/rain_SatSinOscar_grid.mat');
%save(outfile,'satS_in_oscar','-append'); 

%%%% make balance terms
%load(outfile,'u_oscar','v_oscar');

%%% make the balance terms for that lat-lon box
%SSS tendency and no tendency, and wake
%LHS= satS_in_oscar(indloW:indloE,indlaN:indlaS,2:end)-satS_in_oscar(indloW:indloE,indlaN:indlaS,1:end-1);  
%for l=1:22 %p %16
%S_wake(:,:,l)= satS_in_oscar(indloW:indloE,indlaN:indlaS,l)-satS_in_oscar(indloW:indloE,indlaN:indlaS,1); 
%end %l
S_raw=squeeze(nanmean(nanmean(SSS_sat)));

%{
Sc_raw=squeeze(nanmean(nanmean(satSc_in_oscar(indloW:indloE,indlaN:indlaS,:))));

% E-P precipitation surface freshwater flux term
P=rain_in_oscar(indloW:indloE,indlaN:indlaS,:)./(6*86400*1000);  
E=zeros(size(P,1),size(P,2),21); %zeros(1,731)
rho=1029; h=30; 
nm_dep= (E-P).*satSc_in_oscar(indloW:indloE,indlaN:indlaS,:); %(:,:,2:end); 
Fwflux= nm_dep./h;

% horizontal adv term % -u(dS/dx)-v(dS/dy)
S=satSc_in_oscar(indloW-1:indloE+1,indlaN-1:indlaS+1,:);
Fx=diff(S)./(0.3333*106000);
Fy=diff(S,1,2)./(0.3333*111000);
dSx=(Fx(1:end-1,:,:)+Fx(2:end,:,:))./2;
dSy=(Fy(:,1:end-1,:)+Fy(:,2:end,:))./2;
udSxf=-(u_oscar(indloW:indloE,indlaN:indlaS,:).*dSx(:,2:end-1,:)); 
vdSyf=-(v_oscar(indloW:indloE,indlaN:indlaS,:).*dSy(2:end-1,:,:));
Horadv=udSxf+vdSyf;
%}
%save
save(outfile,'-v7.3','S_raw');

clear SSS_sat S_raw outfile

end %jj

ii=ii+1
end %m
end %yr


















