PROGRAM UPWIND
C====================================================================
C        UPWIND FINITE DIFFERENCE WITH ADDITIVE VISCOSITY
C N = NUMBER OF ELEMENTS, NN=NUMBER OF NODES
C====================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( MXN=31 )
      DIMENSION A (MXN,MXN), RHS(MXN)
  99  WRITE (*,*)'ENTER U ='
      READ (*,*) U
      IF ( U .EQ. 0. ) STOP
      N = 10
      TL = 1.D0
      DX = TL/FLOAT(N)
      FK = 1.D0
      FKP = DX*U/2.D0
      GBAR = U / FK
      GAMMA = FKP/FK
      WRITE (*,*)' GAMMA=',GAMMA
      C = FKP * ALPHA(GAMMA) + FK
      DX2 = 2.D0* DX
      DXDX = DX * DX
      A1 = - U/DX2 - C / DXDX
      A2 =   2.D0* C / DXDX
      A3 =   U / DX2 - C / DXDX
C--------- RESET MATIX AND RIGHT HAND SIDE
      NN = N + 1
      DO 10 I = 1 , NN
      RHS (I) = 0.
      DO 10 J = 1 , NN
      A (I,J) = 0.
   10 CONTINUE
      DO 20 I = 2 , N
      A (I,I-1) = A1
      A (I,I  ) = A2
      A (I,I+1) = A3
   20 CONTINUE
C--------- IMPLEMENTATION OF BOUNDARY CONDITION
      A (1,1) = 1.D0
      RHS(1) = 1.D0
      A(NN,NN) = 1.D0
      RHS(NN) = 0.
C--------- EVELUATION FOR UNKNOWNS
      CALL SYSTEM ( MXN , NN , A , RHS )
C--------- PRINTING RESULTS
      X = 0.D0
      I = 1
      AA = 1.D0-GAMMA
      BB = 1.D0+GAMMA
      WRITE (*,200)
  200 FORMAT ( 5X,'======== UPWIND PROGRAM =======' )
      EX = EXACT ( X , FK , U, TL)
      WRITE (*,105)
  105 FORMAT ( 1X,6X,'   X-COORD',6X,'     EXACT',' UPWIND SOLUTION',
     * '  WITHOUT UPWIND')
      WRITE (*,100) X , EX, RHS(1), EXACT1 ( I, AA, BB, NN )
  100 FORMAT ( 1X, 4F16.10 )
      DO 30 I = 2 , NN
      X = X + DX
      EX = EXACT( X , FK , U, TL)
      WRITE (*,100) X , EX, RHS(I), EXACT1 ( I, AA, BB, NN )
   30 CONTINUE
      GO TO 99
      END
C
      FUNCTION ALPHA ( G )
      IMPLICIT REAL*8 ( A-H , O-Z )
      ALPHA = 0.D0
      IF ( G .LE. 0.001D0 ) RETURN
      ALPHA = 1.D0
      IF ( G .GT. 15.D0  ) RETURN
      E = DEXP( G )
      COTHG = ( E + 1.D0/E ) / ( E - 1.D0/E )
      ALPHA = COTHG - 1.D0/ G
      RETURN
      END
C
CCCCCCCCCCFUNCTION EXACT ( X , G , U )
C      IMPLICIT REAL*8 ( A-H , O-Z )
C      EXACT = (1.D0/U)*(X-(1-DEXP(G*X))/(1.-DEXP(G)) )
C      RETURN
C     END
C
C
      FUNCTION EXACT ( X , FK , U, TL )
      IMPLICIT REAL*8 ( A-H , O-Z )
C---------- T(X=0)=1. AND T(X=TL)=0.
      EXACT = (DEXP(U*TL/FK)-DEXP(U*X/FK)) / (DEXP(U*TL/FK)-1.)
      RETURN
      END
C
C
      FUNCTION EXACT1 ( I, A, B, N )
      IMPLICIT REAL*8 ( A-H , O-Z )
C---------- T(1)=1. AND T(N)=0.
C---------- EXACT SOLUTION FOR FINITE DIFFRENCE EQUATION
      ABN = (B/A)**N
      EXACT1 = ( ABN - (B/A)**I ) /( ABN - B/A )
      RETURN
      END
C
      SUBROUTINE SYSTEM ( MXN, N, A, R )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A(MXN,MXN), R(MXN)
      DO 40 K = 1 , N
      P = A(K,K)
      DO 10 J = K , N
      A(K,J) = A(K,J) / P
   10 CONTINUE
      R(K) = R(K) / P
      DO 20 I = 1 , N
      IF ( I .EQ. K ) GO TO 20
      AIK = A (I,K)
      DO 30 J = K , N
      A(I,J) = A(I,J)-AIK*A(K,J)
   30 CONTINUE
      R(I) = R(I) - AIK * R(K)
   20 CONTINUE
   40 CONTINUE
      RETURN
      END