-
Notifications
You must be signed in to change notification settings - Fork 1
/
doric.f90
181 lines (141 loc) · 6.14 KB
/
doric.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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
!>
!! \brief This module contains routines for calculating the ionization evolution of a single grid point.
!!
!! Module for Capreole / C2-Ray (f90)
!!
!! \b Author: Garrelt Mellema
!!
!! \b Date:
!!
!! \b Version: H only
!!
!! History: the doric routine goes back to the early 1990-ies. It is an
!! analytical solution to the problem of time dependent photo-ionization
!! under the assumption of constant electron density. In this case the
!! solution is an exponential which goes towards an equilibrium solution.
!! This idea was taken from Schmidt-Voigt & Koeppen (1987). After the
!! C2-Ray algorithm was conceived, the same time dependent solution is
!! used to find the time-averaged ionization fraction needed by the
!! algorithm.
!! Note that the original version of doric also contained helium, following
!! the algorithm from Schmidt-Voigt & Koeppen (1987). This algorithm is
!! wrong, the proper one is given by Altay et al. (2008) and Friedrich
!! et al. (2010).
module doric_module
! This module contains routines having to do with the calculation of
! the ionization evolution of a single grid point.
! - doric : time dependent solution of the ionization equation
! - coldens : column density of a single cell.
! - coldensh_bndry : Boundary condition for H column density
use precision, only: dp
implicit none
! No public data associated with this module
contains
!=======================================================================
!> calculates time dependent ionization state for hydrogen
subroutine doric (dt,temp0,rhe,rhh,xfh,xfh_av,phi)
! calculates time dependent ionization state for hydrogen
! Author: Garrelt Mellema
! Date: 11-May-2005 (f90) (15 June 2004
! Version:
! Simplified version of the original doric as used in Coral,
! - H only
! Modifications:
! 15-Jun-2004 (GM): test for (1.0d0-ee)/deltht to avoid fp errors.
! 11-May-2005 (GM): f90
use cgsconstants, only: bh00,albpow,colh0,temph0
use material, only: clumping
use radiation, only: photrates
use tped, only: electrondens ! should this really be used inside doric?
use c2ray_parameters, only: epsilon
real(kind=dp),intent(in) :: dt !< time step
real(kind=dp),intent(in) :: temp0 !< local temperature
real(kind=dp),intent(inout) :: rhe !< electron density
real(kind=dp),intent(in) :: rhh !< H density (or total density?)
real(kind=dp),dimension(0:1),intent(out) :: xfh !< H ionization fractions
real(kind=dp),dimension(0:1),intent(inout) :: xfh_av !< H ionization fractions (time-averaged)
type(photrates),intent(in) :: phi !< photo-ionization rates
real(kind=dp) :: brech0,sqrtt0,acolh0
real(kind=dp) :: rhe0,xfh1old,xfh0old,aih0,delth,eqxfh0,eqxfh1
real(kind=dp) :: aphoth0
real(kind=dp) :: deltht,ee
real(kind=dp) :: avg_factor
! find the hydrogen recombination rate at the local temperature
! We assume that the recombination rate can be approximated as
! a single powerlaw.
brech0=clumping*bh00*(temp0/1e4)**albpow
! find the hydrogen collisional ionization rate at the local
! temperature
sqrtt0=sqrt(temp0)
acolh0=colh0*sqrtt0*exp(-temph0/temp0)
! Find the true photo-ionization rate
aphoth0=phi%h
! determine the hydrogen ionization state and
! electron density
! (schmidt-voigt & koeppen 1987)
! The electron density is the time-averaged value.
! This needs to be iterated, as the electron density
! depends on the hydrogen ionization state.
! In this version the iteration
! is assumed to take place outside the doric routine.
! Save old values
xfh1old=xfh(1)
xfh0old=xfh(0)
aih0=aphoth0+rhe*acolh0
delth=aih0+rhe*brech0
eqxfh1=aih0/delth
eqxfh0=rhe*brech0/delth
deltht=delth*dt
ee=exp(-deltht)
xfh(1)=(xfh1old-eqxfh1)*ee+eqxfh1
xfh(0)=(xfh0old-eqxfh0)*ee+eqxfh0
rhe=electrondens(rhh,xfh) ! should this really be used inside doric?
! determine neutral densities (take care of precision fluctuations)
!if (xfh(0) < epsilon .and. abs(xfh(0)).lt.1.0e-10) then
if (xfh(0) < epsilon) then
xfh(0)=epsilon
xfh(1)=1.0_dp-epsilon
endif
! Determine average ionization fraction over the time step
! Mind fp fluctuations. (1.0-ee)/deltht should go to 1.0 for
! small deltht, but finite precision leads to values slightly
! above 1.0 and for very small values even to 0.0.
if (deltht.lt.1.0e-8) then
avg_factor=1.0
else
avg_factor=(1.0-ee)/deltht
endif
! The question here is whether it would be better to calculate
! xfh_av(0) first, and xfh_av(1) from it.
xfh_av(1)=eqxfh1+(xfh1old-eqxfh1)*avg_factor
xfh_av(0)=1.0_dp-xfh_av(1)
! Take care of precision
!if (xfh_av(0).lt.epsilon.and.abs(xfh_av(0)).lt.1.0e-10) xfh_av(0)=epsilon
if (xfh_av(0) < epsilon) xfh_av(0)=epsilon
end subroutine doric
! =======================================================================
!> Calculates the column density (of hydrogen)
!! for a cell of ionization fraction xh, length path,
!! and density ndenstime dependent ionization state for hydrogen
function coldens(path,xh0,ndens)
! Calculates the column density (of hydrogen)
! for a cell of ionization fraction xh, length dr,
! and density ndens
real(kind=dp) :: coldens
real(kind=dp),intent(in) :: path !< path length through cell
real(kind=dp),intent(in) :: xh0 !< neutral H fraction
real(kind=dp),intent(in) :: ndens !< number density of H
! Column density over a distance dr (cell)
coldens=xh0*ndens*path
end function coldens
!=======================================================================
!> Sets the boundary condition for the hydrogen column density
function coldensh_bndry()
! Sets the boundary condition for the hydrogen column density
use cgsphotoconstants, only: sigh
use radiation, only: tauHI
real(kind=dp):: coldensh_bndry
! Column density at the left of the first cell
coldensh_bndry=tauHI/sigh
end function coldensh_bndry
end module doric_module