Bottom_Gravity_Current_Setup6 

% Initialize DO modes.

% extract parameters from app structure

S = app.S ; MC = app.MC ; 

if isfield (app,'DO_IC_flag') 

   DO_IC_flag = app.DO_IC_flag; 

else

   DO_IC_flag = 'SineCosine_Plateau';

end;

if ~isfield (app,'MollifierWidth') 

   app.MollifierWidth = []; 

end; 

ui = zeros(length(u),app.S(1)) ; 

vi = zeros(length(v),app.S(1)) ; 

Pi = zeros(length(P),app.S(1)) ; 

rhoi = zeros(length(rho),app.S(1)) ; 

var = ones(app.S(1),1) ; 

switch DO_IC_flag 

   case {'SineCosine_Plateau' 'SineCosine_Plateau_Sic'}; 

   % Initialize with Sines & Cosines over plateau or only over

   % plateau region with nonzero initial S.

      if strcmp (DO_IC_flag, 'SineCosine_Plateau_Sic')

         a_do_ic = a2;

      else

         a_do_ic = a1;

      end;

      b_do_ic = b-b1;

      initfn1 = @(x,y,M,N) (-(1/b_do_ic)*sin(pi.*x/a_do_ic).*sin(pi.*M.*x/a_do_ic) ...

       .*( cos(pi.*(b-y)/b_do_ic).*sin(pi.*N.*(b-y)/b_do_ic) ...

           + N.*sin(pi.*(b-y)/b_do_ic).*cos(pi.*N.*(b-y)/b_do_ic) )) .* ((x<=a_do_ic) & (y>=b1)) ;

      initfn2 = @(x,y,M,N)  ((1/a_do_ic)*sin(pi.*(b-y)/b_do_ic).*sin(pi.*N.*(b-y)/b_do_ic) ...

       .*( cos(pi.*x/a_do_ic).*sin(pi.*M.*x/a_do_ic) ...

           + M.*sin(pi.*x/a_do_ic).*cos(pi.*M.*x/a_do_ic) )) .* ((x<=a_do_ic) & (y>=b1)) ;

      i = 1 ;

      i_rho = 1;

      for ii = 1:app.S(1)

          for iiy = 1:ceil(ii/2)

              if ( i <= app.S(1) )

                  ix = ii - iiy + 1 ;

                  iy = iiy ;

                  nrm = sqrt(((ix^2 + 1))/16*(3/2)^(iy==1) ...

                        + ((iy^2 + 1))/16*(3/2)^(ix==1)) ;

                  ui(Nbcs+1:length(u),i) = initfn1(XU(idsu), YU(idsu), ix, iy ) / (nrm*sqrt(app.IP(1)));

                  vi(Nbcs+1:length(v),i) = initfn2(XV(idsv), YV(idsv), ix, iy ) / (nrm*sqrt(app.IP(2)));

                  var(i) = 10 * exp((1-i)/2.5) ;

                  i = i + 1;

                  if (i_rho <=app.S(1))

                     rhoi(Nbcs+1:length(rho),i_rho) = initfn1(XP(idsP),YP(idsP),ix,iy)./(nrm.*sqrt(app.IP(3)));

                     i_rho = i_rho + 1;

                  end;

                  if ( i_rho <= app.S(1) )

                     rhoi(Nbcs+1:length(rho),i_rho) = initfn2(XP(idsP),YP(idsP),ix,iy)./(nrm*sqrt(app.IP(3)));

                     i_rho = i_rho + 1;

                  end;

              end

              if ( i <= app.S(1) && ix ~= iy )

                  % ix and iy interchanged

                  iy = ii - iiy + 1 ;

                  ix = iiy ;

                  nrm = sqrt(((ix^2 + 1))/16*(3/2)^(iy==1) ...

                        + ((iy^2 + 1))/16*(3/2)^(ix==1)) ;

                  ui(Nbcs+1:length(u),i) = initfn1(XU(idsu), YU(idsu), ix, iy ) / (nrm*sqrt(app.IP(1)));

                  vi(Nbcs+1:length(v),i) = initfn2(XV(idsv), YV(idsv), ix, iy ) / (nrm*sqrt(app.IP(2)));

                  var(i) = 10 * exp((1-i)/2.5) ;

                  i = i + 1;

                  if (i_rho <=app.S(1))

                     rhoi(Nbcs+1:length(rho),i_rho) = initfn1(XP(idsP),YP(idsP),ix,iy)./(nrm*sqrt(app.IP(3)));

                     i_rho = i_rho + 1;

                  end;

                  if ( i_rho <= app.S(1) )

                     rhoi(Nbcs+1:length(rho),i_rho) = initfn2(XP(idsP),YP(idsP),ix,iy)./(nrm*sqrt(app.IP(3)));

                     i_rho = i_rho + 1;

                  end;

              end

          end

      end

      % orthonormalize DO modes

      uid = Nbcs+1:length(u);

      vid = Nbcs+1:length(v);

      pid = Nbcs+1:length(rho);

      for i = 1:S

          for j = 1:S

              if ( i ~= j )

                  coef = app.IP(1) * InProd(ui(:,i),ui(:,j),dx,dy,uid) ...

                         + app.IP(2) * InProd(vi(:,i),vi(:,j),dx,dy,vid) ...

                         + app.IP(3) * InProd(rhoi(:,i), rhoi(:,i), dx, dy, pid);

                  ui(:,i) = ui(:,i) - coef * ui(:,j) ;

                  vi(:,i) = vi(:,i) - coef * vi(:,j) ;

                  rhoi(:,i) = rhoi(:,i) - coef * rhoi(:,j);

              end

          end

          nrm = sqrt(app.IP(1)*InProd(ui(:,i),ui(:,i),dx,dy,uid) ...

                     + app.IP(2)*InProd(vi(:,i),vi(:,i),dx,dy,vid) ...

                     + app.IP(3)*InProd(rhoi(:,i), rhoi(:,i), dx, dy, pid)) ;

          ui(:,i) = ui(:,i)/nrm ;

          vi(:,i) = vi(:,i)/nrm ;

          rhoi(:,i) = rhoi(:,i) / nrm;

      end

      % Redimensionalize modes

      ui = ui .*app.varscale./sqrt(app.IP(1));

      vi = vi .*app.varscale./sqrt(app.IP(2));

      rhoi = rhoi .*app.varscale./sqrt(app.IP(3));

      %Initialize stochastic coefficients as normal random with variance as calculated

      %in Stoch_Initial_Perturb

      Yi=zeros(app.MC,app.S(1));

      for i=1:app.S(1)

          Yi(:,i)=var(i).*normrnd(0,1,app.MC,1);

      end

   case {'SVD_ICs'}; 

   % This option creates modes by SVD of frequent output from short run.

   [u,v,rho,ui,vi,rhoi,Yi] = BGC_DO_simsvdIC (app); 

   case {'Correlation_Rho' 'Correlation_Rho_Sic' };

      switch (DO_IC_flag)

         case 'Correlation_Rho_Sic'

            a_limit = a2;

         otherwise

            a_limit = a1;

      end;

%keyboard

      [rhoi, Yi] = Stoch_Initial_Perturb_rho (XP,YP,dx,dy,a_limit,b,b1,app,Nbcs,NodeP,idsP);

   case {'Restart'};

   % This option creates modes by SVD of frequent output from short run.

   [u,v,rho,P,ui,vi,rhoi,Pi,Yi] = BGC_DO_restart (app,u,v,rho,P,ui,vi,rhoi,Pi,idsu,idsv,idsP,Nodeu,Nodev,NodeP);

end;