$DEBUG

PROGRAM PT_Truyen_Nhiet
! DU/Dt = K.D^U/DZ^2
! K : He so truyen nhiet
IMPLICIT NONE
INTEGER N, NT
REAL, ALLOCATABLE :: A(:), B(:), C(:),&
                     G(:), U(:,:), UX(:) 
REAL, PARAMETER :: HSO = 3.E-1
INTEGER I, J
REAL H, K, R, T, L, X, PI
PI = 4.*ATAN(1.)
L = 1.   ! Do dai thanh
K = 120 ! Buoc thoi gian (s)
H = 0.1  ! Khoang chia theo x 
R = HSO*K / H ** 2
N = NINT(L/H) ! So khoang chia theo x
NT = 24*3600/K      ! So buoc tich phan
ALLOCATE(A(2:N-1), B(1:N-1), C(1:N-2), &
         G(1:N-1), U(0:N,0:NT), UX(1:N-1))
A = -R        ! Xac dinh 3 duong cheo cua A
B = 2 + 2 * R
C = -R
DO I = 0, N   ! Dieu kien ban dau
  X = I*H
  U(I,0) = 3 * X + 18.5
  IF (X > 0.5) U(I,0) = 20.      
END DO
DO I = 0, NT								! Dieu kien bien
   U(0,I) = 3*cos(2*PI/(24.*3600.)*I*K + 2*PI/3) + 20. 
END DO
U(N,:) = 20.
T = 0
PRINT "(11F7.4)", (I * H, I = 1, N-1)

DO J = 1, NT  ! Thoi diem tich phan
  G = R * (U(0:N-2,0) + U(2:N,0)) + &
     (2 - 2 * R) * U(1:N-1,0)			! Ve phai
  G(1) = G(1) + R * U(0,J)
  G(N-1) = G(N-1) + R * U(N,J)
    
  CALL TriDiag( A, B, C, UX, G )
  U(1:N-1,J) = UX
  U(1:N-1,0) = UX
END DO

DO J=1, NT
   PRINT "(11F7.4)", U(1:N,J)
ENDDO

OPEN (3,FILE='TEST.DAT')
DO I = 0, N
   X = I*H*10
   DO J = 0, NT
      IF (i==0) U(0,J) = 3*cos(2*PI/(24.*3600.)*J*K + 2*PI/3) + 20.
	  T = J*K/3600.
	  WRITE(3, '(3F10.5)') X, T, U(I,J) 
   END DO
END DO

CONTAINS
SUBROUTINE TriDiag( A, B, C, X, G )
IMPLICIT NONE
  REAL B(:)                  ! main diagonal 
  REAL A(2:)                 ! lower diagonal
  REAL C(:)                  ! upper diagonal
  REAL, INTENT(OUT) :: X(:)  ! unknown
  REAL G(:)                  ! RHS
  REAL W( SIZE(B) )          ! working space
  REAL T
  INTEGER I, J, N
  N = SIZE(B)
  W = B
  DO I = 2, N
    T = A(I) / W(I-1)
    W(I) = W(I) - C(I-1) * T
    G(I) = G(I) - G(I-1) * T
  END DO
X(N) = G(N) / W(N)
  DO I = 1, N-1
    J = N-I
    X(J) = (G(J) - C(J) * X(J+1)) / W(J)  
  END DO
END SUBROUTINE TriDiag
				
END