[子程序] 子程序有问题请大家帮忙看下

module mod_fdm_tool_C3
    use iso_fortran_env
    implicit none
    real(real64), parameter :: pi = 4.d0*datan(1.d0)
contains

subroutine Grid_2D_Type(nx, nt, xmin, xmax, tmin, tmax, hx, ht, x, t)
    implicit none
    integer(int32), intent(in)  :: nx, nt
    real(real64),   intent(in)  :: xmin, xmax, tmin, tmax
    real(real64),   intent(out) :: hx, ht, x(0:nx), t(0:nt)
    integer(int32) :: i, j

    hx = (xmax - xmin) / dble(nx)
    ht = (tmax - tmin) / dble(nt)

    do concurrent (i = 0 : nx)
        x(i) = xmin + dble(i) * hx
    end do
    x(nx) = xmax

    do concurrent (j = 0 : nt)
        t(j) = tmin + dble(j) * ht
    end do
    t(nt) = tmax

    return
end subroutine Grid_2D_Type



!===============================================================================
pure elemental function f_rhs(x,t) result(res)
    implicit none
    real(real64), intent(in) :: x, t
    real(real64) :: res

    res = 0

    return
end function f_rhs

!===============================================================================
pure elemental function Solu_Exact(x,t) result(res)
    implicit none
    real(real64), intent(in) :: x, t
    real(real64) :: res

    res = dexp(x+t)

    return
end function Solu_Exact


subroutine Apply_DBC(x, t,u)
    implicit none
    real(real64), intent(in)  :: x(:), t(:)
    real(real64), intent(out) :: u(:,:)
    integer(int32) :: nt

    nt = size(t)

    u(:, 1) = dexp(x)
    u(1,:) = dexp(t)
    u(nt,:) = dexp(1.0+t)



    return
end subroutine Apply_DBC

!======================================================
subroutine Gene_Alg_Sys_DBC(ht,hx,nx,r,coef1,coef2)
    ! Algebraic system, i.e., Au = f
    implicit none
    integer(int32) ,intent(in):: nx
    real(real64), intent(in)  :: hx,ht
    real(real64), intent(out) :: coef1(nx-1,3),coef2(nx-1,3),r

    r = ht/hx/hx


    ! all the following variable are arraies, not scalars
    coef1(:,1) = -r/2.d0
    coef1(:,2) = 1+r
    coef1(:,3) = -r/2.d0

    coef2(:,1) = r/2.d0
    coef2(:,2) = 1-r
    coef2(:,3) = r/2.d0



    coef1(1,1) = 0.d0
    coef1(nx-1,3) = 0.d0

    coef2(1,1) = 0.d0
    coef2(nx-1,3) = 0.d0



    return
end subroutine Gene_Alg_Sys_DBC


!这个子程序有问题======================================================
subroutine Solver_Crank-Nicolson (coef1,coef2,f,u)

    implicit none
    integer(int32) :: nx,nt,i,j
    real(real64), intent(in)  :: coef1(:,:),coef2(:,:)
    real(real64) :: f(:)
    real(real64) :: u(:,:)
    real(real64) ::r
    real(real64),allocatable::g(:),w(:)


    do j = 1, nt-1
      f(1)=coef2(1,2)*u(1,j)
      f(nx-1)=coef2(nx-1,1)*u(nx-2,j)+coef2(nx-1,2)*u(nx-1,j)
      f(2:nx-1)=coef2(2:nx-2,1)*u(1:nx-3,j)+coef2(2:nx-2,2)*u（2:nx-2,j)+coef2(2,nx-2,3)*u(3:nx-1,j)

      allocate(g(nx-1),w(nx-1))

      g(1)=f(1)/coef1(1,2)
      w(1)=coef1(1,3)/coef1(1,2)

      do i=2:nx-1
        w(i)=coef1(i,2)-coef1(i,1)*w(i-1)
        g(i)=(f(i)-coef1(i,1)*g(i-1))/w(i)
        w(i)=coef1(i,3)/w(i)
      end do
      u(nx-1,:)=g(nx-1)

      do i=nx-2,1,-1
        u(i,j+1)=g(i)-w(i)*u(i+1,j+1)
      end do
    end do

      deallocate(g,w)

    return
end subroutine Solver_Crank-Nicolson

end module mod_fdm_tool_C3

!=============================================================
! cLi, 2022-11-24
!
! Solving example 3.1.1 in p.67:
!
!   -u'(t) - u''(x) = 0,  x \in (0, 1),t \in(0,1]
!   u(x,0) =  exp(x),   x \in (0, 1)
!   u(0,t) = exp(t) ,u(1,t) = exp(1+t) ,t \in(0,1]
!
!
!   
!   
!   u(x,t) = exp(x+t)
!=============================================================
program Crank-Nicolson_2D_DBC

    use iso_fortran_env
    use mod_fdm_tool_C3

implicit none

integer(int32) :: nx, nt, m, k, k_max, fileunit,  i,j
real(real64)   :: hx, xmin, xmax
real(real64)   :: ht, tmin, tmax

real(real64)   :: iter_tol, iter_fac, error = 0.d0
real(real64), allocatable :: x(:), t(:), f(:,:), u(:,:), v(:,:)
character(len=:), allocatable :: dirname, filename, output_fmt
logical :: istat, debug = .false.

namelist /nml_para/ nx, nt, xmin, xmax, tmin, tmax, k_max

print '(1X,A)', 'Example_311'

filename = 'nml_para_DBC.txt'
open(newunit=fileunit, file=filename)
read(unit=fileunit, nml=nml_para)
close(fileunit)

dirname = 'Out'
inquire(file=dirname, exist=istat)
if (.not. istat) then
    call Make_Directory(dirname)
end if


filename = dirname//'/results.txt'
open(newunit=fileunit, file=trim(filename))
write(fileunit, '(a)') 'Inf error at x = 0.5, 0.5,  0.5, 0.5,  0.5,'
write(fileunit, '(a)') '         and t = 0.1, 0.2,  0.3, 0.4,  0.5 are'

do k=1,k_max
   m  = 2**(k-1)
    nx = 10*m
    nt = 10*m
    print '(1X, A, I3, A)', 'The spatial step size of 1/', nx, ',the time step size of 1/',nt,' is running...'

    if (allocated(u)) deallocate(x, t, f, u, v)
    allocate(x(0:nx), t(0:nt), f(0:nx,0:nt), u(0:nx,0:nt), v(0:nx,0:nt))

     ! start the FDM solving procedure
    call Grid_2D_Type (nx, nt, xmin, xmax, tmin, tmax, hx, ht, x, t)

    call Apply_DBC(x, t, u)

    call Gene_Alg_Sys_DBC(ht,hx,nx,coef1,coef2)

    do concurrent (i = 0 : nx)

        v(i,:) = Solu_Exact(x(i),t)
    end do


    call Solver_Crank-Nicolson (coef1,coef2,f,u)

     output_fmt = '(A, I3, 2X, 6(ES9.3, 2X), F5.2)'
    write(fileunit, output_fmt) '1/', nx,'1/', nt,       &

            abs(u(nx/2*1,nt/10*1)-v(nx/2*1,nt/10*1), &

            abs(u(nx/2*1,nt/10*2)-v(u(nx/2*1,nt/10*2), &
            abs(u(nx/2*1,nt/10*3)-v(u(nx/2*1,nt/10*3), &
            abs(u(nx/2*1,nt/10*4)-v(nx/2*1,nt/10*4), &
            abs(u(nx/2*1,nt/10*5)-v(nx/2*1,nt/10*5), &
            maxval(abs(u-v)), error/maxval(abs(u-v))

    error = maxval(abs(u-v))
end do

close(fileunit)
print *, ''
print *, 'The program is done.'

return
end program Crank-Nicolson_2D_DBC