* PACKAGE DLNDYN !" 力学 線型項 * *" [HIS] 92/12/05(takepiro) トロイダルポロイダル interaction 微分表現 *" 93/02/15(takepiro) *" 93/06/06(takepiro) * *********************************************************************** SUBROUTINE LNDYN !" 力学 線型項 M ( WTTOR , WTPOR , WTT , I WDTOR , WDPOR , WDT , WDDPOR, C RAYL , TAU , PRND , C NMO , FLAPLA, C ARAD , DRAD ) * *" * #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数，波数 #else #include "zcdim.F" !" 格子点数，波数 #endif * * [MODIFY] REAL WTTOR ( NMDIM , 0:KDIM ) !" トロイダル Ψ REAL WTPOR ( NMDIM , 0:KDIM ) !" ポロイダル Φ REAL WTT ( NMDIM , 0:KDIM ) !" 温度 Ｔ * * [INPUT] REAL WDTOR ( NMDIM , 0:KDIM ) !" トロイダル Ψ REAL WDPOR ( NMDIM , 0:KDIM ) !" ポロイダル Φ REAL WDT ( NMDIM , 0:KDIM ) !" 温度 Ｔ REAL WDDPOR( NMDIM , 0:KDIM ) !" ポロイダル D_l Φ * REAL RAYL !" レイリー数 REAL TAU !" √テイラー数 REAL PRND !" プランドル数 * INTEGER NMO ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番 REAL FLAPLA( NMDIM ) !" ラプラシアンの係数 * REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) * * [INTERNAL WORK] INTEGER NM, K * *" < 1. トロイダル Ψ > * CALL GLOPER M ( WTTOR , I WDTOR , F 'ADD ', C TAU , PRND , C NMO , FLAPLA , C ARAD , DRAD ) * CALL QOPER M ( WTTOR , I WDPOR , F 'SUB' , C TAU , PRND , C NMO , FLAPLA , C ARAD , DRAD ) * DO 1000 K = 0, KDIM DO 1000 NM = 1, NMDIM IF( FLAPLA(NM) .EQ. 0. )THEN WTTOR( NM, K ) = 0. ENDIF 1000 CONTINUE * *" < 2. ポロイダル Φ > * CALL GLOPER M ( WTPOR , I WDDPOR, F 'ADD ', C TAU , PRND , C NMO , FLAPLA, C ARAD , DRAD ) * CALL QOPER M ( WTPOR , I WDTOR , F 'ADD' , C TAU , PRND , C NMO , FLAPLA , C ARAD , DRAD ) * DO 2000 K = 0, KDIM DO 2000 NM = 1, NMDIM IF( FLAPLA(NM) .NE. 0. )THEN WTPOR ( NM, K ) = WTPOR ( NM,K ) & - RAYL * PRND * WDT ( NM,K ) ELSE WTPOR( NM, K ) = 0. ENDIF 2000 CONTINUE * *" < 3. 温度場 Ｔ > * CALL DLOPER M ( WTT , I WDT , F 'ADD ', C FLAPLA, C ARAD , DRAD ) * DO 3000 K = 0, KDIM DO 3000 NM = 1, NMDIM WTT ( NM, K ) = WTT ( NM, K ) & + ( -FLAPLA(NM) ) * WDPOR( NM, K ) 3000 CONTINUE * RETURN END *********************************************************************** SUBROUTINE GLOPER !" 演算子 Pr * G_lm ( 力学:拡散+β )の計算 M ( WODATA, I WDATA , F FUNC , C TAU , PRND , C NMO , FLAPLA, C ARAD , DRAD ) * * #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数，波数 #else #include "zcdim.F" !" 格子点数，波数 #endif * * [MODIFY] REAL WODATA( NMDIM, 0:KDIM ) !" * * [INPUT] REAL WDATA ( NMDIM, 0:KDIM ) !" REAL TAU !" √テイラー数 REAL PRND !" プランドル数 CHARACTER FUNC*4 !" 動作スイッチ INTEGER NMO ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番 REAL FLAPLA( NMDIM ) !" ラプラシアンの係数 REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) * * [INTERNAL WORK] REAL WGLDAT( NMDIM, 0:KDIM ) !" INTEGER M, L, LEND, K, NM REAL AN * *" < 1. ｉｍ をかける > * CALL RESET ( WGLDAT , NMDIM*KMAX ) * DO 1000 M = 0 , MMAX, MINT LEND = MIN( LMAX, NMAX-M ) DO 1000 L = 0 , LEND AN = L + M IF( AN .NE. 0. )THEN DO 1010 K = 0, KDIM WGLDAT( NMO(1,M,L),K ) & = - REAL( M ) * TAU / ( AN*(AN+1) ) & * WDATA( NMO(2,M,L),K ) WGLDAT( NMO(2,M,L),K ) & = REAL( M ) * TAU / ( AN*(AN+1) ) & * WDATA( NMO(1,M,L),K ) 1010 CONTINUE ENDIF 1000 CONTINUE * *" < 2. D_l をかける > * CALL DLOPER M ( WGLDAT, I WDATA , F 'ADD ', C FLAPLA, C ARAD , DRAD ) * *" < 3. 出力 > * IF ( FUNC(1:1) .EQ. 'A' )THEN DO 2100 K = 0, KDIM DO 2100 NM = 1, NMDIM WODATA( NM,K ) = WODATA( NM,K ) + PRND * WGLDAT( NM,K ) 2100 CONTINUE * ELSE IF ( FUNC(1:1) .EQ. 'S' )THEN DO 2200 K = 0, KDIM DO 2200 NM = 1, NMDIM WODATA( NM,K ) = WODATA( NM,K ) - PRND * WGLDAT( NM,K ) 2200 CONTINUE * ELSE IF ( FUNC(1:1) .EQ. 'N' )THEN DO 2300 K = 0, KDIM DO 2300 NM = 1, NMDIM WODATA( NM,K ) = - PRND * WGLDAT( NM,K ) 2300 CONTINUE * ELSE DO 2400 K = 0, KDIM DO 2400 NM = 1, NMDIM WODATA( NM,K ) = PRND * WGLDAT( NM,K ) 2400 CONTINUE ENDIF * RETURN END *********************************************************************** SUBROUTINE DLOPER !" 演算子 D_l の計算 M ( WODATA, I WDATA , F FUNC , C FLAPLA, C ARAD , DRAD ) * * #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数，波数 #else #include "zcdim.F" !" 格子点数，波数 #endif * * [MODIFY] REAL WODATA( NMDIM, 0:KDIM ) !" * * [INPUT] REAL WDATA ( NMDIM, 0:KDIM ) !" CHARACTER FUNC*4 !" 動作スイッチ REAL FLAPLA( NMDIM ) !" ラプラシアンの係数 -n(n+1) REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) * * [INTERNAL WORK] COMMON /COMWRK/ WDLDAT, WDAT1, WDAT2, CDRDAT REAL WDLDAT( NMDIM, 0:KDIM ) !" REAL WDAT2 ( NMDIM, 0:KDIM ) !" REAL WDAT1 ( NMDIM, 0:KDIM ) !" REAL CDRDAT( NMDIM, 0:KDIM ) !" * INTEGER NM, K * *" < 1. 演算子 D_l の計算 > * CALL W2C M ( CDRDAT, I WDATA , I 'DR' , 'POS' , KMAX ) * CALL C2W O ( WDAT2 , I CDRDAT, I 'DR' , 'POS' , KMAX ) * CALL C2W O ( WDAT1 , I CDRDAT, I ' ' , 'POS' , KMAX ) * DO 1000 K = 0, KDIM DO 1000 NM = 1, NMDIM WDLDAT( NM,K ) = WDAT2( NM, K ) & + 2.0/ ARAD( K ) * WDAT1( NM, K ) & - ( -FLAPLA(NM) ) / ARAD(K)**2 * WDATA( NM,K ) 1000 CONTINUE * *" < 2. 出力 > * IF ( FUNC(1:1) .EQ. 'A' )THEN DO 2100 K = 0, KDIM DO 2100 NM = 1, NMDIM WODATA( NM,K ) = WODATA( NM,K ) + WDLDAT( NM,K ) 2100 CONTINUE * ELSE IF ( FUNC(1:1) .EQ. 'S' )THEN DO 2200 K = 0, KDIM DO 2200 NM = 1, NMDIM WODATA( NM,K ) = WODATA( NM,K ) - WDLDAT( NM,K ) 2200 CONTINUE * ELSE IF ( FUNC(1:1) .EQ. 'N' )THEN DO 2300 K = 0, KDIM DO 2300 NM = 1, NMDIM WODATA( NM,K ) = - WDLDAT( NM,K ) 2300 CONTINUE * ELSE DO 2400 K = 0, KDIM DO 2400 NM = 1, NMDIM WODATA( NM,K ) = WDLDAT( NM,K ) 2400 CONTINUE ENDIF * RETURN END *********************************************************************** SUBROUTINE QOPER !" 演算子 Pr * τ * Ｑ/n(n+1) の計算 M ( WODATA, I WDATA , F FUNC , C TAU , PRND , C NMO , FLAPLA , C ARAD , DRAD ) * * #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数，波数 #else #include "zcdim.F" !" 格子点数，波数 #endif * * [MODIFY] REAL WODATA( NMDIM, 0:KDIM ) !" * * [INPUT] REAL WDATA ( NMDIM, 0:KDIM ) !" * CHARACTER FUNC*4 !" 動作スイッチ * REAL TAU !" √テイラー数 REAL PRND !" プランドル数 * INTEGER NMO ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番 REAL FLAPLA( NMDIM ) !" ラプラシアンの係数 REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) * * [INTERNAL WORK] c$$$ COMMON /COMWRK/ WKGRD, WKGL2, WORK REAL WKGRD ( NMDIM, 0:KDIM ) !" REAL WKGL2 ( NMDIM, 0:KDIM ) !" REAL WORK ( NMDIM, 0:KDIM ) !" * INTEGER NM, K * *" < 1. ｋ・▽ の計算 > * CALL KGRAD O ( WKGRD , I WDATA , C NMO , C ARAD , DRAD ) * *" < 2. ｋ・▽ L_2 の計算 > * DO 2000 K = 0, KDIM DO 2000 NM = 1, NMDIM WORK( NM,K ) = - FLAPLA( NM ) * WDATA( NM,K ) 2000 CONTINUE * CALL KGRAD O ( WKGL2 , I WORK , C NMO , C ARAD , DRAD ) * *" < 3. Ｑ の計算 > * DO 3000 K = 0, KDIM DO 3000 NM = 1, NMDIM IF( FLAPLA( NM ) .NE. 0 ) THEN WORK( NM,K ) = TAU / ( -FLAPLA( NM ) ) & * ( WKGRD( NM,K ) & - 0.5 * ( - FLAPLA( NM )* WKGRD( NM,K ) & + WKGL2( NM,K ) ) ) ELSE WORK( NM,K ) = 0. ENDIF 3000 CONTINUE * *" < 4. 出力 > * IF ( FUNC(1:1) .EQ. 'A' )THEN DO 3100 K = 0, KDIM DO 3100 NM = 1, NMDIM WODATA( NM,K ) = WODATA( NM,K ) + PRND * WORK( NM,K ) 3100 CONTINUE * ELSE IF ( FUNC(1:1) .EQ. 'S' )THEN DO 3200 K = 0, KDIM DO 3200 NM = 1, NMDIM WODATA( NM,K ) = WODATA( NM,K ) - PRND * WORK( NM,K ) 3200 CONTINUE * ELSE IF ( FUNC(1:1) .EQ. 'N' )THEN DO 3300 K = 0, KDIM DO 3300 NM = 1, NMDIM WODATA( NM,K ) = - PRND * WORK( NM,K ) 3300 CONTINUE * ELSE DO 3400 K = 0, KDIM DO 3400 NM = 1, NMDIM WODATA( NM,K ) = PRND * WORK( NM,K ) 3400 CONTINUE ENDIF * RETURN END *********************************************************************** SUBROUTINE KGRAD !" 演算子 ｋ・▽ の計算 O ( WODATA, I WDATA , C NMO , C ARAD , DRAD ) * * #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数，波数 #else #include "zcdim.F" !" 格子点数，波数 #endif * * [MODIFY] REAL WODATA( NMDIM, 0:KDIM ) !" * * [INPUT] REAL WDATA ( NMDIM, 0:KDIM ) !" * INTEGER NMO ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番 REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) * * [INTERNAL WORK] COMMON /COMWRK/ WDAT1, WDAT2, WORK, GWORK REAL WDAT2 ( NMDIM, 0:KDIM ) !" REAL WDAT1 ( NMDIM, 0:KDIM ) !" REAL WORK ( NMDIM, 0:KDIM ) !" REAL GWORK ( IDIM*JDIM, 0:KDIM ) !" * INTEGER NM, K * *" < 1. μ d/dr の計算 > * CALL MLTIMU O ( WDAT1 , I WDATA , C NMO ) * CALL W2C M ( WORK , I WDAT1 , I 'DR' , 'POS' , KMAX ) * CALL C2W O ( WDAT1 , I WORK , I ' ' , 'POS' , KMAX ) * *" < 2. (1-μ^2)/r d/dμ の計算 > * CALL W2G O ( GWORK , I WDATA , I 'YGRA', 'POS ', KMAX ) * CALL G2W O ( WDAT2 , I GWORK , I ' ', 'POS ', KMAX ) * *" < 3. ｋ・▽ の計算 > * DO 1000 K = 0, KDIM DO 1000 NM = 1, NMDIM WODATA( NM,K ) = WDAT1( NM,K ) + WDAT2( NM,K ) / ARAD( K ) 1000 CONTINUE * RETURN END *********************************************************************** SUBROUTINE MLTIMU !" μW_n^m Y_n^m の計算 O ( WODATA, I WDATA , C NMO ) * * #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数，波数 #else #include "zcdim.F" !" 格子点数，波数 #endif * * [MODIFY] REAL WODATA( NMDIM, 0:KDIM ) !" * * [INPUT] REAL WDATA ( NMDIM, 0:KDIM ) !" INTEGER NMO ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番 * * [INTERNAL WORK] INTEGER NMOL ( NMDIM ) !" 計算するスペクトルのリスト INTEGER NMLSTL ( NMDIM ) !" INTEGER NMOU ( NMDIM ) !" INTEGER NMLSTU ( NMDIM ) !" * INTEGER NML, NMU * REAL COEFLO( NMDIM ) !" 演算子の係数 REAL COEFUP( NMDIM ) !" 演算子の係数 * SAVE NMOL, NMOU, NMLSTL, NMLSTU, NML, NMU, COEFLO, COEFUP * INTEGER NM, L, M, N, LEND, K, I, IE * LOGICAL OFIRST DATA OFIRST / .TRUE. / * *" < 1. 係数, リスト作成 > * IF ( OFIRST )THEN OFIRST = .FALSE. NML = 0 NMU = 0 * IF ( MMAX .EQ. 0 )THEN !" 東西平均 or not IE = 1 ELSE IE = 2 ENDIF * DO 1000 M = 0 , MMAX, MINT LEND = MIN( LMAX, NMAX-M ) DO 1000 L = 0 , LEND DO 1000 I = 1, IE N = L + M * IF( L-1 .GE. 0 )THEN NML = NML + 1 NMOL ( NML ) = NMO( I,M,L ) NMLSTL( NML ) = NMO( I,M,L-1 ) COEFLO( NML ) = SQRT( REAL( L*(N+M) ) & / REAL( (2*N-1)*(2*N+1) ) ) ENDIF * IF( L+1 .LE. LEND )THEN NMU = NMU + 1 NMOU ( NMU ) = NMO( I,M,L ) NMLSTU( NMU ) = NMO( I,M,L+1 ) COEFUP( NMU ) = SQRT( REAL( (L+1)*(N+M+1) ) & / REAL( (2*N+1)*(2*N+3) ) ) ENDIF 1000 CONTINUE ENDIF * *" < 2. スペクトル係数の計算 > * CALL RESET( WODATA, NMDIM*KMAX ) * DO 2000 K = 0, KDIM DO 2000 NM = 1, NML WODATA( NMOL(NM), K ) = WODATA( NMOL(NM), K ) & + COEFLO( NM ) * WDATA( NMLSTL( NM ), K ) 2000 CONTINUE * DO 2100 K = 0, KDIM DO 2100 NM = 1, NMU WODATA( NMOU(NM), K ) = WODATA( NMOU(NM), K ) & + COEFUP( NM ) * WDATA( NMLSTU( NM ), K ) 2100 CONTINUE * RETURN END