DISPATCH
radau_mod.f90
1 !===============================================================================
2 !===============================================================================
3 MODULE radau_mod
4  implicit none
5  private
6  type, public:: radau_t
7  contains
8  procedure, nopass:: init
9  end type
10  type(radau_t), public:: radau
11 CONTAINS
12 
13 !===============================================================================
14 !===============================================================================
15 SUBROUTINE init (k,xi,wi)
16 !
17 ! Routine to return weights and abscissa points for Gauss-Radau
18 ! integration in the interval [a; b]. The transformation from
19 ! the inherent interval of [-1; 1] is performed via c and d.
20 ! The Radau integration differs from a normal Gauss integration,
21 ! in that one of the end-points are included.
22 ! The number of integration points, k, can be chosen between
23 ! 1 and 10, and the order of the interpolating function is 2*k-1.
24 !
25 ! Update history:
26 ! ---------------
27 ! 29.08.2006 Coded/RT
28 ! 31.08.2006 Fixed single/double precision ambiguity and changed from
29 ! x=0 to x=1 being included as a node/RT
30 ! 15.01.2007 Ported from F77 to F90/RT
31 !
32 ! Known bugs: Hardcoded whether to use upper (coded) or lower end-point
33 ! in the integration.
34 !
35  implicit none
36  integer, parameter:: kmax = 10
37  integer:: k
38  real, dimension(k):: wi, xi
39  real:: a, b
40  real(kind=8), dimension(55):: xn = (/ &
41  -1.00000000000000d0, &! 1
42  -1.00000000000000d0, .33333333333333d0, &! 2
43  -1.00000000000000d0,-.28989794855663d0, .68989794855663d0,&! 3
44  -1.00000000000000d0,-.57531892352169d0, .18106627111853d0,&! 4
45  .82282408097459d0, &
46  -1.00000000000000d0,-.72048027131244d0,-.16718086473783d0,&! 5
47  .44631397272376d0, .88579160777096d0, &
48  -1.00000000000000d0,-.80292982840235d0,-.39092854670727d0,&! 6
49  .12405037950522d0, .60397316425279d0, .92038028589707d0,&
50  -1.00000000000000d0,-.85389134263949d0,-.53846772406011d0,&! 7
51  -.11734303754309d0, .32603061943769d0, .70384280066303d0,&
52  .94136714568042d0, &
53  -1.00000000000000d0,-.88747487892616d0,-.63951861652622d0,&! 8
54  -.29475056577366d0, .09430725266111d0, .46842035443083d0,&
55  .77064189367819d0, .95504122712258d0, &
56  -1.00000000000000d0,-.91073208942006d0,-.71126748591571d0,&! 9
57  -.42635048571114d0,-.09037336960685d0, .25613567083346d0,&
58  .57138304120873d0, .81735278420041d0, .96444016970528d0,&
59  -1.00000000000000d0,-.92748437423359d0,-.76384204242000d0,&! 10
60  -.52564603037008d0,-.23623446939059d0, .07605919783798d0,&
61  .38066484014473d0, .64776668767401d0, .85122522058160d0,&
62  .97117518070224d0/)
63  real(kind=8), dimension(55):: wn = (/ &
64  2.00000000000000d0, &! 1
65  0.50000000000000d0,1.50000000000005d0, &! 2
66  0.22222222222222d0,1.02497165237682d0,0.75280612540100d0,&! 3
67  0.12500000000000d0,0.65768863996010d0,0.77638693768637d0,&! 4
68  0.44092442235367d0, &
69  0.08000000000000d0,0.44620780216715d0,0.62365304595145d0,&! 5
70  0.56271203029887d0,0.28742712158247d0, &
71  0.05555555555556d0,0.31964075322051d0,0.48538718846895d0,&! 6
72  0.52092678318961d0,0.41690133431186d0,0.20158838525329d0,&
73  0.04081632653061d0,0.23922748922532d0,0.38094987364422d0,&! 7
74  0.44710982901453d0,0.42470377900601d0,0.31820423146727d0,&
75  0.14898847111224d0, &
76  0.03125000000000d0,0.18535815480298d0,0.30413062064680d0,&! 8
77  0.37651754538912d0,0.39157216745253d0,0.34701479563444d0,&
78  0.24964790132996d0,0.11450881474417d0, &
79  0.02469135802469d0,0.14765401904631d0,0.24718937820459d0,&! 9
80  0.31684377567046d0,0.34827300277294d0,0.33769396697589d0,&
81  0.28638669635730d0,0.20055329802460d0,0.09071450492309d0,&
82  0.02000000000000d0,0.12029667055749d0,0.20427013187899d0,&! 10
83  0.26819483784117d0,0.30585928772442d0,0.31358245722693d0,&
84  0.29061016483286d0,0.23919343171435d0,0.16437601273700d0,&
85  0.07361700548689d0/)
86  integer, dimension(kmax+1):: k0 = (/1,2,4,7,11,16,22,29,37,46,56/)
87  real c,d
88  integer i
89 
90  a = 0.0
91  b = 1.0
92  if (k.lt.1) then
93  print *,' WARNING(radaui): k=',k,'< 1 does not make sense!'
94  print *,' Setting k=1.'
95  k = 1
96  endif
97  if (k.gt.kmax) then
98  print *,' WARNING(radaui): k=',k,'>',kmax,' is not implemented.'
99  print *,' Setting k=',kmax,'.'
100  k = kmax
101  endif
102  c = .5*(b+a)
103  d = .5*(b-a)
104 !
105 ! Loop through the number of points in the integration
106  do i=1, k
107 !
108 ! Using this form instead of [c + d*xn(k0(k)+i-1)] makes
109 ! x=1 (or b) the predefined node instead of x=-1 (or a).
110  xi(i) = 1.0 - (c + d*xn(k0(k+1)-i))
111  wi(i) = d*wn(k0(k+1)-i)
112  enddo
113 
114  return
115 END SUBROUTINE
116 
117 END MODULE radau_mod