!====================================================================
!  拡散方程式解法ルーチン(C-N法)
!
!  履歴  1999/12/03  小高正嗣
!        2000/02/14  小高正嗣
!====================================================================
      subroutine CLDIFEQ
     I      ( TPOT0, CDTURB_M,
     B        DENS0, 
     C        DTI,
     G        DX   , FDZ, FDZM,
     D        NZ   , NZB     )

!--------------------------------------------------------------------
      implicit real*8 ( a-h,o-z )
!--------------------------------------------------------------------

      dimension TPOT0(-2:NZB), CDTURB_M(-2:NZB), DENS0(-2:NZB)
      dimension FDZ(-2:NZB), FDZM(-2:NZB)

      dimension DENS0_M(-2:NZB)

      dimension ELM_FOR(-2:NZB, -2:NZB), ELM_BCK(-2:NZB, -2:NZB) 

      dimension DL(NZ-1), D(NZ), DU(NZ-1), DU2(NZ-2), B(0:NZ-1)
      dimension IPIV(NZ) 
      character(1) TRANS

!--------------------------------------------------------------------
! 係数行列の設定

      call SETMTRX
     I   ( CDTURB_M, DENS0,
     O     ELM_FOR, ELM_BCK, DL, D, DU,
     G     DX, FDZ, FDZM,
     D     DTI, NZ, NZB      )

!       do iz = 0,9
!          write(*,*) (ELM_FOR(iz,jz),jz=0,9)
!       end do

!       do iz = 0,9
!          write(*,*) (ELM_BCK(iz,jz),jz=0,9)
!       end do

!--------------------------------------------------------------------
! AX=B における係数行列 A の分解(LAPACK)

      call DGTTRF
     \   ( NZ, DL, D, DU, DU2, IPIV, INFO)

!--------------------------------------------------------------------
! AX=B を解いて X を解く(LAPACK)

      do iz = 0, NZ-1
         B(iz) = 0.0D0
         do jz = 0, NZ-1
            B(iz) = B(iz) + ELM_BCK(iz, jz) * TPOT0(jz)
         end do
      end do

      TRANS = 'N'
      call DGTTRS
     \     ( TRANS, NZ, 1, DL, D, DU, DU2, IPIV, B, NZ, INFO)

      do iz = 0, NZ-1
         TPOT0(iz) = B(iz)
      end do

      return
      end 



