!***********************************************************************
!
!
                          SUBROUTINE init_pars
!
!
!***********************************************************************
!
!
!     Subroutine to initialize system parameters that depend on
!     either the input data or defined PARAMETERS.
!
!
!=======================================================================
!
	 USE fly_h 
 	implicit none
!   Declaration of local variables.
!   -------------------------------
        
	REAL(KIND=8) :: xw,xw2,deldrg,xw3,xw4,zero,two,two_o_three,f2
        INTEGER :: i,minustwo


!=======================================================================
!   Initialize misc. useful numbers.
!   --------------------------------
        minustwo = -2.
        tiny=1.e-20
        zero=0.
        one=1.
        two=2.
	alpha=0.5
	two_o_three=2./3.
	one_o_alpha=1./alpha
	three_o_alpha=3./alpha
	f2=2.+one_o_alpha
	alpha2=alpha*alpha

!   Initialize position and velocity times, 1/2 timestep.
!   -----------------------------------------------------
	tnow=1./((1.+znow)**alpha)

        tpos=tnow
        dtime2=.5*dtime
        tol2inv=1./(tol*tol)

!-----------------------------------------------------------------------
!   Initialize variables and arrays for gravitational field smoothing
!   interpolation.  Interpolation performed in distance.
!-----------------------------------------------------------------------
        deldrg=2./ninterp
        
	DO i=0,1+ninterp  
           xw=i*deldrg
           xw2=xw*xw
           xw3=xw2*xw
           xw4=xw2*xw2
     	   IF(xw.LE.one) THEN  
     	    phsmooth(i)=-2.*xw3*(one/3.-3.*xw2/20.+xw3/20.)+7.*xw/5.
            acsmooth(i)=xw3*(4./3.-6.*xw2/5.+0.5*xw3)
 	   ELSE
     	    phsmooth(i)=-one/15.+8.*xw/5.-xw3*(4./3.-xw+3.*xw2/10.-xw3/30.)
            acsmooth(i)=-one/15.+8.*xw3/3.-3.*xw4+6.*xw3*xw2/5.-xw4*xw2/6.
     	   ENDIF      	   
           IF(xw.GE.two) THEN
             phsmooth(i)=one
	     acsmooth(i)=one
 	   ENDIF

        ENDDO
	
	

!--------------------------------------------		
! Ewald definition
!
! Fix L2=Lbox/2.
!--------------------------------------------

	Lbox= rsize
	L2=Lbox/2.

! Fix linv (used in EWALD_R.F)
	
	linv=nknots/Lbox

	om_sum=Omega_cdm+Omega_hdm
	f_ap=(2.*(1.+alpha+Omega_l)-(Omega_cdm+Omega_hdm))/(4.0*alpha)

	
        RETURN
        END