!====================================================================
!  拡散方程式解法ルーチン(C-N法)
!
!  履歴  1999/12/03  小高正嗣
!====================================================================
      program main
!--------------------------------------------------------------------
#include "grid_size_M.f"
!--------------------------------------------------------------------
      implicit real*8 ( a-h,o-z )
!--------------------------------------------------------------------

      dimension TEMP(-2:NZB), CDTURB(-2:NZB), CDTURB_M(-2:NZB)
      dimension ZGRID(-2:NZB), FDZ(-2:NZB), FDZM(-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

      real*4 tplot(nz),zplot(nz)

!--------------------------------------------------------------------
! 格子間隔の設定

      call SETGRD
     O   ( ZGRID, FDZ, FDZM, DX,
     D     NZ, NZB                )

      do iz = 1, nz
         zplot(iz) = REAL(ZGRID(iz-1))
      end do

!--------------------------------------------------------------------
! 初期値の設定

      call INIT
     I   ( ZGRID,
     O     TEMP, CDTURB, 
     D     NZ, NZB )

!--------------------------------------------------------------------
! 係数行列の設定
     
      TIME = 0.0D0
      DTI  = 1.0D2

      do iz = 0, NZ
         CDTURB_M(iz) = ( CDTURB(iz) + CDTURB(iz-1) ) / 2.0D0
      end do

      write(*,*) (CDTURB(iz),iz=0,NZ)
      write(*,*) (CDTURB_M(iz),iz=0,NZ)

      call BOUND( CDTURB_M, NZ, NZB )

      call SETMTRX
     I   ( CDTURB_M,
     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

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

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

      ITIME = 20

      call dclopn

      TEMP(-2)  =0.0D0
      TEMP(-1)  =0.0D0
      TEMP(NZ)  =0.0D0
      TEMP(NZ+1)=0.0D0
      TEMP(NZ+2)=0.0D0

      do iz = 1, NZ
         tplot(iz) = REAL(TEMP(iz-1))
      end do

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

      do it = 1, ITIME
         
         TIME = TIME + DTI

         do iz = 0, NZ-1
            B(iz) = 0.0D0
            do jz = 0, NZ-1
               B(iz) = B(iz) + ELM_BCK(iz, jz) * TEMP(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
            TEMP(iz) = B(iz)
         end do

         call BOUND( TEMP, NZ, NZB )

         TEMP(-2)  =0.0D0
         TEMP(-1)  =0.0D0
         TEMP(NZ)  =0.0D0
         TEMP(NZ+1)=0.0D0
         TEMP(NZ+2)=0.0D0

         do iz = 1, NZ
            tplot(iz) = REAL(TEMP(iz-1))
         end do

       	 call uulinz( nz, 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, 1.0, 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

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



