Interactive session for prepared example 10



        *******************************************************            
        ****  Welcome to the PDE2D 9.6 Interactive Driver  ****            
        *******************************************************            
                                                                           
   PDE2D can solve systems of (linear or nonlinear) steady-state,          
   time-dependent and eigenvalue partial differential equations            
   in 1D intervals, general 2D regions, and in a wide range of simple      
   3D regions.  Ordinary differential equation systems can also be         
   solved.                                                                 
                                                                           
   You will now be asked a series of interactive questions about your      
   problem.  The answers you give will be used to construct a PDE2D        
   FORTRAN program, which can then be compiled and linked with the PDE2D   
   runtime routines to produce an executable program.  The FORTRAN         
   driver program created will be well-documented and highly readable      
   (most of the interactive messages are repeated in the comments),        
   so that minor modifications or corrections can be made directly to      
   the FORTRAN program, without the need to work through a new             
   interactive session.                                                    
                                                                           
   You can alternatively create your PDE2D FORTRAN program using the       
   PDE2D graphical user interface (GUI) ("pde2d_gui [progname]").  It      
                                                                      [RETURN]
   is extraordinarily easy to set up problems using the PDE2D GUI,         
   which handles 0D and 1D problems, and 2D and 3D problems in "a wide     
   range of simple regions".  However, the PDE2D GUI cannot handle         
   complex regions, so if you have a complex 2D region you must use        
   this Interactive Driver.                                                
                                                                           
   If this is your first time to use PDE2D, you may want to work           
   through an example problem before trying one of your own.  Do you       
   want to work through a prepared example?                                
|---- Enter yes or no
yes
   Several prepared examples are available.  Enter:                        
                                                                           
       1 - to see a simple problem: a simply-supported elastic plate       
           equation, with a unit load concentrated at the midpoint of      
           a square.                                                       
       2 - to see a more complex problem: a non-linear, steady-state       
           PDE, solved in an annulus.  Dirichlet (U = ...) boundary        
           conditions are imposed on part of the boundary, and Neumann     
           (dU/dn = ...) conditions are imposed on the other part, in      
           this example.  The initial triangulation is generated           
           automatically, and adaptive grid refinement is illustrated.     
       3 - to see an eigenvalue problem.  The region has a curved          
           interface across which material properties vary abruptly, in    
           this example.                                                   
       4 - to see the first part of a thermal stress problem.  In this     
           part, the temperature distribution in a V-notched block is      
           calculated by solving the time-dependent heat conduction        
           equation, using adaptive time step control.                     
       5 - to see the second part of a thermal stress problem.  In         
           this part, the stresses induced in the V-notched block by       
           thermal expansion are calculated, using the temperature         
           distribution output by example 4.  You must run example 4       
                                                                      [RETURN]
           and save the tabular output before you can run example 5.       
           Examples 4 and 5 illustrate communication between problems.     
       6 - to see a 1D time-dependent integro-differential equation for    
           a financial math application.  In this problem there is a       
           term involving an integral of the solution, which requires      
           that we use PDE2D's feature for interpolating the solution      
           at the last saved time step, for use in the integral term.      
       7 - to see a waveguide problem (an eigenvalue problem in which      
           the eigenvalue appears nonlinearly).  This example shows how    
           to handle boundary conditions of different types on the same    
           arc, and how to produce a plot of a computed integral vs time.  
       8 - to see the Navier-Stokes equations (penalty formulation)        
           solved for a fluid flowing around a bend.                       
       9 - to see a 3D elasticity problem, solved in a torus.  This        
           example illustrates the use of user-defined coordinate          
           transformations to handle more general 3D regions.              
      10 - to see a time-dependent wave equation (reduced to a system      
           of two PDEs), solved in a 3D box.                               
      11 - to see a 3D eigenvalue problem (the Schrodinger equation        
           in a hydrogen atom).  This example illustrates the use of       
           spherical coordinates and periodic boundary conditions.         
      12 - to see a 3D eigenvalue problem, solved in a composite region    
                                                                      [RETURN]
           consisting of two cylinders of different material properties.   
      13 - to see the axisymmetric Navier-Stokes equations solved in a     
           non-rectangular channel, using the collocation FEM.             
      14 - to see a 1D saturated/unsaturated water flow problem.           
      15 - to see a 1D version of the Schrodinger eigenvalue equation      
           of example 11.                                                  
       0 - (no example)                                                    
                                                                           
   If you select one of the examples, the correct answer for each          
   interactive question will be supplied after the question.               
|---- Enter an integer value in the range 0 to 15                 
10                                                                      
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ The example problem is                                               $
   $                                                                      $
   $         Ut = V                  in (-1,1) X (-1,1) X (-1,1)          $
   $         Vt = Uxx + Uyy + Uzz    for 0 < T < 0.6                      $
   $ with boundary conditions                                             $
   $         Ux = 0   Vx = 0         on the faces X=-1 and X=1            $
   $         Uy = 0   Vy = 0         on the faces Y=-1 and Y=1            $
   $         Uz = 0   Vz = 0         on the faces Z=-1 and Z=1            $
   $ and initial conditions                                               $
   $         U = exp(-10*(X**2+Y**2+Z**2))      at T = 0                  $
   $         V = 0                                                        $
   $                                                                      $
   $ (Note: this models the wave equation Utt = Uxx + Uyy + Uzz, reduced  $
   $ to a system of two equations by the introduction of the velocity     $
   $ variable V=Ut.)                                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
                                                                      [RETURN]
   In what follows, when you are told to enter a 'FORTRAN expression',     
   this means any valid FORTRAN expression of 65 characters or less.       
   In this expression, you may include references to FORTRAN function      
   subprograms.  You may define these functions line by line at the end    
   of the interactive session, when prompted, or add them later using an   
   editor.                                                                 
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you enter a "#" in the first column of any input line, this       +
   + instructs the interactive driver to read this and subsequent input   +
   + lines from the file "pde2d.in".  A "#" in the first column of an     +
   + input line in the file "pde2d.in" (or an end-of-file) instructs the  +
   + driver to switch back to interactive input.                          +
   +                                                                      +
   + All lines input during an interactive session are echo printed to    +
   + a file "echo.out".  You may want to modify this file and rename it   +
   + "pde2d.in", and read some or all of your input from this file during +
   + your next interactive session.                                       +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
                                                                           
   PDE2D can solve problems with 0,1,2 or 3 space variables.  Enter the    
   dimension of your problem:                                              
                                                                      [RETURN]
                                                                           
     0 - to solve a time-dependent ordinary differential equation system,  
         or an algebraic or algebraic eigenvalue system                    
     1 - to solve problems in 1D intervals                                 
     2 - to solve problems in general 2D regions                           
     3 - to solve problems in a wide range of simple 3D regions            
                                                                           
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: 3                                                             $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter an integer value in the range 0 to 3                  
3                                                                       
   Is double precision mode to be used?  Double precision is recommended.  
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If double precision mode is used, variables and functions assigned   +
   + names beginning with a letter in the range A-H or O-Z will be DOUBLE +
   + PRECISION, and you should use double precision constants and FORTRAN +
   + expressions throughout; otherwise such variables and functions will  +
   + be of type REAL.  In either case, variables and functions assigned   +
   + names beginning with I,J,K,L,M or N will be of INTEGER type.         +
   +                                                                      +
   + It is possible to convert a single precision PDE2D program to double +
   + precision after it has been created, using an editor.  Just change   +
   + all occurrences of "real" to "double precision"                      +
   +                    " tdp" to "dtdp"  (note leading blank)            +
   + Any user-written code or routines must be converted "by hand", of    +
   + course.  To convert from double to single, reverse the changes.      +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: yes                                                           $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
yes
   If you don't want to read the FINE PRINT, default NPROB.                
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you want to solve several similar problems in the same run, set   +
   + NPROB equal to the number of problems you want to solve.  Then NPROB +
   + loops through the main program will be done, with IPROB=1,...,NPROB, +
   + and you can make the problem parameters vary with IPROB.  NPROB      +
   + defaults to 1.                                                       +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ press [RETURN] to default NPROB                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NPROB =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
   PDE2D solves the time-dependent system (note: U,F,G,U0 may be vectors,  
   C,RHO may be matrices):                                                 
                                                                           
      C(X,Y,Z,T,U,Ux,Uy,Uz)*d(U)/dT =                                      
                          F(X,Y,Z,T,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz)    
                                                                           
   or the steady-state system:                                             
                                                                           
      F(X,Y,Z,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz) = 0                      
                                                                           
   or the linear and homogeneous eigenvalue system:                        
                                                                           
      F(X,Y,Z,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz) = lambda*RHO(X,Y,Z)*U    
                                                                           
   with boundary conditions:                                               
                                                                           
                G(X,Y,Z,[T],U,Ux,Uy,Uz) = 0                                
         (periodic boundary conditions are also permitted)                 
                                                                           
   For time-dependent problems there are also initial conditions:          
                                                                           
          U = U0(X,Y,Z)   at T=T0                                          
                                                                      [RETURN]
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If your PDEs involve the solution at points other than (P1,P2,P3),   +
   + the function                                                         +
   +             (D)OLDSOL3(IDER,IEQ,PP1,PP2,PP3,KDEG)                    +
   + will interpolate (using interpolation of degree KDEG=1,2 or 3) to    +
   + (PP1,PP2,PP3) the function saved in UOUT(*,*,*,IDER,IEQ,ISET) on the +
   + last time step or iteration (ISET) for which it has been saved.      +
   + Thus, for example, if IDER=1, this will return the latest value of   +
   + component IEQ of the solution at (PP1,PP2,PP3), assuming this has    +
   + not been modified using UPRINT... If your equations involve          +
   + integrals of the solution, for example, you can use (D)OLDSOL3 to    +
   + approximate these using the solution from the last time step or      +
   + iteration.                                                           +
   +                                                                      +
   + CAUTION: For a steady-state or eigenvalue problem, you must reset    +
   + NOUT=1 if you want to save the solution each iteration.              +
   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   + A system of NEQN complex partial differential equations must be      +
   + written as a system of 2*NEQN real equations, by separating the      +
   + equations into their real and imaginary parts.  However, note that   +
   + the complex arithmetic abilities of FORTRAN can be used to simplify  +
                                                                      [RETURN]
   + this separation.  For example, the complex PDE:                      +
   +      I*(Uxx+Uyy+Uzz) - 1/(1+U**10) = 0,   where U = UR + UI*I        +
   + would be difficult to split up analytically, but using FORTRAN       +
   + expressions it is easy:                                              +
   +   F1 = -(UIxx+UIyy+UIzz) -  REAL(1.0/(1.0+CMPLX(UR,UI)**10))         +
   +   F2 =  (URxx+URyy+URzz) - AIMAG(1.0/(1.0+CMPLX(UR,UI)**10))         +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   You may now define global parameters, which may be referenced in any    
   of the "FORTRAN expressions" you input throughout the rest of this      
   interactive session.  You will be prompted alternately for parameter    
   names and their values; enter a blank name when you are finished.       
                                                                           
   Parameter names are valid FORTRAN variable names, starting in           
   column 1.  Thus each name consists of 1 to 6 alphanumeric characters,   
   the first of which must be a letter.  If the first letter is in the     
   range I-N, the parameter must be an integer.                            
                                                                           
   Parameter values are either FORTRAN constants or FORTRAN expressions    
   involving only constants and global parameters defined on earlier       
   lines.  They may also be functions of the problem number IPROB, if      
   you are solving several similar problems in one run (NPROB > 1).  Note  
   that you are defining global CONSTANTS, not functions; e.g., parameter  
                                                                      [RETURN]
   values may not reference any of the independent or dependent variables  
   of your problem.                                                        
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you define other parameters here later, using an editor, you must +
   + add them to COMMON block /PARM8Z/ everywhere this block appears, if  +
   + they are to be "global" parameters.                                  +
   +                                                                      +
   + The variable PI is already included as a global parameter, with an   +
   + accurate value 3.14159...                                            +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ press [RETURN] to skip global parameter definitions                  $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 Parameter name = (type blank line to terminate)
       
   If you don't want to read the FINE PRINT, enter 'no'.                   
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + Do you want to be given a chance to write a FORTRAN block before the +
   + definitions of many functions?  If you answer 'no', you will still   +
   + be given a chance to write code before the definition of the PDE     +
   + coefficients, but not other functions.  Of course, you can always    +
   + add code later directly to the resulting program, using an editor.   +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: no                                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
no 
   You can solve problems in your region only if you can describe it by    
                       X = X(P1,P2,P3)                                     
                       Y = Y(P1,P2,P3)                                     
                       Z = Z(P1,P2,P3)                                     
   with constant limits on the parameters P1,P2,P3.  If your region is     
   rectangular, enter ITRANS=0 and the trivial parameterization            
                       X = P1                                              
                       Y = P2                                              
                       Z = P3                                              
   will be used.  Otherwise, you need to read the FINE PRINT below.        
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If P1,P2,P3 represent cylindrical, spherical or other non-Cartesian  +
   + coordinates, you can reference the Cartesian coordinates X,Y,Z       +
   + and derivatives of your unknowns with respect to these coordinates,  +
   + when you define your PDE coefficients, boundary conditions, and      +
   + volume and boundary integrals, if you enter ITRANS .NE. 0.  Enter:   +
   +   ITRANS = 1, if P1,P2,P3 are cylindrical coordinates, that is, if   +
   +               P1=R, P2=Theta, P3=Z, where X = R*cos(Theta)           +
   +                                           Y = R*sin(Theta)           +
   +                                           Z = Z                      +
   +   ITRANS = -1, same as ITRANS=1, but P1=Theta, P2=R, P3=Z            +
                                                                      [RETURN]
   +   ITRANS = 2, if P1,P2,P3 are spherical coordinates, that is, if     +
   +               P1=Rho, P2=Phi, P3=Theta, where                        +
   +                                    X = Rho*sin(Phi)*cos(Theta)       +
   +                                    Y = Rho*sin(Phi)*sin(Theta)       +
   +                                    Z = Rho*cos(Phi)                  +
   +               (Theta is longitude, Phi is measured from north pole)  +
   +   ITRANS = -2, same as ITRANS=2, but P1=Rho, P2=Theta, P3=Phi        +
   +   ITRANS = 3, to define your own coordinate transformation.  In this +
   +               case, you will be prompted to define X,Y,Z and their   +
   +               first and second derivatives in terms of P1,P2,P3.     +
   +               Because of symmetry, you will not be prompted for all  +
   +               of the second derivatives.  If you make a mistake in   +
   +               computing any of these derivatives, PDE2D will usually +
   +               be able to issue a warning message. (X1 = dX/dP1, etc) +
   +   ITRANS = -3, same as ITRANS=3, but you will only be prompted to    +
   +               define X,Y,Z; their first and second derivatives will  +
   +               be approximated using finite differences.              +
   +   When ITRANS = -3 or 3, the first derivatives of X,Y,Z must all be  +
   +   continuous.                                                        +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: ITRANS = 0                                                    $
                                                                      [RETURN]
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ITRANS =
|---- Enter an integer value in the range -3 to 3                 
0                                                                       
   A collocation finite element method is used, with tri-cubic Hermite     
   basis functions on the elements (small boxes) defined by the grid       
   points:                                                                 
             P1GRID(1),...,P1GRID(NP1GRID)                                 
             P2GRID(1),...,P2GRID(NP2GRID)                                 
             P3GRID(1),...,P3GRID(NP3GRID)                                 
   You will first be prompted for NP1GRID, the number of P1-grid points,   
   then for P1GRID(1),...,P1GRID(NP1GRID).  Any points defaulted will be   
   uniformly spaced between the points you define; the first and last      
   points cannot be defaulted.  Then you will be prompted similarly        
   for the number and values of the P2 and P3-grid points.  The limits     
   on the parameters are then:                                             
             P1GRID(1) < P1 < P1GRID(NP1GRID)                              
             P2GRID(1) < P2 < P2GRID(NP2GRID)                              
             P3GRID(1) < P3 < P3GRID(NP3GRID)                              
                                                                           
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: NP1GRID = 7                                                   $
   $        P1GRID(1) = -1                                                $
   $        P1GRID(NP1GRID) = 1                                           $
   $ and default P1GRID(2),P1GRID(3),P1GRID(4),P1GRID(5),P1GRID(6)        $
   $        NP2GRID = 7                                                   $
                                                                      [RETURN]
   $        P2GRID(1) = -1                                                $
   $        P2GRID(NP2GRID) = 1                                           $
   $ and default P2GRID(2),P2GRID(3),P2GRID(4),P2GRID(5),P2GRID(6)        $
   $        NP3GRID = 7                                                   $
   $        P3GRID(1) = -1                                                $
   $        P3GRID(NP3GRID) = 1                                           $
   $ and default P3GRID(2),P3GRID(3),P3GRID(4),P3GRID(5),P3GRID(6)        $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NP1GRID =
|---- Enter an integer value in the range 1 to +INFINITY          
7                                                                       
  P1GRID(1) =
|----Enter constant or FORTRAN expression-----------------------|
-1                                                               
  P1GRID(2) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P1GRID(3) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P1GRID(4) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P1GRID(5) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P1GRID(6) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P1GRID(NP1GRID) =
|----Enter constant or FORTRAN expression-----------------------|
1                                                                
  NP2GRID =
|---- Enter an integer value in the range 1 to +INFINITY          
7                                                                       
  P2GRID(1) =
|----Enter constant or FORTRAN expression-----------------------|
-1                                                               
  P2GRID(2) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P2GRID(3) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P2GRID(4) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P2GRID(5) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P2GRID(6) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P2GRID(NP2GRID) =
|----Enter constant or FORTRAN expression-----------------------|
1                                                                
  NP3GRID =
|---- Enter an integer value in the range 1 to +INFINITY          
7                                                                       
  P3GRID(1) =
|----Enter constant or FORTRAN expression-----------------------|
-1                                                               
  P3GRID(2) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P3GRID(3) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P3GRID(4) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P3GRID(5) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P3GRID(6) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P3GRID(NP3GRID) =
|----Enter constant or FORTRAN expression-----------------------|
1                                                                
   If you don't want to read the FINE PRINT, enter ISOLVE = 1.             
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + The following linear system solvers are available:                   +
   +                                                                      +
   + 1. Sparse direct method                                              +
   +               Harwell Library routine MA27 (used by permission) is   +
   +               used to solve the (positive definite) "normal"         +
   +               equations A**T*A*x = A**T*b.  The normal equations,    +
   +               which are essentially the equations which would result +
   +               if a least squares finite element method were used     +
   +               instead of a collocation method, are substantially     +
   +               more ill-conditioned than the original system Ax = b,  +
   +               so it may be important to use high precision if this   +
   +               option is chosen.                                      +
   + 2. Frontal method                                                    +
   +               This is an out-of-core band solver.  If you want to    +
   +               override the default number of rows in the buffer (11),+
   +               set a new value for NPMX8Z in the main program.        +
   + 3. Jacobi conjugate gradient iterative method                        +
   +               A preconditioned conjugate gradient iterative method   +
   +               is used to solve the (positive definite) normal        +
                                                                      [RETURN]
   +               equations.  High precision is also important if this   +
   +               option is chosen.  (This solver is MPI-enhanced, if    +
   +               MPI is available.)  If you want to override the        +
   +               default convergence tolerance, set a new relative      +
   +               tolerance CGTL8Z in the main program.                  +
   + 4. Local solver (normal equations)                                   +
   + 5. Local solver (original equations)                                 +
   +               Choose these options ONLY if alterative linear system  +
   +               solvers have been installed locally.  See subroutines  +
   +               (D)TD3M, (D)TD3N in file (d)subs.f for instructions    +
   +               on how to add local solvers.                           +
   + 6. MPI-based parallel band solver                                    +
   +               This is a parallel solver which runs efficiently on    +
   +               multiple processor machines, under MPI.  It is a       +
   +               band solver, with the matrix distributed over the      +
   +               available processors.  Choose this option ONLY if the  +
   +               solver has been activated locally.  See subroutine     +
   +               (D)TD3O in file (d)subs.f for instructions on how to   +
   +               activate this solver and the MPI-enhancements to the   +
   +               conjugate gradient solver.                             +
   +                                                                      +
   + Enter ISOLVE = 1,2,3,4,5 or 6 to select a linear system solver.      +
                                                                      [RETURN]
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
                                                                           
   If you don't want to read the FINE PRINT, enter ISOLVE = 1.             
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: ISOLVE = 1                                                    $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ISOLVE =
|---- Enter an integer value in the range 1 to 6                  
1                                                                       
   What type of PDE problem do you want to solve?                          
                                                                           
      1. a steady-state (time-independent) problem                         
      2. a time-dependent problem                                          
      3. a linear, homogeneous eigenvalue problem                          
                                                                           
   Enter 1,2 or 3 to select a problem type.                                
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: 2                                                             $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter an integer value in the range 1 to 3                  
2                                                                       
   Enter the initial time value (T0) and the final time value (TF), for    
   this time-dependent problem.  T0 defaults to 0.                         
                                                                           
   TF is not required to be greater than T0.                               
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: T0 = 0                                                        $
   $        TF = 0.6                                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  T0 =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
  TF =
|----Enter constant or FORTRAN expression-----------------------|
0.6                                                              
   Is this a linear problem? ("linear" means all differential equations    
   and all boundary conditions are linear).  If you aren't sure, it is     
   safer to answer "no".                                                   
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: yes                                                           $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
yes
   Do you want the time step to be chosen adaptively?  If you answer       
   'yes', you will then be prompted to enter a value for TOLER(1), the     
   local relative time discretization error tolerance.  The default is     
   TOLER(1)=0.01.  If you answer 'no', a user-specified constant time step 
   will be used.  We suggest that you answer 'yes' and default TOLER(1)    
   (although for certain linear problems, a constant time step may be much 
   more efficient).                                                        
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If a negative value is specified for TOLER(1), then ABS(TOLER(1)) is +
   + taken to be the "absolute" error tolerance.  If a system of PDEs is  +
   + solved, by default the error tolerance specified in TOLER(1) applies +
   + to all variables, but the error tolerance for the J-th variable can  +
   + be set individually by specifying a value for TOLER(J) using an      +
   + editor, after the end of the interactive session.                    +
   +                                                                      +
   + Each time step, two steps of size dt/2 are taken, and that solution  +
   + is compared with the result when one step of size dt is taken.  If   +
   + the maximum difference between the two answers is less than the      +
   + tolerance (for each variable), the time step dt is accepted (and the +
   + next step dt is doubled, if the agreement is "too" good); otherwise  +
   + dt is halved and the process is repeated.  Note that forcing the     +
                                                                      [RETURN]
   + local (one-step) error to be less than the tolerance does not        +
   + guarantee that the global (cumulative) error is less than that value.+
   + However, as the tolerance is decreased, the global error should      +
   + decrease correspondingly.                                            +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: no                                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
no 
   If you don't want to read the FINE PRINT, it is safe (though possibly   
   very inefficient) to enter 'no'.                                        
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If your time-dependent problem is linear with all PDE and boundary   +
   + condition coefficients independent of time except inhomogeneous      +
   + terms, then a large savings in execution time may be possible if     +
   + this is recognized (the LU decomposition computed on the first step  +
   + can be used on subsequent steps).  Is this the case for your         +
   + problem?  (Caution: if you answer 'yes' when you should not, you     +
   + will get incorrect results with no warning.)                         +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: yes                                                           $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
yes
   The time stepsize will be constant, DT = (TF-T0)/NSTEPS.  Enter a       
   value for NSTEPS, the number of time steps.                             
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you later turn on adaptive step control, the time stepsize will be+
   + chosen adaptively, between an upper limit of DTMAX = (TF-T0)/NSTEPS  +
   + and a lower limit of 0.0001*DTMAX.                                   +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ So that DT = (0.6-0.0)/30 = 0.02,                                    $
   $ enter: NSTEPS = 30                                                   $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NSTEPS =
|----Enter constant or FORTRAN expression-----------------------|
30                                                               
   If you don't want to read the FINE PRINT, enter 'no'.                   
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + Is the Crank-Nicolson scheme to be used to discretize time?  If you  +
   + answer 'no', a backward Euler scheme will be used.                   +
   +                                                                      +
   + If a user-specified constant time step is chosen, the second order   +
   + Crank Nicolson method is recommended only for problems with very     +
   + well-behaved solutions, and the first order backward Euler scheme    +
   + should be used for more difficult problems.  In particular, do not   +
   + use the Crank Nicolson method if the left hand side of any PDE is    +
   + zero, for example, if a mixed elliptic/parabolic problem is solved.  +
   +                                                                      +
   + If adaptive time step control is chosen, however, an extrapolation   +
   + is done between the 1-step and 2-step answers which makes the Euler  +
   + method second order, and the Crank-Nicolson method strongly stable.  +
   + Thus in this case, both methods have second order accuracy, and both +
   + are strongly stable.                                                 +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: yes                                                           $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
                                                                      [RETURN]
|---- Enter yes or no
yes
   How many differential equations (NEQN) are there in your problem?       
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: NEQN = 2                                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NEQN =
|---- Enter an integer value in the range 1 to 99                 
2                                                                       
   You may now choose names for the component(s) of the (possibly vector)  
   solution U.  Each must be an alphanumeric string of one to three        
   characters, beginning with a letter in the range A-H or O-Z.  The       
   variable names X,Y,Z,P and T must not be used.  The name should         
   start in column 1.                                                      
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: U1 = U                                                        $
   $        U2 = V                                                        $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 U1 =
U   
 U2 =
V   
   You may calculate one or more integrals (over the entire region) of     
   some functions of the solution and its derivatives.  How many integrals 
   (NINT), if any, do you want to calculate?                               
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + In the FORTRAN program created by the preprocessor, the computed     +
   + values of the integrals will be returned in the vector SINT8Z.  If   +
   + several iterations or time steps are done, only the last computed    +
   + values are saved in SINT8Z (all values are printed).                 +
   +                                                                      +
   + A limiting value, SLIM8Z(I), for the I-th integral can be set        +
   + below in the main program.  The computations will then stop          +
   + gracefully whenever SINT8Z(I) > SLIM8Z(I), for any I=1...NINT.       +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: NINT = 1                                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NINT =
|---- Enter an integer value in the range 0 to 20                 
1                                                                       
   Enter FORTRAN expressions for the functions whose integrals are to be   
   calculated and printed.  They may be functions of                       
                                                                           
      X,Y,Z,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz                             
            V,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz                             
                                              and (if applicable) T        
                                                                           
   The parameters P1,P2,P3 and derivatives with respect to these may also  
   be referenced (U1 = dU/dP1, etc):                                       
         U1,U2,U3,U11,U22,U33,U12,U13,U23                                  
         V1,V2,V3,V11,V22,V33,V12,V13,V23                                  
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you only want to integrate a function over part of the region,    +
   + define that function to be zero in the rest of the region.           +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want to compute the total energy in the box, the integral of      $
   $ Ut**2 + Ux**2 + Uy**2 + Uz**2 (which should be constant), so         $
   $ enter: INTEGRAL = V**2+Ux**2+Uy**2+Uz**2                             $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 INTEGRAL =
|----Enter constant or FORTRAN expression-----------------------|
V**2+Ux**2+Uy**2+Uz**2                                           
   You may calculate one or more boundary integrals (over the entire       
   boundary) of some functions of the solution and its derivatives.  How   
   many boundary integrals (NBINT), if any, do you want to calculate?      
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + In the FORTRAN program created by the preprocessor, the computed     +
   + values of the integrals will be returned in the vector BINT8Z.  If   +
   + several iterations or time steps are done, only the last computed    +
   + values are saved in BINT8Z (all values are printed).                 +
   +                                                                      +
   + A limiting value, BLIM8Z(I), for the I-th boundary integral can be   +
   + set below in the main program.  The computations will then stop      +
   + gracefully whenever BINT8Z(I) > BLIM8Z(I), for any I=1...NBINT.      +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: NBINT = 0                                                     $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NBINT =
|---- Enter an integer value in the range 0 to 20                 
0                                                                       
   Now enter FORTRAN expressions to define the PDE coefficients, which     
   may be functions of                                                     
                                                                           
      X,Y,Z,T,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz                           
              V,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz                           
                                                                           
   Recall that the PDEs have the form                                      
                                                                           
                C11*d(U)/dT + C12*d(V)/dT = F1                             
                C21*d(U)/dT + C22*d(V)/dT = F2                             
                                                                           
   The parameters P1,P2,P3 and derivatives with respect to these may also  
   be referenced (U1 = dU/dP1, etc):                                       
         U1,U2,U3,U11,U22,U33,U12,U13,U23                                  
         V1,V2,V3,V11,V22,V33,V12,V13,V23                                  
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ The partial differential equations may be written                    $
   $                                                                      $
   $        1.0*Ut = V                                                    $
   $        1.0*Vt = Uxx + Uyy + Uzz                                      $
   $                                                                      $
   $ Thus when asked if you want to write a FORTRAN block,                $
                                                                      [RETURN]
   $ enter: no                                                            $
   $   then, when prompted, enter:                                        $
   $   C11 = 1    C12 = 0    F1 = V                                       $
   $   C21 = 0    C22 = 1    F2 = Uxx + Uyy + Uzz                         $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

  Do you want to write a FORTRAN block to define some parameters to be
   used in the definition of these coefficients?
|---- Enter yes or no
no 
 C(1,1) =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
1                                                                
 C(1,2) =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
 F1 =
|----Enter constant or FORTRAN expression-----------------------|
V                                                                
 C(2,1) =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
 C(2,2) =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
1                                                                
 F2 =
|----Enter constant or FORTRAN expression-----------------------|
Uxx + Uyy + Uzz                                                  
   Now the initial values must be defined using FORTRAN expressions.       
   They may be functions of X,Y and Z (and the parameters P1,P2,P3), and   
   may also reference the initial time T0.                                 
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: U0 = EXP(-10*(X**2+Y**2+Z**2))                                $
   $        V0 = 0                                                        $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 U0 =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
EXP(-10*(X**2+Y**2+Z**2))                                        
 V0 =         (Press [RETURN] to default to 0)
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
   If you don't want to read the FINE PRINT, enter 'no'.                   
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + Normally, interpolation is done to approximate the initial values    +
   + using cubic Hermites.  Since some derivatives must be interpolated,  +
   + if the initial values are not smooth (ie, have large or infinite     +
   + derivatives), the resulting cubic interpolants may have undesired    +
   + noise or large spikes.  Do you want to compute a least squares       +
   + approximation to the initial values, rather than an interpolant?     +
   + The least squares fit is generally much smoother, but requires one   +
   + extra linear system solution.                                        +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: yes                                                           $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
yes
   If you don't want to read the FINE PRINT, enter 'no'.                   
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + Do you want to read the initial conditions from the restart file,    +
   + if it exists (and use the conditions supplied above if it does not   +
   + exist)?                                                              +
   +                                                                      +
   + If so, PDE2D will dump the final solution at the end of each run     +
   + into a restart file "pde2d.res".  Thus the usual procedure for       +
   + using this dump/restart option is to make sure there is no restart   +
   + file in your directory left over from a previous job, then the       +
   + first time you run this job, the initial conditions supplied above   +
   + will be used, but on the second and subsequent runs the restart file +
   + from the previous run will be used to define the initial conditions. +
   +                                                                      +
   + You can do all the "runs" in one program, by setting NPROB > 1.      +
   + Each pass through the DO loop, T0,TF,NSTEPS and possibly other       +
   + parameters may be varied, by making them functions of IPROB.         +
   +                                                                      +
   + If the 2D or 3D collocation method is used, the coordinate           +
   + transformation should not change between dump and restart.           +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
                                                                      [RETURN]
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: no                                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
no 
   If you do not have any periodic boundary conditions, enter IPERDC=0.    
                                                                           
   Enter IPERDC=1 for periodic conditions at P1 = P1GRID(1),P1GRID(NP1GRID)
         IPERDC=2 for periodic conditions at P2 = P2GRID(1),P2GRID(NP2GRID)
         IPERDC=3 for periodic conditions at P3 = P3GRID(1),P3GRID(NP3GRID)
         IPERDC=4 for periodic conditions on both P1 and P2                
         IPERDC=5 for periodic conditions on both P1 and P3                
         IPERDC=6 for periodic conditions on both P2 and P3                
         IPERDC=7 for periodic conditions on P1, P2 and P3.                
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + When periodic boundary conditions are selected, they apply to all    +
   + variables by default.  To turn off periodic boundary conditions on   +
   + the I-th variable, set PERDC(I) to 0 (or another appropriate value   +
   + of IPERDC) below in the main program and set the desired boundary    +
   + conditions in subroutine GB8Z, "by hand".                            +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: IPERDC = 0                                                    $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  IPERDC =
|---- Enter an integer value in the range 0 to 7                  
0                                                                       
   Enter FORTRAN expressions to define the boundary condition functions,   
   which may be functions of                                               
                                                                           
            X,Y,Z,U,Ux,Uy,Uz,                                              
                  V,Vx,Vy,Vz and (if applicable) T                         
                                                                           
   Recall that the boundary conditions have the form                       
                                                                           
                         G1 = 0                                            
                         G2 = 0                                            
                                                                           
   Enter NONE to indicate "no" boundary condition.                         
                                                                           
   The parameters P1,P2,P3 and derivatives with respect to these may also  
   be referenced (U1 = dU/dP1, etc):                                       
                            U1,U2,U3                                       
                            V1,V2,V3                                       
   The components (NORMx,NORMy,NORMz) of the unit outward normal vector    
   may also be referenced, as well as the normal derivatives Unorm,        
   Vnorm.                                                                  
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
                                                                      [RETURN]
   + If "no" boundary condition is specified, the corresponding PDE is    +
   + enforced at points just inside the boundary (exactly on the          +
   + boundary, if EPS8Z is set to 0 in the main program).                 +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
                                                                           
   First define the boundary conditions on the face P1 = P1GRID(1).        
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: G1 = Ux                                                       $
   $        G2 = Vx                                                       $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 G1 =
|----Enter constant or FORTRAN expression-----------------------|
Ux                                                               
 G2 =
|----Enter constant or FORTRAN expression-----------------------|
Vx                                                               
                                                                           
   Now define the boundary conditions on the face P1 = P1GRID(NP1GRID).    
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: G1 = Ux                                                       $
   $        G2 = Vx                                                       $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 G1 =
|----Enter constant or FORTRAN expression-----------------------|
Ux                                                               
 G2 =
|----Enter constant or FORTRAN expression-----------------------|
Vx                                                               
                                                                           
   Now define the boundary conditions on the face P2 = P2GRID(1).          
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: G1 = Uy                                                       $
   $        G2 = Vy                                                       $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 G1 =
|----Enter constant or FORTRAN expression-----------------------|
Uy                                                               
 G2 =
|----Enter constant or FORTRAN expression-----------------------|
Vy                                                               
                                                                           
   Now define the boundary conditions on the face P2 = P2GRID(NP2GRID).    
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: G1 = Uy                                                       $
   $        G2 = Vy                                                       $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 G1 =
|----Enter constant or FORTRAN expression-----------------------|
Uy                                                               
 G2 =
|----Enter constant or FORTRAN expression-----------------------|
Vy                                                               
                                                                           
   Now define the boundary conditions on the face P3 = P3GRID(1).          
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: G1 = Uz                                                       $
   $        G2 = Vz                                                       $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 G1 =
|----Enter constant or FORTRAN expression-----------------------|
Uz                                                               
 G2 =
|----Enter constant or FORTRAN expression-----------------------|
Vz                                                               
                                                                           
   Now define the boundary conditions on the face P3 = P3GRID(NP3GRID).    
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: G1 = Uz                                                       $
   $        G2 = Vz                                                       $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 G1 =
|----Enter constant or FORTRAN expression-----------------------|
Uz                                                               
 G2 =
|----Enter constant or FORTRAN expression-----------------------|
Vz                                                               
   If you don't want to read the FINE PRINT, default all of the following  
   variables.                                                              
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + Normally, PDE2D saves the values of U,Ux,Uy,Uz,V,Vx,                 +
   + Vy,Vz at the output points.  If different variables are to be        +
   + saved (for later printing or plotting) the following functions can   +
   + be used to re-define the output variables:                           +
   +    define UPRINT(1) to replace  U                                    +
   +           UXPRINT(1)            Ux                                   +
   +           UYPRINT(1)            Uy                                   +
   +           UZPRINT(1)            Uz                                   +
   +           UPRINT(2)             V                                    +
   +           UXPRINT(2)            Vx                                   +
   +           UYPRINT(2)            Vy                                   +
   +           UZPRINT(2)            Vz                                   +
   + Each function may be a function of                                   +
   +                                                                      +
   +    X,Y,Z,U,Ux,Uy,Uz,Uxx,Uyy,Uzz,Uxy,Uxz,Uyz                          +
   +          V,Vx,Vy,Vz,Vxx,Vyy,Vzz,Vxy,Vxz,Vyz                          +
   +                                            and (if applicable) T     +
   +                                                                      +
                                                                      [RETURN]
   + Each may also be a function of the integral estimates SINT(1),...,   +
   + BINT(1),...                                                          +
   +                                                                      +
   + The parameters P1,P2,P3 and derivatives with respect to these may    +
   + also be referenced (U1 = dU/dP1, etc):                               +
   +       U1,U2,U3,U11,U22,U33,U12,U13,U23                               +
   +       V1,V2,V3,V11,V22,V33,V12,V13,V23                               +
   +                                                                      +
   + The default for each variable is no change, for example, UPRINT(1)   +
   + defaults to U.  Enter FORTRAN expressions for each of the            +
   + following functions (or default).                                    +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want to plot the magnitude of U rather than U, so                 $
   $ enter: UPRINT(1) = ABS(U)                                            $
   $        and default the other output modification variables           $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      Replace U for postprocessing?
  UPRINT(1) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
ABS(U)                                                           
      Replace Ux for postprocessing?
  UXPRINT(1) =            (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
      Replace Uy for postprocessing?
  UYPRINT(1) =            (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
      Replace Uz for postprocessing?
  UZPRINT(1) =            (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
      Replace V for postprocessing?
  UPRINT(2) =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
      Replace Vx for postprocessing?
  UXPRINT(2) =            (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
      Replace Vy for postprocessing?
  UYPRINT(2) =            (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
      Replace Vz for postprocessing?
  UZPRINT(2) =            (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
   The solution is normally saved on an NP1+1 by NP2+1 by NP3+1            
   rectangular grid of points,                                             
                  P1 = P1A + I*(P1B-P1A)/NP1,    I = 0,...,NP1             
                  P2 = P2A + J*(P2B-P2A)/NP2,    J = 0,...,NP2             
                  P3 = P3A + K*(P3B-P3A)/NP3,    K = 0,...,NP3             
   Enter values for NP1, NP2 and NP3.  Suggested values: NP1=NP2=NP3=16.   
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you want to save the solution at an arbitrary user-specified set  +
   + of points, set NP2=NP3=0 and NP1+1=number of points.  In this case   +
   + you can request tabular output, but no plots can be made.            +
   +                                                                      +
   + If you set NEAR8Z=1 in the main program, the values saved at each    +
   + output point will actually be the solution as evaluated at a nearby  +
   + collocation point.  For most problems this obviously will produce    +
   + less accurate output or plots, but for certain (rare) problems, a    +
   + solution component may be much less noisy when plotted only at       +
   + collocation points.                                                  +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: NP1 = 16                                                      $
   $        NP2 = 16                                                      $
                                                                      [RETURN]
   $        NP3 = 16                                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NP1 =
|---- Enter an integer value in the range 1 to +INFINITY          
16                                                                      
  NP2 =
|---- Enter an integer value in the range 0 to +INFINITY          
16                                                                      
  NP3 =
|---- Enter an integer value in the range 1 to +INFINITY          
16                                                                      
   The solution is saved on an NP1+1 by NP2+1 by NP3+1 rectangular grid    
   covering the box (P1A,P1B) x (P2A,P2B) x (P3A,P3B).  Enter values for   
   P1A,P1B,P2A,P2B,P3A,P3B.  These variables are usually defaulted.        
                                                                           
   The defaults are P1A = P1GRID(1), P1B = P1GRID(NP1GRID)                 
                    P2A = P2GRID(1), P2B = P2GRID(NP2GRID)                 
                    P3A = P3GRID(1), P3B = P3GRID(NP3GRID)                 
                                                                           
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ press [RETURN] to default P1A,P1B,P2A,P2B,P3A,P3B                    $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  P1A =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P1B =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P2A =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P2B =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P3A =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  P3B =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
   The solution will be saved (for possible postprocessing) at the NSAVE+1 
   time points                                                             
                        T0 + K*(TF-T0)/NSAVE                               
   K=0,...,NSAVE.  Enter a value for NSAVE.                                
                                                                           
   If a user-specified constant time step is used, NSTEPS must be an       
   integer multiple of NSAVE.                                              
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ To save the solution at the 31 points T = 0.00, 0.02, ..., 0.60,     $
   $ enter: NSAVE = 30                                                    $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NSAVE =
|---- Enter an integer value in the range 1 to +INFINITY          
30                                                                      
   To generate tabular or graphical output, choose an output option from   
   the list below.                                                         
                                                                           
   0. No further output is desired                                         
   1. Table of values at output points                                     
         The tabulated output is saved in a file.                          
   2. A plot of the contour (level) surfaces of a scalar variable.         
         This plot gives only a coarse overview of the solution, but is    
         unique in that it shows the general form of the solution in a     
         single picture.  The different levels are identified by color.    
         These plots will only reflect the true geometry if ITRANS=0.      
   3. Surface plot of a scalar variable, at P1, P2 or P3 = constant        
         cross-sections.                                                   
   4. Contour plot of a scalar variable, at P1, P2 or P3 = constant        
         cross-sections.  These plots can be made to reflect the true      
         geometry of the cross-section.                                    
   5. Vector field plot, at P1, P2 or P3 = constant cross-sections.        
         A 3D vector field is plotted using arrows to represent the        
         in-plane components, with a contour plot of the out-of-plane      
         component superimposed.  These plots can be made to reflect       
         the true geometry of the cross-section.                           
   6. One dimensional cross-sectional plots (versus P1, P2, P3 or T)       
                                                                      [RETURN]
                                                                           
   Enter 0,1,2,3,4,5 or 6 to select an output option.                      
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you set NP2=NP3=0 and saved the solution at an arbitrary set of   +
   + user-specified output points, you can only request tabular output.   +
   +                                                                      +
   + If you decide later that you want additional types of plots not      +
   + requested during this interactive session, you will have to work     +
   + through a new interactive session, so it is recommended that you     +
   + request all output or plots you think you MIGHT eventually want now, +
   + during this session.                                                 +
   +                                                                      +
   + Regardless of the options you select, a dummy subroutine POSTPR      +
   + will be included in the program created by the interactive driver;   +
   + you can add your own postprocessing code to this subroutine.         +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want three types of plots for ABS(U): level surface plots,        $
   $ surface plots at a constant Z (P3) cross-section, and a couple of    $
   $ time-plots, so                                                       $
   $ enter: 2, the first time you see this message and                    $
                                                                      [RETURN]
   $        3, the second time                                            $
   $        6, the third time                                             $
   $        0, the fourth time                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter an integer value in the range 0 to 6                  
2                                                                       
   Enter a value for IVAR, to select the variable to be plotted or         
   printed:                                                                
       IVAR = 1 means U  (possibly as modified by UPRINT,..)               
              2       Ux                                                   
              3       Uy                                                   
              4       Uz                                                   
              5       V                                                    
              6       Vx                                                   
              7       Vy                                                   
              8       Vz                                                   
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want to plot ABS(U), which has been saved in U, so                $
   $ enter: IVAR = 1                                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  IVAR =
|---- Enter an integer value in the range 1 to 8                  
1                                                                       
   If you don't want to read the FINE PRINT, default ISET1,ISET2,ISINC.    
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + The tabular output or plots will be made at times:                   +
   +        T(K) = T0 + K*(TF-T0)/NSAVE                                   +
   + for    K = ISET1, ISET1+ISINC, ISET1+2*ISINC,..., ISET2              +
   + Enter values for ISET1, ISET2 and ISINC.                             +
   +                                                                      +
   + The default is ISET1=0, ISET2=NSAVE, ISINC=1, that is, the tabular   +
   + output or plots will be made at all time values for which the        +
   + solution has been saved.                                             +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ To get plots at T=0, 0.2, 0.4 and 0.6,                               $
   $ Enter: ISET1=0, ISET2=NSAVE, ISINC=10                                $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ISET1 =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
  ISET2 =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
NSAVE                                                            
  ISINC =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
10                                                               
   Enter the view latitude, VLAT, and the view longitude, VLON, desired    
   for this plot, in degrees.  VLAT and VLON must be between 10 and 80     
   degrees; each defaults to 45 degrees.  VLAT and VLON are usually        
   defaulted.                                                              
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ press [RETURN] to default VLAT and VLON                              $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  VLAT =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  VLON =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
   Enter the number (NC) of contour surfaces to be drawn.  A maximum       
   of 10 are allowed, and more than the default of NC=4 will often         
   be too many.                                                            
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want to draw only one level surface, corresponding to mid-range   $
   $ intensity, so                                                        $
   $ enter: NC = 1                                                        $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NC =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
1                                                                
   Do you want the contour surfaces to be drawn using solid bands?         
   (If you answer 'no' they will be drawn using hollow bands.)             
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: no                                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
no 
   Enter lower (UMIN) and upper (UMAX) bounds for the contour values.      
   UMIN and UMAX are often defaulted.                                      
                                                                           
   Labeled contours will be drawn corresponding to the values              
                                                                           
            UMIN + (I-0.5)/NC*(UMAX-UMIN), for I=1,2,...,NC                
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + By default, UMIN and UMAX are set to the minimum and maximum values  +
   + of the variable to be plotted.  For a common scaling, you may want   +
   + to set UMIN=ALOW, UMAX=AHIGH.  ALOW and AHIGH are the minimum and    +
   + maximum values over all output points and over all saved time steps  +
   + or iterations.                                                       +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ press [RETURN] to default UMIN,UMAX                                  $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  UMIN =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  UMAX =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
   Enter a title, WITHOUT quotation marks.  A maximum of 40 characters     
   are allowed.  The default is no title.                                  
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: 3D wave problem                                               $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  TITLE =             (Press [RETURN] to default)
|----Enter title or name---------------|
3D wave problem                         
   To generate tabular or graphical output, choose an output option from   
   the list below.                                                         
                                                                           
   0. No further output is desired                                         
   1. Table of values at output points                                     
         The tabulated output is saved in a file.                          
   2. A plot of the contour (level) surfaces of a scalar variable.         
         This plot gives only a coarse overview of the solution, but is    
         unique in that it shows the general form of the solution in a     
         single picture.  The different levels are identified by color.    
         These plots will only reflect the true geometry if ITRANS=0.      
   3. Surface plot of a scalar variable, at P1, P2 or P3 = constant        
         cross-sections.                                                   
   4. Contour plot of a scalar variable, at P1, P2 or P3 = constant        
         cross-sections.  These plots can be made to reflect the true      
         geometry of the cross-section.                                    
   5. Vector field plot, at P1, P2 or P3 = constant cross-sections.        
         A 3D vector field is plotted using arrows to represent the        
         in-plane components, with a contour plot of the out-of-plane      
         component superimposed.  These plots can be made to reflect       
         the true geometry of the cross-section.                           
   6. One dimensional cross-sectional plots (versus P1, P2, P3 or T)       
                                                                      [RETURN]
                                                                           
   Enter 0,1,2,3,4,5 or 6 to select an output option.                      
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you set NP2=NP3=0 and saved the solution at an arbitrary set of   +
   + user-specified output points, you can only request tabular output.   +
   +                                                                      +
   + If you decide later that you want additional types of plots not      +
   + requested during this interactive session, you will have to work     +
   + through a new interactive session, so it is recommended that you     +
   + request all output or plots you think you MIGHT eventually want now, +
   + during this session.                                                 +
   +                                                                      +
   + Regardless of the options you select, a dummy subroutine POSTPR      +
   + will be included in the program created by the interactive driver;   +
   + you can add your own postprocessing code to this subroutine.         +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want three types of plots for ABS(U): level surface plots,        $
   $ surface plots at a constant Z (P3) cross-section, and a couple of    $
   $ time-plots, so                                                       $
   $ enter: 2, the first time you see this message and                    $
                                                                      [RETURN]
   $        3, the second time                                            $
   $        6, the third time                                             $
   $        0, the fourth time                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter an integer value in the range 0 to 6                  
3                                                                       
   Enter a value for IVAR, to select the variable to be plotted or         
   printed:                                                                
       IVAR = 1 means U  (possibly as modified by UPRINT,..)               
              2       Ux                                                   
              3       Uy                                                   
              4       Uz                                                   
              5       V                                                    
              6       Vx                                                   
              7       Vy                                                   
              8       Vz                                                   
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want to plot ABS(U), which has been saved in U, so                $
   $ enter: IVAR = 1                                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  IVAR =
|---- Enter an integer value in the range 1 to 8                  
1                                                                       
   Plots can be made at:                                                   
                                                                           
   1. P1 = constant cross-sections                                         
   2. P2 = constant cross-sections                                         
   3. P3 = constant cross-sections                                         
                                                                           
   Enter 1,2 or 3 to select the cross-section type.                        
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: 3                                                             $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter an integer value in the range 1 to 3                  
3                                                                       
   Plots of the variable or vector as a function of P1 and P2 will be      
   made, at the output grid P3-points closest to                           
                    P3 = P3CROSS(1),...,P3CROSS(NP3VALS)                   
                                                                           
   Enter values for NP3VALS and P3CROSS(1),...,P3CROSS(NP3VALS).           
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want one cross-section per time value, at Z(P3)=0, so             $
   $ enter: NP3VALS = 1                                                   $
   $        P3CROSS(1) = 0                                                $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NP3VALS =
|---- Enter an integer value in the range 1 to 100                
1                                                                       
  P3CROSS(1) =         
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
   If you don't want to read the FINE PRINT, default ISET1,ISET2,ISINC.    
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + The tabular output or plots will be made at times:                   +
   +        T(K) = T0 + K*(TF-T0)/NSAVE                                   +
   + for    K = ISET1, ISET1+ISINC, ISET1+2*ISINC,..., ISET2              +
   + Enter values for ISET1, ISET2 and ISINC.                             +
   +                                                                      +
   + The default is ISET1=0, ISET2=NSAVE, ISINC=1, that is, the tabular   +
   + output or plots will be made at all time values for which the        +
   + solution has been saved.                                             +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ To get plots at T=0, 0.2, 0.4 and 0.6,                               $
   $ Enter: ISET1=0, ISET2=NSAVE, ISINC=10                                $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ISET1 =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
  ISET2 =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
NSAVE                                                            
  ISINC =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
10                                                               
   Enter the view latitude, VLAT, and the view longitude, VLON, desired    
   for this plot, in degrees.  VLAT and VLON must be between 10 and 80     
   degrees; each defaults to 45 degrees.  VLAT and VLON are usually        
   defaulted.                                                              
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ press [RETURN] to default VLAT and VLON                              $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  VLAT =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
  VLON =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
                                                                 
   Specify the range (UMIN,UMAX) for the dependent variable axis.  UMIN    
   and UMAX are often defaulted.                                           
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + By default, each plot will be scaled to just fit in the plot area.   +
   + For a common scaling, you may want to set UMIN=ALOW, UMAX=AHIGH.     +
   + ALOW and AHIGH are the minimum and maximum values over all output    +
   + points and over all saved time steps or iterations.                  +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: UMIN = ALOW                                                   $
   $        UMAX = AHIGH                                                  $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  UMIN =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
ALOW                                                             
  UMAX =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
AHIGH                                                            
   Enter a title, WITHOUT quotation marks.  A maximum of 40 characters     
   are allowed.  The default is no title.                                  
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: 3D wave problem                                               $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  TITLE =             (Press [RETURN] to default)
|----Enter title or name---------------|
3D wave problem                         
   To generate tabular or graphical output, choose an output option from   
   the list below.                                                         
                                                                           
   0. No further output is desired                                         
   1. Table of values at output points                                     
         The tabulated output is saved in a file.                          
   2. A plot of the contour (level) surfaces of a scalar variable.         
         This plot gives only a coarse overview of the solution, but is    
         unique in that it shows the general form of the solution in a     
         single picture.  The different levels are identified by color.    
         These plots will only reflect the true geometry if ITRANS=0.      
   3. Surface plot of a scalar variable, at P1, P2 or P3 = constant        
         cross-sections.                                                   
   4. Contour plot of a scalar variable, at P1, P2 or P3 = constant        
         cross-sections.  These plots can be made to reflect the true      
         geometry of the cross-section.                                    
   5. Vector field plot, at P1, P2 or P3 = constant cross-sections.        
         A 3D vector field is plotted using arrows to represent the        
         in-plane components, with a contour plot of the out-of-plane      
         component superimposed.  These plots can be made to reflect       
         the true geometry of the cross-section.                           
   6. One dimensional cross-sectional plots (versus P1, P2, P3 or T)       
                                                                      [RETURN]
                                                                           
   Enter 0,1,2,3,4,5 or 6 to select an output option.                      
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you set NP2=NP3=0 and saved the solution at an arbitrary set of   +
   + user-specified output points, you can only request tabular output.   +
   +                                                                      +
   + If you decide later that you want additional types of plots not      +
   + requested during this interactive session, you will have to work     +
   + through a new interactive session, so it is recommended that you     +
   + request all output or plots you think you MIGHT eventually want now, +
   + during this session.                                                 +
   +                                                                      +
   + Regardless of the options you select, a dummy subroutine POSTPR      +
   + will be included in the program created by the interactive driver;   +
   + you can add your own postprocessing code to this subroutine.         +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want three types of plots for ABS(U): level surface plots,        $
   $ surface plots at a constant Z (P3) cross-section, and a couple of    $
   $ time-plots, so                                                       $
   $ enter: 2, the first time you see this message and                    $
                                                                      [RETURN]
   $        3, the second time                                            $
   $        6, the third time                                             $
   $        0, the fourth time                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter an integer value in the range 0 to 6                  
6                                                                       
   Enter a value for IVAR, to select the variable to be plotted or         
   printed:                                                                
       IVAR = 1 means U  (possibly as modified by UPRINT,..)               
              2       Ux                                                   
              3       Uy                                                   
              4       Uz                                                   
              5       V                                                    
              6       Vx                                                   
              7       Vy                                                   
              8       Vz                                                   
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want to plot ABS(U), which has been saved in U, so                $
   $ enter: IVAR = 1                                                      $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  IVAR =
|---- Enter an integer value in the range 1 to 8                  
1                                                                       
   Which type of cross-sectional plots do you want?                        
                                                                           
   1. Plots of output variable as function of P1 (constant P2,P3 [and T])  
   2. Plots of output variable as function of P2 (constant P1,P3 [and T])  
   3. Plots of output variable as function of P3 (constant P1,P2 [and T])  
      or, if applicable:                                                   
   4. Plots of output variable as function of T (constant P1,P2 and P3)    
                                                                           
   Enter 1,2,3 or 4 to select a plot type.                                 
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: 4                                                             $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter an integer value in the range 1 to 4                  
4                                                                       
   One-dimensional plots of the output variable as a function of T         
   will be made, at the output grid points (P1,P2,P3) closest to           
          (P1CROSS(I),P2CROSS(J),P3CROSS(K)),                              
                I=1,...,NP1VALS, J=1,...,NP2VALS, K=1,...,NP3VALS          
                                                                           
   Enter values for NP1VALS, P1CROSS(1),...,P1CROSS(NP1VALS),              
                and NP2VALS, P2CROSS(1),...,P2CROSS(NP2VALS),              
                and NP3VALS, P3CROSS(1),...,P3CROSS(NP3VALS)               
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want two time plots, at (0,0,0) and (0,0,0.5), so                 $
   $ enter: NP1VALS = 1                                                   $
   $        P1CROSS(1) = 0                                                $
   $        NP2VALS = 1                                                   $
   $        P2CROSS(1) = 0                                                $
   $        NP3VALS = 2                                                   $
   $        P3CROSS(1) = 0                                                $
   $        P3CROSS(2) = 0.5                                              $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  NP1VALS =
|---- Enter an integer value in the range 1 to 100                
1                                                                       
  P1CROSS(1) =         
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
  NP2VALS =
|---- Enter an integer value in the range 1 to 100                
1                                                                       
  P2CROSS(1) =         
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
  NP3VALS =
|---- Enter an integer value in the range 1 to 100                
2                                                                       
  P3CROSS(1) =         
|----Enter constant or FORTRAN expression-----------------------|
0                                                                
  P3CROSS(2) =         
|----Enter constant or FORTRAN expression-----------------------|
0.5                                                              
   Specify the range (UMIN,UMAX) for the dependent variable axis.  UMIN    
   and UMAX are often defaulted.                                           
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + By default, each plot will be scaled to just fit in the plot area.   +
   + For a common scaling, you may want to set UMIN=ALOW, UMAX=AHIGH.     +
   + ALOW and AHIGH are the minimum and maximum values over all output    +
   + points and over all saved time steps or iterations.                  +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: UMIN = ALOW                                                   $
   $        UMAX = AHIGH                                                  $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  UMIN =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
ALOW                                                             
  UMAX =             (Press [RETURN] to default)
|----Enter constant or FORTRAN expression-----------------------|
AHIGH                                                            
   Enter a title, WITHOUT quotation marks.  A maximum of 40 characters     
   are allowed.  The default is no title.                                  
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: 3D wave problem                                               $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  TITLE =             (Press [RETURN] to default)
|----Enter title or name---------------|
3D wave problem                         
   To generate tabular or graphical output, choose an output option from   
   the list below.                                                         
                                                                           
   0. No further output is desired                                         
   1. Table of values at output points                                     
         The tabulated output is saved in a file.                          
   2. A plot of the contour (level) surfaces of a scalar variable.         
         This plot gives only a coarse overview of the solution, but is    
         unique in that it shows the general form of the solution in a     
         single picture.  The different levels are identified by color.    
         These plots will only reflect the true geometry if ITRANS=0.      
   3. Surface plot of a scalar variable, at P1, P2 or P3 = constant        
         cross-sections.                                                   
   4. Contour plot of a scalar variable, at P1, P2 or P3 = constant        
         cross-sections.  These plots can be made to reflect the true      
         geometry of the cross-section.                                    
   5. Vector field plot, at P1, P2 or P3 = constant cross-sections.        
         A 3D vector field is plotted using arrows to represent the        
         in-plane components, with a contour plot of the out-of-plane      
         component superimposed.  These plots can be made to reflect       
         the true geometry of the cross-section.                           
   6. One dimensional cross-sectional plots (versus P1, P2, P3 or T)       
                                                                      [RETURN]
                                                                           
   Enter 0,1,2,3,4,5 or 6 to select an output option.                      
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you set NP2=NP3=0 and saved the solution at an arbitrary set of   +
   + user-specified output points, you can only request tabular output.   +
   +                                                                      +
   + If you decide later that you want additional types of plots not      +
   + requested during this interactive session, you will have to work     +
   + through a new interactive session, so it is recommended that you     +
   + request all output or plots you think you MIGHT eventually want now, +
   + during this session.                                                 +
   +                                                                      +
   + Regardless of the options you select, a dummy subroutine POSTPR      +
   + will be included in the program created by the interactive driver;   +
   + you can add your own postprocessing code to this subroutine.         +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ We want three types of plots for ABS(U): level surface plots,        $
   $ surface plots at a constant Z (P3) cross-section, and a couple of    $
   $ time-plots, so                                                       $
   $ enter: 2, the first time you see this message and                    $
                                                                      [RETURN]
   $        3, the second time                                            $
   $        6, the third time                                             $
   $        0, the fourth time                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter an integer value in the range 0 to 6                  
0                                                                       
   Do you want to define any FORTRAN function subprograms used in any of   
   the FORTRAN 'expressions' entered earlier, entering them line by line?  
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + If you selected double precision accuracy earlier, be sure to        +
   + declare these functions and their arguments DOUBLE PRECISION.        +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: no                                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
no 
   If you don't want to read the FINE PRINT, enter 'no'.                   
                                                                           
   +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++
   + Do you want to define any FORTRAN function subprograms used in any   +
   + of the FORTRAN 'expressions' entered earlier, by interpolating the   +
   + tabular output saved in a file created on an earlier PDE2D run?      +
   ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   $ enter: no                                                            $
   $$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
|---- Enter yes or no
no 
   More detailed information about PDE2D can be found in the book          
   "Solving Partial Differential Equation Applications with PDE2D,"        
   Granville Sewell, John Wiley & Sons, 2018.                              
                                                                           
              *******************************************                  
              ***** Input program has been created  *****                  
              *******************************************