* PACKAGE DMTRX !" 力学 時間積分用行列 GCM5(DMTRX)
*
**********************************************************************
      SUBROUTINE DMTRIX        !" 時間積分 ( スペクトル ) 用のLU行列
     O         ( DMLU  , KPIV  ,
     I           DELT  ,
     C           WHMTX , GCMTX ,
     C           NMO   , ER    , DIFV  , DIFT  )
*
*"              implicit diffusion
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       DMLU  ( NMDIM, KMAX, KMAX ) !" ＳＩ行列のＬＵ分解
      INTEGER    KPIV  ( NMDIM, KMAX       ) !" ＳＩ行列のピボット
*
*   [INPUT] 
      REAL       DELT                        !" 時間刻みΔt
      REAL       WHMTX ( KMAX, KMAX )        !" セミインプリシット用
      REAL       GCMTX ( KMAX, KMAX )        !" セミインプリシット用
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       ER                          !" 地球半径
      REAL       DIFV  ( NMDIM  )            !" 運動量水平拡散係数
      REAL       DIFT  ( NMDIM  )            !" 熱，水水平拡散係数
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           DM
      REAL       DM    ( KMAX, KMAX )        !" ワーク
      INTEGER    KMAXD
      PARAMETER  (KMAXD=100)
      INTEGER    KP    ( KMAXD )             !" ピボット
*
      INTEGER    MXNM, L, M, N, NM, K, KK
      REAL       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 = - Real( N ) * Real( 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
