* PACKAGE DGDYN  !" 力学 非線型項 
*
*"  [HIS]  93/02/01(takepiro)
*"         93/09/16(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       TANG ( IDIM*JDIM )          !" tan(μ)
*
      INTEGER    IJ, I, J, K
      SAVE       TANG
*
*   [INTERNAL ONCE]
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE       OFIRST
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' NONLINEAR DYNAMICS:  DATE=93/09/16'
         OFIRST = .FALSE.
         IJ = 0
         DO 100 J = 1, JDIM
            DO 100 I = 1, IDIM
               IJ = IJ + 1 
               TANG ( IJ ) = TAN ( ALAT(J) )
 100     CONTINUE 
      ENDIF
*
*"         < 1.  ＵＡ  > 
*
      DO 1100 K = 0, KDIM
         DO 1100 IJ = 1, IDIM*JDIM
            FLW( IJ,K ) = GAU( IJ,K ) * GAW( IJ,K ) 
            FLV( IJ,K ) = GAU( IJ,K ) * GAV( IJ,K ) 
            FLU( IJ,K ) = GAU( IJ,K ) * GAU( IJ,K ) 
 1100 CONTINUE
*
      CALL DIVERT
     O     (  GTUA  , 
     I        FLU   , FLV  , FLW , 
     C        UVFACT, TANG , 
     C        ARAD  , DRAD         )
*
      DO 1200 K = 0, KDIM
         DO 1200 IJ = 1, IDIM*JDIM
            GTUA( IJ,K ) = GTUA( IJ,K )
     &                   + (   GAW( IJ,K ) * GAU( IJ,K )
     &                       - GAV( IJ,K ) * GAU( IJ,K ) * TANG( IJ ) )
     &                   / ARAD( K ) 
 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 ) = GAV( IJ,K ) * GAW( IJ,K ) 
            FLV( IJ,K ) = GAV( IJ,K ) * GAV( IJ,K ) 
            FLU( IJ,K ) = GAV( IJ,K ) * GAU( IJ,K ) 
 2100 CONTINUE
*
      CALL DIVERT
     O     (  GTVA  , 
     I        FLU   , FLV  , FLW , 
     C        UVFACT, TANG , 
     C        ARAD  , DRAD         )
*
      DO 2200 K = 0, KDIM
         DO 2200 IJ = 1, IDIM*JDIM
            GTVA( IJ,K ) = GTVA( IJ,K )
     &                   + (   GAW( IJ,K ) * GAV( IJ,K )
     &                       + GAU( IJ,K ) * GAU( IJ,K ) * TANG( IJ ) )
     &                   / ARAD( K ) 
 2200 CONTINUE
*
      DO 2300 K = 0, KDIM
         DO 2300 IJ = 1, IDIM*JDIM
            GTVA( IJ,K ) = GTVA( IJ,K )
     &                   * UVFACT( IJ ) / 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 DIVERG
     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 DIVERG    !" 球座標系の発散(格子→格子)
     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
*********************************************************************
      SUBROUTINE DIVERT   !" 球座標系の発散(格子→格子)
     O     (  GDDIV , 
     I        GDU   , GDV  , GDW , 
     C        UVFACT, TANG ,
     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       TANG ( IDIM*JDIM )          !" tan(μ)
*
      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       GDWW( IDIM*JDIM, 0:KDIM )    !" フラックス W/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 )
     &                                   * UVFACT( IJ ) / ARAD( K )
            GDVV( IJ, K ) = GDV( IJ, K ) * UVFACT( IJ ) / ARAD( K )
     &                                   * UVFACT( IJ ) / ARAD( K )
            GDWW( IJ, K ) = GDW( IJ, K ) * UVFACT( IJ ) / ARAD( K )
 1000 CONTINUE 
*
*"   < 1. スペクトルでの発散 >
*
      CALL DIVERW
     O     (  WDDIV , 
     I        GDUU  , GDVV , GDWW , 
     F       'POS ' ,
     C        ARAD  , DRAD         )
*
*"   < 2. 格子点へ変換 >
*
      CALL W2G
     O     (  GDDIV , 
     I        WDDIV , 
     I        ' '   , 'POS' , KMAX )
*
*"   < 3. 補正 >
*
      DO 3000 K = 0, KDIM
         DO 3000 IJ = 1, IDIM*JDIM
            GDDIV( IJ, K ) 
     &         = GDDIV( IJ, K ) / UVFACT( IJ ) * ARAD( K )
     &            + GDV( IJ, K ) * TANG( IJ ) / ARAD( K )
     &            + GDW( IJ, K ) / ARAD( K ) 
 3000 CONTINUE 
*
      RETURN
      END
