function [Pstore, Ptime, grd]=ADVOFF_main_adj(YEAR, exp_name) 
% ADVOFF_main_adj(YEAR, exp_name)
% clear
% YEAR=1980;
% exp_name='TR';

PLOTTING=0;
%% Load Grid and Data Structure Array
[grd, ctl]=ADVOFF_GetGrid_DataStruct; dx=grd.dx; dy=grd.dy;

%% MAIN Define Integration parameters
K=0;
dt=1200;
day1=24*60*60;
dt=1200;
tbackward=360; % days
tstart=datenum(YEAR,2,25);
tend=datenum(YEAR,10,30);
T=floor((tend-tstart)*day1/dt);
disp(['Start day :  ', datestr(tstart)]);
disp(['End   day :  ', datestr(tend)]);
disp(['Total integration in days = ', num2str(dt*T/day1)]);
disp(['Total integration tsteps = ', num2str(T)]);



%% Load velocity data
disp('Loading velocity data ...');
istr=find(ctl.time < tstart);
iend=find(ctl.time > tend);
in=istr(end):iend(1);
U=rnt_loadvar(ctl,in, 'u');
V=rnt_loadvar(ctl,in, 'v');
UV_TIME=ctl.time(in);
U=U/100;
V=V/100;
U(isnan(U))=0;
V(isnan(V))=0;



%% Initialize Passive tracer
Po=U(:,:,1)*0; 

% EXP tracer boxes inverse

% Transition Region 36-40N 143-150W
if strcmp(exp_name, 'TR');
    xlim=[143 150]-360; ylim=[36 40];
    [Isource,Jsource]=rgrd_FindIJ(grd.lon, grd.lat, xlim, ylim);
    Po(Isource,Jsource)=100;
end


% % Oyashio Region 40-43N 143-150W
if strcmp(exp_name, 'OY');
    xlim=[143 150]-360; ylim=[40 44];
    [Isource,Jsource]=rgrd_FindIJ(grd.lon, grd.lat, xlim, ylim);
    Po(Isource,Jsource)=100;
end

trelax=day1*10; %10 day relaxation for source
tpersistence=30*4; %days
ndecorr=tpersistence*day1/dt;
for ntime=1:T
    tcurrent(ntime)=tstart+ (dt*(T-ntime))/day1;
end
%tfilter=exp(-[1:T]/ndecorr);
tfilter=tcurrent*0;
[year,month]=dates_datenum(tcurrent);
in=find(month>=5);
tfilter(in)=1;

if PLOTTING==1
    mfig(1); clf;
    subplot(2,1,1);
    rnc_map(Po, grd);
    subplot(2,1,2);
    plot(tcurrent, tfilter); datetick
end


%% Configure storing parameters
save_interval=5; % days
nrec_length=save_interval*day1;
nstore=T*dt/nrec_length;
nfrac=nrec_length/dt;
[i,j]=size(Po);
Pstore=zeros(i,j,nstore+1);
Ptime=zeros(nstore+1,1);


%% Run Advection Model (inverse)
P=Po;

% Inverse Calculation

iwrt=0;
for ntime=1:T
    tcurrent=tstart+ (dt*(T-ntime))/day1;
    
    % Load velocity field
    
    istr=find( tcurrent > UV_TIME); istr=istr(end);
    iend=find( tcurrent <= UV_TIME); iend=iend(1);
    c1=(UV_TIME(iend)-tcurrent)/(UV_TIME(iend)-UV_TIME(istr));
    c2=(tcurrent-UV_TIME(istr))/(UV_TIME(iend)-UV_TIME(istr));
    uvel=U(:,:,istr)*c1 + U(:,:,iend)*c2;
    vvel=V(:,:,istr)*c1 + V(:,:,iend)*c2;
    uvel=-uvel;
    vvel=-vvel;
    
    u=P;
    RHS = ADVOFF_rhs(u,uvel,vvel,dx,dy,K,grd);
    
    P =  u + RHS*dt;
    P = ADVOFF_bc( P );
    
    irec=1+floor(dt*ntime/nrec_length);
    Pstore(:,:,irec)=Pstore(:,:,irec) + P/nfrac;
    Ptime(irec)=Ptime(irec)+tcurrent/nfrac;
    
    % Diagnostic and outout storing
    if irec > iwrt        
        maxval=max(P(:)); minval=min(P(:));
        disp(['Running timestep ',num2str(ntime), ' REC= ', num2str(irec),'  TIME= ', datestr(tcurrent) ...
            ,'   MinVal = ',num2str(minval), '    MaxVal = ', num2str(maxval)]);
%         clf; rnc_map(P, grd);
%         caxis([-110 110]); gradsmap4; colorbar 'h'
%         title (datestr(tcurrent));
%         pause(0.1)
        if maxval > 200
            disp('Blowup.');
        end
        iwrt=irec;
        
    end
    % Apply nudging in source region
    P(Isource,Jsource)=P(Isource,Jsource)+(Po(Isource,Jsource)-P(Isource,Jsource))*dt/trelax*tfilter(ntime);
    
end


%% Post proessing tasks
Pstore=Pstore(:,:,1:end-1); Ptime=Ptime(1:end-1);

str=[num2str(YEAR),'_Pstore_', exp_name, '.mat'];
save(str, 'Pstore','Ptime');

return

%% Plotting

for ntime=1:length(Ptime)
    clf; rnc_map(Pstore(:,:,ntime), grd);
    caxis([0 110]); colorbar 'h'
    title (datestr(Ptime(ntime)));
    pause(0.1)
end