C     **************************                                              
C     * PDE2D 6.7 MAIN PROGRAM *                                              
C     **************************                                              
C     *** 1D PROBLEM SOLVED ***                                               
C##############################################################################
C     Is double precision mode to be used?  Double precision is recommended   #
C     on 32-bit computers.                                                    #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + If double precision mode is used, variables and functions assigned   +#
C     + names beginning with a letter in the range A-H or O-Z will be DOUBLE +#
C     + PRECISION, and you should use double precision constants and FORTRAN +#
C     + expressions throughout; otherwise such variables and functions will  +#
C     + be of type REAL.  In either case, variables and functions assigned   +#
C     + names beginning with I,J,K,L,M or N will be of INTEGER type.         +#
C     +                                                                      +#
C     + It is possible to convert a single precision PDE2D program to double +#
C     + precision after it has been created, using an editor.  Just change   +#
C     + all occurrences of "REAL" to "DOUBLE PRECISION"                      +#
C     +                    " TDP" to "DTDP"  (note leading blanks)           +#
C     +                    " TD3" to "DTD3"                                  +#
C     + Any user-written code or routines must be converted "by hand", of    +#
C     + course.  To convert from double to single, reverse the changes.      +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NEQNMX=  99)
C##############################################################################
C     NXGRID = number of X-grid lines                                         #
C##############################################################################
      PARAMETER (NXGRID = 21)        
C##############################################################################
C     How many partial differential equations (NEQN) are there in your        #
C     problem?                                                                #
C##############################################################################
      PARAMETER (NEQN = 4)           
C        DIMENSIONS OF WORK ARRAYS RWRK8Z AND IWRK8Z                          
      PARAMETER (IRWK8Z=       11928)
      PARAMETER (IIWK8Z=       16968)
      PARAMETER (NXP8Z=1001,KDEG8Z=1)
C##############################################################################
C     The solution is saved on a uniform grid of NX+1 points                  #
C                        XA + I*(XB-XA)/NX                                    #
C     I=0,...,NX.  Enter a value for NX (suggested value = 50).               #
C                                                                             #
C##############################################################################
      PARAMETER (NX = 50)            
      PARAMETER (NSAVE = 1)
      COMMON/PARM8Z/ PAR18Z
      DIMENSION XGRID(NXGRID),XOUT8Z(0:NX),XCROSS(100),TOUT8Z(0:NSAVE),U       
     &OUT8Z(0:NX,2*NEQN,0:NSAVE)                                               
      DIMENSION IWRK8Z(IIWK8Z),RWRK8Z(IRWK8Z)                                  
      CHARACTER*40 TITLE                                                       
      LOGICAL LINEAR,CRANKN,NOUPDT,NODIST,FILLIN,EIGEN,EVCMPX,ADAPT,PLOT       
     &,LSQFIT,FDIFF,SOLID,ECON8Z,NCON8Z                                        
      COMMON/DTDP14/ SINT8Z(20),BINT8Z(20),SLIM8Z(20),BLIM8Z(20)               
      COMMON/DTDP15/ EVLR8Z,EV0R,EVLI8Z,EV0I,EVCMPX                            
      COMMON/DTDP16/ NINT,KINT8Z,NBINT                                         
      COMMON/DTDP19/ ADAPT,TOLER(NEQNMX)                                       
      COMMON/DTDP26/ IPERDC,PLOT,LSQFIT,FDIFF,SOLID                            
      COMMON/DTDP28/ EPS8Z                                                     
      COMMON/DTDP30/ ECON8Z,NCON8Z                                             
      COMMON/DTDP42/ NXA8Z,KD8Z                                                
      COMMON/DTDP43/ WORK8Z(NXP8Z+3)                                           
      COMMON/DTDP45/ PERDC(NEQNMX)                                             
      NXA8Z = NXP8Z                                                            
      KD8Z = KDEG8Z                                                            
C##############################################################################
C     If you don't want to read the FINE PRINT, default NPROB.                #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + If you want to solve several similar problems in the same run, set   +#
C     + NPROB equal to the number of problems you want to solve.  Then NPROB +#
C     + loops through the main program will be done, with IPROB=1,...,NPROB, +#
C     + and you can make the problem parameters vary with IPROB.  NPROB      +#
C     + defaults to 1.                                                       +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
      NPROB = 1                                                                
      DO 78755 IPROB=1,NPROB                                                   
C##############################################################################
C     Definitions for global parameters (passed through COMMON /PARM8Z/)      #
C##############################################################################
C##############################################################################
C     A collocation finite element method is used, with cubic Hermite         #
C     basis functions on the subintervals defined by the grid points:         #
C              XGRID(1),XGRID(2),...,XGRID(NXGRID)                            #
C     You will first be prompted for NXGRID, the number of X-grid points,     #
C     then for XGRID(1),...,XGRID(NXGRID).  Any points defaulted will be      #
C     uniformly spaced between the points you define; the first and last      #
C     points cannot be defaulted.  The interval over which the PDE system     #
C     is to be solved is then:                                                #
C               XGRID(1) < X < XGRID(NXGRID)                                  #
C                                                                             #
C##############################################################################
      CALL DTD3W(XGRID,NXGRID,0)                                               
C        XGRID DEFINED                                                         
      XGRID(1) =                                                               
     & 0                                                                       
      XGRID(NXGRID) =                                                          
     & 1                                                                       
      CALL DTD3W(XGRID,NXGRID,1)                                               
C        *******STEADY-STATE PROBLEM                                           
      EIGEN = .FALSE.                                                          
      T0 = 0.0                                                                 
      DT = 1.0                                                                 
      CRANKN = .FALSE.                                                         
      NOUPDT = .FALSE.                                                         
C##############################################################################
C     Are all partial differential equations and all boundary conditions      #
C     linear?                                                                 #
C##############################################################################
      LINEAR = .FALSE.                                                         
C##############################################################################
C     Give an upper limit on the number of Newton's method iterations         #
C     (NSTEPS) to be allowed for this nonlinear problem.  A suggested value   #
C     is NSTEPS = 20.                                                         #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + The iteration will stop if convergence occurs before the upper       +#
C     + limit has been reached.  In the FORTRAN program created by the       +#
C     + preprocessor, NCON8Z will be .TRUE. if Newton's method has converged.+#
C     +                                                                      +#
C     + For highly non-linear problems you may want to construct a one-      +#
C     + parameter family of problems using the variable T, such that for     +#
C     + T=1 the problems is easy (e.g. linear) and for T > N (N is less      +#
C     + NSTEPS), the problem reduces to the original highly nonlinear        +#
C     + problem.  For example, the nonlinear term(s) may be multiplied by    +#
C     + MIN(1.0,(T-1.0)/N).                                                  +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
      NSTEPS =         20                                                      
C##############################################################################
C     If you don't want to read the FINE PRINT, enter 'yes'.                  #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + The partial derivatives of some of the PDE and boundary condition    +#
C     + coefficients are required by PDE2D.  These may be calculated         +#
C     + automatically using a finite difference approximation, or supplied   +#
C     + by the user.  Do you want them to be calculated automatically?       +#
C     +                                                                      +#
C     + If you answer 'yes', you will not be asked to supply the derivatives,+#
C     + but there is a small risk that the inaccuracies introduced by the    +#
C     + finite difference approximation may cause the Newton iteration       +#
C     + to converge more slowly or to diverge, especially if low precision   +#
C     + is used.  This risk is very low, however, and since answering 'no'   +#
C     + means you may have to compute many partial derivatives, it is        +#
C     + recommended you answer 'yes' unless you have some reason to believe  +#
C     + there is a problem with the finite difference approximations.        +#
C     +                                                                      +#
C     + If you supply analytic partial derivatives, PDE2D will do some spot  +#
C     + checking and can usually issue a warning if any are supplied         +#
C     + incorrectly.                                                         +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
      FDIFF = .TRUE.                                                           
C##############################################################################
C     You may calculate one or more integrals (over the entire region) of     #
C     some functions of the solution and its derivatives.  How many integrals #
C     (NINT), if any, do you want to calculate?                               #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + In the FORTRAN program created by the preprocessor, the computed     +#
C     + values of the integrals will be returned in the vector SINT8Z.  If   +#
C     + several iterations or time steps are done, only the last computed    +#
C     + values are saved in SINT8Z (all values are printed).                 +#
C     +                                                                      +#
C     + A limiting value, SLIM8Z(I), for the I-th integral can be set at the +#
C     + beginning of the main program.  The computations will then stop      +#
C     + gracefully whenever SINT8Z(I) > SLIM8Z(I), for any I=1...NINT.       +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
      NINT =          1                                                        
C##############################################################################
C     You may calculate one or more boundary integrals (over the entire       #
C     boundary) of some functions of the solution and its derivatives.  How   #
C     many boundary integrals (NBINT), if any, do you want to calculate?      #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + In the FORTRAN program created by the preprocessor, the computed     +#
C     + values of the integrals will be returned in the vector BINT8Z.  If   +#
C     + several iterations or time steps are done, only the last computed    +#
C     + values are saved in BINT8Z (all values are printed).                 +#
C     +                                                                      +#
C     + A limiting value, BLIM8Z(I), for the I-th boundary integral can be   +#
C     + set at the beginning of the main program.  The computations will then+#
C     + stop gracefully whenever BINT8Z(I) > BLIM8Z(I), for any I=1...NBINT. +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
      NBINT =          0                                                       
      LSQFIT = .FALSE.                                                         
C##############################################################################
C     If you do not have periodic boundary conditions, enter IPERDC=0.        #
C                                                                             #
C     Enter IPERDC=1 for periodic conditions at X = XGRID(1),XGRID(NXGRID)    #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + When periodic boundary conditions are selected, they apply to all    +#
C     + variables by default.  To prevent them from applying to the I-th     +#
C     + variable, set PERDC(I)=0 at the beginning of the main program and    +#
C     + set the desired boundary conditions in subroutine GB8Z, "by hand".   +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
      IPERDC =          0                                                      
C##############################################################################
C     The solution is saved on a uniform grid of NX+1 points, covering the    #
C     interval (XA,XB).  Enter values for XA,XB.  These variables are usually #
C     defaulted.                                                              #
C                                                                             #
C     The defaults are XA = XGRID(1), XB = XGRID(NXGRID)                      #
C                                                                             #
C##############################################################################
C        DEFAULTS FOR XA,XB                                                    
      XA = XGRID(1)                                                            
      XB = XGRID(NXGRID)                                                       
C        DEFINE XA,XB HERE:                                                    
      NPTS8Z = NX+1                                                            
      HX8Z = (XB-XA)/NX                                                        
      DO 78756 I8Z=0,NX                                                        
         XOUT8Z(I8Z) = XA + I8Z*HX8Z                                           
78756 CONTINUE                                                                 
C        SOLUTION SAVED EVERY NOUT ITERATIONS                                  
      NOUT = NSTEPS                                                            
      II8Z = IIWK8Z                                                            
      IR8Z = IRWK8Z                                                            
C        *******DRAW GRID POINTS?                                              
      PLOT = .TRUE.                                                            
C        *******CALL PDE SOLVER                                                
      CALL DTDP1 (XGRID, NXGRID, NEQN, XOUT8Z, UOUT8Z, TOUT8Z, NPTS8Z, T       
     &0, DT, NSTEPS, NOUT, NSAVE+1, CRANKN, NOUPDT, EIGEN, LINEAR, RWRK8       
     &Z, IR8Z, IWRK8Z, II8Z)                                                   
C        *******CALL USER-WRITTEN POSTPROCESSOR                                
      CALL POSTPR(TOUT8Z,NSAVE,XOUT8Z,NX,UOUT8Z,NEQN)                          
C        *******LINE PLOTS                                                     
C##############################################################################
C     Enter a value for IVAR, to select the variable to be plotted or         #
C     printed:                                                                #
C         IVAR = 1 means Q  (possibly as modified by UPRINT,..)               #
C                2       Qx                                                   #
C                3       DQ                                                   #
C                4       DQx                                                  #
C                5       P                                                    #
C                6       Px                                                   #
C                7       DP                                                   #
C                8       DPx                                                  #
C                .        .                                                   #
C                .        .                                                   #
C##############################################################################
      IVAR =          6                                                        
C        X IS VARIABLE                                                         
      ICS8Z = 1                                                                
      ISET1 = 1                                                                
      ISET2 = NSAVE                                                            
      ISINC = 1                                                                
C##############################################################################
C     Specify the range (UMIN,UMAX) for the dependent variable axis.  UMIN    #
C     and UMAX are often defaulted.                                           #
C                                                                             #
C     The default is UMIN=Ulow, UMAX=Uhigh, where Ulow and Uhigh are the      #
C     minimum and maximum values of the variable to be plotted.  (Note:       #
C     you must specify both UMIN and UMAX, or default both.)                  #
C##############################################################################
      UMIN = 0.0                                                               
      UMAX = 0.0                                                               
C##############################################################################
C     Enter a title, WITHOUT quotation marks.  A maximum of 40 characters     #
C     are allowed.  The default is no title.                                  #
C##############################################################################
      TITLE = ' '                                                              
      TITLE = 'Q                                       '                       
      CALL DTD3R(TOUT8Z,NSAVE,ISET1,ISET2,ISINC)                               
      DO 78757 IS8Z=ISET1,ISET2,ISINC                                          
      CALL DTDPZP(ICS8Z,IVAR,TOUT8Z,NSAVE,XOUT8Z,NX,UOUT8Z,NEQN,TITLE,UM       
     &IN,UMAX,IX8Z,IS8Z)                                                       
78757 CONTINUE                                                                 
C        *******LINE PLOTS                                                     
C##############################################################################
C     Enter a value for IVAR, to select the variable to be plotted or         #
C     printed:                                                                #
C         IVAR = 1 means Q  (possibly as modified by UPRINT,..)               #
C                2       Qx                                                   #
C                3       DQ                                                   #
C                4       DQx                                                  #
C                5       P                                                    #
C                6       Px                                                   #
C                7       DP                                                   #
C                8       DPx                                                  #
C                .        .                                                   #
C                .        .                                                   #
C##############################################################################
      IVAR =          5                                                        
C        X IS VARIABLE                                                         
      ICS8Z = 1                                                                
      ISET1 = 1                                                                
      ISET2 = NSAVE                                                            
      ISINC = 1                                                                
C##############################################################################
C     Specify the range (UMIN,UMAX) for the dependent variable axis.  UMIN    #
C     and UMAX are often defaulted.                                           #
C                                                                             #
C     The default is UMIN=Ulow, UMAX=Uhigh, where Ulow and Uhigh are the      #
C     minimum and maximum values of the variable to be plotted.  (Note:       #
C     you must specify both UMIN and UMAX, or default both.)                  #
C##############################################################################
      UMIN = 0.0                                                               
      UMAX = 0.0                                                               
C##############################################################################
C     Enter a title, WITHOUT quotation marks.  A maximum of 40 characters     #
C     are allowed.  The default is no title.                                  #
C##############################################################################
      TITLE = ' '                                                              
      TITLE = 'P                                       '                       
      CALL DTD3R(TOUT8Z,NSAVE,ISET1,ISET2,ISINC)                               
      DO 78758 IS8Z=ISET1,ISET2,ISINC                                          
      CALL DTDPZP(ICS8Z,IVAR,TOUT8Z,NSAVE,XOUT8Z,NX,UOUT8Z,NEQN,TITLE,UM       
     &IN,UMAX,IX8Z,IS8Z)                                                       
78758 CONTINUE                                                                 
78755 CONTINUE                                                                 
      CALL ENDGKS                                                              
      STOP                                                                     
      END                                                                      
                                                                               
                                                                               
      SUBROUTINE PDES8Z(YD8Z,I8Z,X,T,UU8Z)                                     
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      PARAMETER (NEQNMX=  99)                                                  
C        UN8Z(1,I),UN8Z(2,I),... HOLD THE (RARELY USED) VALUES                 
C        OF UI,UIx,... FROM THE PREVIOUS ITERATION OR TIME STEP                
      COMMON /DTD32/UN8Z(3,NEQNMX)                                             
      COMMON /DTDP11/NORMX                                                     
      COMMON /DTDP16/NINT8Z,KINT8Z,NBIN8Z                                      
      COMMON /DTDP44/PDER8Z,J8Z                                                
      DOUBLE PRECISION NORMX                                                   
      DIMENSION UU8Z(3,NEQNMX)                                                 
      COMMON/PARM8Z/ PAR18Z                                                    
      Q  = UU8Z(1, 1)                                                          
      QX = UU8Z(2, 1)                                                          
      QXX= UU8Z(3, 1)                                                          
      DQ  = UU8Z(1, 2)                                                         
      DQX = UU8Z(2, 2)                                                         
      DQXX= UU8Z(3, 2)                                                         
      P  = UU8Z(1, 3)                                                          
      PX = UU8Z(2, 3)                                                          
      PXX= UU8Z(3, 3)                                                          
      DP  = UU8Z(1, 4)                                                         
      DPX = UU8Z(2, 4)                                                         
      DPXX= UU8Z(3, 4)                                                         
                          IF (I8Z.EQ.0) THEN                                   
      YD8Z = 0.0                                                               
C##############################################################################
C     Enter FORTRAN expressions for the functions whose integrals are to be   #
C     calculated and printed.  They may be functions of                       #
C                                                                             #
C       X,Q,Qx,Qxx,DQ,DQx,DQxx,P,Px,Pxx,DP,DPx,DPxx ...                       #
C                                                and (if applicable) T        #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + If you only want to integrate a function over part of the interval,  +#
C     + define that function to be zero on the rest of the interval.         +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
C                                                  INTEGRAL DEFINED            
      IF (KINT8Z.EQ.  1) YD8Z =                                                
     & Q                                                                       
C                                                  BND. INTEGRAL1 DEFINED      
C     IF (KINT8Z.EQ.-1) YD8Z =                                                 
C    & [DEFAULT SELECTED, DEFINITION COMMENTED OUT]                            
                          ELSE                                                 
C##############################################################################
C     Now enter FORTRAN expressions to define the PDE coefficients, which     #
C     may be functions of                                                     #
C                                                                             #
C        X,Q,Qx,Qxx,DQ,DQx,DQxx,P,Px,Pxx,DP,DPx,DPxx,...                      #
C                                                                             #
C     and, in some cases, of the parameter T.                                 #
C                                                                             #
C     Recall that the PDEs have the form                                      #
C                                                                             #
C                        F1 = 0                                               #
C                        F2 = 0                                               #
C                        F3 = 0                                               #
C                        F4 = 0                                               #
C                         .   .                                               #
C                                                                             #
C##############################################################################
                IF (J8Z.EQ.0) THEN                                             
      YD8Z = 0.0                                                               
C                                                  F1 DEFINED                  
      IF (I8Z.EQ.  1) YD8Z =                                                   
     & Qx-DQ                                                                   
C                                                  F2 DEFINED                  
      IF (I8Z.EQ.  2) YD8Z =                                                   
     & DQx-1.0                                                                
C                                                  F3 DEFINED                  
      IF (I8Z.EQ.  3) YD8Z =                                                   
     & Px-DP                                                                   
C                                                  F4 DEFINED                  
      IF (I8Z.EQ.  4) YD8Z =                                                   
     & .06283*DPx-(Qx**2)*P*(.01**2+P*P-P*(.01**2+P*P)**.5)/
     & (.01**2+P*P)**.5       
                ELSE                                                           
                ENDIF                                                          
                          ENDIF                                                
      RETURN                                                                   
      END                                                                      
                                                                               
                                                                               
      FUNCTION U8Z(I8Z,X)                                                      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      LOGICAL RESTRT,GRIDID                                                    
      COMMON/PARM8Z/ PAR18Z                                                    
      U8Z = 0.0                                                                
C##############################################################################
C     Now the initial values for Newton's method must be defined using        #
C     FORTRAN expressions.  They may be functions of X.                       #
C                                                                             #
C     It is important to provide initial values which are at least of the     #
C     correct order of magnitude.                                             #
C##############################################################################
C                                                  Q0 DEFINED                  
      IF (I8Z.EQ.  1) U8Z =                                                    
     & 4.7                                                                   
C                                                  DQ0 DEFINED                 
      IF (I8Z.EQ.  2) U8Z =                                                    
     & 0                                                                       
C                                                  P0 DEFINED                  
      IF (I8Z.EQ.  3) U8Z =                                                    
     & 1                                                                    
C                                                  DP0 DEFINED                 
      IF (I8Z.EQ.  4) U8Z =                                                    
     & 0                                                                       
C##############################################################################
C     If you don't want to read the FINE PRINT, enter 'no'.                   #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + Do you want to read the initial conditions from the restart file,    +#
C     + if it exists (and use the conditions supplied above if it does not   +#
C     + exist)?                                                              +#
C     +                                                                      +#
C     + If so, PDE2D will dump the final solution at the end of each run     +#
C     + into a restart file "pde2d.res".  Thus the usual procedure for       +#
C     + using this dump/restart option is to make sure there is no restart   +#
C     + file in your directory left over from a previous job, then the       +#
C     + first time you run this job, the initial conditions supplied above   +#
C     + will be used, but on the second and subsequent runs the restart file +#
C     + from the previous run will be used to define the initial conditions. +#
C     +                                                                      +#
C     + You can do all the "runs" in one program, by setting NPROB > 1.      +#
C     + Each pass through the DO loop, T0,TF,NSTEPS and possibly other       +#
C     + parameters may be varied, by making them functions of IPROB.         +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
      RESTRT = .FALSE.                                                         
C     GRIDID = .FALSE. IF FINITE ELEMENT GRID CHANGES BETWEEN DUMP, RESTART    
      GRIDID = .TRUE.                                                          
      IF (I8Z.EQ.0 .AND. RESTRT .AND. GRIDID) U8Z = 1.0                        
      IF (I8Z.EQ.0 .AND. RESTRT .AND. (.NOT.GRIDID)) U8Z = -1.0                
      RETURN                                                                   
      END                                                                      
                                                                               
                                                                               
      SUBROUTINE GB8Z(GD8Z,IFAC8Z,I8Z,X,T,UU8Z)                                
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      PARAMETER (NEQNMX=  99)                                                  
      DIMENSION UU8Z(3,NEQNMX)                                                 
C        UN8Z(1,I),UN8Z(2,I),... HOLD THE (RARELY USED) VALUES                 
C        OF UI,UIx,... FROM THE PREVIOUS ITERATION OR TIME STEP                
      COMMON /DTD32/ UN8Z(3,NEQNMX)                                            
      COMMON /DTDP44/PDER8Z,J8Z                                                
      DOUBLE PRECISION NONE                                                    
      COMMON/PARM8Z/ PAR18Z                                                    
      NONE = DTD3L(2)                                                          
      Q  = UU8Z(1, 1)                                                          
      QX = UU8Z(2, 1)                                                          
      DQ  = UU8Z(1, 2)                                                         
      DQX = UU8Z(2, 2)                                                         
      P  = UU8Z(1, 3)                                                          
      PX = UU8Z(2, 3)                                                          
      DP  = UU8Z(1, 4)                                                         
      DPX = UU8Z(2, 4)                                                         
      IF (J8Z.EQ.0) GD8Z = 0.0                                                 
C##############################################################################
C     Enter FORTRAN expressions to define the boundary condition functions,   #
C     which may be functions of                                               #
C                                                                             #
C            X,Q,Qx,DQ,DQx,P,Px,DP,DPx,... and (if applicable) T              #
C                                                                             #
C     Recall that the boundary conditions have the form                       #
C                                                                             #
C                           G1 = 0                                            #
C                           G2 = 0                                            #
C                           G3 = 0                                            #
C                           G4 = 0                                            #
C                            .   .                                            #
C     Enter NONE to indicate "no" boundary condition.                         #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + If "no" boundary condition is specified, the corresponding PDE is    +#
C     + enforced at a point just inside the boundary (exactly on the         +#
C     + boundary, if EPS8Z is set to 0 in the main program).                 +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
            IF (IFAC8Z.EQ. 1) THEN                                             
C##############################################################################
C                                                                             #
C     First define the boundary conditions at the point X = XGRID(1).         #
C##############################################################################
                          IF (J8Z.EQ.0) THEN                                   
C                                                  G1 DEFINED                  
      IF (I8Z.EQ.  1) GD8Z =                                                   
     & Q-4.7                                                                 
C                                                  G2 DEFINED                  
      IF (I8Z.EQ.  2) GD8Z =                                                   
     & DQ                                                                      
C                                                  G3 DEFINED                  
      IF (I8Z.EQ.  3) GD8Z =                                                   
     & P-.1                                                                 
C                                                  G4 DEFINED                  
      IF (I8Z.EQ.  4) GD8Z =                                                   
     & none                                                                   
                          ELSE                                                 
                          ENDIF                                                
            ENDIF                                                              
            IF (IFAC8Z.EQ. 2) THEN                                             
C##############################################################################
C                                                                             #
C     Now define the boundary conditions at the point X = XGRID(NXGRID).      #
C##############################################################################
                          IF (J8Z.EQ.0) THEN                                   
C                                                  G1 DEFINED                  
      IF (I8Z.EQ.  1) GD8Z =                                                   
     & none                                                                    
C                                                  G2 DEFINED                  
      IF (I8Z.EQ.  2) GD8Z =                                                   
     & none                                                                    
C                                                  G3 DEFINED                  
      IF (I8Z.EQ.  3) GD8Z =                                                   
     & none                                                                    
C                                                  G4 DEFINED                  
      IF (I8Z.EQ.  4) GD8Z =                                                   
     & P-1.0                                                                    
                          ELSE                                                 
                          ENDIF                                                
            ENDIF                                                              
      RETURN                                                                   
      END                                                                      
                                                                               
                                                                               
      SUBROUTINE PMOD8Z(X,T,UU8Z,UPRINT,UXPRNT)                                
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      DIMENSION UU8Z(3,*),UPRINT(*),UXPRNT(*)                                  
      COMMON/DTDP14/SINT(20),BINT(20),SLIM8Z(20),BLIM8Z(20)                    
      COMMON/PARM8Z/ PAR18Z                                                    
      Q  = UU8Z(1, 1)                                                          
      QX = UU8Z(2, 1)                                                          
      QXX= UU8Z(3, 1)                                                          
      DQ  = UU8Z(1, 2)                                                         
      DQX = UU8Z(2, 2)                                                         
      DQXX= UU8Z(3, 2)                                                         
      P  = UU8Z(1, 3)                                                          
      PX = UU8Z(2, 3)                                                          
      PXX= UU8Z(3, 3)                                                          
      DP  = UU8Z(1, 4)                                                         
      DPX = UU8Z(2, 4)                                                         
      DPXX= UU8Z(3, 4)                                                         
C##############################################################################
C     If you don't want to read the FINE PRINT, default all of the following  #
C     variables.                                                              #
C                                                                             #
C     +++++++++++++++ THE "FINE PRINT" (CAN USUALLY BE IGNORED) ++++++++++++++#
C     + Normally, PDE2D saves the values of Q,Qx,DQ,DQx,P,Px,DP,             +#
C     + DPx...at the output points.  If different variables are to be saved  +#
C     + (for later printing or plotting) the following functions can be used +#
C     + to re-define the output variables:                                   +#
C     +    define UPRINT(1) to replace  Q                                    +#
C     +           UPRINT(2)             DQ                                   +#
C     +           UPRINT(3)             P                                    +#
C     +           UPRINT(4)             DP                                   +#
C     +                  .               .                                   +#
C     +           UXPRNT(1)             Qx                                   +#
C     +           UXPRNT(2)             DQx                                  +#
C     +           UXPRNT(3)             Px                                   +#
C     +           UXPRNT(4)             DPx                                  +#
C     +                   .              .                                   +#
C     + Each function may be a function of                                   +#
C     +                                                                      +#
C     +    X,Q,Qx,Qxx,DQ,DQx,DQxx,P,Px,Pxx,DP,DPx,DPxx,...                   +#
C     +                  and (if applicable) T                               +#
C     +                                                                      +#
C     + Each may also be a function of the integral estimates SINT(1),...,   +#
C     + BINT(1),...                                                          +#
C     +                                                                      +#
C     + The default for each variable is no change, for example, UPRINT(1)   +#
C     + defaults to Q.  Enter FORTRAN expressions for each of the            +#
C     + following functions (or default).                                    +#
C     ++++++++++++++++++++++++++ END OF "FINE PRINT" +++++++++++++++++++++++++#
C##############################################################################
C        DEFINE UPRINT(*),UXPRNT(*) HERE:                                      
      RETURN                                                                   
      END                                                                      
C        DUMMY ROUTINES                                                        
      FUNCTION XY8Z(I8Z,IARC8Z,S)                                              
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      XY8Z = 0                                                                 
      RETURN                                                                   
      END                                                                      
      FUNCTION DIS8Z(I8Z,X,Y,KTRI)                                             
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      DIS8Z = 0                                                                
      RETURN                                                                   
      END                                                                      
      FUNCTION FB8Z(I8Z,IARC8Z,S,X,Y,T)                                        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      FB8Z = 0                                                                 
      RETURN                                                                   
      END                                                                      
                                                                               
                                                                               
      SUBROUTINE POSTPR(TOUT,NSAVE,XOUT,NX,UOUT,NEQN)                          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                                      
      DIMENSION XOUT(0:NX),TOUT(0:NSAVE)                                       
      DIMENSION UOUT(0:NX,2,NEQN,0:NSAVE)                                      
C     UOUT(I,1,IEQ,L) = U-sub-IEQ                                              
C     UOUT(I,2,IEQ,L) = UX-sub-IEQ                                             
C        (possibly as modified by UPRINT,..)                                   
C        at the point XOUT(I)                                                  
C        at the time TOUT(L).                                                  
C        ******* ADD POSTPROCESSING CODE HERE:                                 
      RETURN                                                                   
      END                                                                      

