function phi = fvm_soln_SS_adv_diff( X_edges, Z_edges, X_P, Z_P, K_P, ...
                 rho, c, u_velo, w_velo, S_C, S_P, phi_u, phi_d, q_e, q_w )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  find steady-state phi(x,z), for combined advection and diffusion
%%  with given heat sources S and conductivity K_P.
%%
%%      div( rho u phi) = div( K grad(phi) ) + S
%%
%%  Uses Patankar Finite Volume method.
%%
%%
%%  THIS CODE INCLUDES ADVECTIVE FLUX WITH CENTERED DIFFERENCE
%%            (NO UPWINDING SCHEME)  
%%    *** DON'T USE THIS AT HOME FOR SERIOUS WORK!!  ***
%%
%%             EXTERNALLY SPECIFIED VELOCITY MUST SATISFY CONTINUITY,
%%            i.e. IS A CONSTANT
%%
%%  Boundary Conditions:
%%  vector phi_u is fixed value of phi on upper boundary
%%  vector phi_d is fixed value of phi on lower boundary
%%  vector q_w   is fixed flux on west boundary
%%  vector q_e   is fixed flux on east boundary
%%
%% In this treatment, flux boundary conditions are applied on edges 
%% of finite volumes. This means that I don't have to use half-volumes.
%% Function-value boundary conditions are applied on nodes, i.e. with
%% half-volumes at Up and Down boundaries.
%%
%%  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
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% 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 where boundaries are at interfaces
       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
% ------------------------------------------------------
    K_E   =  [ K_P(:,2:end) K_P(:, end) ];
    K_W   =  [ K_P(:,1) K_P(:,1:end-1) ];
    K_U   =  [ K_P(1,:)' K_P(1:end-1,:)' ]';
    K_D   =  [ K_P(2:end,:)' K_P(end,:)' ]';
%
%
% set up interface conductivity matrices
% --------------------------------------
    K_e =  1./ ( (1 - f_e)./K_P + f_e./K_E );
    K_w =  1./ ( (1 - f_w)./K_P + f_w./K_W );
    K_u =  1./ ( (1 - f_u)./K_P + f_u./K_U );
    K_d =  1./ ( (1 - f_d)./K_P + f_d./K_D );
%
%
% set up velocities at interfaces
%% (this is quick and dirty because velocity is uniform)
% ------------------------------------------------------
    u_e = u_velo;
    u_w = u_velo;
    w_u = w_velo;
    w_d = w_velo;

%
% set up coefficients at adjacent points
% ---------------------------------------
    a_E = ( K_e ./ dX_e - rho.* c.* u_e/2 ).* dZ;
    a_W = ( K_w ./ dX_w + rho.* c.* u_w/2 ).* dZ; 
    a_U = ( K_u ./ dZ_u + rho.* c.* w_u/2 ).* dX;
    a_D = ( K_d ./ dZ_d - rho.* c.* w_d/2 ).* dX;

         
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%    specified flux boundaries
%    ---------------------------
%     East boundary - outward flux  (q_w = 0)
%           use  q_e = - K_e * (T_E - T_P)/dx   (flux toward E) 
%          i.e.  T_E = T_P - dx * (q_e/K_e) to eliminate T_E
         a_P(:,end) = a_P(:,end) - a_E(:,end); 
         b(:,end) = b(:,end) - dX_e.* (q_e./K_e(:,end)) .* a_E(:,end);
         a_E(:,end) = 0;

%     West boundary - zero inward flux (q_w = 0)
%           use  q_w = - K_w * (T_P - T_W)/dx  (flux from W)
%          i.e.  T_W = T_P + dx * (q_w/K_w) to eliminate T_W
         a_P(:,1)   = a_P(:,1) - a_W(:,1);
         b(:,1) = b(:,1) + dX_w.* (q_w./K_w(:,1)) .* a_W(:,1);
         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 );
%
%
%
%