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 : ! FORTRAN 90 module for calculation of radiative accelerations,
21 : ! based on the Opacity Project (OP) code "OPserver".
22 : ! See CHANGES_HU for changes made to the original code.
23 : !
24 : ! Haili Hu 2010
25 :
26 : module op_load
27 : use math_lib
28 : use op_def
29 : logical :: have_loaded_op = .false.
30 :
31 : contains
32 : ! *****************************************************************
33 0 : subroutine op_dload(path, cache_filename, ierr)
34 : implicit none
35 : character (len=*), intent(in) :: path, cache_filename
36 : integer, intent(out) :: ierr
37 :
38 0 : real :: am,amm,delp,dpack
39 : integer :: ios,it,ite11,ite22,ite33,itt,itte1,itte2,itte3,izz,jne,ite
40 : integer :: jne1,jne22,jne33,jnn,k,n,ncount2,ncount3,ja,jn
41 : integer :: jne11,jne2,nccc,ne,nfff,ntott,nn
42 0 : real :: orss,um,ux,umaxx,uminn,u
43 : real,dimension(nptot):: umesh, semesh
44 : integer :: ntotv
45 : real :: dv,dv1
46 :
47 : integer :: cache_version
48 :
49 : common /mesh/ ntotv,dv,dv1,umesh,semesh
50 : ! common /atomdata/
51 : ! common/atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, &
52 : ! nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1(17,91,25), &
53 : ! ne2(17,91,25),fion(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), &
54 : ! kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), &
55 : ! yy1(33417000),yy2(120000000),nx(19305000),yx(19305000)
56 :
57 : integer,dimension(ipe) :: ifl,iflp
58 : character :: num(0:9)*1,zlab(ipe)*3,tlab*6,zlabp(ipe)*3
59 : DATA NUM/'0','1','2','3','4','5','6','7','8','9'/
60 :
61 : integer :: kz(17)
62 : data kz/1, 2, 6, 7, 8, 10, 11, 12, 13, 14, 16, 18, 20, 24, 25, 26, 28/
63 : save /mesh/ !HH: put common block in static memory
64 :
65 : integer :: nx_temp(nptot)
66 0 : real :: y_temp(nptot)
67 : integer :: nx_index, left_n, right_n, n_index
68 0 : real :: left_val, right_val, cross_section, slope
69 :
70 0 : if (allocated(yy2) .eqv. .false.) then
71 : ! yy2 actually needs 29,563 x 10,000 length
72 0 : ALLOCATE(yy2(30000*10000),nx(19305000),yx(19305000),stat=ierr)
73 0 : if (ierr/=0) return
74 0 : yy2=0.0
75 0 : nx=0.0
76 0 : yx=0.0
77 : ! write(*,*) "ierr",ierr
78 : end if
79 :
80 :
81 0 : ierr=0
82 0 : if (have_loaded_op) return
83 :
84 0 : !$omp critical (critical_do_op_dload)
85 :
86 0 : if (have_loaded_op) GOTO 1001
87 :
88 : !path = '../OP4STARS_1.3'
89 : !call getenv("oppath", path)
90 : !if (len(trim(path)) == 0) then
91 : ! write(6,*) 'Define environmental variable oppath (directory of OP data)'
92 : ! stop
93 : !end if
94 :
95 0 : ios = 0
96 0 : open(1,file=trim(cache_filename),action='read',status='old',iostat=ios,form='unformatted')
97 0 : if (ios == 0) then
98 0 : write(*,*) 'reading OP cache file ' // trim(cache_filename)
99 0 : read(1,iostat=ios) cache_version,ntotv,dv,dv1,umesh,ite1,ite2,ite3,jn1,jn2,jne3,
100 0 : & umin,umax,ntotp,nc,nf,int,epatom,oplnck,ne1p,ne2p,fionp,np,kp1,kp2,kp3,npp,yy2,nx,yx
101 0 : write(*,*) 'done reading OP cache file'
102 0 : close(1)
103 0 : if (cache_version /= op_cache_version) then
104 0 : write(*,*) 'wrong version of OP cache'
105 0 : write(*,*) 'cache file path is set by op_mono_data_cache_filename'
106 0 : write(*,*) 'perhaps cache is shared between different MESA versions'
107 0 : write(*,*) 'please delete cache file and try again'
108 0 : stop 'op_load'
109 : end if
110 0 : if (ios == 0) then
111 0 : have_loaded_op = .true.
112 0 : call IMESH(UMESH,NTOTV)
113 0 : GOTO 1001
114 : end if
115 0 : write(*,*) 'failed in reading cache file'
116 0 : write(*,*) 'please delete cache file and try again'
117 0 : stop 'op_load'
118 : end if
119 :
120 0 : ite1=140
121 0 : ite2=320
122 0 : ite3=2
123 0 : do n=1,ipe
124 0 : ifl(n)=50+n
125 0 : zlab(n)='m'//num(kz(n)/10)//num(kz(n)-10*(kz(n)/10))
126 0 : iflp(n)=70+n
127 0 : zlabp(n)='a'//num(kz(n)/10)//num(kz(n)-10*(kz(n)/10))
128 : end do
129 :
130 0 : write(*,*) 'loading OP mono data...'
131 :
132 : ! READ INDEX FILES
133 : ! FIRST FILE
134 0 : NN=1
135 : ! write(*,*) ' Opening '//'./'//zlab(1)//'.index'
136 0 : OPEN(1,FILE=trim(path)//'/'//ZLAB(1)//'.index',STATUS='OLD',iostat=ios)
137 0 : if (ios /= 0) then
138 0 : write(*,*) 'failed to open ' // trim(path) // '/' // ZLAB(1) // '.index'
139 0 : ierr = -1
140 0 : GOTO 1001
141 : end if
142 0 : read(1,*)IZZ,AMM
143 0 : read(1,*)ITTE1,ITTE2,ITTE3
144 0 : read(1,*)UMIN,UMAX
145 0 : read(1,*)NC,NF
146 0 : read(1,*)DPACK
147 0 : CLOSE(1)
148 0 : if (IZZ /= KZ(1)) then
149 0 : write(6,6001)zlab(1),izz,nn,kz(1)
150 0 : ierr=1
151 0 : GOTO 1001
152 : end if
153 0 : NTOTP=NF
154 0 : if (NTOTP > nptot) then
155 0 : write(6,6002)ntotp,nptot
156 0 : ierr=2
157 0 : GOTO 1001
158 : end if
159 0 : INT(1)=1
160 0 : if (ITTE3 /= ITE3) then
161 0 : write(6,6077)ite3,itte3,nn
162 0 : ierr=3
163 0 : GOTO 1001
164 : end if
165 : ! ITE1=MAX(ITE1,ITTE1)
166 : ! ITE2=MIN(ITE2,ITTE2)
167 : !
168 : ! READ MESH FILES
169 0 : OPEN(1,FILE=trim(path)//'/'//ZLAB(1)//'.mesh',status='old',form='unformatted',iostat=ios)
170 0 : if (ios /= 0) then
171 0 : write(*,*) 'failed to open ' // trim(path) // '/' // ZLAB(1) // '.mesh'
172 0 : ierr = -1
173 0 : GOTO 1001
174 : end if
175 0 : read(1)DV,NTOTV,(UMESH(N),N=1,NTOTV)
176 0 : umin=umesh(1)
177 0 : umax=umesh(ntotv)
178 0 : DV1=DV
179 0 : CLOSE(1)
180 : !
181 : ! GET MESH FOR SCREEN
182 0 : call IMESH(UMESH,NTOTV)
183 : !
184 : ! SUBSEQUENT FILES
185 0 : do N=2,ipe
186 0 : NN=N
187 0 : OPEN(1,FILE=trim(path)//'/'//ZLAB(N)//'.index',STATUS='OLD')
188 0 : read(1,*)IZZ,AMM
189 0 : read(1,*)ITE11,ITE22,ITE33
190 0 : read(1,*)UMINN,UMAXX
191 0 : read(1,*)NC,NF
192 0 : read(1,*)DPACK
193 0 : CLOSE(1)
194 0 : if (ITE33 /= ITE3) then
195 0 : write(6,6077)ite3,ite33,nn
196 0 : ierr=4
197 0 : GOTO 1001
198 : end if
199 : ! ITE1=MAX(ITE1,ITE11)
200 : ! ITE2=MIN(ITE2,ITE22)
201 0 : if (IZZ /= KZ(N)) then
202 0 : write(6,6001)zlab(n),izz,nn,kz(nn)
203 0 : ierr=5
204 0 : GOTO 1001
205 : end if
206 0 : NTOTT=NF
207 0 : if (NTOTT > NTOTP) then
208 0 : write(6,6006)nn,ntott,ntotp
209 0 : ierr=6
210 0 : GOTO 1001
211 : end if
212 : ! !! if (UMIN /= UMINN.OR.UMAX /= UMAXX) GOTO 1003 !!
213 0 : INT(N)=NTOTP/NTOTT
214 0 : if (INT(N)*NTOTT /= NTOTP) then
215 0 : WRITE(6,6009)NN,NTOTT,NTOTP
216 0 : ierr=7
217 0 : GOTO 1001
218 : end if
219 0 : if (INT(N) /= 1)WRITE(6,6007)N,INT(N)
220 : !
221 : ! READ MESH FILES
222 : !
223 0 : OPEN(1,FILE=trim(path)//'/'//ZLAB(N)//'.mesh', status='old',form='unformatted',iostat=ios)
224 0 : if (ios /= 0) then
225 0 : write(*,*) 'failed to open ' // trim(path) // '/' // ZLAB(N) // '.mesh'
226 0 : ierr = -1
227 0 : GOTO 1001
228 : end if
229 0 : read(1)DV
230 0 : if (DV /= DV1) then
231 : ! write(*,*) ' OP: N=',N,', DV=',DV,' NOT EQUAL TO DV1=',DV1
232 0 : ierr=8
233 0 : GOTO 1001
234 : end if
235 0 : CLOSE(1)
236 : end do
237 : !
238 : ! START TEMPERATURE LOOP
239 : !
240 0 : ncount2=0
241 0 : ncount3=0
242 0 : do it=ite1,ite2,ite3
243 : !
244 : ! OPEN FILES
245 : !
246 0 : TLAB='.'//NUM(IT/100)//NUM(IT/10-10*(IT/100))//NUM(IT-10*(IT/10))
247 0 : do n=1,ipe
248 : ! if (SKIP(N))GOTO 70
249 0 : NN=N
250 0 : OPEN(IFL(N),FILE=trim(path)//'/'//ZLAB(N)//TLAB,FORM='UNFORMATTED',STATUS='OLD',iostat=ios)
251 0 : if (ios /= 0) then
252 0 : write(*,*) 'failed to open ' // trim(path)//'/'//ZLAB(N)//TLAB
253 0 : ierr = -1
254 0 : GOTO 1001
255 : end if
256 0 : if (n > 2) then
257 0 : OPEN(IFLP(N),FILE=trim(path)//'/'//ZLABP(N)//TLAB,FORM='UNFORMATTED',STATUS='OLD',iostat=ios)
258 0 : if (ios /= 0) then
259 0 : write(*,*) 'failed to open ' // trim(path)//'/'//ZLABP(N)//TLAB
260 0 : ierr = -1
261 0 : GOTO 1001
262 : end if
263 : end if
264 : end do
265 : ! READ HEADINGS
266 0 : NN=1
267 0 : read(IFL(1))IZZ,ITE,AM,UM,UX,NCCC,NFFF,DelP,JNE1,JNE2,JNE3
268 0 : do n=2,ipe
269 : ! if (SKIP(N))GOTO 80
270 0 : NN=N
271 0 : read(IFL(N))IZZ,ITE,AM,UM,UX,NC,NF,DelP,JNE11,JNE22,JNE33
272 0 : if (n > 2) read(iflp(n))
273 0 : if (JNE33 /= JNE3) then
274 0 : write(6,6099)jne3,jne33,nn
275 0 : ierr=9
276 0 : GOTO 1001
277 : end if
278 0 : JNE1=MAX(JNE1,JNE11)
279 0 : JNE2=MIN(JNE2,JNE22)
280 : end do
281 0 : itt=(it-ite1)/2+1
282 0 : jn1(itt)=jne1
283 0 : jn2(itt)=jne2
284 : !
285 : ! WRITE(98,9802)ITE,JNE1,JNE2,JNE3
286 : !
287 : ! START DENSITY LOOP
288 : !
289 0 : do n=1,ipe
290 0 : do jn=jne1,jne2,jne3
291 0 : jnn=(jn-jne1)/2+1
292 : !
293 : ! START LOOP ON ELEMENTS
294 : !
295 0 : 95 read(IFL(N))JNE,EPATOM(n,itt,jnn),OPLNCK(n,itt,jnn),ORSS,
296 0 : & NE1P(n,itt,jnn),NE2P(n,itt,jnn),
297 0 : & (FIONP(NE,n,itt,jnn),NE=NE1P(n,itt,jnn),NE2P(n,itt,jnn))
298 0 : read(ifl(n))np(n,itt,jnn)
299 0 : if (np(n,itt,jnn) > 0) then
300 0 : read(ifl(n))(nx_temp(k),y_temp(k),k=1,np(n,itt,jnn))
301 0 : do nx_index = 2, np(n,itt,jnn)
302 0 : left_val = y_temp(nx_index-1)
303 0 : right_val = y_temp(nx_index)
304 0 : left_n = nx_temp(nx_index-1)
305 0 : right_n = nx_temp(nx_index)
306 0 : slope = (right_val - left_val)/float(right_n - left_n)
307 :
308 0 : do n_index = left_n, right_n
309 0 : cross_section = left_val + (n_index-left_n)*slope
310 0 : yy2(ncount2 + n_index) = cross_section
311 : end do
312 0 : yy2(ncount2 + left_n) = left_val
313 0 : yy2(ncount2 + right_n) = right_val
314 : end do
315 0 : kp2(n, itt, jnn) = ncount2
316 0 : ncount2 = ncount2 + ntotp
317 : else
318 0 : read(ifl(n))(yy2(k+ncount2),k=1,ntotp)
319 0 : kp2(n,itt,jnn)=ncount2
320 0 : ncount2=ncount2+ntotp
321 : end if
322 0 : if (n > 2) then
323 0 : read(iflp(n))ja,npp(n,itt,jnn)
324 0 : if (npp(n,itt,jnn) > 0) then
325 0 : read(iflp(n))(nx(k+ncount3),yx(k+ncount3),k=1,npp(n,itt,jnn))
326 0 : kp3(n,itt,jnn)=ncount3
327 0 : ncount3=ncount3+npp(n,itt,jnn)
328 : end if
329 : end if
330 : end do
331 : end do
332 :
333 : ! write(6,610)it
334 : ! write(6,*)'ncount1 = ',ncount1
335 : ! write(6,*)'ncount2 = ',ncount2
336 : ! write(6,*)'ncount3 = ',ncount3
337 : !
338 : ! CLOSE FILES
339 :
340 0 : do N=1,ipe
341 0 : close(IFL(N))
342 0 : close(iflp(n))
343 : end do
344 :
345 : end do
346 :
347 0 : write(*,*) 'done loading OP mono data'
348 0 : have_loaded_op = .true.
349 :
350 : !write(*,*)'ncount1 = ',ncount1
351 : !write(6,*)'ncount2 = ',ncount2
352 : !write(6,*)'ncount3 = ',ncount3
353 0 : ios = 0
354 0 : open(1, file=trim(cache_filename), iostat=ios, action='write', form='unformatted')
355 0 : if (ios == 0) then
356 0 : write(*,*) 'write ' // trim(cache_filename)
357 0 : write(1) op_cache_version, ntotv,dv,dv1,umesh,
358 0 : & ite1,ite2,ite3,jn1,jn2,jne3,umin,umax,ntotp,nc,nf,int,epatom,oplnck, ne1p,
359 0 : & ne2p,fionp,np,kp1,kp2,kp3,npp,yy2,nx,yx
360 0 : close(1)
361 : end if
362 :
363 : 1001 continue
364 :
365 : ! pre-calculate semesh
366 0 : do n = 1, nptot
367 0 : u = umesh(n)
368 0 : semesh(n) = 1.d0 - exp(dble(-u))
369 : end do
370 :
371 : !$omp end critical (critical_do_op_dload)
372 :
373 :
374 : return
375 : 610 format(10x,'Done IT= ',i3)
376 : 1004 WRITE(6,6004)ZLAB(NN),TLAB
377 : STOP
378 : 6001 FORMAT(//5X,'*** OP: FILE ',A3,' GIVES IZZ=',I3,'NOT EQUAL TO IZ(',I2,')=',I2,' ***')
379 : 6002 FORMAT(//5X,'*** OP: NTOT=',I7,' GREATER THAN nptot=',I7,' ***')
380 : 6003 FORMAT(//5X,'*** OP: DISCREPANCY BETWEEN DATA ON FILES ',A3,' AND ',A3,' ***')
381 : 6004 FORMAT(//5X,'*** OP: ERROR OPENING FILE ',A3,A6,' ***')
382 : 6006 FORMAT(//5X,'OP: N=',I2,', NTOTT=',I7,', GREATER THAN NTOT=',I7)
383 : 6007 FORMAT(/5X,'OP: N=',I2,', INT(N)=',I4)
384 : 6009 FORMAT(' OP: N=',I5,', NTOTT=',I10,', NTOT=',I10/' NTOT NOT MULTIPLE OF NTOTT')
385 : ! 6012 FORMAT(/10X,'ERROR, SEE WRITE(6,6012)'/10X,'IT=',I3,', JN=',I3,', N=',I3,', JNE=',I3/)
386 : 6077 FORMAT(//5X,'OP: DISCREPANCY IN ITE3'/10X,I5,' READ FROM UNIT 5'/10X,I5,' FROM INDEX FILE ELEMENT',I5)
387 : 6099 FORMAT(//5X,'OP: DISCREPANCY IN JNE3'/10X,I5,' READ FOR N=1'/10X,I5,' READ FOR N=',I5)
388 :
389 : ! 8000 FORMAT(5X,I5,F10.4/5X,3I5/2E10.2/2I10/10X,E10.2)
390 : 1010 write(*,*) ' OP: ERROR OPENING FILE '//'./'//ZLAB(1)//'.index'
391 : stop
392 : 1011 write(*,*) ' OP: ERROR OPENING FILE '//'./'//ZLAB(1)//'.mesh'
393 : stop
394 : end subroutine op_dload
395 :
396 : ! **********************************************************************
397 0 : subroutine IMESH(UMESH,NTOT)
398 :
399 : DIMENSION UMESH(nptot)
400 : COMMON/CIMESH/U(100),AA(nptot),BB(nptot),IN(nptot),ITOT,NN
401 : save /cimesh/
402 :
403 0 : UMIN=UMESH(1)
404 0 : UMAX=UMESH(NTOT)
405 :
406 0 : II=100
407 0 : A=(II*UMIN-UMAX)/REAL(II-1)
408 0 : B=(UMAX-UMIN)/REAL(II-1)
409 0 : do I=1,II
410 0 : U(I)=A+B*I
411 : end do
412 :
413 0 : ib=2
414 0 : ub=u(ib)
415 0 : ua=u(ib-1)
416 0 : d=ub-ua
417 0 : ibb=0
418 0 : do n=2,ntot
419 0 : if (umesh(n) > ub) then
420 0 : ua=ub
421 0 : ib=ib+1
422 0 : ub=u(ib)
423 0 : d=ub-ua
424 0 : if (umesh(n) > ub) then
425 0 : nn=n-1
426 0 : ibb=ib-1
427 0 : GOTO 1
428 : end if
429 : end if
430 0 : in(n)=ib
431 0 : aa(n)=(ub-umesh(n))/d
432 0 : bb(n)=(umesh(n)-ua)/d
433 : end do
434 :
435 0 : 1 ib=ibb
436 0 : do n=nn+1,ntot
437 0 : ib=ib+1
438 0 : in(n)=ib
439 0 : u(ib)=umesh(n)
440 : end do
441 0 : itot=ib
442 :
443 0 : end subroutine IMESH
444 :
445 :
446 0 : subroutine msh(dv, ntot, umesh, semesh, uf, dscat)
447 : implicit none
448 : integer, intent(out) :: ntot
449 : real, intent(out) :: dv, uf(0:100), dscat
450 : real, intent(out) :: umesh(:), semesh(:) ! (nptot)
451 : integer :: i, ntotv
452 0 : real :: dvp, dv1, umin, umax, umeshp(nptot), semeshp(nptot)
453 : common /mesh/ ntotv, dvp, dv1, umeshp, semeshp
454 : save /mesh/
455 :
456 0 : ntot = ntotv
457 0 : dv = dvp
458 0 : do i=1,ntot
459 0 : umesh(i) = umeshp(i)
460 : end do
461 0 : do i=1,ntot
462 0 : semesh(i) = semeshp(i)
463 : end do
464 :
465 0 : umin = umesh(1)
466 0 : umax = umesh(ntot)
467 0 : dscat = (umax - umin)*0.01
468 0 : do i = 0, 100
469 0 : uf(i) = umin + i*dscat
470 : end do
471 :
472 0 : end subroutine msh
473 :
474 :
475 0 : subroutine solve(u,v,z,uz,ierr)
476 : integer, intent(inout) :: ierr
477 : dimension u(4)
478 :
479 : ! If P(R) = u(1) u(2) u(3) u(4)
480 : ! for R = -3 -1 1 3
481 : ! then a cubic fit is:
482 0 : P(R)=(
483 : & 27*(u(3)+u(2))-3*(u(1)+u(4)) +R*(
484 : & 27*(u(3)-u(2))-(u(4)-u(1)) +R*(
485 : & -3*(u(2)+u(3))+3*(u(4)+u(1)) +R*(
486 : & -3*(u(3)-u(2))+(u(4)-u(1)) ))))/48.
487 : ! First derivative is:
488 : PP(R)=(
489 : & 27*(u(3)-u(2))-(u(4)-u(1))+ 2*R*(
490 : & -3*(u(2)+u(3))+3*(u(4)+u(1)) +3*R*(
491 : & -3*(u(3)-u(2))+(u(4)-u(1)) )))/48.
492 :
493 : ! ierr = 0
494 : ! Find value of z giving P(z)=v
495 : ! First estimate
496 0 : z=(2.*v-u(3)-u(2))/(u(3)-u(2))
497 : ! Newton-Raphson iterations
498 0 : do k=1,10
499 0 : uz=pp(z)
500 0 : d=(v-p(z))/uz
501 0 : z=z+d
502 0 : if (abs(d) < 1.e-4) return
503 : end do
504 :
505 : ! write(*,*) ' Not converged after 10 iterations in SOLVE'
506 : ! write(*,*) ' v=',v
507 : ! do N=1,4
508 : ! write(*,*) ' N, U(N)=',N,U(N)
509 : ! end do
510 0 : ierr = 10
511 0 : return
512 : ! stop
513 :
514 : end subroutine solve
515 : ! **********************************************************************
516 :
517 0 : subroutine BRCKR(T,FNE,RION,NION,U,NFREQ,SF, ierr)
518 : integer, intent(inout) :: ierr
519 : !
520 : ! CODE FOR COLLECTIVE EFFECTS ON THOMSON SCATTERING.
521 : ! METHOD OF D.B. BOERCKER, AP. J., 316, L98, 1987.
522 : !
523 : ! INPUT:-
524 : ! T=TEMPERATURTE IN K
525 : ! FNE=ELECTRON DENSITY IN CM**(-3)
526 : ! ARRAY RION (DIMENSIONED FOR 30 IONS).
527 : ! RION(IZ) IS NUMBER OF IONS WITH NET CHARGE IZ.
528 : ! NORMALISATION OF RION IS OF NO CONSEQUENCE.
529 : ! NION=NUMBER OF IONS INCLUDED.
530 : ! ARRAY U (DIMENSIONED FOR 1000). VALUES OF (H*NU/K*T).
531 : ! NFREQ=NUMBER OF FREQUENCY POINTS.
532 : !
533 : ! OUTPUT:-
534 : ! ARRAY SF, GIVING FACTORS BY WHICH THOMSON CROSS SECTION
535 : ! SHOULD BE MULTIPLIED TO ALLOW FOR COLLECTIVE EFFECTS.
536 : !
537 : ! MODIFFICATIONS:-
538 : ! (1) REPLACE (1.-Y) BY EXP(-Y) TO AVOID NEGATIVE FACTORS FOR
539 : ! HIGHLY-DEGENERATE CASES.
540 : ! (2) INCLUDE RELATIVISTIC CORRECTION.
541 : !
542 : parameter (IPZ=28,IPNC=100)
543 : DIMENSION RION(IPZ),U(0:IPNC),SF(0:IPNC)
544 :
545 0 : AUNE=1.48185E-25*FNE
546 0 : AUT=3.16668E-6*T
547 0 : C1=-1.0650E-4*AUT
548 0 : C2=+1.4746E-8*AUT**2
549 0 : C3=-2.0084E-12*AUT*AUT*AUT
550 0 : V=7.8748*AUNE/(AUT*sqrt(AUT))
551 0 : call FDETA(V,ETA, ierr) ! 23.10.93
552 0 : W=exp(dble(ETA)) ! 23.10.93
553 0 : 11 R=FMH(W)/V
554 0 : A=0.
555 0 : B=0.
556 0 : do I=1,NION
557 0 : A=A+I*RION(I)
558 0 : B=B+I**2*RION(I)
559 : end do
560 0 : X=R+B/A
561 :
562 0 : Y=.353553*W
563 0 : C=1.1799E5*X*AUNE/(AUT*AUT*AUT)
564 0 : do N=0,NFREQ
565 0 : D=C/U(N)**2
566 0 : if (D > 5.) then
567 0 : D=-2./D
568 0 : F=2.666667*(1.+D*(.7+D*(.55+.341*D)))
569 : else
570 0 : G=2.*D*(1+D)
571 0 : F=D*((G+D*D*D)*log(dble(D/(2.+D)))+G+2.6666667)
572 : end if
573 0 : DELTA=.375*R*F/X
574 : SF(N)=(1.-R*DELTA-Y*FUNS(W))*
575 0 : & (1.+U(N)*(C1+U(N)*(C2+U(N)*C3))) !SAMPSON CORRECTION
576 : end do
577 :
578 0 : return
579 :
580 : 600 FORMAT(5X,'NOT CONVERGED IN LOOP 10 OF BRCKR'/5X,'T=',1P,E10.2,', FNE=',E10.2)
581 :
582 : end subroutine BRCKR
583 : ! **********************************************************************
584 0 : function FUNS(A)
585 : !
586 0 : if (A <= 0.001) then
587 : FUNS=1.
588 0 : else if (A <= 0.01) then
589 : FUNS=(1.+A*(-1.0886+A*(1.06066+A*1.101193)))/
590 0 : & (1.+A*(0.35355+A*(0.19245+A+0.125)))
591 : else
592 : FUNS=( 1./(1.+0.81230*A)**2+
593 : & 0.92007/(1.+0.31754*A)**2+
594 : & 0.05683/(1.+0.04307*A)**2 )/
595 : & ( 1./(1.+0.65983*A)+
596 : & 0.92007/(1.+0.10083*A)+
597 0 : & 0.05683/(1.+0.00186*A) )
598 : end if
599 :
600 0 : end function FUNS
601 : ! **********************************************************************
602 0 : function FMH(W)
603 : !
604 : ! CALCULATES FD INTEGRAL I_(-1/2)(ETA). INCLUDES FACTOR 1/GAMMA(1/2).
605 : ! ETA=log(W)
606 : !
607 0 : if (W <= 2.718282) then
608 0 : FMH=W*(1+W*(-.7070545+W*(-.3394862-W*6.923481E-4))/(1.+W*(1.2958546+W*.35469431)))
609 0 : else if (W <= 54.59815) then
610 0 : X=log(dble(W))
611 0 : FMH=(.6652309+X*(.7528360+X*.6494319))/(1.+X*(.8975007+X*.1153824))
612 : else
613 0 : X=log(dble(W))
614 0 : Y=1./X**2
615 0 : FMH=sqrt(X)*(1.1283792+(Y*(-.4597911+Y*(2.286168-Y*183.6074)))/(1.+Y*(-10.867628+Y*384.61501)))
616 : end if
617 :
618 0 : end function FMH
619 : ! **********************************************************************
620 0 : subroutine FDETA(X,ETA, ierr)
621 : !
622 : ! GIVEN X=N_e/P_e, CALCULATES FERMI-DIRAC ETA
623 : ! USE CHEBYSHEV FITS OF W.J. CODY AND H.C. THACHER,
624 : ! MATHS. OF COMP., 21, 30, 1967.
625 : !
626 : integer, intent(inout) :: ierr
627 : DIMENSION D(2:12)
628 : DATA D/
629 : & 3.5355339E-01, 5.7549910E-02, 5.7639604E-03, 4.0194942E-04,
630 : & 2.0981899E-05, 8.6021311E-07, 2.8647149E-08, 7.9528315E-10,
631 : & 1.8774422E-11, 3.8247505E-13, 6.8427624E-15/
632 :
633 : integer :: n, k
634 :
635 : ! ierr = 0
636 0 : a=x*0.88622693
637 :
638 0 : if (X < 1) then
639 0 : v=x
640 0 : S=V
641 0 : U=V
642 0 : do N=2,12
643 0 : S=S*V
644 0 : SS=S*D(N)
645 0 : U=U+SS
646 0 : if (ABS(SS) < 1.E-6*U)GOTO 11
647 : end do
648 : ! write(*,*) ' COMPLETED LOOP 10 IN FDETA'
649 0 : ierr = 11
650 0 : return
651 : ! STOP
652 0 : 11 ETA=log(dble(U))
653 :
654 : else
655 0 : if (a < 2) then
656 0 : E=log(dble(X))
657 : else
658 0 : e=pow(1.5d0*a,2d0/3d0)
659 : end if
660 0 : do k=1,10
661 0 : call FDF1F2(E,F1,F2)
662 0 : DE=(A-F2)*2./F1
663 0 : E=E+DE
664 0 : if (abs(dE) < 1.e-4*abs(E))GOTO 21
665 : end do
666 : ! write(*,*) ' completed loop 20 IN FDETA'
667 0 : ierr = 12
668 0 : return
669 : ! stop
670 0 : 21 ETA=E
671 :
672 : end if
673 :
674 : end subroutine FDETA
675 : ! **********************************************************************
676 0 : subroutine FDF1F2(ETA,F1,F2)
677 : !
678 : ! CALCULATES FD INTEGRALS F1, F2=F(-1/2), F(+1/2)
679 : ! USE CHEBYSHEV FITS OF W.J. CODY AND H.C. THACHER,
680 : ! MATHS. OF COMP., 21, 30, 1967.
681 : !
682 0 : if (ETA <= 1) then
683 0 : X=exp(dble(ETA))
684 : F1=X*(1.772454+X*(-1.2532215+X*(-0.60172359-X*0.0012271551))/
685 0 : & (1.+X*(1.2958546+X*0.35469431)))
686 : F2=X*(0.88622693+X*(-0.31329180+X*(-0.14275695-
687 : & X*0.0010090890))/
688 0 : & (1.+X*(0.99882853+X*0.19716967)))
689 0 : else if (ETA <= 4) then
690 0 : X=ETA
691 : F1=(1.17909+X*(1.334367+X*1.151088))/
692 0 : & (1.+X*(0.8975007+X*0.1153824))
693 : F2=(0.6943274+X*(0.4918855+X*0.214556))/
694 0 : & (1.+X*(-0.0005456214+X*0.003648789))
695 : else
696 0 : X=sqrt(ETA)
697 0 : Y=1./ETA**2
698 : F1=X*(2.+Y*(-0.81495847+Y*(4.0521266-Y*325.43565))/
699 0 : & (1.+Y*(-10.867628+Y*384.61501)))
700 : F2=ETA*X*(0.666666667+Y*(0.822713535+Y*(5.27498049+
701 : & Y*290.433403))/
702 0 : & (1.+Y*(5.69335697+Y*322.149800)))
703 : end if
704 :
705 0 : end subroutine FDF1F2
706 :
707 :
708 0 : subroutine screen2(ft,fne,rion,epa,ntot,umin,umax,umesh,p)
709 : parameter(ipz=28)
710 : real :: umesh(:) ! (nptot)
711 : real :: p(:) ! (nptot)
712 0 : dimension rion(ipz),f(100)
713 : dimension x(3),wt(3)
714 : data x/0.415775,2.294280,6.289945/
715 : data wt/0.711093,0.278518,0.0103893/
716 : data twopi/6.283185/
717 : COMMON/CIMESH/U(100),AA(nptot),BB(nptot),IN(nptot),ITOT,NN
718 : save /cimesh/
719 :
720 0 : rydt=ft/157894.
721 0 : aune=1.48185e-25*fne
722 :
723 : ! get alp2=1/(Debye)**2
724 0 : b=0
725 0 : do i=1,ipz
726 0 : b=b+rion(i)*i**2
727 : end do
728 0 : alp2=(5.8804e-19)*fne*b/(epa*ft)
729 0 : if (alp2/ft < 5e-8) return !!!!!!!!!!!
730 :
731 0 : c=1.7337*aune/sqrt(rydt)
732 :
733 0 : do i=1,itot
734 0 : w=u(i)*rydt
735 0 : f(i)=0.
736 0 : do k=1,ipz
737 0 : if (rion(k) <= 0.01) cycle
738 0 : crz=c*rion(k)*k**2
739 0 : ff=0
740 0 : do j=1,3
741 0 : e=x(j)*rydt
742 0 : fk=sqrt(e)
743 0 : fkp=sqrt(e+w)
744 0 : x1=1.+alp2/(fkp+fk)**2
745 0 : x2=1.+alp2/(fkp-fk)**2
746 0 : q=(1./x2-1./x1+log(dble(x1/x2)))*
747 0 : & (fkp*(1.-exp(dble(-twopi*k/fkp))))/(fk*(1.-exp(dble(-twopi*k/fk))))
748 0 : ff=ff+wt(j)*q
749 : end do
750 0 : f(i)=f(i)+crz*ff
751 : end do
752 : end do
753 :
754 0 : p(1)=f(1)
755 0 : do n=2,nn
756 0 : w=umesh(n)*rydt
757 0 : p(n)=p(n)+(aa(n)*f(in(n)-1)+bb(n)*f(in(n)))/(w*w*w)
758 : end do
759 0 : do n=nn+1,ntot
760 0 : w=umesh(n)*rydt
761 0 : p(n)=p(n)+f(in(n))/(w*w*w)
762 : end do
763 :
764 : end subroutine screen2
765 :
766 : end module op_load
|
