* PACKAGE DG2WD  !" 力学 格子→スペクトル
*
*"  [HIS] 90/08/31(numaguti)
*"        92/05/18(masaaki)       VOR,DIV(LMAX=0) の取扱い変更
*"        93/02/05(takepiro)
*
**********************************************************************
      SUBROUTINE TENG2W         !" 時間変化項 ( スペクトル化 )
     O         ( WTTOR , WTPOR , WTT   , 
     I           GTUA  , GTVA  , GTWA  , 
     I           GTUT  , GTVT  , GTWT  , GTH , 
     C           FLAPLA, 
     C           ARAD  , DRAD                  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)
#else
#include        "zcdim.F"
#endif
*
*   [OUTPUT]
      REAL       WTTOR ( NMDIM, KMAX      )  !" トロイダル Ψ 時間変化
      REAL       WTPOR ( NMDIM, KMAX      )  !" ポロイダル Φ 時間変化
      REAL       WTT   ( NMDIM, KMAX      )  !" 温度  Ｔ 時間変化
*
*   [INPUT]
      REAL       GTUA  ( IDIM, JDIM, KMAX )  !" 東西運動量移流項ＵＡ
      REAL       GTVA  ( IDIM, JDIM, KMAX )  !" 南北運動量移流項ＶＡ
      REAL       GTWA  ( IDIM, JDIM, KMAX )  !" 鉛直運動量移流項ｗＡ
      REAL       GTUT  ( IDIM, JDIM, KMAX )  !" 温度東西移流項  ＵＴ
      REAL       GTVT  ( IDIM, JDIM, KMAX )  !" 温度南北移流項  ＶＴ
      REAL       GTWT  ( IDIM, JDIM, KMAX )  !" 温度鉛直移流項  ｗＴ
      REAL       GTH   ( IDIM, JDIM, KMAX )  !" 温度時間変化項  Ｈ
*
      REAL       FLAPLA( NMDIM )             !" ラプラシアンの係数 -n(n+1)
*
      REAL       ARAD  ( KMAX )              !" ｒレベル(整数)
      REAL       DRAD  ( KMAX )              !" Δｒ(整数)
*
*   [INTERNAL WORK]
      COMMON    /COMWRK/
     &           WTWRK
      REAL       WTWRK  ( NMDIM, KMAX )      !" 作業領域
*
      INTEGER    K, NM
*
*"  << 1.  トロイダル Ψ  >>
*
*"      ( - ∂(ＵＡ)/∂ｙ )
*
      CALL G2W
     O         ( WTTOR ,
     I           GTUA  ,
     I           'YGRA', 'POS ',  KMAX  )
*
*"      (  ∂(ＶＡ)/∂ｘ )
*
      IF ( MMAX .GE. 1 ) THEN
         CALL G2W
     O         ( WTTOR ,
     I           GTVA  ,
     I           'XGRA', 'ADD ', KMAX  )
      ENDIF
*
*"      ( r をかける )
*
      DO 1000 K = 1, KMAX
         DO 1000 NM = 1, NMDIM
            IF ( FLAPLA( NM ) .EQ. 0 )THEN
               WTTOR( NM,K ) = 0.
            ELSE
               WTTOR( NM,K ) = - WTTOR( NM,K ) * ARAD(K) / (-FLAPLA(NM))
            ENDIF
 1000 CONTINUE 
*
*
*"  << 2.  ポロイダル Φ  >>
*
*"      (  ∂(ＶＡ)/∂ｙ )
*
      CALL G2W
     O         ( WTWRK ,
     I           GTVA  ,
     I           'YGRA', 'NEG ', KMAX )
*
      CALL RADDIV
     O         ( WTPOR ,
     I           WTWRK ,
     F           'POS ', 
     I           ARAD  , DRAD  )
*
*"      (  ∂(ＵＡ)/∂ｘ )
*
      IF ( MMAX .GE. 1 ) THEN
         CALL G2W
     O         ( WTWRK ,
     I           GTUA  ,
     I           'XGRA', 'POS ', KMAX )
*
         CALL RADDIV
     O         ( WTPOR ,
     I           WTWRK ,
     F           'ADD ', 
     I           ARAD  , DRAD  )
      ENDIF
*
*"      (  Ｌ_2 ＷＴ / r^2  )
*
         CALL G2W
     O         ( WTWRK ,
     I           GTWA  ,
     I           ' '   , 'POS ', KMAX )
*
      DO 2000 K = 1, KMAX
         DO 2000 NM = 1, NMDIM
            WTPOR( NM,K ) = WTPOR( NM,K ) 
     &                + ( -FLAPLA(NM) ) * WTWRK( NM,K ) / ARAD(K)**2
 2000 CONTINUE 
*
*"      ( r をかける )
*
      DO 2100 K = 1, KMAX
         DO 2100 NM = 1, NMDIM
            IF ( FLAPLA( NM ) .EQ. 0 )THEN
               WTPOR( NM,K ) = 0 
            ELSE
               WTPOR( NM,K ) = WTPOR( NM,K ) * ARAD(K) / ( -FLAPLA(NM) )
            ENDIF
 2100 CONTINUE 
*
*
*"   << 3. 温度 Ｔ >>
*
      CALL DIVERW
     O     (  WTT   , 
     I        GTUT  , GTVT  , GTWT , 
     F       'NEG ' ,
     C        ARAD  , DRAD           )
*
*"      ( Ｈ )
*
      CALL G2W
     O         ( WTT   ,
     I           GTH   ,
     I           '    ', 'ADD ', KMAX )
*
      RETURN
      END
*********************************************************************
      SUBROUTINE DIVERW   !" 球座標系の発散(格子→スペクトル)
     O     (  WDDIV , 
     I        GDU   , GDV  , GDW , 
     F        HFUNC , 
     C        ARAD  , DRAD         )
*
*   [PARAMETER]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT]
      REAL       WDDIV( NMDIM , KMAX )        !" 発散(スペクトル)
*
*   [INPUT]
      REAL       GDU( IDIM*JDIM, 0:KDIM )    !" フラックス U / r
      REAL       GDV( IDIM*JDIM, 0:KDIM )    !" フラックス V / r
      REAL       GDW( IDIM*JDIM, 0:KDIM )    !" フラックス w
*
      CHARACTER  HFUNC*4                     !" 符号フラグ
*
      REAL       ARAD  ( 0:KDIM )            !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )            !" Δｒ(整数)
*
*   [WORK]
      REAL       WDW  ( NMDIM , KMAX )        !" 発散(チェビシェフ)
      REAL       WDATA( NMDIM , KMAX )        !" 発散(チェビシェフ)
*
      INTEGER    NM, K
*
*"  < 1. 鉛直微分 >
*
      CALL G2W
     O         ( WDW ,
     I           GDW ,
     I           ' ',  'POSO', KMAX )
*
      CALL RADDIV
     M         ( WDATA ,
     I           WDW   ,
     F           'POS ', 
     I           ARAD  , DRAD  )
*
*"  <  2. 水平発散 > 
*
      CALL G2W
     O         ( WDATA ,
     I           GDU ,
     I           'XGRA', 'ADD ' , KMAX )
*
      CALL G2W
     O         ( WDATA ,
     I           GDV ,
     I           'YGRA', 'SUB ' , KMAX )
*
*"  < 3. 出力 >
*
      IF( HFUNC(1:1) .EQ. 'A' )THEN
         DO 2100 K = 1, KMAX
            DO 2100 NM = 1, NMDIM
               WDDIV ( NM,K ) = WDDIV( NM,K ) + WDATA( NM,K )
 2100    CONTINUE
*
      ELSE IF( HFUNC(1:1) .EQ. 'S' )THEN
         DO 2200 K = 1, KMAX
            DO 2200 NM = 1, NMDIM
               WDDIV ( NM,K ) = WDDIV( NM,K ) - WDATA( NM,K )
 2200    CONTINUE
*
      ELSE IF( HFUNC(1:1) .EQ. 'N' )THEN
         DO 2300 K = 1, KMAX
            DO 2300 NM = 1, NMDIM
               WDDIV ( NM,K ) =  - WDATA( NM,K )
 2300    CONTINUE
*
      ELSE
         DO 2400 K = 1, KMAX
            DO 2400 NM = 1, NMDIM
               WDDIV ( NM,K ) =  WDATA( NM,K )
 2400    CONTINUE
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE RADDIV  !" r^^2 d/dr r^2 (鉛直発散)の計算
     M             ( WODAT ,
     I               WIDAT ,
     F               HFUNC , 
     I               ARAD  , DRAD  )
*
*   [PARAMETER]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT]
      REAL       WODAT( NMDIM , KMAX  )      !" 鉛直発散(スペクトル)
*
*   [INPUT]
      REAL       WIDAT( NMDIM , KMAX  )      !" 鉛直フラックス(スペクトル)
*
      CHARACTER  HFUNC*4                     !" 符号フラグ
*
      REAL       ARAD  ( KMAX )              !" ｒレベル(整数)
      REAL       DRAD  ( KMAX )              !" Δｒ(整数)
*
*   [WORK]
      REAL       WDATA( NMDIM , KMAX )        !" 発散(スペクトル)
      REAL       CDATA( NMDIM , KMAX )        !" 発散(チェビシェフ)
*
      INTEGER    NM, K
*
      DO 1000 K = 1, KMAX
         DO 1000 NM = 1, NMDIM
            WDATA ( NM,K ) = WIDAT ( NM,K ) 
     &                     * ARAD( K ) * ARAD( K )
 1000 CONTINUE
*
      CALL W2C
     M         ( CDATA ,
     I           WDATA ,
     I           ' '   , 'POS ', KMAX )
*
      CALL C2W
     M         ( WDATA ,
     I           CDATA ,
     I           'DR ' , 'POS ', KMAX )
*
      DO 1100 K = 1, KMAX
         DO 1100 NM = 1, NMDIM
            WDATA ( NM,K ) = WDATA( NM,K ) /( ARAD( K ) * ARAD( K ))
 1100 CONTINUE
*
      IF( HFUNC(1:1) .EQ. 'A' )THEN
         DO 2100 K = 1, KMAX
            DO 2100 NM = 1, NMDIM
               WODAT ( NM,K ) = WODAT( NM,K ) + WDATA( NM,K )
 2100    CONTINUE
*
      ELSE IF( HFUNC(1:1) .EQ. 'S' )THEN
         DO 2200 K = 1, KMAX
            DO 2200 NM = 1, NMDIM
               WODAT ( NM,K ) = WODAT( NM,K ) - WDATA( NM,K )
 2200    CONTINUE
*
      ELSE IF( HFUNC(1:1) .EQ. 'N' )THEN
         DO 2300 K = 1, KMAX
            DO 2300 NM = 1, NMDIM
               WODAT ( NM,K ) =  - WDATA( NM,K )
 2300    CONTINUE
*
      ELSE
         DO 2400 K = 1, KMAX
            DO 2400 NM = 1, NMDIM
               WODAT ( NM,K ) =  WDATA( NM,K )
 2400    CONTINUE
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE GD2WD     !" 予報変数 ( スペクトル化 )
     O         ( WDTOR , WDPOR , WDT   , WDDPOR ,
     I           GDTOR , GDPOR , GDT   ,
     C           FLAPLA,
     C           ARAD  , DRAD                     )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)
#else
#include        "zcdim.F"
#endif
*
*   [OUTPUT]
      REAL       WDTOR ( NMDIM , KMAX     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , KMAX     )  !" ポロイダル Φ
      REAL       WDT   ( NMDIM , KMAX     )  !" 温度  Ｔ
      REAL       WDDPOR( NMDIM , KMAX     )  !" ポロイダル D_l Φ
*
*   [INPUT]
      REAL       GDTOR ( IDIM*JDIM, KMAX )   !" トロイダル Ψ
      REAL       GDPOR ( IDIM*JDIM, KMAX )   !" ポロイダル Φ
      REAL       GDT   ( IDIM*JDIM, KMAX )   !" 温度  Ｔ
*
      REAL       FLAPLA( NMDIM )             !" ラプラシアンの係数
      REAL       ARAD  ( KMAX )              !" ｒレベル(整数)
      REAL       DRAD  ( KMAX )              !" Δｒ(整数)
*
*"         < 1. 予報変数 Ψ, Φ, Ｔ の変換 >
*
      CALL G2W
     O         ( WDTOR ,
     I           GDTOR ,
     I           '    ', 'POS ', KMAX )
*
      CALL G2W
     O         ( WDPOR ,
     I           GDPOR ,
     I           '    ', 'POS ', KMAX )
*
      CALL G2W
     O         ( WDT   ,
     I           GDT   ,
     I           '    ', 'POSO', KMAX )
*
      CALL DLOPER
     M           ( WDDPOR, 
     I             WDPOR , 
     F             'POS ',
     C             FLAPLA,
     C             ARAD  , DRAD    )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE GV2POT    !" （Ｕ,Ｖ,Ｗ）格子点 → （Ψ,Φ）スペクトル
     O         ( WDTOR , WDPOR ,
     I           GDU   , GDV   , GDW   ,
     C           UVFACT, FLAPLA,
     C           ARAD  , DRAD              )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE    (ZCDIM)
#else
#include         "zcdim.F"
#endif
*
*   [OUTPUT]
      REAL       WDTOR ( NMDIM , KMAX     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , KMAX     )  !" ポロイダル Φ
*
*   [INPUT]
      REAL       GDU   ( IDIM*JDIM, KMAX )   !" 西風  ｕ
      REAL       GDV   ( IDIM*JDIM, KMAX )   !" 南風  ｖ
      REAL       GDW   ( IDIM*JDIM, KMAX )   !" 鉛直風 Ｗ
*
      REAL       UVFACT( IDIM*JDIM )         !" u→U のファクター
      REAL       FLAPLA( NMDIM )             !" ラプラシアンの係数
*
      REAL       ARAD   ( KMAX )             !" ｒレベル
      REAL       DRAD   ( KMAX )             !" Δｒ
*
*   [INTERNAL WORK]
      COMMON    /COMWRK/
     &           GDUU  , GDVV
      REAL       GDUU  ( IDIM*JDIM, KMAX )   !" 西風  Ｕ
      REAL       GDVV  ( IDIM*JDIM, KMAX )   !" 南風  Ｖ
*
      INTEGER    IJ, K, NM
*
*"         < 1. Ψ のスペクトル >
*
      DO 1100 K = 1, KMAX
         DO 1100 IJ = 1, IDIM*JDIM
            GDUU ( IJ, K ) = GDU ( IJ, K ) * UVFACT( IJ )
            GDVV ( IJ, K ) = GDV ( IJ, K ) * UVFACT( IJ )
 1100 CONTINUE
*
      CALL G2W
     O         ( WDTOR ,
     I           GDUU  ,
     I           'YGRA', 'POS ', KMAX )
*
      IF ( MMAX .GE. 1 ) THEN
         CALL G2W
     O         ( WDTOR ,
     I           GDVV  ,
     I           'XGRA', 'ADD ', KMAX )
      ENDIF
*
      DO 1200 K = 1, KMAX
         DO 1200 NM = 1, NMDIM
            IF ( FLAPLA( NM ) .LT. 0.0 )THEN
               WDTOR( NM,K )= WDTOR( NM,K ) / FLAPLA( NM )
            ENDIF
 1200 CONTINUE 
*
*"         < 2. Φ のスペクトル >
*
      CALL G2W
     O         ( WDPOR ,
     I           GDW   ,
     F           ' '   , 'POS ', KMAX )
*
      DO 2100 K = 1, KMAX
         DO 2100 NM = 1, NMDIM
            IF ( FLAPLA( NM ) .LT. 0. )THEN
               WDPOR( NM,K ) = - WDPOR( NM,K )/FLAPLA( NM )*ARAD( K )
            ENDIF
 2100 CONTINUE 
*
      RETURN
      END
