* PACKAGE DLPFRG    !" 時間積分( Leap Frog Scheme )
*
*"  [HIS] 92/01/09 (momoko)
*"        95/08/14 (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 により次ステップの物理量を計算する >>
*"                                      92/01/13 (石渡正樹)
*
* [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/08/14' )
      ENDIF
*
*"        < 1. 内部領域の渦度の計算 >
*
      CALL DLPFZ
     O      ( ZETA3 , 
     I        PSI2  , ZETA2 , T2  , 
     I        ZETA1 , 
     C        RAYL  , PRND  ,
     C        DX    , DZ    , DT  ,
     W        AJAC  , ALAP   )
*
*"        < 2. 内部領域の温度の計算 >
*
      CALL DLPFT
     O      ( T3    , 
     I        PSI2  , T2    , 
     I        T1    , 
     C        QINT  , 
     C        DX    , DZ    , DT  ,
     W        AJAC  , ALAP    )
*
*"        < 3. 境界条件の適用と流線関数の計算 >
*
      CALL DBNDR
     O       ( PSI3  , 
     M         ZETA3 , T3   ,
     I         DX    , DZ    )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE DLPFZ     !" Leapfrog Scheme による渦度
     O             ( ZETA3 , 
     I               PSI2  , ZETA2 , T2  , 
     I               ZETA1 , 
     C               RAYL  , PRND  ,
     C               DX    , DZ    , DT  ,
     W               AJAC2 , ALAP1   )
*
* [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 ZETA1 ( 0:NX+1 , 0:NZ+1 ) !"渦度(2STEP前の値)
*
      REAL PRND                      !"プラントル数
      REAL RAYL                      !"レイリー数
      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 , ZETA2,
     D         DX   , DZ    )
*
*"        < 2. ラプラシアンの計算 >
*
      CALL DLAPLA
     O       ( ALAP1,
     I         ZETA1,
     C         DX   , DZ   )
*
*"        < 3. 各グリッドの渦度の値の計算 >
*
      DO 3020 IZ = 1 , NZ
        DO 3020 IX = 1 , NX
              ZETA3( IX,IZ ) 
     &    =   ZETA1( IX,IZ )
     &      + 2.0D0*DT
     &        *(   AJAC2(IX,IZ) + PRND*ALAP1(IX,IZ)
     &           - PRND*RAYL*( T2(IX+1,IZ) - T2(IX-1,IZ) )/(2.0D0*DX)
     &         )
 3010   CONTINUE
 3020 CONTINUE
*
      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
