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
| program theta0
implicit none
real*8, dimension(:,:), allocatable :: R,T1,delta,phi,C,A,B,delta1
real*8 :: c0,c1,c2,alpha,beta,r0,dt,dx,Tmax,l,g
integer :: nb_iter,m,i,j
print*, ' m = '
read*, m
!100
print*, ' nb_iter = '
read*, nb_iter
!qu'importe
allocate(R(0:nb_iter,0:m))
allocate(delta(0:nb_iter-1,1:m-1))
allocate(T1(1:nb_iter,1:m-1))
allocate(phi(0:nb_iter,0:m-1))
allocate(C(0:nb_iter,0:m-1))
allocate(B(0:nb_iter,0:m-1))
allocate(delta1(0:nb_iter-1,1:m-1))
! Constantes
c0 = 10.
alpha = 0.01
beta = 0.1
r0 = 1.
Tmax = 0.5
! Conditions initiales
C(0,1:m-1) = c0
phi(0,1:m-1) = alpha*c0 + beta
R(0,1:m-1) = 0.
B(0,1:m-1) = 0.
! Conditions aux limites en 0
R(0:nb_iter,0) = r0
do i = 0,nb_iter
C(i,0) = c0*exp(-r0*i/10)
phi(i,0) = alpha*C(i,0) + beta
B(i,0) = r0/phi(i,0)
end do
! Conditions aux limites en m
R(0:nb_iter,m) = R(0:nb_iter,m-1)
dt = Tmax/nb_iter
dx = 1./m
!B(0:nb_iter,0:m) = R(0:nb_iter,0:m)/phi(0:nb_iter,0:m)
do i = 1,nb_iter
do j = 2,m
C(i,j-1) = C(i-1,j-1)*exp(-dt*R(i-1,j-1))
phi(i,j-1) = alpha*C(i,j-1) + beta
delta(i-1,j-1) = (phi(i-1,j-1)+phi(i-1,j))*(B(i-1,j)-B(i-1,j-1))-(phi(i-1,j-2)+phi(i-1,j-1))*(B(i-1,j-1)-B(i-1,j-2))
delta1(i-1,j-1) = (dt/dx**2)*delta(i-1,j-1)
T1(i,j-1) = 1./(1-dt*C(i,j-1))
R(i,j-1) = T1(i,j-1)*(R(i-1,j-1)+delta1(i-1,j-1))
B(i,j-1) = R(i,j-1)/phi(i,j-1)
end do
end do
open(unit=10,file='theta1.dat')
do i = 0,m-1
write(10,*) i*dx, R(0,i), R(1,i), R(5,i), R(10,i), R(20,i), R(50,i), R(300,i), C(0,i), C(1,i), C(3,i), C(5,i), C(10,i), C(20,i)
end do
end program theta0
end program theta0 |
