* PACKAGE USPHE   !"  球面調和関数変換のテスト
*
*"  [HIS] 95/06/07(takepiro)
*
*
**********************************************************************
      PROGRAM TSPG2W   !" 格子→スペクトル
*
*"   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#endif
*
*   [PARAM] 
c$$$      INTEGER    IMAX
c$$$      INTEGER    JMAX
c$$$      INTEGER    KMAX
c$$$      INTEGER    IDIM
c$$$      INTEGER    JDIM
c$$$      INTEGER    LMAX
c$$$      INTEGER    MMAX
c$$$      INTEGER    NMAX
c$$$      INTEGER    MINT
c$$$      INTEGER    NMDIM
c$$$      INTEGER    JMXHF
*
*   [MODIFY] 
      REAL       WDATA ( NMDIM, KMAX     )    !" スペクトルデータ
*
*   [INPUT] 
      REAL       GDATA ( IDIM, JDIM, KMAX )    !" 格子点データ
*
      REAL       PNM   ( NMDIM, JMXHF )       !" ルジャンドル関数
      REAL       DPNM  ( NMDIM, JMXHF )       !" ルジャンドル関数微分
      INTEGER    NMO   ( 2, 0:MMAX , 0:LMAX ) !" スペクトルの添字順番
      REAL * 8   QSINLA( JDIM )              !" sin(緯度) (倍精度)
      REAL * 8   QGW    ( JDIM )               !" ガウス荷重
      REAL       GW    ( JDIM )               !" ガウス荷重
      REAL       TRIGS ( IDIM*2 )            !" 三角関数表
      INTEGER    IFAX  ( 10 )                !" IMAX の因数分解
*
      CHARACTER  HGRAD*4                      !" 微分フラグ
      CHARACTER  HFUNC*4                      !" 符号フラグ
*
*   [WORK] 
      REAL * 8   QPNM  ( 0:NMAX+1, 0:MMAX )  !" Ｐnm ルジャンドル
      REAL * 8   QDPNM ( 0:NMAX+1, 0:MMAX )  !" Ｐnm μ微分
      REAL       ZDATA ( IDIM*JDIM, KMAX )    !" 東西スペクトル
      REAL       WORK  ( IDIM*JDIM, KMAX )    !" ワーク
*
*   [INTERNAL WORK] 
      INTEGER    I, J, K
*
      CALL RESET( WDATA, NMDIM*KMAX )
*
*"         < 1.1 緯度, ガウス荷重 >
*
      PI = ATAN(1.) *4
      CALL GAUSS
     O         ( QSINLA , QGW    ,
     I           JMAX             )
*
*"         < 1.1 FFTの準備 >
*
      IF ( IMAX .GE. 2 ) THEN
         CALL RFFTIM
     I         ( IMAX  ,
     O           TRIGS , IFAX )
      ENDIF
*
*"         < 1.2 ルジャンドル関数の計算 >
*
      IF ( JMAX .GE. 2 ) THEN
         CALL SETPNM
     O         ( PNM   , DPNM  , NMO   ,
     I           QSINLA,
     D           JMAX  , MMAX  , LMAX  , NMAX  , MINT  ,
     D           NMDIM , JMXHF ,
     W           QPNM  , QDPNM                          )
      ELSE
         NMH = 0
         DO 2200 M = 0, MMAX, MINT
            NMH = NMH + 1
            NMO ( 1, M, 0 ) = 2* NMH - 1
            NMO ( 2, M, 0 ) = 2* NMH
 2200    CONTINUE
      ENDIF
*
*"         < 1.3 ガウス荷重 >
*
      DO 3100 J = 1, JMAX
         GW  ( J ) = QGW( J )
c$$$         GWDEL( J ) = QGW( J ) / ER / ( 1. - QSINLA(J)**2 )
 3100 CONTINUE
*
*"         < 格子データの設定 >
*
      WRITE ( 6, * ) 'N,M?'
      READ (5,*) NN,MM
      DO 4100 K=1, KMAX
         DO 4100 J=1, JDIM
            Y = QSINLA(J)
            DO 4100 I=1, IDIM
               X = 2 * PI * (I-1)/REAL(IDIM-1) 
               GDATA( I,J,K ) = YLMC(NN,MM,X,Y)
 4100 CONTINUE 
*
*"         < 球面調和函数変換 >
*
      CALL 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                                    )
*
*"         < 出力 >
*
      K = KMAX/2
      DO 5100 L=0,LMAX
         DO 5100 M=0,MMAX
            IF ( NMO(1,M,L) .LE. NMDIM ) THEN
               IF( ABS(WDATA( NMO(1,M,L),K ) ) .GT. 0.0001 )THEN
                  WRITE ( 6, * ) '1:', L+M, M, L, WDATA( NMO(1,M,L),K )
               ENDIF
            ENDIF
            IF ( NMO(1,M,L) .LE. NMDIM ) THEN
               IF( ABS(WDATA( NMO(2,M,L),K ) ) .GT. 0.0001 )THEN
                  WRITE ( 6, * ) '2:', L+M, M, L, WDATA( NMO(2,M,L),K )
               ENDIF
            ENDIF
 5100 CONTINUE 
*
      STOP
      END
*======================================================================
      FUNCTION YLMS(L,M,X,Y)
      REAL    YLMS, X, Y
      INTEGER L, M
      YLMS = PLGNDR(L,M,Y) * SIN(X)
      RETURN
      END
*
      FUNCTION YLMC(L,M,X,Y)
      REAL    YLMC, X, Y
      INTEGER L, M
      YLMC = PLGNDR(L,M,Y) * COS(X)
      RETURN
      END
*======================================================================
*"表題 ルジャンドル陪函数(plgndr)
*
*"履歴  93/11/07 竹広真一 Numerical Recipies より
*
*"注意  岩波数学公式と定義が異なる. (-1)^m のファクターがついている. 
*
*
      FUNCTION PLGNDR(L,M,X)
*
      REAL    PLGNDR        !" ルジャンドル函数 関数値
      INTEGER L, M          !" ルジャンドル函数 次数
      REAL    X             !" ルジャンドル函数 引数
*
      REAL    PMM, PMMP1, PLL, FACT, SOMX2
      INTEGER I, LL
*
      IF(M.LT.0.OR.M.GT.L.OR.ABS(X).GT.1.)PAUSE 'bad arguments'
      PMM=1.
      IF(M.GT.0) THEN
        SOMX2=SQRT((1.-X)*(1.+X))
        FACT=1.
        DO 11 I=1,M
          PMM=-PMM*FACT*SOMX2
          FACT=FACT+2.
11      CONTINUE
      ENDIF
      IF(L.EQ.M) THEN
        PLGNDR=PMM
      ELSE
        PMMP1=X*(2*M+1)*PMM
        IF(L.EQ.M+1) THEN
          PLGNDR=PMMP1
        ELSE
          DO 12 LL=M+2,L
            PLL=(X*(2*LL-1)*PMMP1-(LL+M-1)*PMM)/(LL-M)
            PMM=PMMP1
            PMMP1=PLL
12        CONTINUE
          PLGNDR=PLL
        ENDIF
      ENDIF
      RETURN
      END
