Line data Source code
1 : MODULE m_intgr
2 :
3 : !**********************************************************************
4 : ! intgr[0-3]:
5 : ! Integrators of a function y(jri) on a logarithmic mesh with jri
6 : ! mesh points. The output is either a scalar [z] or a field [z(jri)].
7 : ! Either the first meshpoint [r0,r1] and the increment [h] or the
8 : ! array of meshpoints [rmsh(jri)] is supplied:
9 : !
10 : ! intgz0 & intgz1 :
11 : ! (in)definite integrators on linear mesh. tail=.true. will include
12 : ! a tail correction assuming that the function is a simple
13 : ! decaying exponential between the first mesh point and infinity.
14 : ! y contains the nmz function values tabulated at a spacing of h.
15 : !
16 : ! integrator: ---- input ---- output
17 : ! intgr0: definite y r0 h jri | z
18 : ! intgr1: indefinite y r1 h jri | z(jri)
19 : ! intgr2: indefinite y rmsh(jri) h jri | z(jri)
20 : ! intgr3: definite y rmsh(jri) h jri | z
21 : ! intgz0: definite y tail n nmz ! z
22 : ! intgz1: indefinite y tail n nmz ! z(nmz)
23 : !
24 : ! m. weinert
25 : !**********************************************************************
26 :
27 : IMPLICIT NONE
28 :
29 : !INTRINSIC exp,log
30 : INTERFACE
31 : REAL FUNCTION ddot( n, f1, is1, f2, is2 )
32 : INTEGER, INTENT (IN) :: n, is1, is2
33 : REAL, INTENT (IN) :: f1(n), f2(n)
34 : END FUNCTION
35 : END INTERFACE
36 :
37 : ! interface intgz1
38 : ! module procedure intgz1Real, intgz1Complex
39 : ! end interface
40 :
41 : interface intgz1Reverse
42 : module procedure intgz1RealReverse, intgz1ComplexReverse
43 : end interface
44 :
45 : INTEGER, PARAMETER, PRIVATE :: nr = 7 , nr1 = 6
46 : REAL, PARAMETER, PRIVATE :: h0 = 140.
47 :
48 : ! lagrangian integration coefficients (simpson 7 point rule: error h**9)
49 :
50 : INTEGER, DIMENSION(7), PARAMETER, PRIVATE :: ih = (/41,216,27,272,27,216,41/)
51 :
52 : REAL, DIMENSION(7,5), PARAMETER, PRIVATE :: a = RESHAPE( &
53 : (/19087.,65112.,-46461., 37504.,-20211., 6312.,-863., &
54 : -863.,25128., 46989.,-16256., 7299.,-2088., 271., &
55 : 271.,-2760., 30819., 37504., -6771., 1608.,-191., &
56 : -191., 1608., -6771., 37504., 30819.,-2760., 271., &
57 : 271.,-2088., 7299.,-16256., 46989.,25128.,-863./),(/7,5/))
58 :
59 :
60 : CONTAINS
61 :
62 : !**********************************************************************
63 682071 : SUBROUTINE intgr0( y, r0, h, jri, z )
64 : !**********************************************************************
65 : ! ..
66 : ! .. Arguments ..
67 : INTEGER, INTENT (IN) :: jri
68 : REAL, INTENT (IN) :: h, r0
69 : REAL, INTENT (IN) :: y(jri)
70 : REAL, INTENT (OUT) :: z
71 : ! ..
72 : ! .. Locals ..
73 : INTEGER :: m, n0, nsteps
74 : REAL :: dr, r(7)
75 : INTEGER :: i, j
76 : REAL :: alpha, yr(7)
77 :
78 : !
79 : !---> integral from 0 to r1 approximated by leading term in power
80 : !---> series expansion of y(r)
81 : !
82 682071 : z = 0.0
83 682071 : IF (y(1)*y(2).GT.0.0) THEN
84 601490 : alpha = 1.0 + log(y(2)/y(1))/h
85 601490 : IF (alpha.GT.0.0) z = r0*y(1)/alpha
86 : ENDIF
87 : !
88 : !---> determine steps and starting point for simpson
89 : !
90 682071 : nsteps = (jri-1)/nr1
91 682071 : n0 = jri - nr1*nsteps
92 682071 : dr = exp(h)
93 682071 : r(1) = r0
94 4774497 : DO i = 2,7
95 4774497 : r(i) = dr*r(i-1)
96 : ENDDO
97 : !
98 : !---> lagrange integration for points 1<j<n0, error: h**9
99 : !
100 682071 : IF (n0.GT.1) THEN
101 2332912 : DO i = 1,7
102 2332912 : yr(i) = r(i)*y(i)
103 : ENDDO
104 1068177 : DO j = 1,n0 - 1
105 1068177 : z = z + h*ddot(7,a(1,j),1,yr,1)/60480.
106 : ENDDO
107 : ENDIF
108 682071 : r(1) = r(n0)
109 : !
110 : !---> simpson integration
111 : !
112 84442317 : DO m = 1,nsteps
113 586321722 : DO i = 2,nr
114 586321722 : r(i) = dr*r(i-1)
115 : ENDDO
116 670081968 : DO i = 1,nr
117 670081968 : yr(i) = h*ih(i)*r(i)/h0
118 : ENDDO
119 83760246 : z = z + ddot(nr,yr,1,y(n0),1)
120 83760246 : n0 = n0 + nr1
121 84442317 : r(1) = r(nr)
122 : ENDDO
123 :
124 682071 : RETURN
125 : END SUBROUTINE intgr0
126 :
127 : !**********************************************************************
128 48049 : SUBROUTINE intgr1( y, r1, h, jri, z )
129 : !**********************************************************************
130 : ! ..
131 : ! .. Arguments ..
132 : INTEGER, INTENT (IN) :: jri
133 : REAL, INTENT (IN) :: h,r1
134 : REAL, INTENT (IN) :: y(jri)
135 : REAL, INTENT (OUT):: z(jri)
136 : ! ..
137 : ! .. Locals ..
138 : REAL :: dr, rr, r(7)
139 : INTEGER :: i, j
140 : REAL :: alpha, yr(7)
141 : !
142 : !---> integral from 0 to r1 approximated by leading term in power
143 : !---> series expansion of y(r)
144 : !
145 48049 : z(1) = 0.0
146 48049 : IF (y(1)*y(2).GT.0.0) THEN
147 48049 : alpha = 1.0 + log(y(2)/y(1))/h
148 48049 : IF (alpha.GT.0.0) z(1) = r1*y(1)/alpha
149 : ENDIF
150 : !
151 : !---> lagrange integration for points 1<j<nr, error: h**9
152 : !
153 48049 : dr = exp(h)
154 48049 : rr = r1
155 384392 : DO i = 1,7
156 336343 : yr(i) = rr*y(i)
157 336343 : r(i) = rr
158 384392 : rr = dr*rr
159 : ENDDO
160 288294 : DO j = 1,nr - 2
161 288294 : z(j+1) = z(j) + h*ddot(7,a(1,j),1,yr,1)/60480.
162 : ENDDO
163 : !
164 : !---> simpson integration, j>nr-1
165 : !
166 384392 : DO i = 1,nr
167 384392 : r(i) = h*ih(i)*r(i)/h0
168 : ENDDO
169 43113362 : DO j = nr,jri
170 43065313 : z(j) = z(j-nr1) + ddot(nr,r,1,y(j-nr1),1)
171 344570553 : DO i = 1,7
172 344522504 : r(i) = dr*r(i)
173 : ENDDO
174 : ENDDO
175 :
176 48049 : RETURN
177 : END SUBROUTINE intgr1
178 :
179 : !**********************************************************************
180 31338 : SUBROUTINE intgr2( y, rmsh, h, jri, z )
181 : !**********************************************************************
182 : ! ..
183 : ! .. Arguments ..
184 : INTEGER, INTENT (IN) :: jri
185 : REAL, INTENT (IN) :: h
186 : REAL, INTENT (IN) :: rmsh(jri), y(jri)
187 : REAL, INTENT (OUT):: z(jri)
188 : ! ..
189 : ! .. Locals ..
190 : REAL :: dr, r(7)
191 : INTEGER :: i, j
192 : REAL :: alpha, yr(7)
193 : !
194 : !---> integral from 0 to r1 approximated by leading term in power
195 : !---> series expansion of y(r)
196 : !
197 31338 : z(1) = 0.0
198 31338 : IF (y(1)*y(2).GT.0.0) THEN
199 24348 : alpha = 1.0 + log(y(2)/y(1))/h
200 24348 : IF (alpha.GT.0.0) z(1) = rmsh(1)*y(1)/alpha
201 : ENDIF
202 : !
203 : !---> lagrange integration for points 1<j<nr, error: h**9
204 : !
205 31338 : dr = exp(h)
206 250704 : DO i = 1,7
207 219366 : r(i) = rmsh(i)
208 250704 : yr(i) = rmsh(i)*y(i)
209 : ENDDO
210 188028 : DO j = 1,nr - 2
211 188028 : z(j+1) = z(j) + h*ddot(7,a(1,j),1,yr,1)/60480.
212 : ENDDO
213 : !
214 : !---> simpson integration, j>nr-1
215 : !
216 250704 : DO i = 1,nr
217 250704 : r(i) = h*ih(i)*r(i)/h0
218 : ENDDO
219 22213732 : DO j = nr,jri
220 22182394 : z(j) = z(j-nr1) + ddot(nr,r,1,y(j-nr1),1)
221 177490490 : DO i = 1,7
222 177459152 : r(i) = dr*r(i)
223 : ENDDO
224 : ENDDO
225 :
226 31338 : RETURN
227 : END SUBROUTINE intgr2
228 : !**********************************************************************
229 0 : SUBROUTINE intgr3_modern( y, r, h, jri, z )
230 : !**********************************************************************
231 : ! ..
232 : ! .. Arguments ..
233 : INTEGER, INTENT (IN) :: jri
234 : REAL, INTENT (IN) :: h
235 : REAL, INTENT (IN) :: r(:)
236 : REAL, INTENT (IN) :: y(:)
237 : REAL, INTENT (OUT):: z
238 : ! ..
239 : ! .. Locals ..
240 : INTEGER :: m, n0, nsteps
241 : REAL :: tiny, h1, z1, ih1(nr)
242 : INTEGER :: i, j
243 : REAL :: alpha
244 : !
245 : !---> integral from 0 to r1 approximated by leading term in power
246 : !---> series expansion of y(r)
247 : !
248 : ! DO i=1,jri
249 : ! IF (abs(y(i)).LT.tiny) y(i) = tiny
250 : ! ENDDO
251 : !
252 0 : z = 0.0
253 0 : IF (y(1)*y(2).GT.0.0) THEN
254 0 : alpha = 1.0 + log(y(2)/y(1))/h
255 0 : IF (alpha.GT.0.0) z = r(1)*y(1)/alpha
256 : ENDIF
257 : !
258 : !---> determine steps and starting point for simpson
259 : !
260 0 : nsteps = (jri-1)/nr1
261 0 : n0 = jri - nr1*nsteps
262 : !
263 : !---> lagrange integration for points 1<j<n0, error: h**9
264 : !
265 0 : z1 = 0.
266 0 : DO j = 1,n0 - 1
267 0 : z1 = z1 + dot_product(a(:,j), r(1:7)*y(1:7))
268 : ENDDO
269 0 : z = z + z1 * h / 60480.
270 : !
271 : !---> simpson integration
272 : !
273 0 : h1 = h / h0
274 0 : DO i = 1,nr
275 0 : ih1(i) = h1 * ih(i)
276 : ENDDO
277 0 : DO m = 1,nsteps
278 0 : z = z + dot_product(ih1*r(n0:n0+6),y(n0:n0+nr-1))
279 0 : n0 = n0 + nr1
280 : ENDDO
281 :
282 0 : RETURN
283 : END SUBROUTINE intgr3_modern
284 : !**********************************************************************
285 10194211 : SUBROUTINE intgr3( y, r, h, jri, z )
286 : !**********************************************************************
287 : ! ..
288 : ! .. Arguments ..
289 : INTEGER, INTENT (IN) :: jri
290 : REAL, INTENT (IN) :: h
291 : REAL, INTENT (IN) :: r(jri)
292 : REAL, INTENT (IN) :: y(jri)
293 : REAL, INTENT (OUT):: z
294 : ! ..
295 : ! .. Locals ..
296 : INTEGER :: m, n0, nsteps
297 : REAL :: tiny, h1, z1, ih1(nr)
298 : INTEGER :: i, j
299 : REAL :: alpha
300 : !
301 : !---> integral from 0 to r1 approximated by leading term in power
302 : !---> series expansion of y(r)
303 : !
304 : ! DO i=1,jri
305 : ! IF (abs(y(i)).LT.tiny) y(i) = tiny
306 : ! ENDDO
307 : !
308 10194211 : z = 0.0
309 10194211 : IF (y(1)*y(2).GT.0.0) THEN
310 8538929 : alpha = 1.0 + log(y(2)/y(1))/h
311 8538929 : IF (alpha.GT.0.0) z = r(1)*y(1)/alpha
312 : ENDIF
313 : !
314 : !---> determine steps and starting point for simpson
315 : !
316 10194211 : nsteps = (jri-1)/nr1
317 10194211 : n0 = jri - nr1*nsteps
318 : !
319 : !---> lagrange integration for points 1<j<n0, error: h**9
320 : !
321 10194211 : z1 = 0.
322 19511137 : DO j = 1,n0 - 1
323 84729619 : z1 = z1 + dot_product(a(:,j), r(1:7)*y(1:7))
324 : ENDDO
325 10194211 : z = z + z1 * h / 60480.
326 : !
327 : !---> simpson integration
328 : !
329 10194211 : h1 = h / h0
330 81553688 : DO i = 1,nr
331 81553688 : ih1(i) = h1 * ih(i)
332 : ENDDO
333 1257365243 : DO m = 1,nsteps
334 9977368256 : z = z + dot_product(ih1*r(n0:n0+6),y(n0:n0+nr-1))
335 1257365243 : n0 = n0 + nr1
336 : ENDDO
337 :
338 10194211 : RETURN
339 : END SUBROUTINE intgr3
340 :
341 : !**********************************************************************
342 7001169 : SUBROUTINE intgz0( y, h, nmz, z, tail )
343 : !**********************************************************************
344 : ! ..
345 : ! .. Arguments ..
346 : INTEGER, INTENT (IN) :: nmz
347 : REAL, INTENT (IN) :: h
348 : REAL, INTENT (IN) :: y(nmz)
349 : LOGICAL, INTENT (IN) :: tail
350 : REAL, INTENT (OUT) :: z
351 : ! ..
352 : ! .. Locals ..
353 : REAL :: yl, ys
354 : INTEGER :: m, n0, nsteps
355 : INTEGER :: i, j
356 : REAL :: alpha, yr(7)
357 : !
358 : !---> integral from minus infinity to the first mesh point assuming
359 : !---> exponential decay in this region. this contribution is limited
360 : !---> to alpha>0.1 which corresponds to square of wavefunctions
361 : !---> of energy > 0.00125 a.u. below the vacuum level
362 : !
363 7001169 : z = 0.0
364 7001169 : IF (tail) THEN
365 6141124 : IF (y(1)*y(2).GT.0.0) THEN
366 4654546 : alpha = log(y(2)/y(1))/h
367 4654546 : IF (alpha.GT.0.1) z = y(1)/alpha
368 : ENDIF
369 : ENDIF
370 : !
371 : !---> determine steps and starting point for simpson
372 : !
373 7001169 : nsteps = (nmz-1)/nr1
374 7001169 : n0 = nmz - nr1*nsteps
375 : !
376 : !---> lagrange integration for points 1<j<n0, error: h**9
377 : !
378 7001169 : yl = 0.0
379 7001169 : IF (n0.GT.1) THEN
380 27322379 : DO j = 1, n0 - 1
381 27322379 : yl = yl + ddot(7,a(1,j),1,y,1)
382 : ENDDO
383 6774656 : yl = h*yl/60480.
384 : END IF
385 : !
386 : !---> simpson integration
387 : !
388 7001169 : ys = 0.0
389 7001169 : n0 = n0 - 1
390 117807787 : DO m = 1,nsteps
391 886452944 : DO i = 1,nr
392 886452944 : ys = ys + ih(i)*y(n0+i)
393 : ENDDO
394 117807787 : n0 = n0 + nr1
395 : ENDDO
396 7001169 : ys = h*ys/h0
397 7001169 : z = z + yl + ys
398 :
399 7001169 : RETURN
400 : END SUBROUTINE intgz0
401 :
402 : !**********************************************************************
403 37966 : SUBROUTINE intgz1( y, h, nmz, z, tail )
404 : !**********************************************************************
405 : ! ..
406 : ! .. Arguments ..
407 : INTEGER, INTENT (IN) :: nmz
408 : LOGICAL, INTENT (IN) :: tail
409 : REAL, INTENT (IN) :: h
410 : REAL, INTENT (IN) :: y(nmz)
411 : REAL, INTENT (OUT) :: z(nmz)
412 : ! ..
413 : ! .. Locals ..
414 : REAL :: yl, ys
415 : REAL, PARAMETER :: eps = 1.e-38
416 : INTEGER :: i, j
417 : REAL :: alpha, yr(7)
418 : !
419 : !---> integral from minus infinity to the first mesh point assuming
420 : !---> exponential decay in this region. this contribution is limited
421 : !---> to alpha>0.1 which corresponds to square of wavefunctions
422 : !---> of energy > 0.00125 a.u. below the vacuum level
423 : !
424 37966 : z(1) = 0.0
425 37966 : IF (tail) THEN
426 37966 : IF (abs(y(1)).GT.eps) THEN
427 2770 : IF (y(1)*y(2).GT.0.0) THEN
428 2766 : alpha = log(y(2)/y(1))/h
429 2766 : IF (alpha.GT.0.1) z(1) = y(1)/alpha
430 : ENDIF
431 : ENDIF
432 : ENDIF
433 : !
434 : !---> lagrange integration for points 1<j<n0, error: h**9
435 : !
436 227796 : DO j = 1,nr - 2
437 189830 : yl = 0
438 189830 : yl = yl + ddot(7,a(1,j),1,y,1)
439 227796 : z(j+1) = z(j) + h*yl/60480.
440 : ENDDO
441 : !
442 : !---> simpson integration
443 : !
444 3606770 : DO j = nr,nmz
445 : ys = 0.0
446 28550432 : DO i = 1,nr
447 28550432 : ys = ys + ih(i)*y(j-nr+i)
448 : ENDDO
449 3606770 : z(j) = z(j-nr1) + h*ys/h0
450 : ENDDO
451 :
452 37966 : RETURN
453 : END SUBROUTINE intgz1
454 :
455 :
456 :
457 0 : subroutine intgz1Complex( y, h, nmz, z, tail )
458 :
459 : integer, intent(in) :: nmz
460 : logical, intent(in) :: tail
461 : real, intent(in) :: h
462 : complex, intent(in) :: y(nmz)
463 : complex, intent(out) :: z(nmz)
464 :
465 0 : real :: zr(nmz), zi(nmz)
466 :
467 0 : call intgz1( real(y), h, nmz, zr, tail )
468 0 : call intgz1( aimag(y), h, nmz, zi, tail )
469 0 : z = cmplx( zr, zi )
470 :
471 0 : end subroutine intgz1Complex
472 :
473 :
474 :
475 0 : subroutine intgz1RealReverse( y, h, nmz, z, tail )
476 :
477 : integer, intent(in) :: nmz
478 : logical, intent(in) :: tail
479 : real, intent(in) :: h
480 : real, intent(in) :: y(nmz)
481 : real, intent(out) :: z(nmz)
482 :
483 0 : real :: y_reverse(nmz), z_reverse(nmz)
484 : integer :: i
485 :
486 0 : do i = 1, nmz
487 0 : y_reverse(i) = y(nmz+1-i)
488 : end do
489 0 : call intgz1( y_reverse, h, nmz, z_reverse, tail )
490 0 : do i = 1, nmz
491 0 : z(i) = z_reverse(nmz+1-i)
492 : end do
493 :
494 0 : end subroutine intgz1RealReverse
495 :
496 :
497 :
498 0 : subroutine intgz1ComplexReverse( y, h, nmz, z, tail )
499 :
500 : integer, intent(in) :: nmz
501 : logical, intent(in) :: tail
502 : real, intent(in) :: h
503 : complex, intent(in) :: y(nmz)
504 : complex, intent(out) :: z(nmz)
505 :
506 0 : complex :: y_reverse(nmz), z_reverse(nmz)
507 : integer :: i
508 :
509 0 : do i = 1, nmz
510 0 : y_reverse(i) = y(nmz+1-i)
511 : end do
512 0 : call intgz1Complex( y_reverse, h, nmz, z_reverse, tail )
513 0 : do i = 1, nmz
514 0 : z(i) = z_reverse(nmz+1-i)
515 : end do
516 :
517 0 : end subroutine intgz1ComplexReverse
518 :
519 324 : SUBROUTINE sfint(y,x,h,jri,z)
520 : ! Seperate spline interpolation integrator for source-free calculations,
521 : ! as the input here is notoriously unsmooth and intgr2 fails.
522 :
523 : INTEGER, INTENT (IN) :: jri
524 : REAL, INTENT (IN) :: h
525 : REAL, INTENT (IN) :: x(jri), y(jri)
526 : REAL, INTENT (OUT):: z(jri)
527 :
528 : REAL :: alpha
529 : INTEGER :: i
530 :
531 245592 : z = 0.0
532 324 : IF (y(1)*y(2).GT.0.0) THEN
533 148 : alpha = 1.0 + log(y(2)/y(1))/h
534 148 : IF (alpha.GT.0.0) z(1) = x(1)*y(1)/alpha
535 : ENDIF
536 :
537 324 : z(2)=z(1)+h*(x(1)*y(1)+x(2)*y(2))/2
538 :
539 244620 : DO i=3, jri-1
540 244620 : z(i)=z(i-1)-h*(x(i-2)*y(i-2)-13*x(i-1)*y(i-1)-13*x(i)*y(i)+x(i+1)*y(i+1))/24
541 : END DO
542 :
543 324 : z(i)=z(i-1)+h*(x(i-1)*y(i-1)+x(i)*y(i))/2
544 :
545 324 : END SUBROUTINE
546 :
547 : ! For juPhon:
548 0 : SUBROUTINE intgr3LinIntp(y,r,h,jri,z, i1)
549 :
550 : USE m_juDFT_stop, ONLY : juDFT_error
551 :
552 : INTEGER, INTENT (IN) :: jri
553 : INTEGER, INTENT (IN) :: i1
554 : REAL, INTENT (IN) :: h
555 : REAL, INTENT (IN) :: r(jri)
556 : REAL, INTENT (IN) :: y(jri)
557 : REAL, INTENT (OUT):: z
558 :
559 : INTEGER m, n0, nsteps, i, j
560 : REAL tiny, yr(nr), h1, z1, ih1(nr)
561 : real x1, x2, y1, y2, n
562 :
563 : z = 0.0
564 :
565 0 : z = 0.5 * r(i1) * y(i1)
566 0 : IF (i1.NE.1) CALL juDFT_error("intgr3LinIntp buggy for interpolations that do not start at mesh point 1.",calledby="intgr3LinIntp")
567 :
568 0 : nsteps = (jri-1)/nr1
569 0 : n0 = jri - nr1*nsteps
570 :
571 0 : IF (n0.GT.1) THEN
572 0 : DO i = 1, 7
573 0 : yr(i) = r(i)*y(i)
574 : END DO
575 : z1 = 0.
576 0 : DO j = 1, n0 - 1
577 0 : z1 = z1 + ddot(7,a(1,j),1,yr,1)
578 : END DO
579 0 : z = z + z1 * h / 60480.
580 : END IF
581 :
582 0 : h1 = h / h0
583 :
584 0 : DO i = 1,nr
585 0 : ih1(i) = h1 * ih(i)
586 : END DO
587 0 : DO m = 1,nsteps
588 0 : DO i = 1,nr
589 0 : yr(i) = ih1(i)*r(i+n0-1)
590 : END DO
591 0 : z = z + ddot(nr,yr,1,y(n0),1)
592 0 : n0 = n0 + nr1
593 : END DO
594 :
595 0 : RETURN
596 : END SUBROUTINE intgr3LinIntp
597 :
598 0 : SUBROUTINE intgr2LinIntp(y,rmsh,h,jri,z)
599 :
600 : INTEGER, INTENT (IN) :: jri
601 : REAL, INTENT (IN) :: h
602 : REAL, INTENT (IN) :: rmsh(jri),y(jri)
603 : REAL, INTENT (OUT):: z(jri)
604 :
605 : REAL :: dr, r(7), yr(nr)
606 : INTEGER :: i, j
607 :
608 0 : z = 0.0
609 :
610 0 : z(1) = 0.5 * rmsh(1) * y(1)
611 :
612 0 : dr = exp(h)
613 0 : DO i = 1,7
614 0 : r(i) = rmsh(i)
615 0 : yr(i) = rmsh(i)*y(i)
616 : END DO
617 :
618 0 : DO j = 1,nr - 2
619 0 : z(j+1) = z(j) + h*ddot(7,a(1,j),1,yr,1)/60480.
620 : END DO
621 :
622 0 : DO i = 1,nr
623 0 : r(i) = h*ih(i)*r(i)/h0
624 : END DO
625 :
626 0 : DO j = nr,jri
627 0 : z(j) = z(j-nr1) + ddot(nr,r,1,y(j-nr1),1)
628 0 : DO i = 1,7
629 0 : r(i) = dr*r(i)
630 : END DO
631 : END DO
632 :
633 0 : RETURN
634 : END SUBROUTINE intgr2LinIntp
635 :
636 : END MODULE m_intgr
|
