C
C     ..................................................................
C
C        SUBROUTINE DRHARM
C
C        PURPOSE
C           FINDS THE FOURIER COEFFICIENTS OF ONE DIMENSIONAL DOUBLE
C           PRECISION REAL DATA
C
C        USAGE
C           CALL DRHARM(A,M,INV,S,IFERR)
C
C        DESCRIPTION OF PARAMETERS
C           A     - A DOUBLE PRECISION VECTOR
C                   AS INPUT, CONTAINS ONE DIMENSIONAL REAL DATA. A IS
C                   2*N+4 CORE LOCATIONS, WHERE N = 2**M. 2*N REAL
C                   NUMBERS ARE PUT INTO THE FIRST 2*N CORE LOCATIONS
C                   OF A
C                   AS OUTPUT, A CONTAINS THE FOURIER COEFFICIENTS
C                   A0/2,B0=0,A1,B1,A2,B2,...,AN/2,BN=0 RESPECTIVELY IN
C                   THE FIRST 2N+2 CORE LOCATIONS OF A
C           M     - AN INTEGER WHICH DETERMINES THE SIZE OF THE VECTOR
C                   A. THE SIZE OF A IS 2*(2**M) + 4
C           INV   - A VECTOR WORK AREA FOR BIT AND INDEX MANIPULATION OF
C                   DIMENSION ONE EIGHTH THE NUMBER OF REAL INPUT, VIZ.,
C                   (1/8)*2*(2**M)
C           S     - A DOUBLE PRECISION VECTOR WORK AREA FOR SINE TABLES
C                   WITH DIMENSION THE SAME AS INV
C           IFERR - A RETURNED VALUE OF 1 MEANS THAT M IS LESS THAN 3 OR
C                   GREATER THAN 20. OTHERWISE IFERR IS SET = 0
C
C        REMARKS
C           THIS SUBROUTINE GIVES THE FOURIER COEFFICIENTS OF 2*(2**M)
C           REAL POINTS. SEE SUBROUTINE DHARM FOR THREE DIMENSIONAL,
C           DOUBLE PRECISION, COMPLEX FOURIER TRANSFORMS.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           DHARM
C
C        METHOD
C           THE FOURIER COEFFICIENTS A0,B0=0,A1,B1,...,AN,BN=0 ARE
C           OBTAINED FOR INPUT XJ, J=0,1,2,...,2N-1 FOR THE FOLLOWING
C           EQUATION (PI = 3.14159...)
C
C                 N-1                                               J
C     XJ=(1/2)A0+SUM (AK*COS(PI*J*K/N)+BK*SIN(PI*J*K/N))+(1/2)AN(-1)
C                 K=1
C
C           SEE REFERENCE UNDER SUBROUTINE DHARM
C
C     ..................................................................
C
      SUBROUTINE DRHARM(A,M,INV,S,IFERR)
      DIMENSION A(1),L(3),INV(1),S(1)
      DOUBLE PRECISION A,SI,AP1IM,FN,CO,CIRE,AP2IM,S,SS,DEL,CIIM,AP1RE,
     1 CNIRE,SC,SIS,AP2RE,CNIIM
      IFSET=1
      L(1)=M
      L(2)=0
      L(3)=0
      NTOT=2**M
      NTOT2 = 2*NTOT
      FN = NTOT
      DO   3 I = 2,NTOT2,2
   3  A(I) = -A(I)
      DO   6 I = 1,NTOT2
   6  A(I) = A(I)/FN
      CALL DHARM(A,L,INV,S,IFSET,IFERR)
C
C     MOVE LAST HALF OF A(J)S DOWN ONE SLOT AND ADD A(N) AT BOTTOM TO
C     GIVE ARRAY FOR A1PRIME AND A2PRIME CALCULATION
C
   21 DO  52 I=1,NTOT,2
      J0=NTOT2+2-I
      A(J0)=A(J0-2)
   52 A(J0+1)=A(J0-1)
      A(NTOT2+3)=A(1)
      A(NTOT2+4)=A(2)
C
C     CALCULATE A1PRIMES AND STORE IN FIRST N SLOTS
C     CALCULATE A2PRIMES AND STORE IN SECOND N SLOTS IN REVERSE ORDER
      K0=NTOT+1
      DO 104 I=1,K0,2
      K1=NTOT2-I+4
      AP1RE=.5*(A(I)+A(K1))
      AP2RE=-.5*(A(I+1)+A(K1+1))
      AP1IM=.5*(-A(I+1)+A(K1+1))
      AP2IM=-.5*(A(I)-A(K1))
      A(I)=AP1RE
      A(I+1)=AP1IM
      A(K1)=AP2RE
  104 A(K1+1)=AP2IM
      NTO = NTOT/2
  110 NT=NTO+1
      DEL=3.141592653589793/DFLOAT(NTOT)
      SS=DSIN(DEL)
      SC=DCOS(DEL)
      SI=0.0
      CO=1.0
C
C     COMPUTE C(J)S FOR J=0 THRU J=N
  114 DO 116 I=1,NT
      K6=NTOT2-2*I+5
      AP2RE=A(K6)*CO+A(K6+1)*SI
      AP2IM=-A(K6)*SI+A(K6+1)*CO
      CIRE=.5*(A(2*I-1)+AP2RE)
      CIIM=.5*(A(2*I)+AP2IM)
      CNIRE=.5*(A(2*I-1)-AP2RE)
      CNIIM=.5*(A(2*I)-AP2IM)
      A(2*I-1)=CIRE
      A(2*I)=CIIM
      A(K6)=CNIRE
      A(K6+1)=-CNIIM
      SIS=SI
      SI=SI*SC+CO*SS
  116 CO=CO*SC-SIS*SS
C
C     SHIFT C(J)S FOR J=N/2+1 TO J=N UP ONE SLOT
      DO 117 I=1,NTOT,2
      K8=NTOT+4+I
      A(K8-2)=A(K8)
  117 A(K8-1)=A(K8+1)
      DO 500 I=3,NTOT2,2
      A(I) = 2. * A(I)
  500 A(I + 1) = -2. * A(I + 1)
      RETURN
      END