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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
!-----------------------------------------------------------------------------
!
!
SUBROUTINE dt_comp
!
!
!-----------------------------------------------------------------------------
!
!
! Subroutine to compute the DT integration interval, FLY_2.0 Sept. 2001
!
!
!-----------------------------------------------------------------------------
USE fly_h
implicit none
INCLUDE 'mpif.h'
! Declaration of local variables.
! -------------------------------
REAL(KIND=8) :: dt_act,t_new,a_guess,z_new,max_dt_acc,mod_a
REAL(KIND=8) ::max_acc_local,mod_acc,mod_vel,max_mod_acc,max_mod_vel,max_modg_acc,max_modg_vel
REAL(KIND=8) :: ttn2,ap2,a_dot2,bp2
INTEGER :: p,ibod,ibod2,ibod3
REAL(KIND=8) :: dt_guess
!-----------------------------------------------------------------------------
! max acceleration component in all the system
!-----------------------------------------------------------------------------
max_acc_local=0.0
mod_vel=0.
max_mod_acc=0.
max_mod_vel=0.
max_modg_acc=0.
max_modg_vel=0.
IF(dt_var) THEN
ttn2=tnow+dtime2
ap2=f_ap/ttn2
a_dot2 = om_sum+Omega_l*(ttn2**three_o_alpha)
bp2 = 1./(alpha2*(hubble*hubble)*a_dot2*(ttn2**2.))
DO p=1,nb_res_loc(me+1)
mod_acc=abs(acc(1,p)*acc(1,p) + &
acc(2,p)*acc(2,p)+acc(3,p)*acc(3,p))
mod_vel=abs(vel(1,p)*vel(1,p) + vel(2,p)*vel(2,p) + &
vel(3,p)*vel(3,p))
mod_a=abs(-2.*ap2*(vel(1,p)*vel(1,p) + vel(2,p)*vel(2,p) + &
vel(3,p)*vel(3,p)) + bp2*(acc(1,p)*acc(1,p) + &
acc(2,p)*acc(2,p)+acc(3,p)*acc(3,p)))
IF(mod_a.ge.max_acc_local)THEN
max_acc_local=mod_a
ibod=p
ENDIF
IF(mod_acc.ge.max_mod_acc)THEN
max_mod_acc=mod_acc
ibod2=p
ENDIF
IF(mod_vel.ge.max_mod_vel) THEN
max_mod_vel=mod_vel
ibod3=p
ENDIF
ENDDO
max_acc_local=sqrt(max_acc_local)
NLONG=1
CALL MPI_ALLREDUCE(max_acc_local, max_dt_acc, NLONG, MPI_REAL8, &
MPI_MAX, MPI_COMM_WORLD, ierror)
CALL MPI_ALLREDUCE(max_mod_acc, max_modg_acc, NLONG, MPI_REAL8, &
MPI_MAX, MPI_COMM_WORLD, ierror)
CALL MPI_ALLREDUCE(max_mod_vel, max_modg_vel, NLONG, MPI_REAL8, &
MPI_MAX, MPI_COMM_WORLD, ierror)
dt_guess=sqrt(0.5*eps/max_dt_acc)
ELSE ! if dt_var
dt_guess=dtime
ENDIF ! if dt_var
!-----------------------------------------------------------------------------
! znow determination for output
!-----------------------------------------------------------------------------
wr_out32=.FALSE.
! We now use the new dt
t_new=tnow+dt_guess
a_guess=(t_new)**(1.0/alpha)
z_new=1.0/a_guess-1.0
dt_act=dt_guess
IF(z_new.le.z32_ou(next_out)) THEN
dt_act=(1/(1+z32_ou(next_out)))**alpha - tnow
IF(dt_act.eq.0.0) dt_act=dt_guess
!
wr_out32=.TRUE.
IF(next_out.eq.z32_last)THEN
halt_sim=1 !stop simulation at the end of this time-step
IF(me.eq.0) THEN
write(uterm,*)'SET STOP SIMULATION: at the end of the current time-step'
call flush(uterm)
ENDIF
ENDIF
ENDIF
dtime=dt_act
dtime2=dtime/2
RETURN
END