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