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 0 : jhmax = min(jhmax, jm)
100 0 : cycle
101 : end if
102 : end do
103 : end do
104 :
105 : ! Get flrh
106 0 : do jh=jhmin,jhmax
107 0 : do i=1,4
108 0 : flrh(i, jh) = flmu + 0.25*i3*jh - log10(dble(efa(i, jh)))
109 : end do
110 : end do
111 :
112 : ! Find flrmin and flrmax
113 0 : flrmin = -1000
114 0 : flrmax = 1000
115 0 : do i = 1, 4
116 0 : flrmin = max(flrmin, flrh(i,jhmin))
117 0 : flrmax = min(flrmax, flrh(i,jhmax))
118 : end do
119 :
120 : ! Check range of flrho
121 0 : if(flrho < flrmin .or. flrho > flrmax)then
122 0 : write(*,*) "findne failed because density is out of range for logT, logRho", flt, flrho
123 0 : write(*,*) "Allowed range for logRho is",flrmin," to ",flrmax
124 0 : ierr = 101
125 0 : return
126 : end if
127 :
128 : ! Interpolations in j for flne
129 0 : do jh = jhmin, jhmax
130 0 : if (flrh(2,jh) > flrho) then
131 0 : jm = jh - 1
132 : cycle
133 : end if
134 0 : write(*,*) ' Interpolations in j for flne'
135 0 : write(*,*) ' Not found, i=',i
136 0 : stop
137 : end do
138 0 : jm=max(jm,jhmin+1)
139 0 : jm=min(jm,jhmax-2)
140 :
141 0 : do i = 1, 4
142 0 : do j = 1, 4
143 0 : u(j) = flrh(i, jm+j-2)
144 0 : flr(i,j) = flrh(i, jm+j-2)
145 : end do
146 0 : call solve(u, flrho, zeta, uyi(i), ierr)
147 0 : if (ierr /= 0) return
148 0 : y = jm + 0.5*(zeta+1)
149 0 : flnei(i) = .25*i3*y
150 : end do
151 :
152 : ! Interpolations in i
153 0 : flne = fint(flnei, xi)
154 0 : uy = fint(uyi, xi)
155 : ! Get epa
156 0 : epa = exp10(dble(flne + flmu - flrho))
157 :
158 0 : return
159 :
160 : ! 601 format(' For flt=',1p,e11.3,', flrho=',e11.3,' is out of range. Allowed range for flrho is ',e11.3,' to ',e11.3)
161 0 : end subroutine findne
162 :
163 : ! **********************************************************************
164 0 : subroutine yindex(jhmin, jhmax, flne, jh, i3, eta)
165 : implicit none
166 : integer, intent(in) :: jhmin, jhmax, i3
167 : real, intent(in) :: flne
168 : integer, intent(out) :: jh(4)
169 : real, intent(out) :: eta
170 : integer :: j, k
171 0 : real :: y
172 :
173 0 : y = 4.*flne/real(i3)
174 0 : j = y
175 0 : j = max(j,jhmin+2)
176 0 : j = min(j,jhmax-3)
177 0 : do k = 1, 4
178 0 : jh(k) = j + k - 2
179 : end do
180 0 : eta = 2.*(y-j)-1
181 :
182 0 : end subroutine yindex
183 :
184 : ! **********************************************************************
185 0 : subroutine findux(flr, xi, eta, ux)
186 : implicit none
187 : real, intent(in) :: flr(4, 4), xi, eta
188 : real, intent(out) :: ux
189 : integer :: i, j
190 0 : real :: uxj(4), u(4)
191 :
192 0 : do j = 1, 4
193 0 : do i = 1, 4
194 0 : u(i) = flr(i, j)
195 : end do
196 0 : uxj(j) = fintp(u, xi)
197 : end do
198 0 : ux = fint(uxj, eta)
199 :
200 0 : end subroutine findux
201 :
202 : ! *************************************
203 0 : function fint(u,r)
204 : dimension u(4)
205 :
206 : ! If P(R) = u(1) u(2) u(3) u(4)
207 : ! for R = -3 -1 1 3
208 : ! then a cubic fit is:
209 : P(R)=(27*(u(3)+u(2))-3*(u(1)+u(4)) +R*(27*(u(3)-u(2))-(u(4)-u(1))
210 : & + R*(-3*(u(2)+u(3))+3*(u(4)+u(1)) +R*(-3*(u(3)-u(2))+(u(4)-u(1)) ))))/48.
211 :
212 0 : fint=p(r)
213 :
214 0 : end function fint
215 : ! **********************************************************************
216 0 : function fintp(u,r)
217 : dimension u(4)
218 :
219 : ! If P(R) = u(1) u(2) u(3) u(4)
220 : ! for R = -3 -1 1 3
221 : ! then a cubic fit to the derivative is:
222 : 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.
223 :
224 0 : fintp=pp(r)
225 :
226 0 : end function fintp
227 :
228 : ! **********************************************************************
229 0 : subroutine scatt(ih, jh, rion, uf, f, umesh, semesh, dscat, ntot, epa, ierr)
230 : use op_load, only: BRCKR
231 : integer, intent(inout) :: ierr
232 : real :: umesh(:), semesh(:) ! (nptot)
233 : real :: f(:,:,:) ! (nptot,4,4)
234 0 : dimension rion(28, 4, 4),uf(0:100),fscat(0:100),p(nptot),rr(28),ih(4),jh(4)
235 : integer :: i,j,k,n
236 : ! HH: always use meshtype q='m'
237 0 : ite3=2
238 0 : umin=umesh(1)
239 0 : CSCAT=EPA*2.37567E-8
240 0 : ierr = 0
241 :
242 0 : do i = 1, 4
243 0 : ft = real(exp10(ITE3*dble(ih(i))/40d0))
244 0 : do j = 1, 4
245 0 : fne = real(exp10(ITE3*dble(jh(j))/4d0))
246 0 : do k = 1, ntot
247 0 : p(k) = f(k, i, j)
248 : end do
249 0 : do m = 1, 28
250 0 : rr(m) = rion(m, i, j)
251 : end do
252 0 : call BRCKR(FT,FNE,RR,28,UF,100,FSCAT,ierr)
253 0 : if (ierr /= 0) return
254 0 : do n = 0, 100
255 0 : u = uf(n)
256 0 : fscat(n) = cscat*(fscat(n)-1)
257 : end do
258 0 : do n = 2, ntot-1
259 0 : u=umesh(n)
260 0 : se=semesh(n)
261 0 : m=(u-umin)/dscat
262 0 : ua=umin+dscat*m
263 0 : ub=ua+dscat
264 0 : p(n)=p(n)+((ub-u)*fscat(m)+(u-ua)*fscat(m+1))/(dscat*se)
265 : end do
266 0 : u=umesh(ntot)
267 0 : se=semesh(ntot)
268 0 : p(ntot)=p(ntot)+fscat(100)/se
269 0 : p(1)=p(1)+fscat(1)/(1.-.5*umin)
270 0 : do k=1,ntot
271 0 : f(k,i,j)=p(k)
272 : end do
273 : end do
274 : end do
275 :
276 0 : end subroutine scatt
277 : ! **********************************************************************
278 0 : subroutine screen1(ih,jh,rion,umesh,ntot,epa,f)
279 0 : use op_load, only: screen2
280 : real :: umesh(:) ! (nptot)
281 : real, pointer :: f(:,:,:) ! (nptot,4,4)
282 : dimension uf(0:100),fscat(0:100), ih(4), jh(4)
283 : real, target :: rion(28, 4, 4)
284 : integer :: i, j, k, m
285 0 : real, pointer :: p(:), rr(:)
286 :
287 0 : ite3=2
288 0 : umin=umesh(1)
289 0 : umax=umesh(ntot)
290 :
291 0 : do i = 1, 4
292 0 : ft = real(exp10(ITE3*dble(ih(i))/40d0))
293 0 : do j = 1, 4
294 0 : fne = real(exp10(ITE3*dble(jh(j))/4d0))
295 : ! do k=1,ntot
296 0 : p => f(1:ntot,i,j)
297 : ! end do
298 : ! do m=1,28
299 0 : rr => rion(1:28,i,j)
300 : ! end do
301 0 : call screen2(ft,fne,rr,epa,ntot,umin,umax,umesh,p)
302 : ! do k=1,ntot
303 : ! f(k,i,j)=p(k)
304 : ! end do
305 : end do
306 : end do
307 :
308 0 : end subroutine screen1
309 :
310 : end module op_common
|
