* PACKAGE DHEUN !" 時間積分( Heun Scheme ) * *" [HIS] 92/01/13 (momoko) *" 95/08/14 (takepiro) *" 96/03/16 (takepiro) *********************************************************************** SUBROUTINE DHEUN O ( PSI3 , ZETA3 , T3 , I PSI2 , ZETA2 , T2 , I RAYL , PRND , QINT, C DX , DZ , DT , W TNDZ , TNDT , W AJAC , ALAP ) * *" << Heun 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 TNDZ ( 0:NX+1 , 0:NZ+1 ) !" 第 1 段階での渦度の時間変化 REAL TNDT ( 0:NX+1 , 0:NZ+1 ) !" 第 1 段階での温度の時間変化 REAL AJAC ( 0:NX+1 , 0:NZ+1 ) !"ヤコビアン REAL ALAP ( 0:NX+1 , 0:NZ+1 ) !"ラプラシアン * * [INTERNAL WORK] INTEGER IX !"ループカウンタ INTEGER IZ !"ループカウンタ * LOGICAL OFIRST DATA OFIRST / .TRUE. / * SAVE OFIRST * IF( OFIRST ) THEN OFIRST = .FALSE. CALL MSGDMP('M','DHEUN', & 'HEUN SCHEME TIME INTEGRATE : 95/08/14' ) ENDIF * *" < 1. Euler scheme : dummyの値を計算 > * CALL DEULTZ O ( TNDZ , I PSI2 , ZETA2 , T2 , C RAYL , PRND , C DX , DZ , W AJAC , ALAP ) * DO 1020 IZ = 1 , NZ DO 1010 IX = 1 , NX ZETA3( IX,IZ ) & = ZETA2( IX,IZ ) + DT * TNDZ( IX, IZ ) 1010 CONTINUE 1020 CONTINUE * CALL DEULTT O ( TNDT , I PSI2 , T2 , C QINT , C DX , DZ , W AJAC , ALAP ) * DO 1040 IX = 1 , NX DO 1030 IZ = 1 , NZ T3(IX,IZ) & = T2(IX,IZ) + DT * TNDT( IX, IZ ) 1030 CONTINUE 1040 CONTINUE * CALL DBNDR !" 境界条件の適用と流線関数の計算 O ( PSI3 , M ZETA3 , T3 , I DX , DZ ) * *" < 2. 反復する > * CALL DHEUZZ M ( ZETA3 , I PSI3 , T3 , I ZETA2 , TNDZ , C RAYL , PRND , C DX , DZ , DT , W AJAC , ALAP ) * CALL DBNDRZ !" 運動学/力学境界条件の適用と流線関数の計算 O ( PSI3 , M ZETA3, I DX , DZ ) * CALL DHEUZT M ( T3 , I PSI3 , I T2 , TNDT , C QINT , C DX , DZ , DT , W AJAC , ALAP ) * CALL DBNDRT !" 熱的境界条件の適用 M ( T3 , I DX , DZ ) * RETURN END *********************************************************************** SUBROUTINE DHEUZZ !" Heun scheme による渦度(第 2 段階) M ( ZETA3 , I PSI3 , T3 , I ZETA2 , TNDZ , C RAYL , PRND , C DX , DZ , DT , W AJAC3 , ALAP3 ) * * [PARAM] #include "zcdim.F" !" 格子点数, 波数 * * [MODIFY] REAL ZETA3 ( 0:NX+1 , 0:NZ+1 ) !" 渦度(現STEPの値) * * [INPUT] REAL PSI3 ( 0:NX+1 , 0:NZ+1 ) !" 流線関数(現STEPの値, 予想値) REAL T3 ( 0:NX+1 , 0:NZ+1 ) !" 温度(現STEP前の値) REAL ZETA2 ( 0:NX+1 , 0:NZ+1 ) !" 渦度(1STEP前の値) REAL TNDZ ( 0:NX+1 , 0:NZ+1 ) !" 第 1 段階での渦度の時間変化 * REAL PRND !" プラントル数 REAL RAYL !" レイリー数 REAL DX !" X 軸上の格子点間隔 REAL DZ !" Z 軸上の格子点間隔 REAL DT !" 時間間隔 * * [WORK] REAL AJAC3 ( 0:NX+1 , 0:NZ+1 ) !"ヤコビアン REAL ALAP3 ( 0:NX+1 , 0:NZ+1 ) !"ラプラシアン * * [INTERNAL WORK] INTEGER IX !"ループカウンタ! INTEGER IZ !"ループカウンタ * *" < 1. ヤコビアンの計算 > * CALL DARJAC O ( AJAC3, I PSI3 , ZETA3, C DX , DZ ) * *" < 2. ラプラシアンの計算 > * CALL DLAPLA O ( ALAP3, I ZETA3, C DX , DZ ) * *" < 3. 各グリッドの渦度の値の計算 > * DO 3020 IZ = 1 , NZ DO 3010 IX = 1 , NX ZETA3( IX,IZ ) & = ZETA2( IX,IZ ) & + DT * ( & 0.5*TNDZ(IX,IZ) & + 0.5*( AJAC3(IX,IZ) + PRND * ALAP3(IX,IZ) & - PRND*RAYL*( T3(IX+1,IZ)-T3(IX-1,IZ) )/(2.*DX) ) & ) 3010 CONTINUE 3020 CONTINUE * RETURN END ***************************************************************** SUBROUTINE DHEUZT !" Heun scheme による温度(第 2 段階) M ( T3 , I PSI3 , I T2 , TNDT , C QINT , C DX , DZ , DT , W AJAC3 , ALAP3 ) * * [PARAM] #include "zcdim.F" !" 格子点数, 波数 * * [OUTPUT] REAL T3 ( 0:NX+1 , 0:NZ+1 ) !" 温度(現STEPの値, 予想値) * * [INPUT] REAL PSI3 ( 0:NX+1 , 0:NZ+1 ) !" 流線関数(現STEPの値) REAL T2 ( 0:NX+1 , 0:NZ+1 ) !" 温度(1STEP前の値) REAL TNDT ( 0:NX+1 , 0:NZ+1 ) !" 第 1 段階での温度の時間変化 * REAL QINT !" 加熱率 REAL DX !" X 軸上の格子点間隔 REAL DZ !" Z 軸上の格子点間隔 REAL DT !" 時間間隔 * * [WORK] REAL AJAC3 ( 0:NX+1 , 0:NZ+1 ) !" ヤコビアン REAL ALAP3 ( 0:NX+1 , 0:NZ+1 ) !" ラプラシアン * * [INTERNAL WORK] INTEGER IX !" ループカウンタ INTEGER IZ !" ループカウンタ * *" < 1. ヤコビアンの計算 > * CALL DARJAC O ( AJAC3, I PSI3 , T3 , C DX , DZ ) * *" < 2. ラプラシアンの計算 > * CALL DLAPLA O ( ALAP3, I T3 , C DX , DZ ) * *" < 3. 各グリッドの温度の値の計算 > * DO 3020 IX = 1 , NX DO 3010 IZ = 1 , NZ T3( IX,IZ ) & = T2( IX,IZ ) & + DT*( 0.5 * TNDT( IX,IZ ) & + 0.5 * ( AJAC3(IX,IZ) + ALAP3(IX,IZ) + QINT ) ) 3010 CONTINUE 3020 CONTINUE * RETURN END