function phi = fvm_soln_transient_bdry_nodes( X_edges, Z_edges, X_P, Z_P, ...
                                        nt, dt, rho_c_P, K_P, S_C, S_P, ...
                                        phi_0, phi_u, q_d, q_e, q_w )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  This m-file finds transient phi(x,z), given sources S and
%%  conductivity K_P.  S(x,z) and K_P(x,z) are constant in time.
%%
%%  Main program sets and delivers:
%%    - Mesh geometry (X_P, Z_P, X_edges, Z_edges)
%%    - time-marching info (nt, dt)
%%    - constitutive properties (K_P)
%%    - sources (S_C, S_P)
%%    - initial conditions (phi_0) 
%%    - boundary conditions (phi_u, q_d, q_e, q_w)
%%
%%  Uses Patankar Finite-Volume method.
%%
%%  THIS CODE USES PATANKAR HALF-VOLUMES
%%
%%  Boundary Conditions:
%%  --------------------
%%  At each time step -
%%  vector phi_u is specified value of phi on upper boundary
%%  vector q_w   is specified flux on west boundary
%%  vector q_e   is specified flux on east boundary
%%  vector q_d   is specified flux on down boundary
%%
%%
%%  Initial conditions:
%%  -------------------
%%      phi = phi_0
%%
%%  a_P = coefficient of phi_P at point P = (i,j)
%%  a_E = coefficient of phi_E at point east of (i,j)
%%  a_W, a_U, a_D for points West, Up, Down
%%  a_P_0 = coefficient of phi_P_0, known phi at previous time step
%%
%%  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,
%           initial conditions,and BC terms
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   set time-step partition (for demonstration purposes only!)
%   always use fully implicit scheme (theta=1)
%      theta = 1  fully implicit
%      theta =0.5 Crank-Nicolson
%      theta = 0  explicit
%   ------------------------------------------
    theta = 0.1;


%% INITIALIZING ...
%%%%%%%%%%%%%%%%%%%%

% Geometry
% --------
%
%  get dX, dZ - spacings between finite-volume edges
%  -------------------------------------------------
    dX = diff(X_edges')';
    dZ = diff(Z_edges);
%
%
%  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 [nz, nx]
    dX_e = [ dX_e dX_e(:,end) ];
%
    dX_w = diff(X_P')';
 % add dummy row at western edge to maintain size [nz, nx]
    dX_w = [ dX_w(:,1) dX_w ];
%
    dZ_u = diff(Z_P);
 % add dummy row at top edge to maintain size [nz, nx]
    dZ_u = [ dZ_u(1,:)' dZ_u' ]';
%
    dZ_d = diff(Z_P);
 % add dummy row at bottom edge to maintain size [nz, nx]
    dZ_d = [ dZ_d' dZ_d(end,:)' ]';
%
%
%  get f_e, f_w, f_u, f_d 
%          - fractional distances between volume interfaces
%            and adjacent volume centers
%  ---------------------------------------------------
    f_e = ones(size(X_P)) - (X_edges(:,2:end) - X_P)./ dX_e;
%
    f_w = ones(size(X_P)) - (X_P - X_edges(:,1:(end-1)) )./ dX_w;
%
    f_u = ones(size(Z_P)) - (Z_P - Z_edges(1:(end-1),: ))./ dZ_u;
%
    f_d = ones(size(Z_P)) - (Z_edges(2:end,: ) - Z_P)./ dZ_d;
%

%  Conductivities
%%%%%%%%%%%%%%%%%%%
%    If K = K(t), then conductivities need to be calculated 
%    inside time-stepping loop.
%
%  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 coefficients at adjacent points
% ---------------------------------------
%    If rho_c = rho_c(t), then a_P_0 needs to be calculated 
%    inside time-stepping loop.

    a_E   = (K_e ./ dX_e).* dZ;
    a_W   = (K_w ./ dX_w).* dZ;
    a_U   = (K_u ./ dZ_u).* dX;
    a_D   = (K_d ./ dZ_d).* dX;
    a_P_0 = (rho_c_P./dt).*dX.*dZ;

%
%
  %  plot current solution
        figure(5)
        surf( X_P, Z_P, phi_0 )
        colorbar
        set(gca,'zlim', [ min(min(phi_u)), 0 ] )
        xlabel('Distance  (m)')
        ylabel('Depth  (m)')
        zlabel('Temperature  (^oC)')

 pause


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Time-stepping Loop starts here 
%    (sources and B.C. may be time-dependent)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
     for i_time=1:nt
         
   % set up b matrix - includes constant part of source matrix S 
   %         S = S_C + S_P * phi_P    (units are W m^-3)
   % Boundary-Condition contribution is added below.
   % Initial-Condition contribution must be added at each time step
   %   i.e. b = S_C dX dZ + a_P_0.* T_P_0
   %   Calculated inside time step in case S_C = S_C(t)
%  ---------------------------------------------------------------
    %     b_0 = S_C.* dX.* dZ;
%
%  when using theta < 1, we need old solution at all neighbors
         phi_0_E = [ phi_0(:,2:end) phi_0(:, end) ];
         phi_0_W = [ phi_0(:,1) phi_0(:,1:end-1) ];
         phi_0_U = [ phi_0(1,:)' phi_0(1:end-1,:)' ]';
         phi_0_D = [ phi_0(2:end,:)' phi_0(end,:)' ]';

%         
  %  put initial condition (for this time step) into b
     %%    b = S_C.* dX.* dZ + a_P_0 .* phi_0;
         b = S_C.* dX.* dZ + ( 1 - theta ) * ...
            (a_E.*phi_0_E + a_W.*phi_0_W + a_U.*phi_0_U + a_D.*phi_0_D) ...
           + (a_P_0 - (1 - theta) * (a_E + a_W + a_U + a_D) ).* phi_0;
%

%
%
   % set up coefficient at central point P 
   %   (includes S_P for source dependence on phi(t) )
   %   Calculated inside time step in case S_P = S_P(t)
   % ----------------------------------------------------
      %%   a_P = a_E + a_W + a_U + a_D + a_P_0 - S_P.* dX.* dZ; 
         a_P = theta * (a_E + a_W + a_U + a_D) + a_P_0 - S_P.* dX.* dZ; 
%
%
   %  set Boundary Conditions
   %%%%%%%%%%%%%%%%%%%%%%%%%%%
     %  Since all coefficients have been initialized with
     %  matrix operations assuming all points are interior points,
     %  the coefficients for boundary points must now be updated to 
     %  correctly incorporate B.C.s.
%
   %  specified flux boundaries
   %  ---------------------------
   %   Boundary treatment recommended by Patankar (p. 50) 
   %   half-volumes at domain edges
   %        boundaries are at solution points X_B = X_P
   %   ------------------------------------------------
%
   %  Specified phi
   %  -------------
   %   Boundaries are at solution points X_P (Patankar)
   %   half-volumes at domain edges
   %    Set boundary-point coefficient to unity, b to specified value,
   %    zero out coefficients on neighboring points
   %
%
%

   %   East boundary - outward flux  (x is positive outward)
   %   =====================================================
     %  include specified boundary flux q_e in b
     %  remove a_E from a_P
    %%     a_P(:,end) = a_P(:,end) - a_E(:,end); 
         a_P(:,end) = a_P(:,end) - theta * a_E(:,end); 
         b(:,end) = b(:,end) - q_e(:,i_time).* dZ(:,end);
%
   %  eliminate fictitious East point (outside domain)
         a_E(:,end) = 0;
%

%   West boundary - zero inward flux (x is positive inward)
%   =======================================================
   %  include specified boundary flux q_w in b
   %  remove a_W from a_P
     %%    a_P(:,1) = a_P(:,1) - a_W(:,1); 
         a_P(:,1) = a_P(:,1) - theta * a_W(:,1); 
%
         b(:,1) = b(:,1) + q_w(:,i_time).* dZ(:,1);
%
   %  eliminate fictitious West point (outside domain)
         a_W(:,1) = 0;
%


%   Bottom - geothermal flux
%   ========================  
   %  include specified boundary flux q_d in b
   %  remove a_D from a_P
      %%   a_P(end,:) = a_P(end,:) - a_D(end,:);
         a_P(end,:) = a_P(end,:) - theta * a_D(end,:);
   %
         b(end,:) = b(end,:) - q_d(i_time,:).* dX(end,:);
%
   %  eliminate fictitious Down point (outside domain)
         a_D(end,:) = 0;
%

%
   %   Upper boundary - specified phi = phi_u(t) 
   %   =========================================
         a_P(1,:) = 1; 
         b(1,:) = phi_u(i_time,:);
         a_U(1,:) = 0; 
         a_E(1,:) = 0;
         a_W(1,:) = 0;
         a_D(1,:) = 0;
%
  %  Now call implicit solver
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %%   phi = solver( a_E, a_W, a_U, a_D, a_P, b );
        phi = solver( theta*a_E, theta*a_W, theta*a_U, theta*a_D, a_P, b );
%
%
  %  Save solution as phi_0 for next time step
        phi_0 = phi;
%
  %  plot current solution
        figure(5)
        if( size(X_P,2) <= 2)
            plot( Z_P(:,1), phi(:,1) )
            grid on;box on;
       %%     set(gca,'ylim', [ min(min(phi_u)), 0 ] )
            set(gca,'ylim', [ min( min(min(phi_u)),-max(max(abs(phi))) ), ...
                             max( 0, max(max(abs(phi))) ) ] )
            xlabel('Depth  (m)')
            ylabel('Temperature  (^oC)')

        else
            surf( X_P, Z_P, phi )
            colorbar
    %%        set(gca,'zlim', [ -30, -10 ] )
            set(gca,'zlim', [ min(min(phi_u)), 0 ] )
            xlabel('Distance  (m)')
            ylabel('Depth  (m)')
            zlabel('Temperature  (^oC)')
        end

 pause
        
        
     end
%