! 間隙水圧分布u(z,t)を計算

! gfortran -o u u.f90
! ./u>u.d
! gnuplot> set zrange [] reverse
! gnuplot> splot 'u.d' with pm3d

implicit none
!---手入力---
real(8) :: cv=200.d0          ! cm^2/day : 圧密係数
real(8) :: h=1000.d0          ! cm       : 層厚
real(8) :: p=100.d0           ! kN/m^2   : 圧密応力
integer :: n=1000             !          : 級数の最大値
integer :: t=5000             ! day      : 経過日数
integer :: dt=30              ! day      : 刻み日数
real(8) :: dz=10.d0           ! cm       : 刻み深さ
!---
real(8),allocatable :: u(:,:) ! kN/m^2   : 間隙水圧分布 u(z,t)
real(8) :: pi=2.*asin(1.)     !          : 円周率
integer i,j,k
allocate(u(0:int(h/dz),0:int(t/dt)))

!---間隙水圧分布の計算---
u(:,:)=0. !初期化
do k=0,int(t/dt)
   do j=0,int(h/dz)
      do i=0,n 
         u(j,k)=u(j,k)+4.*p/((2.*i+1.)*pi)*sin((2.*i+1.)/2*pi*(j*dz)/h)&
                &*exp(-((2.*i+1.)**2.)/4*pi**2.*cv*(k*dt)/(h**2.))
      enddo
   enddo
enddo

!---gnuplot用に出力---
do k=0,int(t/dt)
   do j=0,int(h/dz)
      print*,u(j,k),k*dt,j*dz/100
   enddo
   print*,''
enddo
end
