* PACKAGE USPHXZ   !" 球面調和関数変換 (XZ2次元)
*
*"  [HIS]  90/08/31(numaguti)
*"         92/05/18(masaaki)       d/dy(n=0) の追加
*
**********************************************************************
      SUBROUTINE SPW2G      !" 球面調和関数変換 (スペクトル→格子)
     M         ( GDATA ,
     I           WDATA ,
     C           PNM   , NMO   , TRIGS , IFAX ,
     F           HGRAD , HFUNC ,
     D           IMAX  , JMAX  , KMAX  , IDIM  , JDIM  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF, 
     W           ZDATA , WORK                           )
*
*   [PARAM] 
      INTEGER    IMAX
      INTEGER    JMAX
      INTEGER    KMAX
      INTEGER    IDIM
      INTEGER    JDIM
      INTEGER    LMAX
      INTEGER    MMAX
      INTEGER    NMAX
      INTEGER    MINT
      INTEGER    NMDIM
      INTEGER    JMXHF
*
*   [MODIFY] 
      REAL       GDATA ( IDIM*JDIM, KMAX )    !" 格子点データ
*
*   [INPUT] 
      REAL       WDATA ( NMDIM, KMAX     )    !" スペクトルデータ
*
      REAL       PNM   ( * )                  !" ルジャンドル関数
      INTEGER    NMO   ( 2, 0:MMAX , 0:LMAX ) !" スペクトルの添字順番
      REAL       TRIGS ( * )                  !" 三角関数表
      INTEGER    IFAX  ( * )                  !" IMAX の因数分解
*
      CHARACTER  HGRAD*4                      !" 微分フラグ
      CHARACTER  HFUNC*4                      !" 符号フラグ
*
*   [WORK] 
      REAL       ZDATA ( IDIM*JDIM, KMAX )    !" 東西スペクトル
      REAL       WORK  ( IDIM*JDIM, KMAX )    !" ワーク
*
*   [INTERNAL WORK] 
      LOGICAL    LOFFS                        !" オフセットフラグ
      INTEGER    KMAXD
      PARAMETER (KMAXD=100)
      REAL       DOFFS ( KMAXD )              !" オフセット値
      INTEGER    IJ, K, M
*
      IF ( JMAX .NE. 1 ) THEN
         CALL MSGDMP( 'E','SPW2G','THIS ROUTINE IS FOR X-Z 2 DIM.')
      ENDIF
*
*"         < 1. LOFFS フラグ >
*
      IF ( HFUNC(4:4) .EQ. 'O' ) THEN
         LOFFS = .TRUE.
         IF ( KMAXD .LT. KMAX ) THEN
            CALL MSGDMP( 'E','SPW2G','INSUFFICIENT WORK STRAGE')
          ENDIF
      ELSE
         LOFFS = .FALSE.
      ENDIF
*
*"         < 2. 詰め替え >
*
      CALL RESET ( ZDATA, IDIM*JDIM*KMAX )
      CALL RESET ( WORK , IDIM*JDIM*KMAX )
      DO 2100 K = 1, KMAX
         DO 2110 M = 0, MMAX, MINT
            ZDATA( M*2+1, K ) = WDATA( NMO(1,M,0), K )
            ZDATA( M*2+2, K ) = WDATA( NMO(2,M,0), K )
 2110    CONTINUE
*
         IF ( LOFFS ) THEN
            ZDATA( 1, K ) = 0.
            DOFFS( K )    = WDATA( NMO(1,0,0), K )
         ENDIF
 2100 CONTINUE
*
      IF ( HGRAD(1:1) .EQ. 'X' ) THEN
*
*"         < 3. ｘ微分 >
*
         CALL GRADX
     M         ( ZDATA,
     D           IDIM , JDIM , KMAX  , MMAX , MINT ,
     W           WORK                                 )
         LOFFS = .FALSE.
      ENDIF
*
      IF ( HGRAD(1:1) .NE. 'Y' ) THEN
*
*"         < 4. zonal wave → 格子点 >
*
         CALL FFT99X
     O         ( WORK  ,
     M           ZDATA ,
     C           TRIGS , IFAX  ,
     C           1     , IDIM  , IMAX  , JDIM*KMAX , 1    )
*
      ELSE
        DO  4200 K = 1, KMAX
            DO 4210 IJ = 1, IDIM*JDIM
               WORK ( IJ,K ) = WDATA( NMO(1,0,0),K )
 4210       CONTINUE
 4200    CONTINUE
      ENDIF
*
      IF ( LOFFS ) THEN
         DO 4100 K = 1, KMAX
#ifdef SYS_SX3
*vdir noloopchg
#endif
            DO 4110 IJ = 1, IDIM*JDIM
               WORK ( IJ,K ) = WORK ( IJ,K ) + DOFFS( K )
 4110       CONTINUE
 4100    CONTINUE
      ENDIF
*
*"         < 5. 出力データ >
*
      IF      ( HFUNC(1:1) .EQ. 'A' ) THEN
*
*"                               ( add )
        DO 5000 K = 1, KMAX
           DO 5000 IJ = 1, IDIM*JDIM
              GDATA( IJ,K ) = GDATA( IJ,K ) + WORK( IJ,K )
 5000   CONTINUE
*
      ELSE IF ( HFUNC(1:1) .EQ. 'S' ) THEN
*
*"                               ( sub )
        DO 5100 K = 1, KMAX
           DO 5100 IJ = 1, IDIM*JDIM
              GDATA( IJ,K ) = GDATA( IJ,K ) - WORK( IJ,K )
 5100   CONTINUE
*
      ELSE IF ( HFUNC(1:1) .EQ. 'N' ) THEN
*
*"                               ( negative )
        DO 5200 K = 1, KMAX
           DO 5200 IJ = 1, IDIM*JDIM
              GDATA( IJ,K ) = - WORK( IJ,K )
 5200   CONTINUE
*
      ELSE
*"                               ( positive )
        DO 5300 K = 1, KMAX
           DO 5300 IJ = 1, IDIM*JDIM
              GDATA( IJ,K ) = WORK( IJ,K )
 5300   CONTINUE
      ENDIF
*
      RETURN
      END
**********************************************************************
      SUBROUTINE SPG2W     !" 球面調和関数変換 (格子→スペクトル)
     M         ( WDATA ,
     I           GDATA ,
     C           PNM   , NMO   , TRIGS , IFAX , GW   ,
     F           HGRAD , HFUNC ,
     D           IMAX  , JMAX  , KMAX  , IDIM  , JDIM  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF, 
     W           ZDATA , WORK                           )
*
*   [PARAM] 
      INTEGER    IMAX
      INTEGER    JMAX
      INTEGER    KMAX
      INTEGER    IDIM
      INTEGER    JDIM
      INTEGER    LMAX
      INTEGER    MMAX
      INTEGER    NMAX
      INTEGER    MINT
      INTEGER    NMDIM
      INTEGER    JMXHF
*
*   [MODIFY] 
      REAL       WDATA ( NMDIM, KMAX     )    !" スペクトルデータ
*
*   [INPUT] 
      REAL       GDATA ( IDIM*JDIM, KMAX )    !" 格子点データ
*
      REAL       PNM   ( * )                  !" ルジャンドル関数
      INTEGER    NMO   ( 2, 0:MMAX , 0:LMAX ) !" スペクトルの添字順番
      REAL       GW    ( JDIM )               !" ガウス荷重
      REAL       TRIGS ( * )                  !" 三角関数表
      INTEGER    IFAX  ( * )                  !" IMAX の因数分解
*
      CHARACTER  HGRAD*4                      !" 微分フラグ
      CHARACTER  HFUNC*4                      !" 符号フラグ
*
*   [WORK] 
      REAL       ZDATA ( IDIM*JDIM, KMAX )    !" 東西スペクトル
      REAL       WORK  ( IDIM*JDIM, KMAX )    !" ワーク
*
*   [INTERNAL WORK] 
      LOGICAL    LOFFS                        !" オフセットフラグ
      INTEGER    KMAXD
      PARAMETER (KMAXD=100)
      REAL       DOFFS ( KMAXD )              !" オフセット値
      INTEGER    IJ, K, M, I
      REAL       WORKAV
*
      IF ( JMAX .NE. 1 ) THEN
        CALL MSGDMP( 'E','SPG2W','THIS ROUTINE IS FOR X-Z 2 DIM.')
      ENDIF
*
*"         < 1. LOFFS フラグ >
*
      IF ( HFUNC(4:4) .EQ. 'O' ) THEN
         LOFFS = .TRUE.
         IF ( KMAXD .LT. KMAX ) THEN
            CALL MSGDMP( 'E','SPG2W','INSUFFICIENT WORK STRAGE')
         ENDIF
      ELSE
         LOFFS = .FALSE.
      ENDIF
*
*"         < 2. 入力複製 >
*
      IF ( LOFFS ) THEN
         DO 2000 K = 1, KMAX
            DOFFS ( K ) = GDATA( 1,K )
#ifdef SYS_SX3
*vdir noloopchg
#endif            
            DO 2010 IJ = 1, IDIM*JDIM
               WORK ( IJ,K ) = GDATA( IJ,K ) - DOFFS( K )
 2010       CONTINUE
 2000    CONTINUE
      ELSE
         CALL COPY  ( WORK , GDATA , IDIM*JDIM*KMAX )
      ENDIF
*
      IF ( HGRAD(1:1) .EQ. 'Y' ) THEN
*
         CALL RESET ( ZDATA, IDIM*JDIM*KMAX )
*
         DO 2100 K = 1, KMAX
            WORKAV = 0.
            DO 2110 I = 1, IMAX
               WORKAV = WORKAV + WORK( I,K )
 2110       CONTINUE
            ZDATA( NMO(1,0,0),K ) = - WORKAV
 2100    CONTINUE
*
         LOFFS = .FALSE.
*
      ELSE
*
*"         < 3. 格子 → zonal wave >
*
         CALL FFT99X
     M         ( WORK  ,
     O           ZDATA ,
     C           TRIGS , IFAX  ,
     C           1     , IDIM  , IMAX  , JDIM*KMAX , 0    )
*
         IF ( HGRAD(1:1) .EQ. 'X' ) THEN
*
*"                 < 4. ｘ微分 >
*
            CALL GRADX
     M         ( ZDATA,
     D           IDIM , JDIM , KMAX , MMAX , MINT  ,
     W           WORK )
*
            LOFFS = .FALSE.
         ENDIF
*
      ENDIF
*
*"         < 5. 詰め替え >
*
      IF      ( HFUNC(1:1) .EQ. 'A' ) THEN
*
*"                               ( add )
         DO 5100 K = 1, KMAX
            DO 5110 M = 0, MMAX, MINT
               WDATA( NMO(1,M,0), K )
     &       = WDATA( NMO(1,M,0), K ) + ZDATA( M*2+1, K ) * GW( 1 )
               WDATA( NMO(2,M,0), K )
     &       = WDATA( NMO(2,M,0), K ) + ZDATA( M*2+2, K ) * GW( 1 )
 5110       CONTINUE
 5100    CONTINUE
*
      ELSE IF ( HFUNC(1:1) .EQ. 'S' ) THEN
*
*"                               ( sub )
         DO 5200 K = 1, KMAX
            DO 5210 M = 0, MMAX, MINT
               WDATA( NMO(1,M,0), K )
     &       = WDATA( NMO(1,M,0), K ) - ZDATA( M*2+1, K ) * GW( 1 )
               WDATA( NMO(2,M,0), K )
     &       = WDATA( NMO(2,M,0), K ) - ZDATA( M*2+2, K ) * GW( 1 )
 5210       CONTINUE
 5200    CONTINUE
*
      ELSE IF ( HFUNC(1:1) .EQ. 'N' ) THEN
*
*"                               ( negative )
         DO 5300 K = 1, KMAX
            DO 5310 M = 0, MMAX, MINT
               WDATA( NMO(1,M,0), K ) = - ZDATA( M*2+1, K ) * GW( 1 )
               WDATA( NMO(2,M,0), K ) = - ZDATA( M*2+2, K ) * GW( 1 )
 5310       CONTINUE
 5300    CONTINUE
*
      ELSE
*"                               ( positive )
         DO 5400 K = 1, KMAX
            DO 5410 M = 0, MMAX, MINT
               WDATA( NMO(1,M,0), K ) = ZDATA( M*2+1, K ) * GW( 1 )
               WDATA( NMO(2,M,0), K ) = ZDATA( M*2+2, K ) * GW( 1 )
 5410       CONTINUE
 5400    CONTINUE
      ENDIF
*
      IF ( LOFFS ) THEN
         IF (      ( HFUNC(1:1) .EQ. 'N' )
     &        .OR. ( HFUNC(1:1) .EQ. 'S' )    ) THEN
*
            DO 6100 K = 1, KMAX
               WDATA( NMO(1,0,0),K )
     &       = WDATA( NMO(1,0,0),K ) - DOFFS( K )
 6100       CONTINUE
         ELSE
            DO 6200 K = 1, KMAX
               WDATA( NMO(1,0,0),K )
     &       = WDATA( NMO(1,0,0),K ) + DOFFS( K )
 6200       CONTINUE
*
         ENDIF
      ENDIF
*
      RETURN
      END
************************************************************************
      SUBROUTINE GRADX     !" ｘ微分のスペクトル成分
     M         ( ZDATA ,
     D           IDIM  , JDIM  , KMAX  , MMAX  , MINT ,
     W           ZDW    )
*
*   [PARAM] 
      INTEGER    IDIM
      INTEGER    JDIM
      INTEGER    KMAX
      INTEGER    MMAX
      INTEGER    MINT
*
*   [MODIFY] 
      REAL       ZDATA ( 2, 0:IDIM/2-1, JDIM, KMAX ) !" データ
*
*   [WORK] 
      REAL       ZDW   ( 2, 0:IDIM/2-1, JDIM, KMAX ) !" ワーク
*
*   [INTERNAL WORK] 
      INTEGER    M, J, K, MM
*
*"         < 1. 作業領域 >
*
      CALL COPY ( ZDW , ZDATA , IDIM*JDIM*KMAX )
*
*"         < 2. ｉｍ をかける >
*
      DO 2100 M = 0, MMAX, MINT
         MM = M / MINT
         DO 2100 K = 1, KMAX
            DO 2100 J = 1, JDIM
               ZDATA( 1, MM, J, K ) = - REAL( M ) * ZDW  ( 2, MM, J, K )
               ZDATA( 2, MM, J, K ) =   REAL( M ) * ZDW  ( 1, MM, J, K )
 2100 CONTINUE
*
      RETURN
      END
