Line data Source code
1 : module op_common
2 :
3 : use math_lib
4 : use op_def
5 :
6 : contains
7 : ! **********************************************************************
8 0 : subroutine xindex(flt, ilab, xi, ih, i3, ierr)
9 : implicit none
10 : integer, intent(in) :: i3
11 : real, intent(in) :: flt
12 : integer, intent(out) :: ih(4), ilab(4)
13 : real, intent(out) :: xi
14 : integer, intent(out) :: ierr
15 : integer :: i, ih2
16 0 : real :: x
17 :
18 0 : ierr = 0
19 0 : if(flt<3.5) then
20 0 : ierr = 102
21 0 : return
22 0 : else if(flt>8.) then
23 0 : ierr = 102
24 0 : return
25 : end if
26 :
27 0 : x = 40.*flt/real(i3)
28 0 : ih2 = x
29 0 : ih2 = max(ih2, 140/i3+2)
30 0 : ih2 = min(ih2, 320/i3-3)
31 0 : do i = 1, 4
32 0 : ih(i) = ih2 + i - 2
33 0 : ilab(i) = i3*ih(i)
34 : end do
35 0 : xi = 2.*(x-ih2) - 1
36 :
37 : end subroutine xindex
38 :
39 : ! *********************************************************************
40 0 : subroutine jrange(ih, jhmin, jhmax, i3)
41 : implicit none
42 : integer, intent(in) :: ih(4), i3
43 : integer, intent(out) :: jhmin, jhmax
44 : integer :: i
45 :
46 0 : jhmin = 0
47 0 : jhmax = 1000
48 0 : do i = 1, 4
49 0 : jhmin = max(jhmin, js(ih(i)*i3)/i3)
50 0 : jhmax = min(jhmax, je(ih(i)*i3)/i3)
51 : end do
52 :
53 0 : end subroutine jrange
54 : ! *********************************************************************
55 :
56 0 : subroutine findne(ilab, fa, nel, nkz, jhmin, jhmax, ih, flrho, flt, xi, flne, flmu, flr, epa, uy, i3, ierr)
57 : use op_load, only : solve
58 : implicit none
59 : integer, intent(in) :: ilab(4), nel, nkz(ipe), jhmin, ih(4), i3
60 : integer, intent(inout) :: jhmax
61 : integer, intent(out) :: ierr
62 : real, intent(in) :: fa(ipe), flt, xi, flmu
63 : real,intent(out) :: flne, uy, epa
64 : real, intent(inout) :: flrho
65 : integer :: i, j, n, jh, jm, itt, jne, jnn
66 0 : real :: flrmin, flrmax, flr(4,4), uyi(4), efa(4, 7:118), flrh(4, 7:118), u(4), flnei(4), y, zeta, efa_temp
67 : ! declare variables in common block, by default: real (a-h, o-z), integer (i-n)
68 : ! integer :: ite1, ite2, ite3, jn1, jn2, jne3, ntot, nc, nf, int, ne1, ne2, np, kp1, kp2, kp3, npp, mx, nx
69 : ! real :: umin, umax, epatom, oplnck, fion, yy1, yy2, yx
70 : ! common /atomdata/ ite1,ite2,ite3,jn1(91),jn2(91),jne3,umin,umax,ntot, &
71 : ! nc,nf,int(17),epatom(17,91,25),oplnck(17,91,25),ne1(17,91,25), &
72 : ! ne2(17,91,25),fion(-1:28,28,91,25),np(17,91,25),kp1(17,91,25), &
73 : ! kp2(17,91,25),kp3(17,91,25),npp(17,91,25),mx(33417000), &
74 : ! yy1(33417000),yy2(120000000),nx(19305000),yx(19305000)
75 : ! save /atomdata/
76 :
77 : ! efa(i,jh)=sum_n epa(i,jh,n)*fa(n)
78 : ! flrh(i,jh)=log10(rho(i,jh))
79 :
80 : ! Get efa
81 0 : do i = 1, 4
82 0 : itt = (ilab(i)-ite1)/2 + 1
83 0 : do jne = jn1(itt), jn2(itt), i3
84 0 : jnn = (jne-jn1(itt))/2 + 1
85 0 : jh = jne/i3
86 0 : efa_temp = 0.
87 0 : do n = 1, nel
88 0 : efa_temp = efa_temp + epatom(nkz(n), itt, jnn)*fa(n)
89 : end do !n
90 0 : efa(i, jh) = efa_temp
91 : end do !jne
92 : end do !i
93 :
94 : ! Get range for efa > 0
95 0 : do i = 1, 4
96 0 : do jh = jhmin, jhmax
97 0 : if (efa(i, jh) <= 0.) then
98 0 : jm = jh - 1
99 : GOTO 3
100 : end if
101 : end do
102 : GOTO 4
103 0 : 3 jhmax = min(jhmax, jm)
104 0 : 4 continue
105 : end do
106 :
107 : ! Get flrh
108 0 : do jh=jhmin,jhmax
109 0 : do i=1,4
110 0 : flrh(i, jh) = flmu + 0.25*i3*jh - log10(dble(efa(i, jh)))
111 : end do
112 : end do
113 :
114 : ! Find flrmin and flrmax
115 0 : flrmin = -1000
116 0 : flrmax = 1000
117 0 : do i = 1, 4
118 0 : flrmin = max(flrmin, flrh(i,jhmin))
119 0 : flrmax = min(flrmax, flrh(i,jhmax))
120 : end do
121 :
122 : ! Check range of flrho
123 0 : if (flrho < flrmin .or. flrho > flrmax) then
124 0 : write(*,*) "findne failed because density is out of range for logT, logRho", flt, flrho
125 0 : write(*,*) "Allowed range for logRho is",flrmin," to ",flrmax
126 0 : ierr = 101
127 0 : return
128 : end if
129 :
130 : ! Interpolations in j for flne
131 0 : do jh = jhmin, jhmax
132 0 : if (flrh(2,jh) > flrho) then
133 0 : jm = jh - 1
134 : GOTO 5
135 : end if
136 : end do
137 0 : write(*,*) ' Interpolations in j for flne'
138 0 : write(*,*) ' Not found, i=',i
139 0 : stop
140 0 : 5 jm=max(jm,jhmin+1)
141 0 : jm=min(jm,jhmax-2)
142 :
143 0 : do i = 1, 4
144 0 : do j = 1, 4
145 0 : u(j) = flrh(i, jm+j-2)
146 0 : flr(i,j) = flrh(i, jm+j-2)
147 : end do
148 0 : call solve(u, flrho, zeta, uyi(i), ierr)
149 0 : if (ierr /= 0) return
150 0 : y = jm + 0.5*(zeta+1)
151 0 : flnei(i) = .25*i3*y
152 : end do
153 :
154 : ! Interpolations in i
155 0 : flne = fint(flnei, xi)
156 0 : uy = fint(uyi, xi)
157 : ! Get epa
158 0 : epa = exp10(dble(flne + flmu - flrho))
159 :
160 0 : return
161 :
162 : ! 601 format(' For flt=',1p,e11.3,', flrho=',e11.3,' is out of range. Allowed range for flrho is ',e11.3,' to ',e11.3)
163 0 : end subroutine findne
164 :
165 : ! **********************************************************************
166 0 : subroutine yindex(jhmin, jhmax, flne, jh, i3, eta)
167 : implicit none
168 : integer, intent(in) :: jhmin, jhmax, i3
169 : real, intent(in) :: flne
170 : integer, intent(out) :: jh(4)
171 : real, intent(out) :: eta
172 : integer :: j, k
173 0 : real :: y
174 :
175 0 : y = 4.*flne/real(i3)
176 0 : j = y
177 0 : j = max(j,jhmin+2)
178 0 : j = min(j,jhmax-3)
179 0 : do k = 1, 4
180 0 : jh(k) = j + k - 2
181 : end do
182 0 : eta = 2.*(y-j)-1
183 :
184 0 : end subroutine yindex
185 :
186 : ! **********************************************************************
187 0 : subroutine findux(flr, xi, eta, ux)
188 : implicit none
189 : real, intent(in) :: flr(4, 4), xi, eta
190 : real, intent(out) :: ux
191 : integer :: i, j
192 0 : real :: uxj(4), u(4)
193 :
194 0 : do j = 1, 4
195 0 : do i = 1, 4
196 0 : u(i) = flr(i, j)
197 : end do
198 0 : uxj(j) = fintp(u, xi)
199 : end do
200 0 : ux = fint(uxj, eta)
201 :
202 0 : end subroutine findux
203 :
204 : ! *************************************
205 0 : function fint(u,r)
206 : dimension u(4)
207 :
208 : ! If P(R) = u(1) u(2) u(3) u(4)
209 : ! for R = -3 -1 1 3
210 : ! then a cubic fit is:
211 : P(R)=(27*(u(3)+u(2))-3*(u(1)+u(4)) +R*(27*(u(3)-u(2))-(u(4)-u(1))
212 : & + R*(-3*(u(2)+u(3))+3*(u(4)+u(1)) +R*(-3*(u(3)-u(2))+(u(4)-u(1)) ))))/48.
213 :
214 0 : fint=p(r)
215 :
216 0 : end function fint
217 : ! **********************************************************************
218 0 : function fintp(u,r)
219 : dimension u(4)
220 :
221 : ! If P(R) = u(1) u(2) u(3) u(4)
222 : ! for R = -3 -1 1 3
223 : ! then a cubic fit to the derivative is:
224 : PP(R)=(27*(u(3)-u(2))-(u(4)-u(1)) +2.*R*(-3*(u(2)+u(3))+3*(u(4)+u(1)) +3.*R*(-3*(u(3)-u(2))+(u(4)-u(1)) )))/48.
225 :
226 0 : fintp=pp(r)
227 :
228 0 : end function fintp
229 :
230 : ! **********************************************************************
231 0 : subroutine scatt(ih, jh, rion, uf, f, umesh, semesh, dscat, ntot, epa, ierr)
232 : use op_load, only: BRCKR
233 : integer, intent(inout) :: ierr
234 : real :: umesh(:), semesh(:) ! (nptot)
235 : real :: f(:,:,:) ! (nptot,4,4)
236 0 : dimension rion(28, 4, 4),uf(0:100),fscat(0:100),p(nptot),rr(28),ih(4),jh(4)
237 : integer :: i,j,k,n
238 : ! HH: always use meshtype q='m'
239 0 : ite3=2
240 0 : umin=umesh(1)
241 0 : CSCAT=EPA*2.37567E-8
242 0 : ierr = 0
243 :
244 0 : do i = 1, 4
245 0 : ft = real(exp10(ITE3*dble(ih(i))/40d0))
246 0 : do j = 1, 4
247 0 : fne = real(exp10(ITE3*dble(jh(j))/4d0))
248 0 : do k = 1, ntot
249 0 : p(k) = f(k, i, j)
250 : end do
251 0 : do m = 1, 28
252 0 : rr(m) = rion(m, i, j)
253 : end do
254 0 : call BRCKR(FT,FNE,RR,28,UF,100,FSCAT,ierr)
255 0 : if (ierr /= 0) return
256 0 : do n = 0, 100
257 0 : u = uf(n)
258 0 : fscat(n) = cscat*(fscat(n)-1)
259 : end do
260 0 : do n = 2, ntot-1
261 0 : u=umesh(n)
262 0 : se=semesh(n)
263 0 : m=(u-umin)/dscat
264 0 : ua=umin+dscat*m
265 0 : ub=ua+dscat
266 0 : p(n)=p(n)+((ub-u)*fscat(m)+(u-ua)*fscat(m+1))/(dscat*se)
267 : end do
268 0 : u=umesh(ntot)
269 0 : se=semesh(ntot)
270 0 : p(ntot)=p(ntot)+fscat(100)/se
271 0 : p(1)=p(1)+fscat(1)/(1.-.5*umin)
272 0 : do k=1,ntot
273 0 : f(k,i,j)=p(k)
274 : end do
275 : end do
276 : end do
277 :
278 0 : end subroutine scatt
279 : ! **********************************************************************
280 0 : subroutine screen1(ih,jh,rion,umesh,ntot,epa,f)
281 0 : use op_load, only: screen2
282 : real :: umesh(:) ! (nptot)
283 : real, pointer :: f(:,:,:) ! (nptot,4,4)
284 : dimension uf(0:100),fscat(0:100), ih(4), jh(4)
285 : real, target :: rion(28, 4, 4)
286 : integer :: i, j, k, m
287 0 : real, pointer :: p(:), rr(:)
288 :
289 0 : ite3=2
290 0 : umin=umesh(1)
291 0 : umax=umesh(ntot)
292 :
293 0 : do i = 1, 4
294 0 : ft = real(exp10(ITE3*dble(ih(i))/40d0))
295 0 : do j = 1, 4
296 0 : fne = real(exp10(ITE3*dble(jh(j))/4d0))
297 : ! do k=1,ntot
298 0 : p => f(1:ntot,i,j)
299 : ! end do
300 : ! do m=1,28
301 0 : rr => rion(1:28,i,j)
302 : ! end do
303 0 : call screen2(ft,fne,rr,epa,ntot,umin,umax,umesh,p)
304 : ! do k=1,ntot
305 : ! f(k,i,j)=p(k)
306 : ! end do
307 : end do
308 : end do
309 :
310 0 : end subroutine screen1
311 :
312 : end module op_common
|
