      SUBROUTINE GAUSS(NRACP,RACP,PG,SIA,RAD,PGSSIN2)
C
C     * REVISION:      NOV 3/98 - ML. ADD IMPLICIT REAL*I DECLARATION
C     *                               SO THAT ACCURACY WILL BE MAINTAINED
C     *                               EVEN WHEN RUNNING IN 32-BIT MODE.
C     ****************************************************************
C     AUTHOR: Evhen Yakimiw,RPN,Dorval
C             September 1995
C     ****************************************************************
C
C     COMMENTS: This subroutine uses a highly efficient and accurate 
C       method developed by the author for evaluating the Gaussian 
C       latitudes and the Gaussian weights required in the Gauss-Legendre 
C       quadrature rule. It is especially well suited for extremely high
C       order quadrature rules.
C       
C     ****************************************************************
C
C     *****************************************************************
C     NRACP        : NUMBER OF POSITIVE ROOTS FOR THE LEGENDRE POLYNOMIALS
C                    OF DEGREE 2*NRACP.
C     RACP(I)      : POSITIVE ROOTS FOR THE LEGENDRE POLYNOMIALS.
C                    THESE ARE XO, THE SINES OF THE GAUSSIAN LATITUDES.
C     PG(I)        : GAUSSIAN WEIGHTS AT THE GAUSSIAN LATITUDES.
C     SIA(I)       : SIN(COLATITUDE) = COS(LATITUDE).
C     RAD(I)       : COLATITUDE IN RADIANS.
C     PGSSIN2(I)   : PG(I)/(SIA(I)**2
C     SINM1(I)     : 1./SIA(I)      **** removed for CCRD use ****
C     SINM2(I)     : 1./SIA(I)**2   **** removed for CCRD use ****
C     SIN2(I)      : SIA(I)**2      **** removed for CCRD use ****
C     *****************************************************************
C
      IMPLICIT REAL*8 (A-H,O-Z),
     +INTEGER (I-N)
      REAL*8 RACP(NRACP),PG(NRACP),SIA(NRACP),RAD(NRACP)
      REAL*8 PGSSIN2(NRACP)
C 
      INTEGER LG
C======================================================================
C
C     Defining the constants needed
C
      N=2*NRACP
C---------------------------------
C     The following code may be used for even or odd values of N
C     ranging from 2 and up.
C
      M=(N+1)/2
      FN=FLOAT(N)
      F2N=FN+FN
      FN4=1./(384.*FN**4)
      FNN=FLOAT(N*(N+1))
      FNN2=FNN-2.
      FNN3=FLOAT(N*(N+1)-6)
      FNN4=FNN-12.
      FM=FLOAT(N-1)
      FP=FLOAT(N+1)
      LG=3.3*LOG(FN)
C
      AC1=SQRT(2.D0)
      DO 1 K=1,N
        FK=FLOAT(K)
        F2K=2.*FK
        F2KSQ=F2K*F2K
        AC1=AC1*SQRT(1.-1./F2KSQ)
  1   CONTINUE
      AC1=AC1/SQRT(FN+.5)
C
C---------------------------------
      X1=ACOS(-1.D0)/FLOAT(4*N+2)
      X2=1.-.125*FM/FN**3
      A=-1.
C**********************************
      DO 2 I=1,M
        A=A+4.
        T=A*X1
C---------------------------------
C     A high order accurate initial guess of the Gaussian latitudes is used.
C     Ref: F.G.Lether in J.Comput.App.Math.4,pp.47-52(1978). Precision O(n-5).
C---------------------------------
        SN2=1.-COS(T)*COS(T)
        XO=(X2-FN4*(39.-28./SN2))*COS(T)
C---------------------------------
C     To accurately determine the Gaussian latitudes, two methods are combined.
C     The corrections dx to the Gaussian sin(latitude)=xo are first evaluated 
C     by reversing the Taylor expansion of pn, the Legendre polynomials, near 
C     its zeros,followed by a Newton's rule. The procedure can be repeated.
C     In this code, a 5th order scheme is implemented, leading to Gaussian
C     latitudes accurate to machine precision for all values of n. 
C     NO NEED TO ITERATE.
C---------------------------------
        PM=1.
        PN=XO
        DO 3 K=2,N
          T1=XO*PN
          PP=T1-PM
          PP=PP-PP/FLOAT(K)+T1
          PM=PN
          PN=PP
   3    CONTINUE
        C2=(1.+XO)*(1.-XO)
        CM2=1./C2
        DP1=FN*(PM-XO*PN)
        DPN=DP1*CM2
        F1=PN/DPN
        F2=(2.*XO-FNN*F1)*CM2
        F3=(4.*XO*F2-FNN2)*CM2
        F4=(6.*XO*F3-FNN3*F2)*CM2
        F5=(8.*XO*F4-FNN4*F3)*CM2
        F22=F2*F2
C---------------------------------
C     Reversion of Taylor series (Handbook of Mathematical
C     Functions, Abramowitz and Stegun, p.16,Eq.3.6.25)
C
        DX=-F1*(1.+.5*F1*(F2+F1*(F22-F3/3.-F1*(F2*(10.*F3-15.*F22)-F4
     +     -F1*(3.*F2*F4+2.*F3*F3+21.*F22*(F22-F3)-.2*F5))/12.)))
C
C---------------------------------
C     Newton's rule (5th order)(Abramowitz and Stegun, p.18,Eq.3.9.5)
C
        P=F1+DX*(1.+.5*DX*(F2+DX*(F3+.25*DX*(F4+.2*DX*F5))/3.))
        DP=1.+DX*(F2+.5*DX*(F3+DX*(F4+.25*DX*F5)/3.))
C
C---------------------------------
        DX=DX-P/DP
        XO=XO+DX
C---------------------------------
C     This completes the evaluation of the Gaussian latitudes.
C---------------------------------
C
C     The next step calculates the Gaussian weights. To proceed, one needs
C     the values of the Legendre polynomials and its first derivative at
C     the Gaussian latitudes. Two algorithms are used. The first algorithm
C     gives the best accuracy at very high latitudes. The second algorithm 
C     is more accurate at lower latitudes as it is also much more effficient.
C
        THETA=ACOS(XO)
        C2=(1.+XO)*(1.-XO)
        C1=SIN(THETA)
        CM1=1./C1
        CM2=1./C2
        IF(I.LE.LG) THEN
C
          FKM=0.
          S1=0.
          S2=0.
          AC=1.
C
          DO 4 K=1,N-1,2
            FK=FLOAT(K)
            FKP=FK+1.
            AC0=FN-FKM
            ANG=AC0*THETA
            S1=S1+AC*COS(ANG)
            S2=S2-AC*AC0*SIN(ANG)
            AC=(FK*(F2N-FKM))/(FKP*(F2N-FK))*AC
            FKM=FKP
   4      CONTINUE
          IF((K-1).EQ.N) AC=0.5*AC
          AC0=FN-FKM
          ANG=AC0*THETA
C
          S1=S1+AC*COS(ANG)
          S2=S2-AC*AC0*SIN(ANG)
          PN=AC1*S1
          DP1=AC1*S2
C
          DPN=-CM1*DP1
        ELSE
          PM=1.
          PN=XO
          DO 5 K=2,N
            T1=XO*PN
            PP=T1-PM
            PP=PP-PP/FLOAT(K)+T1
            PM=PN
            PN=PP
   5      CONTINUE
          DP1=FN*(PM-XO*PN)
          DPN=DP1*CM2
C-----------------------------------
        ENDIF
        F1=PN/DPN
C---------------------------------
        GP=DPN*(C2-XO*F1+.5*(FNN+CM2)*F1**2)
        PGS=2./(GP*GP)
        PG(I)=C2*PGS
        RACP(I)=XO
C---------------------------------
C       Calculate other quantities 
C
        RAD(I)=THETA
        SIA(I)=C1
        PGSSIN2(I)=PGS
 2    CONTINUE
C ***********************************
      RETURN
      END
