* PACKAGE USPHE   !"  球面調和関数変換
*
*"  [HIS] 90/08/31(numaguti)
*
*
**********************************************************************
      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   ( NMDIM, JMXHF )       !" ルジャンドル関数
      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    LDPNM                        !" ｙ微分フラグ
      LOGICAL    LOFFS                        !" オフセットフラグ
      INTEGER    KMAXD
      PARAMETER (KMAXD=100)
      REAL       DOFFS ( KMAXD )              !" オフセット値
      INTEGER    IJ, K
*
*"         < 1. LOFFS, LDPNM : フラグ >
*
      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
*
      IF ( HGRAD(1:1) .EQ. 'Y' ) THEN
         LDPNM = .TRUE.
         LOFFS = .FALSE.
      ELSE
         LDPNM = .FALSE.
      ENDIF
*
*"         < 2. スペクトル → zonal wave >
*
      IF ( LOFFS ) THEN
         DO 2000 K = 1, KMAX
            DOFFS( K ) = WDATA( NMO(1,0,0), K )
            WDATA( NMO(1,0,0), K ) = 0.  
 2000    CONTINUE
      ENDIF
*
      CALL SPW2Z
     O         ( ZDATA ,
     I           WDATA ,
     C           PNM   , NMO   ,
     F           LDPNM ,
     D           JMAX  , KMAX  , IDIM  , JDIM  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF ,
     W           WORK                                           )
*
      IF ( LOFFS ) THEN
         DO 2100 K = 1, KMAX
            WDATA( NMO(1,0,0), K ) = DOFFS( K )
 2100    CONTINUE
      ENDIF
*
      IF ( HGRAD(1:1) .EQ. 'X' ) THEN
*
*"         < 3. ｘ微分 >
*
         CALL GRADX
     M         ( ZDATA,
     D           IDIM , JDIM , KMAX  , MMAX , MINT ,
     W           WORK                                 )
         LOFFS = .FALSE.
      ENDIF
*
*"         < 4. zonal wave → 格子点 >
*
      CALL FFT99X
     O         ( WORK  ,
     M           ZDATA ,
     C           TRIGS , IFAX  ,
     C           1     , IDIM  , IMAX  , JDIM*KMAX , 1    )
*
      IF ( LOFFS ) THEN
         DO 4100 K = 1, KMAX
#ifdef SYS_SX3
*vdir noloopchg
#endif
            DO 4100 IJ = 1, IDIM*JDIM
               WORK ( IJ,K ) = WORK ( IJ,K ) + DOFFS( K )
 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   ( NMDIM, JMXHF )       !" ルジャンドル関数
      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    LDPNM                        !" ｙ微分フラグ
      LOGICAL    LOFFS                        !" オフセットフラグ
      INTEGER    KMAXD
      PARAMETER (KMAXD=100)
      REAL       DOFFS ( KMAXD )              !" オフセット値
      INTEGER    IJ, K
*
*"         < 1. LOFFS, LDPNM : フラグ >
*
      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
*
      IF ( HGRAD(1:1) .EQ. 'Y' ) THEN
         LDPNM = .TRUE.
         LOFFS = .FALSE.
      ELSE
         LDPNM = .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
*
*"         < 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
*
*"         < 5. zonal wave → スペクトル >
*
      CALL SPZ2W
     M         ( WDATA ,
     I           ZDATA ,
     C           PNM   , NMO   , GW    ,
     F           LDPNM , HFUNC ,
     D           JMAX  , KMAX  , IDIM  , JDIM  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF ,
     W           WORK                                           )
*
      IF ( LOFFS ) THEN
         IF (      ( HFUNC(1:1) .EQ. 'N' )
     &        .OR. ( HFUNC(1:1) .EQ. 'S' )    ) THEN
            DO 5100 K = 1, KMAX
               WDATA( NMO(1,0,0),K ) = WDATA( NMO(1,0,0),K )
     &                               - DOFFS( K )
 5100       CONTINUE
         ELSE
            DO 5200 K = 1, KMAX
               WDATA( NMO(1,0,0),K ) = WDATA( NMO(1,0,0),K )
     &                               + DOFFS( K )
 5200       CONTINUE
         ENDIF
      ENDIF
*
      RETURN
      END
**********************************************************************
      SUBROUTINE SPW2Z  !" ルジャンドル変換 (スペクトル→格子)
     O         ( ZDATA   ,
     I           WDATA ,
     C           PNM   , NMO   ,
     F           LDPNM ,
     D           JMAX  , KMAX  , IDIM  , JDIM  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF ,
     W           ZDW                                            )
*
*   [PARAM] 
      INTEGER    JMAX
      INTEGER    KMAX
      INTEGER    IDIM
      INTEGER    JDIM
      INTEGER    LMAX
      INTEGER    MMAX
      INTEGER    NMAX
      INTEGER    MINT
      INTEGER    NMDIM
      INTEGER    JMXHF
*
*   [OUTPUT] 
      REAL       ZDATA ( IDIM, JDIM, KMAX )   !" 東西スペクトル
*
*   [INPUT] 
      REAL       WDATA ( NMDIM, KMAX      )   !" スペクトル
*
      REAL       PNM   ( NMDIM, JMXHF )       !" ルジャンドル関数
      INTEGER    NMO   ( 2, 0:MMAX , 0:LMAX ) !" スペクトルの添字順番
*
      LOGICAL    LDPNM                        !" ｙ微分フラグ
*
*   [WORK] 
      REAL       ZDW   ( IDIM, JDIM, KMAX )   !" ワーク
*
*   [INTERNAL ONCE] 
      INTEGER    NMDIMD
      PARAMETER (NMDIMD=12000)
      INTEGER    MLIST ( NMDIMD )             !" リストベクトル
      INTEGER    JLIST ( NMDIMD )             !" リストベクトル
      LOGICAL    OLSET
      DATA       OLSET / .FALSE. /
      SAVE       MLIST, JLIST, OLSET
*
*   [INTERNAL WORK] 
      INTEGER    L, M, NM, K
      INTEGER    LEND, MM, IM
      INTEGER    J, JP, JE, JO, JN, JS
*
*
*"         < 0. リストベクトル作成 >
*
*   [ONCE] 
      IF ( .NOT. OLSET ) THEN
         OLSET = .TRUE.
*
         IF ( NMDIMD .LT. NMDIM ) THEN
            CALL MSGDMP( 'E','SPW2Z','INSUFFICIENT WORK STRAGE')
         ENDIF
*
         DO 200 M = 0, MMAX, MINT
            LEND = MIN ( LMAX, NMAX-M )
            MM = M / MINT
*
            DO 100 L = 0, LEND
               MLIST (  NMO ( 1,M,L )  ) = 2 * MM + 1
               MLIST (  NMO ( 2,M,L )  ) = 2 * MM + 2
*
               IF ( MOD( L,2 ) .EQ. 0 ) THEN
                  JLIST (  NMO ( 1,M,L )  ) = 0
                  JLIST (  NMO ( 2,M,L )  ) = 0
               ELSE
                  JLIST (  NMO ( 1,M,L )  ) = JMAX / 2
                  JLIST (  NMO ( 2,M,L )  ) = JMAX / 2
               ENDIF
*
  100       CONTINUE
  200    CONTINUE
*
      ENDIF
*
*"         < 1. 作業領域リセット >
*
      CALL RESET ( ZDW  , IDIM*JDIM*KMAX )
      CALL RESET ( ZDATA, IDIM*JDIM*KMAX )
*
*"         < 2.  ルジャンドル変換 >
*
      DO 2300 NM = 1 , NMDIM
         IM = MLIST( NM )
*
         DO 2200 K = 1 , KMAX
#ifdef SYS_HITAC
*VOPTION INDEP(ZDW)
#endif
#ifdef SYS_SX3
*vdir nodep(ZDW)
#endif
            DO 2100 J = 1 , JMAX/2
               JP = JLIST( NM ) + J
*
               ZDW( IM,JP,K   ) = ZDW( IM,JP,K   )
     &                        +  PNM  ( NM,J ) * WDATA( NM,K   )
 2100       CONTINUE
 2200    CONTINUE
*
 2300 CONTINUE
*
*"         < 3. 対称・反対称成分 → 通常の並び >
*
      DO 3200 J = 1 , JMAX/2
         JN = J
         JS = JMAX+1 - J
         IF ( .NOT. LDPNM ) THEN
*"                                (  PNM を用いる場合 )
            JE = J
            JO = JMAX/2 + J
         ELSE
*"                                ( DPNM を用いる場合 )
            JE = JMAX/2 + J
            JO = JN
         ENDIF
*
         DO 3100 K = 1 , KMAX
            DO 3100 IM = 1 , IDIM
               ZDATA( IM,JN,K ) = ZDW( IM,JE,K ) + ZDW( IM,JO,K )
               ZDATA( IM,JS,K ) = ZDW( IM,JE,K ) - ZDW( IM,JO,K )
 3100    CONTINUE
*
 3200 CONTINUE
*
*
      RETURN
      END
**********************************************************************
      SUBROUTINE SPZ2W     !" ルジャンドル変換 (格子→スペクトル)
     M         ( WDATA ,
     I           ZDATA   ,
     C           PNM   , NMO   , GW    ,
     F           LDPNM , HFUNC ,
     D           JMAX  , KMAX  , IDIM  , JDIM  ,
     D           LMAX  , MMAX  , NMAX  , MINT  , NMDIM , JMXHF ,
     W           ZDW                                            )
*
*   [PARAM] 
      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       ZDATA ( IDIM, JDIM, KMAX )   !" 東西スペクトル
*
      REAL       PNM   ( NMDIM, JMXHF )       !" ルジャンドル関数
      INTEGER    NMO   ( 2, 0:MMAX , 0:LMAX ) !" スペクトルの添字順番
      REAL       GW    ( JDIM )               !" ガウス荷重
*
      LOGICAL    LDPNM                        !" ｙ微分フラグ
      CHARACTER  HFUNC*4                      !" 符号フラグ
*
*   [WORK] 
      REAL       ZDW   ( IDIM, JDIM, KMAX )   !" ワーク
*
*   [INTERNAL ONCE] 
      INTEGER    NMDIMD
      PARAMETER (NMDIMD=12000)
      INTEGER    MLIST ( NMDIMD )             !" リストベクトル
      INTEGER    JLIST ( NMDIMD )             !" リストベクトル
      LOGICAL    OLSET
      DATA       OLSET / .FALSE. /
      SAVE       MLIST, JLIST, OLSET
*
*   [INTERNAL WORK] 
      INTEGER    L, M, NM, K
      INTEGER    LEND, MM, IM
      INTEGER    J, JP, JE, JO, JN, JS
*
*
*"         < 0. リストベクトル作成 >
*
*   [ONCE] 
      IF ( .NOT. OLSET ) THEN
         OLSET = .TRUE.
*
         IF ( NMDIMD .LT. NMDIM ) THEN
            CALL MSGDMP( 'E','SPZ2W','INSUFFICIENT WORK STRAGE')
         ENDIF
*
         DO 200 M = 0, MMAX, MINT
            LEND = MIN ( LMAX, NMAX-M )
            MM   = M / MINT
            DO 100 L = 0, LEND
               MLIST (  NMO ( 1,M,L )  ) = 2* MM + 1
               MLIST (  NMO ( 2,M,L )  ) = 2* MM + 2
*
               IF ( MOD( L,2 ) .EQ. 0 ) THEN
                  JLIST (  NMO ( 1,M,L )  ) = 0
                  JLIST (  NMO ( 2,M,L )  ) = 0
               ELSE
                  JLIST (  NMO ( 1,M,L )  ) = JMAX / 2
                  JLIST (  NMO ( 2,M,L )  ) = JMAX / 2
               ENDIF
*
  100       CONTINUE
  200    CONTINUE
*
      ENDIF
*
*
      IF (      ( HFUNC(1:1) .NE. 'A' )
     &     .AND.( HFUNC(1:1) .NE. 'S' )    ) THEN
*
*"                 < 1. 和差をとらないときリセット >
*
         CALL RESET ( WDATA, NMDIM*KMAX )
      ENDIF
*
*"         < 2. ガウス荷重をかける >
*"              ( 通常の並び → 対称・反対称成分 )
*
      DO 2200 J = 1 , JMAX/2
*
         JN = J
         JS = JMAX+1 - J
         IF ( .NOT. LDPNM ) THEN
*"                                (  PNM を用いる場合 )
            JE = J
            JO = JMAX/2 + J
         ELSE
*"                                ( DPNM を用いる場合 )
            JE = JMAX/2 + J
            JO = JN
         ENDIF
*
         DO 2100 K = 1 , KMAX
            DO 2100 IM = 1 , IDIM
               ZDW( IM,JE,K ) = GW ( J ) *
     &                      ( ZDATA( IM,JN,K ) + ZDATA( IM,JS,K ) )
               ZDW( IM,JO,K ) = GW ( J ) *
     &                      ( ZDATA( IM,JN,K ) - ZDATA( IM,JS,K ) )
 2100    CONTINUE
*
 2200 CONTINUE
*
*
      IF (      ( HFUNC(1:1) .EQ. 'N' )
     &     .OR. ( HFUNC(1:1) .EQ. 'S' )    ) THEN
*
*"         < 3. ルジャンドル変換−負のとき   >
*
         DO 3300 J = 1 , JMAX/2
*
            DO 3200 K = 1 , KMAX
               DO 3100 NM = 1 , NMDIM
*
                  IM = MLIST( NM )
                  JP = JLIST( NM ) + J
*
                  WDATA( NM,K   )  = WDATA( NM,K   )
     &                             -  PNM ( NM,J ) * ZDW( IM,JP,K   )
 3100          CONTINUE
 3200       CONTINUE
*
 3300    CONTINUE
*
      ELSE
*
*"         < 4. ルジャンドル変換−正のとき   >
*
         DO 4300 J = 1 , JMAX/2
*
            DO 4200 K = 1 , KMAX
               DO 4100 NM = 1 , NMDIM
*
                  IM = MLIST( NM )
                  JP = JLIST( NM ) + J
*
                  WDATA( NM,K   )  = WDATA( NM,K   )
     &                             +  PNM ( NM,J ) * ZDW( IM,JP,K   )
 4100          CONTINUE
 4200       CONTINUE
*
 4300    CONTINUE
*
      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 ( IDIM, JDIM, KMAX ) !" データ
*
*   [WORK] 
      REAL       ZDW   ( IDIM, JDIM, KMAX ) !" ワーク
*
*   [INTERNAL WORK] 
      INTEGER    M, MR, MI, MM
      INTEGER    J, K
*
*
*"         < 1. 作業領域 >
*
      CALL COPY ( ZDW , ZDATA , IDIM*JDIM*KMAX )
*
*"         < 2. ｉｍ をかける >
*
      DO 2100 M = 0, MMAX, MINT
         MM = M / MINT
         MR = 2*MM + 1
         MI = 2*MM + 2
         DO 2110 K = 1, KMAX
            DO 2110 J = 1, JDIM
               ZDATA( MR, J, K ) = - REAL( M ) * ZDW  ( MI, J, K )
               ZDATA( MI, J, K ) =   REAL( M ) * ZDW  ( MR, J, K )
 2110    CONTINUE      
 2100 CONTINUE
*
      RETURN
      END