program main
  implicit none
  
  integer(kind=4),parameter:: neq=3         ！解く方程式の数
  integer(kind=4),parameter:: nmax=10000
  
  real(kind=8),parameter:: eps=1.0d-10
  
  external ex2,ex1
  common S,S1,P0,PA,M0,L,H,RHO,x0,g,M    ！ex1,ex2と共通の変数

  real(kind=8),dimension(neq):: y0,yn
  real(kind=8) t0,tn,dt
  real(kind=8) S,S1,P0,PA,M0,L,H,RHO,T_0       ！新たな変数の定義
  real(kind=8) x0,v,P,F,M,A,g,E,W,T            ！新たな変数の定義
  
  real(kind=8),dimension(neq,12):: work
  
  integer(kind=4) ierr,i,n
  
  character(len=12) fname1,fname2    ！ファイル名
  
  
  write(*,'(a,$)') 'filename1 ?'         ！ファイル名の入力
  read(*,*) fname1

  write(*,'(a,$)') 'filename2 ?'
  read(*,*) fname2
  
  open(9,file=fname1)
  open(10,file=fname2)
  
  n=0
  t0=0.0d0
  dt=0.01d0
 
  S=0.00031/2.                
  S1=0.0064
  P0=200000
  PA=100000                     ！初期条件などの入力（MKS単位系）
  M0=0.06                       
  L=0.25
  x0=0.2
  H=1.4
  RHO=1000
  T_0=300
  g=9.8
  
  y0(1) = x0               ！ペットボトル内の水面の位置
  y0(2) = 0                ！ペットボトルの速度
  y0(3) = 0.               ！ペットボトルの高さ
  
  V=sqrt(2.*(P0-PA)/RHO)    ！（１）式
  F=RHO*S*V*V               ！（６）式
  M=M0+RHO*S1*(L-X0)        ！（７）式
  A=F/M-g                   ！（８）式
  E=(P0-PA)*S*V             ！ 噴出する水の運動エネルギー
  W=PA*S*V                  ！ 水が大気を押しのける仕事率
  
  write(9,'(i5,1x,f7.3,10d23.15)') n,t0,x0,V,P0,F,M,A,T_0,E,W
  write(10,'(i5,1x,f7.3,3d23.15)') n,t0,y0(2),y0(3)

 do n=1,nmax
     if (M0 < M) then        ！水を噴出中      
        tn=t0+dt
        call runge(neq,ex2,t0,tn,y0,yn,eps,ierr,work)
        t0=tn
        do i=1,neq
           y0(i)=yn(i)
        end do
        
        P=P0*(X0/y0(1))**H      ！（３）式
        V=sqrt(2.*(P-PA)/RHO)   ！（１）式
        F=RHO*S*V*V             ！（６）式
        M=M0+RHO*S1*(L-y0(1))   ！（７）式
        A=F/M-g                 ！（８）式
        T=T_0*(X0/y0(1))**(H-1.)！（４）式
        E=(P-PA)*S*V            ！ 噴出する水の運動エネルギー
        W=PA*S*V                ！ 水が大気を押しのける仕事率

        write(9,'(i5,1x,f7.3,10d23.15)') n,t0,y0(1),V,P,F,M,A,T,E,W
        write(10,'(i5,1x,f7.3,3d23.15)') n,t0,y0(2),y0(3)     
        
     else                   ！水を出し切った後
        tn=t0+dt
        call runge(neq,ex1,t0,tn,y0,yn,eps,ierr,work)
        t0=tn
        do i=1,neq
           y0(i)=yn(i)
        end do
        
        write(10,'(i5,1x,f7.3,4d23.15)') n,t0,y0(2),y0(3)     
     end if
        
     if(y0(3)<0.)exit   ！地上に戻ったら計算終了
  end do
  
  
  close(9)
  close(10)
  
end program main

real(kind=8) function ex2(i,t,y)     ！水を噴出中
  implicit none
  
  integer(kind=4),parameter:: neq=3  ！解く方程式の数
  
  common S,S1,P0,PA,M0,L,H,RHO,x0,g,M  ！mainと共通の変数
  integer(kind=4) i
  real(kind=8) t
  real(kind=8),dimension(neq):: y
  real(kind=8) S,S1,P0,PA,M0,L,H,RHO,x0,g,M  ！変数の定義
       
  if ((L-y(1))<0.) y(1)=L         ！水を最後に噴出したら水の位置はX＝L
  if (i==1) then
     ex2= S/S1*sqrt(2./RHO*(P0*(X0/y(1))**H-PA))  ！（５）式
  else if (i==2) then
     ex2= 2.*S*(P0*(X0/y(1))**H-PA)/(M0+RHO*S1*(L-y(1)))-g  ！ （８）式
  else if (i==3) then
     ex2 = y(2)                      ！（９）式
  else
     ex2=0.0d0
  end if
  
  return
end function ex2


real(kind=8) function ex1(i,t,y)    ！水を出し切った後
  implicit none
  
  integer(kind=4),parameter:: neq=3  ！解く方程式の数
  
   common S,S1,P0,PA,M0,L,H,RHO,x0,g,M  ！mainと共通の変数
  integer(kind=4) i
  real(kind=8) t
  real(kind=8),dimension(neq):: y
  real(kind=8)  S,S1,P0,PA,M0,L,H,RHO,x0,g,M   ！変数の定義
  
     if (i==1) then
        ex1=0.     
     else if (i==2) then
        ex1= -g
     else if (i==3) then
        ex1 = y(2)
     else
        ex1=0.0d0
     end if
  
     return
end function ex1