PROGRAM ODE1
  
  IMPLICIT NONE
  INTEGER(8),PARAMETER :: N=100
  INTEGER(8) :: I
  REAL(8) :: T,H,DN2,X(0:1),TEND

  TEND=1000 ! 計算終了の t の値
  DN2=0.0001_8 ! ブラントバイサラ振動数の2乗
  H=TEND/N ! t 方向の刻み幅の設定

  OPEN(1,FILE='buoyancy.dat',FORM='UNFORMATTED',ACCESS='STREAM')
  X(0)=1
  X(1)=0
  T=0
  WRITE(1) N
  WRITE(1) H  
  WRITE(1) X
  DO I=0,N-1
     X=X+H*F(X,T) ! オイラー法の例
     T=T+H
     PRINT*,T
     WRITE(1) X     
  END DO
  CLOSE(1)
  
CONTAINS ! 以下の関数副プログラムを内部関数とする
  !----------------    
  FUNCTION F(X,T) ! 関数副プログラムの定義
    REAL(8) :: F(0:1)
    REAL(8) :: T,X(0:1) ! X(0)が z, X(1)が dz/dt
    F(0)=X(1)
    F(1)=-DN2*X(0)
  END FUNCTION F
  !----------------
END PROGRAM ODE1
