C=========================================================================
C     INTEGRATE ADVECTION EQUATION BY EXPLICIT-FTCS METHOD
C     dc/dt + u*dc/dx = 0
C          INITIAL-CONDITION  : c(x,0) = 0.
C          BOUNDARY-CONDITION : c(0,t) = 1.
C=========================================================================
C------------------ MAIN PROGRAM -----------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXNX=1000,MXNT=1000)
      DOUBLE PRECISION C(0:MXNX,0:MXNT)
      CHARACTER CHARTRIAL*2

      CALL INPUT(XX,DX,NX,TT,DT,NT,U,SIGMA,CHARTRIAL)
      CALL SOLVER(MXNX,MXNT,NX,NT,SIGMA,C)
      CALL OUTPUT(MXNX,MXNT,XX,DX,NX,TT,DT,NT,U,SIGMA,C,CHARTRIAL)

      STOP 'NORMAL TERMINATION'
      END
C-------------------------------------------------------------------------
C------------------ SUB ROUTINE ------------------------------------------
      SUBROUTINE INPUT(XX,DX,NX,TT,DT,NT,U,SIGMA,CHARTRIAL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER CHARTRIAL*2
   10 WRITE(*,*) 'Input simulation length XX and simulation time TT -->'
      WRITE(*,*) 'Input simulation length XX'
      READ(*,*) XX
      WRITE(*,*) 'Input simulation time TT'
      READ(*,*) TT
      WRITE(*,*) '     '
      WRITE(*,*) 'Input length increment DX and time increment DT -->'
      WRITE(*,*) 'length increment DX'
      READ(*,*) DX
      WRITE(*,*) 'time increment DT'
      READ(*,*) DT
      WRITE(*,*) '     '
      WRITE(*,*) 'Input advection velocity U -->'
      WRITE(*,*) 'advection velocity U'
      READ(*,*) U
      WRITE(*,*) '     '
      WRITE(*,*) 'Input 2 characters to append to filename -->'
      WRITE(*,*) 'appendix to filename CHARACTER'
      READ(*,'(A2)') CHARTRIAL
      NX = INT(XX/DX)
      NT = INT(TT/DT)
      SIGMA = U*DT/DX
      WRITE(*,*) '     '
      WRITE(*,*) '==========  INPUT PARAMETER  ====================='
      WRITE(*,'(A15,F10.4)') 'XX        = ',XX
      WRITE(*,'(A15,F10.4)') 'DX        = ',DX
      WRITE(*,'(A15,I10)')   'NX        = ',NX
      WRITE(*,'(A15,F10.4)') 'TT        = ',TT
      WRITE(*,'(A15,F10.4)') 'DT        = ',DT
      WRITE(*,'(A15,I10)')   'NT        = ',NT
      WRITE(*,'(A15,F10.4)') 'U         = ',U
      WRITE(*,'(A15,F10.4)') 'SIGMA     = ',SIGMA
      WRITE(*,'(A15,A2)')    'CHARACTER = ',CHARTRIAL
      WRITE(*,*) '=================================================='
      WRITE(*,*) '     '
      WRITE(*,*) 'Is these parameters OK ? ----->'
      WRITE(*,*) 'Yes --- enter 1, No --- enter 0'
      READ(*,*) ID
      IF(ID.EQ.0) GO TO 10
      RETURN
      END



      SUBROUTINE SOLVER(MXNX,MXNT,NX,NT,SIGMA,C)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION C(0:MXNX,0:MXNT)
      IF(NX.GT.MXNX) STOP 'NX > MXNX'
      IF(NT.GT.MXNT) STOP 'NT > MXNT'
C---- RESET
      DO N = 0,NT
        DO J = 0,NX
          C(N,J) = 0.
        END DO
      END DO
C---- ENTRY OF INITIAL-CONDITION AND BOUDARY-CONDITION
      DO J = 0,NX
        C(J,0) = 0.
      END DO
      DO N = 1,NT
        C(0,N) = 1.
      END DO
C---- COMPUTATION OF NUMERICAL SOLUTION BY EXPLICIT-FTCS METHOD
      DO N = 0,NT-1
        DO J = 1,NX-1
          C(J,N+1) = C(J,N) - SIGMA*(C(J+1,N)-C(J-1,N))/2.
        END DO
        C(NX,N+1) = C(NX,N) - SIGMA*(C(NX,N)-C(NX-1,N))
      END DO
      RETURN
      END



      SUBROUTINE OUTPUT(MXNX,MXNT,XX,DX,NX,TT,DT,NT,U,SIGMA,C,CHARTRIAL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION C(0:MXNX,0:MXNT)
      CHARACTER CHARTRIAL*2,INPDATA*13,OUTDATA*14
      INPDATA = './input'//CHARTRIAL//'.dat'
      OUTDATA = './output'//CHARTRIAL//'.dat'
C---- OUTPUT OF INPUT-PARAMETER
      OPEN(1,FILE=INPDATA,STATUS='REPLACE')
      WRITE(1,'(A30,F10.4)') 'SIMULATION LENGTH         = ',XX
      WRITE(1,'(A30,F10.4)') 'LENGTH INCREMENT          = ',DX
      WRITE(1,'(A30,I10)')   'TOTAL STEPS FOR LENGTH    = ',NX
      WRITE(1,'(A30,F10.4)') 'SIMULATION TIME           = ',TT
      WRITE(1,'(A30,F10.4)') 'TIME INCREMENT            = ',DT
      WRITE(1,'(A30,I10)')   'TOTAL STEPS FOR TIME      = ',NT
      WRITE(1,'(A30,F10.4)') 'ADVECTION VELOCITY        = ',U
      WRITE(1,'(A30,F10.4)') 'NON-DIMENSIONAL PARAMETER = ',SIGMA
      CLOSE(1)
C---- OUTPUT OF NUMERICAL SOLUTION
      OPEN(2,FILE=OUTDATA,STATUS='REPLACE')
      DO J = 0,NX
        WRITE(2,'(2F20.8)') J*DX, C(J,NT)
      END DO
      CLOSE(2)
      RETURN
      END