Line data Source code
1 : #undef _GAUSS_TABLE
2 : module quadrature_mod
3 : use shr_kind_mod, only: r8=>shr_kind_r8
4 :
5 : implicit none
6 : private
7 :
8 : type, public :: quadrature_t
9 : real (kind=r8), dimension(:), pointer :: points
10 : real (kind=r8), dimension(:), pointer :: weights
11 : end type quadrature_t
12 :
13 : public :: gausslobatto
14 : public :: test_gausslobatto
15 : public :: gauss
16 : public :: test_gauss
17 : public :: legendre
18 : public :: jacobi
19 : public :: quad_norm
20 :
21 : public :: trapezoid
22 : private :: trapN
23 : public :: simpsons
24 : public :: gaussian_int
25 :
26 : private :: gausslobatto_pts
27 : private :: gausslobatto_wts
28 : private :: gauss_pts
29 : private :: gauss_wts
30 : private :: jacobi_polynomials
31 : private :: jacobi_derivatives
32 :
33 :
34 : contains
35 :
36 : ! ==============================================================
37 : ! gauss:
38 : !
39 : ! Find the Gauss collocation points and the corresponding weights.
40 : !
41 : ! ==============================================================
42 :
43 0 : function gauss(npts) result(gs)
44 : integer, intent(in) :: npts
45 : type (quadrature_t) :: gs
46 :
47 0 : allocate(gs%points(npts))
48 0 : allocate(gs%weights(npts))
49 :
50 0 : gs%points=gauss_pts(npts)
51 0 : gs%weights=gauss_wts(npts,gs%points)
52 :
53 0 : end function gauss
54 :
55 : #if defined(_GAUSS_TABLE)
56 : function gauss_pts(npts) result(pts)
57 :
58 : integer, intent(in) :: npts
59 : real (kind=r8) :: pts(npts)
60 :
61 : pts(1) = -0.93246951420315202781_r8
62 : pts(2) = -0.66120938646626451366_r8
63 : pts(3) = -0.23861918608319690863_r8
64 : pts(4) = -pts(3)
65 : pts(5) = -pts(2)
66 : pts(6) = -pts(1)
67 :
68 : end function gauss_pts
69 :
70 :
71 : function gauss_wts(npts,pts) result(wts)
72 :
73 : integer, intent(in) :: npts
74 : real (kind=r8) :: pts(npts)
75 : real (kind=r8) :: wts(npts)
76 :
77 : wts(1) = 0.17132449237917034504_r8
78 : wts(2) = 0.36076157304813860756_r8
79 : wts(3) = 0.46791393457269104738_r8
80 : wts(4) = wts(3)
81 : wts(5) = wts(2)
82 : wts(6) = wts(1)
83 :
84 : end function gauss_wts
85 : #else
86 :
87 : ! ==============================================================
88 : ! gauss_pts:
89 : !
90 : ! Compute the Gauss Collocation points
91 : ! for Jacobi Polynomials
92 : !
93 : ! ==============================================================
94 :
95 0 : function gauss_pts(np1) result(pts)
96 : use physconst, only: pi
97 :
98 : integer, intent(in) :: np1 ! Number of velocity grid points
99 : real (kind=r8) :: pts(np1)
100 :
101 : ! Local variables
102 :
103 : real (kind=r8) :: alpha,beta
104 0 : real (kind=r8) :: xjac(0:np1-1)
105 0 : real (kind=r8) :: jac(0:np1)
106 0 : real (kind=r8) :: djac(0:np1)
107 :
108 : integer prec ! number of mantissa bits
109 : real (kind=r8) eps ! machine epsilon
110 : real (kind=r8), parameter :: convthresh = 10 ! convergence threshold relative\
111 :
112 : ! to machine epsilon
113 : integer, parameter :: kstop = 30 ! max iterations for polynomial deflation
114 :
115 : real (kind=r8) :: poly
116 : real (kind=r8) :: pder
117 : real (kind=r8) :: recsum,thresh
118 : real (kind=r8) :: dth
119 :
120 : real (kind=r8) :: x
121 : real (kind=r8) :: delx
122 : real (kind=r8) :: c0,c1,c2,c10
123 :
124 : integer i,j,k
125 : integer n, nh
126 :
127 0 : n = np1 - 1
128 0 : c0 = 0.0_r8
129 0 : c1 = 1.0_r8
130 0 : c2 = 2.0_r8
131 0 : c10 = 10.0_r8
132 0 : alpha = c0
133 0 : beta = c0
134 :
135 : ! =========================================================
136 : ! compute machine precision and set the convergence
137 : ! threshold thresh to 10 times that level
138 : ! =========================================================
139 :
140 0 : prec = precision(c10)
141 0 : eps = c10**(-prec)
142 0 : thresh = convthresh*eps
143 :
144 : ! ============================================================
145 : ! Compute first half of the roots by "polynomial deflation".
146 : ! ============================================================
147 :
148 0 : dth = PI/(2*n+2)
149 :
150 0 : nh = (n+1)/2
151 :
152 0 : do j=0,nh-1
153 0 : x=COS((c2*j+1)*dth) ! first guess at root
154 0 : k=0
155 0 : delx=c1
156 0 : do while(k<kstop .and. ABS(delx) > thresh)
157 0 : call jacobi(n+1,x,alpha,beta,jac(0:n+1),djac(0:n+1))
158 0 : poly = jac(n+1)
159 0 : pder = djac(n+1)
160 0 : recsum=c0
161 0 : do i=0,j-1
162 0 : recsum = recsum + c1/(x-xjac(i))
163 : end do
164 0 : delx = -poly/(pder-recsum*poly)
165 0 : x = x + delx
166 0 : k = k + 1
167 : end do
168 :
169 0 : xjac(j)=x
170 :
171 : end do
172 :
173 : ! ================================================
174 : ! compute the second half of the roots by symmetry
175 : ! ================================================
176 :
177 0 : do j=0,nh
178 0 : xjac(n-j) = -xjac(j)
179 : end do
180 :
181 0 : if (MODULO(n,2)==0) xjac(nh)=c0
182 :
183 : ! ====================================================
184 : ! Reverse the sign of everything so that indexing
185 : ! increases with position
186 : ! ====================================================
187 :
188 0 : do j=0,n
189 0 : pts(j+1) = -xjac(j)
190 : end do
191 :
192 0 : end function gauss_pts
193 :
194 : ! ================================================
195 : ! gauss_wts:
196 : !
197 : ! Gauss Legendre Weights
198 : ! ================================================
199 :
200 0 : function gauss_wts(np1, gpts) result(wts)
201 :
202 : integer, intent(in) :: np1
203 : real (kind=r8), intent(in) :: gpts(np1) ! Gauss-Legendre points
204 : real (kind=r8) :: wts(np1) ! Gauss-Legendre weights
205 :
206 : ! Local variables
207 :
208 : real (kind=r8) :: c0,c1,c2
209 : real (kind=r8) :: alpha
210 : real (kind=r8) :: beta
211 0 : real (kind=r8) :: djac(np1)
212 : integer i,n
213 :
214 0 : c0 = 0.0_r8
215 0 : c1 = 1.0_r8
216 0 : c2 = 2.0_r8
217 :
218 0 : alpha = c0
219 0 : beta = c0
220 0 : n = np1-1
221 :
222 0 : djac=jacobi_derivatives(np1,alpha,beta,np1,gpts)
223 :
224 0 : do i=1,np1
225 0 : wts(i)=c2/((c1-gpts(i)**2)*djac(i)*djac(i))
226 : end do
227 :
228 0 : end function gauss_wts
229 :
230 : #endif
231 :
232 : ! ==============================================================
233 : ! test_gauss:
234 : !
235 : ! Unit Tester for Gaussian Points, Weights
236 : ! ==============================================================
237 :
238 0 : subroutine test_gauss(npts)
239 :
240 : integer, intent(in) :: npts
241 : type (quadrature_t) :: gs
242 :
243 : integer i
244 : real (kind=r8) :: gssum
245 0 : gs=gauss(npts)
246 :
247 0 : print *
248 0 : print *,"============================================"
249 0 : print *," Testing Gaussian Quadrature..."
250 0 : print *
251 0 : print *," points weights"
252 0 : print *,"============================================"
253 0 : do i=1,npts
254 0 : print *,i,gs%points(i),gs%weights(i)
255 : end do
256 0 : print *,"============================================"
257 0 : gssum=SUM(gs%weights(:))
258 0 : print *,"sum of Gaussian weights=",gssum
259 0 : print *,"============================================"
260 :
261 0 : deallocate(gs%points)
262 0 : deallocate(gs%weights)
263 :
264 0 : end subroutine test_gauss
265 :
266 : ! ==============================================================
267 : ! gausslobatto:
268 : !
269 : ! Find the Gauss-Lobatto Legendre collocation points xgl(i) and the
270 : ! corresponding weights.
271 : !
272 : ! ==============================================================
273 :
274 2622024 : function gausslobatto(npts) result(gll)
275 :
276 : integer, intent(in) :: npts
277 : type (quadrature_t) :: gll
278 :
279 7866072 : allocate(gll%points(npts))
280 5244048 : allocate(gll%weights(npts))
281 :
282 13110120 : gll%points=gausslobatto_pts(npts)
283 15732144 : gll%weights=gausslobatto_wts(npts,gll%points)
284 :
285 2622024 : end function gausslobatto
286 :
287 : ! ==============================================================
288 : ! gausslobatto_pts:
289 : !
290 : ! Compute the Gauss-Lobatto Collocation points
291 : ! for Jacobi Polynomials
292 : !
293 : ! ==============================================================
294 :
295 5244048 : function gausslobatto_pts(np1) result(pts)
296 : use physconst, only: pi
297 :
298 : integer, intent(in) :: np1 ! Number of velocity grid points
299 : real (kind=r8) :: pts(np1)
300 :
301 : ! Local variables
302 :
303 : real (kind=r8) :: alpha,beta
304 5244048 : real (kind=r8) :: xjac(0:np1-1)
305 5244048 : real (kind=r8) :: jac(0:np1)
306 5244048 : real (kind=r8) :: jacm1(0:np1)
307 2622024 : real (kind=r8) :: djac(0:np1)
308 :
309 : integer prec ! number of mantissa bits
310 : real (kind=r8) eps ! machine epsilon
311 : real (kind=r8), parameter :: convthresh = 10 ! convergence threshold relative
312 : ! to machine epsilon
313 : integer, parameter :: kstop = 30 ! max iterations for polynomial deflation
314 :
315 : real (kind=r8) :: a,b,det
316 : real (kind=r8) :: poly
317 : real (kind=r8) :: pder
318 : real (kind=r8) :: recsum,thresh
319 : real (kind=r8) :: dth,cd,sd,cs,ss,cstmp
320 :
321 : real (kind=r8) :: x
322 : real (kind=r8) :: delx
323 : real (kind=r8) :: c0,c1,c2,c10
324 :
325 : integer i,j,k
326 : integer n, nh
327 :
328 2622024 : n = np1 - 1
329 2622024 : c0 = 0.0_r8
330 2622024 : c1 = 1.0_r8
331 2622024 : c2 = 2.0_r8
332 2622024 : c10 = 10.0_r8
333 :
334 2622024 : alpha = c0
335 2622024 : beta = c0
336 :
337 : ! =========================================================
338 : ! compute machine precision and set the convergence
339 : ! threshold thresh to 10 times that level
340 : ! =========================================================
341 :
342 2622024 : prec = PRECISION(c10)
343 2622024 : eps = c10**(-prec)
344 2622024 : thresh = convthresh*eps
345 :
346 : ! =====================================================
347 : ! initialize the end points
348 : ! =====================================================
349 :
350 2622024 : xjac(0) = c1
351 2622024 : xjac(n) = -c1
352 :
353 : ! ============================================================
354 : ! Compute first half of the roots by "polynomial deflation".
355 : ! ============================================================
356 :
357 : ! ============================================================
358 : ! compute the parameters in the polynomial whose
359 : ! roots are desired...
360 : ! ============================================================
361 :
362 2622024 : call jacobi(n+1, c1,alpha,beta,jac(0:n+1),djac(0:n+1))
363 2622024 : call jacobi(n+1,-c1,alpha,beta,jacm1(0:n+1),djac(0:n+1))
364 :
365 2622024 : det = jac(n )*jacm1(n-1)-jacm1(n )*jac(n-1)
366 2622024 : a = -(jac(n+1)*jacm1(n-1)-jacm1(n+1)*jac(n-1))/det
367 2622024 : b = -(jac(n )*jacm1(n+1)-jacm1(n )*jac(n+1))/det
368 :
369 2622024 : dth = PI/(2*n+1)
370 2622024 : cd = COS(c2*dth)
371 2622024 : sd = SIN(c2*dth)
372 2622024 : cs = COS(dth)
373 2622024 : ss = SIN(dth)
374 :
375 2622024 : nh = (n+1)/2
376 :
377 5244048 : do j=1,nh-1
378 2622024 : x=cs ! first guess at root
379 2622024 : k=0
380 2622024 : delx=c1
381 20976192 : do while(k<kstop .and. ABS(delx) > thresh)
382 18354168 : call jacobi(n+1,x,alpha,beta,jac(0:n+1),djac(0:n+1))
383 18354168 : poly = jac(n+1)+a* jac(n)+b* jac(n-1)
384 18354168 : pder = djac(n+1)+a*djac(n)+b*djac(n-1)
385 18354168 : recsum=c0
386 36708336 : do i=0,j-1
387 36708336 : recsum = recsum + c1/(x-xjac(i))
388 : end do
389 18354168 : delx = -poly/(pder-recsum*poly)
390 18354168 : x = x + delx
391 18354168 : k = k + 1
392 : end do
393 :
394 2622024 : xjac(j)=x
395 :
396 : ! =====================================================
397 : ! compute the guesses for the roots
398 : ! for the next points, i.e :
399 : !
400 : ! ss = sn(theta) => sin(theta+2*dth)
401 : ! cs = cs(theta) => cs(theta+2*dth)
402 : ! =====================================================
403 :
404 2622024 : cstmp=cs*cd-ss*sd
405 2622024 : ss=cs*sd+ss*cd
406 5244048 : cs=cstmp
407 : end do
408 :
409 : ! ================================================
410 : ! compute the second half of the roots by symmetry
411 : ! ================================================
412 :
413 7866072 : do j=1,nh
414 7866072 : xjac(n-j) = -xjac(j)
415 : end do
416 :
417 2622024 : if (MODULO(n,2)==0) xjac(nh)=c0
418 :
419 : ! ====================================================
420 : ! Reverse the sign of everything so that indexing
421 : ! increases with position
422 : ! ====================================================
423 :
424 13110120 : do j=0,n
425 13110120 : pts(j+1) = -xjac(j)
426 : end do
427 :
428 2622024 : end function gausslobatto_pts
429 :
430 : ! ================================================
431 : ! Gauss Lobatto Legendre Weights
432 : ! ================================================
433 :
434 5244048 : function gausslobatto_wts(np1, glpts) result(wts)
435 :
436 : integer, intent(in) :: np1
437 : real (kind=r8), intent(in) :: glpts(np1)
438 : real (kind=r8) :: wts(np1)
439 :
440 : ! Local variables
441 :
442 : real (kind=r8) :: c0,c2
443 : real (kind=r8) :: alpha
444 : real (kind=r8) :: beta
445 2622024 : real (kind=r8) :: jac(np1)
446 : integer i,n
447 :
448 2622024 : c0 = 0.0_r8
449 2622024 : c2 = 2.0_r8
450 2622024 : alpha = c0
451 2622024 : beta = c0
452 2622024 : n = np1-1
453 :
454 2622024 : jac=jacobi_polynomials(n,alpha,beta,np1,glpts)
455 :
456 13110120 : do i=1,np1
457 13110120 : wts(i)=c2/(n*(n+1)*jac(i)*jac(i))
458 : end do
459 :
460 2622024 : end function gausslobatto_wts
461 :
462 : ! ==============================================================
463 : ! test_gausslobatto:
464 : !
465 : ! Unit Tester for Gaussian Lobatto Quadrature...
466 : ! ==============================================================
467 :
468 0 : subroutine test_gausslobatto(npts)
469 : integer, intent(in) :: npts
470 : type (quadrature_t) :: gll
471 :
472 : integer i
473 : real (kind=r8) :: gllsum
474 0 : gll=gausslobatto(npts)
475 :
476 0 : print *
477 0 : print *,"============================================"
478 0 : print *," Testing Gauss-Lobatto Quadrature..."
479 0 : print *
480 0 : print *," points weights"
481 0 : print *,"============================================"
482 0 : do i=1,npts
483 0 : print *,i,gll%points(i),gll%weights(i)
484 : end do
485 0 : print *,"============================================"
486 0 : gllsum=SUM(gll%weights(:))
487 0 : print *,"sum of Gauss-Lobatto weights=",gllsum
488 0 : print *,"============================================"
489 :
490 0 : deallocate(gll%points)
491 0 : deallocate(gll%weights)
492 :
493 0 : end subroutine test_gausslobatto
494 :
495 : ! ================================================
496 : !
497 : ! subroutine jacobi:
498 : !
499 : ! Computes the Jacobi Polynomials (jac) and their
500 : ! first derivatives up to and including degree n
501 : ! at point x on the interval (-1,1).
502 : !
503 : ! See for example the recurrence relations
504 : ! in equation 2.5.4 (page 70) in
505 : !
506 : ! "Spectral Methods in Fluid Dynamics",
507 : ! by C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A.Zang
508 : ! Springer-Verlag, 1988.
509 : ! ================================================
510 :
511 23598216 : subroutine jacobi(n, x, alpha, beta, jac, djac)
512 :
513 : integer, intent(in) :: n
514 : real (kind=r8), intent(in) :: x
515 : real (kind=r8), intent(in) :: alpha
516 : real (kind=r8), intent(in) :: beta
517 : real (kind=r8) :: jac(0:n)
518 : real (kind=r8) :: djac(0:n)
519 :
520 : ! Local variables
521 :
522 : real (kind=r8) :: a1k
523 : real (kind=r8) :: a2k
524 : real (kind=r8) :: a3k
525 : real (kind=r8) :: da2kdx
526 :
527 : real (kind=r8) :: c2,c1,c0
528 :
529 : integer :: k
530 :
531 23598216 : c0 = 0.0_r8
532 23598216 : c1 = 1.0_r8
533 23598216 : c2 = 2.0_r8
534 :
535 23598216 : jac(0)=c1
536 23598216 : jac(1)=(c1 + alpha)*x
537 :
538 23598216 : djac(0)=c0
539 23598216 : djac(1)=(c1 + alpha)
540 :
541 94392864 : do k=1,n-1
542 70794648 : a1k = c2*( k + c1 )*( k + alpha + beta + c1 )*( c2*k + alpha + beta )
543 70794648 : da2kdx = ( c2*( k + c1 ) + alpha + beta )*( c2*k + alpha + beta + c1 )*( c2*k + alpha + beta )
544 70794648 : a2k = ( c2*k + alpha + beta + c1 )*( alpha*alpha - beta*beta ) + x*da2kdx
545 70794648 : a3k = c2*(k + alpha)*( k + beta )*( c2*k + alpha + beta + c2 )
546 70794648 : jac(k+1) = ( a2k*jac(k)-a3k*jac(k-1) )/a1k
547 94392864 : djac(k+1)= ( a2k*djac(k) + da2kdx*jac(k) - a3k*djac(k-1) )/a1k
548 : end do
549 :
550 23598216 : end subroutine jacobi
551 :
552 :
553 : ! ==========================================================
554 : ! This routine computes the Nth order Jacobi Polynomials
555 : ! (jac) for a vector of positions x on the interval (-1,1),
556 : ! of length npoints.
557 : !
558 : ! See for example the recurrence relations
559 : ! in equation 2.5.4 (page 70) in
560 : !
561 : ! "Spectral Methods in Fluid Dynamics",
562 : ! by C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A.Zang
563 : ! Springer-Verlag, 1988.
564 : !
565 : ! ===========================================================
566 :
567 2622024 : function jacobi_polynomials(n, alpha, beta, npoints, x) result(jac)
568 :
569 : integer, intent(in) :: n ! order of the Jacobi Polynomial
570 : real (kind=r8) :: alpha
571 : real (kind=r8) :: beta
572 : integer, intent(in) :: npoints
573 : real (kind=r8) :: x(npoints)
574 : real (kind=r8) :: jac(npoints)
575 :
576 : ! Local variables
577 :
578 : real (kind=r8) :: a1k
579 : real (kind=r8) :: a2k
580 : real (kind=r8) :: a3k
581 : real (kind=r8) :: da2kdx
582 :
583 : real (kind=r8) :: jacp1
584 : real (kind=r8) :: jacm1
585 : real (kind=r8) :: jac0
586 : real (kind=r8) :: xtmp
587 :
588 : real (kind=r8) :: c2,c1,c0
589 : integer j,k
590 :
591 : c0 = 0.0_r8
592 : c1 = 1.0_r8
593 : c2 = 2.0_r8
594 :
595 13110120 : do j = 1,npoints
596 :
597 10488096 : xtmp=x(j)
598 :
599 10488096 : jacm1=c1
600 10488096 : jac0 =(c1+alpha)*xtmp
601 :
602 31464288 : do k=1,n-1
603 20976192 : a1k=c2*(k+c1)*(k+alpha+beta+c1)*(c2*k+alpha+beta)
604 20976192 : da2kdx=(c2*k+alpha+beta+c2)*(c2*k+alpha+beta+c1)*(c2*k+alpha+beta)
605 20976192 : a2k=(c2*k+alpha+beta+c1)*(alpha*alpha-beta*beta) + xtmp*da2kdx
606 20976192 : a3k=c2*(k+alpha)*(k+beta)*(c2*k+alpha+beta+c2)
607 20976192 : jacp1=(a2k*jac0-a3k*jacm1)/a1k
608 20976192 : jacm1=jac0
609 31464288 : jac0 =jacp1
610 : end do
611 :
612 10488096 : if (n==0)jac0=jacm1
613 13110120 : jac(j)=jac0
614 : end do
615 :
616 2622024 : end function jacobi_polynomials
617 :
618 : ! ================================================
619 : ! This routine computes the first derivatives of Nth
620 : ! order Jacobi Polynomials (djac) for a vector of
621 : ! positions x on the interval (-1,1), of length npoints.
622 : !
623 : ! See for example the recurrence relations
624 : ! in equation 2.5.4 (page 70) in
625 : !
626 : ! "Spectral Methods in Fluid Dynamics",
627 : ! by C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A.Zang
628 : ! Springer-Verlag, 1988.
629 : !
630 : ! ================================================
631 :
632 0 : function jacobi_derivatives(n, alpha, beta, npoints, x) result(djac)
633 :
634 : integer , intent(in) :: n ! order of the Jacobi Polynomial
635 : real (kind=r8), intent(in) :: alpha
636 : real (kind=r8), intent(in) :: beta
637 : integer , intent(in) :: npoints
638 : real (kind=r8), intent(in) :: x(npoints)
639 :
640 : real (kind=r8) :: djac(npoints)
641 :
642 : ! Local variables
643 :
644 : ! Local variables
645 :
646 : real (kind=r8) :: a1k
647 : real (kind=r8) :: a2k
648 : real (kind=r8) :: a3k
649 : real (kind=r8) :: da2kdx
650 :
651 : real (kind=r8) :: jacp1
652 : real (kind=r8) :: jacm1
653 : real (kind=r8) :: jac0
654 : real (kind=r8) :: djacp1
655 : real (kind=r8) :: djacm1
656 : real (kind=r8) :: djac0
657 :
658 : real (kind=r8) :: xtmp
659 :
660 : real (kind=r8) :: c2,c1,c0
661 : integer j,k
662 :
663 : c0 = 0.0_r8
664 : c1 = 1.0_r8
665 : c2 = 2.0_r8
666 :
667 0 : do j = 1,npoints
668 :
669 0 : xtmp=x(j)
670 :
671 0 : jacm1=c1
672 0 : jac0 =(c1+alpha)*xtmp
673 :
674 0 : djacm1 = c0
675 0 : djac0 = (c1+alpha)
676 :
677 0 : do k=1,n-1
678 0 : a1k=c2*(k+c1)*(k+alpha+beta+c1)*(c2*k+alpha+beta)
679 0 : da2kdx=(c2*k+alpha+beta+c2)*(c2*k+alpha+beta+c1)*(c2*k+alpha+beta)
680 0 : a2k=(c2*k+alpha+beta+c1)*(alpha*alpha-beta*beta) + xtmp*da2kdx
681 0 : a3k=c2*(k+alpha)*(k+beta)*(c2*k+alpha+beta+c2)
682 :
683 0 : jacp1=(a2k*jac0-a3k*jacm1)/a1k
684 0 : djacp1=(a2k*djac0+da2kdx*jac0-a3k*djacm1)/a1k
685 :
686 0 : jacm1=jac0
687 0 : jac0=jacp1
688 :
689 0 : djacm1=djac0
690 0 : djac0=djacp1
691 :
692 : end do
693 :
694 0 : if (n==0)djac0=djacm1
695 0 : djac(j)=djac0
696 :
697 : end do
698 :
699 0 : end function jacobi_derivatives
700 :
701 : ! ===================================================
702 : !
703 : ! legendre:
704 : !
705 : ! Compute the legendre polynomials using
706 : ! the recurrence relationship.
707 : ! return leg(m+1) = P_N(x) for m=0..N
708 : ! p_3 = Legendre polynomial of degree N
709 : ! p_2 = Legendre polynomial of degree N-1 at x
710 : ! p_1 = Legendre polynomial of degree N-2 at x
711 : !
712 : ! ===================================================
713 :
714 20894784 : function legendre(x,N) result(leg)
715 :
716 : integer :: N
717 : real (kind=r8) :: x
718 : real (kind=r8) :: leg(N+1)
719 :
720 : real (kind=r8) :: p_1, p_2, p_3
721 : integer :: k
722 :
723 20894784 : p_3 = 1.0_r8
724 20894784 : leg(1)=p_3
725 20894784 : if (n.ne.0) then
726 20894784 : p_2 = p_3
727 20894784 : p_3 = x
728 20894784 : leg(2)=p_3
729 62684352 : do k = 2,N
730 41789568 : p_1 = p_2
731 41789568 : p_2 = p_3
732 41789568 : p_3 = ( (2*k-1)*x*p_2 - (k-1)*p_1 ) / k
733 62684352 : leg(k+1)=p_3
734 : end do
735 : end if
736 :
737 20894784 : end function legendre
738 :
739 :
740 : ! ===========================================
741 : ! quad_norm:
742 : !
743 : ! compute normalization constants
744 : ! for k=1,N order Legendre polynomials
745 : !
746 : ! e.g. gamma(k) in Canuto, page 58.
747 : !
748 : ! ===========================================
749 :
750 2608200 : function quad_norm(gquad,N) result(gamma)
751 : type (quadrature_t), intent(in) :: gquad
752 : integer , intent(in) :: N
753 :
754 : real (kind=r8) :: gamma(N)
755 :
756 : ! Local variables
757 2608200 : real (kind=r8) :: leg(N)
758 : integer :: i,k
759 :
760 13041000 : gamma(:)=0.0_r8
761 :
762 13041000 : do i=1,N
763 10432800 : leg=legendre(gquad%points(i),N-1)
764 54772200 : do k=1,N
765 52164000 : gamma(k)= gamma(k)+leg(k)*leg(k)*gquad%weights(i)
766 : end do
767 : end do
768 :
769 2608200 : end function quad_norm
770 :
771 : ! =======================
772 : ! TrapN:
773 : ! Numerical recipes
774 : ! =======================
775 :
776 0 : subroutine trapN(f,a,b,N,it,s)
777 : INTERFACE
778 : FUNCTION f(x) RESULT(f_x) ! Function to be integrated
779 : use shr_kind_mod, only: r8=>shr_kind_r8
780 : real(kind=r8), INTENT(IN) :: x
781 : real(kind=r8) :: f_x
782 : END FUNCTION f
783 : END INTERFACE
784 :
785 : real(kind=r8),intent(in) :: a,b
786 : integer, intent(in) :: N
787 : integer, intent(inout) :: it
788 : real(kind=r8), intent(inout) :: s
789 :
790 : real(kind=r8) :: ssum
791 : real(kind=r8) :: del
792 : real(kind=r8) :: rtnm
793 : real(kind=r8) :: x
794 :
795 : integer :: j
796 :
797 0 : if (N==1) then
798 0 : s = 0.5_r8*(b-a)*(f(a) + f(b))
799 0 : it =1
800 : else
801 0 : ssum = 0.0_r8
802 0 : rtnm =1.0_r8/it
803 0 : del = (b-a)*rtnm
804 0 : x=a+0.5_r8*del
805 0 : do j=1,it
806 0 : ssum = ssum + f(x)
807 0 : x=x+del
808 : end do
809 0 : s=0.5_r8*(s + del*ssum)
810 0 : it=2*it
811 : end if
812 :
813 0 : end subroutine trapN
814 :
815 : ! ==========================================
816 : ! Trapezoid Rule for integrating functions
817 : ! from a to b with residual error eps
818 : ! ==========================================
819 :
820 0 : function trapezoid(f,a,b,eps) result(Integral)
821 :
822 : integer, parameter :: Nmax = 25 ! At most 2^Nmax + 1 points in integral
823 :
824 : INTERFACE
825 : FUNCTION f(x) RESULT(f_x) ! Function to be integrated
826 : use shr_kind_mod, only: r8=>shr_kind_r8
827 : real(kind=r8), INTENT(IN) :: x
828 : real(kind=r8) :: f_x
829 : END FUNCTION f
830 : END INTERFACE
831 :
832 : real(kind=r8), intent(in) :: a,b ! The integral bounds
833 : real(kind=r8), intent(in) :: eps ! relative error bound for integral
834 : real(kind=r8) :: Integral ! the integral result (within eps)
835 : real(kind=r8) :: s ! Integral approximation
836 : real(kind=r8) :: sold ! previous integral approx
837 :
838 : integer :: N
839 : integer :: it
840 :
841 : ! ==============================================================
842 : ! Calculate I here using trapezoid rule using f and a DO loop...
843 : ! ==============================================================
844 :
845 0 : s = 1.0e30_r8
846 0 : sold = 0.0_r8
847 0 : N=1
848 0 : it=0
849 0 : do while(N<=Nmax .and. ABS(s-sold)>eps*ABS(sold))
850 0 : sold=s
851 0 : call trapN(f,a,b,N,it,s)
852 0 : N=N+1
853 : end do
854 :
855 0 : Integral = s
856 :
857 0 : end function trapezoid
858 :
859 : ! ==========================================
860 : ! Simpsons Rule for integrating functions
861 : ! from a to b with residual error eps
862 : ! ==========================================
863 :
864 0 : function simpsons(f,a,b,eps) result(Integral)
865 :
866 : integer, parameter :: Nmax = 25 ! At most 2^Nmax + 1 points in integral
867 :
868 : INTERFACE
869 : FUNCTION f(x) RESULT(f_x) ! Function to be integrated
870 : use shr_kind_mod, only: r8=>shr_kind_r8
871 : real(kind=r8), INTENT(IN) :: x
872 : real(kind=r8) :: f_x
873 : END FUNCTION f
874 : END INTERFACE
875 :
876 : real(kind=r8), intent(in) :: a,b ! The integral bounds
877 : real(kind=r8), intent(in) :: eps ! relative error bound for integral
878 : real(kind=r8) :: Integral ! the integral result (within eps)
879 : real(kind=r8) :: s ! Integral approximation
880 : real(kind=r8) :: os ! previous integral approx
881 : real(kind=r8) :: st ! Integral approximation
882 : real(kind=r8) :: ost ! previous integral approx
883 :
884 : integer :: N
885 : integer :: it
886 :
887 : ! ==============================================================
888 : ! Calculate I here using trapezoid rule using f and a DO loop...
889 : ! ==============================================================
890 :
891 0 : ost= 0.0_r8
892 0 : s = 1.0e30_r8
893 0 : os = 0.0_r8
894 :
895 0 : N=1
896 0 : it=0
897 0 : do while ((N<=Nmax .and. ABS(s-os)>eps*ABS(os) ) .or. N<=2)
898 0 : os = s
899 0 : call trapN(f,a,b,N,it,st)
900 0 : s=(4.0_r8*st-ost)/3.0_r8
901 0 : ost=st
902 0 : N=N+1
903 : end do
904 :
905 0 : Integral = s
906 :
907 0 : end function simpsons
908 :
909 :
910 : ! ==========================================
911 : ! gaussian_int:
912 : !
913 : ! Gaussian Quadrature Rule for integrating
914 : ! function f from a to b with gs weights and
915 : ! points with precomputed gaussian quadrature
916 : ! and weights.
917 : ! ==========================================
918 :
919 0 : function gaussian_int(f,a,b,gs) result(Integral)
920 :
921 : integer, parameter :: Nmax = 10 ! At most 2^Nmax + 1 points in integral
922 :
923 : INTERFACE
924 : FUNCTION f(x) RESULT(f_x) ! Function to be integrated
925 : use shr_kind_mod, only: r8=>shr_kind_r8
926 : real(kind=r8), INTENT(IN) :: x
927 : real(kind=r8) :: f_x
928 : END FUNCTION f
929 : END INTERFACE
930 :
931 : real(kind=r8), intent(in) :: a,b ! The integral bounds
932 : type(quadrature_t), intent(in) :: gs ! gaussian points/wts
933 : real(kind=r8) :: Integral ! the integral result (within eps)
934 :
935 : integer :: i
936 : real (kind=r8) :: s,x
937 : ! ==============================================================
938 : ! Calculate I = S f(x)dx here using gaussian quadrature
939 : ! ==============================================================
940 :
941 0 : s = 0.0_r8
942 0 : do i=1,SIZE(gs%points)
943 0 : x = 0.50_r8*((b-a)*gs%points(i) + (b+a))
944 0 : s = s + gs%weights(i)*f(x)
945 : end do
946 0 : Integral = s*(0.5_r8*(b-a))
947 :
948 0 : end function gaussian_int
949 :
950 0 : end module quadrature_mod
951 :
952 :
953 :
954 :
955 :
|
