* PACKAGE DINTG  !" 力学 時間積分
*
*"   [HIS] 92/11/26(takepiro)
*
***********************************************************************
      SUBROUTINE TINTGR         !" 時間積分 ( スペクトル )
     M         ( WDTOR , WDPOR , WDDPOR, WDT   , 
     I           WTTOR , WTPOR , WTT   , 
     I           DELT  , ITRUNG, 
     C           FLAPLA, 
     C           ARAD  , DRAD             )
*
*"                   4-th order Runge-Kutta scheme
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [MODIFY] 
      REAL       WDTOR ( NMDIM , 0:KDIM     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , 0:KDIM     )  !" ポロイダル Φ
      REAL       WDDPOR( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
      REAL       WDT   ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
*
*   [INPUT] 
      REAL       WTTOR ( NMDIM, 0:KDIM      )  !" トロイダル Ψ
      REAL       WTPOR ( NMDIM, 0:KDIM      )  !" ポロイダル Φ
      REAL       WTT   ( NMDIM, 0:KDIM      )  !" 温度  Ｔ
*
      INTEGER    ITRUNG                        !" Runge-Kutta 制御変数
      REAL       DELT                          !" 時間刻みΔt
*
      REAL       FLAPLA( NMDIM  )              !" ラプラシアンの係数
*
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )              !" Δｒ(整数)
*
*   [INTERNAL WORK] 
      REAL       WDTOR0( NMDIM , 0:KDIM     )  !" トロイダル Ψ
      REAL       WDDPO0( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
      REAL       WDT0  ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
*"       : 昔の値をセーブ 
*
      REAL       WDTOR1( NMDIM , 0:KDIM     )  !" トロイダル Ψ
      REAL       WDDPO1( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
      REAL       WDT1  ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
*"       : 新しい値をつみ重ねる
*
      REAL       DELTA1, DELTA2, DELTA3
      REAL       WEIGH1, WEIGH2, WEIGH3, WEIGH4
      INTEGER    NM, K
*
*"         < 0. 時間制御 >
*
      IF ( ITRUNG . EQ. 1 ) THEN
         DELTA1  =  DELT * 0.5
         DELTA2  =  DELT * 0.5
         DELTA3  =  DELT
*
         WEIGH1  =       DELT / 6.
         WEIGH2  =  2. * DELT / 6. 
         WEIGH3  =  2. * DELT / 6. 
         WEIGH4  =       DELT / 6. 
*
         CALL COPY ( WDTOR0, WDTOR , NMDIM*KMAX )
         CALL COPY ( WDDPO0, WDDPOR, NMDIM*KMAX )
         CALL COPY ( WDT0,   WDT   , NMDIM*KMAX )
         CALL COPY ( WDTOR1, WDTOR , NMDIM*KMAX )
         CALL COPY ( WDDPO1, WDDPOR, NMDIM*KMAX )
         CALL COPY ( WDT1,   WDT   , NMDIM*KMAX )
      ENDIF
*
*"         < 1. 積分 >
*
*
      IF ( ITRUNG . EQ. 1 )THEN
         DO 1100 K = 0 , KDIM
            DO 1100 NM = 1 , NMDIM
*"    < トロイダル > 
               WDTOR ( NM,K ) = WDTOR0( NM,K )
     &                     + WTTOR ( NM,K )* DELTA1
               WDTOR1( NM,K ) = WDTOR1( NM,K )
     &                     + WTTOR ( NM,K )* WEIGH1
*"    < 温度 >
               WDT ( NM,K ) = WDT0 ( NM,K )
     &                     + WTT ( NM,K )* DELTA1
               WDT1( NM,K ) = WDT1 ( NM,K )
     &                     + WTT ( NM,K )* WEIGH1
 1100    CONTINUE
         DO 1150 K = 0 , KDIM
            DO 1150 NM = 1 , NMDIM
*"    < ポロイダル >
               WDDPOR ( NM,K ) = WDDPO0 ( NM,K )
     &                     + WTPOR ( NM,K )* DELTA1
               WDDPO1 ( NM,K ) = WDDPO1 ( NM,K )
     &                     + WTPOR ( NM,K )* WEIGH1
 1150    CONTINUE 
*
      ELSE IF ( ITRUNG . EQ. 2 )THEN
         DO 1200 K = 0 , KDIM
            DO 1200 NM = 1 , NMDIM
*"    < トロイダル > 
               WDTOR ( NM,K ) = WDTOR0( NM,K )
     &                     + WTTOR ( NM,K )* DELTA2
               WDTOR1( NM,K ) = WDTOR1( NM,K )
     &                     + WTTOR ( NM,K )* WEIGH2
*"    < 温度 >
               WDT ( NM,K ) = WDT0 ( NM,K )
     &                     + WTT ( NM,K )* DELTA2
               WDT1( NM,K ) = WDT1 ( NM,K )
     &                     + WTT ( NM,K )* WEIGH2
 1200    CONTINUE
         DO 1250 K = 0 , KDIM
            DO 1250 NM = 1 , NMDIM
*"    < ポロイダル >
               WDDPOR ( NM,K ) = WDDPO0 ( NM,K )
     &                     + WTPOR ( NM,K )* DELTA2
               WDDPO1 ( NM,K ) = WDDPO1 ( NM,K )
     &                     + WTPOR ( NM,K )* WEIGH2
 1250    CONTINUE 
*
      ELSE IF ( ITRUNG . EQ. 3 )THEN
         DO 1300 K = 0 , KDIM
            DO 1300 NM = 1 , NMDIM
*"    < トロイダル > 
               WDTOR ( NM,K ) = WDTOR0( NM,K )
     &                     + WTTOR ( NM,K )* DELTA3
               WDTOR1( NM,K ) = WDTOR1( NM,K )
     &                     + WTTOR ( NM,K )* WEIGH3
*"    < 温度 >
               WDT ( NM,K ) = WDT0 ( NM,K )
     &                     + WTT ( NM,K )* DELTA3
               WDT1( NM,K ) = WDT1 ( NM,K )
     &                     + WTT ( NM,K )* WEIGH3
 1300    CONTINUE
         DO 1350 K = 0, KDIM
            DO 1350 NM = 1 , NMDIM
*"    < ポロイダル >
               WDDPOR ( NM,K ) = WDDPO0 ( NM,K )
     &                     + WTPOR ( NM,K )* DELTA3
               WDDPO1 ( NM,K ) = WDDPO1 ( NM,K )
     &                     + WTPOR ( NM,K )* WEIGH3
 1350    CONTINUE 
*
      ELSE IF ( ITRUNG . EQ. 4 ) THEN
         DO 1400 K = 0, KDIM
            DO 1400 NM = 1, NMDIM
*"    < トロイダル >
               WDTOR( NM,K ) = WDTOR1( NM,K )
     &                     + WTTOR ( NM,K )* WEIGH4
*"    < 温度 >
               WDT ( NM,K ) = WDT1 ( NM,K )
     &                     + WTT ( NM,K )* WEIGH4
 1400    CONTINUE
         DO 1450 K = 0 , KDIM
            DO 1450 NM = 1 , NMDIM
*"    < ポロイダル > 
               WDDPOR ( NM,K ) = WDDPO1 ( NM,K )
     &                     + WTPOR ( NM,K )* WEIGH4
 1450    CONTINUE 
      ELSE
         CALL MSGDMP ( 'E', 'TINTGR', 
     &                 'THE ORDER OF RUNGE-KUTTA SCHEME INVALID' )
      ENDIF
*
*"         < 2. 境界値 >
*
      CALL BNDRY
     M       ( WDTOR , WDPOR , WDT  , WDDPOR, 
     C         FLAPLA,
     C         ARAD  , DRAD                   )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE MOVGD        !"  格子点値の移動
     O         ( GBU   , GBV   , GBW   , GBT   , 
     O           GBTOR , GBPOR ,
     I           GAU   , GAV   , GAW   , GAT   , 
     I           GATOR , GAPOR                          )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       GBU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ
      REAL       GBV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ
      REAL       GBW   ( IDIM, JDIM, KMAX )  !" 鉛直風 ｗ
      REAL       GBT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
      REAL       GBTOR ( IDIM, JDIM, KMAX )  !" トロイダル Ψ
      REAL       GBPOR ( IDIM, JDIM, KMAX )  !"ポロイダル Φ
*"       : 格子点データ(t-Δt) ->  格子点データ(t)
*
*   [INPUT] 
      REAL       GAU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ
      REAL       GAV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ
      REAL       GAW   ( IDIM, JDIM, KMAX )  !" 鉛直風 ｗ
      REAL       GAT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
      REAL       GATOR ( IDIM, JDIM, KMAX )  !" トロイダル Ψ
      REAL       GAPOR ( IDIM, JDIM, KMAX )  !" ポロイダル Φ
*"       : 格子点データ(t)
*
      CALL COPY ( GBU   , GAU   , IDIM*JDIM*KMAX )
      CALL COPY ( GBV   , GAV   , IDIM*JDIM*KMAX )
      CALL COPY ( GBW   , GAW   , IDIM*JDIM*KMAX )
      CALL COPY ( GBT   , GAT   , IDIM*JDIM*KMAX )
      CALL COPY ( GBTOR , GATOR , IDIM*JDIM*KMAX )
      CALL COPY ( GBPOR , GAPOR , IDIM*JDIM*KMAX )
*
      RETURN
      END
