* PACKAGE USPHEZ  !" 球面調和関数変換 (南北1次元) GCM5(USPHEZ)
*
*"  [HIS] !" 90/08/31(numaguti)
*
**********************************************************************
      SUBROUTINE SPW2GM    !" 球面調和関数変換 (スペクトル→格子)
     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    J, K
*
*
      IF ( IMAX .NE. 1 ) THEN
         CALL MSGDMP( 'E','SPW2G','THIS ROUTINE IS FOR AXSYMMETRIC')
      ENDIF
*
      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. 'X' ) THEN
*
*"         < 1. ｘ微分 = 0 >
*
          CALL RESET ( WORK , IDIM*JDIM*KMAX )
          LOFFS = .FALSE.
*
      ELSE
*
*"         < 2. LDPNM : ｙ微分のフラグ >
*
         IF ( HGRAD(1:1) .EQ. 'Y' ) THEN
            LDPNM = .TRUE.
            LOFFS = .FALSE.
         ELSE
            LDPNM = .FALSE.
         ENDIF
*
*"         < 3. スペクトル → [ zonal mean ] value >
*
*
         IF ( LOFFS ) THEN
            DO 3000 K = 1, KMAX
               DOFFS( K ) = WDATA( NMO(1,0,0), K )
               WDATA( NMO(1,0,0), K ) = 0.  
 3000       CONTINUE
         ENDIF
*
         CALL SPW2ZM
     O           ( WORK  ,
     I             WDATA ,
     C             PNM   , NMO   ,
     F             LDPNM ,
     D             JMAX  , KMAX  , JDIM  ,
     D             LMAX  , NMAX  , NMDIM , JMXHF ,
     W             ZDATA                           )
*
      ENDIF
*
      IF ( LOFFS ) THEN
         DO 3100 K = 1, KMAX
            WDATA( NMO(1,0,0), K ) = DOFFS( K )
 3100    CONTINUE
*
         DO 3200 K = 1, KMAX
#ifdef SYS_SX3
*vdir noloopchg
#endif
            DO 3200 J = 1, JDIM
               WORK ( J,K ) = WORK ( J,K ) + DOFFS( K )
 3200    CONTINUE
      ENDIF
*
*"         < 5. 出力データ >
*
      IF      ( HFUNC(1:1) .EQ. 'A' ) THEN
*
*"                               ( add )
        DO 5000 K = 1, KMAX
           DO 5000 J = 1, JDIM
                  GDATA ( J,K ) = GDATA( J,K ) + WORK( J,K )
 5000   CONTINUE
*
      ELSE IF ( HFUNC(1:1) .EQ. 'S' ) THEN
*
*"                               ( sub )
        DO 5100 K = 1, KMAX
          DO 5100 J = 1, JDIM
                  GDATA ( J,K ) = GDATA( J,K ) - WORK( J,K )
 5100   CONTINUE
*
      ELSE IF ( HFUNC(1:1) .EQ. 'N' ) THEN
*
*"                               ( negative )
        DO 5200 K = 1, KMAX
          DO 5200 J = 1, JDIM
                  GDATA ( J,K ) = - WORK( J,K )
 5200   CONTINUE
*
      ELSE
*"                               ( positive )
        DO 5300 K = 1, KMAX
          DO 5300 J = 1, JDIM
                  GDATA ( J,K ) = WORK( J,K )
 5300   CONTINUE
*
      ENDIF
*
      RETURN
      END
**********************************************************************
      SUBROUTINE SPG2WM    !" 球面調和関数変換 (格子→スペクトル)
     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    J, K
*
      IF ( IMAX .NE. 1 ) THEN
         CALL MSGDMP( 'E','SPG2W','THIS ROUTINE IS FOR AXSYMMETRIC')
      ENDIF
*
      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. 'X' ) THEN
*
*"         < 1. ｘ微分 = 0 >
*
         IF (      ( HFUNC(1:1) .NE. 'A' )
     &        .AND.( HFUNC(1:1) .NE. 'S' )    ) THEN
            CALL RESET ( WDATA , NMDIM*KMAX )
*
         ENDIF
*
      ELSE
*
*"         < 2. LDPNM : ｙ微分のフラグ >
*
         IF ( HGRAD(1:1) .EQ. 'Y' ) THEN
            LDPNM = .TRUE.
            LOFFS = .FALSE.
         ELSE
            LDPNM = .FALSE.
         ENDIF
*
*
*"         < 3. [ zonal mean ] wave → スペクトル >
*
         IF ( LOFFS ) THEN
            DO 3000 K = 1, KMAX
               DOFFS ( K ) = GDATA( 1,K )
#ifdef SYS_SX3
*vdir noloopchg
#endif
               DO 3020 J = 1, JDIM
                  GDATA( J,K ) = GDATA( J,K ) - DOFFS( K )
 3020          CONTINUE
 3000       CONTINUE
         ENDIF
*
         CALL SPZ2WM
     M           ( WDATA ,
     I             GDATA ,
     C             PNM   , NMO   , GW    ,
     F             LDPNM , HFUNC ,
     D             JMAX  , KMAX  , JDIM  ,
     D             LMAX  , NMAX  , NMDIM , JMXHF ,
     W             WORK                            )
*
         IF ( LOFFS ) THEN
            IF (      ( HFUNC(1:1) .EQ. 'N' )
     &           .OR. ( HFUNC(1:1) .EQ. 'S' )    ) THEN
               DO 3100 K = 1, KMAX
                  WDATA( NMO(1,0,0),K ) = WDATA( NMO(1,0,0),K )
     &                                  - DOFFS( K )
 3100          CONTINUE
            ELSE
               DO 3200 K = 1, KMAX
                  WDATA( NMO(1,0,0),K ) = WDATA( NMO(1,0,0),K )
     &                                  + DOFFS( K )
 3200          CONTINUE
            ENDIF
*
            DO 3300 K = 1, KMAX
               DO 3300 J = 1, JDIM
                  GDATA( J,K ) = GDATA( J,K ) + DOFFS( K )
 3300       CONTINUE
         ENDIF
*
      ENDIF
*
      RETURN
      END
**********************************************************************
      SUBROUTINE SPW2ZM    !" ルジャンドル変換 (スペクトル→格子)
     O         ( ZDATA   ,
     I           WDATA ,
     C           PNM   , NMO   ,
     F           LDPNM ,
     D           JMAX  , KMAX  , JDIM  ,
     D           LMAX  , NMAX  , NMDIM , JMXHF ,
     W           ZDW                            )
*
*   [PARAM] 
      INTEGER    JMAX
      INTEGER    KMAX
      INTEGER    JDIM
      INTEGER    LMAX
      INTEGER    NMAX
      INTEGER    NMDIM
      INTEGER    JMXHF
*
*   [OUTPUT] 
      REAL       ZDATA ( JDIM, KMAX    )      !" 東西スペクトル
*
*   [INPUT] 
      REAL       WDATA ( NMDIM, KMAX   )      !" スペクトル
*
      REAL       PNM   ( NMDIM, JMXHF )       !" ルジャンドル関数
      INTEGER    NMO   ( 2, 0:LMAX )          !" スペクトルの添字順番
*
      LOGICAL    LDPNM                        !" ｙ微分フラグ
*
*   [WORK] 
      REAL       ZDW   ( JDIM, KMAX )         !" ワーク
*
*   [INTERNAL ONCE] 
      INTEGER    NMDIMD
      PARAMETER (NMDIMD=500)
      INTEGER    JLIST ( NMDIMD )             !" リストベクトル
      LOGICAL    OLSET
      DATA       OLSET / .FALSE. /
      SAVE       JLIST, OLSET
*
*   [INTERNAL WORK] 
      INTEGER    NM, K, L, LEND
      INTEGER    J, JP, JN, JS, JE, JO
*
*"         < 0. リストベクトル作成 >
*
*   [ONCE] 
      IF ( .NOT. OLSET ) THEN
        OLSET = .TRUE.
*
        IF ( NMDIMD .LT. NMDIM ) THEN
            CALL MSGDMP( 'E','SPW2ZM','INSUFFICIENT WORK STRAGE')
        ENDIF
*
        LEND = MIN ( LMAX, NMAX )
*
        DO 100 L = 0, LEND
*
          IF ( MOD( L,2 ) .EQ. 0 ) THEN
            JLIST (  NMO ( 1,L )  ) = 0
          ELSE
            JLIST (  NMO ( 1,L )  ) = JMAX / 2
          ENDIF
*
  100   CONTINUE
*
      ENDIF
*
*"         < 1. 作業領域リセット >
*
      CALL RESET ( ZDW  , JDIM*KMAX )
      CALL RESET ( ZDATA, JDIM*KMAX )
*
*"         < 2.  ルジャンドル変換 >
*
      DO 2300 NM = 1 , NMDIM
         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( JP,K ) = ZDW( 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 = J
        ENDIF
*
        DO 3100 K = 1 , KMAX
          ZDATA( JN,K ) = ZDW( JE,K ) + ZDW( JO,K )
          ZDATA( JS,K ) = ZDW( JE,K ) - ZDW( JO,K )
 3100   CONTINUE
 3200 CONTINUE
*
*
      RETURN
      END
**********************************************************************
      SUBROUTINE SPZ2WM    !" ルジャンドル変換 (格子→スペクトル)
     M         ( WDATA ,
     I           ZDATA   ,
     C           PNM   , NMO   , GW    ,
     F           LDPNM , HFUNC ,
     D           JMAX  , KMAX  , JDIM  ,
     D           LMAX  , NMAX  , NMDIM , JMXHF ,
     W           ZDW                             )
*
*   [PARAM] 
      INTEGER    JMAX
      INTEGER    KMAX
      INTEGER    JDIM
      INTEGER    LMAX
      INTEGER    NMAX
      INTEGER    NMDIM
      INTEGER    JMXHF
*
*   [MODIFY] 
      REAL       WDATA ( NMDIM, KMAX   )      !" スペクトル
*
*   [INPUT] 
      REAL       ZDATA ( JDIM, KMAX    )      !" 東西スペクトル
*
      REAL       PNM   ( NMDIM, JMXHF )       !" ルジャンドル関数
      INTEGER    NMO   ( 2, 0:LMAX )          !" スペクトルの添字順番
      REAL       GW    ( JDIM )               !" ガウス荷重
*
      LOGICAL    LDPNM                        !" ｙ微分フラグ
      CHARACTER  HFUNC*4                      !" 符号フラグ
*
*   [WORK] 
      REAL       ZDW   ( JDIM, KMAX )         !" ワーク
*
*   [INTERNAL ONCE] 
      INTEGER    NMDIMD
      PARAMETER (NMDIMD=4000)
      INTEGER    JLIST ( NMDIMD )             !" リストベクトル
      LOGICAL    OLSET
      DATA       OLSET / .FALSE. /
      SAVE       JLIST, OLSET
*
*   [INTERNAL WORK] 
      INTEGER    NM, K, L, LEND
      INTEGER    J, JP, JN, JS, JE, JO
*
*"         < 0. リストベクトル作成 >
*
*   [ONCE] 
      IF ( .NOT. OLSET ) THEN
        OLSET = .TRUE.
*
        IF ( NMDIMD .LT. NMDIM ) THEN
           CALL MSGDMP( 'E','SPZ2WM','INSUFFICIENT WORK STRAGE')
        ENDIF
*
        LEND = MIN ( LMAX, NMAX )
*
        DO 100 L = 0, LEND
*
          IF ( MOD( L,2 ) .EQ. 0 ) THEN
            JLIST (  NMO ( 1,L )  ) = 0
          ELSE
            JLIST (  NMO ( 1,L )  ) = JMAX / 2
          ENDIF
*
  100   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 = J
        ENDIF
*
        DO 2100 K = 1 , KMAX
          ZDW( JE,K ) = GW ( J ) *
     &                      ( ZDATA( JN,K ) + ZDATA( JS,K ) )
          ZDW( JO,K ) = GW ( J ) *
     &                      ( ZDATA( JN,K ) - ZDATA( 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
*
              JP = JLIST( NM ) + J
*
              WDATA( NM,K )  = WDATA( NM,K   )
     &                       -  PNM ( NM,J   )
     &                           * ZDW( 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
*
              JP = JLIST( NM ) + J
*
              WDATA( NM,K )  = WDATA( NM,K   )
     &                       +  PNM ( NM,J   )
     &                           * ZDW( JP,K )
 4100       CONTINUE
 4200     CONTINUE
 4300   CONTINUE
*
      ENDIF
*
      RETURN
      END
