* 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