      PROGRAM FTCS
C
C     ANOTHER PROGRAM TO TEST TIME SCHEMES; PARTICULARLY
C     THE CLANCY TIME SCHEME AS OPPOSED TO MINE.  IT IS BEING WRITTEN
C     WITH THE IDEA OF EXTENSIVE EXPERIMENTATION.
C
C     THIS IS THE STANDARD VERSION ON A FORTY POINT GRID WITH 
C     PERIODIC BOUNDARY CONDITIONS
C
C     I AM GOING TO CALCULATE DURING EACH RUN WITH DIFFERENT VARIABLES
C     THE TRANSPORT USING
C
C     1) FORWARD TIME SECOND ORDER CENTERED IN SPACE
C     2) CENTERED TIME SECOND ORDER CENTERED IN SPACE
C     3) CENTERED TIME FOURTH ORDER CENTERED IN SPACE
C
C
      PARAMETER I$=40
      REAL*4 S0(I$), S1(I$)
     1,      F0(I$), F1(I$), F2(I$)
     2,      G0(I$), G1(I$), G2(I$),XX(I$)
      LOGICAL*1 CHNG, SCHEME(3)
C
C     SET UP THE PROBLEM PARAMETERS AND THE INITIAL CONDITIONS
C
      T=0.
      TCOUNT=0.
      TOUT=1.
      TSTOP=50.
      DELTX=1.
      DELTT=.5
      read (5,*) deltt
      U=1.
      A=0.
      SCHEME(1)=.FALSE.
      SCHEME(2)=.TRUE.
      SCHEME(3)=.FALSE.
      PRINT 11
   11 FORMAT(' CHANGE THE INITIAL DATA T/F')
      READ 30, CHNG
   30 FORMAT (3L2)
      IF (.NOT.CHNG) GO TO 2000
      PRINT 12, DELTT, U, A, TOUT, TSTOP
   12 FORMAT(' ENTER THE NEW VALUES UNDER THE OLD VALUES'/
     1 ' DT             U (VELOCITY)   A (DIFFUSION)  '
     2, 'TOUT           TSTOP          '/1X,5G15.5)
      READ 31, DELTT, U, A, TOUT, TSTOP
   31 FORMAT (5G15.5)
      PRINT 15
   15 FORMAT(' WHICH SCHEMES??? FTCS CTCS2 CTCS4  (3L2)')
      READ 30, SCHEME(1), SCHEME(2), SCHEME(3)      
 2000 CONTINUE
      WRITE(6,13) DELTT,U,A,DELTX,TOUT,TSTOP
   13 FORMAT(/' PARAMETER VALUES (CGS UNITS)'/
     1' DT             U (VELOCITY)   A (DIFFUSION)  DX             '
     2,/1X,4G15.5//
     3,' TOUT           TSTOP          '/1X,2G15.5)
C
C     CALCULATE THE ADVECTIVE AND DIFFUSIVE TIME CONSTANTS
C
      IF (U.EQ.0.) THEN
      ADVECTIVE=1.E25
      ELSE
      ADVECTIVE=DELTX/U
      END IF
      IF (A.EQ.0.) THEN
      DIFFUSIVE=1.E25
      ELSE
      DIFFUSIVE=DELTX**2/A
      END IF
      RATIO=ADVECTIVE/DIFFUSIVE
      WRITE(6,14) ADVECTIVE, DIFFUSIVE, RATIO
   14 FORMAT (//
     1' TIME CONSTANTS IF RATIO < 1 THEN ADVECTION IS DOMINANT'/
     1,' ADVECTIVE      DIFFUSIVE      RATIO'/1X,3G15.5)
C
C     INITIAL CONDITIONS
C
      DO I=1,I$
      S0(I)=0.
      F0(I)=0.
      G0(I)=0.
      XX(I)=FLOAT(I)
      END DO
      DO I=5,15
      S0(I)=.5
C      F0(I)=.5
C      G0(I)=.5
      F0(I)=.5-.4/5.*ABS(10.-I)
      G0(I)=.5-.4/5.*ABS(10.-I)
      END DO
C
C     DO THE FIRST STEP FOR ALL SCHEMES
C
      T=T+DELTT
      TCOUNT=TCOUNT+DELTT
      I$1=I$-1
C
C     CLANCY
C
      DO I=2,I$1
      S1(I)=S0(I)
     1     -DELTT*U/(2.*DELTX)*(S0(I+1)-S0(I-1))           
     2     +DELTT*A/DELTX**2*(S0(I+1)-2.*S0(I)+S0(I-1))
      END DO
      I=1
      S1(I)=S0(I)
     1     -DELTT*U/(2.*DELTX)*(S0(I+1)-S0(I$))           
     2     +DELTT*A/DELTX**2*(S0(I+1)-2.*S0(I)+S0(I$))
      I=I$
      S1(I)=S0(I)
     1     -DELTT*U/(2.*DELTX)*(S0(1)-S0(I-1))           
     2     +DELTT*A/DELTX**2*(S0(1)-2.*S0(I)+S0(I-1))
C
C     SCHEME (2)
C
      DO I=2,I$1
      F1(I)=F0(I)
     1     -DELTT*U/(2.*DELTX)*(F0(I+1)-F0(I-1))           
     2     +DELTT*A/DELTX**2*(F0(I+1)-2.*F0(I)+F0(I-1))
      END DO
      I=1
      F1(I)=F0(I)
     1     -DELTT*U/(2.*DELTX)*(F0(I+1)-F0(I$))           
     2     +DELTT*A/DELTX**2*(F0(I+1)-2.*F0(I)+F0(I$))
      I=I$
      F1(I)=F0(I)
     1     -DELTT*U/(2.*DELTX)*(F0(1)-F0(I-1))           
     2     +DELTT*A/DELTX**2*(F0(1)-2.*F0(I)+F0(I-1))
C
C     SCHEME 3
C
      I$2=I$-2
      DO I=3,I$2
      G1(I)=G0(I)
     1     -DELTT*U/(DELTX)*2./3.*(G0(I+1)-G0(I-1))           
     1     +DELTT*U/(DELTX)*1./12.*(G0(I+2)-G0(I-2))           
     2     +DELTT*A/DELTX**2*4./3.*(G0(I+1)-2.*G0(I)+G0(I-1))
     2     -DELTT*A/DELTX**2*1./12.*(G0(I+2)-2.*G0(I)+G0(I-2))
      END DO
      I=1
      G1(I)=G0(I)
     1     -DELTT*U/(DELTX)*2./3.*(G0(I+1)-G0(I$))           
     1     +DELTT*U/(DELTX)*1./12.*(G0(I+2)-G0(I$1))           
     2     +DELTT*A/DELTX**2*4./3.*(G0(I+1)-2.*G0(I)+G0(I$))
     2     -DELTT*A/DELTX**2*1./12.*(G0(I+2)-2.*G0(I)+G0(I$1))
      I=2
      G1(I)=G0(I)
     1     -DELTT*U/(DELTX)*2./3.*(G0(I+1)-G0(I-1))           
     1     +DELTT*U/(DELTX)*1./12.*(G0(I+2)-G0(I$))           
     2     +DELTT*A/DELTX**2*4./3.*(G0(I+1)-2.*G0(I)+G0(I-1))
     2     -DELTT*A/DELTX**2*1./12.*(G0(I+2)-2.*G0(I)+G0(I$))
      I=I$1
      G1(I)=G0(I)
     1     -DELTT*U/(DELTX)*2./3.*(G0(I+1)-G0(I-1))           
     1     +DELTT*U/(DELTX)*1./12.*(G0(1)-G0(I-2))           
     2     +DELTT*A/DELTX**2*4./3.*(G0(I+1)-2.*G0(I)+G0(I-1))
     2     -DELTT*A/DELTX**2*1./12.*(G0(1)-2.*G0(I)+G0(I-2))
      I=I$
      G1(I)=G0(I)
     1     -DELTT*U/(DELTX)*2./3.*(G0(1)-G0(I-1))           
     1     +DELTT*U/(DELTX)*1./12.*(G0(2)-G0(I-2))           
     2     +DELTT*A/DELTX**2*4./3.*(G0(1)-2.*G0(I)+G0(I-1))
     2     -DELTT*A/DELTX**2*1./12.*(G0(2)-2.*G0(I)+G0(I-2))
C
C     LETS DO THE CLANCY PROBLEM (1)
C
C     THE BOUNDARIES ARE DONE ANALYTICALLY
C
C
C     MOVE THE CLANCY FROM THE FIRST TIME STEP
C
      DO I=1,I$
      S0(I)=S1(I)
      END DO
 1000 CONTINUE
      T=T+DELTT
      TCOUNT=TCOUNT+DELTT
      IF (.NOT.SCHEME(1)) GO TO 1100
      DO I=2,I$1
      S1(I)=S0(I)
     1     -DELTT*U/(2.*DELTX)*(S0(I+1)-S0(I-1))           
     2     +DELTT*A/DELTX**2*(S0(I+1)-2.*S0(I)+S0(I-1))
      IF (S1(I).GE.2.) THEN
      WRITE (6,42) T,I
   42 FORMAT(' SCHEME 1 IS UNSTABLE TIME (DAYS)',F10.3,I5)
      SCHEME(1)=.FALSE.
      END IF
      END DO
      I=1
      S1(I)=S0(I)
     1     -DELTT*U/(2.*DELTX)*(S0(I+1)-S0(I$))           
     2     +DELTT*A/DELTX**2*(S0(I+1)-2.*S0(I)+S0(I$))
      I=I$
      S1(I)=S0(I)
     1     -DELTT*U/(2.*DELTX)*(S0(1)-S0(I-1))           
     2     +DELTT*A/DELTX**2*(S0(1)-2.*S0(I)+S0(I-1))
C
C     END CLANCY
C
C
C     SCHEME (2)
C
 1100 CONTINUE
      IF (.NOT.SCHEME(2)) GO TO 1200
      DO I=2,I$1
      F2(I)=F0(I)
     1     -2.*DELTT*U/(2.*DELTX)*(F1(I+1)-F1(I-1))           
     2     +2.*DELTT*A/DELTX**2*(F0(I+1)-2.*F0(I)+F0(I-1))
      IF (f2(I).GE.2.) THEN
      WRITE (6,43) T,I
   43 FORMAT(' SCHEME 2 IS UNSTABLE TIME (DAYS)',F10.3,I5)
      end if
      END DO
      I=1
      F2(I)=F0(I)
     1     -2.*DELTT*U/(2.*DELTX)*(F1(I+1)-F1(I$))           
     2     +2.*DELTT*A/DELTX**2*(F0(I+1)-2.*F0(I)+F0(I$))
      I=I$
      F2(I)=F0(I)
     1     -2.*DELTT*U/(2.*DELTX)*(F1(1)-F1(I-1))           
     2     +2.*DELTT*A/DELTX**2*(F0(1)-2.*F0(I)+F0(I-1))
C
C     SCHEME 3
C
 1200 CONTINUE
      IF (.NOT.SCHEME(3)) GO TO 1300
      I$2=I$-2
      DO I=3,I$2
      G2(I)=G0(I)
     1     -2.*DELTT*U/(DELTX)*2./3.*(G1(I+1)-G1(I-1))           
     1     +2.*DELTT*U/(DELTX)*1./12.*(G1(I+2)-G1(I-2))           
     2     +2.*DELTT*A/DELTX**2*4./3.*(G0(I+1)-2.*G0(I)+G0(I-1))
     2     -2.*DELTT*A/DELTX**2*1./12.*(G0(I+2)-2.*G0(I)+G0(I-2))
      END DO
      I=1
      G2(I)=G0(I)
     1     -2.*DELTT*U/(DELTX)*2./3.*(G1(I+1)-G1(I$))           
     1     +2.*DELTT*U/(DELTX)*1./12.*(G1(I+2)-G1(I$1))           
     2     +2.*DELTT*A/DELTX**2*4./3.*(G0(I+1)-2.*G0(I)+G0(I$))
     2     -2.*DELTT*A/DELTX**2*1./12.*(G0(I+2)-2.*G0(I)+G0(I$1))
      I=2
      G2(I)=G0(I)
     1     -2.*DELTT*U/(DELTX)*2./3.*(G1(I+1)-G1(I-1))           
     1     +2.*DELTT*U/(DELTX)*1./12.*(G1(I+2)-G1(I$))           
     2     +2.*DELTT*A/DELTX**2*4./3.*(G0(I+1)-2.*G0(I)+G0(I-1))
     2     -2.*DELTT*A/DELTX**2*1./12.*(G0(I+2)-2.*G0(I)+G0(I$))
      I=I$1
      G2(I)=G0(I)
     1     -2.*DELTT*U/(DELTX)*2./3.*(G1(I+1)-G1(I-1))           
     1     +2.*DELTT*U/(DELTX)*1./12.*(G1(1)-G1(I-2))           
     2     +2.*DELTT*A/DELTX**2*4./3.*(G0(I+1)-2.*G0(I)+G0(I-1))
     2     -2.*DELTT*A/DELTX**2*1./12.*(G0(1)-2.*G0(I)+G0(I-2))
      I=I$
      G2(I)=G0(I)
     1     -2.*DELTT*U/(DELTX)*2./3.*(G1(1)-G1(I-1))           
     1     +2.*DELTT*U/(DELTX)*1./12.*(G1(2)-G1(I-2))           
     2     +2.*DELTT*A/DELTX**2*4./3.*(G0(1)-2.*G0(I)+G0(I-1))
     2     -2.*DELTT*A/DELTX**2*1./12.*(G0(2)-2.*G0(I)+G0(I-2))
C
C     ALL SCHEMES HAVE BEEN ADVANCED THE SAME
C
C     DO MOVES FOR NEXT TIME STEP
C
 1300 CONTINUE
      DO I=1,I$
      S0(I)=S1(I)
      F0(I)=F1(I)
      G0(I)=G1(I)
      F1(I)=F2(I)
      G1(I)=G2(I)
      END DO
C
C     CHECK FOR PLOTS AND GO TO NEXT TIME STEP
C
      IF (TCOUNT.LT.TOUT) GO TO 1000
      WRITE (6,10) T
      TCOUNT=0.
C      
C     PLOT THE ANSWER 
C
C      CALL USXUSY(1.,40.,.FALSE.,-1.0,1.0,.FALSE.
C     1, 'TRANSPORTED QUANTITY$')
C      IF (SCHEME(1)) CALL IPLOT(XX,S1,I$,'1',.TRUE.)
C      IF (SCHEME(2)) CALL IPLOT(XX,F2,I$,'2',.TRUE.)
C      IF (SCHEME(3)) CALL IPLOT(XX,G2,I$,'3',.TRUE.)
C      CALL IPLOT(XX,ANAY,I$,'A',.FALSE.)
      CALL RSCALE(60,20)
      IF (SCHEME(1)) THEN
      CALL USXUSY(1.,40.,.FALSE.,-.6,1.0,.FALSE.
     1, 'SCHEME 1$')
C      CALL IPLOT(XX,ANAY,I$,'A',.FALSE.)
      CALL IPLOT(XX,S1,I$,'1',.FALSE.)
      END IF
      IF (SCHEME(2)) THEN
      CALL USXUSY(1.,40.,.FALSE.,-.6,1.0,.FALSE.
     1, 'SCHEME 2$')
C      CALL IPLOT(XX,ANAY,I$,'A',.FALSE.)
      CALL IPLOT(XX,F2,I$,'2',.FALSE.)
      END IF
      IF (SCHEME(3)) THEN
      CALL USXUSY(1.,40.,.FALSE.,-.6,1.0,.FALSE.
     1, 'SCHEME 3$')
C      CALL IPLOT(XX,ANAY,I$,'A',.FALSE.)
      CALL IPLOT(XX,G2,I$,'3',.FALSE.)
      END IF
      IF (T.LT.TSTOP) GO TO 1000
      STOP
   10 FORMAT(' THE FOLLOWING PLOTS ARE AT (DAYS)', F10.4)
      END
