-
Notifications
You must be signed in to change notification settings - Fork 5
/
photonstatistics.F90
295 lines (235 loc) · 9.63 KB
/
photonstatistics.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
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
!>
!! \brief This module contains data and routines for calculating the photon statistics
!!
!! Module for C2-Ray
!!
!! \b Author: Garrelt Mellema
!!
!! \b Date: 26-Feb-2008
!!
!! \b Version: 3D
module photonstatistics
! Photon statistics
! photon_loss: this is kept here, but calculated in the evolve module.
use precision, only: dp
use my_mpi, only: rank
use file_admin, only: logf
use cgsconstants, only: albpow,bh00,colh0,temph0
use cgsphotoconstants, only: sigh => sigma_HI_at_ion_freq
use sizes, only: mesh
use grid, only: vol
!use material, only: ndens, temper, clumping, clumping_point
use density_module, only: ndens
use temperature_module, only: temper
use temperature_module, only: temperature_states_dbl
use temperature_module, only: get_temperature_point
use clumping_module, only: clumping, clumping_point
use tped, only: electrondens
use sourceprops, only: NormFlux_stellar, NumSrc
use radiation_sed_parameters, only: S_star
use radiation_sizes, only: NumFreqBnd
use c2ray_parameters, only: type_of_clumping
implicit none
!> true if checking photonstatistics
logical,parameter :: do_photonstatistics=.true.
!> Total number of recombinations
real(kind=dp) :: totrec
!> Total number of collisional ionizations
real(kind=dp) :: totcollisions
!> Change in number of neutral H atoms
real(kind=dp) :: dh0
!> Total number of ionizing photons used
real(kind=dp) :: total_ion
!> Total number of photons lost in LLSs
real(kind=dp) :: LLS_loss
!> Grand total number of ionizing photons used
real(kind=dp) :: grtotal_ion
!> Grand total number of ionizing photons produced
real(kind=dp) :: grtotal_src
!> Number of photons leaving the grid
real(kind=dp) :: photon_loss(NumFreqBnd)
real(kind=dp),private :: h0_before !< number of H atoms at start of time step
real(kind=dp),private :: h0_after !< number of H atoms at end of time step
real(kind=dp),private :: h1_before !< number of H ions at start of time step
real(kind=dp),private :: h1_after !< number of H ions at end of time step
integer,private :: i !< mesh loop index (x)
integer,private :: j !< mesh loop index (y)
integer,private :: k !< mesh loop index (z)
contains
!----------------------------------------------------------------------------
!> Initialize the photon statistics
subroutine initialize_photonstatistics ()
! set total number of ionizing photons used to zero
grtotal_ion=0.0
grtotal_src=0.0
end subroutine initialize_photonstatistics
!----------------------------------------------------------------------------
!> Call the individual routines needed for photon statistics calculation
subroutine calculate_photon_statistics (dt,xh_l,xh_r)
real(kind=dp),intent(in) :: dt !< time step
#ifdef ALLFRAC
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3),0:1),intent(in) :: xh_l
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3),0:1),intent(in) :: xh_r
#else
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3)),intent(in) :: xh_l
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3)),intent(in) :: xh_r
#endif
! Call the individual routines needed for this calculation
call state_after (xh_l) ! number of neutrals after integration
call total_rates (dt,xh_r) ! total photons used in balancing recombinations etc.
call total_ionizations () ! final statistics
end subroutine calculate_photon_statistics
!----------------------------------------------------------------------------
!> Calculates the number of neutrals and ions at the start of the time step
subroutine state_before (xh_l)
#ifdef ALLFRAC
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3),0:1),intent(in) :: xh_l
#else
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3)),intent(in) :: xh_l
#endif
! Photon statistics: calculate the number of neutrals before integration
h0_before=0.0
h1_before=0.0
do k=1,mesh(3)
do j=1,mesh(2)
do i=1,mesh(1)
#ifdef ALLFRAC
h0_before=h0_before+ndens(i,j,k)*xh_l(i,j,k,0)
h1_before=h1_before+ndens(i,j,k)*xh_l(i,j,k,1)
#else
h0_before=h0_before+ndens(i,j,k)*(1.0_dp-xh_l(i,j,k))
h1_before=h1_before+ndens(i,j,k)*xh_l(i,j,k)
#endif
enddo
enddo
enddo
h0_before=h0_before*vol
h1_before=h1_before*vol
end subroutine state_before
!----------------------------------------------------------------------------
!> Calculates total number of recombinations and collisions
subroutine total_rates(dt,xh_l)
real(kind=dp),intent(in) :: dt !< time step
#ifdef ALLFRAC
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3),0:1),intent(in) :: xh_l
#else
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3)),intent(in) :: xh_l
#endif
real(kind=dp),dimension(0:1) :: yh
real(kind=dp) :: ndens_p ! needed because ndens may be single precision
type(temperature_states_dbl) :: temperature
! Photon statistics: Determine total number of recombinations/collisions
! Should match the code in doric_module
totrec=0.0
totcollisions=0.0
do k=1,mesh(3)
do j=1,mesh(2)
do i=1,mesh(1)
#ifdef ALLFRAC
yh(0)=xh_l(i,j,k,0)
yh(1)=xh_l(i,j,k,1)
#else
yh(0)=1.0_dp-xh_l(i,j,k)
yh(1)=xh_l(i,j,k)
#endif
ndens_p=ndens(i,j,k)
call get_temperature_point (i,j,k,temperature)
! Set clumping to local value if we have a clumping grid
if (type_of_clumping == 3 .or. type_of_clumping == 4 .or. type_of_clumping == 5) &
call clumping_point (i,j,k)
totrec=totrec+ndens_p*yh(1)* &
electrondens(ndens_p,yh)* &
clumping*bh00*(temperature%average/1e4)**albpow
totcollisions=totcollisions+ndens_p* &
yh(0)*electrondens(ndens_p,yh)* &
colh0*sqrt(temperature%average)*&
exp(-temph0/temperature%average)
enddo
enddo
enddo
totrec=totrec*vol*dt
totcollisions=totcollisions*vol*dt
end subroutine total_rates
!----------------------------------------------------------------------------
!> Calculates the number of neutrals and ions at the end of the time step
subroutine state_after(xh_l)
#ifdef ALLFRAC
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3),0:1),intent(in) :: xh_l
#else
real(kind=dp),dimension(mesh(1),mesh(2),mesh(3)),intent(in) :: xh_l
#endif
! Photon statistics: Calculate the number of neutrals after the integration
h0_after=0.0
h1_after=0.0
do k=1,mesh(3)
do j=1,mesh(2)
do i=1,mesh(1)
#ifdef ALLFRAC
h0_after=h0_after+ndens(i,j,k)*xh_l(i,j,k,0)
h1_after=h1_after+ndens(i,j,k)*xh_l(i,j,k,1)
#else
h0_after=h0_after+ndens(i,j,k)*(1.0_dp-xh_l(i,j,k))
h1_after=h1_after+ndens(i,j,k)*xh_l(i,j,k)
#endif
enddo
enddo
enddo
h0_after=h0_after*vol
h1_after=h1_after*vol
end subroutine state_after
!----------------------------------------------------------------------------
!> Calculate the total number of ionizing photons used
subroutine total_ionizations ()
! Photon statistics: Total number of new ionizations
dh0=(h0_before-h0_after)
total_ion=totrec+dh0
end subroutine total_ionizations
!----------------------------------------------------------------------------
!> Calculate the total number of photons lost in LLSs
subroutine total_LLS_loss (phi_out,coldensh_LLS)
real(kind=dp),intent(in) :: phi_out !< Photon flux out off cell
real(kind=dp),intent(in) :: coldensh_LLS !< local (cell) column density of LLSs
! Note LLS_loss is set to zero at the start of a time step in
! evolve:pass_all_sources
real(kind=dp) :: tau_LLS
tau_LLS=sigh*coldensh_LLS
! Note: This expression is only valid for grey opacities!!
! GM (110131): Why do we divide by exp(-tau_LLS)? The numbers
! generated by this expression seem to be wrong. Comment out
! exp(-tau_LLS)
LLS_loss = LLS_loss + phi_out*(1.0-exp(-tau_LLS))!/exp(-tau_LLS)
end subroutine total_LLS_loss
!----------------------------------------------------------------------------
!> Calculate the total number of ionizing photons used
subroutine report_photonstatistics (dt)
real(kind=dp),intent(in) :: dt
real(kind=dp) :: totalsrc,photcons,total_photon_loss,total_LLS_loss
! By this time photon_loss has been renormalized to photon_loss per cell
! in the simulation. Multiply by the number of cells to restore the
! total photon loss.
total_photon_loss=sum(photon_loss)*dt* &
real(mesh(1))*real(mesh(2))*real(mesh(3))
total_LLS_loss = LLS_loss*dt
totalsrc=sum(NormFlux_stellar(1:NumSrc))*s_star*dt
!photcons=(total_ion-totcollisions)/totalsrc
photcons=(total_ion+LLS_loss-totcollisions)/totalsrc
if (rank == 0) then
write(logf,"(8(1pe10.3))") &
total_ion, totalsrc, &
photcons, &
dh0/total_ion, &
totrec/total_ion, &
total_LLS_loss/totalsrc, &
total_photon_loss/totalsrc, &
totcollisions/total_ion
write(logf,*) h1_before,h1_after
endif
end subroutine report_photonstatistics
!----------------------------------------------------------------------------
!> Calculate the total number of ionizing photons produced
subroutine update_grandtotal_photonstatistics (dt)
real(kind=dp),intent(in) :: dt !< time step
grtotal_src=grtotal_src+sum(NormFlux_stellar(1:NumSrc))*S_star*dt
grtotal_ion=grtotal_ion+total_ion-totcollisions
end subroutine update_grandtotal_photonstatistics
end module photonstatistics