% Double Gyre Setup
% Forcing Function
Fu = @(time,p) ...
             -1000/2/pi*cos(2.*pi.*p(:,2))+0.*time+0.*p(:,1);
Fv = @(time,p) 0;
% Corriolis Coefficient f=f0+b0*y
app.Fcor = @(time,p)...
            0+1000.*p(:,2)+0.*p(:,2)+0.*time;

%# of interior CV's
a    = 1;          b    = 1;
dx   = a./(app.Nx);    dy   = b./(app.Ny); 
 bcsDP= [];        bcsNP= [0 0 0 0];     
 bcsDrho=[1 1 1 1]  ; bcsNrho = [];
 bcsDu= [0 0]  ;  bcsNu= [0 0]; 
 bcsDv= [0 0 0 0]  ;  bcsNv= [];        
 
Nt   = ceil(app.T/app.dt);

%% create geometry mask
Nbcs = length(bcsDP) + length(bcsNP);
bnd = [1 1 1 1];
Mask = zeros(app.Ny, app.Nx);

%create Node matrices
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,bnd,1);

%% Initialize fields
xp = linspace(dx/2, a-dx/2, app.Nx  );yp = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xu = linspace( dx   ,a-dx , app.Nx-1);yu = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xv = linspace(dx/2, a-dx/2, app.Nx  );yv = fliplr(linspace(dy    , b-dy, app.Ny-1));
[XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);
P = [bcsDP';bcsNP'; 0*idsP];
u = [bcsDu';bcsNu';idsu*0];
v = [bcsDv';bcsNv'; 0*idsv];
ucentid = Nodeu(:,floor(app.Nx/2+1));
vcentid = Nodev(:,floor(app.Nx/2+1));
