      SUBROUTINE INITAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=61)
      COMMON/T/        T,     NSTOP,    NORUN
     1      /Y/   CA(NX),    CB(NX)
     2      /F/  CAT(NX),   CBT(NX)
     3      /S/  CAX(NX),  CAXX(NX),  CBX(NX),  CBXX(NX),
     4      /R/   XA(NX),    XB(NX),       XL,        DX,
     5                DA,        DB,      CA0,       CB0,
     6                 H
     7      /I/       IP
C...
C...  HALF LENGTH
      XL=60.D0
C...
C...  GRID INCREMENT
      DX=XL/DFLOAT(NX-1)
C...
C...  GRID POINTS
      DO 1 I=1,NX
         XA(I)=-XL+DFLOAT(I-1)*DX
         XB(I)=    DFLOAT(I-1)*DX
1     CONTINUE
C...
C...  DIFFUSIVITIES
      DA=1.0D0
      DB=1.0D0
C...
C...  TRANSFER COEFFICIENT
      H=0.1D0
C...
C...  INITIAL CONDITION
      CA0=0.D0
      CB0=1.0D0
      DO 2 I=1,NX
         CA(I)=CA0
         CB(I)=CB0
2     CONTINUE
C...
C...  INITIALIZE COUNTER FOR PRINTED AND PLOTTED SOLUTION
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=61)
      COMMON/T/        T,     NSTOP,    NORUN
     1      /Y/   CA(NX),    CB(NX)
     2      /F/  CAT(NX),   CBT(NX)
     3      /S/  CAX(NX),  CAXX(NX),  CBX(NX),  CBXX(NX),
     4      /R/   XA(NX),    XB(NX),       XL,        DX,
     5                DA,        DB,      CA0,       CB0,
     6                 H
     7      /I/       IP
C...
C...  BOUNDARY CONDITION AT X = -XL
      NL=1
      CA(1)=CA0
C...
C...  BOUNDARY CONDITION AT X = 0 FOR CA
      NU=2
      CAX(NX)=-(H/DA)*(CA(NX)-CB(1))
C...
C...  DERIVATIVE CA
C...               XX
      CALL DSS044(-XL,0.D0,NX,CA,CAX,CAXX,NL,NU)
C...
C...  BOUNDARY AT X = 0 FOR CB
      NL=2
      CBX(1)=(DA/DB)*CAX(NX)
C...
C...  BOUNDARY CONDITION AT X = XL
      NU=1
      CB(NX)=CB0
C...
C...  DERIVATIVE CB
C...               XX
      CALL DSS044(0.D0,XL,NX,CB,CBX,CBXX,NL,NU)
C...
C...  PDE FOR CA
      DO 1 I=1,NX
         CAT(I)=DA*CAXX(I)
1     CONTINUE
C...
C...  PDE FOR CA
      DO 2 I=1,NX
         CBT(I)=DB*CBXX(I)
2     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=61)
      COMMON/T/        T,     NSTOP,    NORUN
     1      /Y/   CA(NX),    CB(NX)
     2      /F/  CAT(NX),   CBT(NX)
     3      /S/  CAX(NX),  CAXX(NX),  CBX(NX),  CBXX(NX),
     4      /R/   XA(NX),    XB(NX),       XL,        DX,
     5                DA,        DB,      CA0,       CB0,
     6                 H
     7      /I/       IP
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IP=IP+1
      IF(IP.EQ.1)THEN
         WRITE(NO,1)
1        FORMAT(7X,'t'
     1   ,1X,'Ca(-60,t)',1X,'Ca(-30,t)',1X,'  Ca(0,t)'
     2   ,1X,'  Cb(0,t)',1X,' Cb(30,t)',1X,' Cb(60,t)')
      END IF            
C...
C...  PRINT THE NUMERICAL SOLUTION
      NIN=30
      WRITE(NO,3)T,(CA(I),I=1,NX,NIN),(CB(I),I=1,NX,NIN)
3     FORMAT(F8.3,6F10.5,/)
C...
C...  PREPARE A FILE FOR PLOTTING
      CALL PLOT1
      RETURN
      END
      SUBROUTINE PLOT1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NX=61)
      COMMON/T/        T,     NSTOP,    NORUN
     1      /Y/   CA(NX),    CB(NX)
     2      /F/  CAT(NX),   CBT(NX)
     3      /S/  CAX(NX),  CAXX(NX),  CBX(NX),  CBXX(NX),
     4      /R/   XA(NX),    XB(NX),       XL,        DX,
     5                DA,        DB,      CA0,       CB0,
     6                 H
     7      /I/       IP
C...
C...  OPEN FILE FOR TOP DRAWER PLOTTING      
      IF(IP.EQ.1)THEN
         OPEN(4,FILE='T.TOP',STATUS='NEW')
C...
C...     WRITE TOP DRAW FILE FOR PLOTTING CONCENTRATION PROFILES
         WRITE(4,17)
17       FORMAT(' SET LIMITS X FROM -60 TO 60 Y FROM 0 TO 1',/,
     1          ' SET FONT DUPLEX')
C...
C...    SET THE CASE NUMBER FOR THE PLOT LABEL
        NCASE=1
      END IF
C...
C...  WRITE THE SOLUTION AT THE PRESCRIBED TIMES, T = 0, 1, 10, 50, 100
      IF((IP.EQ.  1).OR.(IP.EQ.  2).OR.(IP.EQ. 11).OR.
     1   (IP.EQ. 51).OR.(IP.EQ.101))THEN
         WRITE(4,14)(XA(I),CA(I),I=1,NX)
14       FORMAT(2F10.3)
         WRITE(4,14)(XB(I),CB(I),I=1,NX)
         WRITE(4,16)
16       FORMAT(' JOIN 1')
      END IF
C...
C...  LABELS FOR THE PLOT
      IF(IP.EQ.101)WRITE(4,18)NCASE
18    FORMAT(
     1' TITLE 4.75 9.5 "Figure ',I2,
     1': Ca(x,t), Cb(x,t) vs x"'
     2,/,' TITLE LEFT "Ca(x,t), Cb(x,t)"'
     3,/,' TITLE BOTTOM "x"'
     4,/,' TITLE 5.5 0.75 '
     5   '"t = 0, 1, 10, 50, 100"'
     6,/,' TITLE 5.25 0.35 '
     7   '"Da = 1, Db = 1, h = 0.1"')
      RETURN
      END
