-
Notifications
You must be signed in to change notification settings - Fork 0
/
Mad_He2.f90
129 lines (101 loc) · 2.85 KB
/
Mad_He2.f90
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
117
118
119
120
121
122
123
124
125
126
127
128
module m_mad_he2
contains
subroutine Mad_He2 (time,temperature,ntime,Apatite_age,iflag)
use m_tridag
! iflag=1 He in Apatite
! iflag=2 He in Zircon
! iflag=3 Ar in K-feldspar
! iflag=4 Ar in Biotite
! iflag=5 Ar in Muscovite
! iflag=6 Ar in Hornblende
implicit none
integer(4) :: ntime, iflag, i, istep, itime, n, nstep
real(8) :: time(ntime),temperature(ntime), Apatite_age
real(8), dimension(:), allocatable :: age,diag,sup,inf,f
real(8) :: agei, alpha, beta, D0a2, Da2now, Da2then
real(8) :: dr, dt, dt0, Ea, f1, f2, fact, fstep, R
real(8) :: temp, tempf, temps
! D0a2=10.**7.7*3600.*24.*365.25e6
! Ea=36.2e3*4.184
agei = 0.0
if (iflag.eq.1) then
! He in Ap
! 2014.06.10, WK: changes from Byron Adams
D0a2=3.6e6*3600.*24.*365.25e6
Ea=142.0e3
elseif (iflag.eq.2) then
! He in Zi
D0a2=4.6e3*3600.*24.*365.25e6
Ea=169.0e3
elseif (iflag.eq.3) then
! Ar in Ks
D0a2=5.6*3600.*24.*365.25e6
Ea=120.e3
elseif (iflag.eq.4) then
! Ar in Bi 500mic
D0a2=160.*3600.*24.*365.25e6
Ea=211.e3
elseif (iflag.eq.5) then
! Ar in Mu 500mic
D0a2=13.2*3600.*24.*365.25e6
Ea=183.e3
elseif (iflag.eq.6) then
! Ar in Ho 500mic
! 2014.06.10, WK: changes from Byron Adams
D0a2=9.6*3600.*24.*365.25e6
Ea=268.0e3
else
stop 'option not supported in Mad_He'
endif
R=8.314
dt0=0.1
n=100
allocate (age(n),diag(n),sup(n),inf(n),f(n))
age=0.
do itime=1,ntime-1
dt=dt0
nstep=max(1,int((time(itime)-time(itime+1)+tiny(dt))/dt))
dt=(time(itime)-time(itime+1))/nstep
alpha=0.5
dr=1./(n-1)
! beta determines the geometry: 2 = spherical
! 1 = cylindrical
beta=2.
temps=temperature(itime)
tempf=temperature(itime+1)
Da2now=D0a2*exp(-Ea/R/(temps+273.))
do istep=1,nstep
fstep=float(istep)/(nstep)
temp=temps+(tempf-temps)*fstep
Da2then=Da2now
Da2now=D0a2*exp(-Ea/R/(temp+273.))
f1=alpha*dt*Da2now/dr**2
f2=(1.-alpha)*dt*Da2then/dr**2
do i=2,n-1
diag(i)=1.+2.*f1
sup(i)=-f1*(1.+beta/(i-1)/2.)
inf(i)=-f1*(1.-beta/(i-1)/2.)
f(i)=age(i)+f2* &
((age(i+1)-2.*age(i)+age(i-1)) + &
beta*(age(i+1)-age(i-1))/(i-1)/2.)+dt
enddo
diag(1)=1.
sup(1)=-1.
f(1)=0.
diag(n)=1.
inf(n)=0.
f(n)=0.
call tridag (inf,diag,sup,f,age,n)
agei=0.
do i=1,n
fact=1.
if (i.eq.1.or.i.eq.n) fact=0.5
agei=agei+age(i)*fact*dr**3*(i-1)*(i-1)
enddo
agei=3.*agei
enddo
enddo
Apatite_age=agei
return
end subroutine Mad_He2
end module m_mad_he2