function [time, P]=ADVOFF_Model (Po, uvel, vvel, grd)

%% MAIN Define Integration parameters
dx=grd.dx;
dy=grd.dy;
K=0;
dt=1200*6;
day1=24*60*60;
trelax=day1*10;     %10 day relaxation for source
tdecay=day1*30*6;   % decay timescale set to 6 month

tstart=datenum(2000,1,1);
tend=datenum(2000,1,10);
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)]);


%% Configure storing parameters
save_interval=1; % days
nrec_length=save_interval*day1;
nstore=T*dt/nrec_length;
nfrac=nrec_length/dt;
[i,j]=size(uvel);
nrec=10;
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;
    
    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);
    
    % 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);
%         maxp=max(P(:));
%         caxis([-maxp maxp]); gradsmap4; colorbar 'h'
%         title (datestr(tcurrent));
%         pause(0.1)
%        
        myrec=myrec+1;
        PSTORE(:,:,myrec)=P;
        PTIME(myrec)=tcurrent;
         if maxval > 200
            disp('Blowup.');
        end
        iwrt=irec;            
    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');
% 

P=PSTORE;
time=PTIME;
return

%% Plotting

for ntime=1:10
    clf; rnc_map(PSTORE(:,:,ntime), grd);
    tmp=PSTORE(:,:,ntime);
     maxp=max(tmp(:));
     caxis([-maxp maxp]); gradsmap4; 
     colorbar 'h'
    title (datestr(PTIME(ntime)));
    pause(0.1)
end

%% Code used to generate reduced velocity.
%% Load velocity and grid information
clear
load VELOCITY_GRID.mat
% i=81:122;
% j=18:48;
load IJ.mat
uvel=uvel(i,j); vvel=vvel(i,j);
grd2.lon=grd.lon(i,j); grd2.lat=grd.lat(i,j); grd2.mask=grd.mask(i,j);
[I, J]=size(uvel);
grd2.I=I; grd2.J=J;
grd2.dx=grd.dx; grd2.dy=grd.dy;
grd=grd2;
save ADVOFF_data.mat grd uvel vvel