*" PACKAGE AKEANL !" 運動エネルギー生成散逸計算 *" ver 1.00 94/01/25 takepiro *" ver 1.01 95/06/12 takepiro Bug fix *********************************************************************** PROGRAM AEDYFL !" メインルーチン * *" Boussinesq 全球モデル *" 非線形項の計算, 出力 * * [PARAM] #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数, 波数 INCLUDE (ZHDIM) !" 文字列文字数 INCLUDE (ZCCOM) !" 標準物理定数 #else #include "zcdim.F" !" 格子点数, 波数 #include "zhdim.F" !" 文字列文字数 #include "zccom.F" !" 標準物理定数 #endif * * [VAR] REAL GAU ( IDIM, JDIM, 0:KDIM ) !" 西風 ｕ REAL GAV ( IDIM, JDIM, 0:KDIM ) !" 南風 ｖ REAL GAW ( IDIM, JDIM, 0:KDIM ) !" 鉛直風 ｗ REAL GAT ( IDIM, JDIM, 0:KDIM ) !" 温度 Ｔ REAL GATOR ( IDIM, JDIM, 0:KDIM ) !" トロイダル Ψ REAL GAPOR ( IDIM, JDIM, 0:KDIM ) !" ポロイダル Φ *" : 格子点データ(t) 格子点データ(t+Δt) * REAL GEXX ( IDIM, JDIM, 0:KDIM ) !" 変形テンソル e_xx REAL GEYY ( IDIM, JDIM, 0:KDIM ) !" 変形テンソル e_yy REAL GEZZ ( IDIM, JDIM, 0:KDIM ) !" 変形テンソル e_zz REAL GEXY ( IDIM, JDIM, 0:KDIM ) !" 変形テンソル e_xy REAL GEYZ ( IDIM, JDIM, 0:KDIM ) !" 変形テンソル e_yz REAL GEZX ( IDIM, JDIM, 0:KDIM ) !" 変形テンソル e_zx * REAL GKEGE ( IDIM, JDIM, 0:KDIM ) !" 運動エネルギー生成 REAL GKEDS ( IDIM, JDIM, 0:KDIM ) !" 運動エネルギー散逸 * *" : 格子点変化項(D) * INTEGER ISTEP !" 通しステップ数 INTEGER ITA !" 通し時間(t), 標準時間単位 c$$$ INTEGER IDATEA( 6 ) !" 時刻年月日時分秒(ダミー) *" : 時刻等 * * [ONCE] INTEGER ITSTRT !" 計算開始時刻(標準時間単位) INTEGER ITEND !" 計算終了時刻(標準時間単位) INTEGER IORSTR !" 出力間隔：再出発(標準時間単位) *" : 実験管理パラメータ * INTEGER IDELT !" 標準時間刻み(SEC) *" : 時間ステップパラメータ * *" COMMON /COMCON/ ( include ZCCOM ) * REAL ERI !" 内径 * REAL ERO !" 外径 * REAL ETA !" 内径/外径 *" COMMON /COMCON/ end *" : 物理パラメータ * REAL ALAT ( JDIM ) !" 緯度 REAL DLAT ( JDIM ) !" 緯度荷重 REAL ALON ( IDIM ) !" 経度 REAL DLON ( IDIM ) !" 経度荷重 REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) *" : 座標値 * * [WORK] COMMON /COMWRK/ WORK REAL WORK ( NWORK ) !" ワーク領域 * INTEGER IUNIT !" 標準出力装置番号 * c$$$ DATA IDATEA / 1992, 12, 18, 0, 0, 0/ * #ifdef DEBUG external common_handler integer i, ieee_handler i = ieee_handler("set","common",common_handler ) if( i.ne. 0 ) print *, "Could not establish fp signal handler" #endif * *" << SETPUP : 初期設定 >> * CALL GLPGET( 'MSGUNIT', IUNIT ) WRITE ( IUNIT, * )' KINETIC ENERGY GENERATE/DISSIPATION AKEANL' & //NAME & //' ver.1.01, 95/06/12' * CALL CLCSTR ( 'SETUP' ) CALL YPREP !" システム前処理 * CALL SETPAR !" 実験パラメータ O ( ITSTRT , ITEND , IORSTR , O IDELT ) * CALL SETCOR !" 座標値 O ( ALON , DLON , O ALAT , DLAT , O ARAD , DRAD ) * CALL RDSTRT !" 初期値の読み込みと再生成 O ( GAU , GAV , GAW , GAT , O GATOR , GAPOR , O ITA , ISTEP , C ITSTRT, ALAT , DLAT , ARAD , DRAD ) * CALL ADMXMN I ( GAU , GAV , GAW , GAT , I GATOR , GAPOR , ITA , 'start GA' ) * CALL AHSTRG !" 出力設定 * CALL SETTIM !" 時刻を記憶 I ( ITA-IDELT, ISTEP-1 ) * CALL HISTOU * CALL SETTIM !" 時刻を記憶 I ( ITA , ISTEP ) * * *" << 運動エネルギー散逸計算 >> * CALL DKEANL O ( GKEGE, GKEDS, O GEXX , GEYY , GEZZ , O GEXY , GEYZ , GEZX , I GAU , GAV , GAW , GAT , C ALON , DLON , ALAT , DLAT , C ARAD , DRAD ) * *" << 散逸項の出力 >> * CALL HISTIN ( GKEGE , 'KEGE' ) !" 出力領域へ CALL HISTIN ( GKEDS , 'KEDS' ) CALL HISTIN ( GEXX , 'EXX' ) CALL HISTIN ( GEYY , 'EYY' ) CALL HISTIN ( GEZZ , 'EZZ' ) CALL HISTIN ( GEXY , 'EXY' ) CALL HISTIN ( GEYZ , 'EYZ' ) CALL HISTIN ( GEZX , 'EZX' ) * CALL HISTOU * CALL YFINE !" システム後処理 * STOP END *********************************************************************** SUBROUTINE DKEANL !" 運動エネルギー生成散逸項計算 O ( GKEGE, GKEDS, O GEXX , GEYY , GEZZ , O GEXY , GEYZ , GEZX , I GDU , GDV , GDW , GDT , C ALON , DLON , ALAT , DLAT , C ARAD , DRAD ) * * [PARAM] #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数, 波数 INCLUDE (ZHDIM) !" 文字列文字数 #else #include "zcdim.F" !" 格子点数, 波数 #include "zhdim.F" !" 文字列文字数 #endif * *" [OUTPUT] REAL GKEGE ( IDIM*JDIM, 0:KDIM ) !" 運動エネルギー生成 REAL GKEDS ( IDIM*JDIM, 0:KDIM ) !" 運動エネルギー散逸 * REAL GEXX ( IDIM*JDIM, 0:KDIM ) !" 変形テンソル e_xx REAL GEYY ( IDIM*JDIM, 0:KDIM ) !" 変形テンソル e_yy REAL GEZZ ( IDIM*JDIM, 0:KDIM ) !" 変形テンソル e_zz REAL GEXY ( IDIM*JDIM, 0:KDIM ) !" 変形テンソル e_xy REAL GEYZ ( IDIM*JDIM, 0:KDIM ) !" 変形テンソル e_yz REAL GEZX ( IDIM*JDIM, 0:KDIM ) !" 変形テンソル e_zx * *" [INPUT] REAL GDU ( IDIM*JDIM, 0:KDIM ) !" 西風 ｕ REAL GDV ( IDIM*JDIM, 0:KDIM ) !" 南風 ｖ REAL GDW ( IDIM*JDIM, 0:KDIM ) !" 鉛直風 ｗ REAL GDT ( IDIM*JDIM, 0:KDIM ) !" 温度 Ｔ * REAL ALAT ( JDIM ) !" 緯度 REAL DLAT ( JDIM ) !" 緯度荷重 REAL ALON ( IDIM ) !" 経度 REAL DLON ( IDIM ) !" 経度荷重 REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) * * [INTERNAL WORK] REAL GDVDX ( IDIM*JDIM, 0:KDIM ) !" 速度場の X 微分 REAL GDVDY ( IDIM*JDIM, 0:KDIM ) !" 速度場の Y 微分 REAL GDVDZ ( IDIM*JDIM, 0:KDIM ) !" 速度場の Z 微分 * * [INTERNAL ONCE] INTEGER NMO ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番 REAL FLAPLA( NMDIM ) !" ラプラシアンの係数 REAL EDEL ( NMDIM ) !" ζ，Ｄ→Ｕ，Ｖ REAL UVFACT ( IDIM, JDIM ) !" u→U のファクター REAL TANG ( IDIM*JDIM ) !" tan(φ) * REAL RAYL !" レイリー数 REAL TAU !" √テイラー数 REAL PRND !" プランドル数 LOGICAL ONONLN !" 非線形計算スイッチ(ダミー) * INTEGER I, J, K, IJ * LOGICAL OFIRST DATA OFIRST / .TRUE. / SAVE OFIRST , & NMO , FLAPLA , EDEL , UVFACT, TANG , & RAYL , TAU , PRND * *" << 0. 定数設定 >> * IF ( OFIRST ) THEN OFIRST = .FALSE. * *" < 0.1 スペクトル定数設定 > * CALL DSETC O ( NMO , FLAPLA, EDEL , UVFACT, C ALAT , DLAT ) * IJ = 0 DO 100 J = 1, JDIM DO 100 I = 1, IDIM IJ = IJ + 1 TANG( IJ ) = TAN( ALAT( J ) ) 100 CONTINUE * *" < 0.2 力学定数設定 > * CALL DPARM !" パラメター設定 O ( RAYL , TAU , PRND , O ONONLN ) * CALL HSPARM I ( 'RAYLEIGH NUMBER', RAYL ) CALL HSPARM I ( 'TAYLOR NUMBER' , TAU**2 ) CALL HSPARM I ( 'PRANDTL NUMBER' , PRND ) * ENDIF * *" << 1. 変形テンソルの計算 >> * CALL RESET( GEXX , IDIM*JDIM*KMAX ) CALL RESET( GEYY , IDIM*JDIM*KMAX ) CALL RESET( GEZZ , IDIM*JDIM*KMAX ) CALL RESET( GEXY , IDIM*JDIM*KMAX ) CALL RESET( GEYZ , IDIM*JDIM*KMAX ) CALL RESET( GEZX , IDIM*JDIM*KMAX ) * *" < 1.1 ｕ の微分 > * CALL VDERIV O ( GDVDX , GDVDY , GDVDZ , I GDU , D ARAD , DRAD , C UVFACT ) * DO 1100 K = 0, KDIM DO 1100 IJ = 1, IDIM*JDIM GEXX( IJ, K ) = GEXX( IJ, K ) & + 2 * GDVDX( IJ, K ) & - 2 * GDV ( IJ, K ) * TANG( IJ ) / ARAD( K ) & + 2 * GDW ( IJ, K ) / ARAD( K ) GEXY( IJ, K ) = GEXY( IJ, K ) & + GDVDY( IJ, K ) & + GDU( IJ, K ) * TANG( IJ ) / ARAD( K ) GEZX( IJ, K ) = GEZX( IJ, K ) & + GDVDZ( IJ, K ) & - GDU( IJ,K ) / ARAD( K ) 1100 CONTINUE * *" < 1.2 ｖ の微分 > * CALL VDERIV O ( GDVDX , GDVDY , GDVDZ , I GDV , D ARAD , DRAD , C UVFACT ) * DO 1200 K = 0, KDIM DO 1200 IJ = 1, IDIM*JDIM GEXY( IJ, K ) = GEXY( IJ, K ) & + GDVDX( IJ, K ) GEYY( IJ, K ) = GEYY( IJ, K ) & + 2 * GDVDY( IJ, K ) & + 2 * GDW( IJ, K ) / ARAD( K ) GEYZ( IJ, K ) = GEYZ( IJ, K ) & + GDVDZ( IJ, K ) & - GDV( IJ,K ) / ARAD( K ) 1200 CONTINUE * *" < 1.3 ｗ の微分 > * CALL VDERIV O ( GDVDX , GDVDY , GDVDZ , I GDW , D ARAD , DRAD , C UVFACT ) * DO 1300 K = 0, KDIM DO 1300 IJ = 1, IDIM*JDIM GEZX( IJ, K ) = GEZX( IJ, K ) & + GDVDX( IJ, K ) GEYZ( IJ, K ) = GEYZ( IJ, K ) & + GDVDY( IJ, K ) GEZZ( IJ, K ) = GEZZ( IJ, K ) & + 2 * GDVDZ( IJ, K ) 1300 CONTINUE * *" << 2. 運動エネルギー散逸の計算 >> * DO 2000 K = 0, KDIM DO 2000 IJ = 1, IDIM*JDIM GKEDS( IJ, K ) = & 0.5 * PRND * & ( GEXX( IJ,K )**2 & + GEYY( IJ,K )**2 & + GEZZ( IJ,K )**2 & + 2 * GEXY( IJ,K )**2 & + 2 * GEYZ( IJ,K )**2 & + 2 * GEZX( IJ,K )**2 ) 2000 CONTINUE * *" << 3. 運動エネルギー生成の計算 >> * DO 3000 K = 0, KDIM DO 3000 IJ = 1, IDIM*JDIM GKEGE( IJ,K ) = RAYL * PRND & * ARAD( K ) * GDT( IJ,K ) * GDW( IJ,K ) 3000 CONTINUE * RETURN END *********************************************************************** SUBROUTINE VDERIV !" 速度場の微分計算 O ( GDVDX , GDVDY , GDVDZ , I GDV , D ARAD , DRAD , C UVFACT ) * * [PARAM] #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数, 波数 INCLUDE (ZHDIM) !" 文字列文字数 #else #include "zcdim.F" !" 格子点数, 波数 #include "zhdim.F" !" 文字列文字数 #endif * * [INPUT] REAL GDV ( IDIM*JDIM, 0:KDIM ) !" 速度場 * REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) REAL UVFACT ( IDIM*JDIM ) !" u→U のファクター * * [OUTPUT] REAL GDVDX ( IDIM*JDIM, 0:KDIM ) !" 速度場の X 微分 REAL GDVDY ( IDIM*JDIM, 0:KDIM ) !" 速度場の Y 微分 REAL GDVDZ ( IDIM*JDIM, 0:KDIM ) !" 速度場の Z 微分 * * [INTERNAL WORK] REAL WDV ( NMDIM, 0:KDIM ) !" 速度場スペクトルデータ REAL WORK( NMDIM, 0:KDIM ) !" INTEGER IJ , K * *" < 0. スペクトル変換 > * CALL G2W O ( WDV , I GDV , I ' ' , ' ' , KMAX ) * *" < 1. X 微分 > * CALL W2G O ( GDVDX , I WDV , I 'X' , 'POS' , KMAX ) * DO 1100 K = 0, KDIM DO 1100 IJ = 1, IDIM*JDIM GDVDX( IJ, K ) = GDVDX( IJ, K ) & / UVFACT( IJ ) / ARAD( K ) 1100 CONTINUE * *" < 2. Y 微分 > * CALL W2G O ( GDVDY , I WDV , I 'Y' , 'POS' , KMAX ) * DO 1200 K = 0, KDIM DO 1200 IJ = 1, IDIM*JDIM GDVDY( IJ, K ) = GDVDY( IJ, K ) & / UVFACT( IJ ) / ARAD( K ) 1200 CONTINUE * *" < 3. Z 微分 > * CALL W2C O ( WORK , I WDV , I ' ' , 'POS' , KMAX ) * CALL C2G O ( GDVDZ , I WORK , I 'DR ' , 'POS' , KMAX ) * RETURN END *********************************************************************** SUBROUTINE AHSTRG !" 標準時間平均出力の登録 * * [INTERN PARAM] INTEGER ISTYPL DATA ISTYPL / 1 / REAL VMISS * CALL GZDBGT ( 'MISS', VMISS ) !" 欠損値 * CALL HISTRG !" 出力の登録 I ( 'KEGE ', 'k.energy generate ' ,'W/m^3 ', 'ALEV', I VMISS , VMISS , 0.1 , 0.5 , ISTYPL , I ' ' , ' ' , 0 , 0 , ' ' ,'(F12.3)' ) * CALL HISTRG I ( 'KEDS ', 'k.energy dissip ' ,'W/m^3 ', 'ALEV', I VMISS , VMISS , 0.1 , 0.5 , ISTYPL , I ' ' , ' ' , 0 , 0 , ' ' ,'(F12.3)' ) * CALL HISTRG I ( 'EXX ', 'xx-strain ' ,'1/s ', 'ALEV', I VMISS , VMISS , 0.1 , 0.5 , ISTYPL , I 'CON' , ' ' , 0 , 0 , ' ' ,'(F12.3)' ) * CALL HISTRG !" 出力の登録 I ( 'EYY ', 'yy-strain ' ,'1/s ', 'ALEV', I VMISS , VMISS , 0.1 , 0.5 , ISTYPL , I 'CON' , ' ' , 0 , 0 , ' ' ,'(F12.3)' ) * CALL HISTRG !" 出力の登録 I ( 'EZZ ', 'zz-strain ' ,'1/s ', 'ALEV', I VMISS , VMISS , 0.1 , 0.5 , ISTYPL , I 'CON' , ' ' , 0 , 0 , ' ' ,'(F12.3)' ) * CALL HISTRG !" 出力の登録 I ( 'EXY ', 'xy-strain ' ,'1/s ', 'ALEV', I VMISS , VMISS , 0.1 , 0.5 , ISTYPL , I 'CON' , ' ' , 0 , 0 , ' ' ,'(F12.3)' ) * CALL HISTRG !" 出力の登録 I ( 'EYZ ', 'yz-strain ' ,'1/s ', 'ALEV', I VMISS , VMISS , 0.1 , 0.5 , ISTYPL , I 'CON' , ' ' , 0 , 0 , ' ' ,'(F12.3)' ) * CALL HISTRG !" 出力の登録 I ( 'EZX ', 'zx-strain ' ,'1/s ', 'ALEV', I VMISS , VMISS , 0.1 , 0.5 , ISTYPL , I 'CON' , ' ' , 0 , 0 , ' ' ,'(F12.3)' ) * RETURN END ************************************************************************* SUBROUTINE ADMXMN !" debug monitor I ( GDU , GDV , GDW , GDT , I GDTOR , GDPOR , IT , HLABEL ) * * [PARAM] #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数, 波数 #else #include "zcdim.F" !" 格子点数, 波数 #endif * * [INPUT] REAL GDU ( IDIM, JDIM, KMAX ) !" 西風 ｕ REAL GDV ( IDIM, JDIM, KMAX ) !" 南風 ｖ REAL GDW ( IDIM, JDIM, KMAX ) !" 鉛直風 ｗ REAL GDT ( IDIM, JDIM, KMAX ) !" 温度 Ｔ REAL GDTOR ( IDIM, JDIM, KMAX ) !" トロイダル Ψ REAL GDPOR ( IDIM, JDIM, KMAX ) !" ポロイダル Φ * INTEGER IT CHARACTER HLABEL *(*) * WRITE ( 6,* ) HLABEL, ' IT= ', IT CALL MAXMIN( GDU, IDIM, JDIM, KMAX, IDIM, JDIM, 'U' ) CALL MAXMIN( GDV, IDIM, JDIM, KMAX, IDIM, JDIM, 'V' ) CALL MAXMIN( GDW, IDIM, JDIM, KMAX, IDIM, JDIM, 'W' ) CALL MAXMIN( GDT, IDIM, JDIM, KMAX, IDIM, JDIM, 'T' ) CALL MAXMIN( GDTOR, IDIM, JDIM, KMAX, IDIM, JDIM, 'TOR' ) CALL MAXMIN( GDPOR, IDIM, JDIM, KMAX, IDIM, JDIM, 'POR' ) * RETURN END