* PACKAGE DGENG  !" 力学 スペクトル→格子点値
*
************************************************************************
      SUBROUTINE GENGD          !" 新格子点値の生成
     O         ( GDU   , GDV   , GDW   , GDT   , 
     O           GDTOR , GDPOR ,
     I           WDTOR , WDPOR , WDT  , 
     C           FLAPLA, UVFACT,
     C           ARAD  , DRAD                    )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)
#else
#include        "zcdim.F"
#endif
*
*   [OUTPUT]
      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 )   !" ポロイダル Φ
*
*   [INPUT]
      REAL       WDTOR ( NMDIM , KMAX     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , KMAX     )  !" ポロイダル Φ
      REAL       WDT   ( NMDIM , KMAX     )  !" 温度  Ｔ
*
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       UVFACT( IDIM*JDIM )         !" u→U のファクター
      REAL       ARAD   ( KMAX )             !" ｒレベル
      REAL       DRAD   ( KMAX )             !" Δｒ
*
*
*"         < 1. Ψ 格子点値 >
*
      CALL W2G
     O         ( GDTOR ,
     I           WDTOR,
     I           '    ', 'POS', KMAX )
*
*"         < 2. Φ 格子点値 >
*
      CALL W2G
     O         ( GDPOR ,
     I           WDPOR,
     I           '    ', 'POS', KMAX )
*
*"         < 2. Ｔ 格子点値 >
*
      CALL W2G
     O         ( GDT   ,
     I           WDT   ,
     I           '    ', 'POSO', KMAX )
*
*"         < 3. Ｕ, Ｖ, Ｗ の格子点値 >
*
      CALL POTVEL
     O         ( GDU   , GDV   , GDW   , 
     I           WDTOR , WDPOR , 
     C           FLAPLA, UVFACT        ,
     C           ARAD  , DRAD               )
*
      RETURN
      END
*
************************************************************************
      SUBROUTINE POTVEL          !" 速度場の計算
     O         ( GDU   , GDV   , GDW   , 
     I           WDTOR , WDPOR , 
     C           FLAPLA, UVFACT        ,
     C           ARAD  , DRAD            )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)
#else
#include        "zcdim.F"
#endif
*
*   [OUTPUT]
      REAL       GDU   ( IDIM*JDIM, KMAX )   !" 西風  u
      REAL       GDV   ( IDIM*JDIM, KMAX )   !" 南風  v
      REAL       GDW   ( IDIM*JDIM, KMAX )   !" 鉛直風 w
*
*   [INPUT]
      REAL       WDTOR ( NMDIM , KMAX     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , KMAX     )  !" ポロイダル Φ
*
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       UVFACT( IDIM*JDIM )         !" u→U のファクター
      REAL       ARAD   ( KMAX )             !" ｒレベル
      REAL       DRAD   ( KMAX )             !" Δｒ
*
*   [INTERNAL WORK]
      COMMON    /COMWRK/
     &           WDPHI , CDPOR
      REAL       CDPOR ( NMDIM , KMAX  )  !" ポロイダル Φ
      REAL       WDPHI ( NMDIM , KMAX  )  !" rD(rΦ) スペクトル(鉛直)格子点
*
      INTEGER    IJ, K, NM
*
*"         < 1.1 トロイダルポテンシャル >
*
*"         < 1.2  ポロイダルポテンシャルのｒ微分 >
*
      DO 3000 K = 1, KMAX
         DO 3000 NM = 1, NMDIM
            WDPHI ( NM, K ) = WDPOR ( NM, K ) / ARAD( K )
 3000 CONTINUE 
*
      CALL W2C
     M         ( CDPOR ,
     I           WDPOR ,
     I           ' '   , 'POS ', KMAX )
*
      CALL C2W
     M         ( WDPHI ,
     I           CDPOR ,
     I           'DR  ', 'ADD ', KMAX )
*
*"         < 2.  Ｕ の計算 >
*
      CALL W2G
     O         ( GDU   ,
     I           WDTOR ,
     I           'YGRA', 'POS ', KMAX )
*
      IF ( MMAX .GE. 1 ) THEN
         CALL W2G
     O         ( GDU   ,
     I           WDPHI ,
     I           'XGRA', 'ADD ', KMAX )
      ENDIF
*
*"         < 3. Ｖ の計算 >
*
      CALL W2G
     O         ( GDV   ,
     I           WDPHI ,
     I           'YGRA', 'POS ', KMAX )
*
      IF ( MMAX .GE. 1 ) THEN
         CALL W2G
     O         ( GDV   ,
     I           WDTOR ,
     I           'XGRA', 'SUB ', KMAX )
      ENDIF
*
      DO 3100 K = 1, KMAX
         DO 3100 IJ = 1, IDIM*JDIM
            GDU ( IJ,K ) = GDU ( IJ,K ) / UVFACT( IJ )
            GDV ( IJ,K ) = GDV ( IJ,K ) / UVFACT( IJ )
 3100 CONTINUE
*
*"         < 4. w の計算 >
*
      DO 3200 K = 1, KMAX
         DO 3200 NM = 1, NMDIM
            WDPHI ( NM, K ) = - FLAPLA( NM ) * WDPOR( NM, K )
 3200 CONTINUE 
*
      CALL W2G
     O         ( GDW   ,
     I           WDPHI ,
     I           '    ', 'POS ', KMAX )
*
      DO 3300 K = 1, KMAX
         DO 3300 IJ = 1, IDIM*JDIM
            GDW ( IJ,K ) = GDW ( IJ,K ) / ARAD( K )
 3300 CONTINUE
*
      RETURN
      END
