forked from lanl/SuperNu
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtbxsmod.f
251 lines (250 loc) · 8.45 KB
/
tbxsmod.f
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
* © 2023. Triad National Security, LLC. All rights reserved.
* This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National
* Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of
* Energy/National Nuclear Security Administration. All rights in the program are reserved by Triad
* National Security, LLC, and the U.S. Department of Energy/National Nuclear Security Administration.
* The Government is granted for itself and others acting on its behalf a nonexclusive, paid-up,
* irrevocable worldwide license in this material to reproduce, prepare. derivative works, distribute
* copies to the public, perform publicly and display publicly, and to permit others to do so.
module tbxsmod
c ---------------
implicit none
************************************************************************
* Opacity tables by C. Fontes
************************************************************************
c-- number of tabulated elements, densities, temperatures
integer,parameter :: tb_nrho=17
integer,parameter :: tb_ntemp=27
integer :: tb_nelem=0
c-- density, temperature table points
real*8 :: tb_rho(tb_nrho) !(nrho)
real*8 :: tb_temp(tb_ntemp) !(ntemp)
c-- raw table input
real*4,allocatable,private :: tb_raw(:,:,:,:,:) !(ncol,ngr,ntemp,nrho,nelem)
c-- grey scattering opacity
real*8,allocatable :: tb_sig(:,:,:)
c-- group-coarsened table
real*4,allocatable :: tb_cap(:,:,:,:) !(ng,ntemp,nrho,nelem)
c-- number of energy points
integer,parameter,private :: ngr=14900
c-- indicies of structure elements (from input.str)
integer, allocatable :: tb_ielem(:) !(nelem)
c
save
public
c
contains
c
c
c
subroutine tbxs_dealloc
c -----------------------
implicit none
if(allocated(tb_ielem)) deallocate(tb_ielem)
if(allocated(tb_sig)) deallocate(tb_sig)
if(allocated(tb_cap)) deallocate(tb_cap)
end subroutine tbxs_dealloc
c
c
c
subroutine read_tbxs
c --------------------
use miscmod, only:lcase
use elemdatamod, only: elem_neldata, elem_data
use inputstrmod, only: str_nabund,str_iabund
implicit none
************************************************************************
* Read Chris Fontes opacity data tables to total opacity array.
************************************************************************
integer,parameter :: ncol=7
character(3),parameter :: word1='op_'
character(4),parameter :: word2='_1Em'
character(9),parameter :: word3='gcc.table'
integer :: istat,ierr
integer :: ielem,itemp,irho,iirho,l
real*8 :: dmy
character(6) :: sdmy
character(2) :: fid,fnum
character(20) :: fname
c
c-- sanity check
if(str_nabund==0) stop 'read_tbxs: str_nabund=0'
c
c-- find number of tabulated elements from iabund
do ielem=1,elem_neldata
if(any(str_iabund==ielem)) tb_nelem=tb_nelem+1
enddo
c-- need at least one tabulated element in input.str
if(tb_nelem<=0) stop 'read_tbxs: tb_nelem<=0'
c-- store correct indices from elem_data
allocate(tb_ielem(tb_nelem))
l = 0
do ielem=1,elem_neldata
if(any(str_iabund==ielem)) then
l=l+1
tb_ielem(l)=ielem
endif
enddo
c
c-- raw table
allocate(tb_raw(ncol,ngr,tb_ntemp,tb_nrho,tb_nelem))
c
c-- elements
do l=1,tb_nelem
c-- skip if not in input structure
ielem=tb_ielem(l)
c-- file element
fid=lcase(trim(elem_data(ielem)%sym))
do irho=tb_nrho,1,-1
c-- density value
iirho=irho+3
tb_rho(tb_nrho-irho+1)=10d0**(-iirho)
c-- file density id
if(iirho/10==0) then
write(fnum,'("0"i1)') iirho
else
write(fnum,'(i2)') iirho
endif
c-- file name
if(fid(2:2) == ' ') then
fname=trim(word1//fid(1:1)//word2//fnum//word3)
else
fname=trim(word1//fid//word2//fnum//word3)
endif
open(4,file='Table/'//adjustl(fname),status='old',
& action='read',iostat=istat)
c-- require all possible data (for now)
if(istat/=0) then
write(6,*) 'file: ',fname
stop 'read_tbxs: missing file'
endif
do itemp=1,tb_ntemp
c-- temperature value
read(4,*) sdmy, dmy, tb_temp(itemp)
read(4,*)
c-- all data at temp-rho(-elem) point
read(4,*,iostat=ierr) tb_raw(:,:,itemp,tb_nrho-irho+1,l)
if(ierr/=0) stop 'read_tbxs format err: body'
enddo
c-- ensure no residual file data
read(4,*,iostat=ierr) sdmy
if(ierr/=-1) then
write(6,*) 'sdmy: ',sdmy
write(6,*) 'file: ',fname
stop 'read_tbxs: body too long'
endif
close(4)
enddo
enddo
c
end subroutine read_tbxs
c
c
c
subroutine coarsen_tbxs(lopac,ng,wl)
c ------------------------------------------------------------------
use miscmod, only:binsrch
use physconstmod
implicit none
logical, intent(in) :: lopac(4)
integer, intent(in) :: ng
real*8, intent(in) :: wl(ng+1)
************************************************************************
* Coursen Chris Fontes opacity tables.
************************************************************************
integer :: ig,igr,itemp,irho,ielem,n
integer :: ig1,ig2
real*4 :: cap,cap1,cap2,capl,capr
real*8 :: help1,help2,help3,help4,wll,wlr
real*8,allocatable :: evinvarr(:)
c
c-- coarse energy array
allocate(evinvarr(ng+1))
c-- final opacity tables
allocate(tb_sig(tb_ntemp,tb_nrho,tb_nelem))
allocate(tb_cap(ng,tb_ntemp,tb_nrho,tb_nelem))
c-- initialize
evinvarr=pc_ev*wl/(pc_h*pc_c)
tb_sig=0d0
tb_cap=0.
c-- quick exit
if(all(lopac)) return
c
do ielem=1,tb_nelem
do irho=1,tb_nrho
do itemp=1,tb_ntemp
c-- scattering opacity
if(.not.lopac(4)) tb_sig(itemp,irho,ielem)=
& dble(sum(tb_raw(7,:,itemp,irho,ielem))/ngr)
c-- absorption opacity
do igr=ngr-1,1,-1
c-- binary search group containing igr point
wll=1d0/dble(tb_raw(1,igr+1,itemp,irho,ielem))
wlr=1d0/dble(tb_raw(1,igr,itemp,irho,ielem))
ig1=binsrch(wll,evinvarr,ng+1,.true.)
ig1=max(ig1,1)
ig2=binsrch(wlr,evinvarr,ng+1,.true.)
ig2=min(ig2,ng)
c-- exclude points out of wl domain
if(ig2==0 .or. ig1==ng+1) cycle
do ig=ig1,ig2
c-- interpolate trapezoid
help1=max(wll,evinvarr(ig))
help2=min(wlr,evinvarr(ig+1))
help3=wlr*(help1-wll)/(help1*(wlr-wll)) ! high energy, low wavelength basis
help4=wll*(wlr-help2)/(help2*(wlr-wll)) ! low energy, high wavelength basis
c-- add opacity contributions
cap=0.
c-- bound-bound
if(.not.lopac(1)) then
capl=tb_raw(4,igr+1,itemp,irho,ielem)
capr=tb_raw(4,igr,itemp,irho,ielem)
cap1=capl*sngl(1d0-help3)+capr*sngl(help3)
cap2=capl*sngl(help4)+capr*sngl(1d0-help4)
cap=cap+.5*sngl((help2-help1)/(help1*help2))*(cap1+cap2)
endif
c-- bound-free
if(.not.lopac(2)) then
capl=tb_raw(5,igr+1,itemp,irho,ielem)
capr=tb_raw(5,igr,itemp,irho,ielem)
cap1=capl*sngl(1d0-help3)+capr*sngl(help3)
cap2=capl*sngl(help4)+capr*sngl(1d0-help4)
cap=cap+.5*sngl((help2-help1)/(help1*help2))*(cap1+cap2)
endif
c-- free-free
if(.not.lopac(3)) then
capl=tb_raw(6,igr+1,itemp,irho,ielem)
capr=tb_raw(6,igr,itemp,irho,ielem)
cap1=capl*sngl(1d0-help3)+capr*sngl(help3)
cap2=capl*sngl(help4)+capr*sngl(1d0-help4)
cap=cap+.5*sngl((help2-help1)/(help1*help2))*(cap1+cap2)
endif
c-- tally absorption opacity in group
tb_cap(ig,itemp,irho,ielem) =
& tb_cap(ig,itemp,irho,ielem)+cap
enddo
enddo
c-- average opacity
tb_cap(:,itemp,irho,ielem) =
& tb_cap(:,itemp,irho,ielem) *
& sngl(evinvarr(2:)*evinvarr(:ng) /
& (evinvarr(2:)-evinvarr(:ng)))
enddo
enddo
enddo
c
c-- deallocate large raw table
deallocate(evinvarr)
deallocate(tb_raw)
c
c-- convert temperature from eV to K
tb_temp=tb_temp*pc_ev/pc_kb
c
c-- print alloc size (keep this updated)
c---------------------------------------
n=tb_nelem*tb_nrho*tb_ntemp*(8+ng*4)/1024 !kB
write(6,*) 'ALLOC tbxs :',n,"kB",n/1024,"MB",n/1024**2,"GB"
c
end subroutine coarsen_tbxs
c
end module tbxsmod