[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode
program main
  implicit none
  integer,parameter::N=11
  integer t
 real*8::A(N),B(N),C(N),D(N),H(N)
 integer,parameter::Dx=10,Dt=1
  !空间和时间的步长
 real,parameter:: k=0.1
 real,parameter:: lemad=real(k*Dt)/real(Dx**2)
  !含水层参数
 integer i,j
  !声明
    do i=2,N-1
     A(i)=-lemad
     B(i)=1+2*lemad
     C(i)=-lemad
     end do
    write(*,20)lemad
     20	format(/,5X,'lemad=',F9.6)
    do i =1,N
    H(i)=10
    end do
    write (*,*) H(2:N-1)
    H(1)=10
    H(N)=12
    write(*,15)
     15   format(5X,'TIME',5X,'HEAD')!赋初值和边界值
     t = 0
    do j=1,101
    t=t+Dt 
    do  i=2,N-1
       D(i)=H(i)
    D(2)=D(2)+H(1)*lemad
    D(N-1)=D(N-1)+H(N)*lemad
    end do
    call  chase (N,A,B,C,D,H)
     write(*,30) t,H(i)
    30	format(1X,I3,5X,11F9.6)
    end do
    read(*,*)t
    end program
    subroutine chase (N,A,B,C,D,H)
     !追赶法计算三对角线方程组
     implicit none
      integer:: N
     real*8::A(N),B(N),C(N),D(N),H(N)   
     real*8::L(N),U(N)
    integer::i,j
     L(2)=B(2)
     U(2)=D(2)/L(2)
     do  i = 3,N-1
      L(i) = B(i)-A(i)*C(i-1)/L(i-1)
      U(i) =(D(i)-A(i)*U(i-1))/L(i)
   end do
    !回代求H
    H(N-1)=U(N-1)
     do j= 2,N-3
     i=N-1-j
     H(i)=U(i)-C(i)*H(i+1)/L(i)
    end do
    return
    end subroutine chase