* PACKAGE DHEUN    !" P \infty 近似 時間積分( Heun Scheme )
*
*"  [HIS] 96/03/16 (takepiro) 
*
***********************************************************************
      SUBROUTINE DHEUN
     O             ( PSI3  , ZETA3 , T3  , 
     I               PSI2  , ZETA2 , T2  , 
     I               RAYL  , PRND  , QINT, 
     C               DX    , DZ    , DT  ,
     W               TNDZ  , TNDT  ,
     W               AJAC  , ALAP   ) 
*
*"    << Heun schemeにより次ステップの物理量を計算する (P \infty 近似) >>
*" 
* [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 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]
      INTEGER IX                   !"ループカウンタ
      INTEGER IZ                   !"ループカウンタ
*
      LOGICAL OFIRST
      DATA    OFIRST  / .TRUE. /
*    
      SAVE  OFIRST 
*
      IF( OFIRST ) THEN
         OFIRST = .FALSE.
         CALL MSGDMP('M','DHEUN',
     &               'HEUN SCHEME TIME INTEGRATE : 96/03/16' )
         CALL MSGDMP('M','DHEUN',
     &               '***** INFINITE PRNDTL NUMBER APPROX. *****' )
      ENDIF
*
*"        < 1. Euler scheme : dummyの値を計算 >
*
      CALL DEULTT
     O      ( TNDT  , 
     I        PSI2  , T2  , 
     C        QINT  , 
     C        DX    , DZ    , 
     W        AJAC  , ALAP    )
*
      DO 1040 IX = 1 , NX
        DO 1030 IZ = 1 , NZ
              T3(IX,IZ)
     &    =   T2(IX,IZ) + DT * TNDT( IX, IZ )
 1030   CONTINUE
 1040 CONTINUE
*
      CALL DTMP2Z               !"  温度から渦度の計算(P inf 近似)
     O       ( ZETA3 , 
     I         T3    , 
     C         RAYL  , 
     C         DX    , DZ    , DT  ,
     W         ALAP    )
*
      CALL DBNDR                !" 境界条件の適用と流線関数の計算
     O      ( PSI3  , 
     M        ZETA3 , T3   ,
     I        DX    , DZ    )
*
*"        < 2. 反復する >
*
      CALL DHEUZT
     M       ( T3    , 
     I         PSI3  , 
     I         T2    , TNDT  , 
     C         QINT  , 
     C         DX    , DZ    , DT  ,
     W         AJAC  , ALAP    )
*
      CALL DTMP2Z               !"  温度から渦度の計算(P inf 近似)
     O       ( ZETA3 , 
     I         T3    , 
     C         RAYL  , 
     C         DX    , DZ    , DT  ,
     W         ALAP    )
*
      CALL DBNDR                !" 境界条件の適用と流線関数の計算
     O      ( PSI3  , 
     M        ZETA3 , T3   ,
     I        DX    , DZ    )
*
      RETURN
      END
*****************************************************************
      SUBROUTINE DHEUZT        !" Heun scheme による温度(第 2 段階)
     M             ( T3    , 
     I               PSI3  , 
     I               T2    , TNDT  , 
     C               QINT  , 
     C               DX    , DZ    , DT  ,
     W               AJAC3 , ALAP3   )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [OUTPUT]
      REAL T3 ( 0:NX+1 , 0:NZ+1 )    !" 温度(現STEPの値, 予想値)
*
* [INPUT]
      REAL PSI3  ( 0:NX+1 , 0:NZ+1 ) !" 流線関数(現STEPの値)
      REAL T2    ( 0:NX+1 , 0:NZ+1 ) !" 温度(1STEP前の値)
      REAL TNDT  ( 0:NX+1 , 0:NZ+1 ) !" 第 1 段階での温度の時間変化
*
      REAL QINT                      !" 加熱率
      REAL DX                        !" X 軸上の格子点間隔
      REAL DZ                        !" Z 軸上の格子点間隔
      REAL DT                        !" 時間間隔
*
* [WORK]
      REAL AJAC3 ( 0:NX+1 , 0:NZ+1 ) !" ヤコビアン
      REAL ALAP3 ( 0:NX+1 , 0:NZ+1 ) !" ラプラシアン
*
*  [INTERNAL WORK]
      INTEGER IX                     !" ループカウンタ
      INTEGER IZ                     !" ループカウンタ
*
*"        < 1. ヤコビアンの計算 >
*
      CALL DARJAC
     O       ( AJAC3,
     I         PSI3 , T3  ,
     C         DX   , DZ   )
*
*"        < 2. ラプラシアンの計算 >
*
      CALL DLAPLA
     O       ( ALAP3,
     I         T3   ,
     C         DX   , DZ   )
*
*"        < 3. 各グリッドの温度の値の計算 >
*
      DO 3020 IX = 1 , NX
        DO 3010 IZ = 1 , NZ
             T3( IX,IZ ) 
     &    =  T2( IX,IZ ) 
     &     + DT*(   0.5 * TNDT( IX,IZ )  
     &            + 0.5 * ( AJAC3(IX,IZ) + ALAP3(IX,IZ) + QINT ) )
 3010   CONTINUE
 3020 CONTINUE
*
      RETURN
      END