Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-2019 Haili Hu & 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_osc
21 : use math_lib
22 : use op_def
23 : use const_def
24 : use kap_def, only: kap_test_partials, kap_test_partials_val, kap_test_partials_dval_dx
25 :
26 : contains
27 : ! **********************************************************************
28 0 : subroutine abund(nel, izz, fa, flmu, nkz)
29 : implicit none
30 : integer, intent(in) :: nel, izz(ipe)
31 : real, intent(in) :: fa(ipe)
32 : real, intent(out) :: flmu
33 : integer, intent(out) :: nkz(ipe)
34 : integer :: k, k1, k2, m
35 0 : real :: amamu(ipe), fmu
36 :
37 : ! Get k1,get amamu(k)
38 0 : do k = 1, nel
39 0 : do m = 1, ipe
40 0 : if (izz(k) == kz(m)) then
41 0 : amamu(k) = amass(m)
42 0 : nkz(k) = m
43 : GOTO 1
44 : end if
45 : end do
46 0 : write(*,*) ' k=',k,', izz(k)=',izz(k)
47 0 : write(*,*) ' kz(m) not found'
48 0 : stop
49 0 : 1 continue
50 : end do
51 :
52 : ! Mean atomic weight = fmu
53 0 : fmu = 0.
54 0 : do k = 1, nel
55 0 : fmu = fmu + fa(k)*amamu(k)
56 : end do
57 :
58 0 : fmu = fmu*1.660531e-24 ! Convert to cgs
59 0 : flmu = log10(dble(fmu))
60 :
61 0 : end subroutine abund
62 : ! *********************************************************************
63 0 : subroutine xindex(flt, ilab, xi, ih, i3, ierr)
64 : implicit none
65 : integer, intent(in) :: i3
66 : real, intent(in) :: flt
67 : integer, intent(out) :: ih(0:5), ilab(0:5)
68 : real, intent(out) :: xi
69 : integer, intent(out) :: ierr
70 : integer :: i, ih2
71 0 : real :: x
72 :
73 0 : ierr = 0
74 0 : if (flt < 3.5) then
75 0 : ierr = 102
76 0 : return
77 0 : else if (flt > 8.) then
78 0 : ierr = 102
79 0 : return
80 : end if
81 :
82 0 : x = 40.*flt/real(i3)
83 0 : ih2 = x
84 0 : ih2 = max(ih2, 140/i3+2)
85 0 : ih2 = min(ih2, 320/i3-3)
86 0 : do i = 0, 5
87 0 : ih(i) = ih2 + i - 2
88 0 : ilab(i) = i3*ih(i)
89 : end do
90 0 : xi = 2.*(x-ih2) - 1
91 :
92 : end subroutine xindex
93 : ! *********************************************************************
94 0 : subroutine jrange(ih, jhmin, jhmax, i3)
95 : implicit none
96 : integer, intent(in) :: ih(0:5), i3
97 : integer, intent(out) :: jhmin, jhmax
98 : integer :: i
99 :
100 0 : jhmin = 0
101 0 : jhmax = 1000
102 0 : do i = 0, 5
103 0 : jhmin = max(jhmin, js(ih(i)*i3)/i3)
104 0 : jhmax = min(jhmax, je(ih(i)*i3)/i3)
105 : end do
106 :
107 0 : end subroutine jrange
108 : ! *********************************************************************
109 0 : subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, flrho, flt, xi, flne, flmu, flr, epa, uy, i3, ierr)
110 : use op_load, only : solve
111 : implicit none
112 : integer, intent(in) :: ilab(0:5), nel, nkz(ipe), jhmin, ih(0:5), i3
113 : integer, intent(inout) :: jhmax
114 : integer, intent(out) :: ierr
115 : real, intent(in) :: fa(ipe), flt, xi, flmu
116 : real,intent(out) :: flne, uy, epa
117 : real, intent(inout) :: flrho
118 : integer :: i, j, n, jh, jm, itt, jne, jnn
119 0 : real :: flrmin, flrmax, flr(4,4), uyi(4), efa(0:5, 7:118), flrh(0:5, 7:118), u(4), flnei(4), y, zeta, efa_temp
120 : ! declare variables in common block, by default: real (a-h, o-z), integer (i-n)
121 : ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, ne1, ne2, np, kp1, kp2, kp3, npp, mx, nx
122 : ! real :: umin, umax, epatom, oplnck, fion, yy1, yy2, yx
123 : ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, &
124 : ! nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1(17,91,25), &
125 : ! ne2(17,91,25),fion(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), &
126 : ! kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), &
127 : ! yy1(33417000),yy2(120000000),nx(19305000),yx(19305000)
128 : ! save /atomdata/
129 : !
130 : ! efa(i,jh)=sum_n epa(i,jh,n)*fa(n)
131 : ! flrh(i,jh)=log10(rho(i,jh))
132 : !
133 : ! Get efa
134 0 : do i = 0, 5
135 0 : itt = (ilab(i)-ite1)/2 + 1
136 0 : do jne = jn1(itt), jn2(itt), i3
137 0 : jnn = (jne-jn1(itt))/2 + 1
138 0 : jh = jne/i3
139 0 : efa_temp = 0.
140 0 : do n = 1, nel
141 0 : efa_temp = efa_temp + epatom(nkz(n), itt, jnn)*fa(n)
142 : end do !n
143 0 : efa(i, jh) = efa_temp
144 : end do !jne
145 : end do !i
146 :
147 : ! Get range for efa > 0
148 0 : do i = 0, 5
149 0 : do jh = jhmin, jhmax
150 0 : if (efa(i, jh) <= 0.) then
151 0 : jm = jh - 1
152 : GOTO 3
153 : end if
154 : end do
155 : GOTO 4
156 0 : 3 jhmax = min(jhmax, jm)
157 0 : 4 continue
158 : end do
159 :
160 : ! Get flrh
161 0 : do jh = jhmin,jhmax
162 0 : do i = 0,5
163 0 : flrh(i, jh) = flmu + 0.25*i3*jh - log10(dble(efa(i,jh)))
164 : end do
165 : end do
166 :
167 : ! Find flrmin and flrmax
168 0 : flrmin = -1000
169 0 : flrmax = 1000
170 0 : do i = 0, 5
171 0 : flrmin = max(flrmin, flrh(i,jhmin))
172 0 : flrmax = min(flrmax, flrh(i,jhmax))
173 : end do
174 :
175 : ! Check range of flrho
176 0 : if (flrho < flrmin .or. flrho > flrmax) then
177 0 : ierr = 101
178 0 : return
179 : end if
180 :
181 : ! Interpolations in j for flne
182 0 : do jh = jhmin, jhmax
183 0 : if (flrh(2,jh) > flrho) then
184 0 : jm = jh - 1
185 : GOTO 5
186 : end if
187 : end do
188 0 : write(*,*) ' Interpolations in j for flne'
189 0 : write(*,*) ' Not found, i=',i
190 0 : stop
191 0 : 5 jm=max(jm,jhmin+1)
192 0 : jm=min(jm,jhmax-2)
193 :
194 0 : do i = 1, 4
195 0 : do j = 1, 4
196 0 : u(j) = flrh(i, jm+j-2)
197 0 : flr(i,j) = flrh(i, jm+j-2)
198 : end do
199 0 : call solve(u, flrho, zeta, uyi(i), ierr)
200 0 : if (ierr /= 0) return
201 0 : y = jm + 0.5*(zeta+1)
202 0 : flnei(i) = .25*i3*y
203 : end do
204 :
205 : ! Interpolations in i
206 0 : flne = fint(flnei, xi)
207 0 : uy = fint(uyi, xi)
208 : ! Get epa
209 0 : epa = exp10(dble(flne + flmu - flrho))
210 :
211 0 : return
212 :
213 : ! 601 format(' For flt=',1p,e11.3,', flrho=',e11.3,' is out of range. Allowed range for flrho is ',e11.3,' to ',e11.3)
214 0 : end subroutine findne
215 :
216 : ! **********************************************************************
217 0 : subroutine yindex(jhmin, jhmax, flne, jh, i3, eta)
218 : implicit none
219 : integer, intent(in) :: jhmin, jhmax, i3
220 : real, intent(in) :: flne
221 : integer, intent(out) :: jh(0:5)
222 : real, intent(out) :: eta
223 : integer :: j, k
224 0 : real :: y
225 :
226 0 : y = 4.*flne/real(i3)
227 0 : j = y
228 0 : j = max(j,jhmin+2)
229 0 : j = min(j,jhmax-3)
230 0 : do k = 0, 5
231 0 : jh(k) = j + k - 2
232 : end do
233 0 : eta = 2.*(y-j)-1
234 :
235 0 : end subroutine yindex
236 : ! **********************************************************************
237 0 : subroutine findux(flr, xi, eta, ux)
238 : implicit none
239 : real, intent(in) :: flr(4, 4), xi, eta
240 : real, intent(out) :: ux
241 : integer :: i, j
242 0 : real :: uxj(4), u(4)
243 :
244 0 : do j = 1, 4
245 0 : do i = 1, 4
246 0 : u(i) = flr(i, j)
247 : end do
248 0 : uxj(j) = fintp(u, xi)
249 : end do
250 0 : ux = fint(uxj, eta)
251 :
252 0 : end subroutine findux
253 : ! *********************************************************************
254 0 : subroutine rd(nel, nkz, izz, ilab, jh, n_tot, ff, rr, i3, umesh, fac)
255 : implicit none
256 : integer, intent(in) :: nel, nkz(ipe), izz(ipe), ilab(0:5), jh(0:5), n_tot, i3
257 : real, intent(in) :: umesh(nptot)
258 : real(dp), intent(in) :: fac(nel)
259 : real, intent(out) :: ff(:,:,0:,0:) ! (nptot, ipe, 6, 6)
260 : real, intent(out) :: rr(28, ipe, 0:5, 0:5)
261 : integer :: i, j, k, l, m, n, itt, jnn, izp, ne1, ne2, ne, ib, ia
262 0 : real :: fion(-1:28), yb, ya, d
263 : ! declare variables in common block (instead of by default: real (a-h, o-z), integer (i-n))
264 : ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, ne1p, ne2p, np, kp1, kp2, kp3, npp, mx, nx
265 : ! real :: umin, umax, epatom, oplnck, fionp, yy1, yy2, yx
266 : ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, &
267 : ! nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1p(17,91,25), &
268 : ! ne2p(17,91,25),fionp(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), &
269 : ! kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), &
270 : ! yy1(33417000),yy2(120000000),nx(19305000),yx(19305000)
271 : ! save /atomdata/
272 : !
273 : ! i=temperature index
274 : ! j=density index
275 : ! k=frequency index
276 : ! n=element index
277 : ! Get:
278 : ! mono opacity cross-section ff(k,n,i,j)
279 : ! modified cross-section for selected element, ta(k,i,j)
280 : !
281 : ! Initialisations
282 0 : rr=0.
283 0 : ff=0.
284 :
285 : ! Start loop on i (temperature index)
286 0 : do i = 0, 5
287 0 : itt = (ilab(i) - ite1)/2 + 1
288 0 : do j = 0, 5
289 0 : jnn = (jh(j)*i3 - jn1(itt))/2 + 1
290 : ! Read mono opacities
291 0 : do n = 1, nel
292 0 : izp = izz(n)
293 0 : ne1 = ne1p(nkz(n), itt, jnn)
294 0 : ne2 = ne2p(nkz(n), itt, jnn)
295 0 : do ne = ne1, ne2
296 0 : fion(ne) = fionp(ne, nkz(n), itt, jnn)
297 : end do
298 0 : do ne = ne1, min(ne2, izp-2)
299 0 : rr(izp-1-ne, n, i, j) = fion(ne)
300 : end do
301 :
302 0 : do k = 1, n_tot
303 0 : ff(k, n, i, j) = yy2(k+kp2(nkz(n), itt, jnn))
304 : end do
305 :
306 0 : if (fac(n) /= 1d0) then
307 0 : do k = 1, size(ff,dim=1)
308 0 : ff(k,n,i,j) = fac(n)*ff(k,n,i,j)
309 : end do
310 : end if
311 : end do !n
312 : end do !j
313 : end do !i
314 :
315 0 : return
316 :
317 : end subroutine rd
318 : ! **********************************************************************
319 0 : subroutine ross(flmu, dv, ntot,rs, rossl)
320 : implicit none
321 : integer, intent(in) :: ntot
322 : real, intent(in) :: flmu, dv, rs(nptot, 0:5, 0:5)
323 : real, intent(out) :: rossl(0:5, 0:5)
324 : integer :: i, j, n
325 0 : real(dp) :: drs, dd, oross
326 : real :: fmu, tt
327 :
328 : ! oross=cross-section in a.u.
329 : ! rossl=log10(ROSS in cgs)
330 0 : do i = 0, 5
331 0 : do j = 0, 5
332 0 : drs = 0.d0
333 0 : do n = 1, ntot
334 0 : dd = 1.d0/rs(n, i, j)
335 0 : drs = drs + dd
336 : end do
337 0 : oross = 1.d0/(drs*dv)
338 0 : rossl(i, j) = log10(oross) - 16.55280d0 - flmu !log10(fmu)
339 : end do !j
340 : end do !i
341 :
342 0 : end subroutine ross
343 : ! **********************************************************************
344 0 : subroutine mix(ntot, nel, fa, ff, rs, rr, rion)
345 : implicit none
346 : integer, intent(in) :: ntot, nel
347 : real, intent(in) :: ff(nptot, ipe, 0:5, 0:5), fa(ipe), rr(28, 17, 0:5, 0:5)
348 : real, intent(out) :: rs(nptot, 0:5, 0:5), rion(28, 0:5, 0:5)
349 : integer :: i, j, k, n, m
350 : real :: rs_temp, rion_temp
351 :
352 0 : do i = 0, 5
353 0 : do j = 0, 5
354 0 : do n = 1, ntot
355 : !rs_temp = ff(n,1,i,j)*fa(1)
356 : !do k = 2, nel
357 : ! rs_temp = rs_temp + ff(n,k,i,j)*fa(k)
358 : !end do
359 : !rs(n,i,j) = rs_temp
360 0 : rs(n, i, j) = dot_product(ff(n,1:nel,i,j),fa(1:nel))
361 : end do
362 0 : do m = 1, 28
363 : !rion_temp = rr(m, 1, i, j)*fa(1)
364 : !do k = 2, nel
365 : ! rion_temp = rion_temp + rr(m,k,i,j)*fa(k)
366 : !end do
367 : !rion(m,i,j) = rion_temp
368 0 : rion(m,i,j) = dot_product(rr(m,1:nel,i,j),fa(1:nel))
369 : end do
370 : end do
371 : end do
372 :
373 0 : end subroutine mix
374 : ! **********************************************************************
375 0 : subroutine interp(nel, rossl, xi, eta, g, i3, ux, uy, gx, gy)
376 : implicit none
377 : integer, intent(in) :: nel, i3
378 : real, intent(in) :: ux, uy, xi, eta
379 : real, intent(out) :: gx, gy, g
380 : integer :: i, j, l
381 : real :: V(4), U(4), vyi(4)
382 0 : real :: x3(3), fx3!, fxxy(0:5, 0:5), fyxy(0:5, 0:5)
383 : ! pointers and targets
384 0 : real, target :: fx(0:5, 0:5), fy(0:5, 0:5), fxy(0:5, 0:5), fyx(0:5, 0:5), fxx(0:5, 0:5), fyy(0:5, 0:5), rossl(0:5, 0:5)
385 0 : real, pointer :: f3(:), fin(:, :), finx(:, :), finy(:, :)
386 : !
387 : ! interpolation of g (=rosseland mean opacity)
388 : ! Use refined techniques of bi-cubic spline interpolation (Seaton 1993):
389 : ! call deriv(rossl, fx, fy, fxy)
390 : ! call interp2(rossl, fx, fy, fxy, xi, eta, g, gx, gy)
391 : ! gy = 0.5*gy/uy
392 : ! gx = (80./real(i3))*(0.5*gx-gy*ux)
393 : ! Alternatively, use interpolation techniques by M.-A. Dupret to ensure smooothness required
394 : ! for pulsation studies:
395 0 : do i = 0, 3
396 0 : x3(1) = i
397 0 : x3(2) = i+1
398 0 : x3(3) = i+2
399 0 : do j = 1, 4
400 0 : f3 => rossl(i:i+2, j)
401 0 : call deriv3(f3, x3, fx3)
402 0 : fx(i+1, j) = fx3
403 : end do
404 : end do
405 0 : do i = 0, 3
406 0 : x3(1) = i
407 0 : x3(2) = i+1
408 0 : x3(3) = i+2
409 0 : do j = 1, 4
410 0 : f3 => rossl(j, i:i+2)
411 0 : call deriv3(f3, x3, fx3)
412 0 : fy(j, i+1) = fx3
413 : end do
414 : end do
415 0 : do i = 1, 2
416 0 : x3(1) = i
417 0 : x3(2) = i+1
418 0 : x3(3) = i+2
419 0 : do j = 2, 3
420 0 : f3 => fx(i:i+2, j)
421 0 : call deriv3(f3, x3, fx3)
422 0 : fxx(i+1, j) = fx3
423 : end do
424 : end do
425 0 : do i = 1, 2
426 0 : x3(1) = i
427 0 : x3(2) = i+1
428 0 : x3(3) = i+2
429 0 : do j = 2, 3
430 0 : f3 => fx(j, i:i+2)
431 0 : call deriv3(f3, x3, fx3)
432 0 : fxy(j, i+1) = fx3
433 : end do
434 : end do
435 0 : do i = 1, 2
436 0 : x3(1) = i
437 0 : x3(2) = i+1
438 0 : x3(3) = i+2
439 0 : do j = 2, 3
440 0 : f3 => fy(i:i+2, j)
441 0 : call deriv3(f3, x3, fx3)
442 0 : fyx(i+1, j) = fx3
443 : end do
444 : end do
445 0 : do i = 1, 2
446 0 : x3(1) = i
447 0 : x3(2) = i+1
448 0 : x3(3) = i+2
449 0 : do j = 2, 3
450 0 : f3 => fy(j, i:i+2)
451 0 : call deriv3(f3, x3, fx3)
452 0 : fyy(j, i+1) = fx3
453 : end do
454 : end do
455 : ! call deriv(rossl, fx, fy, fxy)
456 : ! call deriv(fx, fxx, fxy, fxxy)
457 : ! call deriv(fy, fyx, fyy, fyxy)
458 0 : fin => rossl(2:3, 2:3)
459 0 : finx => fx(2:3, 2:3)
460 0 : finy => fy(2:3, 2:3)
461 0 : call interp3(fin, finx, finy, xi, eta, g)
462 0 : fin => fx(2:3, 2:3)
463 0 : finx => fxx(2:3, 2:3)
464 0 : finy => fxy(2:3, 2:3)
465 0 : call interp3(fin, finx, finy, xi, eta, gx)
466 0 : fin => fy(2:3, 2:3)
467 0 : finx => fyx(2:3, 2:3)
468 0 : finy => fyy(2:3, 2:3)
469 0 : call interp3(fin, finx, finy, xi, eta, gy)
470 0 : gy = 0.5*gy/uy
471 0 : gx = (80./real(i3))*(0.5*gx-gy*ux)
472 :
473 0 : end subroutine interp
474 : ! *************************************
475 0 : function fint(u,r)
476 : dimension u(4)
477 :
478 : ! If P(R) = u(1) u(2) u(3) u(4)
479 : ! for R = -3 -1 1 3
480 : ! then a cubic fit is:
481 : P(R)=(
482 : & 27*(u(3)+u(2))-3*(u(1)+u(4)) +R*(
483 : & 27*(u(3)-u(2))-(u(4)-u(1)) +R*(
484 : & -3*(u(2)+u(3))+3*(u(4)+u(1)) +R*(
485 : & -3*(u(3)-u(2))+(u(4)-u(1)) ))))/48.
486 :
487 0 : fint=p(r)
488 :
489 0 : end function fint
490 : ! **********************************************************************
491 0 : function fintp(u,r)
492 : dimension u(4)
493 :
494 : ! If P(R) = u(1) u(2) u(3) u(4)
495 : ! for R = -3 -1 1 3
496 : ! then a cubic fit to the derivative is:
497 : PP(R)=(
498 : & 27*(u(3)-u(2))-(u(4)-u(1)) +2.*R*(
499 : & -3*(u(2)+u(3))+3*(u(4)+u(1)) +3.*R*(
500 : & -3*(u(3)-u(2))+(u(4)-u(1)) )))/48.
501 :
502 0 : fintp=pp(r)
503 :
504 0 : end function fintp
505 :
506 : ! **********************************************************************
507 0 : subroutine DERIV(f, fx, fy, fxy)
508 :
509 : real, intent(in) :: f(0:5, 0:5)
510 : real, intent(out) :: fx(0:5, 0:5), fy(0:5, 0:5), fxy(0:5, 0:5)
511 0 : real :: C(6)
512 :
513 : ! GET FX
514 0 : do J = 0, 5
515 0 : L=0
516 0 : do I = 0, 5
517 0 : L=L+1
518 0 : C(L)=F(I,J)
519 : end do
520 0 : call GET(C,L)
521 0 : L=0
522 0 : do I = 0, 5
523 0 : L = L + 1
524 0 : FX(I, J) = C(L)
525 : end do
526 : end do
527 :
528 : ! GET FY
529 0 : do I = 0, 5
530 0 : L=0
531 0 : do J = 0, 5
532 0 : L = L + 1
533 0 : C(L) = F(I, J)
534 : end do
535 0 : call GET(C,L)
536 0 : L=0
537 0 : do J = 0, 5
538 0 : L = L + 1
539 0 : FY(I,J) = C(L)
540 : end do
541 : end do
542 :
543 : ! GET FXY
544 0 : do I = 0, 5
545 0 : L = 0
546 0 : do J = 0, 5
547 0 : L = L + 1
548 0 : C(L) = FX(I, J)
549 : end do
550 0 : call GET(C,L)
551 0 : L=0
552 0 : do J = 0, 5
553 0 : L = L + 1
554 0 : FXY(I,J) = C(L)
555 : end do
556 : end do
557 :
558 0 : end subroutine DERIV
559 : ! *****************************************************************
560 : !
561 0 : subroutine GET(F,N)
562 : !
563 : ! SIMPLIFIED CODE FOR SPLINE COEFFICIENTS, FOR CASE OF INTERVALS
564 : ! OF UNITY.
565 : ! RETURNS DERIVATIVES OF ORIGINAL F IN LOCATION F
566 : !
567 : ! REVISED 5.5.95
568 : !
569 : parameter (IPI=6)
570 0 : DIMENSION F(IPI),D(IPI),T(IPI)
571 :
572 0 : if (N <= 0) then
573 0 : WRITE(6,*)' Error in subroutine GET: N=',N
574 0 : STOP
575 0 : else if (N == 1) then
576 0 : F(1)=0.
577 0 : return
578 0 : else if (N == 2) then
579 0 : F(1)=F(2)-F(1)
580 0 : F(2)=F(1)
581 0 : return
582 0 : else if (N == 3) then
583 0 : FP1=.5*(-3.*F(1)+4.*F(2)-F(3))
584 0 : FPN=.5*(F(1)-4.*F(2)+3.*F(3))
585 : else
586 0 : FP1=(-11.*F(1)+18.*F(2)-9.*F(3)+2.*F(4))/6.
587 0 : FPN=(11.*F(N)-18.*F(N-1)+9.*F(N-2)-2.*F(N-3))/6.
588 : end if
589 :
590 0 : D(1)=-.5
591 0 : T(1)=.5*(-F(1)+F(2)-FP1)
592 :
593 0 : do J=2,N-1
594 0 : D(J)=-1./(4.+D(J-1))
595 0 : T(J)=-D(J)*(F(J-1)-2.*F(J)+F(J+1)-T(J-1))
596 : end do
597 :
598 0 : D(N)=(FPN+F(N-1)-F(N)-T(N-1))/(2.+D(N-1))
599 :
600 0 : do J=N-1,1,-1
601 0 : D(J)=D(J)*D(J+1)+T(J)
602 : end do
603 :
604 0 : do J=2,N-1
605 0 : F(J)=-F(J)+F(J+1)-2.*D(J)-D(J+1)
606 : end do
607 0 : F(1)=FP1
608 0 : F(N)=FPN
609 :
610 : end subroutine GET
611 :
612 : ! **********************************************************************
613 0 : subroutine INTERP2(f, fx, fy, fxy, xi, eta, g, gx, gy)
614 : real, intent(in) :: eta, xi, f(0:5, 0:5), fx(0:5, 0:5), fy(0:5, 0:5), fxy(0:5, 0:5)
615 : real, intent(out) :: g, gx, gy
616 : integer :: i , j
617 0 : real :: x, y, B(16)
618 : !
619 : ! function DEFINITIONS FOR CUBIC EXPANSION
620 : !
621 0 : FF(S,T)= B( 1)+T*(B( 2)+T*(B( 3)+T*B( 4)))
622 : & +S*( B( 5)+T*(B( 6)+T*(B( 7)+T*B( 8)))
623 : & +S*( B( 9)+T*(B(10)+T*(B(11)+T*B(12)))
624 : & +S*( B(13)+T*(B(14)+T*(B(15)+T*B(16))) )))
625 :
626 : FFX(S,T)= B( 5)+T*(B( 6)+T*(B( 7)+T*B( 8)))
627 : & +S*( 2*(B( 9)+T*(B(10)+T*(B(11)+T*B(12))))
628 : & +S*( 3*(B(13)+T*(B(14)+T*(B(15)+T*B(16)))) ))
629 :
630 : FFY(S,T)= B( 2)+S*(B( 6)+S*(B(10)+S*B(14)))
631 : & +T*( 2*(B( 3)+S*(B( 7)+S*(B(11)+S*B(15))))
632 : & +T*( 3*(B( 4)+S*(B( 8)+S*(B(12)+S*B(16)))) ))
633 :
634 0 : Y = (eta + 5.)/2.
635 0 : X = (xi + 5.)/2.
636 : ! i = floor(x)
637 : ! j = floor(y)
638 0 : I = X + 1.E-5
639 0 : if (ABS(X-I) <= 1.E-5) X = I
640 0 : J = Y + 1.E-5
641 0 : if (ABS(Y-J) <= 1.E-5) Y = J
642 : !
643 : ! INTERPOLATE
644 : !
645 : ! GIVEN FUNCTIONS AND DERIVATIVES AT GRID POINTS, COMPUTE COEFFICIENTS.
646 0 : B(1)=F(I,J)
647 0 : B(2)=FY(I,J)
648 0 : B(3)=3*(-F(I,J)+F(I,J+1))-2*FY(I,J)-FY(I,J+1)
649 0 : B(4)=2*(F(I,J)-F(I,J+1))+FY(I,J)+FY(I,J+1)
650 :
651 0 : B(5)=FX(I,J)
652 0 : B(6)=FXY(I,J)
653 0 : B(7)=3*(-FX(I,J)+FX(I,J+1))-2*FXY(I,J)-FXY(I,J+1)
654 0 : B(8)=2*(FX(I,J)-FX(I,J+1))+FXY(I,J)+FXY(I,J+1)
655 :
656 0 : B(9)=3*(-F(I,J)+F(I+1,J))-2*FX(I,J)-FX(I+1,J)
657 0 : B(10)=3*(-FY(I,J)+FY(I+1,J))-2*FXY(I,J)-FXY(I+1,J)
658 : B(11)=9*(F(I,J)-F(I+1,J)+F(I+1,J+1)-F(I,J+1))
659 : & +6*(FX(I,J)-FX(I,J+1)+FY(I,J)-FY(I+1,J))
660 : & +4*FXY(I,J)
661 : & +3*(FX(I+1,J)-FX(I+1,J+1)-FY(I+1,J+1)+FY(I,J+1))
662 : & +2*(FXY(I,J+1)+FXY(I+1,J))
663 0 : & +FXY(I+1,J+1)
664 : B(12)=6*(-F(I,J)+F(I+1,J)-F(I+1,J+1)+F(I,J+1))
665 : & +4*(-FX(I,J)+FX(I,J+1))
666 : & +3*(-FY(I,J)+FY(I+1,J)+FY(I+1,J+1)-FY(I,J+1))
667 : & +2*(-FX(I+1,J)+FX(I+1,J+1)-FXY(I,J)-FXY(I,J+1))
668 0 : & -FXY(I+1,J)-FXY(I+1,J+1)
669 :
670 0 : B(13)=2*(F(I,J)-F(I+1,J))+FX(I,J)+FX(I+1,J)
671 0 : B(14)=2*(FY(I,J)-FY(I+1,J))+FXY(I,J)+FXY(I+1,J)
672 : B(15)=6*(-F(I,J)+F(I+1,J)-F(I+1,J+1)+F(I,J+1))
673 : & +4*(-FY(I,J)+FY(I+1,J))
674 : & +3*(-FX(I,J)-FX(I+1,J)+FX(I+1,J+1)+FX(I,J+1))
675 : & +2*(FY(I+1,J+1)-FY(I,J+1)-FXY(I,J)-FXY(I+1,J))
676 0 : & -FXY(I+1,J+1)-FXY(I,J+1)
677 : B(16)=4*(F(I,J)-F(I+1,J)+F(I+1,J+1)-F(I,J+1))
678 : & +2*(FX(I,J)+FX(I+1,J)-FX(I+1,J+1)-FX(I,J+1)
679 : & +FY(I,J)-FY(I+1,J)-FY(I+1,J+1)+FY(I,J+1))
680 0 : & +FXY(I,J)+FXY(I+1,J)+FXY(I+1,J+1)+FXY(I,J+1)
681 : !
682 : ! GET G=LOG10(ROSS), DGDT=d LOG10(ROSS)/d LOG10(T),
683 : ! DGDRHO=d LOG10(ROSS)/d LOG10(RHO)
684 0 : U = X - I
685 0 : V = Y - J
686 0 : G = FF(U, V)
687 0 : gx = FFX(U, V)
688 0 : gy = FFY(U, V)
689 : ! DGDT=(1./CT)*FFX(U,V)-(3./CN)*FFY(U,V)
690 : ! DGDRHO=(1./CN)*FFY(U,V)
691 :
692 0 : end subroutine INTERP2
693 : ! ************************************************************
694 : ! This subroutine estimates the partial derivative of a function
695 : ! as described in the PhD thesis by M.-A. Dupret.
696 0 : subroutine deriv3(f, x, fx)
697 : real, intent(in) :: f(3), x(3)
698 : real, intent(out) :: fx
699 0 : real :: a, a1, a2, b, x1, x2
700 :
701 0 : x1 = (f(3)- f(2))/(x(3)-x(2))
702 0 : x2 = (f(2)-f(1))/(x(2)-x(1))
703 0 : b = abs(x1/x2)
704 0 : a = (2.*x(2) - x(1) - x(3))/((x(2)-x(1))*(x(2)-x(3)))
705 0 : a1 = (x(2)-x(3))/((x(1)-x(3))*(x(1)-x(2)))
706 0 : a2 = (x(2)-x(1))/((x(3)-x(2))*(x(3)-x(1)))
707 0 : if (b >= 0.2 .and. b <= 5.) then
708 0 : fx = a*f(2) + a1*f(1) + a2*f(3)
709 : else
710 0 : if (abs(x1)<abs(x2)) fx = sign(1.0,x1)*min(abs(x1), abs(x2))
711 0 : if (abs(x2)<abs(x1)) fx = sign(1.0,x2)*min(abs(x1), abs(x2))
712 : end if
713 :
714 0 : end subroutine deriv3
715 : ! **********************************************************************
716 0 : subroutine interp3(f, fx, fy, xi, eta, g)
717 : real, intent(in) :: eta, xi, f(2:3, 2:3), fx(2:3, 2:3), fy(2:3, 2:3)
718 : real, intent(out) :: g
719 : real :: x, y
720 :
721 0 : x = (xi + 5.)/2.
722 0 : y = (eta + 5.)/2.
723 :
724 0 : P11 = P1(x,2.,3.)*P1(y,2.,3.)
725 0 : P21 = P2(x,2.,3.)*P1(y,2.,3.)
726 0 : P12 = P1(x,2.,3.)*P2(y,2.,3.)
727 0 : P22 = P2(x,2.,3.)*P2(y,2.,3.)
728 0 : Px11 = Px1(x,2.,3.)*P1(y,2.,3.)
729 0 : Px21 = Px2(x,2.,3.)*P1(y,2.,3.)
730 0 : Px12 = Px1(x,2.,3.)*P2(y,2.,3.)
731 0 : Px22 = Px2(x,2.,3.)*P2(y,2.,3.)
732 0 : Py11 = P1(x,2.,3.)*Px1(y,2.,3.)
733 0 : Py21 = P2(x,2.,3.)*Px1(y,2.,3.)
734 0 : Py12 = P1(x,2.,3.)*Px2(y,2.,3.)
735 0 : Py22 = P2(x,2.,3.)*Px2(y,2.,3.)
736 :
737 : g = f(2,2)*P11 + f(3,2)*P21 + f(2,3)*P12 + f(3,3)*P22
738 : & + fx(2,2)*Px11 + fx(3,2)*Px21 + fx(2,3)*Px12 + fx(3,3)*Px22
739 0 : & + fy(2,2)*Py11 + fy(3,2)*Py21 + fy(2,3)*Py12 + fy(3,3)*Py22
740 :
741 0 : end subroutine interp3
742 : ! **********************************************************************
743 0 : function P1(x, xi, xj)
744 :
745 0 : P1 = (x - xj)*(x - xj) * (2*x - 3*xi + xj)/((xj - xi)*(xj - xi)*(xj - xi))
746 :
747 0 : end function P1
748 : ! **********************************************************************
749 0 : function P2(x, xi, xj)
750 : real, intent(in) :: x, xi, xj
751 :
752 0 : P2 = -(x - xi)*(x - xi) * (2*x - 3*xj + xi)/((xj - xi)*(xj - xi)*(xj - xi))
753 :
754 0 : end function P2
755 : ! **********************************************************************
756 0 : function Px1(x, xi, xj)
757 :
758 0 : Px1 = (x - xi) * (x - xj)*(x - xj)/((xj - xi)*(xj - xi))
759 :
760 0 : end function Px1
761 : ! **********************************************************************
762 0 : function Px2(x, xi, xj)
763 :
764 0 : Px2 = (x - xi)*(x - xi) * (x - xj)/((xj - xi)*(xj - xi))
765 :
766 0 : end function Px2
767 : ! **********************************************************************
768 0 : subroutine scatt(ih, jh, rion, uf, f, umesh, semesh, dscat, ntot, epa, ierr)
769 : use op_load, only: BRCKR
770 : integer, intent(inout) :: ierr
771 : dimension rion(28, 0:5, 0:5),uf(0:100),f(nptot,0:5,0:5),umesh(nptot),
772 0 : & semesh(nptot),fscat(0:100),p(nptot),rr(28),ih(0:5),jh(0:5)
773 : integer :: i,j,k,n
774 : ! HH: always use meshtype q='m'
775 0 : ite3=2
776 0 : umin=umesh(1)
777 0 : CSCAT=EPA*2.37567E-8
778 0 : ierr = 0
779 :
780 0 : do i = 0, 5
781 0 : ft = exp10(ITE3*dble(ih(i))/40d0)
782 0 : do j = 0, 5
783 0 : fne = exp10(ITE3*dble(jh(j))/4d0)
784 0 : do k = 1, ntot
785 0 : p(k) = f(k, i, j)
786 : end do
787 0 : do m = 1, 28
788 0 : rr(m) = rion(m, i, j)
789 : end do
790 0 : call BRCKR(FT,FNE,RR,28,UF,100,FSCAT,ierr)
791 0 : if (ierr /= 0) return
792 0 : do n = 0, 100
793 0 : u = uf(n)
794 0 : fscat(n) = cscat*(fscat(n)-1)
795 : end do
796 0 : do n = 2, ntot-1
797 0 : u = umesh(n)
798 0 : se = semesh(n)
799 0 : m=(u-umin)/dscat
800 0 : ua=umin+dscat*m
801 0 : ub=ua+dscat
802 0 : p(n)=p(n)+((ub-u)*fscat(m)+(u-ua)*fscat(m+1))/(dscat*se)
803 : end do
804 0 : u=umesh(ntot)
805 0 : se=semesh(ntot)
806 0 : p(ntot)=p(ntot)+fscat(100)/se
807 0 : p(1)=p(1)+fscat(1)/(1.-.5*umin)
808 0 : do k=1,ntot
809 0 : f(k,i,j)=p(k)
810 : end do
811 : end do
812 : end do
813 :
814 0 : end subroutine scatt
815 : ! **********************************************************************
816 0 : subroutine screen1(ih,jh,rion,umesh,ntot,epa,f)
817 0 : use op_load, only: screen2
818 : dimension uf(0:100), umesh(nptot), fscat(0:100), ih(0:5), jh(0:5)
819 : real, target :: f(nptot, 0:5, 0:5), rion(28, 0:5, 0:5)
820 : integer :: i, j, k, m
821 0 : real, pointer :: p(:), rr(:)
822 :
823 0 : ite3=2
824 0 : umin=umesh(1)
825 0 : umax=umesh(ntot)
826 :
827 0 : do i = 0, 5
828 0 : ft=exp10(ITE3*dble(ih(i))/40d0)
829 0 : do j = 0, 5
830 0 : fne=exp10(ITE3*dble(jh(j))/4d0)
831 : ! do k=1,ntot
832 0 : p => f(1:ntot,i,j)
833 : ! end do
834 : ! do m=1,28
835 0 : rr => rion(1:28,i,j)
836 : ! end do
837 0 : call screen2(ft,fne,rr,epa,ntot,umin,umax,umesh,p)
838 : ! do k=1,ntot
839 : ! f(k,i,j)=p(k)
840 : ! end do
841 : end do
842 : end do
843 :
844 0 : end subroutine screen1
845 :
846 : end module op_osc
|
