本プログラムは、数値計算ハンドブック（オーム社）に基づき
作成させていただきました。

C --- Scalar equation, Ut + A*Ux = 0 ---
      PROGRAM SCALAR
      PARAMETER (N=50)
      DIMENSION A(N),F(N),U(N),UB(N),DF(N)
      OPEN(5,FILE='SCALAR.DAT',STATUS='NEW',ACCESS='SEQUENTIAL')
      DATA NSTEP, NPLOT, ALAM, DX/
     &       40 ,   10 ,  0.75 , 0.05/
C --- ALAM : Courant Number Cdt/dx
C     INITIAL VALUE
      CALL INITVL(U,N,AMAXI)
C --- AKAPPA : dt/dx
      AKAPPA= ALAM/AMAXI
C --- DLT : dt
      DLT=AKAPPA*DX
      WRITE(5,*) 'Dx  = ',DX
      WRITE(5,*) 'Dt  = ',DLT
      WRITE(5,*) 'CFL = ',ALAM
C --- time advance
C
      TIME=0.0
      DO 10 ISTEP=1,NSTEP
         TIME=ISTEP*DLT
           CALL FIRST(A,F,U,DF,N,AKAPPA)
           IF (MOD(ISTEP,NPLOT).EQ.0) THEN
               WRITE(5,*) 'ISTEP=',ISTEP 
               CALL GRAPH(U,N,TIME)
           ENDIF
   10 CONTINUE
      CLOSE(5)
      STOP
      END
C ---------------------------------------
      SUBROUTINE INITVL(U,N,AMAXI)
      DIMENSION U(N)
      DATA    N1,    N2,     UR,    UL/
     &        10,    30,    0.0,   1.0/
C
C --- initial value
C
       DU= (UR-UL)/FLOAT(N2-N1)
         DO 10 J=1, N1
           U(J)=UL
   10 CONTINUE
      DO 20 J=N1+1,N2-1
         U(J)=UL + (J-N1)*DU
   20 CONTINUE
      DO 30 J=N2,N
       U(J)= UR
   30 CONTINUE
      AMAXI=MAX(UL,UR)
      RETURN
      END
C --------------------------------------
      SUBROUTINE BOUND1(U,N)
      DIMENSION  U(N)
        U(1)=2.*U(2) - U(3)
        U(N)=2.*U(N-1)-U(N-2)
      RETURN
      END
C --------------------------------------
      SUBROUTINE BOUND2(U,N)
      DIMENSION  U(N)
        U(2)=2.*U(3)-U(4)
        U(1)=2.*U(2)-U(3)
        U(N-1)=2.*U(N-2)-U(N-3)
        U(N)=2.*U(N-1)-U(N-2)
       RETURN
      END
C --------------------------------------
      SUBROUTINE GRAPH (U,N,TIME)
      DIMENSION U(N)
      WRITE(5,*) 'TIME =' , TIME
      DO 10 J=1,N
      WRITE(5,*) 'N=',J,'     U=',U(J)
   10 CONTINUE
      RETURN
      END

C --------------------------------------
      SUBROUTINE FIRST(A,F,U,DF,N,AKAPPA)
C --------------------------------------
C     First-order upwind scheme
C --------------------------------------
      DIMENSION A(N),F(N),U(N),DF(N)
      DO 10 J=1,N
      A(J)=U(J)
      F(J)=U(J)*U(J)/2.
   10 CONTINUE    
      DO 20 J=2,N-1
      AP=(1+SIGN(1.,A(J)))/2.
      AM=(1-SIGN(1.,A(J)))/2.
      DF(J)=AKAPPA*(AP*(F(J)-F(J-1))
     &             +AM*(F(J+1)-F(J)))
   20 CONTINUE
      DO 30 J=2,N-1
      U(J)=U(J)-DF(J)
   30 CONTINUE
      CALL BOUND1(U,N)
      RETURN
      END
C --------------------------------------
      SUBROUTINE SECOND(F,U,DF,N,AKAPPA)
C --------------------------------------
C second-order central difference scheme
C --------------------------------------
      DIMENSION F(N),U(N),DF(N)
      DO 10 J=1,N
      F(J)=U(J)*U(J)/2.
   10 CONTINUE
      DO 20 J=2,N-1
         DF(J)=AKAPPA*(F(J+1)-F(J-1))/2.
   20 CONTINUE
       DO 30 J=2,N-1
       U(J)=U(J)-DF(J)
   30 CONTINUE
      CALL BOUND1(U,N)
      RETURN
      END
C --------------------------------------
      SUBROUTINE LAXWEND(A,F,U,DF,N,AKAPPA)
C --------------------------------------
C     Lax-Wendroff scheme
C --------------------------------------
      DIMENSION A(N),F(N),U(N),DF(N)
      DO 10 J=1,N
      A(J)=U(J)
      F(J)=U(J)*U(J)/2.
   10 CONTINUE
      ALMDA=AKAPPA/2.
      DO 20 J=2,N-1
      AJP2=(A(J+1)+A(J))/2.
      AJM2=(A(J)+A(J-1))/2.
      DF(J)=ALMDA*(F(J+1)-F(J-1)
     &    -AKAPPA*(AJP2*(F(J+1)-F(J))
     &            -AJM2*(F(J)-F(J-1))))
   20 CONTINUE
      DO 30 J=2,N-1
      U(J)=U(J)-DF(J)
   30 CONTINUE
      CALL BOUND1(U,N)
      RETURN
      END
C --------------------------------------
      SUBROUTINE MACORMAK(F,U,UB,DF,N,AKAPPA)
C --------------------------------------
C     MacCormack scheme
C --------------------------------------
      DIMENSION F(N),U(N),UB(N),DF(N)
C --- Predictor step -------------------
      DO 10 J=1,N
      F(J)=U(J)*U(J)/2.
   10 CONTINUE
      DO 20 J=2,N-1
         DF(J)=AKAPPA*(F(J+1)-F(J))
   20 CONTINUE
       DO 30 J=2,N-1
       UB(J)=U(J)-DF(J)
   30 CONTINUE
      CALL BOUND1(UB,N)
C --- Corrector step -------------------
      DO 11 J=1,N
      F(J)=UB(J)*UB(J)/2.
   11 CONTINUE
       DO 21 J=2,N-1
         DF(J)=AKAPPA*(F(J)-F(J-1))
   21 CONTINUE
       DO 31 J=2,N-1
       U(J)=(U(J)+UB(J)-DF(J))/2.
   31 CONTINUE
      CALL BOUND1(U,N)
      RETURN
      END