!====================================================================
!  拡散方程式解法ルーチン(C-N法)
!  係数行列の計算
!
!  履歴  2000/01/07  小高正嗣
!====================================================================
      subroutine SETMTRX
     I   ( CDTURB_M, DENS0,
     O     ELM_FOR, ELM_BCK, DL, D, DU,
     G     DX, FDZ, FDZM,
     D     DTI, NZ, NZB      )
!--------------------------------------------------------------------
      implicit real*8 ( a-h,o-z )
!--------------------------------------------------------------------

      dimension CDTURB_M(-2:NZB), DENS0(-2:NZB)
      dimension 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)

      dimension DENS0_M(-2:NZB)

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

      do iz = 0, NZ+1
         DENS0_M(iz) = ( DENS0(iz) + DENS0(iz-1) ) * 0.5D0
      end do

      do jz = -2, NZB
         do iz = -2, NZB
            ELM_FOR(iz,jz) = 0.0D0
            ELM_BCK(iz,jz) = 0.0D0
         end do
      end do

      ELM_FOR(0,0) = 1.0D0 + DTI / FDZ(0) / DENS0(0) / 2.0D0 *
     \     ( DENS0_M(1) * CDTURB_M(1) / FDZM(1) +  
     \       DENS0_M(0) * CDTURB_M(0) / FDZM(0) * 2.0D0 ) / DX ** 2

      ELM_BCK(0,0) = 1.0D0 - DTI / FDZ(0) / DENS0(0) / 2.0D0 *
     \     ( DENS0_M(1) * CDTURB_M(1) / FDZM(1) + 
     \       DENS0_M(0) * CDTURB_M(0) / FDZM(0) * 2.0D0 ) / DX ** 2

      do iz = 1, NZ-2
         ELM_FOR(iz,iz) = 1.0D0 + DTI / FDZ(iz) / DENS0(iz) / 2.0D0 *
     \        ( DENS0_M(iz+1) * CDTURB_M(iz+1) / FDZM(iz+1) + 
     \          DENS0_M(iz) * CDTURB_M(iz) / FDZM(iz) ) / DX ** 2

         ELM_BCK(iz,iz) = 1.0D0 - DTI / FDZ(iz) / DENS0(iz) / 2.0D0 *
     \        ( DENS0_M(iz+1) * CDTURB_M(iz+1) / FDZM(iz+1) +  
     \          DENS0_M(iz) * CDTURB_M(iz) / FDZM(iz) ) / DX ** 2
      end do

      ELM_FOR(NZ-1,NZ-1) = 1.0D0 
     \     + DTI / FDZ(NZ-1) / DENS0(NZ-1) / 2.0D0 *
     \     ( DENS0_M(NZ) * CDTURB_M(NZ) / FDZM(NZ) * 2.0D0 +  
     \       DENS0_M(NZ-1) * CDTURB_M(NZ-1) / FDZM(NZ-1) ) / DX ** 2

      ELM_BCK(NZ-1,NZ-1) = 1.0D0 
     \     - DTI / FDZ(NZ-1) / DENS0(NZ-1) / 2.0D0 *
     \     ( DENS0_M(NZ) * CDTURB_M(NZ) / FDZM(NZ) * 2.0D0 +  
     \       DENS0_M(NZ-1) * CDTURB_M(NZ-1) / FDZM(NZ-1) ) / DX ** 2

      do iz = 0, NZ-2
         ELM_FOR(iz,iz+1) = - DTI / FDZ(iz) / DENS0(iz) / 2.0D0 * 
     \        ( DENS0_M(iz+1) * CDTURB_M(iz+1) / FDZM(iz+1) ) / DX ** 2
         ELM_BCK(iz,iz+1) = - ELM_FOR(iz,iz+1) 
      end do

      do iz = 1, NZ-1
         ELM_FOR(iz,iz-1) = - DTI / FDZ(iz) / DENS0(iz) / 2.0D0 * 
     \        ( DENS0_M(iz) * CDTURB_M(iz) / FDZM(iz) ) / DX ** 2
         ELM_BCK(iz,iz-1) = - ELM_FOR(iz,iz-1) 
      end do

!--------------------------------------------------------------------
      do iz = 1, NZ
         D(iz)  = ELM_FOR(iz-1,iz-1)
      end do
         
      do iz = 1, NZ-1
         DU(iz) = ELM_FOR(iz-1,iz)
         DL(iz) = ELM_FOR(iz,iz-1)
      end do
!--------------------------------------------------------------------

      return
      end
