c enable d-lines compliation if you wish detailed printout of process
      implicit real*8(a-h,o-z)
      parameter (n=33 )
      dimension   z(n), zout(n)
      dimension   bs(n)  ! work vector

! two g-exp filter parameters: beta and L
      data beta/0.5d0/
      data L/4/           ! L must be greater than zero

      dsd=0.729432d1    ! seed; will be changed by repeated calls to randno()

      pie2=4.0d0*dasin(1.0d0)

      do i=1,n
      r=randno(dsd)
         z(i)=1.0d0*dsin(pie2*dble(i-1)/dble(n-1)) + (r - 0.5d0)*0.1d0   ! arbitrary signal
      enddo

      call       gexp(z,n,beta,L, bs, zout)

      write(6,*) '  i    z(i)        zout(i)'
      do i=1,n
        write(6,1000) i,z(i),zout(i)
 1000 format(1x,i3,1p,2g13.5)
      enddo

      stop
      end

      subroutine gexp(z,n,beta,L, bs, zout)
      implicit real*8(a-h,o-z)
c     z(n) - input
c     zout(n) - output
c     n - number of points to be smoothed
c     L == lrw - length of window (can be zero but not recommended)
c     beta == alfa - is the exponent  (very large beta effectively causes no smoothing, small beta flattens the points)
c     bs - work vector
      dimension z(n), zout(n)
      dimension bs(n)             ! work vector (declared outside)
      dimension g(1000)           ! work vector  - size of the g filter (declared in sub())

      Lsave=L
      L=iabs(L)
      if(L.eq.0) stop 'L=0; program terminated.'

      do L1=1,L

d     write(6,*) ' '
d     write(6,*) 'L1 = ',L1

      ng1=2*L1+1
d     write(6,*) 'Construct the filter as in Eqs. (22) '
      sum=0.0
      do i=1,ng1
        x=dble(i)-dble(ng1+1)/dble(2)
        if(L1.gt.0) x=x/dble(L1)
        g(i)=exp(-beta*x**2)
        sum=sum+g(i)
d       write(6,*) i,'  ',x, '  ', g(i)
      enddo

d     write(6,*) ' '
d     write(6,*) 'Normalize the filter as in Eq. (23)'
      do i=1,ng1
        g(i)=g(i)/sum
d       write(6,*) i,'  ',g(i)
      enddo

d     write(6,*) ' '
d     write(6,*) 'Execute the convolution as in Eq. (24) or (25)'
      do is=1,n-2*L1
        bs(is)=0.0
          do ir=1,2*L1+1
            bs(is)=bs(is) + g(ir)*z(is+ir-1)
          enddo
        zout(is+L1)=bs(is)
d       write(6,*) is,'  '       , z(is+L1), '   ',zout(is+L1)
      enddo

      enddo  ! loop L1

      zout(1)=z(1)
      zout(n)=z(n)

d     write(6,*) ' '
d     write(6,*) 'z(i) and zout(i)                           '
      do i=1,n
d       write(6,*) i,'  ',z(i), '  ', zout(i)
!       z(i)=zout(i)       ! overwrite the input if you wish
      enddo
      L=Lsave

      return
      end

      real*8 function randno (dsd)
      implicit real*8(a-h,o-z)
      real*4 ratio
      data    d2p31m/2147483647.d0/
      data    d2p31 /2147483648.d0/
      dsd = dmod(16807.d0*dsd,d2p31m)
      ratio = dsd / d2p31
      randno = dble(ratio)
      return
      end

