subroutine recfour(a,b,n,f,k)

        implicit none
        integer n,k
        real    a(n/2+1),b(n/2),f(n)

c ***** THIS SUBROUTINE RECOVERS THE FOURIER SERIES FROM THE
c ***** COEFFICIENTS COMPUTED BY n-POINT HARMONIC ANALYSIS

        integer i,j
        real    pi,ang,suma,sumb
        data    pi      / 3.1415926 /

c ***** COMPUTE SERIES FOR WAVENUMBER k=0
        if (k.eq.0) then
           do i = 1,n
              f(i) = a(1)/2.
           enddo

c ***** COMPUTE SERIES FOR WAVENUMBER k
        elseif ((k.gt.0).and.(k.lt.n/2)) then
           do i = 1,n
              ang  = 2.*pi*float(i-1)/float(n)
              f(i) = a(k+1)*cos(float(k)*ang)+b(k+1)*sin(float(k)*ang)
c              print*,'i= ',i
           enddo

c ***** COMPUTE FOR WAVENUMBER k=n/2
        elseif (k.eq.n/2) then
           do i = 1,n
              ang  = 2.*pi*float(i-1)/float(n)
              f(i) = 0.5*a(k+1)*cos(float(k)*ang)
           enddo

c ***** COMPUTE SERIES FOR k=-1 (FULL SERIES)
       elseif (k.eq.-1) then
           do i = 1,n
              ang  = 2.*pi*float(i-1)/float(n)
              suma = 0.
              sumb = 0.
              do j = 1,n/2-1
                 suma = suma + a(j+1)*cos(float(j)*ang)
              enddo
              do j = 1,n/2-1
                 sumb = sumb + b(j+1)*sin(float(j)*ang)
              enddo
              j    = n/2
              suma = suma + 0.5*a(1) + 0.5*a(j+1)*cos(float(j)*ang)
              f(i) = suma + sumb
           enddo
        elseif (k.lt.-1) then
           write(6,10) 'INVALID WAVENUMBER'
        else
           write(6,10) 'MAXIMUM OF WAVENUMBER',n/2,'ALLOWED.'
        endif
 10     format(i2)

        return
        end