C
C       ******************************************************
C      *                                                      *
C      *                  SUBROUTINE EXERFC                   *
C      *                                                      *
C      *            VERSION CURRENT AS OF 10/01/87            *
C      *                                                      *
C       ******************************************************
C
      SUBROUTINE EXERFC (X,YY,Z)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION P1(5),Q1(5),P2(9),Q2(9),P3(6),Q3(6)
C
C       THIS ROUTINE USES RATIONAL CHEBYSHEV APPROXIMATIONS
C       FOR EVALUATING THE ERROR FUNCTION AND COMPLEMENTARY
C       ERROR FUNCTION IN ORDER TO EVALUATE THE PRODUCT OF
C       EXP(X) AND ERFC(Y)
C
      DATA P1/3.209377589138469472562D03,3.774852376853020208137D02,
     1        1.138641541510501556495D02,3.161123743870565596947D0,
     2        1.857777061846031526730D-01/
      DATA Q1/2.844236833439170622273D03,1.282616526077372275645D03,
     1        2.440246379344441733056D02,2.360129095234412093499D01,
     2        1.0D0 /
      DATA P2/1.23033935479799725272D03,2.05107837782607146532D03,
     1        1.71204761263407058314D03,8.81952221241769090411D02,
     2        2.98635138197400131132D02,6.61191906371416294775D01,
     3       8.88314979438837594118D00,5.64188496988670089180D-01,
     4       2.15311535474403846343D-08/
      DATA Q2/1.23033935480374942043D03,3.43936767414372163696D03,
     1        4.36261909014324715820D03,3.29079923573345962678D03,
     2        1.62138957456669018874D03,5.37181101862009857509D02,
     3        1.17693950891312499305D02,1.57449261107098347253D01,
     4        1.0D0 /
      DATA P3/-6.58749161529837803157D-04,-1.60837851487422766278D-02,
     1        -1.25781726111229246204D-01,-3.60344899949804439429D-01,
     2        -3.05326634961232344035D-01,-1.63153871373020978498D-02/
      DATA Q3/2.33520497626869185443D-03,6.05183413124413191178D-02,
     1        5.27905102951428412248D-01,1.87295284992346047209D00,
     2        2.56852019228982242072D00,1.0D0/
C
      IF(YY.EQ.0.0D0) Z=DEXP(X)
      IF(YY.EQ.0.0D0) RETURN
      Y=DABS(YY)
C
C        FOR 0.0 < Y < .46875
      IF (Y.GT.0.46875D0) GO TO 20
      SUMP=0.0D0
      SUMQ=0.0D0
      DO 10 I=1,5
      Y2I=Y**(2*(I-1))
      SUMP=SUMP+P1(I)*Y2I
      SUMQ=SUMQ+Q1(I)*Y2I
10    CONTINUE
      ERF=Y*SUMP/SUMQ
      IF(YY.LT.0.0) ERF=-ERF
      ERFCY=1.0D0-ERF
      Z=DEXP(X)*ERFCY
      RETURN
C
C        FOR 0.0 < Y < .46875
20    IF (Y.GT.4.0D0) GO TO 40
      SUMP=0.0D0
      SUMQ=0.0D0
      DO 30 I=1,9
      YI=Y**(I-1)
      SUMP=SUMP+P2(I)*YI
      SUMQ=SUMQ+Q2(I)*YI
30    CONTINUE
      Z=DEXP(X-Y*Y)*SUMP/SUMQ
      IF(YY.LT.0.0D0) Z=2.0D0*DEXP(X)-Z
      RETURN
40    SUMP=0.0D0
      SUMQ=0.0D0
      DO 50 I=1,6
      Y2I=Y**(-2*(I-1))
      SUMP=SUMP+P3(I)*Y2I
      SUMQ=SUMQ+Q3(I)*Y2I
50    CONTINUE
      SQRTPI=0.5641895835477562869481D0
      Z=SQRTPI+SUMP/(Y*Y*SUMQ)
      Z=DEXP(X-Y*Y)*Z/Y
      IF(YY.LT.0.0D0) Z=2.0D0*DEXP(X)-Z
      RETURN
      END
C
C       ******************************************************
C      *                                                      *
C      *                  SUBROUTINE GLQPTS                   *
C      *                                                      *
C      *            VERSION CURRENT AS OF 10/01/87            *
C      *                                                      *
C       ******************************************************
C
      SUBROUTINE GLQPTS (N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CHARACTER*1 SKIP
      COMMON /GLPTS/ WN(256),ZN(256)
      COMMON /IOUNIT/ IN,IO
C
C       THIS ROUTINE READS THE NORMALIZED ROOTS ZN(I) AND WEIGHTS WN(I)
C       OF THE LEGENDRE POLYNOMIALS FROM THE DATA FILE 'GLQ.PTS'
C
C       N IS THE NUMBER OF INTEGRATION POINTS AND CAN ONLY HAVE A
C       VALUE OF EITHER 4,20,60,104,OR 256
C
      IN2=77
      OPEN(IN2,FILE='GLQ.PTS',STATUS='OLD')
C
C       SKIP LINES IN FILE UNTIL CORRECT COEFFICIENTS ARE REACHED
      ISKIP=-1
      IF(N.EQ.4) ISKIP=7
      IF(N.EQ.20) ISKIP=9
      IF(N.EQ.60) ISKIP=15
      IF(N.EQ.104) ISKIP=31
      IF(N.EQ.256) ISKIP=57
      IF (ISKIP.EQ.-1) WRITE(IO,201)
      IF (ISKIP.EQ.-1) STOP
      DO 60 I=1,ISKIP
60    READ(IN2,101) SKIP
C
C       READ IN ZN(I) AND WN(I), FOUR VALUES PER LINE
      NC=N/8
      IF (MOD(N,8).NE.0) NC=NC+1
      DO 80 I=1,NC
      K=(I-1)*8-1
      READ(IN2,102) (ZN(K+J*2),J=1,4)
80    CONTINUE
      DO 100 I=1,NC
      K=(I-1)*8-1
      READ(IN2,102) (WN(K+J*2),J=1,4)
100   CONTINUE
C
C       FILL IN THE SYMMETRIC TERMS
      DO 120 J=2,N,2
      J1=J-1
      ZN(J)=-ZN(J1)
120   WN(J)=WN(J1)
      CLOSE(IN2)
      RETURN
C
C       FORMAT STATEMENTS
101   FORMAT(A1)
102   FORMAT(4D20.0)
201   FORMAT(1H0,20X,'*****  ERROR IN ROUTINE GLQPTS  *****'/
     1 1H ,20X,'NO. OF ROOTS SPECIFIED MUST EQUAL 4,20,60,104 OR 256')
      END
C
C       ******************************************************
C      *                                                      *
C      *                  SUBROUTINE TITLE4                   *
C      *                                                      *
C      *            VERSION CURRENT AS OF 10/01/87            *
C      *                                                      *
C       ******************************************************
C
      SUBROUTINE TITLE4
      CHARACTER*1 LINE1(60),EQUAL,BLANK
      CHARACTER DATE*16,TIME*8
      CHARACTER*60 LINE
      CHARACTER*61 TITLE,T1
      COMMON /IOUNIT/ IN,IO
      COMMON /TITLES/ TITLE(4)
      DATA EQUAL/'='/,BLANK/' '/
C
C        THIS ROUTINE CREATES A TITLE BOX ON THE FIRST PAGE OF
C        PROGRAM OUTPUT. THE ROUTINE READS AND PRINTS ALL DATA
C        CARDS UNTIL IT ENCOUNTERS AN '=' IN COLUMN 1. THE FIRST 4
C        LINES READ IN ARE ALSO USED AS TITLES ON PLOTS.
C
CCC   CALL TIME$A (TIME)
CCC   CALL DATE$A(DATE)
      WRITE(IO,201)
      DO 10 L=1,60
      READ(IN,101,END=20) LINE
      IF (LINE(1:1).EQ.EQUAL) GOTO 60
      T1=LINE
C       STRIP OFF TRAILING BLANKS AND CENTER LINE
      DO 15 N=1,60
      NN=61-N
15    IF(LINE(NN:NN).NE.BLANK) GOTO 20
20    NN1=NN+1
      T1(NN1:NN1)='$'
      IF(L.LT.5) TITLE(L)=T1
      NS=(60-NN)/2
      IF(NS.EQ.0) GO TO 35
      DO 30 I=1,60
30    LINE1(I)=BLANK
35    NS1=NS+1
      DO 40 I=1,NN
40    LINE1(NS+I)=LINE(I:I)
10    WRITE(IO,202) (LINE1(I),I=1,60)
CCC   WRITE(IO,203) DATE,TIME
      WRITE(IO,603)
60    CONTINUE
      RETURN
C
C       FORMAT STATEMENTS
101   FORMAT (A60)
201   FORMAT(1H1//////////1H ,16X,68(1H*))
202   FORMAT(1H ,16X,1H*,66X,1H*/1H ,16X,1H*,3X,60A1,3X,1H*)
203   FORMAT(1H ,16X,1H*,66X,1H*/1H ,16X,1H*,12X,'PROGRAM RUN ON ',
     1 A16,' AT ',A8,11X,1H*/1H ,16X,1H*,66X,1H*/1H ,16X,68(1H*)
     2 /1H1)
603   FORMAT(1H ,16X,1H*,66X,1H*/1H ,16X,1H*,66X,1H*/1H ,16X,68(1H*)
     2 /1H1)
      END
C
C       ******************************************************
C      *                                                      *
C      *                  SUBROUTINE OFILE                    *
C      *                                                      *
C      *            VERSION CURRENT AS OF 10/01/87            *
C      *                                                      *
C       ******************************************************
C
      SUBROUTINE OFILE
      CHARACTER*50 IFNAME,OFNAME
      CHARACTER*1 STAR
      COMMON /IOUNIT/ IN,IO
      DATA STAR/'*'/
      IN=15
      IO=16
      WRITE(*,5)
      READ(*,7) IFNAME
      WRITE(*,6)
      READ(*,7) OFNAME
      OPEN (IN,FILE=IFNAME,STATUS='OLD')
      OPEN (IO,FILE=OFNAME)
      RETURN
C
C       FORMAT STATEMENTS
5     FORMAT(5X,'TYPE IN INPUT FILE NAME')
6     FORMAT(5X,'TYPE IN OUTPUT FILE NAME')
7     FORMAT(A50)
      END
C
C       ******************************************************
C      *                                                      *
C      *                  SUBROUTINE PLOT1D                   *
C      *                                                      *
C      *            VERSION CURRENT AS OF 10/01/87            *
C      *                                                      *
C       ******************************************************
C
      SUBROUTINE PLOT1D (X,C,NX,T,IT,NT,TUNITS,LUNITS,XSCLP)
      DIMENSION X(NX),C(NX)
      CHARACTER*10 LUNITS,TUNITS
C
C       THIS ROUTINE PLOTS CONCENTRATION VS. DISTANCE AT EACH OF THE
C       TIMES SPECIFIED IN THE INPUT DATA.
C
C       *** PLOT88 VERSION WILL BE AVAILABLE SOON. ***
C
      WRITE(*,*) '*** PLOT ROUTINES FOR PC WILL BE AVAILABLE SOON ***'
      RETURN
      END
C
C       ******************************************************
C      *                                                      *
C      *                  SUBROUTINE PLOT2D                   *
C      *                                                      *
C      *            VERSION CURRENT AS OF 10/01/87            *
C      *                                                      *
C       ******************************************************
C
      SUBROUTINE PLOT2D (XP,YP,CP,TP,DELTA,NX,NY,NXY,NXY2,IT,NT,IPLT,
     1 TUNITS,LUNITS,XSCLP,YSCLP,XPC,YPC,IFLAG)
      CHARACTER*10 TUNITS,LUNITS
      DIMENSION XP(NX),YP(NY),CP(NXY),XPC(50),YPC(50),IFLAG(NXY2)
C
C       THIS ROUTINE INITIALIZES A CONTOUR PLOT ON THE RECTANGULAR GRID
C       DEFINED IN THE X-Y PLANE BY THE X AND Y VALUES READ IN. ONE
C       SUBPLOT IS GENERATED FOR EACH TIME VALUE.  THE ROUTINE USES
C       DISSPLA (ISCO) SOFTWARE SUBROUTINE CALLS.
C
C       *** PLOT88 VERSION WILL BE AVAILABLE SOON. ***
C
      RETURN
      END
C
C       ******************************************************
C      *                                                      *
C      *                  SUBROUTINE PLOT3D                   *
C      *                                                      *
C      *            VERSION CURRENT AS OF 10/01/87            *
C      *                                                      *
C       ******************************************************
C
      SUBROUTINE PLOT3D (XP,YP,ZP,CP,TP,DELTA,NX,NY,NXY,NXY2,IZ,NZ,IPLT,
     1 TUNITS,LUNITS,XSCLP,YSCLP,XPC,YPC,IFLAG)
      CHARACTER*10 TUNITS,LUNITS
      DIMENSION XP(NX),YP(NY),CP(NXY),XPC(50),YPC(50),IFLAG(NXY2)
C
C       THIS ROUTINE INITIALIZES A CONTOUR PLOT ON THE RECTANGULAR GRID
C       DEFINED IN THE X-Y PLANE BY THE X AND Y VALUES READ IN. ONE
C       SUBPLOT IS GENERATED FOR EACH Z VALUE AND A NEW PLOT IS
C       GENERATED FOR EACH TIME VALUE.  THE ROUTINE USES DISSPLA (ISCO)
C       SOFTWARE SUBROUTINE CALLS.
C
C       *** PLOT88 VERSION WILL BE AVAILABLE SOON. ***
C
      WRITE(*,*) '*** PLOT ROUTINES FOR PC WILL BE AVAILABLE SOON ***'
      RETURN
      END
      BLOCK DATA
      CHARACTER*61 TITLE,T1
      COMMON /IOUNIT/ IN,IO
      COMMON /TITLES/ TITLE(4)
      DATA TITLE/4*'
     1         '/
      END
