function phi = fvm_SS_adv_diff_2D( x_edges_vec, z_edges_vec, ...
                     x_P_vec, z_P_vec, 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.
%%    - Advective flux included through PATANKAR POWER-LAW SCHEME
%%    - Gamma_P is initial guess at diffusion coefficient
%%      Gamma is updated at each iteration using value of phi from
%%      previous iteration
%%
%%
%%  EXTERNALLY SPECIFIED VELOCITY MUST SATISFY CONTINUITY
%%   i.e. 
%%     div(rho u) = 0
%%
%%  Boundary Conditions:
%%
%%  Boundary Conditions:
%%   All boundaries are finite-volume interfaces
%%   (Patankar Practice B, page 69)
%%   However, there is an additional "point" X_B located on boundary.
%%   Its value and coefficients are determined by phi value or gradient
%%   boundary conditions there.
%%   Adjacent interior point is X_I
%%   (X_I - X_B) = dX_B
%%
%%    phi_B = phi_0
%%        a_B = 1,  all other a_J = 0
%%        b   = phi_0
%%
%%    d phi/dx = dphi_dx_0
%%        use gradient between X_B and X_I (interior point)
%%        d phi/dx|B = Gamma_I (phi_I - phi_B)/dx_B
%%    or
%%        phi_B = dphi_dx_0 * dx_B + phi_I
%%
%%  For each boundary there is a matrix bc_x
%%     with either 2 rows (Up,Down) or 2 columns (East,West).
%%  e.g. East  
%%   bc_e is [nz_P x 2] where nz_P is number of nodes on boundary 
%%     - first column contains value of BC at that node on boundary
%%     - second column specifies type of BC
%%            - 1 -> specified phi
%%            - 2 -> specified phi gradient across boundary
%%
%%   Similarly bc_w, bc_u, bc_d on other boundaries
%%
%%
%%  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 phi_ref Gamma_ref Gamma_grad  %  Gamma(phi) parameterization

    global  conv_test   max_iter
    global  which_Gamma

%
%
% get number of Finite-Volume nodes
% ----------------------------------
    nx_P = length(x_P_vec);
    nz_P = length(z_P_vec);
%
% get number of finite-volume faces
     nz_fv = nz_P - 2;
     nx_fv = nx_P - 2;
%
% set up arrays defining nodes locations
     X_P = repmat(x_P_vec, nz_P, 1);
     Z_P = repmat(z_P_vec, 1, nx_P);


%  get dX, dZ - spacings between volume edges
%  -------------------------------------------
   %%    dx_vec = [ 0  diff(x_edges_vec)  0];
   %%    dz_vec = [ 0  diff(z_edges_vec)'  0]';
   %% This leads to division by zero for Peclet numbers
      dx_vec = [ 1  diff(x_edges_vec)  1 ];
      dz_vec = [ 1  diff(z_edges_vec)'  1 ]';
%
      dX = repmat( dx_vec, nz_P, 1 );
      dZ = repmat( dz_vec, 1, nx_P );
%

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

%

%  get f_e, f_w, f_u, f_d 
%          - fractional distances to interface in adjacent volume
%  --------------------------------------------------------------
      f_e = [(1 -(x_edges_vec -x_P_vec(1:end-1))./dX_e(1,1:end-1)) 0];
      f_e = repmat( f_e, nz_P, 1 );
       
      f_w = [0 (1-(x_P_vec(2:end) -x_edges_vec )./dX_w(1,2:end) )];
      f_w = repmat( f_w, nz_P, 1 );
       
      f_u = [0 (1 -(z_P_vec(2:end) -z_edges_vec)./dZ_u(2:end,1) )']';
      f_u = repmat( f_u, 1, nx_P );
%
      f_d = [(1 -(z_edges_vec -z_P_vec(1:end-1) )./dZ_d(1:end-1,1))' 0]';
      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;
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% iteration loop on phi-dependent coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      phi      = repmat(phi_ref, size(X_P) ); % initial guess
      dphi_max = 2*conv_test;        % to get into convergence loop
%
%  initialize matrix to save intermediate solutions
      phi_save = [];
      iter = 0;      % iteration number
%
     while ( (dphi_max > conv_test) & (iter < max_iter) )
         
        phi_old = phi;
        iter = iter+1;
        
%
%   set coefficients based on old solution phi_old
       if( which_Gamma == 0)
          Gamma_P = Gamma_phi(phi_old);
       elseif(which_Gamma == 1 )
          Gamma_P = Gamma_phi_bad(phi_old);
       end
%
%
% 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 positions X and Z where velocities u and w
%  are calculated 
%  - u on East and West interfaces
%  - w on Up and Down interfaces
      X_u_velo = repmat( x_edges_vec, nz_P, 1 );
      Z_u_velo = repmat( z_P_vec, 1, nx_P-1 );
      X_w_velo = repmat( x_P_vec, nz_P-1, 1 );
      Z_w_velo = repmat( z_edges_vec, 1, nx_P );
%
%  find the velocities
      u_edges = u(X_u_velo, Z_u_velo);
      w_edges = w(X_w_velo, Z_w_velo);
%
%  split them out into East and West, Up and Down
      u_e = [ u_edges  u_edges(:, end) ];
      u_w = [ u_edges(:, 1)  u_edges ];
%
      w_u = [ w_edges(1,:)'  w_edges' ]';
      w_d = [ w_edges'  w_edges(end,:)' ]';
%
%  clean up work space
      clear u_edges w_edges
%
%
% 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
% ---------------------------------------
      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 );
 %

%         
% 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 values phi_P = phi_value at x_P
%     set a_P = 1, all other a_X = 0
%           b = phi_value = bc_x(index, 1);
%
%  specified gradients
%     gradient = (phi_I - phi_P)/dX_I  where dX_I = (X_I - X_P)
%  or
%     phi_P = phi_I + dX_I * gradient 
%  so
%    a_P = 1
%    a_I = 1
%    b   = dX_I * gradient
%  -------------------------------------
%
%  East Boundary
%  --------------
      a_P(:,end) = 1; 
      a_U(:,end) = 0;
      a_D(:,end) = 0;
      a_W(:,end) = 0;
      a_E(:,end) = 0;
%     
%  recover info about type of BC on each volume
      bc_type = bc_e(:,2);
%
  % index identifies boundary nodes where phi is specified
      index = find(bc_type == 1);
      if( ~isempty(index) )
        b(index, end) = bc_e(index, 1);
      end
%
  % index identifies boundary nodes where gradient is specified
      index = find( bc_type == 2 );
      if( ~isempty(index) )
         a_W(index, end) = 1;
         b(index, end)   = dX_e(index, end).* bc_e(index, 1);
      end
%
  % index identifies boundary nodes where diffusive flux is specified
      index = find(bc_type == -2);
      if( ~isempty(index) )
         a_W(index, end) = 1;
         b(end, index)   = - dZ_w(index, end) ...
                            .* bc_e(index,1)./Gamma_P(index, end);
      end

%
%  West Boundary
%  --------------
      a_P(:,1) = 1; 
      a_U(:,1) = 0;
      a_D(:,1) = 0;
      a_W(:,1) = 0;
      a_E(:,1) = 0;
%
      bc_type = bc_w(:,2);
%
  % index identifies boundary nodes where phi is specified
      index = find( bc_type == 1 );
      if( ~isempty(index) )
         b(index, 1) = bc_w(index, 1);
      end
%
  % index identifies boundary nodes where gradient is specified
      index = find( bc_type == 2 );
      if( ~isempty(index) )
         a_E(index, 1) = 1;
         b(index, 1)   = - dX_w(index,1).* bc_w(index, 1);
      end
%
  % index identifies boundary nodes where diffusive flux is specified
      index = find(bc_type == -2);
      if( ~isempty(index) )
         a_E(1, index) = 1;
         b(index, 1)   = - dX_w(index, 1) ...
                            .* bc_w(index, 1)./Gamma_P(1, index);
      end

%
%
%  Up Boundary
%  --------------
      a_P(1,:) = 1; 
      a_U(1,:) = 0;
      a_D(1,:) = 0;
      a_W(1,:) = 0;
      a_E(1,:) = 0;
%
      bc_type = bc_u(2,:);
%
  % index identifies boundary nodes where phi is specified
      index = find( bc_type == 1 );
      if( ~isempty(index) )
         b(1, index) = bc_u(1, index);
      end
%
  % index identifies boundary nodes where gradient is specified
      index = find( bc_type == 2 );
      if( ~isempty(index) )
         a_D(1, index) = 1;
         b(1, index)   = - dZ_d(1, index).* bc_u(1, index);
      end
%
  % index identifies boundary nodes where diffusive flux is specified
      index = find(bc_type == -2);
      if( ~isempty(index) )
         a_D(1, index) = 1;
         b(end, index)   = - dZ_d(1, index) ...
                            .* bc_u(1, index)./Gamma_P(1, index);
      end

      %
%
%  Down Boundary
%  --------------
      a_P(end,:) = 1; 
      a_U(end,:) = 0;
      a_D(end,:) = 0;
      a_W(end,:) = 0;
      a_E(end,:) = 0;
%
      bc_type = bc_d(2,:);
%
  % index identifies boundary nodes where phi is specified
      index = find( bc_type == 1 );
      if( ~isempty(index) )
         b(end, index) = bc_d(1, index);
      end
%
  % index identifies boundary nodes where gradient is specified
      index = find(bc_type == 2);
      if( ~isempty(index) )
         a_U(end, index) = 1;
         b(end, index)   = dZ_u(end, index).* bc_d(1, index);
      end
%
  % index identifies boundary nodes where diffusive flux is specified
      index = find(bc_type == -2);
      if( ~isempty(index) )
         a_U(end, index) = 1;
    %%     b(end, index)   = dZ_u(end, index).* bc_d(1, index);
         b(end, index)   = - dZ_u(end, index) ...
                            .* bc_d(1, index)./Gamma_P(end,index);
      end
%%

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Now call matrix solver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      phi = solver( a_E, a_W, a_U, a_D, a_P, b );
%
      dphi_max = max(max(abs( phi - phi_old) ) );
%

%  plot current solution
%  -----------------------
      if( which_Gamma == 0 )
          figure(2)
      elseif( which_Gamma == 1 )
          figure(3)
      end

      plot(Z_P(:,1), phi(:,1) )
      hold on;grid on; box on
      xlabel('Depth  (m)')
      ylabel('phi')
      
      iter
      dphi_max
      min_Gamma = min(min(Gamma_P))

pause


   end   %        while (dphi_max > conv_test)



%
%