デバッグ終わらず

 波動方程式を計算するプログラムを作成しましたが、デバッグが終了しません。


!module globals

module globals

implicit none

integer,save:: n,itrmax,pintv,i

real(8),save:: dt,xl,fr,gr,h0,dx,dh

real(8),allocatable,save:: x(:),u(:),v(:),w(:),h(:),p(:),q(:),pn(:),qn(:)

!public n,dt,xl,dx,fr,gr,itrmax,pintv,h0,dh,x,u,v,w,h,p,q,pn,qn

end module globals



!module subprograms

module subprogs

use globals

implicit none

contains

subroutine set_init(n,dt,xl,dx,fr,gr,itrmax,pintv,h0,dh,&

x,u,v,w,h,p,q,pn,qn)

!implicit none

real(8),allocatable:: x(:),u(:),v(:),w(:),h(:),p(:),q(:),pn(:),qn(:)

integer n,itrmax,pintv

real(8) dt,xl,fr,gr,h0,dx,dh

open(10,file='data.d')

read(10,*) n,itrmax,pintv

read(10,*) dt,xl,fr,gr,h0

close(10)

dx=xl/dble(n-1)

dh=0.1d0*h0

allocate (x(n),u(n),v(n),w(n),h(n),p(n),q(n),pn(n),qn(n))

x(:)=(/(dx*dble(i-1),i=1,n)/)

u(:)=sqrt(gr*h0)*fr

h(:)=h0+dh*exp(-(x(:)-0.5d0*xl)**2/1.0d2)

call uh2pqvw(pn,qn,gr,u,h,v,w,n)

end subroutine set_init


subroutine print_uh(x,h,u,fo)

!implicit none

real(8),allocatable::x(:),h(:),u(:)

allocate (x(n),h(n),u(n))

!real(8) gr

!integer fopen

!fopen = 0

if(fo==1)then

open(20,file='uh.d')

else if (fo==-1) then

close(20)

return

endif

do i=1,n

write(20,'(10e16.8)')x(i),h(i),u(i),u(i)/sqrt(gr*h(i))

enddo

write(20,*)''

end subroutine print_uh



subroutine chk_cno(dx,dt,v,w)

!implicit none

real(8) cno1,cno2,cno,dt,dx

real(8) v(:),w(:)

cno1=maxval(abs(v(:))*dt/dx)

cno2=maxval(abs(w(:))*dt/dx)

cno=max(cno1,cno2)

if(cno>=1)then

write(*,*) 'stop, cno >= 1, cno = ',cno

stop

endif

end subroutine chk_cno



subroutine cm1d(pn,p,i,v,dt,dx)

!implicit none

real(8),allocatable:: v(:),p(:),q(:),w(:),pn(:),qn(:)

real(8) dt,dx,cno

integer i,n

do i=2,n-1

cno=v(i)*dt/dx

if(cno>=0.0d0)then

pn(i)=p(i) - cno * (p(i ) - p(i-1))

else

pn(i)=p(i) - cno * (p(i+1) - p(i ))

endif

enddo

end subroutine cm1d



subroutine bc_thru(pn,qn,p,q,v,w,n,dt,dx)

!implicit none

integer i,n

real(8),allocatable::p(:),q(:),pn(:),qn(:),v(:),w(:)

real(8) cno1,cno2,dt,dx

do i=1,n,n-1

pn(i)=p(i)

qn(i)=q(i)

cno1=v(i)*dt/dx

cno2=w(i)*dt/dx

if(i==1)then

if (cno1 < 0.0d0) pn(i) = p(i) - cno1 * (p(i+1) - p(i))

if (cno2 < 0.0d0) qn(i) = q(i) - cno2 * (q(i+1) - q(i))

endif

enddo

end subroutine bc_thru



subroutine pq2uhvw(pn,qn,gr,u,h,v,w,n)

!implicit none

real(8) ,allocatable:: u(:),p(:),q(:),v(:),w(:),h(:),pn(:),qn(:)

real(8) gr,cno

integer i,n

allocate (x(n),u(n),v(n),w(n),h(n),p(n),q(n),pn(n),qn(n))

do i=1,n

u(i) = p(i)-q(i)

cno = 1.0/2.0*(p(i)+q(i))

h(i)=cno**2/gr

v(i)=u(i)+cno

w(i)=u(i)-cno

enddo

end subroutine pq2uhvw


end module subprogs



!main program
program main

use globals

use subprogs

implicit none

!real(8),allocatable :: x(:),u(:),v(:),w(:),h(:)

!real(8),allocatable :: p(:),q(:),pn(:),qn(:)

!real(8) dt,xl,dx,fr,gr,h0,dh

integer itr,fopen

call set_init(n,dt,xl,dx,fr,gr,itrmax,pintv,h0,dh,&

x,u,v,w,h,p,q,pn,qn)

call print_uh(x,h,u,1)

do itr=1 ,itrmax

call chk_cno(dx,dt,v,w)

do i=2,n-1

call cm1d(pn,p,i,v,dt,dx)

call cm1d(qn,q,i,w,dt,dx)

enddo

call bc_thru(pn,qn,p,q,v,w,n,dt,dx)

call pq2uhvw(pn,qn,gr,u,h,v,w,n)

p(:)=pn(:)
q(:)=qn(:)

if(mod(itr,pintv)==0) call print_uh(x,h,u,0)

enddo

call print_uh(x,h,u,-1)

deallocate(x,u,v,w,h,p,q,pn,qn)

end program main


!data.d

100 200 20

1.0d-1 1.0d2 0.5d0 9.8d0 1.0d-1

ブログ気持玉

クリックして気持ちを伝えよう!

ログインしてクリックすれば、自分のブログへのリンクが付きます。

→ログインへ

なるほど(納得、参考になった、ヘー)
驚いた
面白い
ナイス
ガッツ(がんばれ!)
かわいい

気持玉数 : 0

この記事へのコメント

この記事へのトラックバック