!====================================================================
!  拡散方程式解法ルーチン(C-N法)
!  係数行列の計算
!
!  履歴  2000/01/07  小高正嗣
!  履歴  2000/09/17  小高正嗣
!====================================================================
      subroutine SETMTRXG
     I   ( CDTURBG,
     O     ELMG_FOR, ELMG_BCK, DGL, DG, DGU,
     G     FDZ_G, FDZM_G,
     D     DTI, NZG, NZGB      )
!--------------------------------------------------------------------
      implicit real*8 ( a-h,o-z )
!--------------------------------------------------------------------

      dimension FDZ_G(0:NZGB), FDZM_G(0:NZG)
      dimension ELMG_FOR(NZG, NZG), ELMG_BCK(NZG,NZG) 
      dimension DGL(NZG-1), DG(NZG), DGU(NZG-1)

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

      do jz = 1, NZG
         do iz = 1, NZG
            ELMG_FOR(iz,jz) = 0.0D0
            ELMG_BCK(iz,jz) = 0.0D0
         end do
      end do

      ELMG_FOR(1,1) = 1.0D0 + 
     \     CDTURBG * DTI / FDZ_G(1) / FDZM_G(1) * 0.25D0
      ELMG_BCK(1,1) = 1.0D0 - 
     \     CDTURBG * DTI / FDZ_G(1) / FDZM_G(1) * 0.25D0

      do iz = 2, NZG-1
         ELMG_FOR(iz,iz) = 1.0D0 + CDTURBG * DTI / FDZ_G(iz) *
     \        ( 1.0D0 / FDZM_G(iz) + 1.0D0  / FDZM_G(iz-1) ) * 0.25D0
         ELMG_BCK(iz,iz) = 1.0D0 - CDTURBG * DTI / FDZ_G(iz) *
     \        ( 1.0D0 / FDZM_G(iz) + 1.0D0  / FDZM_G(iz-1) ) * 0.25D0
      end do

      ELMG_FOR(NZG,NZG) = 1.0D0 + 
     \     CDTURBG * DTI / FDZ_G(NZG) / FDZM_G(NZG-1) * 0.25D0
      ELMG_BCK(NZG,NZG) = 1.0D0 - 
     \     CDTURBG * DTI / FDZ_G(NZG) / FDZM_G(NZG-1) * 0.25D0

      do iz = 1, NZG-1
         ELMG_FOR(iz,iz+1) = 
     \        - CDTURBG * DTI / FDZ_G(iz) / FDZM_G(iz) * 0.25D0
         ELMG_BCK(iz,iz+1) = - ELMG_FOR(iz,iz+1) 
      end do

      do iz = 1, NZG-1
         ELMG_FOR(iz+1,iz) = 
     \        - CDTURBG * DTI / FDZ_G(iz+1) / FDZM_G(iz) * 0.25D0
         ELMG_BCK(iz+1,iz) = - ELMG_FOR(iz+1,iz) 
      end do

!--------------------------------------------------------------------
      do iz = 1, NZG
         DG(iz)  = ELMG_FOR(iz,iz)
      end do
         
      do iz = 1, NZG-1
         DGU(iz) = ELMG_FOR(iz,iz+1)
         DGL(iz) = ELMG_FOR(iz+1,iz)
      end do
!--------------------------------------------------------------------

      return
      end
