17 integer iupdte,nvar,md,mtable,njon
18 real(kind=4):: eosxmin,ur,ul,ut,eps,tff,grv,abnd,dbox
19 real(kind=4),
allocatable,
dimension(:):: tmean,tamp,rhm,xcorr,thmin,thmax, &
20 dth,eemin,eemax,deetab,tab
21 integer,
allocatable,
dimension(:):: itab,mtab
22 integer,
save:: ncall=0, mbox
23 character(len=64):: top, tablefile, file
24 logical,
save:: do_eos
25 logical,
save:: first_time=.true.
29 integer,
parameter:: n_lambda=4
31 integer :: n_lambda = n_lambda
32 real :: w_lambda(n_lambda) = 1.0
35 logical :: on = .false.
39 procedure:: lookup_table
41 procedure:: temperature
43 integer,
save:: verbose=0
44 type(
eos_t),
public:: eos
50 SUBROUTINE init (self)
52 real,
dimension(:,:,:),
pointer:: d, ee, pg
53 real,
save :: kappa = 1.0
55 namelist /eos_params/ do_eos, top, tablefile, kappa, verbose
57 call trace%begin (
'eos_t%init')
63 tablefile =
'data/eos/stagger/table.dat' 64 rewind(io%input);
read (io%input,eos_params)
65 if (io%master)
write (*,eos_params)
69 file = trim(top)//
'/'//trim(tablefile)
70 open (io_unit%data, file=trim(file), form=
'unformatted', status=
'old', action=
'read')
71 read (io_unit%data) md,iupdte,nvar,mbox,eosxmin,dbox,ul,ut,ur,eps,tff,grv,abnd
73 allocate (tmean(md),tamp(md),rhm(md),xcorr(md),thmin(md),thmax(md))
74 allocate (dth(md),eemin(md),eemax(md),deetab(md),itab(md),mtab(md))
75 if (io%do_trace.and.io%master) print *,mpi%rank,
'reading 12 times', md,
' values' 76 read (io_unit%data) tmean,tamp,rhm,xcorr,thmin,thmax,dth,eemin,eemax,deetab,itab,mtab
77 read (io_unit%data) mtable
78 allocate (tab(mtable))
79 read (io_unit%data) tab
82 print *, trim(tablefile)//
': md, mtable =', md, mtable
83 allocate (d(1,1,1), ee(1,1,1), pg(1,1,1))
87 print *,
'd, ee, pg =', d, ee, pg
88 deallocate (d, ee, pg)
90 print *,
' density min(E) max(E) min(T) max(T)' 92 print
'(i4,1p,e12.4,2(2x,2e12.4))', ir, rhm(ir), eemin(ir), eemax(ir), &
93 tab(itab(ir)+2*3*mtab(ir)), tab(itab(ir)+(2*3+1)*mtab(ir)-1)
107 FUNCTION pressure (self, d, e)
RESULT (pg)
109 real(kind=KindScalarVar),
dimension(:,:,:),
pointer:: d, e
110 real(kind=KindScalarVar),
dimension(size(d,1),size(d,2),size(d,3)):: pg
112 pg = (self%gamma-1.0)*d*e
117 FUNCTION temperature (self, d, e)
RESULT (tt)
119 real,
dimension(:,:,:),
pointer:: d, e
120 real,
dimension(size(d,1),size(d,2),size(d,3)):: tt
122 tt = (self%gamma-1.0)*e
128 SUBROUTINE lookup (self, dim, lnd, ee, lnx, x, lny, y, pg, tt, ss, rk, src, gamma)
131 real,
optional:: gamma
132 real,
dimension(:,:,:),
intent(in),
pointer,
optional :: lnx, x, lny, y, lnd, ee
133 real,
dimension(:,:,:),
optional :: pg, tt, ss
134 real,
dimension(:,:,:,:),
optional :: src, rk
135 real,
dimension(:,:,:),
pointer :: x_loc, y_loc, lnd_l, ee_l
137 integer,
save:: itimer=0
140 call trace%begin (
'eos_t%lookup', itimer=itimer)
141 if (
present(gamma))
then 147 else if (
present(lny))
then 148 allocate (y_loc(dim(1),dim(2),dim(3)))
151 call mpi%abort (
'eos_t: missing 2nd argument')
156 else if (
present(lnx))
then 157 allocate (x_loc(dim(1),dim(2),dim(3)))
160 call mpi%abort (
'eos_t: missing 1st argument')
163 if (
present(lnd).and.
present(ee))
then 164 if (
present(src))
then 165 call lookup_table (self, dim, ee=ee, lnd=lnd, rk=rk, src=src)
167 call lookup_table (self, dim, ee=ee, lnd=lnd, rk=rk)
169 else if (
present(src).and.
present(rk))
then 170 allocate (lnd_l(dim(1),dim(2),dim(3)))
171 allocate (ee_l(dim(1),dim(2),dim(3)))
174 call lookup_table (self, dim, ee=ee_l, lnd=lnd_l, src=src, rk=rk)
175 deallocate (lnd_l, ee_l)
177 if (
present(tt))
then 178 tt = (self%gamma-1.0)*x_loc
180 if (
present(src))
then 181 stefan = cgs%stefan/(scaling%p*scaling%u)*scaling%temp**4
183 src(:,:,:,i) = stefan*
real(((self%gamma-1.0)*x_loc),kind=8)**4
186 if (
present(pg))
then 187 pg = (self%gamma-1.0)*y_loc*x_loc
189 if (
present(rk))
then 191 rk(:,:,:,i) = y_loc*self%kappa
195 if (
present(lnx))
deallocate(x_loc)
196 if (
present(lny))
deallocate(y_loc)
197 call trace%end(itimer)
198 END SUBROUTINE lookup
203 SUBROUTINE lookup_table (self, dim, e, ee, d, lnd, lne, pg, tt, ne, rk, src, &
207 real,
dimension(:,:,:) ,
pointer,
optional :: lnd,d,lne,e,ee,tt,ne,ss
208 real,
dimension(:,:,:,:) ,
optional :: src,rk
209 real,
dimension(:,:,:) ,
optional :: pg
210 real,
optional:: gamma
212 real,
dimension(:,:,:),
pointer :: lnd_loc, ee_loc
213 integer,
dimension(dim(1)) :: np,np1,ik,ntab,ik1,ntab1
214 real(kind=4),
dimension(dim(1)) :: px,py,f00,f01,f10,f11, &
215 fx00,fx01,fx10,fx11,fy00,fy01,fy10,fy11
216 real :: rhm1,rhm2,drhm,algrk,eek,eek1,py1
217 real :: qx,qy,pxqx,pxpx,qxqx,pypy,pyqy,qyqy,pxpy,qxqy,pxqy,qxpy
219 integer :: ix,iy,iz,kee,kee1,j
220 integer :: mx,my,mz,mbox,nbelow_d, nbelow_e
221 integer,
save:: itimer=0
225 call trace%begin (
'eos_t%lookup_table', itimer=itimer)
226 if (io%verbose>2)
then 227 print *,
'present(lnd) =',
present(lnd)
228 print *,
'present(ee) =',
present(ee)
229 print *,
'present(pg) =',
present(pg)
230 print *,
'present(src) =',
present(src)
231 print *,
'present(rk) =',
present(rk)
238 if (
present(lnd))
then 240 else if (
present(d))
then 241 allocate (lnd_loc(dim(1),dim(2),dim(3)))
244 call mpi%abort (
'eos_t: missing density')
246 if (
present(ee))
then 248 else if (
present(e))
then 249 allocate (ee_loc(dim(1),dim(2),dim(3)))
253 ee_loc = e/exp(lnd_loc)
255 else if (
present(lne))
then 256 allocate (ee_loc(dim(1),dim(2),dim(3)))
259 call mpi%abort (
'eos_t: missing internal energy')
264 drhm=(rhm2-rhm1)/(md-1)
276 algrk=1.+(lnd_loc(ix,iy,iz)-rhm1)/drhm
277 nbelow_d = nbelow_d + merge(1,0,lnd_loc(ix,iy,iz)<rhm1)
278 np(ix)=max0(1,min0(md-1,int(algrk)))
285 ntab(ix) = mtab(np(ix))
286 eek = 1.+(ee_loc(ix,iy,iz)-eemin(np(ix)))/deetab(np(ix))
287 nbelow_e = nbelow_e + merge(1,0,eek<1.0)
289 kee = min(ntab(ix)-1,max(1,int(eek)))
291 ik(ix) = itab(np(ix))+kee-1
292 ntab1(ix) = mtab(np(ix)+1)
295 kee1 = kee - nint((eemin(np(ix)+1)-eemin(np(ix))) &
297 kee1 = min0(ntab1(ix)-1,max0(1,kee1))
298 ik1(ix) = itab(np(ix)+1)+kee1-1
307 pxpx = px(ix) * px(ix)
312 pypy = py(ix) * py(ix)
316 pxpy = px(ix) * py(ix)
320 f00(ix) = qxqy * (1. + pxqx - pxpx + pyqy - pypy)
321 f01(ix) = qxpy * (1. + pxqx - pxpx + pyqy - qyqy)
322 f10(ix) = pxqy * (1. - qxqx + pxqx + pyqy - pypy)
323 f11(ix) = pxpy * (1. - qxqx + pxqx + pyqy - qyqy)
325 fx00(ix) = qxqy * pxqx
326 fx01(ix) = qxpy * pxqx
327 fx10(ix) = - pxqy * pxqx
328 fx11(ix) = - pxpy * pxqx
330 fy00(ix) = qxqy * pyqy
331 fy01(ix) = - qxpy * pyqy
332 fy10(ix) = pxqy * pyqy
333 fy11(ix) = - pxpy * pyqy
339 if (
present(pg))
then 341 pg(ix,iy,iz) = exp( &
342 f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
343 + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
344 + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
345 + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
346 + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
347 + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
348 + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
349 + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
350 + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
351 + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
358 if (
present(rk))
then 360 rk(ix,iy,iz,1) = exp( &
361 f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
362 + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
363 + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
364 + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
365 + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
366 + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
367 + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
368 + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
369 + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
370 + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
373 if (
size(rk,4)>1)
then 375 rk(:,iy,iz,j) = rk(:,iy,iz,j-1)*10.0
382 if (
present(tt))
then 385 f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
386 + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
387 + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
388 + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
389 + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
390 + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
391 + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
392 + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
393 + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
394 + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
401 if (
present(ne))
then 403 ne(ix,iy,iz) = exp( &
404 f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
405 + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
406 + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
407 + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
408 + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
409 + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
410 + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
411 + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
412 + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
413 + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
418 if (
present(src))
then 421 ik(ix) = ik(ix) + 3*ntab(ix)
422 ik1(ix) = ik1(ix) + 3*ntab1(ix)
423 src(ix,iy,iz,j) = exp( &
424 f00(ix) * tab(ik(ix))+ f01(ix) * tab(ik(ix) + 1) &
425 + f10(ix) * tab(ik1(ix))+ f11(ix) * tab(ik1(ix) + 1) &
426 + fx00(ix) * tab(ik(ix) + 2 * ntab(ix) ) &
427 + fx01(ix) * tab(ik(ix) + 2 * ntab(ix) + 1) &
428 + fx10(ix) * tab(ik1(ix) + 2 * ntab1(ix) ) &
429 + fx11(ix) * tab(ik1(ix) + 2 * ntab1(ix) + 1) &
430 + fy00(ix) * tab(ik(ix) + ntab(ix) ) &
431 + fy01(ix) * tab(ik(ix) + ntab(ix) + 1) &
432 + fy10(ix) * tab(ik1(ix) + ntab1(ix) ) &
433 + fy11(ix) * tab(ik1(ix) + ntab1(ix) + 1) &
440 if (nbelow_d+nbelow_e > 0)
then 442 'WARNING: number of points below table d & e:', nbelow_d, nbelow_e
448 if (
present(e).or.
present(lne))
then 452 call trace%end (itimer)
453 END SUBROUTINE lookup_table
Define code units, in terms of (here) CGS units.
Equation of state module for any sort of tables, provided by a reader.
Fundamental constants in CGS and SI units.