Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module op_eval_mombarg
21 :
22 : use const_def, only: dp, pi
23 :
24 : implicit none
25 :
26 : private
27 : public :: compute_grad
28 : public :: compute_grad_fast
29 : public :: compute_gamma_grid
30 : public :: compute_kappa
31 : public :: compute_kappa_fast
32 : public :: interpolate_kappa
33 : public :: compute_kappa_grid
34 :
35 : real(dp), parameter :: log_c = 10.476820702927927d0 !log10_cr(dble(299792458e2)) c = speed of light
36 : real(dp), parameter :: log10_bohr_radius_sqr = -16.55280d0
37 : real(dp), parameter :: logNa = 23.779750912481397d0 !log10_cr(6.0221409d23) Avogadro's number
38 :
39 : real(dp), parameter :: diff_v = (1.0552976117319748d0 - 0.00010565516589892675d0) !v(u(1)) - v(u(nptot))
40 :
41 : contains
42 :
43 :
44 : !!! Compute g_rad for a single cell. (Currently not used.)
45 0 : subroutine compute_grad(k, fk, logT_face, logRho_face,l, r,lkap_ross_cell, lgrad_cell, ierr,ite,jne,epatom,amamu,sig,eumesh)
46 : ! OP mono data for: H, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, and Ni.
47 : use chem_def, only: chem_isos, ih1, ihe3, ihe4, ic12, in14, io16, ine20, ina23, &
48 : img24, ial27, isi28, is32, iar40, ica40, icr52, imn55, ife56, ini58
49 : use interp_2d_lib_sg
50 : use cubic_interpolator, only: interpolator
51 :
52 : integer, intent(inout) :: ierr
53 : integer, intent(in) :: k
54 : real(dp), intent(in) :: fk(:), logT_face, logRho_face, l, r
55 : integer :: nel, nptot
56 : parameter(nel = 17, nptot = 10000) !number of elements and number of u-mesh points.
57 : real(dp), intent(out) :: lkap_ross_cell, lgrad_cell(nel)
58 :
59 : integer , pointer :: ite(:),jne(:)
60 : real(dp), pointer :: sig(:,:,:)
61 : real(dp), pointer :: epatom(:,:),amamu(:),eumesh(:,:,:)
62 :
63 :
64 : ! ite: temperature index, logT = 0.025*ite, (1648)
65 : ! jne: electron density index, logNe = 0.25*jne, (1648)
66 : ! epatom: electrons per atom, (nel,nptot).
67 : ! amamu: mu = mean atomic weight (nel)
68 : ! sig: monochromatic cross sections in (a0^2), (nel,1648,nptot)
69 : ! eumesh: sigma_i*(1 - exp(-u(v))) - a_k,i, where u=h*nu/(kT). OP mono grid is equally spaced in variable v. a_k,i are the correction factors. (nel,1648,nptot)
70 :
71 :
72 : integer :: ke, id, m, ik, i
73 :
74 0 : real(dp):: epa_mix_cell(1648), amu_mix_cell, logRho(1648),logT(1648) ! Number of electrons per atom, mean molecular weight, density and temperature as a function of ite (temp index) and jne (density index) from the OP mono data.
75 0 : real(dp) :: delta(1648)
76 0 : real(dp) :: lgamm_cell(nel) !interpolated log kappa Rossland and log gamma_k in each cell (k = element index).
77 : real(dp), parameter :: dv = diff_v/nptot !v(u(1)) - v(u(nptot))/nptot
78 0 : real(dp) :: mH
79 :
80 : integer :: delta_min_idx
81 0 : real(dp) :: lkap_Ross(4,4),gamma_k(4,4,nel),lgamm(4,4,nel),sig_Ross(4,4) ! Rossland cross section, log kappa, gamma_k, log gamma_k for interpolation.
82 0 : real(dp) :: sf, flux !'scale factor' for g_rad, local flux.
83 0 : real, allocatable :: sig_mix_cell(:,:,:),sig_int(:)
84 : integer :: ii, jj, ite_min, jne_min, ii_min, jj_min, ite_step, jne_step
85 : integer :: ite_i, jne_i, dite, djne, i_grid(4,4)
86 0 : real(dp) :: logT_min, logRho_min, logT_grid(4,4), logRho_grid(4,4)
87 : integer :: offset1, offset2, tries, missing_point(4,4)
88 0 : real(dp) :: log_amu_mix_cell, gam
89 : integer :: imin, imax
90 : logical :: retry
91 :
92 : type(interpolator) :: rossl_interpolator
93 : type(interpolator) :: gaml_interpolators(nel)
94 :
95 0 : ierr = 0
96 :
97 :
98 0 : imin = 1 !-1
99 0 : imax = 1648 !-1
100 0 : ite_step = 2
101 0 : jne_step = 2
102 :
103 : !!! Compute the number of electrons per atom for the local mixture.
104 0 : epa_mix_cell = 0d0
105 0 : do i=imin,imax
106 0 : epa_mix_cell(i) = dot_product(fk,epatom(1:nel,i))
107 : end do
108 :
109 0 : amu_mix_cell = dot_product(fk,amamu)
110 :
111 0 : mH = chem_isos% W(ih1) * 1.660538782d-24
112 0 : log_amu_mix_cell = log10(amu_mix_cell * mH)
113 : !logRho = 0.25d0*jne + log_amu_mix_cell - log10(epa_mix_cell)
114 0 : logRho = 0.25d0*jne + log10(amu_mix_cell) - log10(epa_mix_cell) -logNa
115 0 : logT = 0.025d0*ite
116 : !write(*,*) 'rho computed'
117 :
118 : !!! Select nearest points in logT,logRho for interpolation.
119 : !!! First, find the nearest OP data point, and check which of the four possible positionings this minimum has wrt to (T,Rho)_cell.
120 : !!! Acquire the remaining 15 points of the 4x4 grid, where (T,Rho)_cell is located in the inner square.
121 :
122 0 : delta = sqrt((logRho - logRho_face)*(logRho - logRho_face)/0.25d0 + (logT-logT_face)*(logT-logT_face)/0.025d0)
123 :
124 0 : delta_min_idx = MINLOC(delta, DIM=1)
125 0 : ite_min = ite(delta_min_idx)
126 0 : jne_min = jne(delta_min_idx)
127 0 : logT_min = logT(delta_min_idx)
128 0 : logRho_min = logRho(delta_min_idx)
129 0 : if (logT_min <= logT_face .and. logRho_min <= logRho_face) then
130 0 : ii_min = 2
131 0 : jj_min = 2
132 0 : else if (logT_min < logT_face .and. logRho_min > logRho_face) then
133 0 : ii_min = 2 !3
134 0 : jj_min = 3 !2
135 0 : else if (logT_min > logT_face .and. logRho_min < logRho_face) then
136 0 : ii_min = 3
137 0 : jj_min = 2
138 0 : else if (logT_min > logT_face .and. logRho_min > logRho_face) then
139 0 : ii_min = 3
140 0 : jj_min = 3
141 : end if
142 : !write(*,*) 'ii, jj min', ii_min, jj_min, logT_min, XC, logRho_min, logRho_cell
143 :
144 0 : offset1 = 0
145 0 : offset2 = 0
146 0 : missing_point = 1
147 : tries = 0
148 : retry = .true.
149 : !do while (point_not_found .ne. 0) !If (T,Rho)_cell lies to close to the edge of the OP mono grid (in Rho), shift the 4x4 grid up/down by 1 in jne.
150 0 : do while (retry) !If (T,Rho)_cell lies to close to the edge of the OP mono grid (in Rho), shift the 4x4 grid up/down by 1 in jne.
151 0 : missing_point = 1
152 0 : retry = .false.
153 0 : do jj=1,4
154 0 : do ii=1,4
155 0 : dite = (ii - ii_min + offset1)*ite_step !(ii - ii_min)*ite_step + offset2
156 0 : djne = (jj - jj_min + offset2)*jne_step !(jj - jj_min)*jne_step + offset1
157 0 : do i =imin,imax
158 0 : ite_i = ite(i)
159 0 : jne_i = jne(i)
160 :
161 0 : if (ite_i == ite_min + dite .and. jne_i == jne_min + djne) THEN
162 0 : logT_grid(ii,jj) = logT(i)
163 0 : logRho_grid(ii,jj) = logRho(i)
164 0 : i_grid(ii,jj) = i
165 : !write(*,*) ite_i, jne_i, ii, jj, i_grid(ii,jj), ',',logT_grid(ii,jj),&
166 : !',', logRho_grid(ii,jj)
167 0 : missing_point(ii,jj) = 0
168 : end if
169 : end do
170 : end do
171 : end do
172 0 : if (SUM(missing_point) > 0) then
173 0 : retry = .true.
174 0 : if (SUM(missing_point(2:4,1:3)) == 0) then
175 0 : offset1 = offset1 + 1
176 0 : offset2 = offset2 - 1
177 0 : else if (SUM(missing_point(1:3,1:3)) == 0) then
178 0 : offset1 = offset1 - 1
179 0 : offset2 = offset2 - 1
180 0 : else if (SUM(missing_point(2:4,2:4)) == 0) then
181 0 : offset1 = offset1 + 1
182 0 : offset2 = offset2 + 1
183 0 : else if (SUM(missing_point(1:3,2:4)) == 0) then
184 0 : offset1 = offset1 - 1
185 0 : offset2 = offset2 + 1
186 : else
187 : !write(*,*) 'warning interpolation quality in cell grad', k,ite_min,jne_min, &
188 : !logT_face, logRho_face, ii_min, jj_min
189 0 : if (ii_min == 3 .and. jj_min == 2) then
190 0 : offset1 = offset1 + 2
191 0 : offset2 = offset2 - 2
192 0 : else if (ii_min == 2 .and. jj_min == 3) then
193 0 : offset1 = offset1 - 2
194 0 : offset2 = offset2 + 2
195 0 : else if (ii_min == 2 .and. jj_min == 2) then
196 0 : offset1 = offset1 - 2
197 0 : offset2 = offset2 - 2
198 0 : else if (ii_min == 3 .and. jj_min == 3) then
199 0 : offset1 = offset1 + 2
200 0 : offset2 = offset2 - 1
201 : end if
202 : end if
203 : end if
204 :
205 :
206 0 : if (tries > 2) THEN ! To prevent loop from getting stuck.
207 0 : write(*,*) 'Cannot find points for interpolation compute_grad', &
208 0 : k, ite_min, jne_min, logT_min, logRho_min, logT_face, &
209 0 : logRho_face, missing_point, ii_min, jj_min, offset1, offset2, &
210 0 : imin, imax, log_amu_mix_cell
211 0 : ierr = 1
212 : return
213 : end if
214 0 : tries = tries + 1
215 : end do
216 : !write(*,*) 'Points for interpolation selected', k
217 :
218 : !!! Compute the monochromatic cross-section for the local mixture.
219 0 : allocate(sig_mix_cell(4,4,nptot),sig_int(nptot-1), stat=ierr)
220 0 : sig_mix_cell = 0d0
221 0 : do jj=1,4
222 0 : do ii=1,4
223 0 : ik = i_grid(ii,jj) !+ 1648*(ke-1)
224 0 : do m=1,nptot
225 0 : sig_mix_cell(ii,jj,m) = dot_product(fk,sig(:,ik,m))
226 : end do
227 : end do
228 : end do
229 :
230 : !!! Compute the Rossland mean cross-section by integrating over variable v (mesh equally spaced in v).
231 0 : sig_Ross = 0d0
232 0 : do jj=1,4
233 0 : do ii=1,4
234 0 : sig_int(1:nptot-1) = (1/sig_mix_cell(ii,jj,1:nptot-1) + 1/sig_mix_cell(ii,jj,2:nptot))/2.0d0 * dv
235 0 : do m=1, nptot-1
236 : !sig_Ross(ii,jj) = sig_Ross(ii,jj) + (1/sig_mix_cell(ii,jj,m) + 1/sig_mix_cell(ii,jj,m+1))/2.0d0 * dv
237 0 : sig_Ross(ii,jj) = sig_Ross(ii,jj) + sig_int(m)
238 : end do
239 : end do
240 : end do
241 :
242 0 : lkap_Ross = log10_bohr_radius_sqr - log_amu_mix_cell - log10(sig_Ross)
243 :
244 0 : call rossl_interpolator% initialize()
245 0 : do jj = 1, 4
246 0 : do ii = 1, 4
247 0 : call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_Ross(ii,jj))
248 : end do
249 : end do
250 :
251 0 : lkap_ross_cell = rossl_interpolator% evaluate(logT_face,logRho_face)
252 :
253 :
254 : !!! Compute gamma for each element for interpolation.
255 0 : gamma_k = 0d0
256 :
257 0 : do jj=1,4
258 0 : do ii=1,4
259 0 : ik = i_grid(ii,jj) ! + 1648*(ke-1)
260 0 : do ke=3,nel
261 : gam = 0
262 0 : do m=1, nptot-1
263 0 : gam = gam + ((eumesh(ke,ik,m)/ sig_mix_cell(ii,jj,m)) +(eumesh(ke,ik,m+1)/ sig_mix_cell(ii,jj,m+1)))
264 : end do
265 0 : gamma_k(ii,jj,ke) = gam/2.0d0 * dv
266 : end do
267 : end do
268 : end do
269 :
270 0 : deallocate(sig_mix_cell,sig_int)
271 :
272 0 : where (gamma_k < 0.0d0)
273 : gamma_k = 10d-30
274 : endwhere
275 0 : lgamm = log10(gamma_k)
276 :
277 0 : do ke = 3, nel
278 0 : call gaml_interpolators(ke)% initialize()
279 0 : do jj = 1, 4
280 0 : do ii = 1, 4
281 0 : call gaml_interpolators(ke)% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lgamm(ii,jj,ke))
282 : end do
283 : end do
284 :
285 0 : lgamm_cell(ke) = gaml_interpolators(ke)% evaluate(logT_face, logRho_face) !cell
286 : end do
287 0 : lgamm_cell(1) = -30.
288 0 : lgamm_cell(2) = -30.
289 : !write(*,*) 'lgamm_cell computed'!, ke, lgamm_cell(ke)
290 :
291 : ! Compute log g_rad.
292 0 : flux = l / (4d0*pi*r*r)
293 0 : sf = lkap_ross_cell - log_c + log10(flux) + log10(amu_mix_cell)
294 0 : lgrad_cell = sf + lgamm_cell - log10(amamu)
295 0 : lgrad_cell(1) = -30.
296 0 : lgrad_cell(2) = -30.
297 : !write(*,*) 'log kappa + log g_rad Fe cell', k, lkap_ross_cell, lgrad_cell(16)
298 : !write(*,*) 'compute_grad done', k
299 :
300 0 : end subroutine compute_grad
301 :
302 :
303 : !!! Compute gamma factors and kappa_Ross for all OP mono data points for a given mixture.
304 0 : subroutine compute_gamma_grid(ngp, fk_all,lgamm_pcg, lkap_face_pcg, &
305 : logT_pcg, logRho_pcg, ierr,ite,jne,epatom,amamu,sig,eumesh)
306 : ! OP mono data for: H, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, and Ni.
307 0 : use chem_def, only: chem_isos, ih1, ihe3, ihe4, ic12, in14, io16, ine20, ina23, &
308 : img24, ial27, isi28, is32, iar40, ica40, icr52, imn55, ife56, ini58
309 : use interp_2d_lib_sg
310 : use cubic_interpolator, only: interpolator
311 :
312 : integer, intent(inout) :: ierr
313 : integer, intent(in) :: ngp
314 : real(dp), intent(in) :: fk_all(:)
315 : integer :: nel, nptot
316 : parameter(nel = 17, nptot = 10000) !number of elements and number of u-mesh points.
317 : !real(dp), pointer, intent(out)::lgamm_pcg(:,:,:),lkap_face_pcg(:,:)
318 : real(dp), intent(out)::lgamm_pcg(nel,1648),lkap_face_pcg(1648)
319 :
320 : real(dp), intent(out) :: logT_pcg(1648), logRho_pcg(1648)
321 :
322 : integer , pointer :: ite(:),jne(:)
323 : real(dp), pointer :: sig(:,:,:) !, ak_f(:,:)
324 : real(dp), pointer :: epatom(:,:),amamu(:),eumesh(:,:,:)
325 : integer :: imin, imax
326 :
327 : integer :: n, ke, nz, id, m, ik, i, j
328 :
329 : !real(dp) :: fk_norm_fac !Local fractional abundance per element and normalization factor.
330 0 : real(dp):: epa_mix_cell(1648), amu_mix_cell, fk(nel) ! Number of electrons per atom, mean molecular weight, density and temperature as a function of ite (temp index) and jne (density index) from the OP mono data.
331 : !integer :: eid(nel)
332 : real(dp), parameter :: dv = diff_v/nptot !v(u(1)) - v(u(nptot))/nptot
333 0 : real(dp) :: mH
334 :
335 : !!!! For interpolator.
336 :
337 :
338 0 : real(dp) :: sig_Ross(1648), gamma_k(nel,1648) !, lgamm(nel,1648) ! Rossland cross section, log kappa, gamma_k, log gamma_k for interpolation.
339 :
340 :
341 0 : real, allocatable :: sig_mix_cell(:,:),sig_int(:)
342 0 : real(dp) :: log_amu_mix_cell, gam
343 :
344 :
345 0 : ierr = 0
346 :
347 : !!! Compute an estimated temperature range.
348 0 : imin = 1 !-1
349 0 : imax = 1648 !-1
350 :
351 0 : fk = fk_all !(j,:)
352 : !!! Compute the number of electrons per atom for the local mixture.
353 0 : epa_mix_cell = 0d0
354 0 : do i=imin,imax
355 0 : epa_mix_cell(i) = dot_product(fk,epatom(1:nel,i))
356 : end do
357 :
358 0 : amu_mix_cell = dot_product(fk,amamu)
359 :
360 0 : mH = chem_isos% W(ih1) * 1.660538782d-24
361 0 : log_amu_mix_cell = log10(amu_mix_cell * mH)
362 :
363 :
364 : !!! Compute the monochromatic cross-section for the local mixture.
365 0 : allocate(sig_mix_cell(1648,nptot),sig_int(nptot-1), stat=ierr)
366 0 : sig_mix_cell = 0d0
367 0 : !$OMP PARALLEL DO PRIVATE(i,m) SCHEDULE(guided)
368 : do i=1,1648
369 : do m=1,nptot
370 : sig_mix_cell(i,m) = dot_product(fk,sig(:,i,m))
371 : end do
372 : end do
373 : !$OMP END PARALLEL DO
374 :
375 : !!! Compute the Rossland mean cross-section by integrating over variable v (mesh equally spaced in v).
376 0 : sig_Ross = 0d0
377 0 : !$OMP PARALLEL DO PRIVATE(i,m,sig_int) SCHEDULE(guided)
378 : do i=1,1648
379 : sig_int(1:nptot-1) = (1/sig_mix_cell(i,1:nptot-1) + 1/sig_mix_cell(i,2:nptot))/2.0d0 * dv !inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot)
380 : do m=1, nptot-1
381 : sig_Ross(i) = sig_Ross(i) + sig_int(m)
382 : end do
383 : end do
384 : !$OMP END PARALLEL DO
385 :
386 :
387 0 : logT_pcg = 0.025d0*ite
388 0 : logRho_pcg = 0.25d0*jne + log10(amu_mix_cell) - log10(epa_mix_cell) - logNa
389 : !if (j==1) allocate(lkap_face_pcg(ngp,1648),stat=ierr)
390 0 : lkap_face_pcg = log10_bohr_radius_sqr - log_amu_mix_cell - log10(sig_Ross)
391 0 : !$OMP PARALLEL DO PRIVATE(i,ke,m,gam) SCHEDULE(guided)
392 : do i=imin,imax
393 : do ke=3,nel
394 : gam = 0d0
395 : do m=1, nptot-1
396 : gam = gam + ((eumesh(ke,i,m)/ sig_mix_cell(i,m))+(eumesh(ke,i,m+1)/ sig_mix_cell(i,m+1)))
397 : end do
398 : gamma_k(ke,i) = gam/2.0d0 * dv
399 : end do
400 : end do
401 : !$OMP END PARALLEL DO
402 :
403 0 : where (gamma_k < 0.0d0)
404 : gamma_k = 10d-30
405 : endwhere
406 :
407 0 : deallocate(sig_mix_cell,sig_int)
408 :
409 0 : lgamm_pcg = log10(gamma_k)
410 0 : end subroutine compute_gamma_grid
411 :
412 :
413 : !!! Compute g_rad from precomputeds for grid for gamma and kappa_Ross for the entire OP mono data.
414 0 : subroutine compute_grad_fast(k, fk, logT_face, logRho_face, l, r,lgrad_cell, &
415 : ierr,ite,jne,epatom,amamu,logT_pcg,logRho_pcg,lgamm_pcg,lkap_face_pcg)
416 :
417 0 : use cubic_interpolator, only: interpolator
418 :
419 : integer, intent(inout) :: ierr
420 : integer, intent(in) :: k
421 : real(dp), intent(in) :: fk(:), logT_face, logRho_face, l, r
422 : real(dp), intent(in) :: logT_pcg(1648), logRho_pcg(1648) !, lgamm_pcg(:,:), lkap_face_pcg(:)
423 : integer :: nel, nptot
424 : parameter(nel = 17, nptot = 10000) !number of elements and number of u-mesh points.
425 : real(dp), intent(in) :: lgamm_pcg(nel,1648), lkap_face_pcg(1648)
426 : real(dp), intent(out) :: lgrad_cell(nel)
427 :
428 : integer , pointer :: ite(:),jne(:)
429 : real(dp), pointer :: epatom(:,:),amamu(:)
430 :
431 : integer :: n, ke, nz, id, m, ik, i
432 :
433 0 : real(dp):: epa_mix_cell(1648), amu_mix_cell ! Number of electrons per atom, mean molecular weight, density and temperature as a function of ite (temp index) and jne (density index) from the OP mono data.
434 0 : real(dp) :: delta(1648)
435 : integer :: eid(nel)
436 0 : real(dp) :: lgamm_cell(nel) !interpolated log kappa Rossland and log gamma_k in each cell (k = element index).
437 : real(dp), parameter :: dv = diff_v/nptot !v(u(1)) - v(u(nptot))/nptot
438 : real(dp) :: mH
439 :
440 : !!!! For interpolator.
441 : integer :: delta_min_idx
442 0 : real(dp) :: lkap_Ross(4,4),sig_Ross(4,4) ! Rossland cross section, log kappa, gamma_k, log gamma_k for interpolation.
443 0 : real(dp) :: sf, flux !'scale factor' for g_rad, local flux.
444 :
445 : integer :: ii, jj, ite_min, jne_min, ii_min, jj_min, ite_step, jne_step
446 : integer :: ite_i, jne_i, dite, djne, i_grid(4,4)
447 0 : real(dp) :: logT_min, logRho_min, logT_grid(4,4),logRho_grid(4,4),lgamm(4,4,nel)
448 : integer :: offset1, offset2, tries, missing_point(4,4)
449 0 : real(dp) :: log_amu_mix_cell, lkap_ross_cell
450 : integer :: imin, imax
451 0 : real(dp) :: logT(1648), logRho(1648) !, ite_grid(4,4), jne_grid(4,4)
452 : logical :: retry, do_difficult_point
453 : type(interpolator) :: rossl_interpolator
454 : type(interpolator) :: gaml_interpolators(nel)
455 :
456 0 : ierr = 0
457 :
458 0 : imin = 1
459 0 : imax = 1648
460 :
461 :
462 :
463 0 : amu_mix_cell = dot_product(fk,amamu)
464 :
465 0 : logT = logT_pcg !0.025d0*ite
466 0 : logRho = logRho_pcg !0.25d0*jne + log10_cr(amu_mix_cell) - log10(epa_mix_cell) - logNa
467 :
468 : !!! Compute an estimated temperature range.
469 : imin = 1
470 : imax = 1648
471 0 : ite_step = 2
472 0 : jne_step = 2
473 :
474 0 : delta = sqrt((logRho - logRho_face)*(logRho - logRho_face)/0.25d0 +(logT-logT_face)*(logT-logT_face)/0.025d0)
475 :
476 0 : delta_min_idx = MINLOC(delta, DIM=1)
477 0 : ite_min = ite(delta_min_idx)
478 0 : jne_min = jne(delta_min_idx)
479 0 : logT_min = logT(delta_min_idx)
480 0 : logRho_min = logRho(delta_min_idx)
481 0 : if (logT_min <= logT_face .and. logRho_min <= logRho_face) then
482 0 : ii_min = 2
483 0 : jj_min = 2
484 0 : else if (logT_min < logT_face .and. logRho_min > logRho_face) then
485 0 : ii_min = 2
486 0 : jj_min = 3
487 0 : else if (logT_min > logT_face .and. logRho_min < logRho_face) then
488 0 : ii_min = 3
489 0 : jj_min = 2
490 0 : else if (logT_min > logT_face .and. logRho_min > logRho_face) then
491 0 : ii_min = 3
492 0 : jj_min = 3
493 : end if
494 :
495 :
496 0 : offset1 = 0
497 0 : offset2 = 0
498 0 : missing_point = 1
499 : tries = 0
500 : retry = .true.
501 0 : do while (retry) !If (T,Rho)_cell lies to close to the edge of the OP mono grid (in Rho), shift the 4x4 grid up/down by 1 in jne.
502 0 : missing_point = 1
503 0 : retry = .false.
504 0 : do jj=1,4
505 0 : do ii=1,4
506 0 : dite = (ii - ii_min + offset1)*ite_step
507 0 : djne = (jj - jj_min + offset2)*jne_step
508 0 : do i =imin,imax
509 0 : ite_i = ite(i)
510 0 : jne_i = jne(i)
511 :
512 0 : if (ite_i == ite_min + dite .and. jne_i == jne_min + djne) THEN
513 0 : logT_grid(ii,jj) = logT(i)
514 0 : logRho_grid(ii,jj) = logRho(i)
515 0 : i_grid(ii,jj) = i
516 0 : missing_point(ii,jj) = 0
517 :
518 : end if
519 : end do
520 : end do
521 : end do
522 :
523 0 : if (SUM(missing_point) > 0) then
524 0 : retry = .true.
525 0 : if (SUM(missing_point(2:4,1:3)) == 0) then
526 0 : offset1 = offset1 + 1
527 0 : offset2 = offset2 - 1
528 0 : else if (SUM(missing_point(1:3,1:3)) == 0) then
529 0 : offset1 = offset1 - 1
530 0 : offset2 = offset2 - 1
531 0 : else if (SUM(missing_point(2:4,2:4)) == 0) then
532 0 : offset1 = offset1 + 1
533 0 : offset2 = offset2 + 1
534 0 : else if (SUM(missing_point(1:3,2:4)) == 0) then
535 0 : offset1 = offset1 - 1
536 0 : offset2 = offset2 + 1
537 : else
538 0 : if (ii_min == 3 .and. jj_min == 2) then
539 0 : offset1 = offset1 + 2
540 0 : offset2 = offset2 - 2
541 0 : else if (ii_min == 2 .and. jj_min == 3) then
542 0 : offset1 = offset1 - 2
543 0 : offset2 = offset2 + 2
544 0 : else if (ii_min == 2 .and. jj_min == 2) then
545 0 : offset1 = offset1 - 2
546 0 : offset2 = offset2 - 2
547 0 : else if (ii_min == 3 .and. jj_min == 3) then
548 0 : offset1 = offset1 + 2
549 0 : offset2 = offset2 - 1
550 : end if
551 : end if
552 : end if
553 :
554 0 : if (tries > 2) THEN ! To prevent loop from getting stuck.
555 0 : do_difficult_point = .false.
556 0 : if (ite_min == 226 .and. jne_min == 90 .and. ii_min == 2 .and. jj_min == 3) do_difficult_point = .true.
557 0 : if (ite_min == 226 .and. jne_min == 88 .and. ii_min == 2 .and. jj_min == 2) do_difficult_point = .true.
558 0 : if (ite_min == 230 .and. jne_min == 88 .and. ii_min == 3 .and. jj_min == 2) do_difficult_point = .true.
559 0 : if (ite_min == 230 .and. jne_min == 90 .and. ii_min == 3 .and. jj_min == 3) do_difficult_point = .true.
560 :
561 0 : if (do_difficult_point) then
562 0 : do jj=1,4
563 0 : do ii =1,4
564 0 : dite = (ii - ii_min)*2*ite_step
565 0 : djne = (jj - jj_min - 1)*jne_step
566 0 : do i=imin,imax
567 0 : ite_i = ite(i)
568 0 : jne_i = jne(i)
569 0 : if (ite_i == ite_min + dite .and. jne_i == jne_min + djne) then
570 0 : logT_grid(ii,jj) = logT(i)
571 0 : logRho_grid(ii,jj) = logRho(i)
572 0 : i_grid(ii,jj) = i
573 :
574 : end if
575 : end do
576 : end do
577 : end do
578 : retry = .false.
579 : else
580 :
581 0 : write(*,*) 'Cannot find points for interpolation compute_grad_fast', &
582 0 : k, ite_min, jne_min, logT_min, logRho_min, logT_face, logRho_face, &
583 0 : missing_point, ii_min, jj_min, offset1, offset2,imin, imax
584 0 : ierr = 1
585 0 : return
586 : end if
587 : end if
588 0 : tries = tries + 1
589 : end do
590 :
591 0 : do jj=1,4
592 0 : do ii=1,4
593 0 : ik = i_grid(ii,jj)
594 0 : do ke=3,nel
595 0 : lgamm(ii,jj,ke) = lgamm_pcg(ke,ik)
596 : end do
597 0 : lkap_Ross(ii,jj) = lkap_face_pcg(ik)
598 : end do
599 : end do
600 :
601 :
602 0 : call rossl_interpolator% initialize()
603 0 : do jj = 1, 4
604 0 : do ii = 1, 4
605 0 : call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_Ross(ii,jj))
606 : end do
607 : end do
608 :
609 0 : lkap_ross_cell = rossl_interpolator% evaluate(logT_face,logRho_face)
610 :
611 0 : do ke = 3, nel
612 0 : call gaml_interpolators(ke)% initialize()
613 0 : do jj = 1, 4
614 0 : do ii = 1, 4
615 0 : call gaml_interpolators(ke)% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lgamm(ii,jj,ke))
616 : end do
617 : end do
618 :
619 0 : lgamm_cell(ke) = gaml_interpolators(ke)% evaluate(logT_face, logRho_face) !cell
620 : end do
621 0 : lgamm_cell(1) = -30d0
622 0 : lgamm_cell(2) = -30d0
623 :
624 0 : flux = l / (4d0*pi*r*r)
625 0 : sf = lkap_ross_cell - log_c + log10(flux) + log10(amu_mix_cell)
626 0 : lgrad_cell = sf + lgamm_cell - log10(amamu)
627 0 : lgrad_cell(1) = -30.
628 0 : lgrad_cell(2) = -30.
629 :
630 : end subroutine compute_grad_fast
631 :
632 :
633 : !!! Compute Rossland opacity and derivatives for a single cell.
634 0 : subroutine compute_kappa(k,fk,logT_cntr,logRho_cntr,lkap_ross_cell, &
635 : dlnkap_rad_dlnT,dlnkap_rad_dlnRho,ierr,ite,jne,epatom,amamu,sig,logT_grid,logRho_grid,lkap_Ross)
636 : ! OP mono data for: H, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, and Ni.
637 : use chem_def, only: chem_isos, ih1, ihe3, ihe4, ic12, in14, io16, ine20, ina23, &
638 : img24, ial27, isi28, is32, iar40, ica40, icr52, imn55, ife56, ini58
639 : use interp_2d_lib_sg
640 : use cubic_interpolator, only: interpolator
641 :
642 : integer, intent(inout) :: ierr
643 : integer, intent(in) :: k
644 : real(dp), intent(in) :: fk(:), logT_cntr, logRho_cntr
645 : integer :: nel, nptot
646 : parameter(nel = 17, nptot = 10000) !number of elements and number of u-mesh points.
647 : real(dp), intent(out) :: lkap_ross_cell, dlnkap_rad_dlnT, dlnkap_rad_dlnRho
648 : real(dp), intent(out), dimension(4,4) :: logT_grid, logRho_grid, lkap_Ross
649 :
650 :
651 : integer , pointer :: ite(:),jne(:)
652 : real(dp), pointer :: sig(:,:,:)
653 : real(dp), pointer :: epatom(:,:),amamu(:)
654 :
655 : integer :: n, ke, nz, id, m, ik, i
656 :
657 0 : real(dp):: epa_mix_cell(1648), amu_mix_cell, logRho(1648), logT(1648) ! Number of electrons per atom, mean molecular weight, density and temperature as a function of ite (temp index) and jne (density index) from the OP mono data.
658 0 : real(dp) :: delta(1648)
659 : real(dp) :: lgamm_cell(nel) !interpolated log kappa Rossland and log gamma_k in each cell (k = element index).
660 : real(dp), parameter :: dv = diff_v/nptot !v(u(1)) - v(u(nptot))/nptot
661 0 : real(dp) :: mH
662 :
663 :
664 : integer :: delta_min_idx
665 0 : real(dp) :: sig_Ross(4,4) !,lkap_Ross(4,4), ! Rossland cross section, log kappa, gamma_k, log gamma_k for interpolation.
666 : real(dp) :: sf, flux !'scale factor' for g_rad, local flux.
667 :
668 0 : real, allocatable :: sig_mix_cell(:,:,:),sig_int(:)
669 : integer :: ii, jj, ite_min, jne_min, ii_min, jj_min, ite_step, jne_step
670 : integer :: ite_i, jne_i, dite, djne, i_grid(4,4)
671 : real(dp) :: logT_min, logRho_min !, logT_grid(4,4), logRho_grid(4,4)
672 : integer :: offset1, offset2, tries, missing_point(4,4)
673 : real(dp) :: log_amu_mix_cell
674 : integer :: imin, imax
675 : real(dp) :: sig_grid(nel,4,4,nptot)
676 : logical :: retry
677 :
678 : type(interpolator) :: rossl_interpolator
679 :
680 0 : ierr = 0
681 :
682 0 : imin = 1 !-1
683 0 : imax = 1648 !-1
684 0 : ite_step = 2
685 0 : jne_step = 2
686 :
687 : !!! Compute the number of electrons per atom for the local mixture.
688 0 : epa_mix_cell = 0d0
689 0 : do i=imin,imax
690 0 : epa_mix_cell(i) = dot_product(fk,epatom(1:nel,i))
691 : end do
692 :
693 :
694 0 : amu_mix_cell = dot_product(fk,amamu)
695 :
696 0 : mH = chem_isos% W(ih1) * 1.660538782d-24
697 0 : log_amu_mix_cell = log10(amu_mix_cell * mH)
698 : !logRho = 0.25d0*jne + log_amu_mix_cell - log10(epa_mix_cell)
699 0 : logRho = 0.25d0*jne + log10(amu_mix_cell) - log10(epa_mix_cell) -logNa
700 :
701 0 : logT = 0.025d0*ite
702 :
703 : !!! Select nearest points in logT,logRho for interpolation.
704 : !!! First, find the nearest OP data point, and check which of the four possible positionings this minimum has wrt to (T,Rho)_cell.
705 : !!! Acquire the remaining 15 points of the 4x4 grid, where (T,Rho)_cell is located in the inner square.
706 :
707 :
708 0 : delta = sqrt((logRho - logRho_cntr)*(logRho - logRho_cntr)/0.25d0 +(logT-logT_cntr)*(logT-logT_cntr)/0.025d0)
709 :
710 0 : delta_min_idx = MINLOC(delta, DIM=1)
711 0 : ite_min = ite(delta_min_idx)
712 0 : jne_min = jne(delta_min_idx)
713 0 : logT_min = logT(delta_min_idx)
714 0 : logRho_min = logRho(delta_min_idx)
715 0 : if (logT_min <= logT_cntr .and. logRho_min <= logRho_cntr) then
716 0 : ii_min = 2
717 0 : jj_min = 2
718 0 : else if (logT_min < logT_cntr .and. logRho_min > logRho_cntr) then
719 0 : ii_min = 2 !3
720 0 : jj_min = 3 !2
721 0 : else if (logT_min > logT_cntr .and. logRho_min < logRho_cntr) then
722 0 : ii_min = 3
723 0 : jj_min = 2
724 0 : else if (logT_min > logT_cntr .and. logRho_min > logRho_cntr) then
725 0 : ii_min = 3
726 0 : jj_min = 3
727 : end if
728 :
729 0 : offset1 = 0
730 0 : offset2 = 0
731 0 : missing_point = 1
732 : tries = 0
733 : retry = .true.
734 0 : do while (retry) !If (T,Rho)_cell lies to close to the edge of the OP mono grid (in Rho), shift the 4x4 grid up/down by 1 in jne.
735 0 : missing_point = 1
736 0 : retry = .false.
737 0 : do jj=1,4
738 0 : do ii=1,4
739 0 : dite = (ii - ii_min + offset1)*ite_step
740 0 : djne = (jj - jj_min + offset2)*jne_step
741 0 : do i =imin,imax
742 0 : ite_i = ite(i)
743 0 : jne_i = jne(i)
744 :
745 0 : if (ite_i == ite_min + dite .and. jne_i == jne_min + djne) THEN
746 0 : logT_grid(ii,jj) = logT(i)
747 0 : logRho_grid(ii,jj) = logRho(i)
748 0 : i_grid(ii,jj) = i
749 :
750 0 : missing_point(ii,jj) = 0
751 : end if
752 : end do
753 : end do
754 : end do
755 :
756 0 : if (SUM(missing_point) > 0) then
757 0 : retry = .true.
758 0 : if (SUM(missing_point(2:4,1:3)) == 0) then
759 0 : offset1 = offset1 + 1
760 0 : offset2 = offset2 - 1
761 0 : else if (SUM(missing_point(1:3,1:3)) == 0) then
762 0 : offset1 = offset1 - 1
763 0 : offset2 = offset2 - 1
764 0 : else if (SUM(missing_point(2:4,2:4)) == 0) then
765 0 : offset1 = offset1 + 1
766 0 : offset2 = offset2 + 1
767 0 : else if (SUM(missing_point(1:3,2:4)) == 0) then
768 0 : offset1 = offset1 - 1
769 0 : offset2 = offset2 + 1
770 : else
771 0 : if (ii_min == 3 .and. jj_min == 2) then
772 0 : offset1 = offset1 + 2
773 0 : offset2 = offset2 - 2
774 0 : else if (ii_min == 2 .and. jj_min == 3) then
775 0 : offset1 = offset1 - 2
776 0 : offset2 = offset2 + 2
777 0 : else if (ii_min == 2 .and. jj_min == 2) then
778 0 : offset1 = offset1 - 2
779 0 : offset2 = offset2 - 2
780 0 : else if (ii_min == 3 .and. jj_min == 3) then
781 0 : offset1 = offset1 + 2
782 0 : offset2 = offset2 - 1
783 : end if
784 : end if
785 : end if
786 :
787 :
788 :
789 0 : if (tries > 4) THEN ! To prevent loop from getting stuck.
790 0 : write(*,*) 'Cannot find points for interpolation compute_kappa', &
791 0 : k, ite_min, jne_min, logT_min, logRho_min,logT_cntr, logRho_cntr, &
792 0 : missing_point, ii_min, jj_min, offset1, offset2,imin, imax, log_amu_mix_cell
793 0 : ierr = 1
794 : return
795 : end if
796 0 : tries = tries + 1
797 : end do
798 :
799 :
800 : !!! Compute the monochromatic cross-section for the local mixture.
801 0 : allocate(sig_mix_cell(4,4,nptot),sig_int(nptot-1), stat=ierr)
802 :
803 :
804 0 : sig_mix_cell = 0d0
805 0 : do jj=1,4
806 0 : do ii=1,4
807 0 : ik = i_grid(ii,jj)
808 0 : do m=1,nptot
809 0 : sig_mix_cell(ii,jj,m) = dot_product(fk,sig(:,ik,m))
810 : end do
811 : end do
812 : end do
813 :
814 : !!! Compute the Rossland mean cross-section by integrating over variable v (mesh equally spaced in v).
815 0 : sig_Ross = 0d0
816 0 : do jj=1,4
817 0 : do ii=1,4
818 0 : sig_int = (1/sig_mix_cell(ii,jj,1:nptot-1) + 1/sig_mix_cell(ii,jj,2:nptot))/2.0d0 * dv !(1:nptot-1) inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot)
819 0 : do m=1, nptot-1
820 0 : sig_Ross(ii,jj) = sig_Ross(ii,jj) + sig_int(m)
821 : end do
822 : end do
823 : end do
824 :
825 0 : lkap_Ross = log10_bohr_radius_sqr - log_amu_mix_cell - log10(sig_Ross)
826 :
827 :
828 0 : call rossl_interpolator% initialize()
829 0 : do jj = 1, 4
830 0 : do ii = 1, 4
831 0 : call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_Ross(ii,jj))
832 : end do
833 : end do
834 :
835 0 : lkap_ross_cell = rossl_interpolator% evaluate(logT_cntr,logRho_cntr)
836 0 : dlnkap_rad_dlnT = rossl_interpolator% evaluate_deriv(logT_cntr,logRho_cntr, .true., .false.)
837 0 : dlnkap_rad_dlnRho = rossl_interpolator% evaluate_deriv(logT_cntr, logRho_cntr, .false., .true.)
838 :
839 :
840 0 : deallocate(sig_mix_cell,sig_int)
841 :
842 0 : end subroutine compute_kappa
843 :
844 :
845 :
846 : !!! Compute the Rossland opacity for the entire OP mono data set given a mixture.
847 0 : subroutine compute_kappa_grid(fk,lkap_ross_pcg, ierr,ite,jne,epatom,amamu,sig)
848 : ! OP mono data for: H, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, and Ni.
849 0 : use chem_def, only: chem_isos, ih1, ihe3, ihe4, ic12, in14, io16, ine20, ina23, &
850 : img24, ial27, isi28, is32, iar40, ica40, icr52, imn55, ife56, ini58
851 : use interp_2d_lib_sg
852 : use cubic_interpolator, only: interpolator
853 :
854 : integer, intent(inout) :: ierr
855 : real(dp), intent(in) :: fk(:)
856 : integer :: nel, nptot
857 : parameter(nel = 17, nptot = 10000) !number of elements and number of u-mesh points.
858 : real(dp), pointer, intent(out) :: lkap_ross_pcg(:) !, fk_pcg(:)
859 : real(dp) :: logT_pcg(1648), logRho_pcg(1648)
860 :
861 : integer , pointer :: ite(:),jne(:)
862 : real(dp), pointer :: sig(:,:,:) !, ak_f(:,:)
863 : real(dp), pointer :: epatom(:,:),amamu(:)
864 : integer :: imin, imax
865 :
866 : integer :: n, ke, nz, id, m, ik, i
867 :
868 0 : real(dp):: epa_mix_cell(1648), amu_mix_cell ! Number of electrons per atom, mean molecular weight, density and temperature as a function of ite (temp index) and jne (density index) from the OP mono data.
869 : real(dp) :: lgamm_cell(nel) !interpolated log kappa Rossland and log gamma_k in each cell (k = element index).
870 : real(dp), parameter :: dv = diff_v/nptot !v(u(1)) - v(u(nptot))/nptot
871 0 : real(dp) :: mH
872 :
873 : !!!! For interpolator.
874 : integer :: delta_min_idx
875 :
876 0 : real(dp) :: sig_Ross(1648) ! Rossland cross section, log kappa, gamma_k, log gamma_k for interpolation.
877 :
878 0 : real, allocatable :: sig_mix_cell(:,:),sig_int(:)
879 0 : real(dp) :: log_amu_mix_cell
880 :
881 : type(interpolator) :: rossl_interpolator
882 :
883 0 : ierr = 0
884 :
885 : !!! Compute an estimated temperature range.
886 0 : imin = 1
887 0 : imax = 1648
888 :
889 : !!! Compute the number of electrons per atom for the local mixture.
890 : epa_mix_cell = 0d0
891 0 : do i=imin,imax
892 : epa_mix_cell(i) = dot_product(fk,epatom(1:nel,i))
893 : end do
894 :
895 0 : amu_mix_cell = dot_product(fk,amamu)
896 :
897 0 : mH = chem_isos% W(ih1) * 1.660538782d-24
898 0 : log_amu_mix_cell = log10(amu_mix_cell * mH)
899 :
900 : !!! Compute the monochromatic cross-section for the local mixture.
901 0 : allocate(sig_mix_cell(1648,nptot),sig_int(nptot-1), stat=ierr)
902 0 : sig_mix_cell = 0d0
903 0 : !$OMP PARALLEL DO PRIVATE(i,m) SCHEDULE(guided)
904 :
905 : do i=1,1648
906 : do m=1,nptot
907 : sig_mix_cell(i,m) = dot_product(fk,sig(:,i,m))
908 : end do
909 : end do
910 : !$OMP END PARALLEL DO
911 :
912 : !!! Compute the Rossland mean cross-section by integrating over variable v (mesh equally spaced in v).
913 0 : sig_Ross = 0
914 0 : !$OMP PARALLEL DO PRIVATE(i,m,sig_int) SCHEDULE(guided)
915 : do i=1,1648
916 : sig_int(1:nptot-1) = (1/sig_mix_cell(i,1:nptot-1) + 1/sig_mix_cell(i,2:nptot))/2.0d0 * dv
917 : ! inv_sig_mix_cell(ii,jj,1:nptot-1) + inv_sig_mix_cell(ii,jj,2:nptot)
918 : do m=1, nptot-1
919 : sig_Ross(i) = sig_Ross(i) + sig_int(m)
920 : end do
921 : end do
922 : !$OMP END PARALLEL DO
923 :
924 0 : deallocate(sig_mix_cell,sig_int)
925 :
926 : logT_pcg = 0.025d0*ite
927 : logRho_pcg = 0.25d0*jne + log10(amu_mix_cell) - log10(epa_mix_cell) - logNa
928 0 : allocate(lkap_ross_pcg(1648),stat=ierr)
929 0 : lkap_ross_pcg = log10_bohr_radius_sqr - log_amu_mix_cell - log10(sig_Ross)
930 :
931 0 : end subroutine compute_kappa_grid
932 :
933 :
934 : !!! Routine to compute Rossland opacity from a precomputed grid of the entire OP mono data.
935 0 : subroutine compute_kappa_fast(k, fk, logT_cntr, logRho_cntr, lkap_ross_cell, dlnkap_rad_dlnT, dlnkap_rad_dlnRho, &
936 : ierr,ite,jne,epatom,amamu,sig,lkap_ross_pcg)
937 : ! OP mono data for: H, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Cr, Mn, Fe, and Ni.
938 0 : use chem_def, only: chem_isos, ih1, ihe3, ihe4, ic12, in14, io16, ine20, ina23, &
939 : img24, ial27, isi28, is32, iar40, ica40, icr52, imn55, ife56, ini58
940 : use cubic_interpolator, only: interpolator
941 :
942 : integer, intent(inout) :: ierr
943 : integer, intent(in) :: k
944 : real(dp), intent(in) :: fk(:), logT_cntr, logRho_cntr
945 : real(dp), pointer, intent(in) :: lkap_ross_pcg(:)
946 : integer :: nel, nptot
947 : parameter(nel = 17, nptot = 10000) !number of elements and number of u-mesh points.
948 : real(dp), intent(out) :: lkap_ross_cell, dlnkap_rad_dlnT, dlnkap_rad_dlnRho
949 : integer , pointer :: ite(:),jne(:)
950 : real(dp), pointer :: sig(:,:,:)
951 : real(dp), pointer :: epatom(:,:),amamu(:)
952 :
953 : integer :: n, ke, nz, id, m, ik, i
954 :
955 0 : real(dp):: epa_mix_cell(1648), amu_mix_cell ! Number of electrons per atom, mean molecular weight, density and temperature as a function of ite (temp index) and jne (density index) from the OP mono data.
956 0 : real(dp) :: delta(1648)
957 : real(dp) :: lgamm_cell(nel) !interpolated log kappa Rossland and log gamma_k in each cell (k = element index).
958 : real(dp), parameter :: dv = diff_v/nptot !v(u(1)) - v(u(nptot))/nptot
959 : real(dp) :: mH
960 :
961 : !!!! For interpolator.
962 : integer :: delta_min_idx
963 0 : real(dp) :: lkap_Ross(4,4),sig_Ross(4,4) ! Rossland cross section, log kappa, gamma_k, log gamma_k for interpolation.
964 : real(dp) :: sf, flux !'scale factor' for g_rad, local flux.
965 :
966 : integer :: ii, jj, ite_min, jne_min, ii_min, jj_min, ite_step, jne_step
967 : integer :: ite_i, jne_i, dite, djne, i_grid(4,4)
968 0 : real(dp) :: logT_min, logRho_min, logT_grid(4,4), logRho_grid(4,4)
969 : integer :: offset1, offset2, tries, missing_point(4,4)
970 : real(dp) :: log_amu_mix_cell
971 : integer :: imin, imax
972 0 : real(dp) :: logT(1648), logRho(1648) !, ite_grid(4,4), jne_grid(4,4)
973 : logical :: retry, do_difficult_point
974 :
975 : type(interpolator) :: rossl_interpolator
976 :
977 0 : ierr = 0
978 :
979 0 : imin = 1
980 0 : imax = 1648
981 :
982 : !!! Compute the number of electrons per atom for the local mixture.
983 0 : epa_mix_cell = 0d0
984 0 : do i=imin,imax
985 0 : epa_mix_cell(i) = dot_product(fk,epatom(:,i))
986 : end do
987 :
988 0 : amu_mix_cell = dot_product(fk,amamu)
989 :
990 0 : logT = 0.025d0*ite
991 0 : logRho = 0.25d0*jne + log10(amu_mix_cell) - log10(epa_mix_cell) - logNa
992 :
993 : imin = 1
994 : imax = 1648
995 0 : ite_step = 2
996 0 : jne_step = 2
997 :
998 : !!! Select nearest points in logT,logRho for interpolation.
999 : !!! First, find the nearest OP data point, and check which of the four possible positionings this minimum has wrt to (T,Rho)_cell.
1000 : !!! Acquire the remaining 15 points of the 4x4 grid, where (T,Rho)_cell is located in the inner square.
1001 :
1002 0 : delta = sqrt((logRho - logRho_cntr)*(logRho - logRho_cntr)/0.25d0 +(logT-logT_cntr)*(logT-logT_cntr)/0.025d0)
1003 :
1004 0 : delta_min_idx = MINLOC(delta, DIM=1)
1005 : !delta_min(1) = MINVAL(delta)(1)
1006 0 : ite_min = ite(delta_min_idx)
1007 0 : jne_min = jne(delta_min_idx)
1008 0 : logT_min = logT(delta_min_idx)
1009 0 : logRho_min = logRho(delta_min_idx)
1010 0 : if (logT_min <= logT_cntr .and. logRho_min <= logRho_cntr) then
1011 0 : ii_min = 2
1012 0 : jj_min = 2
1013 0 : else if (logT_min < logT_cntr .and. logRho_min > logRho_cntr) then
1014 0 : ii_min = 2 !3
1015 0 : jj_min = 3 !2
1016 0 : else if (logT_min > logT_cntr .and. logRho_min < logRho_cntr) then
1017 0 : ii_min = 3
1018 0 : jj_min = 2
1019 0 : else if (logT_min > logT_cntr .and. logRho_min > logRho_cntr) then
1020 0 : ii_min = 3
1021 0 : jj_min = 3
1022 : end if
1023 :
1024 :
1025 0 : offset1 = 0
1026 0 : offset2 = 0
1027 0 : missing_point = 1
1028 : tries = 0
1029 : retry = .true.
1030 0 : do while (retry) !If (T,Rho)_cell lies to close to the edge of the OP mono grid (in Rho), shift the 4x4 grid up/down by 1 in jne.
1031 0 : missing_point = 1
1032 0 : retry = .false.
1033 0 : do jj=1,4
1034 0 : do ii=1,4
1035 0 : dite = (ii - ii_min + offset1)*ite_step !(ii - ii_min)*ite_step + offset2
1036 0 : djne = (jj - jj_min + offset2)*jne_step !(jj - jj_min)*jne_step + offset1
1037 0 : do i =imin,imax
1038 0 : ite_i = ite(i)
1039 0 : jne_i = jne(i)
1040 :
1041 0 : if (ite_i == ite_min + dite .and. jne_i == jne_min + djne) THEN
1042 0 : logT_grid(ii,jj) = logT(i)
1043 0 : logRho_grid(ii,jj) = logRho(i)
1044 0 : i_grid(ii,jj) = i
1045 0 : missing_point(ii,jj) = 0
1046 :
1047 : end if
1048 : end do
1049 : end do
1050 : end do
1051 :
1052 0 : if (SUM(missing_point) > 0) then
1053 0 : retry = .true.
1054 0 : if (SUM(missing_point(2:4,1:3)) == 0) then
1055 0 : offset1 = offset1 + 1
1056 0 : offset2 = offset2 - 1
1057 0 : else if (SUM(missing_point(1:3,1:3)) == 0) then
1058 0 : offset1 = offset1 - 1
1059 0 : offset2 = offset2 - 1
1060 0 : else if (SUM(missing_point(2:4,2:4)) == 0) then
1061 0 : offset1 = offset1 + 1
1062 0 : offset2 = offset2 + 1
1063 0 : else if (SUM(missing_point(1:3,2:4)) == 0) then
1064 0 : offset1 = offset1 - 1
1065 0 : offset2 = offset2 + 1
1066 : else
1067 0 : if (ii_min == 3 .and. jj_min == 2) then
1068 0 : offset1 = offset1 + 2
1069 0 : offset2 = offset2 - 2
1070 0 : else if (ii_min == 2 .and. jj_min == 3) then
1071 0 : offset1 = offset1 - 2
1072 0 : offset2 = offset2 + 2
1073 0 : else if (ii_min == 2 .and. jj_min == 2) then
1074 0 : offset1 = offset1 - 2
1075 0 : offset2 = offset2 - 2
1076 0 : else if (ii_min == 3 .and. jj_min == 3) then
1077 0 : offset1 = offset1 + 2
1078 0 : offset2 = offset2 - 1
1079 : end if
1080 :
1081 : end if
1082 : end if
1083 :
1084 0 : if (tries > 2) THEN ! To prevent loop from getting stuck.
1085 0 : do_difficult_point = .false.
1086 0 : if (ite_min == 226 .and. jne_min == 90 .and. ii_min == 2 .and. jj_min == 3) do_difficult_point = .true.
1087 0 : if (ite_min == 226 .and. jne_min == 88 .and. ii_min == 2 .and. jj_min == 2) do_difficult_point = .true.
1088 0 : if (ite_min == 230 .and. jne_min == 88 .and. ii_min == 3 .and. jj_min == 2) do_difficult_point = .true.
1089 0 : if (ite_min == 230 .and. jne_min == 90 .and. ii_min == 3 .and. jj_min == 3) do_difficult_point = .true.
1090 :
1091 0 : if (do_difficult_point) then
1092 0 : do jj=1,4
1093 0 : do ii =1,4
1094 0 : dite = (ii - ii_min)*2*ite_step
1095 0 : djne = (jj - jj_min - 1)*jne_step
1096 0 : do i=imin,imax
1097 0 : ite_i = ite(i)
1098 0 : jne_i = jne(i)
1099 0 : if (ite_i == ite_min + dite .and. jne_i == jne_min + djne) then
1100 0 : logT_grid(ii,jj) = logT(i)
1101 0 : logRho_grid(ii,jj) = logRho(i)
1102 0 : i_grid(ii,jj) = i
1103 : end if
1104 : end do
1105 : end do
1106 : end do
1107 : retry = .false.
1108 : else
1109 0 : write(*,*) 'Cannot find points for interpolation compute_kappa_fast', &
1110 0 : k, ite_min, jne_min, logT_min, logRho_min,logT_cntr, logRho_cntr, &
1111 0 : missing_point, ii_min, jj_min, offset1, offset2,imin, imax
1112 0 : ierr = 1
1113 : return
1114 : !stop
1115 : end if
1116 : end if
1117 0 : tries = tries + 1
1118 : end do
1119 :
1120 0 : do jj=1,4
1121 0 : do ii=1,4
1122 0 : ik = i_grid(ii,jj)
1123 0 : lkap_Ross(ii,jj) = lkap_ross_pcg(ik)
1124 : end do
1125 : end do
1126 :
1127 0 : call rossl_interpolator% initialize()
1128 0 : do jj = 1, 4
1129 0 : do ii = 1, 4
1130 0 : call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_Ross(ii,jj))
1131 : end do
1132 : end do
1133 :
1134 0 : lkap_ross_cell = rossl_interpolator% evaluate(logT_cntr,logRho_cntr)
1135 0 : dlnkap_rad_dlnT = rossl_interpolator% evaluate_deriv(logT_cntr, logRho_cntr, .true., .false.)
1136 0 : dlnkap_rad_dlnRho = rossl_interpolator% evaluate_deriv(logT_cntr, logRho_cntr, .false., .true.)
1137 :
1138 0 : end subroutine compute_kappa_fast
1139 :
1140 :
1141 : !!! Interpolate Rossland opacity and derivatives from a 4x4 grid computed last time step.
1142 0 : subroutine interpolate_kappa(k, logT_cntr, logRho_cntr,lkap_ross_cell, dlnkap_rad_dlnT, dlnkap_rad_dlnRho, &
1143 : ierr,logT_grid, logRho_grid, lkap_grid)
1144 :
1145 0 : use cubic_interpolator, only: interpolator
1146 :
1147 : integer, intent(inout) :: ierr
1148 : integer, intent(in) :: k
1149 : real(dp), intent(in) :: logT_cntr, logRho_cntr
1150 : real(dp), intent(in) :: logT_grid(4,4), logRho_grid(4,4), lkap_grid(4,4)
1151 : real(dp), intent(out) :: lkap_ross_cell, dlnkap_rad_dlnT, dlnkap_rad_dlnRho
1152 :
1153 : integer :: ii,jj
1154 :
1155 : type(interpolator) :: rossl_interpolator
1156 :
1157 0 : ierr = 0
1158 :
1159 0 : call rossl_interpolator% initialize()
1160 0 : do jj = 1, 4
1161 0 : do ii = 1, 4
1162 0 : call rossl_interpolator% add_point(logT_grid(ii,jj), logRho_grid(ii,jj), lkap_grid(ii,jj))
1163 : end do
1164 : end do
1165 :
1166 0 : lkap_ross_cell = rossl_interpolator% evaluate(logT_cntr,logRho_cntr)
1167 0 : dlnkap_rad_dlnT = rossl_interpolator% evaluate_deriv(logT_cntr,logRho_cntr, .true., .false.)
1168 0 : dlnkap_rad_dlnRho = rossl_interpolator% evaluate_deriv(logT_cntr,logRho_cntr, .false., .true.)
1169 :
1170 0 : end subroutine interpolate_kappa
1171 : end module op_eval_mombarg
|
