!====================================================================
!  拡散方程式解法ルーチン(C-N法)
!
!  履歴  2000/01/18  小高正嗣
!====================================================================
      program main
!--------------------------------------------------------------------
      implicit real*8 ( a-h,o-z )
!--------------------------------------------------------------------
      parameter ( NZ = 5 )
      parameter ( NZB = NZ + 1 )
!--------------------------------------------------------------------

      dimension T(NZ)
      dimension ZGRID(NZB), FDZ(NZ)

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

      real*4 Tplot(NZB), zplot(nzb)

!--------------------------------------------------------------------

      do iz = 1, NZB
         ZGRID(iz) = 0.1D0 * ( iz - 1 )
      end do

      do iz = 1, NZ
         FDZ(iz) = ZGRID(iz+1) - ZGRID(iz)
      end do

      do iz = 1, NZ
         T(iz) = 2.0D0 * ZGRID(iz+1) 
      end do

      DTI = 1.0D-2

      call SETMTRX
     O   ( ELM_FOR, ELM_BCK, DL, D, DU,
     G     FDZ, 
     D     DTI, NZ, NZB      )

       do iz = 1,NZ
          write(*,*) (ELM_FOR(iz,jz),jz=1,NZ)
       end do

       do iz = 1,NZ
          write(*,*) (ELM_BCK(iz,jz),jz=1,NZ)
       end do

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

!--------------------------------------------------------------------
! 時間積分

      ITIME = 5

      call dclopn

      do iz = 1,NZB
         zplot(iz) = REAL(ZGRID(iz))
      end do

      Tplot(1) = 0.0
      do iz = 2,NZB
         Tplot(iz) = REAL(T(iz-1))
      end do

      call uulinz( nzb, zplot, Tplot, 1, 2 ) !-DCLでお絵書き


      do it = 1, ITIME

         TIME = TIME + DTI

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

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

         do iz = 1, NZ
            T(iz) = B(iz)
         end do

         Tplot(1) = 0.0
         do iz = 2,NZB
            Tplot(iz) = REAL(T(iz-1))
         end do

       	 call uulinz( nzb, zplot, Tplot, 1, 2 )		!-DCLでお絵書き

      end do

      call grcls

      end 
!------------------------------------------------------------------------
      subroutine dclopn

      call sglset( 'LCORNER', .FALSE. )
      write(*,*) 'WORKSTATION ID (I) ?;'
      call sgpwsn
      read(*,*) iws
      call gropn( iws )
      
      call grfrm
      call grswnd( 0.0, 0.5, 0.0, 1.0 )
      call grsvpt( 0.2, 0.8, 0.2, 0.8 )
      call grstrn( 1 )
      call grstrf
      call usdaxs
      call uxsttl( 'B', 'X', 0. )
      call uysttl( 'L', 'T', 0. )

      return
      end

!------------------------------------------------------------------------

