* PACKAGE DINTG  !" 力学 時間積分
*
*"[履歴]  97/05/19 堀之内武 bug fix
***********************************************************************
      SUBROUTINE TINTGR         !" 時間積分 ( スペクトル )
     M         ( WDVOR , WDDIV , WDT   , WDPI  , WDQ   ,
     I           WTVOR , WTDIV , WTT   , WTPI  , WTQ   ,
     I           DMLU  , KPIV  , DELT  ,
     C           DIFV  , DIFT  , FLAPLA, WPHIS ,
     C           GG    , TMTX  , DSIG                   )
*
*"                 Arakawa & Suarez 鉛直差分 σ座標
*"                   semi implicit scheme
*"                   implicit horizontal diffusion
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [MODIFY] 
      REAL       WDVOR ( NMDIM , KMAX     )  !" 渦度  ζ
      REAL       WDDIV ( NMDIM , KMAX     )  !" 発散  Ｄ
      REAL       WDT   ( NMDIM , KMAX     )  !" 温度  Ｔ
      REAL       WDPI  ( NMDIM            )  !" log Ps ＝π
      REAL       WDQ   ( NMDIM , KMAX     )  !" 比湿  ｑ
*
*   [INPUT] 
      REAL       WTVOR ( NMDIM, KMAX      )  !" 渦度  ζ 時間変化
      REAL       WTDIV ( NMDIM, KMAX      )  !" 発散  Ｄ 時間変化
      REAL       WTT   ( NMDIM, KMAX      )  !" 温度  Ｔ 時間変化
      REAL       WTPI  ( NMDIM            )  !" logPs π 時間変化
      REAL       WTQ   ( NMDIM, KMAX      )  !" 比湿  ｑ 時間変化
*
      REAL       DMLU  ( NMDIM, KMAX, KMAX ) !" ＳＩ行列のＬＵ分解
      INTEGER    KPIV  ( NMDIM, KMAX       ) !" ＳＩ行列のピボット
      REAL       DELT                        !" 時間刻みΔt
*
      REAL       DIFV  ( NMDIM  )            !" 運動量水平拡散係数
      REAL       DIFT  ( NMDIM  )            !" 熱，水水平拡散係数
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       WPHIS ( NMDIM )             !" 地表ジオポテンシャル
*
      REAL       GG    ( KMAX )              !" セミインプリシット用
      REAL       TMTX  ( KMAX, KMAX )        !" セミインプリシット用
      REAL       DSIG  ( KMAX )              !" Δσ
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           WDTX  , WDPIX ,
     &           WDPHIX, WBRDIV
      REAL       WDTX     ( NMDIM, KMAX )    !" 温度(仮積分)
      REAL       WDPIX    ( NMDIM )          !" π  (仮積分)
      REAL       WDPHIX   ( NMDIM, KMAX )    !" Φ  (仮積分)
      REAL       WBRDIV   ( NMDIM, KMAX )    !" 時間平均のＤ
*
      REAL       DELT2
      INTEGER    NM, K, KK
*
      DELT2  = DELT*2.  
*
*"         < 1. 非重力波項の仮積分 >
*
      DO 1100 K = 1 , KMAX
         DO 1100 NM=1,NMDIM !"97/05/19堀之内：下の行で DELT2*.. -> DELT*..
            WDTX  ( NM,K ) = WDT  ( NM,K )*( 1. - DELT*DIFT( NM ) )
     &                     + WTT  ( NM,K )* DELT
 1100 CONTINUE
*
      DO 1200 NM = 1 , NMDIM
            WDPIX ( NM )   = WDPI ( NM ) + WTPI  ( NM )* DELT
 1200 CONTINUE
*
*"         < 2. ジオポテンシャル >
*
      CALL WHYDRO
     O         ( WDPHIX,
     I           WDTX   )
*
*"         < 3. 発散方程式の右辺  >
*
      DO 3100 K = 1 , KMAX
*
         DO 3110 NM = 1 , NMDIM
*
            WBRDIV( NM,K )
     &       =  WDDIV( NM,K ) * ( 1. + DELT *DIFV( NM ) )
     &                        * ( 1. + DELT2*DIFT( NM ) )
     &           + DELT*(   WTDIV ( NM,K )
     &                    - FLAPLA( NM )
     &                      * (  WDPHIX( NM,K )
     &                           + ( 1. + DELT2*DIFT( NM ) )
     &                             * (  WPHIS ( NM )
     &                                + WDPIX ( NM )*GG( K ) )  )   )
*
 3110    CONTINUE
 3100 CONTINUE
*
*"         < 4. 時間平均のＤをLU行列で解く >
*
      CALL LUSOLV
     M         ( WBRDIV,
     I           DMLU  , KPIV  ,
     D           1     , NMDIM , KMAX   )
*
*"         < 5. 熱力, 地表気圧の時間変化項 >
*
      DO 5100 K = 1, KMAX
         DO 5100 KK = 1 , KMAX
            DO 5100 NM = 1, NMDIM
               WTT ( NM,K ) = WTT   ( NM,K )
     &                       - TMTX ( K,KK ) * WBRDIV( NM,KK )
 5100 CONTINUE
*
      DO 5200 KK = 1 , KMAX
         DO 5200 NM = 1, NMDIM
              WTPI ( NM   ) = WTPI  ( NM   )
     &                        - DSIG(   KK ) * WBRDIV( NM,KK )
 5200 CONTINUE
*
*"         < 7. 時間積分  拡散は前方差分 >
*
      DO 7100  K = 1 , KMAX
         DO 7100  NM = 1 , NMDIM
*
*"        ------------------  < 渦度 > ---------------------------
*
            WDVOR( NM,K ) = ( WDVOR( NM,K ) + WTVOR( NM,K )*DELT2 )
     &                     /( 1. + DELT2*DIFV( NM )               )
*
*"        ------------------  < 発散 >  -------------------------
*
            WDDIV ( NM,K ) =  WBRDIV( NM,K )*2. - WDDIV ( NM,K )
*
*"        ------------------  < 温度 >  -------------------------
*
            WDT  ( NM,K ) = ( WDT  ( NM,K ) + WTT  ( NM,K )*DELT2 )
     &                    / ( 1. + DELT2*DIFT( NM )               )
*
*"        ------------------  < 比湿 >  -------------------------
*
            WDQ  ( NM,K ) = ( WDQ  ( NM,K ) + WTQ  ( NM,K )*DELT2 )
     &                    / ( 1. + DELT2*DIFT( NM )               )
*
 7100 CONTINUE
*
      DO 7200 NM = 1 , NMDIM
*
*"        ------------------  < ln(地表気圧) >  --------------------
*
            WDPI ( NM   ) =   WDPI ( NM   ) + WTPI ( NM   )*DELT2
*
 7200 CONTINUE
*
      RETURN
      END
***********************************************************************
      SUBROUTINE MOVGD        !"  格子点値の移動
     O         ( GBU   , GBV   , GBT   , GBPS  , GBQ   ,
     O           GBVOR , GBDIV ,
     I           GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     I           GAVOR , GADIV                          )
*

#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       GBU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ
      REAL       GBV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ
      REAL       GBT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
      REAL       GBPS  ( IDIM, JDIM       )  !" 地表気圧
      REAL       GBQ   ( IDIM, JDIM, KMAX )  !" 比湿  ｑ
      REAL       GBVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ
      REAL       GBDIV ( IDIM, JDIM, KMAX )  !" 発散  Ｄ
*"       : 格子点データ(t-Δt) ->  格子点データ(t)
*
*   [INPUT] 
      REAL       GAU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ
      REAL       GAV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ
      REAL       GAT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
      REAL       GAPS  ( IDIM, JDIM       )  !" 地表気圧
      REAL       GAQ   ( IDIM, JDIM, KMAX )  !" 比湿  ｑ
      REAL       GAVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ
      REAL       GADIV ( IDIM, JDIM, KMAX )  !" 発散  Ｄ
*"       : 格子点データ(t)
*
      CALL COPY ( GBU   , GAU   , IDIM*JDIM*KMAX )
      CALL COPY ( GBV   , GAV   , IDIM*JDIM*KMAX )
      CALL COPY ( GBT   , GAT   , IDIM*JDIM*KMAX )
      CALL COPY ( GBPS  , GAPS  , IDIM*JDIM      )
      CALL COPY ( GBQ   , GAQ   , IDIM*JDIM*KMAX )
      CALL COPY ( GBVOR , GAVOR , IDIM*JDIM*KMAX )
      CALL COPY ( GBDIV , GADIV , IDIM*JDIM*KMAX )
*
      RETURN
      END
