Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
!***********************************************************************
!
!
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