function [P u v rho Tsw Ssw NodeP Nodeu Nodev] =...
    BoussinesqTS(app,SetupScript,PlotScript)
% function [P u v rho Tsw Ssw NodeP Nodeu Nodev] =
% BoussinesqTS(app,SetupScript,PlotScript)
% Solves the 2D Boussinesq problem using a simple Finite Volume Incremental
% Pressure Correction in Rotational form projection method.
% This function uses a first-order implicit Euler integration scheme, and
% uses LU factorization to solve the matrices. Uses QUICK scheme (mexed) for 
%advection. 
%
%EXAMPLE: 
%>>app.Nx = 30; app.Ny=30; app.T = 20; app.dt = 1/200; app.nu = 0.01;
%>>app.kappaT = 0.01; app.kappaS = 0.01; app.PlotIntrvl = 100;
%>>Thermohaline(app,'ThermohalineSetup','ThermohalinePlots')
%       
%INPUTS:
%       app:	    The application data structured. At a minimum it should 
%		     contain the following fields (which may be added through
%                    the SetupScript).
%          app.Nx:         Number of interior nodes in x-direction (integer number)
%          app.Ny:         Number of interior nodes in y-direction (integer number)
%          app.T:          Final time
%          app.dt:         Timestep size
%          app.nu:         Viscosity
%          app.kappaT:     Thermal conductivity (Optional)
%          app.kappaS:     Saline conductivity
%          app.PlotIntrvl: Number of time integration intervals between of calling
%                           "plotscript."
%       PlotScript: (String input) of filename used to plot specific problem.
%                   Can also be used to save outputs and perform other tasks.
%                       Note: Variables u, v, and P, will be available for
%                       plotting.
%       SetupScript: String input of filename used to setup the problem.
%           This file needs to correctly create the following variables:
%	    Nt:     Number of timesteps (integer number)
%           dx:     Size of control volume in x-direction
%           dy:     Size of control volume in y-direction
%           app.dt:     Size of timestep (TotalTime/Nt)
%           bcsDT:  Temperature Dirichlet boundary conditions
%           bcsDS:  Salinity Dirichlet boundary conditions
%           bcsDP:  Pressure Dirichlet boundary conditions
%           bcsDu:  U-velocity Dirichlet boundary conditions
%           bcsDv:  V-velocity Dirichlet boundary conditions
%           bcsNT:  Temperature Neumann boundary conditions
%           bcsNS:  Salinity Neumann boundary conditions
%           bcsNP:  Pressure Neumann boundary conditions
%           bcsNu:  U-velocity Neumann boundary conditions
%           bcsNv:  V-velocity Neumann boundary conditions
%           NodeT:  Temperature node-numbering matrix (see NodePad.m)
%           NodeS:  Salinity node-numbering matrix (see NodePad.m)
%           NodeP:  Pressure node-numbering matrix (see NodePad.m)
%           Nodeu:  U-velocity node-numbering matrix (see NodePad.m)
%           Nodev:  V-velocity node-numbering matrix (see NodePad.m)
%           idsP:   Interior Pressure node id's (see NodePad.m)
%           idsu:   Interior U-velocity node id's (see NodePad.m)
%           idsv:   Interior V-velocity node id's (see NodePad.m)
%           FT(time,p):   Forcing function for Tsw
%                         p =[ XP(idsP),YP(idsP) ]
%           FS(time,p):   Forcing function for Ssw
%                         p =[ XP(idsP),YP(idsP) ]
%           Fu(time,p):   Forcing function for U-velocity
%                         p =[ XU(idsu),YU(idsu) ]
%           Fv(time,p):   Forcing function for V-velocity
%                         p =[ XV(idsv),YV(idsv) ]
%           XP:     Matrix of x-coordinates corresponding P, Tsw & Ssw CV
%                    centers. XP(idsP) used as 'x' input to FT & FS. (see
%                    meshgrid)
%           XU:     Matrix of x-coordinates corresponding to U-velocity CV
%                    centers. XU(idsu) used as 'x' input to Fu. (see meshgrid)
%           XV:     Matrix of x-coordinates corresponding to V-velocity CV
%                    centers. XV(idsv) used as 'x' input to Fv. (see meshgrid)
%           YP:     Matrix of y-coordinates corresponding to P, Tsw & Ssw CV
%                    centers. YP(idsP) used as 'y' input to FT & FS. (see
%                    meshgrid)
%           YU:     Matrix of y-coordinates corresponding to U-velocity CV
%                    centers. YU(idsu) used as 'y' input to Fu. (see meshgrid)
%           YV:     Matrix of y-coordinates corresponding to V-velocity CV
%                    centers. YV(idsv) used as 'y' input to Fv. (see meshgrid)
%           Tsw:    Initial Temperature vector. The first #Nbc terms contain
%                   the boundary condition values.
%           Ssw:      Initial Salinity vector. The first #Nbc terms contain
%                   the boundary condition values.
%           P:      Initial Pressure vector. The first #Nbc terms contain
%                   the boundary condition values.
%           u:      Initial U-velocity vector. The first #Nbc terms contain
%                   the boundary condition values.
%           v:      Initial V-velocity vector. The first #Nbc terms contain
%                   the boundary condition values.
%
%OUTPUTS:
%           P:      Final Pressure vector, including the boundary
%                   conditions.
%           u:      Final U-velocity vector, including the boundary
%                   conditions.
%           v:      Final V-velocity vector, including the boundary
%                   conditions.
%           rho:    Final density vector, including the boundary
%                   conditions.
%           Tsw:    Final Temperature vector, including the boundary
%                   conditions.
%           Ssw:      Final Salinity vector, including the boundary
%                   conditions.
%           NodeP:  Pressure node-numbering matrix (see NodePad.m)
%           Nodeu:  U-velocity node-numbering matrix (see NodePad.m)
%           Nodev:  V-velocity node-numbering matrix (see NodePad.m)
%
%
%
%   Uses Matlab functions: symamd, lu
%   See also: mk_DiffusionOperator, FixDiffuvDbcs, Grad2D, Div2D
%
%
%ADAPTED by Theo Touvet in July 2010 (to introduce Tsw and Ssw)
%FROM the "NavierStokesInc.m" file by Matt Ueckermann, Pierre Lermusiaux, 
%Pat Haley for MIT course 2.29

                    
%% Run setup script
eval(SetupScript)
if ~exist('bcsDT','var'),bcsDT=zeros(size(bcsDP));end
if ~exist('bcsDS','var'),bcsDS=zeros(size(bcsDP));end
if ~exist('bcsNT','var'),bcsNT=zeros(size(bcsNP));end
if ~exist('bcsNS','var'),bcsNS=zeros(size(bcsNP));end
if ~exist('Tsw','var'),Tsw=zeros(size(P));end
if ~exist('Ssw','var'),Ssw=zeros(size(P));end
if ~exist('FT','var'),FT=@(time,p) 0;end
if ~exist('FS','var'),FS=@(time,p) 0;end
if ~isfield(app,'kappaT'),app.kappaT=0;end
if ~isfield(app,'kappaS'),app.kappaS=0;end
if ~exist('bcsOP','var'),bcsOP=[];end
if ~exist('bcsOu','var'),bcsOu=[];end
if ~exist('bcsOv','var'),bcsOv=[];end

%% Find number of boundary conditions
NbcsP = length(bcsDP) + length(bcsNP);
Nbcsu = length(bcsDu) + length(bcsNu);
Nbcsv = length(bcsDv) + length(bcsNv);
if Nbcsu ~= Nbcsv
    error('For this code we require the same number of boundary conditions on u and v');
end
NbcsT = length(bcsDT) + length(bcsNT);
NbcsS = length(bcsDS) + length(bcsNS);

%% Constructing A and RHS
[AT rhsTbc] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDT,bcsNT, NodeP,'P');
[AS rhsSbc] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDS,bcsNS, NodeP,'P');
[AP rhsPbc] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDP,bcsNP, NodeP,'P');
[Au rhsubc] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDu,bcsNu, Nodeu,'u');
[Av rhsvbc] =  mk_DiffusionOperator(app.Nx, app.Ny, dx, dy, bcsDv,bcsNv, Nodev,'v');
[AP rhsPbc] = FixDiffPObcs(AP,bcsOP,NbcsP,...
    rhsPbc,NodeP,dx,dy);bcsDP(bcsOP)=[];
[Au rhsubc] = FixDiffPObcs(Au,bcsOu,Nbcsu,...
    rhsubc,Nodeu,dx,dy);bcsDu(bcsOu)=[];
[Av rhsvbc] = FixDiffPObcs(Av,bcsOv,Nbcsv,...
    rhsvbc,Nodev,dx,dy);bcsDv(bcsOv)=[];
[Au Av rhsubc rhsvbc] = FixDiffuvDbcs(Au,Av,bcsDu,bcsDv,...
    rhsubc,rhsvbc,Nodeu, Nodev, dx,dy,Nbcsu);

AT = app.kappaT*AT + 1/app.dt*speye(size(AT));
AS = app.kappaS*AS + 1/app.dt*speye(size(AS));
Au = app.nu*Au + 1/app.dt*speye(size(Au));
Av = app.nu*Av + 1/app.dt*speye(size(Av));


%% LU factorization of matrices
imposemean = false;
if isempty(bcsDP) 
    imposemean = true;
end
%First order matrix to give best spartsity pattern
Tperm = symamd(AT);
Sperm = symamd(AS);
Pperm = symamd(AP);
uperm = symamd(Au);
vperm = symamd(Av);

% % LU factorization
[LT UT] = lu(AT(Tperm,Tperm));
[LS US] = lu(AS(Sperm,Sperm));
[LP UP] = lu(AP(Pperm,Pperm));
[Lu Uu] = lu(Au(uperm,uperm));
[Lv Uv] = lu(Av(vperm,vperm));

% Get interior ids
pid = NbcsP+1:length(P);
uid = Nbcsu+1:length(u);
vid = Nbcsv+1:length(v);

srcid     = 2:app.Nx+1;
srcTopid  = 2:app.Ny;
srcBotid  = 3:app.Ny+1;

Neu = [length(bcsDu)+1,Nbcsu];
Nev = [length(bcsDv)+1,Nbcsv];

NbcsDu = length(bcsDu); %Needed for advection
NbcsDv = length(bcsDv); %Needed for advection

%% Solving system
q = zeros(size(P));
T=0; k=0; eval(PlotScript);drawnow
display('Paused at initial condition. Hit any key to continue');drawnow;pause
for k=1:Nt
    T=app.dt*k;
    [dpdx dpdy]   = Grad2D(P,NodeP,dx,dy,idsu,idsv);
    % Calculate RHS vectors for u and v
%     [fu fv] = UVadvect_DO_mex(app.Nx, app.Ny, dx, dy, u, v, u, v,Nodeu, Nodev, idsu, idsv, Nbcsu, NbcsDu, NbcsDv);
    [fu fv] = UVadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, u, v, u, v,Nodeu, Nodev, idsu, idsv, Nbcs, NbcsDu, NbcsDv);

    % NON-LINEAR equation of state using SW_DENS0: density of sea water AT
    % ATMOSPHERIC PRESSURE
    rho_0=1028; g=9.81;
    rho = g*(sw_dens0(Ssw,Tsw)-rho_0*ones(size(P)))/rho_0; %perturbation density divided by rho_0 and multiplied by g. Here Ssw has to be in psu and Tsw in Celsius degree!

    src = 0.5*(rho(NodeP(srcTopid,srcid)) + rho(NodeP(srcBotid,srcid)) );
    src = src(idsv);
    RHSu = 1/app.dt*u(uid) + fu + rhsubc*app.nu - dpdx + ...
        Fu(T,[XU(idsu) YU(idsu)]);
    RHSv = 1/app.dt*v(vid) + fv + rhsvbc*app.nu - dpdy -src+...
        Fv(T,[XV(idsv) YV(idsv)]);
    u(uid(uperm))=Uu\(Lu\RHSu(uperm));
    v(vid(vperm))=Uv\(Lv\RHSv(vperm));

    %Calculate RHS vector for P
    divU = Div2D(u,v,Nodeu,Nodev,dx,dy,idsP,Neu,Nev);
    RHSq = rhsPbc - 1/app.dt*divU;
    q(pid(Pperm))=UP\(LP\RHSq(Pperm));

    %Correct Velocity
    [dqdx dqdy]   = Grad2D(q,NodeP,dx,dy,idsu,idsv);
    u(uid) = u(uid) - app.dt*dqdx;
    v(vid) = v(vid) - app.dt*dqdy;

    %correct pressure
    P(pid) = P(pid) + q(pid) - app.nu*divU;
    %remove mean pressure
    if imposemean 
        P(pid)=P(pid)-mean(P(pid));        
    else
        P([1:length(bcsDP),pid])=...
            P([1:length(bcsDP),pid])-mean(P(pid));        
    end
    %Advect Tracers (re-use RHSq to save a little memory)
    RHSq=1/app.dt*Tsw(pid) + app.kappaT*rhsTbc...
        +SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, Tsw, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,length(bcsDT),NbcsDu,NbcsDv)+...        
        FT(T + app.dt,[XP(idsP) YP(idsP)]);
        %         +SCAadvect_mex(app.Nx, app.Ny, dx, dy, Tsw, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,length(bcsDT),NbcsDu,NbcsDv)+...
    Tsw(pid(Tperm)) = UT\(LT\RHSq(Tperm));
    
    RHSq=1/app.dt*Ssw(pid) + app.kappaS*rhsSbc+...
        SCAadvect_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, Ssw, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,length(bcsDS),NbcsDu,NbcsDv)+...
        FS(T + app.dt,[XP(idsP) YP(idsP)]);
%         SCAadvect_mex(app.Nx, app.Ny, dx, dy, Ssw, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,length(bcsDS),NbcsDu,NbcsDv)+...        
    Ssw(pid(Sperm)) = US\(LS\RHSq(Sperm));

    %Plot solution
    if ~mod(k,app.PlotIntrvl)
        eval(PlotScript); drawnow
    end
end
