function phi = fvm_soln_SS_adv_diff( X_edges, Z_edges, X_P, Z_P, Gamma_P, ...
                 rho, S_C, S_P, bc_u, bc_d, bc_e, bc_w )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  find steady-state phi(x,z), for combined advection and diffusion
%%  with given heat sources S and conductivity Gamma_P.
%%
%%      div( rho u phi) = div( K grad(phi) ) + S
%%
%%  Uses Patankar Finite Volume method.
%%
%%  IF upwind = 0,    (in global block)
%%  THIS CODE INCLUDES ADVECTIVE FLUX WITH CENTERED DIFFERENCE
%%            (NO UPWINDING SCHEME)  
%%    *** DON'T USE THIS AT HOME FOR SERIOUS WORK!!  ***
%%
%%  IF upwind = 1,
%%  THIS CODE INCLUDES ADVECTIVE FLUX WITH PATANKAR POWER-LAW SCHEME
%%
%%
%%  EXTERNALLY SPECIFIED VELOCITY MUST SATISFY CONTINUITY,
%%   i.e. IS A CONSTANT IN 1-D
%%
%%  Boundary Conditions:
%%  vector bc_u = phi_u      set value of phi on upper boundary
%%  vector bc_d = phi_d      set value of phi on lower boundary
%%  vector bc_w = dphi_dx_w  set gradient on west boundary
%%  vector bc_e = dphi_dx_e  set gradient on east boundary
%%
%%
%%  THIS CODE APPLIES FLUX BOUNDARY CONDITIONS ON VOLUME INTERFACES
%%  on East and West
%%  APPLIES FIXED VALUE BOUNDARY CONDITIONS ON VOLUME CENTERS
%%  in Up and Down 
%%
%%  a_P = coefficient of phi_P at point P = (i,j)
%%  a_E = coefficient of phi_E at point east of (i,j)
%%  a_E, a_W, a_U, a_D for points East, West, Up, Down
%%
%%  System of equations to be solved for phi
%%
%%       big_A * phi = rhs
%%
%%  big_A = matrix with coefficients a_E etc, modified to include BC
%%  phi   = unknown values in Control-volumes
%%  rhs   = right-hand-side column vector containing sources and BC terms
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
    global upwind  Peclet

%
% get Finite-Volume dimensions
% ------------------------------
       [nz_P, nx_P] = size(X_P);


%  get dX, dZ - spacings between volume edges
%  -------------------------------------------
       dX = diff(X_edges')';
       dZ = diff(Z_edges);
%
%   reduce dimensions in vertical, where boundaries are at nodes
       dX = dX(1:end-1,:);


%  get dX_E, dx_W, dZ_U, dZ_D - 
%          - distances to adjacent volume centers
%  -----------------------------------------------
%
      if( size(X_P,2) == 1)
          dX_e = ones( size(X_P) );
          dX_w = ones( size(X_P) );
      else
          dX_e = diff(X_P')';
       % add dummy row at eastern edge to maintain size [nx, nz]
          dX_e = [ dX_e dX_e(:,end) ];
%
          dX_w = diff(X_P')';
       % add dummy row at western edge to maintain size [nx, nz]
          dX_w = [ dX_w(:,1) dX_w ];
      end
%
       dZ_u = diff(Z_P);
    % add dummy row at top edge to maintain size [nx, nz]
       dZ_u = [ dZ_u(1,:)' dZ_u' ]';
%
       dZ_d = diff(Z_P);
    %  add dummy row at bottom edge to maintain size [nx, nz]
       dZ_d = [ dZ_d' dZ_d(end,:)' ]';

%

%  get f_e, f_w, f_u, f_d 
%          - fractional distances to volume interfaces
%  ---------------------------------------------------
       f_e = 1 - (X_edges(1,2:end) - X_P(1,:))./ dX_e(1,:);
       f_e = repmat( f_e, nz_P, 1 );
        
       f_w = 1 - (X_P(1,:) - X_edges(1,1:(end-1)) ) ./ dX_w(1,:);
       f_w = repmat( f_w, nz_P, 1 );
       
       f_u = 1 - (Z_P(:,1) - Z_edges(1:(end-1),1 ))./ dZ_u(:,1);
       f_u = repmat( f_u, 1, nx_P );
%
       f_d = 1 - (Z_edges(2:end,1 ) - Z_P(:,1) )./ dZ_d(:,1);
       f_d = repmat( f_d, 1, nx_P );
%



% set up b matrix - includes constant part of source matrix S 
%         S = S_C + S_P * phi_P    (units are W m^-3)
% ----------------------------------------------------------
         b = S_C.* dX.* dZ;
%
%
% set up matrices for conductivities in adjacent volumes
% ------------------------------------------------------
    Gamma_E   =  [ Gamma_P(:,2:end) Gamma_P(:, end) ];
    Gamma_W   =  [ Gamma_P(:,1) Gamma_P(:,1:end-1) ];
    Gamma_U   =  [ Gamma_P(1,:)' Gamma_P(1:end-1,:)' ]';
    Gamma_D   =  [ Gamma_P(2:end,:)' Gamma_P(end,:)' ]';
%
%
% set up interface conductivity matrices
% --------------------------------------
    Gamma_e =  1./ ( (1 - f_e)./Gamma_P + f_e./Gamma_E );
    Gamma_w =  1./ ( (1 - f_w)./Gamma_P + f_w./Gamma_W );
    Gamma_u =  1./ ( (1 - f_u)./Gamma_P + f_u./Gamma_U );
    Gamma_d =  1./ ( (1 - f_d)./Gamma_P + f_d./Gamma_D );
%
%
%  advection
%%%%%%%%%%%%%%%%%
%  set up uniform velocity for 1-D
     u_0 = 0;             % m s^-1
     w_0 = 2;            % m s^-1
 %
     u_velo = repmat( u_0, size(X_P) ); 
     w_velo = repmat( w_0, size(X_P) ); 
%

% set up velocities at interfaces
%% (this is quick and dirty because velocity is uniform)
%%  In a real problem, you would call a subroutine that
%%  provided u velocities on the east-west interfaces,
%%  and w velocities on the Up-Down interfaces
% ------------------------------------------------------
    u_e = u_velo;
    u_w = u_velo;
    w_u = w_velo;
    w_d = w_velo;

%
% set up Flow and Diffusion coefficients
% ---------------------------------------
    D_e = (Gamma_e ./ dX_e).* dZ;
    D_w = (Gamma_w ./ dX_w).* dZ; 
    D_u = (Gamma_u ./ dZ_u).* dX;
    D_d = (Gamma_d ./ dZ_d).* dX;
%
    F_e = ( rho.* u_e ).* dZ;
    F_w = ( rho.* u_w ).* dZ; 
    F_u = ( rho.* w_u ).* dX;
    F_d = ( rho.* w_d ).* dX;
%

% set up Peclet numbers
% ----------------------
    P_e = F_e./ D_e;
    P_w = F_w./ D_w;
    P_u = F_u./ D_u;
    P_d = F_d./ D_d;

%
% set up coefficients at adjacent points
% ---------------------------------------

   if(upwind == 0 )
      a_E = ( Gamma_e ./ dX_e - rho.* u_e/2 ).* dZ;
      a_W = ( Gamma_w ./ dX_w + rho.* u_w/2 ).* dZ; 
      a_U = ( Gamma_u ./ dZ_u + rho.* w_u/2 ).* dX;
      a_D = ( Gamma_d ./ dZ_d - rho.* w_d/2 ).* dX;
%
   else
%
     a_E = D_e .* A( P_e ) + F_upwind( -F_e );
     a_W = D_w .* A( P_w ) + F_upwind(  F_w );
     a_U = D_u .* A( P_u ) + F_upwind(  F_u );
     a_D = D_d .* A( P_d ) + F_upwind( -F_d );
 %
   end

         
% set up coefficient at central point P (incl source dependent on phi)
% --------------------------------------------------------------------
    a_P = a_E + a_W + a_U + a_D - S_P.* dX.* dZ; 
%

%
%  set Boundary Conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
     phi_u     = bc_u;
     phi_d     = bc_d;
     dphi_dx_e = bc_e;
     dphi_dx_w = bc_w;

%    specified gradient boundaries
%    ------------------------------
%     East boundary  (dphi/dx = dphi_dx_e)
%           use  dphi/dx = (phi_E - phi_P)/dx  
%          i.e.    phi_E = phi_P + dx * dphi_dx_e
         a_P(:,end) = a_P(:,end) - a_E(:,end); 
         b(:,end)   = b(:,end) + a_E(:,end).* dX_e(:,end).* dphi_dx_e;
         a_E(:,end) = 0;

%     West boundary  (dphi/dx = dphi_dx_w)
%           use  dphi/dx = (phi_P - phi_W)/dx
%          i.e.    phi_W = phi_P - dx * dphi_dx_w
         a_P(:,1) = a_P(:,1) - a_W(:,1);
         b(:,1)   = b(:,1) -  a_W(:,1).* dX_w(:,1).* dphi_dx_w ;
         a_W(:,1) = 0;
%
    
%
%    specified phi = phi_u on upper boundary
%    -------------------------------------------
         a_P(1,:) = 1;
         b(1,:)   = phi_u;
         a_U(1,:) = 0;
         a_D(1,:) = 0;
         a_E(1,:) = 0;
         a_W(1,:) = 0;
%

%    specified phi = phi_d on lower boundary
%    -------------------------------------------
         a_P(end,:) = 1;
         b(end,:)   = phi_d;
         a_D(end,:) = 0;            
         a_U(end,:) = 0;
         a_E(end,:) = 0;
         a_W(end,:) = 0;
%


%%%%%%%%%%%%%%

%
%    Now call matrix solver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        phi = solver( a_E, a_W, a_U, a_D, a_P, b );
%
%


%%Peclet = ( rho(1,1) * w(1,1) * dz_u(1,1) ) / Gamma_P(1,1)
  Peclet = P_u(1,1);
%
%