% Thermohaline circulation setup script
% WRITTEN by Theo Touvet in July 2010

%Forcing Function
% presently NONE
FT = @(time,p) 0;
FS = @(time,p) 0;
Fu = @(time,p) 0;
Fv = @(time,p) 0;

%Domain definition and discretisation
a      = 2;          b      = 1;
dx     = a./(app.Nx);    dy     = b./(app.Ny); 

%Boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Change this:
%bcsDT= [0 0 0 zeros(1,app.Nx)];  bcsNT= [];
bcsDT= [(0) (0) (0) (10+10*cos(pi*linspace(-a/2 + dx/2,a/2 - dx/2,app.Nx)/(a/2)))]';  bcsNT= [];   % bcsDT has to be in Celsius degree with the non linear equation of state used!
bcsDS= [];  bcsNS= [0 0 0 0.05*cos(pi*linspace(-a/2 + dx/2,a/2 - dx/2,app.Nx)/(a/2))]';                     % bcsDS has to be in       psu      with the non linear equation of state used!
% bcsDT= [(273+4-10) (273+4-10) (273+4-10) ((273+4)+10*cos(pi*linspace(-a/2 + dx/2,a/2 - dx/2,app.Nx)/(a/2)))];  bcsNT= [];
% bcsDS= [];  bcsNS= [0 0 0 0.5*cos(pi*linspace(-a/2 + dx/2,a/2 - dx/2,app.Nx)/(a/2))];
%%%%% ### eg: bcsDS= [];  bcsNS= [0.1 -0.1 0 0];% ~ Salinity flux: positive input at the equator (x=0) due to evaporation, release of fresh water at the pole (x=1) due to iceberg melting!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bcsDP  = [];         bcsNP  = [0 0 0 zeros(1,app.Nx)]';        
bcsDu  = [0 0 0 0]';  bcsNu  = [];        
bcsDv  = [0 0 0 0]';  bcsNv  = [];        

%Set the number of timesteps.
Nt = ceil(app.T/app.dt);

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

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

%% Change pressure node number so that a non uniform Temperature AND Salinity boundary condition can be added on top of the domain.
%For instance, this gives the following result for NodeP with app.Nx=5 and app.Ny=3:
% 4  4  4  4  4  4  4          1  4  5  6  7  8  2
% 1  5  8 11 14 17  2          1  9 12 15 18 21  2
% 1  6  9 12 15 18  2    =>    1 10 13 16 19 22  2
% 1  7 10 13 16 19  2          1 11 14 17 20 23  2
% 3  3  3  3  3  3  3          1  3  3  3  3  3  2

%first, re-define the corners (so that the function applied on top does not
%jut out from  the actual domain of length app.Nx):
NodeP(1,1)=1;       NodeP(1,end)=2;
NodeP(end,1)=1;     NodeP(end,end)=2; %this line is not very useful here, but just to be coherent with the former line...

%then, increment every "interior" number of NodeP by the number of unique
%boundary entries that I am adding. In this case (top BCs only): (app.Nx-1). %% app.Nx-1=(size(NodeP,2)-2)-1 = one for every column minus 2 to keep the corner unchanged and minus 1
NodeP(2:(end-1),(2:end-1)) = NodeP(2:(end-1),(2:end-1)) + app.Nx - 1;

%finally, add the unique node numbers. That is in the example above allocate 4, 5,
%6, 7, 8 to NodeP(1,2)...NodeP(1,6).
NodeP(1,2:(end-1)) = 4:(4+(app.Nx-1));

%As far as I can tell, there is no need for the moment to repeat this for
%other node matrices (ie just let Nodeu and Nodev unchanged).

%% Initialize fields
xp = linspace(-a/2 + dx/2,a/2 - dx/2,app.Nx);yp = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xu = linspace(-a/2 + dx,a/2 - dx,app.Nx-1);  yu = fliplr(linspace(dy/2, b-dy/2, app.Ny  ));
xv = linspace(-a/2 + dx/2,a/2 - 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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Change this:
Tswinit=@(X,Y) (Y>-1)*(4);%Average Temperature in the ocean. %Has to be in Celsius degree with the non linear equation of state used!
Sinit=@(X,Y) (Y>-1)*34.7;%Average Salinity in the ocean. %Has to be in       psu      with the non linear equation of state used!
%%%%% ### eg: Tinit=@(X,Y) ((X>0.5) - (X<0.5))*(-10); %% ~ Thermaly driven circulation from equator to North pole at the beginning with warm pool at the equator region and cool water at th Pole!!!!!!!!
%%%%% ### eg: Sinit=@(X,Y) ((X>0.5) - (X<0.5))*(-1); %% ~ Salinity distribution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pinit = @(X,Y) zeros(size(X));
uinit = @(X,Y) zeros(size(X));
vinit = @(X,Y) zeros(size(X));
% vinit = @(X,Y) rhoinit(X,Y);
Tsw = [bcsDT;bcsNT;  Tswinit(XP(idsP), YP(idsP))];
Ssw = [bcsDS;bcsNS;  Sinit(XP(idsP), YP(idsP))];
P = [bcsDP;bcsNP;  pinit(XP(idsP), YP(idsP))];
u = [bcsDu;bcsNu;  uinit(XU(idsu), YU(idsu))];
v = [bcsDv;bcsNv;  vinit(XV(idsv), YV(idsv))];
