function free_fem_navier_stokes (  node_file_name, element_file_name )

%% MAIN is the main routine of FREE_FEM_NAVIER_STOKES.
%
%  Discussion:
%
%    This program solves the steady incompressible Navier Stokes equations for 
%    velocity vector W and scalar pressure P:
%
%      -nu * Laplacian W(X,Y) + W dot Grad W + Grad P(X,Y) = F(X,Y)
%
%                                               Div W(X,Y) = G(X,Y)
%
%    in an arbitrary triangulated region in the plane.  
%
%    Let U and V denote the scalar components of the velocity vector W.
%
%    Along the boundary of the region, the user controls the type of
%    boundary condition to be imposed, if any.  Currently, these
%    conditions may be of Dirichlet form:
%
%      U(X,Y) = U_BC(X,Y)
%      V(X,Y) = V_BC(X,Y)
%      P(X,Y) = P_BC(X,Y)
%
%    or Neumann form with ZERO right hand side:
%
%      dU/dn(X,Y) = 0
%      dV/dn(X,Y) = 0
%      dP/dn(X,Y) = 0
%
%    The code uses the finite element method.  The Taylor-Hood element
%    is used, in which a single reference triangle is used to define
%    both a piecewise quadratic representation of velocity, and a piecewise
%    linear representation of pressure.
%
%  Geometry specification:
%
%    The user defines the geometry by supplying two data files
%    which list the node coordinates, and list the nodes that make up
%    each element.
%
%  Equation specification:
%
%    The user specifies 
%
%    * the kinematic viscosity NU;
%
%       function nu = constants ( 'DUMMY' )
%
%    * the right hand side of the Stokes equations:
%
%      function [ u_rhs, v_rhs, p_rhs ] = rhs ( node_num, node_xy )
%
%    * the type of boundary conditions imposed:
%
%      [ node_u_condition, node_v_condition, node_p_condition ] =
%        boundary_type ( node_num, node_xy, node_boundary, node_type,
%        node_u_condition, node_v_condition, node_p_condition )
%
%    * the right hand side of any Dirichlet boundary conditions:
%
%      function [ u_bc, v_bc, p_bc ] = dirichlet_condition ( node_num, node_xy )
%
%    * the right hand side of any Neumann boundary conditions:
%      (currently, nonzero values will be ignored)
%
%      function [ u_bc, v_bc, p_bc ] = neumann_condition ( node_num, node_xy )
%
%  Usage:
%
%    free_fem_navier_stokes ( node_file, element_file )
%
%    invokes the program:
%
%    * "node_file" contains the coordinates of the nodes;
%    * "element_file" contains the indices of nodes that make up each element.
%
%    Graphics files created include:
%
%    * "nodes6.eps", an image of the nodes;
%    * "triangles6.eps", an image of the elements;
%
%    Data files created include:
%
%    * "nodes3.txt", the nodes associated with pressure;
%    * "triangles3.txt", the linear triangles associated with pressure;
%    * "stokes_pressure3.txt", the Stokes pressure at the pressure nodes;
%    * "stokes_velocity6.txt", the Stokes velocity at the velocity nodes.
%    * "pressure3.txt", the pressure at the pressure nodes;
%    * "velocity6.txt", the velocity at the velocity nodes;
%
%  Modified:
%
%    26 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string NODE_FILE_NAME, the name of the node file.  If
%    this argument is not supplied, it will be requested.
%
%    Input, string ELEMENT_FILE_NAME, the name of the element file.  If
%    this argument is not supplied, it will be requested.
%
%  Local parameters:
%
%    Local, real A(3*IB+1,VARIABLE_NUM), the VARIABLE_NUM by VARIABLE_NUM
%    coefficient matrix (for the Stokes problem) or jacobian matrix
%    (for the Navier Stokes problem).
%
%    Local, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); 
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Local, integer ELEMENT_NUM, the number of elements.
%
%    Local, integer ELEMENT_ORDER, the element order.
%
%    Local, real F(VARIABLE_NUM), the right hand side.
%
%    Local, integer IB, the half-bandwidth of the matrix.
%
%    Local, integer IT_MAX, the maximum number of Newton iterations allowed.
%    0 is a legal value.  It simply solves the Stokes problem, computes
%    the Navier-Stokes residual, and stops.
%    1 does the above, and then takes one Newton step, and so on.
%
%    Local, real NODE_C(VARIABLE_NUM), the finite element
%    coefficients of the current solution estimate.
%
%    Local, real NODE_C_DEL(VARIABLE_NUM), the correction to
%    the finite element coefficients of the current solution estimate.
%
%    Local, real NODE_C_OLD(VARIABLE_NUM), the finite element
%    coefficients computed by the prefious step of the iteration.
%
%    Local, integer NODE_NUM, the number of nodes.
%
%    Local, integer NODE_P_CONDITION(NODE_NUM), 
%    indicates the condition used to determine pressure at a node.
%    0, there is no condition at this node.
%    1, a finite element equation is used;
%    2, a Dirichlet condition is used. 
%    3, a Neumann condition is used. 
%
%    Local, integer NODE_P_VARIABLE(NODE_NUM), 
%    is the index of the pressure variable associated with the node,
%    or -1 if there is no associated pressure variable.
%
%    Local, integer NODE_TYPE(NODE_NUM), determines if the node is a 
%    vertex or midside node.
%    1, the node is a vertex (P, U, V variables are associated with it).
%    2, the node is a midside node (only U and V variables are associated.)
%
%    Local, integer NODE_U_CONDITION(NODE_NUM), 
%    indicates the condition used to determine horizontal velocity at a node.
%    0, there is no condition at this node.
%    1, a finite element equation is used;
%    2, a Dirichlet condition is used. 
%    3, a Neumann condition is used. 
%
%    Local, integer NODE_U_VARIABLE(NODE_NUM), 
%    is the index of the horizontal velocity variable associated with the node.
%
%    Local, integer NODE_V_CONDITION(NODE_NUM), 
%    indicates the condition used to determine vertical velocity at a node.
%    0, there is no condition at this node.
%    1, a finite element equation is used;
%    2, a Dirichlet condition is used.
%    3, a Neumann condition is used.  
%
%    Local, integer NODE_V_VARIABLE(NODE_NUM), 
%    is the index of the vertical velocity variable associated with the node.
%
%    Local, real NODE_XY(2,NODE_NUM), the coordinates of nodes.
%
%    Local, integer NODE3_NUM, the number of pressure nodes.
%
%    Local, integer NODE3_LABEL(NODE_NUM), contains the renumbered
%    labels of pressures nodes, and -1 for nodes that are not pressure nodes.
%
%    Local, real NU, the kinematic viscosity.
%
%    Local, integer QUAD_NUM, the number of quadrature points used for 
%    assembly.  This is currently set to 3, the lowest reasonable value.
%    Legal values are 1, 3, 4, 6, 7, 9, 13, and for some problems, a value 
%    of QUAD_NUM greater than 3 may be appropriate.
%
%    Local, integer VARIABLE_NUM, the number of variables.
%
  debugging = 0;
  it_max = 5;

  nu = constants ( 'DUMMY' );

  quad_num = 7;
  res_tol = 0.0001;

  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FREE_FEM_NAVIER_STOKES\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Finite element solution of the \n' );
  fprintf ( 1, '  steady incompressible Navier Stokes equations\n' );
  fprintf ( 1, '  on a triangulated region in 2 dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  - nu * ( Uxx + Uyy ) + UUx + VUy + dPdx = F1(x,y)\n' );
  fprintf ( 1, '  - nu * ( Vxx + Vyy ) + UVx + VVy + dPdy = F2(x,y)\n' );
  fprintf ( 1, '                          Ux +  Vy        = F3(x,y).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Boundary conditions may be of Dirichlet type:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    U(x,y) = U_BC(x,y)\n' );
  fprintf ( 1, '    V(x,y) = V_BC(x,y)\n' );
  fprintf ( 1, '    P(x,y) = P_BC(x,y)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  or of Neumann type with zero right hand side:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    dU/dn(x,y) = 0\n' );
  fprintf ( 1, '    dV/dn(x,y) = 0\n' );
  fprintf ( 1, '    dP/dn(x,y) = 0\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The finite element method uses Taylor-Hood\n' );
  fprintf ( 1, '  triangular elements which are linear for pressure\n' );
  fprintf ( 1, '  and quadratic for velocity.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum number of Newton iterations IT_MAX = %d:\n', it_max );
  fprintf ( 1, '  Quadrature order =       %d\n', quad_num );
  fprintf ( 1, '  Kinematic viscosity NU = %f\n', nu );
%
%  Get the file names, if not specified on command line.
%
  if ( nargin < 2 )
    element_file_name = [];
  end
  
  if ( nargin < 1 )
    node_file_name = [];
  end
  
  [ node_file_name, element_file_name ] = file_name_specification ( ...
    node_file_name, element_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Node file is "%s".\n', node_file_name );
  fprintf ( 1, '  Element file is "%s".\n', element_file_name );
%
%  Read the node coordinate file.
%
  [ dim_num, node_num ] = dtable_header_read ( node_file_name );

  fprintf ( 1, '  Number of nodes =          %d\n', node_num );

  node_xy = dtable_data_read ( node_file_name, dim_num, node_num );

  r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, 2, 10, ...
    '  First 10 nodes' );
%
%  Read the element file.
%
  [ element_order, element_num ] = itable_header_read ( element_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Element order =            %d\n', element_order );
  fprintf ( 1, '  Number of elements =       %d\n', element_num );

  if ( element_order ~= 6 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FREE_FEM_NAVIER_STOKES - Fatal error!\n' );
    fprintf ( 1, '  The input triangulation has order %d.\n', element_order );
    fprintf ( 1, '  However, a triangulation of order 6 is required.' );
    error ( 'FREE_FEM_NAVIER_STOKES - Fatal error!' );
  end

  element_node = itable_data_read ( element_file_name, element_order, ...
    element_num );

  i4mat_transpose_print_some ( 6, element_num, ...
    element_node, 1, 1, 6, 10, '  First 10 elements' );
%
%  Determine which nodes are boundary nodes.
%
  node_boundary = triangulation_order6_boundary_node ( node_num, ...
    element_num, element_node );

  if ( 0 == 1 )
    lvec_print ( node_num, node_boundary, '    Node  Boundary?' );
  end
%
%  Determine the node type (vertex or midside):
%
  node_type(1:node_num) = 1;

  for element = 1 : element_num
    for j = 4 : 6
      node = element_node(j,element);
      node_type(node) = 2;
    end
  end
%
%  Determine the node conditions:
%
%  All conditions begin as finite element conditions.
%
  node_p_condition(1:node_num) = 1;
  node_u_condition(1:node_num) = 1;
  node_v_condition(1:node_num) = 1;
%
%  Conditions on velocities associated with a boundary node are Dirichlet
%  conditions.
%
  for node = 1 : node_num
    if ( node_boundary(node) )
      node_u_condition(node) = 2;
      node_v_condition(node) = 2;
    end
  end
%
%  Midside nodes have no associated pressure variable.
%
  for node = 1 : node_num
    if ( node_type(node) == 2 )
      node_p_condition(node) = 0;
    end
  end
%
%  Replace a single finite element pressure condition by a Dirichlet
%  condition.
%
  p_node = -1;
  for node = 1 : node_num
    if ( node_p_condition(node) == 1 )
      node_p_condition(node) = 2;
      p_node = node;
      break
    end
  end

  if ( p_node == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FREE_FEM_NAVIER_STOKES - Fatal error!\n' );
    fprintf ( 1, '  Unable to find a finite element pressure condition\n' );
    fprintf ( 1, '  suitable for replacement by a Dirichlet condition.\n' );
    error ( 'FREE_FEM_NAVIER_STOKES - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Dirichlet boundary condition on pressure\n' );
  fprintf ( 1, '  will be applied at node %d\n', p_node );
%
%  Number the variables.
%
  variable_num = 0;

  for node = 1 : node_num

    variable_num = variable_num + 1;
    node_u_variable(node) = variable_num;

    variable_num = variable_num + 1;
    node_v_variable(node) = variable_num;

    if ( node_type(node) == 1 )
      variable_num = variable_num + 1;
      node_p_variable(node) = variable_num;
    else
      node_p_variable(node) = -1;
    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Total number of variables is %d\n', variable_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Variable indices per node:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Node  U_index  V_index  P_index\n' );
  fprintf ( 1, '\n' );
  for node = 1 : node_num
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', ...
      node, node_u_variable(node), node_v_variable(node), ...
      node_p_variable(node) );
  end
%
%  Determine the bandwidth of the coefficient matrix.
%
  ib = bandwidth ( element_num, element_node, ...
    node_num, node_p_variable, node_u_variable, node_v_variable );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix half bandwidth is %d\n', ib );
  fprintf ( 1, '  The matrix bandwidth is      %d\n', 2 * ib + 1 );
  fprintf ( 1, '  The storage bandwidth is     %d\n', 3 * ib + 1 );
%
%  Plot the nodes.
%
  if ( node_num <= 100 )

    file_name = 'nodes6.eps';
    node_label = 1;

    points_plot ( file_name, node_num, node_xy, node_label );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Order 6 nodes plotted in "%s".\n', file_name );

  end
%
%  Plot the elements.
%
  if ( node_num <= 100 )

    file_name = 'triangles6.eps';
    node_show = 2;
    triangle_show = 2;

    triangulation_order6_plot ( file_name, node_num, ...
      node_xy, element_num, element_node, node_show, triangle_show );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Order 6 triangles plotted in "%s".\n', ...
      file_name );

  end
%
%  Get an initial condition, by assembling the Stokes coefficient matrix A 
%  and the right-hand side F, and solving.
%
  [ a, f ] = assemble_stokes ( node_num, element_num, quad_num, ...
    variable_num, node_xy, node_p_variable, node_u_variable, ...
    node_v_variable, element_node, nu, ib );

  if ( debugging )

    dgb_print_some ( variable_num, variable_num, ib, ib, a, 1, 1, 10, 10, ...
      '  Initial block of Stokes matrix A:' );

    r8vec_print_some ( variable_num, f, 1, 10, ...
      '  Part of Stokes right hand side vector:' );

  end
%
%  Adjust the linear system to account for Dirichlet boundary conditions.
%
  [ a, f ] = dirichlet_apply ( node_num, node_xy, node_p_variable, ...
    node_u_variable, node_v_variable, node_p_condition, ...
    node_u_condition, node_v_condition, variable_num, ib, a, f );

  if ( debugging )

    dgb_print_some ( variable_num, variable_num, ib, ib, a, 1, 1, 10, 10, ...
      '  Stokes matrix after boundary condition adjustments:' );

    r8vec_print_some ( variable_num, f, 1, 10, ...
      '  Part of adjusted Stokes right hand:' );

  end
%
%  Solve the linear system using a banded solver.
%
  [ a_lu, pivot, info ] = dgb_fa ( variable_num, ib, ib, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FREE_FEM_NAVIER_STOKES - Fatal error!\n' );
    fprintf ( 1, '  DGB_FA returned an error condition.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The linear system was not factored, and the\n' );
    fprintf ( 1, '  algorithm cannot proceed.\n' );
    error ( 'FREE_FEM_NAVIER_STOKES - Fatal error!' );
  end

  job = 0;
  node_c = dgb_sl ( variable_num, ib, ib, a_lu, pivot, f, job );

  if ( debugging )
    r8vec_print_some ( variable_num, node_c, 1, 10, ...
      '  Part of the solution vector:' );
  end
%
%  Print the Stokes solution vector based at nodes.
%
  if ( debugging )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Solution to the STOKES equations:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '    Node        U        V        P\n' );
    fprintf ( 1, '\n' );

    for node = 1 : node_num

      iu = node_u_variable(node);
      iv = node_v_variable(node);
      ip = node_p_variable(node);

      if ( 0 < ip )
        fprintf ( 1, '  %6d  %14f  %14f  %14f\n', ...
          node, node_c(iu), node_c(iv), node_c(ip) );
      else
        fprintf ( 1, '  %6d  %14f  %14f\n', ...
          node, node_c(iu), node_c(iv) );
      end

    end

  end
%
%  Compute a renumbering of the pressure nodes.
%
  node3_num = 0;

  for node = 1 : node_num
    if ( node_type(node) == 1 )
      node3_num = node3_num + 1;
    end
  end

  node3_num = 0;

  for node = 1 : node_num
    if ( node_type(node) == 1 )
      node3_num = node3_num + 1;
      node3_label(node) = node3_num;
    else
      node3_label(node) = -1;
    end
  end
%
%  Write the pressure nodes to a file.
%
  file_name = 'nodes3.txt';

  nodes3_write ( file_name, node_num, node_xy, node_type );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Pressure nodes written to "%s".\n', file_name );
%
%  Write the pressure triangles to a file.
%
  file_name = 'triangles3.txt';

  triangles3_write ( file_name, element_num, element_node, ...
    node_num, node3_label );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Pressure triangles written to "%s".\n', file_name );
%
%  Write the pressures to a file.
%
  if ( 0 )

    file_name = 'stokes_pressure3.txt';

    pressure3_write ( file_name, node_num, node_p_variable, ...
      variable_num, node_c );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Stokes pressures written to "%s".\n', file_name );

  end
%
%  Write an ASCII file that can be read into MATLAB.
%
  if ( 0 )

    file_name = 'stokes_velocity6.txt';

    velocity6_write ( file_name, node_num, node_u_variable, ...
      node_v_variable, variable_num, node_c );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Stokes velocities written to "%s".\n', file_name );

  end
%
%  The Stokes solution provides an initial guess for the Newton iteration
%  applied to the Navier Stokes equations.
%
  it_num = 0;

  while ( 1 )

    node_r = residual_fem ( node_num, node_xy, element_num, ...
      element_node, quad_num, node_u_variable, ...
      node_v_variable, node_p_variable, variable_num, nu, node_c );

    if ( debugging )
      r8vec_print_some ( variable_num, node_r, 1, 10, ...
        '  Part of Navier-Stokes FEM residual:' );
    end

    res_norm = sqrt ( sum ( node_r(1:node_num).^2 ) );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  l2-norm of FEM residual = %f\n', res_norm );

    node_r = residual_adjust_dirichlet ( node_num, node_xy, node_p_variable, ...
      node_u_variable, node_v_variable, node_p_condition, ...
      node_u_condition, node_v_condition, variable_num, node_c, node_r );

    if ( debugging )
      r8vec_print_some ( variable_num, node_r, 1, 10, ...
        '  Part of Navier-Stokes FEM residual adjusted for BC:' );
    end

    node_r_norm = sqrt ( sum ( node_r(1:node_num).^2 ) );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  l2-norm of adjusted FEM residual = %f\n', node_r_norm );

    if ( node_r_norm <= res_tol )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Convergence.\n' );
      break;
    end

    if ( it_max <= it_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  The iteration limit has been reached.\n' );
      break
    end

    it_num = it_num + 1;
 
    a = jacobian_fem ( node_num, node_xy, element_num, ...
      element_node, quad_num, node_u_variable, ...
      node_v_variable, node_p_variable, variable_num, nu, node_c, ib ); 
%
%  Adjust the jacobian for boundary conditions.
%
    a = jacobian_adjust_dirichlet ( node_num, node_xy, ...
      node_p_variable, node_u_variable, node_v_variable, ...
      node_p_condition, node_u_condition, node_v_condition, ...
      variable_num, ib, a );
%
%  Factor the jacobian.
%
    [ a, pivot, ierr ] = dgb_fa ( variable_num, ib, ib, a );

    if ( ierr ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'FREE_FEM_NAVIER_STOKES - Fatal error!\n' );
      fprintf ( 1, '  The Jacobian matrix is singular.\n' );
      error ( 'FREE_FEM_NAVIER_STOKES - Fatal error!' );
    end
%
%  Set up and solve the Newton system J * dX = - F(x)
%
    node_c_del(1:variable_num) = -node_r(1:variable_num);

    node_c_del = dgb_sl ( variable_num, ib, ib, a, pivot, node_c_del, job );

    node_c_old(1:variable_num) = node_c(1:variable_num);
    node_c(1:variable_num) = node_c(1:variable_num) + node_c_del(1:variable_num);

    if ( debugging )
      r8vec_print_some ( variable_num, node_c_del, 1, 10, ...
        '  Part of Newton correction vector:' );
    end

    node_c_del_norm = sqrt ( sum ( node_c_del(1:variable_num).^2 ) );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  l2-norm of Newton correction = %f\n', node_c_del_norm );

  end
%
%  Print the Navier Stokes solution vector based at nodes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution to the NAVIER STOKES equations:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Node        U        V        P\n' );
  fprintf ( 1, '\n' );

  for node = 1 : node_num

    iu = node_u_variable(node);
    iv = node_v_variable(node);
    ip = node_p_variable(node);

    if ( 0 < ip )
      fprintf ( 1, '  %6d  %14f  %14f  %14f\n', ...
       node, node_c(iu), node_c(iv), node_c(ip) );
    else
      fprintf ( 1, '  %6d  %14f  %14f\n', ...
       node, node_c(iu), node_c(iv) );
    end

  end
%
%  Write the pressures to a file.
%
  file_name = 'pressure3.txt';
  
  pressure3_write ( file_name, node_num, node_p_variable, ...
    variable_num, node_c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Navier Stokes pressures written to "%s".\n', file_name );
%
%  Write the velocities to a file.
%
  file_name = 'velocity6.txt';

  velocity6_write ( file_name, node_num, node_u_variable, ...
    node_v_variable, variable_num, node_c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Navier Stokes velocities written to "%s".\n', file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FREE_FEM_NAVIER_STOKES:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp;