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