* PACKAGE LSTEDY !" 線形力学定常解(非線形項に対する応答) *" ver 1.02 93/09/27 takepiro ********************************************************************** SUBROUTINE LSTEDY !" 線形力学定常解 M ( GAU , GAV , GAW , GAT , M GATOR , GAPOR , I MFIX , C ALON , DLON , ALAT , DLAT , C ARAD , DRAD ) * * [PARAM] #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数, 波数 INCLUDE (ZLDIM) !" パラメター数, 係数行列 INCLUDE (ZHDIM) !" 文字列文字数 INCLUDE (ZCCOM) !" 標準物理定数 #else #include "zcdim.F" !" 格子点数, 波数 #include "zldim.F" !" パラメター数, 係数行列 #include "zhdim.F" !" 文字列文字数 #include "zccom.F" !" 標準物理定数 #endif * * [MODIFY] 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) 格子点データ *" 入力-非線形項を計算, 出力-非線形項に対する応答 * * [INPUT] INTEGER MFIX !" 定常解の東西波数 * REAL ALON ( IDIM ) !" 経度 REAL DLON ( IDIM ) !" 経度荷重 REAL ALAT ( JDIM ) !" 緯度 REAL DLAT ( JDIM ) !" 緯度荷重 REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) *" : 座標格子 * * [INTERNAL ONCE] INTEGER NMO ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番 REAL FLAPLA( NMDIM ) !" ラプラシアンの係数 REAL EDEL ( NMDIM ) !" ζ，Ｄ→Ｕ，Ｖ * REAL UVFACT ( IDIM, JDIM ) !" u→U のファクター * REAL RAYL !" レイリー数 REAL TAU !" √テイラー数 REAL PRND !" プランドル数 * SAVE NMO , FLAPLA, EDEL , & UVFACT * LOGICAL OSETC , OFIRST DATA OSETC , OFIRST / .FALSE., .TRUE. / SAVE OSETC , OFIRST * LOGICAL ONONLN !" 非線形力学スイッチ * LOGICAL OATOR , OAPOR , OAT !" 解析対象物理量 SAVE OATOR , OAPOR , OAT * * [INTERNAL WORK] REAL GTUA ( IDIM, JDIM, 0:KDIM ) !" 東西運動量変化項ＵＡ REAL GTVA ( IDIM, JDIM, 0:KDIM ) !" 南北運動量変化項ＶＡ REAL GTWA ( IDIM, JDIM, 0:KDIM ) !" 鉛直運動量変化項ＶＡ REAL GTH ( IDIM, JDIM, 0:KDIM ) !" 温度時間変化項 Ｈ *" : 格子点変化項を足し込む * REAL GTUT ( IDIM, JDIM, 0:KDIM ) !" 温度東西移流項 ＵＴ REAL GTVT ( IDIM, JDIM, 0:KDIM ) !" 温度南北移流項 ＶＴ REAL GTWT ( IDIM, JDIM, 0:KDIM ) !" 温度鉛直移流項 ＷＴ *" : 時間変化項 * REAL WTTOR ( NMDIM, 0:KDIM ) !" トロイダル Ψ REAL WTPOR ( NMDIM, 0:KDIM ) !" ポロイダル Φ REAL WTT ( NMDIM, 0:KDIM ) !" 温度 Ｔ *" : スペクトル時間変化項 * REAL COFMTX( MATDMX , MATDMX ) !" 係数行列 REAL VSTEDY( MATDMX ) !" 非線形項 / 定常解 * INTEGER NMKMTX ( NMDIM, 0:KDIM, 3 ) !" 係数行列の添え字 INTEGER NMKMAX !" 係数行列の最大添え字 INTEGER NMTLST ( NMDIM ) !" MFIX 成分リスト(トロイダル) INTEGER NMPLST ( NMDIM ) !" MFIX 成分リスト(ポロイダル) INTEGER NMQLST ( NMDIM ) !" MFIX 成分リスト(温度) INTEGER NMTMAX !" NMTLST の最大添え字 INTEGER NMPMAX !" NMPLST の最大添え字 INTEGER NMQMAX !" NMQLST の最大添え字 * INTEGER IERR !" サブルーチンエラーコード REAL EPS !" 特異性の判定基準 REAL EPS0 !" 特異性の判定基準 DATA EPS0 / 3.52D-15 / *" : 固有値計算用行列 * INTEGER IUNIT !" 標準出力装置番号 INTEGER I , J * *" << 0.SETCON : 定数設定 >> * * [ONCE] IF ( .NOT. OSETC ) THEN OSETC = .TRUE. CALL DSETC !" スペクトル定数設定 O ( NMO , FLAPLA, EDEL , UVFACT, C ALAT , DLAT ) ENDIF * * [ONCE] IF ( OFIRST ) THEN OFIRST = .FALSE. CALL DPARM !" パラメター設定 O ( RAYL , TAU , PRND , O ONONLN ) * IF( ONONLN )THEN CALL MSGDMP('W','LSTDY', & 'N-LIN. DYNAMICS SWITCH WAS NEGLECTED') ENDIF * CALL HSPARM I ( 'RAYLEIGH NUMBER', RAYL ) CALL HSPARM I ( 'TAYLOR NUMBER' , TAU**2 ) CALL HSPARM I ( 'PRANDTL NUMBER' , PRND ) * c$$$ CALL DHSTRG !" 出力設定 * CALL LSTPHS !" 解析対象物理量の設定 O ( OATOR , OAPOR , OAT ) ENDIF * *" << 1. 係数行列の計算 : 線型力学 >> * CALL CLCSTR( 'LSTMTX' ) * CALL MKMLST O ( NMKMTX , NMKMAX , O NMTLST , NMTMAX , O NMPLST , NMPMAX , O NMQLST , NMQMAX , F OATOR , OAPOR , OAT , I NMO , MFIX ) * CALL LSTMTX O ( COFMTX , I NMKMTX , NMKMAX , I NMTLST , NMTMAX , I NMPLST , NMPMAX , I NMQLST , NMQMAX , I RAYL , TAU , PRND , F OATOR , OAPOR , OAT , C NMO , FLAPLA , C ARAD , DRAD ) * CALL CLCEND( 'LSTMTX' ) * *" << 2. 非線形項の計算 >> * CALL CLCSTR( 'W2VECT' ) * CALL GRDDYN !" 非線型力学項 O ( GTUA , GTVA , GTWA , O GTUT , GTVT , GTWT , GTH , I GAU , GAV , GAW , GAT , C UVFACT, FLAPLA, C ALAT , DLAT , ARAD , DRAD ) * CALL GRDFRC !" 力学外部強制項(格子点) M ( GTUA , GTVA , GTWA , M GTUT , GTVT , GTWT , GTH , C ALON , DLON , ALAT , DLAT , C ARAD , DRAD ) * CALL TENG2W !" 時間変化項をスペクトルに O ( WTTOR , WTPOR , WTT , I GTUA , GTVA , GTWA , I GTUT , GTVT , GTWT , GTH , C FLAPLA, C ARAD , DRAD ) * CALL W2VECT !" 非線形項をベクトルに O ( VSTEDY , I WTTOR , WTPOR , WTT , I NMKMTX , NMKMAX , I NMTLST , NMTMAX , I NMPLST , NMPMAX , I NMQLST , NMQMAX , F OATOR , OAPOR , OAT ) * DO 2000 I = 1, NMKMAX VSTEDY( I ) = - VSTEDY( I ) 2000 CONTINUE * CALL CLCSTR( 'W2VECT' ) * *" << 3. 定常解計算 >> * CALL CLCSTR( 'LGLU' ) * EPS = EPS0 * PRND !" 特異性の基準値を設定 * CALL LEIGLU !" 連立 1 次方程式を解く c$$$ CALL LSVDSL !" 特異連立 1 次方程式を解く M ( COFMTX , M VSTEDY , O IERR , I EPS , I NMKMAX , .TRUE. ) * IF ( IERR .EQ. 0 )THEN CALL LEIVEC O ( GAU , GAV , GAW , GAT , O GATOR , GAPOR , I VSTEDY , VSTEDY , I NMKMTX , I NMTLST , NMTMAX , I NMPLST , NMPMAX , I NMQLST , NMQMAX , F OATOR , OAPOR , OAT , C NMO , FLAPLA , UVFACT , C ARAD , DRAD ) ENDIF CALL CLCSTR( 'LGLU' ) * * *" << OUTPUT : 出力 >> * CALL HISTIN ( GTUA , 'UA' ) !" 出力領域へ CALL HISTIN ( GTVA , 'VA' ) CALL HISTIN ( GTWA , 'WA' ) CALL HISTIN ( GTUT , 'UT' ) !" 出力領域へ CALL HISTIN ( GTVT , 'VT' ) CALL HISTIN ( GTWT , 'WT' ) CALL HISTIN ( GTH , 'H' ) * RETURN END ********************************************************************** SUBROUTINE W2VECT !" 非線形項をベクトルに O ( EIGVEC , I WDTOR , WDDPOR , WDT , I NMKMTX , NMKMAX , I NMTLST , NMTMAX , I NMPLST , NMPMAX , I NMQLST , NMQMAX , F OATOR , OAPOR , OAT ) * * [PARAM] #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数, 波数 INCLUDE (ZLDIM) !" 格子点数, 波数 #else #include "zcdim.F" !" 格子点数, 波数 #include "zldim.F" !" 係数行列の大きさ #endif * *" [OUTPUT] REAL EIGVEC( MATDMX ) !" 固有ベクトル * *" [INPUT] REAL WDTOR ( NMDIM , 0:KDIM ) !" トロイダル Ψ REAL WDDPOR( NMDIM , 0:KDIM ) !" ポロイダル D_l Φ REAL WDT ( NMDIM , 0:KDIM ) !" 温度 Ｔ * INTEGER NMKMTX ( NMDIM, 0:KDIM, 3 ) !" 係数行列の添え字 INTEGER NMKMAX !" 係数行列の最大添え字 INTEGER NMTLST ( NMDIM ) !" MFIX 成分リスト(トロイダル) INTEGER NMPLST ( NMDIM ) !" MFIX 成分リスト(ポロイダル) INTEGER NMQLST ( NMDIM ) !" MFIX 成分リスト(温度) INTEGER NMTMAX !" NMTLST の最大添え字 INTEGER NMPMAX !" NMPLST の最大添え字 INTEGER NMQMAX !" NMQLST の最大添え字 * LOGICAL OATOR , OAPOR , OAT !" 解析対象物理量 * INTEGER K, NML * *" < 1. 固有ベクトルの計算(スペクトル) > * IF ( OATOR ) THEN DO 1100 K = 1, KDIM-1 DO 1100 NML = 1, NMTMAX EIGVEC( NMKMTX(NML,K,1) ) & = WDTOR( NMTLST(NML), K ) 1100 CONTINUE ENDIF * IF ( OAPOR )THEN DO 1200 K = 2, KDIM-2 DO 1200 NML = 1, NMPMAX EIGVEC( NMKMTX(NML,K,2) ) & = WDDPOR( NMPLST(NML), K ) 1200 CONTINUE ENDIF * IF ( OAT ) THEN DO 1300 K = 1, KDIM-1 DO 1300 NML = 1, NMQMAX EIGVEC( NMKMTX(NML,K,3) ) & = WDT( NMQLST(NML), K ) 1300 CONTINUE ENDIF * RETURN END ********************************************************************** SUBROUTINE LSTPHS !" 解析対象物理量の設定 O ( OATOR , OAPOR , OAT ) * *" [OUTPUT] LOGICAL OATOR !" トロイダル解析スイッチ LOGICAL OAPOR !" ポロイダル解析スイッチ LOGiCAL OAT !" 温度解析スイッチ * *" [INTERNAL ONCE] LOGICAL OTOR !" トロイダル解析スイッチ LOGICAL OPOR !" ポロイダル解析スイッチ LOGiCAL OTEMP !" 温度解析スイッチ DATA OTOR , OPOR , OTEMP & / .TRUE. , .TRUE. , .TRUE. / * NAMELIST /NMLPHS/ OTOR , OPOR , OTEMP * * [INTERNAL WORK] INTEGER IFPAR, JFPAR * *" < 1. 解析対象物理量の設定 > * C WRITE ( JFPAR, NMLPHS ) CALL REWNML ( IFPAR , JFPAR ) READ ( IFPAR, NMLPHS , END=1190 ) 1190 WRITE ( JFPAR, NMLPHS ) * OATOR = OTOR OAPOR = OPOR OAT = OTEMP * IF ( (.NOT.OATOR ) .OR. (.NOT. OAPOR ) .OR. (.NOT. OAT ) )THEN CALL MSGDMP( 'W', 'LSTPHS', 'ANALIZING RESTRICTED DYNAMICS' ) ENDIF * RETURN END