dynapw.f90

fortran程序

program DynaPw

implicit none

integer n,ntime,npoin,resunit
integer ipoin,itime,i
real u0,dens,pi,lamta,ypoin,temp1,temp2,h,delta,time,omega,c,iomega,omega1,deltaT
real*8 t
real,allocatable::pw(:,:)
character*20 text,probn

probn='11'
resunit=1

delta=6.86
ntime=500
deltaT=0.01
h=103.
npoin=h/delta   +1
n=100
pi=3.1416
u0=1.
dens=1000.
omega=2.
c=1440.e0        !把e后的0改为10可以认为水体不可压缩
omega1=pi*c/2./h
open(resunit,file=trim(probn)//'.res')

do omega=5.,5.,0.1        !此处值可以改
!write(text,'(f4.1)')omega
!open(resunit,file=trim(probn)//trim(text)//'.res')

allocate(pw(npoin,ntime))
pw=0.

ypoin=0.

do ipoin=1,npoin
	temp1=0.;temp2=0.
	do i=1,n
		if(i/2*2==i)cycle
		lamta=i*pi/2./h
		t=lamta*lamta-(omega/c)**2

		if(t<-1.e-8)then
		temp1=temp1+sin(lamta*ypoin)/i/(sqrt(-t))
		elseif(t>1.e-8)then	
		temp2=temp2+sin(lamta*ypoin)/i/(sqrt(t))
		else
		temp1=1e2;temp2=1e2
		exit
		end if
	end do

	time=0.
	do itime=1,ntime
	time=time+deltaT
	pw(ipoin,itime)=4*u0*dens/pi*(cos(omega*time)*temp1+sin(omega*time)*temp2)
	end do

	ypoin=ypoin+delta
end do

ypoin=0.
do ipoin=1,npoin	
	time=0.
	write(resunit,*)'ipoin=',ipoin,'   ypoin=',ypoin
	do itime=1,ntime
	time=time+deltaT
	write(resunit,'(10e20.6)')time,-1.*pw(ipoin,itime)
	end do
	ypoin=ypoin+delta
end do

write(resunit,*)'pressure alone height'
ypoin=0.
do ipoin=1,npoin	
	write(resunit,'(10e20.6)')maxval(pw(ipoin,:)),ypoin
	ypoin=ypoin+delta
end do
!close(resunit)
!write(resunit,'(10e20.6)')omega,maxval(pw(npoin,:)),ypoin
deallocate(pw)
end do
close(resunit)
end program DynaPw