PROGRAM SHALLOW
  
  IMPLICIT NONE
  INTEGER(8),PARAMETER :: N=100,M=32
  INTEGER(8) :: I
  REAL(8),PARAMETER :: PI=3.1415926535897932384_8
  REAL(8) :: T,H,G,DEPTH,X,DX,DL,ETA0
  REAL(8) :: VAR(0:M,0:1) ! VAR(:,0) は η, VAR(:1) は u
  REAL(8) :: DK1(0:M,0:1),DK2(0:M,0:1),DK3(0:M,0:1),DK4(0:M,0:1),TEND
  
  TEND=2 ! 計算終了の t の値
  DL=1
  DX=DL/M
  H=TEND/N ! t 方向の刻み幅の設定
  ETA0=0.01_8
  G=10
  DEPTH=0.1_8

  VAR=0  
  DO I=0,M-1
     X=(I+0.5_8)*DX
     VAR(I,0)=ETA0*COS(PI*X/DL)
  END DO
  
  OPEN(1,FILE='shallow.dat',FORM='UNFORMATTED',ACCESS='STREAM')
  WRITE(1) M,N
  WRITE(1) DX,H 
  WRITE(1) VAR
  DO I=0,N-1
     DK1=H*DVAR(VAR,T)
     DK2=H*DVAR(VAR+DK1/2,T+H/2)
     DK3=H*DVAR(VAR+DK2/2,T+H/2)
     DK4=H*DVAR(VAR+DK3,T+H)
     VAR=VAR+(DK1+2*DK2+2*DK3+DK4)/6
     T=T+H
     WRITE(1) VAR
  END DO
  CLOSE(1)
  
CONTAINS ! 以下の関数副プログラムを内部関数とする
  !----------------    
  FUNCTION DVAR(VAR,T)
    INTEGER(8) :: I
    REAL(8) :: DVAR(0:M,0:1)
    REAL(8) :: T,VAR(0:M,0:1)

    DVAR(M,0)=0    
    DVAR(0,1)=0
    DVAR(M,1)=0
    
    DO I=0,M-1
       DVAR(I,0)=-DEPTH*(VAR(I+1,1)-VAR(I,1))/DX
    END DO
    DO I=1,M-1
       DVAR(I,1)=-G*(VAR(I,0)-VAR(I-1,0))/DX
    END DO
    
  END FUNCTION DVAR
  !----------------
END PROGRAM SHALLOW
