* PACKAGE DMTRX 
*
**********************************************************************
      SUBROUTINE DMTRIX        
     O         ( DMLU  , KPIV  ,
     I           DELT  ,
     C           WHMTX , GCMTX ,
     C           NMO   , ER    , DIFV  , DIFT  )
*
*"              implicit diffusion
*




* PACKAGE  ZXT10L32 !" 格子点数，波数 X-Z
*
      INTEGER     IMAX           !" 経度格子数
      INTEGER     JMAX           !" 緯度格子数
      INTEGER     IDIM           !" 経度配列サイズ
      INTEGER     JDIM           !" 緯度配列サイズ
      INTEGER     MMAX           !" 最大東西波数
      INTEGER     LMAX           !" 最大南北波数
      INTEGER     NMAX           !" 最大全波数
      INTEGER     MINT           !" 東西波数間隔
      INTEGER     KMAX           !" 鉛直層数
      INTEGER     KDIM           !" 鉛直配列サイズ
      INTEGER     NMDIM          !" 水平波数成分数
      INTEGER     JMXHF          !" JMAX/2+1
      INTEGER     NWORK          !" ワーク領域大きさ
*
      PARAMETER ( IMAX  = 32,     JMAX =  1,
     &            IDIM  = IMAX+1, JDIM = JMAX,
     &            MMAX  = 10,     LMAX =  0,    NMAX = MMAX, MINT = 1,
     &            KMAX  = 32,     KDIM = KMAX,
     &            NMDIM = (MMAX/MINT+1)*(2*(NMAX+1)-MMAX)
     &                  - (NMAX-LMAX)/MINT*(NMAX-LMAX+1),
     &            JMXHF = JMAX/2+1,
     &            NWORK = IDIM*JDIM*KMAX*12                       )

*
*   [OUTPUT] 
      ReaL*8       DMLU  ( NMDIM, KMAX, KMAX ) 
      INTEGER    KPIV  ( NMDIM, KMAX       ) 
*
*   [INPUT] 
      ReaL*8       DELT                        
      ReaL*8       WHMTX ( KMAX, KMAX )        
      ReaL*8       GCMTX ( KMAX, KMAX )        
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) 
      ReaL*8       ER                          
      ReaL*8       DIFV  ( NMDIM  )            
      ReaL*8       DIFT  ( NMDIM  )            
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           DM
      ReaL*8       DM    ( KMAX, KMAX )        
      INTEGER    KMAXD
      PARAMETER  (KMAXD=100)
      INTEGER    KP    ( KMAXD )             
*
      INTEGER    MXNM, L, M, N, NM, K, KK
      ReaL*8       FLAPLA
*
      IF  ( KMAX .GT. KMAXD ) THEN
         CALL MSGDMP( 'E','DMTRIX','INSUFFICIENT WORK STRAGE.')
      ENDIF
*
      MXNM = 0
*
      DO 3300 N = 0, NMAX
*"                                                 ↓=== LOOP N ===
*
*"         < 1. semi-implicit の行列 >
*
         FLAPLA = - DBLE( N ) * DBLE( N+1 ) / ER**2
*
         DO 1000 M = 0, MIN(N,MMAX), MINT
            IF ( N-M .LE. LMAX ) THEN
               NM = NMO( 1, M, N-M )
               GOTO 1010
            ENDIF
 1000    CONTINUE
*"                   no n(m,l) ----->  continue
         WRITE (6,*) ' WARNING [DMTRX]: N=', N, ' DOES NOT EXIST '
         GOTO 3300
*
 1010    CONTINUE
*
         DO 1100 K = 1, KMAX
            DO 1100 KK = 1, KMAX
               DM ( K,KK ) = - DELT**2 * FLAPLA
     &                         * (   WHMTX( K,KK )
     &                             + GCMTX( K,KK )
     &                               * ( 1. + DELT*2.*DIFT(NM) )  )
               IF ( K .EQ. KK) THEN
                  DM ( K,KK ) = DM ( K,KK )
     &                        +   ( 1. + DELT*2.* DIFV(NM) )
     &                          * ( 1. + DELT*2.* DIFT(NM) )
               ENDIF
 1100    CONTINUE
*
*"         < 2. LU行列 >
*
         CALL LUMAKE
     M         ( DM    ,
     O           KP    ,
     D           1     , KMAX   )
*
*"         < 3. 配列の詰め替え >
*
         DO 3200 M = 0 , MMAX, MINT
            L = N - M
*
            IF (   ( L .GE. 0 ) .AND. ( L .LE. LMAX )   )THEN
               DO 3100 K = 1, KMAX
*
                  DO 3110 KK = 1, KMAX
                     DMLU ( NMO(1,M,L),K,KK ) = DM ( K,KK )
                     DMLU ( NMO(2,M,L),K,KK ) = DM ( K,KK )
 3110             CONTINUE
                  KPIV ( NMO(1,M,L),K ) = KP ( K )
                  KPIV ( NMO(2,M,L),K ) = KP ( K )
*
                  MXNM = MAX( MXNM, NMO(1,M,L) )
                  MXNM = MAX( MXNM, NMO(2,M,L) )
 3100          CONTINUE
            ENDIF
*
 3200    CONTINUE
*
*"                                                 ↑=== LOOP N ===
 3300 CONTINUE
*
      DO 4100 NM = MXNM+1, NMDIM
         DO 4100 K = 1, KMAX
             DO 4110 KK = 1, KMAX
                DMLU ( NM,K,KK ) = DM ( K,KK )
 4110        CONTINUE
             KPIV ( NM,K ) = KP ( K )
 4100 CONTINUE
            
      RETURN
      END
