function [Pstore, Ptime, grd]=ADVOFF_main_fwd_long(exp_name)
% ADVOFF_main_adj(YEAR, exp_name)
% clear
% YEAR=1980;
%clear
%exp_name='SB';
%
PLOTTING=0;
%% Load Grid and Data Structure Array
[grd, ctl]=ADVOFF_GetGrid_DataStruct; dx=grd.dx; dy=grd.dy;
[L,P]=size(grd.lon);
grd.L=L;
grd.P=P;

%% MAIN Define Integration parameters
K=0;
dt=1200;
day1=24*60*60;
dt=1200;
tstart=datenum(1966,2,25);
tend=datenum(2000,9,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 ...');
vincrement=30;
istr=find(ctl.time < tstart);
iend=find(ctl.time > tend);
in=iend(1)-vincrement: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

% South Region
if strcmp(exp_name, 'SB');
    xlim=[-228.5265 -221.0984]; ylim=[28.7471   34.1533];
    [Isource,Jsource]=rgrd_FindIJ(grd.lon, grd.lat, xlim, ylim);
    Po(Isource,Jsource)=100;
end

if strcmp(exp_name, 'NB');
    xlim=[-214.6790 -207.8012]; ylim=[42.4514   46.7569];
    [Isource,Jsource]=rgrd_FindIJ(grd.lon, grd.lat, xlim, ylim);
    Po(Isource,Jsource)=100;
end

if strcmp(exp_name, 'KOE');
    xlim=[ -221.6486 -217.3385]; ylim=[32.9792   35.9069];
    [Isource,Jsource]=rgrd_FindIJ(grd.lon, grd.lat, xlim, ylim);
    Po(Isource,Jsource)=100;
end

% 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

% Transition Region 36-40N 143-150W
xlim1=[143 150]-360; ylim1=[36 40];

% Oyashio Region 40-43N 143-150W
xlim2=[143 150]-360; ylim2=[40 44];



trelax=day1*10; %10 day relaxation for source
tdecay=day1*30*6; % decay timescale set to 6 month

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);
nrec=30;
Pstore=zeros(i,j);
Ptime=zeros(1);
PSTORE=zeros(i,j, nrec);
PTIME=zeros(nrec,1);
myrec=0;

ncount=0;


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

% Forward Calculation

iwrt=0;
for ntime=1:T
    %tcurrent=tstart+ (dt*(ntime-1))/day1;
    tcurrent=tstart+ (dt*(T-ntime))/day1;
    
    if tcurrent-5 < UV_TIME(1)
        tstart1=tcurrent+5;
        % Load velocity data
        disp('Loading velocity data ...');
        %pause
        istr=find(ctl.time > tstart1);
        in=istr(1)-vincrement:istr(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;
    end
    
    
    % 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;  % Adjoint
    vvel=-vvel;
    
    u=P;
    RHS = ADVOFF_rhs(u,uvel,vvel,dx,dy,K,grd);
    
    P =  u + RHS*dt - u*dt/tdecay;
    P = ADVOFF_bc( P );
    
    irec=1+floor(dt*ntime/nrec_length);
    ncount=ncount+1;
    Pstore=Pstore + P;
    Ptime=Ptime+tcurrent;
    
    % 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));
        %         % TR box
        %         xlim=xlim1; ylim=ylim1;
        %         x=[xlim(1) xlim(1) xlim(2) xlim(2) xlim(1)];
        %         y=[ylim(1) ylim(2) ylim(2) ylim(1) ylim(1)];
        %         plot(x,y,'color','k', 'linewidth', 2);
        %         % OY box
        %         xlim=xlim2; ylim=ylim2;
        %         x=[xlim(1) xlim(1) xlim(2) xlim(2) xlim(1)];
        %         y=[ylim(1) ylim(2) ylim(2) ylim(1) ylim(1)];
        %         plot(x,y,'color','k', 'linewidth', 2);
        %
        %
        %         pause(0.1)
        Pstore=Pstore/ncount;
        Ptime=Ptime/ncount;
        % store average
        myrec=myrec+1;
        PSTORE(:,:,myrec)=Pstore;
        PTIME(myrec)=Ptime;
        
        if myrec==nrec
            filename=[datestr(tcurrent, 29),'_Pstore_LONG_ADJ', exp_name, '.nc'];
            grd.T=length(PTIME);                                  
            Convert2NC_tracer(grd,filename);
            nc=netcdf(filename, 'w');
            nc{'time'}(:)=PTIME;
            nc{'P'}(:,:,:)=perm(PSTORE);
            close(nc)
            %save(str, 'PSTORE','PTIME');
            disp(['----- Writing to disk ',num2str(ntime), ' REC= ', num2str(irec),'  TIME= ', datestr(Ptime)]);
            PSTORE(:)=0;
            PTIME(:)=0;
            myrec=0;
        end
        
        if maxval > 200
            disp('Blowup.');
        end
        iwrt=irec;
        ncount=0;
        Pstore(:)=0;
        Ptime=0;
        
        
    end
    
    
    % Apply nudging in source region
    %if str2num(datestr(tcurrent,'mm'))<4
    P(Isource,Jsource)=P(Isource,Jsource)+(Po(Isource,Jsource)-P(Isource,Jsource))*dt/trelax;
    %end
    
    
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