d0 = fileparts([pwd, filesep]); 

d0 = d0(1 : strfind(d0,'trunk')+4); 

if ~isfield(app, 'savemat'), app.savemat = 'Save/mat/Bottom_Gravity_Current'; end; 

if ~exist(sprintf('%s/%s', d0, app.savemat), 'dir') 

    mkdir(sprintf('%s/%s', d0, app.savemat)); 

end; 

if ~isfield(app, 'NoPlot'), app.NoPlot = false; end; 

if ~app.NoPlot 

    if ~isfield(app, 'saveimg'), app.saveimg = 'Save/img/Bottom_Gravity_Current'; end;

    if ~exist(sprintf('%s/%s', d0, app.saveimg), 'dir')

        mkdir(sprintf('%s/%s', d0, app.saveimg));

    end;

end;

%keyboard

%% Set up primitive parameters of the bottom gravity current.

if isfield (app, 'pass_params') 

   A = app.pass_params.A; % Horizontal size of solution domain (unit: m) 

   B = app.pass_params.B; % Vertical size of solution domain (unit: m) 

   A1 = app.pass_params.A1; % Length of the initial plateau (unit: m) 

   B1 = app.pass_params.B1; % Height of the initial plateau (unit: m) 

   A2 = app.pass_params.A2; % Length of salinity forcing region (unit: m) 

   h = app.pass_params.h; % Length scale of thickness of grivity current (unit: m) 

%    (The paper does not directly provide this value.)

   rhoinitflag = app.pass_params.rhoinitflag; 

else

   A = 10e3; % Horizontal size of solution domain (unit: m)

   B = 1e3; % Vertical size of solution domain (unit: m)

   A1 = 1e3; % Length of the initial plateau (unit: m)

   B1 = 0.8e3; % Height of the initial plateau (unit: m)

   A2 = 0.6e3; % Length of salinity forcing region (unit: m)

   h = 500; % Length scale of thickness of grivity current (unit: m)

%    (The paper does not directly provide this value.)

   rhoinitflag = 'cosine';

end;

% Maximum variation of salinity (unit: psu(nondimensional))

if isfield(app, 'DelS'), DelS = app.DelS; else DelS = 1.5; end; 

% Slope angle of the inclined plane (unit: degree, not radian!)

if isfield(app, 'theta'), theta = app.theta; else theta = 2.5; end; 

nu_x = 7; % Eddy viscosity in horizontal direction (unit: m^2/s) 

r = 1e-5; % Ratio of vertical to horizontal diffusivity (nu_z / nu_x) 

Pr = 7; % Turbulent Prandtl number (nu_x / K_x) 

k = 1/21600; % Salinity forcing relaxation frequency (unit: Hz) 

Ra = 506 * DelS; % Rayleigh number 

kappa = k*h^2/nu_x; % Nondimensional relaxation frequency 

%% Set up physical properties and forcing terms

% Define momentum forcing terms

% 'XY' is a 2-D array.

% XY(:, 1) contains x coordinates, XY(:, 2) contains y coordinates.

% 't' is just the current time (initial time is t = 0).

% Fu and Fv are required functions that need to be set for the solver

Fu = @(t, XY) 0; 

Fv = @(t, XY) 0; 

if ~isfield(app, 'nu'), app.nu = 1; end; 

if ~isfield(app, 'nu2'), app.nu2 = r; end; 

if ~isfield(app, 'kappa'), app.kappa = 1/Pr; end; 

if ~isfield(app, 'kappa2'), app.kappa2 = r/Pr; end; 

if ~isfield(app, 'g'), app.g = Ra; end; 

%% Set up mesh and solution domain

% Set size of domain

a = A/h; 

b = B/h; 

a1 = A1/h; 

b1 = B1/h; 

a2 = A2/h; 

tan_theta = tan(theta*pi/180); 

if ~isfield(app, 'Nx'), app.Nx = 2000; end; 

if ~isfield(app, 'Ny'), app.Ny = 200; end; 

% Set number of timesteps.

if ~isfield(app, 'T'), app.T = 16200*nu_x/h^2; end; 

if ~isfield(app, 'dt'), app.dt = 0.1*nu_x/h^2; end; 

Nt = ceil(app.T/app.dt); 

if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = round((app.T/app.dt)/270); end; 

% Mesh the domain

% This is a standard call and should not require much modification.

% [0, a] gives the min/max x coordinates, and [0, b] gives the min/max y coordinates.

[XP, YP, XU, YU, XV, YV, dx, dy] = meshdomain(app, [0, a], [0, b]); 

app.IP = 1./(dx*dy*(app.sigma_nd.*app.sigma_nd)); 

 % want norm to be 1, i.e. (dx*dy*var_nd_rho*rhoi*rhoi)=1

 % rescale = sqrt(app.IP(3))

%% Set up the geometry mask

% This mask is user-created. Zeros are unmasked, and non-zero integers

% indicate a masked region. The domain sides should not be considered here.

% Each interior masked region could have a different number to allow

% different boundary conditions.

% Create mask of correct size. This is a standard call and should not

% require modification.

Mask = zeros(app.Ny, app.Nx); 

% User-defined true-false (boolean) functions that define the mask.

% True if masked, false if not masked.

forcing = @(x, y) (x <= a2) & (y >= b1); 

plateau = @(x, y) (x <= a1) & (y <= b1); 

slope = @(x, y) (x > a1) & (y <= max(0, (a1-x)*tan_theta + b1)); 

% Assign masked region.

Mask(plateau(XP, YP)) = 1; 

Mask(slope(XP, YP)) = 1; 

% Logical array add_bnd signifies whether to add (add if =true) boundaries for the 

% [west, east, south, north] sides of the solution domain.

% This is almost a standard procedure except when periodic boundaries are used.

add_bnd = [true, true, true, true]; 

% Create preliminary node numbering matrices using function "NodePad".

% Note: This is a standard procedure and should not be modified.

[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask, add_bnd); 

% Plot the current boundary IDs to facilitate setting up boundary conditions.

if ~app.NoPlot 

    plot_bc_ids(NodeP, Mask);

    fprintf('Paused at boundary id plot. Hit any key to continue.\n\n');

    pause;

end;

%% Set up boundary conditions (BCs)

% User assigns BC types to each boundary ID for u, v, rho in three cell arrays.

% This requires user input. Use the above plot as reference to identify BC IDs.

% e.g. Element 1 in cell array bcid2type_u represents the BC type of the boundary

% with ID 1 for velocity u.

% Note: Only the following four types of BCs can be set for velocities and density:

% Single-value Dirichlet BC   (='D')

% Multiple-value Dirichlet BC (='Dmv')

% Single-value Neumann BC     (='N')

% Multiple-value Neumann BC   (='Nmv')

% Warning: Open boundary conditions (='O') are not allowed for velocities and density!

bcid2type_u = {'D', 'N', 'N', 'D', 'N'};       

bcid2type_v = {'D', 'N', 'N', 'D', 'D'}; 

bcid2type_rho = {'N', 'N', 'N', 'N', 'N'}; 

% Automatically assign BC types to each boundary ID for q which is an intermediate variable

% for which we need to solve a Poisson equation in the projection methods

% used by the Navier-Stokes solver.

% The BC types of q depends on those of velocities and can be either 'N' or 'O'.

% Note: This is a standard procedure and should not be modified.

if ~isfield (app,'qBC') 

   app.qBC='Default'; 

end; 

Nbcs = size(bcid2type_u); Nbcs = Nbcs(2); 

bcid2type_q = cell(1, Nbcs); 

for i = 1 : Nbcs 

    if strcmp(bcid2type_u(i), 'D') || strcmp(bcid2type_u(i), 'Dmv')... 

            || strcmp(bcid2type_v(i), 'D') || strcmp(bcid2type_v(i), 'Dmv')

        bcid2type_q{i} = 'N'; 

    else 

        bcid2type_q{i} = 'O'; 

    end; 

end; 

% Allow over-ride of automatically defined q BC

switch(app.qBC) 

   case 'Dinlet' 

      bcid2type_q{1} = 'D';

   case 'Dinout' 

      bcid2type_q{1} = 'D';

      bcid2type_q{2} = 'D';

   case 'Dirichlet' 

      for i=1:Nbcs

         bcid2type_q{i} = 'D';

      end;

end;

% Set up the final node numbering matrix using function "set_bcs".

% For each single-value BC ('D', 'N', 'O'), all its boundary cells share one ID,

% like in the preliminary node numbering matrices.

% For each multiple-value BC ('Dmv', 'Nmv'), now each boundary cell will be assigned

% a different ID.

% Besides, the IDs are arranged to satisfy: 

%               Dirichlet_BC_IDs < Open_BC_IDs < Neumann_BC_IDs

% Note: This is a standard procedure and should not be modified.

% Without density equation:

% [NodeP, nbcP, Nodeu, nbcu, Nodev, nbcv] = set_bcs(...

%     bcid2type_P, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev);

% With density equation:

[NodeP, nbcq, Nodeu, nbcu, Nodev, nbcv, NodeRho, nbcrho] = set_bcs(...  

    bcid2type_q, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev, bcid2type_rho, NodeP);

% Set preliminary BC values to zero since most BCs are homogeneous.

% Some of them can be modified to be non-homogeneous BCs later.

% Note: This is a standard procedure and should not be modified.

[bcsDq, bcsOq, bcsNq] = init_bcs(nbcq); 

[bcsDu, bcsOu, bcsNu] = init_bcs(nbcu); 

[bcsDv, bcsOv, bcsNv] = init_bcs(nbcv); 

[bcsDrho, bcsOrho, bcsNrho] = init_bcs(nbcrho); 

% User sets non-zero boundary conditions.

% The coordinates of the boundary points on each grid can be obtained...

% by function "get_bc_coord".

% Each of the three XY arrays contains 2 columns with x coordinates in the first and 

% y coordinates in the second stored in order of BC IDs.

% The three XY arrays correspond to boundary points of type Dirichlet ('D' and 'Dmv'),

% open ('O'), and Neumann ('N' and 'Nmv').

% [XYD, XYO, XYN] = get_bc_coord('u', XU, YU, Nodeu, bcsDu, bcsOu, bcsNu);

% ids = (abs(XYD(:, 1)) < eps);

% bcsDu(ids) = 1*h/nu_x;

% [XYD, XYO, XYN] = get_bc_coord('P', XP, YP, NodeRho, bcsDrho, bcsOrho, bcsNrho);

% ids = (abs(XYD(:, 1)) < eps);

% bcsDrho(ids) = 0.5*(1 - cos(pi*(b-XYD(ids, 2)/(b-b1))));

%% Set up initial conditions

% User specifies the functions that provide initial values at given coordinates.

[rhoinit_choice,rhoinit_parms] = strtok (rhoinitflag); 

switch rhoinit_choice 

   case 'cosine' 

      rho_init = @(X, Y) 0.5*(1 - cos(pi*(b-Y)/(b-b1))).*forcing(X,Y);

      P_init   = @(X, Y) (app.g*0.5)*((b-Y) - ((b-b1)/pi).*sin(pi*(b-Y)/(b-b1))).*forcing(X,Y);

   case 'linear' 

      rho_init = @(X, Y) (b-Y)./(b-b1).*forcing(X,Y);

      P_init   = @(X, Y) app.g*((b-Y).^2)./(2.*(b-b1)).*forcing(X,Y);

   case 'exponential' 

      rhodk = str2num(rhoinit_parms); 

      rho_init = @(X, Y) (1-exp(rhodk.*(Y-b)))./(1-exp(rhodk.*(b1-b))).*forcing(X,Y); 

      P_init   = @(X, Y) app.g*((b-Y)-(1-exp(rhodk.*(Y-b)))./rhodk)./(1-exp(rhodk.*(b1-b))).*forcing(X,Y); 

   case 'cosine+exponential'

      rhodk = str2num(rhoinit_parms);

      rho_init = @(X, Y) 0.5.*( 0.5*(1 - cos(pi*(b-Y)/(b-b1))) + (1-exp(rhodk.*(Y-b)))./(1-exp(rhodk.*(b1-b))) ).*forcing(X,Y);

      P_init   = @(X, Y) (app.g*0.5)*( 0.5*((b-Y) - ((b-b1)/pi).*sin(pi*(b-Y)/(b-b1))) + ((b-Y)-(1-exp(rhodk.*(Y-b)))./rhodk)./(1-exp(rhodk.*(b1-b))) ).*forcing(X,Y);

   case 'cosine+linear'

      rho_init = @(X, Y) 0.5.*( 0.5*(1 - cos(pi*(b-Y)/(b-b1))) + (b-Y)./(b-b1) ).*forcing(X,Y);

      P_init   = @(X, Y) (app.g*0.5)*( 0.5*((b-Y) - ((b-b1)/pi).*sin(pi*(b-Y)/(b-b1))) + ((b-Y).^2)./(2.*(b-b1)) ).*forcing(X,Y);

   case 'linear+exponential'

      rhodk = str2num(rhoinit_parms);

      rho_init = @(X, Y) 0.5.*( (b-Y)./(b-b1) + (1-exp(rhodk.*(Y-b)))./(1-exp(rhodk.*(b1-b))) ).*forcing(X,Y);

      P_init   = @(X, Y) (app.g*0.5)*( ((b-Y).^2)./(2.*(b-b1)) + ((b-Y)-(1-exp(rhodk.*(Y-b)))./rhodk)./(1-exp(rhodk.*(b1-b))) ).*forcing(X,Y);

end;

u_init   = @(X, Y) zeros(size(X)); 

v_init   = @(X, Y) zeros(size(X)); 

% Initialize the state variables with values from both boundary and internal cells.

% Note: This is a standard procedure and should not be modified.

rho = [bcsDrho; bcsNrho; rho_init(XP(idsP), YP(idsP))]; 

P   = [bcsDq; bcsNq; P_init(XP(idsP), YP(idsP))]; 

u   = [bcsDu; bcsNu; u_init(XU(idsu), YU(idsu))]; 

v   = [bcsDv; bcsNv; v_init(XV(idsv), YV(idsv))]; 

%keyboard

%% Set up density forcing term

% Define density forcing term

if ~isfield (app,'rhofrc_type') 

   app.rhofrc_type = 'fullvert'; 

end; 

switch app.rhofrc_type 

   case 'fullvert' 

      app.Frho = @(t, XY, rho) -kappa*(rho - 1).*((XY(:,1) <= a1)&(XY(:,2) >= b1)); 

   case 'halfvert'

      app.Frho = @(t, XY, rho) -kappa*(rho - (XY(:,2)<=(b+b1)*0.5)).*((XY(:,1) <= a1)&(XY(:,2) >= b1));

end;

%% Save all the setup data in a mat-file.

save(sprintf('%s/%s/setup', d0, app.savemat));