C     gfortran e:\1DCOMPRE\1DCOMPRE20170116_SEMIRAGU.f -O3 e:\1DCOMPRE\1DCOMPRE20170114.exe -O3
C######################################################################
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      PARAMETER(NX=401,ISTEP=200000/1)

C-------------------------- 配 列 の 定 義 ----------------------------

      DIMENSION X(NX),A(NX),R(NX),U(NX),P(NX),T(NX),PX(NX),PO(NX)
      DIMENSION F(NX,3,3),G(NX,3),FD(NX,2,2),FD1(NX,2,2),FD2(NX,2,2)
      DIMENSION RES(NX,2),WK2(NX,1),FF(NX,2,1)
      DIMENSION DD(NX,1),D(NX,2,2),AA(NX)

C-------------------------- γ の 値 ----------------------------------

      GMA=1.4
      GM1=0.5*(GMA-1.0)
      GM2=GMA/(GMA-1.0)

C-------------------------- Re,Pt,DX,DELT の 設 定 -------------------

      RE=100 !1.   !20.    !5.  ! 1.E12   !5.D0
      PT=0.7   !0.01 !0.5    !0.8 !0.5 !0.4
C
      DX=20.0/FLOAT(NX-1)
C
      DELT=0.01*DX  !0.01*DX

C-------------------------- 初 期 値 の 設 定 ----------------------

      DO I=1,NX

      X(I)=-10.+DX*FLOAT(I-1)
      A(I)= 21.-SQRT(20**2-X(I)**2)
      U(I)=0.1

      ENDDO

C=============== 入 口 以 外 タ ン ク 内 圧 に す る ==================

      T(1)=1.0-GM1*U(1)**2
      P(1)=T(1)**GM2
      R(1)=GMA*P(1)/T(1)
C
      DO I=2,NX

      T(I)=1.0-GM1*U(I)**2
      P(I)=PT*T(I)**GM2
      R(I)=GMA*P(I)/T(I)

      ENDDO

C====================== F1=RA, F2=RAU ==============================

      DO I=1,NX

      F(I,1,1)=R(I)*A(I)
      F(I,2,1)=F(I,1,1)*U(I)
      ENDDO

C===================== IT STEP を 回 す ============================

      DO IT=1, ISTEP

C==================== 値 の 更 新 ==================================
      DO I=1,NX

      F(I,1,2)=F(I,1,1)
      F(I,2,2)=F(I,2,1)

      F(I,1,3)=F(I,1,1)
      F(I,2,3)=F(I,2,1)

      ENDDO
      

C      DO I=2,NX-1
C      WK2(I,1)=0.5D0*(F(I+1,2,2)+F(I-1,2,2))
C      ENDDO
C      DO I=2,NX-1
C      F(I,2,2)=WEIGHT*WK2(I,1)+(1.0D0-WEIGHT)*F(I,2,2)
C      ENDDO


C==============時間微分は TVD 3次ランゲクッタ法 を使う==================

      CALL WENO(NX,DX,DELT,D,FD,U,F,  FF)
C      CALL KAWAMURA(NX,NY,DX,DELT,D,FD,U,  F)
C      CALL KAWAKAMI1(NX,NY,DX,DELT,D,FD,U,F,  FF)
C
      DOI=2,NX-1
      F(I,1,1)=F(I,1,3)+FF(I,1,1)
      F(I,2,1)=F(I,2,3)+FF(I,2,1)
      ENDDO

      DO I=2,NX-1
      F(I,1,2)=F(I,1,1)
      F(I,2,2)=F(I,2,1)
      ENDDO

      CALL WENO(NX,DX,DELT,D,FD,U,F,  FF)
C      CALL KAWAKAMI1(NX,NY,DX,DELT,D,FD,U,F,  FF)
      
      DO I=2,NX-1
      F(I,1,1)=3.D0*F(I,1,3)/4.D0+(F(I,1,2)+FF(I,1,1))/4.D0
      F(I,2,1)=3.D0*F(I,2,3)/4.D0+(F(I,2,2)+FF(I,2,1))/4.D0
      ENDDO

      DO I=2,NX-1
      F(I,1,2)=F(I,1,1)
      F(I,2,2)=F(I,2,1)
      ENDDO
      
      CALL WENO(NX,DX,DELT,D,FD,U,F,  FF)
C      CALL KAWAKAMI1(NX,NY,DX,DELT,D,FD,U,F,  FF)
      
      DO I=2,NX-1
      F(I,1,1)=F(I,1,3)/3.D0+2.D0*(F(I,1,2)+FF(I,1,1))/3.D0
      F(I,2,1)=F(I,2,3)/3.D0+2.D0*(F(I,2,2)+FF(I,2,1))/3.D0
      ENDDO

c      CALL WENO(NX,DX,DELT,D,FD,U,F,  FF)
C      CALL KAWAMURA(NX,NY,DX,DELT,D,FD,U,  F)
C      CALL KAWAKAMI1(NX,NY,DX,DELT,D,FD,U,  F)
C      CALL KAWAKAMI2(NX,NY,DX,DELT,D,FD,U,  F)

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%                   圧力についての差分                   %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C================== 2 次 精 度 中 心 差 分 ======================
C
C      DO I=2,NX-1
C      PO(I)=0.5D0*(P(I)-P(I-1))/DX
C      F(I,2,1)=F(I,2,1)-DELT*PO(I)*A(I)
C      ENDDO
C      

C================== 4 次 精 度 中 心 差 分 ======================

C      DO I=3,NX-2
C
C      PO(I)=(-P(I+2)+P(I-2)+8*P(I+1)-8*P(I-1))/(12.D0*DX)
C      F(I,2,1)=F(I,2,1)-DELT*PO(I)*A(I)
C
C      ENDDO

C      PO(2)=0.5D0*(P(3)-P(1))/DX
C      F(2,2,1)=F(2,2,1)-DELT*PO(2)*A(2)
C      PO(NX-1)=0.5D0*(P(NX)-P(NX-2))/DX
C      F(NX-1,2,1)=F(NX-1,2,1)-DELT*PO(NX-1)*A(NX-1)

C================== 6 次 精 度 中 心 差 分 ======================
CC
C      DO I=4,NX-3
CC
C      PO(I)=(45.D0*(P(I+1)-P(I-1))-9.D0*(P(I+2)-P(I-2))
C     >            + P(I+3)-P(I-3)                    )/(60.D0*DX) 
C      AA(I)=DELT*PO(I)*A(I)
C      F(I,2,1)=F(I,2,1)-DELT*PO(I)*A(I)
C
C      ENDDO
CC      
C      PO(3   )=(8.D0*(P(4   )-P(2   ))-P(5   )+P(1   ))/(12.D0*DX)
C      PO(2   )=0.5D0*(P(3   )-P(1   ))/DX
C      PO(NX-2)=(8.D0*(P(NX-1)-P(NX-3))-P(NX  )+P(NX-4))/(12.D0*DX)
C      PO(NX-1)=0.5D0*(P(NX  )-P(NX-2))/DX
C
C      AA(3)=DELT*PO(3)*A(3)
C      AA(2)=DELT*PO(2)*A(2)
C      AA(NX-2)=DELT*PO(NX-2)*A(NX-2)
C      AA(NX-1)=DELT*PO(NX-1)*A(NX-1)
CC
C      F(3,2,1)=F(3,2,1)-DELT*PO(3)*A(3)
C      F(2,2,1)=F(2,2,1)-DELT*PO(2)*A(3)
C      F(NX-1,2,1)=F(NX-1,2,1)-DELT*PO(NX-1)*A(NX-1)
C      F(NX-2,2,1)=F(NX-2,2,1)-DELT*PO(NX-2)*A(NX-2)
C
C##################  境  界  条  件  ##########################

      U(2)   =F(2   ,2,1)/F(2   ,1,1)
      U(NX-1)=F(NX-1,2,1)/F(NX-1,1,1)
      T(2)=1.0-GM1*U(2)**2
      T(NX-1)=1.0-GM1*U(NX-1)**2
C
      U(1)=U(2)
      T(1)=T(2)
      P(1)=T(1)**GM2
      R(1)=GMA*P(1)/T(1)
      F(1,1,1)=R(1)*A(1)
      F(1,2,1)=F(1,1,1)*U(1)
C
      U(NX)=U(NX-1)
      T(NX)=T(NX-1)
      P(NX)=PT*T(NX)**GM2
      R(NX)=GMA*P(NX)/T(NX)
      F(NX,1,1)=R(NX)*A(NX)
      F(NX,2,1)=F(NX-1,2,1)

C      F(NX,2,1)=F(NX,1,1)*U(NX)

C=========================== 値 の 更 新 ===========================

      DO I=1,NX
      F(I,2,2)=F(I,2,1)
      ENDDO

C####################### 粘 性 項 の 計 算 ##########################

      DO KK=1,3
C-------------------------------------------------------------------

      DO I=2,NX-1

      DIV=1.+DELT/F(I,1,1)/A(I)/RE
      F(I,2,1)=F(I,2,2)/DIV

      ENDDO

C##################  境  界  条  件  ##########################

      U(2)   =F(2   ,2,1)/F(2   ,1,1)
      U(NX-1)=F(NX-1,2,1)/F(NX-1,1,1)
      T(2)=1.0-GM1*U(2)**2
      T(NX-1)=1.0-GM1*U(NX-1)**2
C
      U(1)=U(2)
      T(1)=T(2)
      P(1)=T(1)**GM2
      R(1)=GMA*P(1)/T(1)
      F(1,1,1)=R(1)*A(1)
      F(1,2,1)=F(1,1,1)*U(1)

      U(NX)=U(NX-1)
      T(NX)=T(NX-1)
      P(NX)=PT*T(NX)**GM2
      R(NX)=GMA*P(NX)/T(NX)
      F(NX,1,1)=R(NX)*A(NX)
C      F(NX,2,1)=F(NX,1,1)*U(NX)
      F(NX,2,1)=F(NX-1,2,1)

C####### KK STEPの終わり##########

      ENDDO

C$$$$$$$$$$$$$$ F1 , F2 を そ れ ぞ れ の 値 に 分 配 $$$$$$$$$$$$$$$$$$

      DO I=2,NX-1
      R(I)=F(I,1,1)/A(I)
      U(I)=F(I,2,1)/F(I,1,1)
      T(I)=1.0-GM1*U(I)**2
      P(I)=MIN(R(I)*T(I)/GMA,1.0)
C      P(I)=R(I)*T(I)/GMA
      ENDDO

C$$$$$$$$$$$$$$ F1 F2 PO を コマンドプロンプトに出力 $$$$$$$$$$$$$$$$$$$

      WRITE(*,*)IT,F(1,1,1),F(1,2,1)

C####### I STEPの終わり##########

      ENDDO

C******************** EXCEL と TECPLOT に 出 力 **********************
      OPEN(10, FILE='COMPRE1d2.CSV')
      
      DO I=1,NX
C      WRITE(10,*)I,',',X(I),',',F(I,1,1),',',F(I,2,1),',',
C     >            PO(I),',',AA(I)
      WRITE(10,*) I,',',X(I),',',U(I),',',PO(I),',',F(I,2,1)
     >            ,',',F(I,1,1),',',R(I),',',T(I),',',U(I)/T(I)
      ENDDO
      
      CLOSE(10)
C=====================================================================
      OPEN(25, FILE='COMPRE1d.DAT')

      WRITE (25,*) 'VARIABLES = "X","R","U","P","T","M",'
      WRITE (25,*) 'ZONE F=POINT, I=', NX

      DO I=1,NX
      WRITE(25,*) X(I),',',R(I),',',U(I),',',P(I),',',T(I),',',
     >            U(I)/T(I)
      ENDDO
      
      CLOSE(25)
C*********************************************************************
      STOP
      END


C#####################################################################
      SUBROUTINE WENO(NX,DX,DELT,D,FD,U,F,  FF)
C#####################################################################
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(NX,3,3),FD(NX,2,2),D(NX,2,2),U(NX),FF(NX,2,1)

C========================== W E N O 法 == RA =======================

C##############################後退差分###############################
      DO I=2,NX
      D(I,1,1)=(F(I,1,2)-F(I-1,1,2))
C      DO I=2,NX-1
C      UP=0.5D0*(U(I  )+U(I+1))
C      UM=0.5D0*(U(I-1)+U(I  ))
C      D(I,1,1)=(UP*F(I,1,2)-UM*F(I-1,1,2))/DX
      ENDDO
C#####################################################################
C      DO I=4,NX-2
      DO I=4,NX-3
C
      FUI1= D(I-2,1,1)/3.D0-7.D0*D(I-1,1,1)/6.D0+11.D0*D(I  ,1,1)/6.D0
      FUI2=-D(I-1,1,1)/6.D0+5.D0*D(I  ,1,1)/6.D0+      D(I+1,1,1)/3.D0
      FUI3= D(I  ,1,1)/3.D0+5.D0*D(I+1,1,1)/6.D0-      D(I+2,1,1)/6.D0
C
      S1=13.D0*(  D(I-2,1,1)-2.D0*D(I-1,1,1)+     D(I  ,1,1))**2/12.D0
     >        +(  D(I-2,1,1)-4.D0*D(I-1,1,1)+3.D0*D(I  ,1,1))**2/ 4.D0
C
      S2=13.D0*(  D(I-1,1,1)-2.D0*D(I  ,1,1)+     D(I+1,1,1))**2/12.D0
     >        +(  D(I-1,1,1)                     -D(I+1,1,1))**2/ 4.D0
C
      S3=13.D0*(  D(I  ,1,1)-2.D0*D(I+1,1,1)+     D(I+2,1,1))**2/12.D0
     >     +(3.D0*D(I  ,1,1)-4.D0*D(I+1,1,1)+     D(I+2,1,1))**2/ 4.D0
C
      E=1.D-15

      A1=0.1/(S1+E)**2
      A2=0.6/(S2+E)**2
      A3=0.3/(S3+E)**2
C
      W1=A1/(A1+A2+A3)
      W2=A2/(A1+A2+A3)
      W3=A3/(A1+A2+A3)
C
      FD(I,1,1)=(W1*FUI1+W2*FUI2+W3*FUI3)

      ENDDO
C#################################前進差分############################
      DO I=1,NX-1
      D(I,2,1)=(F(I+1,1,2)-F(I,1,2))

C      DO I=2,NX-1
C      UP=0.5D0*(U(I  )+U(I+1))
C      UM=0.5D0*(U(I+1)+U(I  ))
C      D(I,2,1)=(UP*F(I+1,1,2)-UM*F(I,1,2))/DX
      ENDDO
C#####################################################################
      DO I=4,NX-3
C      DO I=3,NX-3
C
      FUI1=-   D(I-2,2,1)/6.D0+5.D0*D(I-1,2,1)/6.D0+D(I  ,2,1)/3.D0
      FUI2=    D(I-1,2,1)/3.D0+5.D0*D(I  ,2,1)/6.D0-D(I+1,2,1)/6.D0
      FUI3=11.*D(I  ,2,1)/6.D0-7.D0*D(I+1,2,1)/6.D0+D(I+2,2,1)/3.D0

      S1=13.D0*(  D(I-2,2,1)-2.D0*D(I-1,2,1)+     D(I  ,2,1))**2/12.D0
     >        +(  D(I-2,2,1)-4.D0*D(I-1,2,1)+3.D0*D(I  ,2,1))**2/4.D0
C
      S2=13.D0*(  D(I-1,2,1)-2.D0*D(I  ,2,1)+     D(I+1,2,1))**2/12.D0
     >        +(  D(I-1,2,1)                -     D(I+1,2,1))**2/4.D0
C
      S3=13.D0*(  D(I  ,2,1)-2.D0*D(I+1,2,1)+     D(I+2,2,1))**2/12.D0
     >     +(3.D0*D(I  ,2,1)-4.D0*D(I+1,2,1)+     D(I+2,2,1))**2/4.D0
C
      E=1.D-15

      A1=0.3/(S1+E)**2
      A2=0.6/(S2+E)**2
      A3=0.1/(S3+E)**2
C
      W1=A1/(A1+A2+A3)
      W2=A2/(A1+A2+A3)
      W3=A3/(A1+A2+A3)
C
      FD(I,2,1)=(W1*FUI1+W2*FUI2+W3*FUI3)

      ENDDO

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C$$$$                  4からNX-3までをまとめる                     $$$$
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      DO I=4,NX-3
      FF(I,1,1)=-DELT*U(I)*(FD(I,1,1)+FD(I,2,1))/DX
      ENDDO

C================= I=3    の 部 分 の 差 分 ========================

      UP13=0.5D0*(U(3)+U(4))
      UM13=0.5D0*(U(2)+U(3))
      DF13=(UP13*F(3,1,2)-UM13*F(2,1,2))/DX
C      F(3,1,1)=F(3,1,2)-DELT*DF13
      FF(3,1,1)=-DELT*DF13

C================= I=2    の 部 分 の 差 分 ========================

      UP12=0.5D0*(U(2)+U(3))
      UM12=0.5D0*(U(1)+U(2))
      DF12=(UP12*F(2,1,2)-UM12*F(1,1,2))/DX
C      F(2,1,1)=F(2,1,2)-DELT*DF12
      FF(2,1,1)=-DELT*DF12

CC================= I=NX-3 の 部 分 の 差 分 ========================
C
C      UP2NX=0.5D0*(U(NX-3)+U(NX-2))
C      UM2NX=0.5D0*(U(NX-4)+U(NX-3))
C      DF2NX=(UP2NX*F(NX-3,1,2)-UM2NX*F(NX-4,1,2))/DX
CC      F(NX-2,1,1)=F(NX-2,1,2)-DELT*DF2NX
C      FF(NX-3,1,1)=-DELT*DF2NX
C
C================= I=NX-2 の 部 分 の 差 分 ========================

      UP1NX=0.5D0*(U(NX-2)+U(NX-1))
      UM1NX=0.5D0*(U(NX-3)+U(NX-2))
      DF1NX=(UP1NX*F(NX-2,1,2)-UM1NX*F(NX-3,1,2))/DX
C      F(NX-2,1,1)=F(NX-2,1,2)-DELT*DF1NX
      FF(NX-2,1,1)=-DELT*DF1NX

C================= I=NX-1 の 部 分 の 差 分 ========================

      UP1NX1=0.5D0*(U(NX-1)+U(NX  ))
      UM1NX1=0.5D0*(U(NX-2)+U(NX-1))
      DF1NX1=(UP1NX1*F(NX-1,1,2)-UM1NX1*F(NX-2,1,2))/DX
C      F(NX-1,1,1)=F(NX-1,1,2)-DELT*DF1NX1
      FF(NX-1,1,1)=-DELT*DF1NX1
      
C========================= W E N O 法 == RAU =========================

C########################## 後 退 差 分 ###############################

      DO I=2,NX
      D(I,1,2)=(F(I,2,2)-F(I-1,2,2))      
C      DO I=2,NX-1
C      UP=0.5D0*(U(I  )+U(I+1))
C      UM=0.5D0*(U(I-1)+U(I  ))
C      D(I,1,2)=(UP*F(I,2,2)-UM*F(I-1,2,2))/DX      

      ENDDO

C###################################################################
      DO I=4,NX-3
C      DO I=4,NX-2
C######################################################################

      FUI1=    D(I-2,1,2)/3.D0-7.D0*D(I-1,1,2)/6.D0+11.*D(I  ,1,2)/6.D0
      FUI2=   -D(I-1,1,2)/6.D0+5.D0*D(I  ,1,2)/6.D0+    D(I+1,1,2)/3.D0
      FUI3=    D(I  ,1,2)/3.D0+5.D0*D(I+1,1,2)/6.D0-    D(I+2,1,2)/6.D0
C
      S1=13.D0*(  D(I-2,1,2)-2.D0*D(I-1,1,2)+     D(I  ,1,2))**2/12.D0
     >        +(  D(I-2,1,2)-4.D0*D(I-1,1,2)+3.D0*D(I  ,1,2))**2/ 4.D0
C                                                          
      S2=13.D0*(  D(I-1,1,2)-2.D0*D(I  ,1,2)+     D(I+1,1,2))**2/12.D0
     >        +(  D(I-1,1,2)                -     D(I+1,1,2))**2/ 4.D0
C                                                          
      S3=13.D0*(  D(I  ,1,2)-2.D0*D(I+1,1,2)+     D(I+2,1,2))**2/12.D0
     >     +(3.D0*D(I  ,1,2)-4.D0*D(I+1,1,2)+     D(I+2,1,2))**2/ 4.D0
C
      E=1.D-15

      A1=0.1/(S1+E)**2
      A2=0.6/(S2+E)**2
      A3=0.3/(S3+E)**2
C
      W1=A1/(A1+A2+A3)
      W2=A2/(A1+A2+A3)
      W3=A3/(A1+A2+A3)
C
      FD(I,1,2)=(W1*FUI1+W2*FUI2+W3*FUI3)

      ENDDO

C######################### 前 進 差 分 ###############################
      DO I=1,NX-1
      D(I,2,2)=(F(I+1,2,2)-F(I,2,2))

C      DO I=2,NX-1
C      UP=0.5D0*(U(I  )+U(I+1))
C      UM=0.5D0*(U(I+1)+U(I  ))
C      D(I,2,2)=(UP*F(I+1,2,2)-UM*F(I,2,2))/DX
      ENDDO
C######################################################################
      DO I=4,NX-3
C      DO I=3,NX-3
C######################################################################

      FUI1=-   D(I-2,2,2)/6.D0+5.D0*D(I-1,2,2)/6.D0+D(I  ,2,2)/3.D0
      FUI2=    D(I-1,2,2)/3.D0+5.D0*D(I  ,2,2)/6.D0-D(I+1,2,2)/6.D0
      FUI3=11.*D(I  ,2,2)/6.D0-7.D0*D(I+1,2,2)/6.D0+D(I+2,2,2)/3.D0
C
      S1=13.D0*(  D(I-2,2,2)-2.D0*D(I-1,2,2)+     D(I  ,2,2))**2/12.D0
     >        +(  D(I-2,2,2)-4.D0*D(I-1,2,2)+3.D0*D(I  ,2,2))**2/ 4.D0
C                                                                  
      S2=13.D0*(  D(I-1,2,2)-2.D0*D(I  ,2,2)+     D(I+1,2,2))**2/12.D0
     >        +(  D(I-1,2,2)                -     D(I+1,2,2))**2/ 4.D0
C                                                                  
      S3=13.D0*(  D(I  ,2,2)-2.D0*D(I+1,2,2)+     D(I+2,2,2))**2/12.D0
     >     +(3.D0*D(I  ,2,2)-4.D0*D(I+1,2,2)+     D(I+2,2,2))**2/ 4.D0

      E=1.D-15

      A1=0.3/(S1+E)**2
      A2=0.6/(S2+E)**2
      A3=0.1/(S3+E)**2

      W1=A1/(A1+A2+A3)
      W2=A2/(A1+A2+A3)
      W3=A3/(A1+A2+A3)

      FD(I,2,2)=(W1*FUI1+W2*FUI2+W3*FUI3)

      ENDDO
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
C$$$$                  4からNX-3までをまとめる                     $$$$
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      DO I=4,NX-3
      FF(I,2,1)=-DELT*U(I)*(FD(I,1,2)+FD(I,2,2))/DX
      ENDDO

C================= I=3    の 部 分 の 差 分 ========================

      UP23=0.5D0*(U(3)+U(4))
      UM23=0.5D0*(U(2)+U(3))
      DF23=(UP23*F(3,2,2)-UM23*F(2,2,2))/DX
C      F(3,2,1)=F(3,2,2)-DELT*DF23
      FF(3,2,1)=-DELT*DF23

C================= I=2    の 部 分 の 差 分 ========================

      UP22=0.5D0*(U(2)+U(3))
      UM22=0.5D0*(U(1)+U(2))
      DF22=(UP22*F(2,2,2)-UM22*F(1,2,2))/DX
C      F(2,2,1)=F(2,2,2)-DELT*DF22
      FF(2,2,1)=-DELT*DF22

CC================= I=NX-3 の 部 分 の 差 分 ========================
C
C      UP2NX=0.5D0*(U(NX-3)+U(NX-2))
C      UM2NX=0.5D0*(U(NX-4)+U(NX-3))
C      DF2NX=(UP2NX*F(NX-3,2,2)-UM2NX*F(NX-4,2,2))/DX
CC      F(NX-2,2,1)=F(NX-2,2,2)-DELT*DF2NX
C      FF(NX-3,2,1)=-DELT*DF2NX
C
C================= I=NX-2 の 部 分 の 差 分 ========================

      UP2NX=0.5D0*(U(NX-2)+U(NX-1))
      UM2NX=0.5D0*(U(NX-3)+U(NX-2))
      DF2NX=(UP2NX*F(NX-2,2,2)-UM2NX*F(NX-3,2,2))/DX
C      F(NX-2,2,1)=F(NX-2,2,2)-DELT*DF2NX
      FF(NX-2,2,1)=-DELT*DF2NX

C================= I=NX-1 の 部 分 の 差 分 ========================

      UP2NX1=0.5D0*(U(NX-1)+U(NX  ))
      UM2NX1=0.5D0*(U(NX-2)+U(NX-1))
      DF2NX1=(UP2NX1*F(NX-1,2,2)-UM2NX1*F(NX-2,2,2))/DX
C      F(NX-1,2,1)=F(NX-1,2,2)-DELT*DF2NX
      FF(NX-1,2,1)=-DELT*DF2NX1

      RETURN
      END
C#####################################################################
      SUBROUTINE KAWAKAMI1(NX,NY,DX,DELT,D,FD,U,F,  FF)
C#####################################################################
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(NX,3,3),FD(NX,2,2),D(NX,2,2),U(NX),FF(NX,2,1)

      DO I=2,NX-1
C
      UP=0.5D0*(U(I  )+U(I+1))
      UM=0.5D0*(U(I-1)+U(I  ))

      DF1=(UP*F(I,1,2)-UM*F(I-1,1,2))/DX
      DF2=(UP*F(I,2,2)-UM*F(I-1,2,2))/DX
C
c      F(I,1,1)=F(I,1,2)-DELT*DF1
c      F(I,2,1)=F(I,2,2)-DELT*DF2
      FF(I,1,1)=-DELT*DF1
      FF(I,2,1)=-DELT*DF2

      ENDDO

      RETURN
      END

C#####################################################################
      SUBROUTINE KAWAKAMI2(NX,NY,DX,DELT,D,FD,U,  F)
C#####################################################################
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(NX,3,3),FD(NX,2,2),D(NX,2,2),U(NX)

      DO I=2,NX-1
C
      UP=0.5D0*(U(I  )+ABS(U(I)))
      UM=0.5D0*(U(I  )-ABS(U(I)))

      DF1=(UP*(F(I,1,2)-F(I-1,1,2))+UM*(F(I+1,1,2)-F(I,1,2)))/DX
      DF2=(UP*(F(I,2,2)-F(I-1,2,2))+UM*(F(I+1,2,2)-F(I,2,2)))/DX
C
      F(I,1,1)=F(I,1,2)-DELT*DF1
      F(I,2,1)=F(I,2,2)-DELT*DF2

      ENDDO

      RETURN
      END
C#####################################################################
      SUBROUTINE KAWAMURA(NX,NY,DX,DELT,D,FD,U,  F)
C#####################################################################
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(NX,3,3),FD(NX,2,2),D(NX,2,2),U(NX)

      DO I=3,NX-2
C
      DF13=(      U(I+2)*F(I+2,1,2)+U(I-2)*F(I-2,1,2)
     >      -4.*( U(I+1)*F(I+1,1,2)+U(I-1)*F(I-1,1,2))
     >      +6.*( U(I  )*F(I  ,1,2))             )/(4.D0 *DX)
      DF14=(     -U(I+2)*F(I+2,1,2)+U(I-2)*F(I-2,1,2)
     >      +8.*(-U(I-1)*F(I-1,1,2)+U(I+1)*F(I+1,1,2))  )/(12.D0*DX)
     >      +DF13
      DF23=(      U(I+2)*F(I+2,2,2)+U(I-2)*F(I-2,2,2)
     >      -4.*( U(I+1)*F(I+1,2,2)+U(I-1)*F(I-1,2,2))
     >      +6.*( U(I  )*F(I  ,2,2))             )/(4.D0 *DX)
      DF24=(     -U(I+2)*F(I+2,2,2)+U(I-2)*F(I-2,2,2)
     >      +8.*(-U(I-1)*F(I-1,2,2)+U(I+1)*F(I+1,2,2))  )/(12.D0*DX)
     >      +DF23
C
      F(I,1,1)=F(I,1,2)-DELT*DF14
      F(I,2,1)=F(I,2,2)-DELT*DF24

      ENDDO

      UP12=0.5D0*(U(2)+U(3))
      UM12=0.5D0*(U(1)+U(2))
      DF12=(UP12*F(2,1,2)-UM12*F(1,1,2))/DX
      F(2,1,1)=F(2,1,2)-DELT*DF12

      UP1NX=0.5D0*(U(NX-1)+U(NX  ))
      UM1NX=0.5D0*(U(NX-2)+U(NX-1))
      DF1NX=(UP1NX*F(NX-1,1,2)-UM1NX*F(NX-2,1,2))/DX
      F(NX-1,1,1)=F(NX-1,1,2)-DELT*DF1NX

      UP22=0.5D0*(U(2)+U(3))
      UM22=0.5D0*(U(1)+U(2))
      DF22=(UP22*F(2,2,2)-UM22*F(1,2,2))/DX
      F(2,2,1)=F(2,2,2)-DELT*DF22

      UP2NX=0.5D0*(U(NX-1)+U(NX  ))
      UM2NX=0.5D0*(U(NX-2)+U(NX-1))
      DF2NX=(UP2NX*F(NX-1,2,2)-UM2NX*F(NX-2,2,2))/DX
      F(NX-1,2,1)=F(NX-1,2,2)-DELT*DF2NX

      RETURN
      END