* PACKAGE DGDYN !" 力学 非線型項 * *" [HIS] 93/02/01(takepiro) *" 93/02/16(takepiro) *" 93/09/07(takepiro) 東西運動量の角運動量フラックスでの見積 * ********************************************************************* SUBROUTINE GRDDYN !" 格子点上での力学項 O ( GTUA , GTVA , GTWA , O GTUT , GTVT , GTWT , GTH , I GAU , GAV , GAW , GAT , C UVFACT, FLAPLA, C ALAT , DLAT , ARAD , DRAD ) * * * [PARAMETER] #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数，波数 #else #include "zcdim.F" !" 格子点数，波数 #endif * * [OUTPUT] 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 ) !" 温度鉛直移流項 ｗＴ * * [INPUT] REAL GAU ( IDIM*JDIM, 0:KDIM ) !" 西風 u(t) REAL GAV ( IDIM*JDIM, 0:KDIM ) !" 南風 v(t) REAL GAW ( IDIM*JDIM, 0:KDIM ) !" 鉛直風 w(t) REAL GAT ( IDIM*JDIM, 0:KDIM ) !" 温度 Ｔ(t) * REAL UVFACT( IDIM*JDIM ) !" u→U のファクター REAL FLAPLA( NMDIM ) !" ラプラシアンの係数 * REAL ALAT ( JDIM ) !" 緯度 REAL DLAT ( JDIM ) !" 緯度荷重 REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) * * [INTERNAL WORK] COMMON /COMWRK/ FLU, FLV, FLW REAL FLU( IDIM*JDIM, 0:KDIM ) !" 運動量フラックス u REAL FLV( IDIM*JDIM, 0:KDIM ) !" 運動量フラックス v REAL FLW( IDIM*JDIM, 0:KDIM ) !" 運動量フラックス w * REAL SINE ( IDIM*JDIM ) !" tan(μ) * INTEGER IJ, I, J, K SAVE SINE * * [INTERNAL ONCE] LOGICAL OFIRST DATA OFIRST / .TRUE. / SAVE OFIRST IF ( OFIRST ) THEN WRITE ( 6,* ) ' NONLINEAR DYNAMICS: DATE=93/09/07' OFIRST = .FALSE. IJ = 0 DO 100 J = 1, JDIM DO 100 I = 1, IDIM IJ = IJ + 1 SINE ( IJ ) = SIN ( ALAT(J) ) 100 CONTINUE ENDIF * *" < 1. ＵＡ > * DO 1100 K = 0, KDIM DO 1100 IJ = 1, IDIM*JDIM FLW( IJ,K ) = ARAD( K ) * UVFACT( IJ ) * GAU( IJ,K ) & * GAW( IJ,K ) FLV( IJ,K ) = ARAD( K ) * UVFACT( IJ ) * GAU( IJ,K ) & * GAV( IJ,K ) FLU( IJ,K ) = ARAD( K ) * UVFACT( IJ ) * GAU( IJ,K ) & * GAU( IJ,K ) 1100 CONTINUE * CALL DIVER O ( GTUA , I FLU , FLV , FLW , C UVFACT, C ARAD , DRAD ) * DO 1200 K = 0, KDIM DO 1200 IJ = 1, IDIM*JDIM GTUA( IJ,K ) = GTUA( IJ,K ) & / ( ARAD(K) * UVFACT(IJ) ) 1200 CONTINUE * DO 1300 K = 0, KDIM DO 1300 IJ = 1, IDIM*JDIM GTUA( IJ,K ) = GTUA( IJ,K ) & * UVFACT( IJ ) / ARAD( K ) 1300 CONTINUE * *" < 2. ＶＡ > * DO 2100 K = 0, KDIM DO 2100 IJ = 1, IDIM*JDIM FLW( IJ,K ) = ARAD( K ) * UVFACT( IJ ) * GAV( IJ,K ) & * GAW( IJ,K ) FLV( IJ,K ) = ARAD( K ) * UVFACT( IJ ) * GAV( IJ,K ) & * GAV( IJ,K ) FLU( IJ,K ) = ARAD( K ) * UVFACT( IJ ) * GAV( IJ,K ) & * GAU( IJ,K ) 2100 CONTINUE * CALL DIVER O ( GTVA , I FLU , FLV , FLW , C UVFACT, C ARAD , DRAD ) * DO 2200 K = 0, KDIM DO 2200 IJ = 1, IDIM*JDIM GTVA( IJ,K ) = GTVA( IJ,K ) & + ( GAU( IJ,K ) * GAU( IJ,K ) & + GAV( IJ,K ) * GAV( IJ,K ) ) * SINE( IJ ) ) 2200 CONTINUE * DO 2300 K = 0, KDIM DO 2300 IJ = 1, IDIM*JDIM GTVA( IJ,K ) = GTVA( IJ,K ) & / ( ARAD( K ) * ARAD( K ) ) 2300 CONTINUE * *" < 3. ｗＡ > * DO 3100 K = 0, KDIM DO 3100 IJ = 1, IDIM*JDIM FLW( IJ,K ) = GAW( IJ,K ) * GAW( IJ,K ) FLV( IJ,K ) = GAW( IJ,K ) * GAV( IJ,K ) FLU( IJ,K ) = GAW( IJ,K ) * GAU( IJ,K ) 3100 CONTINUE * CALL DIVER O ( GTWA , I FLU , FLV , FLW , C UVFACT, C ARAD , DRAD ) * DO 3200 K = 0, KDIM DO 3200 IJ = 1, IDIM*JDIM GTWA( IJ,K ) = GTWA( IJ,K ) & - ( GAU( IJ,K ) * GAU( IJ,K ) & + GAV( IJ,K ) * GAV( IJ,K ) ) & / ARAD( K ) 3200 CONTINUE * *" < 4. ＵＴ ，ＶＴ , ｗＴ > * DO 4100 K = 0, KDIM DO 4100 IJ = 1, IDIM*JDIM GTUT( IJ,K ) = GAU ( IJ,K ) * GAT ( IJ,K ) & * UVFACT ( IJ ) / ARAD( K ) GTVT( IJ,K ) = GAV ( IJ,K ) * GAT ( IJ,K ) & * UVFACT ( IJ ) / ARAD( K ) GTWT( IJ,K ) = GAW ( IJ,K ) * GAT ( IJ,K ) 4100 CONTINUE * *" < 5. Ｈ > * CALL RESET( GTH , IDIM*JDIM*KMAX ) * RETURN END ********************************************************************* SUBROUTINE DIVER !" 球座標系の発散(格子→格子) O ( GDDIV , I GDU , GDV , GDW , C UVFACT, C ARAD , DRAD ) * * [PARAMETER] #ifdef SYS_IBMS INCLUDE (ZCDIM) !" 格子点数，波数 #else #include "zcdim.F" !" 格子点数，波数 #endif * * [OUTPUT] REAL GDDIV( IDIM*JDIM, 0:KDIM ) !" 発散(格子) * * [INPUT] REAL GDU( IDIM*JDIM, 0:KDIM ) !" フラックス u REAL GDV( IDIM*JDIM, 0:KDIM ) !" フラックス v REAL GDW( IDIM*JDIM, 0:KDIM ) !" フラックス w * REAL UVFACT( IDIM*JDIM ) !" u→U のファクター * REAL ARAD ( 0:KDIM ) !" ｒレベル(整数) REAL DRAD ( 0:KDIM ) !" Δｒ(整数) * * [WORK] REAL GDUU( IDIM*JDIM, 0:KDIM ) !" フラックス U/r REAL GDVV( IDIM*JDIM, 0:KDIM ) !" フラックス V/r REAL WDDIV( NMDIM , KMAX ) !" 発散(スペクトル) * INTEGER IJ , K * *" < 0. 変数変換 > * DO 1000 K = 0, KDIM DO 1000 IJ = 1, IDIM*JDIM GDUU( IJ, K ) = GDU( IJ, K ) * UVFACT( IJ ) / ARAD( K ) GDVV( IJ, K ) = GDV( IJ, K ) * UVFACT( IJ ) / ARAD( K ) 1000 CONTINUE * *" < 1. スペクトルでの発散 > * CALL DIVERW O ( WDDIV , I GDUU , GDVV , GDW , F 'POS ' , C ARAD , DRAD ) * *" < 2. 格子点へ変換 > * CALL W2G O ( GDDIV , I WDDIV , I ' ' , 'POS' , KMAX ) * RETURN END