* PACKAGE DRUNGE    !" P \infty 近似 時間積分( 4th order Runge-Kutta Scheme )
*
*"  [HIS] 96/03/14 (takepiro) 
*
***********************************************************************
      SUBROUTINE DRUNGE
     O             ( PSI3  , ZETA3 , T3  ,
     I               PSI2  , ZETA2 , T2  , 
     I               RAYL  , PRND  , QINT,
     C               DX    , DZ    , DT  ,
     W               PSI   , ZETA  , T   , 
     W               TNDZ  , TNDT  ,
     W               AJAC  , ALAP          )
*
*"      << 4th order Runge-Kutta scheme により時間積分する >>
*"                     96/03/14 (竹広真一)
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [OUTPUT]
      REAL PSI3  ( 0:NX+1 , 0:NZ+1 ) !" 流線関数(現STEPの値)
      REAL ZETA3 ( 0:NX+1 , 0:NZ+1 ) !" 渦度(現STEPの値)
      REAL T3    ( 0:NX+1 , 0:NZ+1 ) !" 温度(現STEPの値)
*
* [INPUT]
      REAL PSI2  ( 0:NX+1 , 0:NZ+1 ) !" 流線関数(1STEP前の値)
      REAL ZETA2 ( 0:NX+1 , 0:NZ+1 ) !" 渦度(1STEP前の値)
      REAL T2    ( 0:NX+1 , 0:NZ+1 ) !" 温度(1STEP前の値)
*
      REAL PRND                 !" プラントル数(ダミー)
      REAL RAYL                 !" レイリー数
      REAL QINT                 !" 加熱率
*
      REAL DX                   !" X 軸上の格子点間隔
      REAL DZ                   !" Z 軸上の格子点間隔
      REAL DT                   !" 時間間隔
*
* [WORK]
      REAL PSI  ( 0:NX+1 , 0:NZ+1 ) !" 流線関数(仮積分用)
      REAL ZETA ( 0:NX+1 , 0:NZ+1 ) !" 渦度(仮積分用)
      REAL T    ( 0:NX+1 , 0:NZ+1 ) !" 温度(仮積分用)
      REAL TNDZ ( 0:NX+1 , 0:NZ+1 ) !" 第 1 段階での渦度の時間変化
      REAL TNDT ( 0:NX+1 , 0:NZ+1 ) !" 第 1 段階での温度の時間変化
      REAL AJAC ( 0:NX+1 , 0:NZ+1 ) !"ヤコビアン
      REAL ALAP ( 0:NX+1 , 0:NZ+1 ) !"ラプラシアン
*
*  [INTERNAL WORK]
      REAL DELT(3)              !" 仮積分用時間間隔
      REAL WEIGH(4)             !" 出力積分用時間間隔
*
      INTEGER IRUNG                !"ループカウンタ
      INTEGER IX                   !"ループカウンタ
      INTEGER IZ                   !"ループカウンタ
*
      LOGICAL OFIRST
      DATA    OFIRST  / .TRUE. /
      SAVE  OFIRST 
*
      IF( OFIRST ) THEN
         OFIRST = .FALSE.
         CALL MSGDMP('M','DRUNGE',
     &               '4TH RUNGE-KUTTA SCHEME : 96/03/14' )
         CALL MSGDMP('M','DRUNGE',
     &               '***** INFINITE PRNDTL NUMBER APPROX. *****' )
      ENDIF
*
*"         < 0. 時間制御 >
*
      DELT(1)  =  DT * 0.5
      DELT(2)  =  DT * 0.5
      DELT(3)  =  DT
*
      WEIGH(1)  =       DT / 6.
      WEIGH(2)  =  2. * DT / 6. 
      WEIGH(3)  =  2. * DT / 6. 
      WEIGH(4)  =       DT / 6. 
*
*"             出力用積分
      CALL COPY ( PSI3  , PSI2  , (NX+2)*(NZ+2) )
      CALL COPY ( ZETA3 , ZETA2 , (NX+2)*(NZ+2) )
      CALL COPY ( T3    , T2    , (NX+2)*(NZ+2) )
*
*"             仮積分
      CALL COPY ( PSI   , PSI2  , (NX+2)*(NZ+2) )
      CALL COPY ( ZETA  , ZETA2 , (NX+2)*(NZ+2) )
      CALL COPY ( T     , T2    , (NX+2)*(NZ+2) )
*
*"         < 1. 積分 >
*
      DO 1000 IRUNG = 1, 4
         CALL DEULTT            !" 温度の式の積分
     O      ( TNDT  , 
     I        PSI   , T   , 
     C        QINT  , 
     C        DX    , DZ  ,
     W        AJAC  , ALAP    )
*
         IF ( IRUNG .LT. 4 ) THEN       !" 仮積分
            DO 1100 IZ = 1, NZ 
               DO 1110 IX = 1 , NX
                  T    ( IX,IZ ) = T2( IX,IZ )
     &                 + TNDT ( IX,IZ )* DELT(IRUNG)
 1110          CONTINUE 
 1100       CONTINUE 
*
            CALL DTMP2Z        !"  温度から渦度の計算(P inf 近似)
     O            ( ZETA  , 
     I              T     , 
     C              RAYL  , 
     C              DX    , DZ    , DT  ,
     W              ALAP   )
*
             CALL DBNDR
     O            ( PSI   , 
     M              ZETA  , T   ,
     I              DX    , DZ    )
          ENDIF                         !" 仮積分ここまで
*
          DO 1200 IZ = 1, NZ            !"   出力用積分
             DO 1210 IX = 1, NX
                T3   ( IX,IZ ) = T3( IX,IZ )
     &               + TNDT ( IX,IZ ) * WEIGH(IRUNG)
 1210        CONTINUE
 1200     CONTINUE
*
 1000 CONTINUE 
*
      CALL DTMP2Z               !"  温度から渦度の計算(P inf 近似)
     O       ( ZETA3 , 
     I         T3    , 
     C         RAYL  , 
     C         DX    , DZ    , DT  ,
     W         ALAP   )
*
*"         < 2. 境界値 >
*
      CALL DBNDR
     O      ( PSI3  , 
     M        ZETA3 , T3  ,
     I        DX    , DZ    )
*
      RETURN
      END