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