C
C     Test Harmon
C
      Parameter  (IDIM=50)
      Dimension  A(0:IDIM,0:IDIM), PNM(0:IDIM,0:IDIM)
C
      PI = 1
      PI = 4*Atan(PI)
      ONE = 1
C
      Read(5,1)  NMAX
    1   Format(I5)
C
C     Harmonic Coefficients
C
      RAN = PI/13
      Do  2  M=0, NMAX
      Do  2  N=0, NMAX
        A(N,M) = Sin((N+M+1)*RAN)
    2 Continue
C
      Read(5,3)  DTHETA, DPHI, R
    3   Format(3F10.0)
      THETA = (PI/180)*DTHETA
      PHI   = (PI/180)*DPHI
      X     = Cos(THETA)
      STH   = Sin(THETA)
C
C     Legendre Function Table
C
      Do  5  M=0, NMAX
        Call  PNM1(X, M, NMAX, PNM(M,M), LEX, IER)
        SCL = Exp(M*Log(STH) + LEX)
        Do  4  N=M, NMAX
          PNM(N,M) = PNM(N,M)*SCL
    4   Continue
    5 Continue
C
C     Straight Sum
C
      SN = 0
      SR = 0
      SP = 0
      ST = 0
      Do  7  N=NMAX, 1,-1
        TN = 0
        TP = 0
        TT = 0
        Do  6  M=N, 1,-1
          TN = TN + (A(N,M)*Cos(M*PHI)+A(M-1,N)*Sin(M*PHI))*PNM(N,M)
          TP = TP + M*(-A(N,M)*Sin(M*PHI)+A(M-1,N)*Cos(M*PHI))*PNM(N,M)
          DP = M*X*PNM(N,M)/STH
          If( M.lt.N )  DP = DP - PNM(N,M+1)*Sqrt((N+M+ONE)*(N-M))
          TT = TT + (A(N,M)*Cos(M*PHI)+A(M-1,N)*Sin(M*PHI))*DP
    6   Continue
        SN = R*SN + TN + A(N,0)*PNM(N,0)
        SR = R*SR + (N+1)*(TN+A(N,0)*PNM(N,0))
        SP = R*SP + TP
        ST = R*ST + (TT-A(N,0)*PNM(N,1)*Sqrt(N*(N+ONE)/2))
    7 Continue
      SN = R*SN + A(0,0)
      SR = R*SR + A(0,0)
      SP = R*SP/STH
      ST = R*ST
C
C     Harmon
C
      Call  Harmon(R, THETA, PHI, NMAX, A, IDIM, VN, LN, IER)
      FC = LN
      VN = VN*Exp(FC)
      Call  Harmor(R, THETA, PHI, NMAX, A, IDIM, VR, LR, IER)
      FC = LR
      VR = VR*Exp(FC)
      Call  Harmot(R, THETA, PHI, NMAX, A, IDIM, VT, LT, IER)
      FC = LT
      VT = VT*Exp(FC)
      Call  Harmop(R, THETA, PHI, NMAX, A, IDIM, VP, LP, IER)
      FC = LP
      VP = VP*Exp(FC)
      EN = SN - VN
      ER = SR - VR
      ET = ST - VT
      EP = SP - VP
C
      Write(6,10)  NMAX, DTHETA, DPHI, R
   10   Format(/3X,'Test Harmon'//5X,'Nmax  =',I4/
     *         5X,'Theta =',F11.6/5X,'Phi   =',F11.6/
     *         5X,'R     =',F11.6)
      Write(6,11)  SN, VN, EN, SR, VR, ER, ST, VT, ET, SP, VP, EP
   11   Format(/12X,'Naive',10X,'Harmo*',9X,'Diff'/
     *         5X,'Vn',1P3E15.6/5X,'Vr',1P3E15.6/5X,'Vt',1P3E15.6/
     *         5X,'Vp',1P3E15.6)
C
      Stop
      End
C
C     Spherical Harmonics Expansion
C
C       M. Saito,  May 18, 1993
C
      Subroutine  Harmon(R, THETA, PHI, NMAX, A, IDIM, VN, LEX, IER)
C
C     Input
C       R    : r
C       THETA: Colatitude
C       PHI  : Longitude
C       NMAX : Max degree
C       A(N,M): Harmonic coefficients a(n,m) and b(n,m)
C       IDIM : 1-st dimension of A
C     Output
C       VN   : VN(r,theta,phi)
C       LEX  : Exponent of scale factor
C       IER  : Error code
C              = 0 ; No error
C
C     Subroutine : SPNM1
C
C     This program calculates the harmonic function
C                         nmax      n
C       VN(r,theta,phi) = Sum r**n Sum [a(n,m)cos(m*phi)
C                         n=0      m=0
C         + b(n,m)sin(m*phi)]Pnm(cos(theta))
C
C     where Pnm is the normalized associated Legendre function.  Pnm
C     is normalized to
C       +1
C       Int [Pnm(x)]**2 dx = 2  (m=0)  or  = 4 (m>0)
C       -1
C
C     a(n,m) and b(n,m) are in
C       a(n,m) in A(N,M)
C       b(n,m) in A(M-1,N)
C
      Dimension  A(0:IDIM,0:*)
      Data  KEX/20/
C
      If( NMAX.le.0 .or. IDIM.lt.NMAX )  Go to  90
C
      EPS = KEX
      EPS = Exp(-EPS)
      K   = 0
C
      CT = Cos(THETA)
      ST = Sin(THETA)
      CP = Cos(PHI)
      SP = Sin(PHI)
      X2 = CP*2
      RT = R*ST
      C0 = 0
      C1 = 0
      S0 = 0
      S1 = 0
      U1 = 1
C
      Do  1  M=NMAX, 0, -1
        Call  SPNM1(R, CT, M, NMAX, A, IDIM, CQ, SQ, L, IER)
        SCL = L - K
        SCL = Exp(-SCL)
        If( L.gt.K )  then
          C0 = C0*SCL
          C1 = C1*SCL
          S0 = S0*SCL
          S1 = S1*SCL
          K  = L
        else If( L.lt.K )  then
          CQ = CQ/SCL
          SQ = SQ/SCL
        Endif
        If( M.eq.0 )  Go to  1
        C2 = C1
        C1 = C0
        S2 = S1
        S1 = S0
        U2 = U1
        U1 = (M+1)*2
        U1 = RT*Sqrt(1 + 1/U1)
        C0 = CQ + U1*(X2*C1 - U2*C2)
        S0 = SQ + U1*(X2*S1 - U2*S2)
        If( Abs(C0)*EPS.gt.1 .or. Abs(S0)*EPS.gt.1 )  then
          C0 = C0*EPS
          C1 = C1*EPS
          S0 = S0*EPS
          S1 = S1*EPS
          K  = K + KEX
        else If( Abs(C0).lt.EPS .or. Abs(S0).lt.EPS )  then
          C0 = C0/EPS
          C1 = C1/EPS
          S0 = S0/EPS
          S1 = S1/EPS
          K  = K - KEX
        Endif
    1 Continue
C
      THR = 3
      FIV = 5
      VN  = CQ + Sqrt(THR)*RT*(CP*C0 + SP*S0 - Sqrt(FIV/4)*RT*C1)
C
      LEX = K
      IER = 0
      Return
C
   90 Write(6,91)
   91   Format(/3X,5('?'),3X,'Invalid Input in HARMON'/)
      IER =-1
      Return
      End