* PACKAGE DEULER    !" 時間積分( Euler Scheme )
*
*"  [HIS] 92/01/13 (momoko)
*"        95/08/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 により次ステップの物理量を計算する >>
*"                     92/01/13 (石渡正樹)
*
* [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 : 95/08/14' )
      ENDIF
*
*"        < 1. 内部領域の渦度の計算 >
*
      CALL DEULRZ
     O       ( ZETA3 , 
     I         PSI2  , ZETA2 , T2  , 
     C         RAYL  , PRND  ,
     C         DX    , DZ    , DT  ,
     W         AJAC2 , ALAP2   )
*
*"        < 2. 内部領域の温度の計算 >
*
      CALL DEULRT
     O       ( T3    , 
     I         PSI2  , T2  , 
     C         QINT  , 
     C         DX    , DZ    , DT  ,
     W         AJAC2 , ALAP2   )
*
*"        < 3. 境界条件の適用と流線関数の計算 >
*
      CALL DBNDR
     O       ( PSI3 , 
     M         ZETA3, T3  ,
     I         DX   , DZ   )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE DEULRZ          !" Euler scheme による渦度
     O             ( ZETA3 , 
     I               PSI2  , ZETA2 , T2  , 
     C               RAYL  , PRND  ,
     C               DX    , DZ    , DT  ,
     W               AJAC2 , ALAP2   )
*
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [OUTPUT]
      REAL ZETA3 ( 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 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 DEULTZ
     O       ( ZETA3 , 
     I         PSI2  , ZETA2 , T2  , 
     C         RAYL  , PRND  ,
     C         DX    , DZ    , 
     W         AJAC2 , ALAP2   )
*
*"        < 2. 各グリッドの渦度の値の計算 >
*
      DO 3020 IZ = 1 , NZ
        DO 3010 IX = 1 , NX
              ZETA3( IX,IZ ) 
     &    =   ZETA2( IX,IZ ) + DT * ZETA3( IX, IZ )
 3010   CONTINUE
 3020 CONTINUE
*
      RETURN
      END
**********************************************************************
      SUBROUTINE DEULTZ          !" Euler scheme による渦度の時間変化
     O             ( TNDZ, 
     I               PSI   , ZETA  , T   , 
     C               RAYL  , PRND  ,
     C               DX    , DZ    , 
     W               AJAC  , ALAP    )
*
* [PARAM]
#include        "zcdim.F"            !" 格子点数, 波数
*
* [OUTPUT]
      REAL TNDZ( 0:NX+1 , 0:NZ+1 ) !" 渦度の時間変化
*
* [INPUT]
      REAL PSI   ( 0:NX+1 , 0:NZ+1 ) !" 流線関数(1STEP前の値)
      REAL ZETA  ( 0:NX+1 , 0:NZ+1 ) !" 渦度(1STEP前の値)
      REAL T     ( 0:NX+1 , 0:NZ+1 ) !" 温度(1STEP前の値)
*
      REAL PRND                      !" プラントル数
      REAL RAYL                      !" レイリー数
      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  , ZETA ,
     C         DX   , DZ    )
*
*"        < 2. ラプラシアンの計算 >
*
      CALL DLAPLA
     O       ( ALAP ,
     I         ZETA ,
     C         DX   , DZ   )
*
*"        < 3. 各グリッドの渦度の時間変化の計算 >
*
      DO 3020 IZ = 1 , NZ
        DO 3010 IX = 1 , NX
               TNDZ( IX,IZ ) 
     &     =   AJAC(IX,IZ) + PRND * ALAP(IX,IZ)
     &       - PRND*RAYL*( T(IX+1,IZ)-T(IX-1,IZ) )/(2.*DX) 
 3010   CONTINUE
 3020 CONTINUE
*
      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
