Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : !--------------------------------------------------------------------------------------------------!
9 : ! Copyright (c) 2008, Joost VandeVondele and Manuel Guidon !
10 : ! All rights reserved. !
11 : ! !
12 : ! Redistribution and use in source and binary forms, with or without !
13 : ! modification, are permitted provided that the following conditions are met: !
14 : ! * Redistributions of source code must retain the above copyright !
15 : ! notice, this list of conditions and the following disclaimer. !
16 : ! * Redistributions in binary form must reproduce the above copyright !
17 : ! notice, this list of conditions and the following disclaimer in the !
18 : ! documentation and/or other materials provided with the distribution. !
19 : ! !
20 : ! THIS SOFTWARE IS PROVIDED BY Joost VandeVondele and Manuel Guidon AS IS AND ANY !
21 : ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED !
22 : ! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE !
23 : ! DISCLAIMED. IN NO EVENT SHALL Joost VandeVondele or Manuel Guidon BE LIABLE FOR ANY !
24 : ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES !
25 : ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; !
26 : ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND !
27 : ! ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT !
28 : ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS !
29 : ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. !
30 : !--------------------------------------------------------------------------------------------------!
31 :
32 : ! **************************************************************************************************
33 : !> \brief ...
34 : !> \author Joost VandeVondele and Manuel Guidon
35 : !> \par History
36 : !> Nov 2008 Joost VandeVondele and Manuel Guidon
37 : !> Sep 2019 A.Bussy: Moved the file to common (together with gamma.F and t_c_g0.F)
38 : ! **************************************************************************************************
39 : MODULE t_sh_p_s_c
40 : USE kinds, ONLY: dp
41 : USE message_passing, ONLY: mp_comm_type
42 : #include "../base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 : PRIVATE
46 :
47 : PUBLIC :: trunc_CS_poly_n20, init, free_C0
48 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE :: C0
49 : INTEGER, PARAMETER :: degree = 15
50 : REAL(KIND=dp), PARAMETER :: target_error = 0.100000E-08
51 : INTEGER, PARAMETER :: nderiv_max = 21
52 : INTEGER, SAVE :: nderiv_init = -1
53 : INTEGER, SAVE :: patches = -1
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief ...
59 : !> \param RES ...
60 : !> \param R ...
61 : !> \param T ...
62 : !> \param NDERIV ...
63 : ! **************************************************************************************************
64 1160608 : SUBROUTINE trunc_CS_poly_n20(RES, R, T, NDERIV)
65 : REAL(KIND=dp), INTENT(OUT) :: RES(*)
66 : REAL(KIND=dp), INTENT(IN) :: R, T
67 : INTEGER, INTENT(IN) :: NDERIV
68 :
69 : REAL(KIND=dp), PARAMETER :: eps = 34.53877639491068526026987182_dp, n = 2.0_dp
70 :
71 : INTEGER :: i
72 : REAL(KIND=dp) :: lower, TG1, TG2, X1, X2
73 :
74 1160608 : IF (T <= 81.0_dp) THEN
75 851492 : X2 = 1/(1 + R)
76 851492 : X1 = 1/(1 + SQRT(T))
77 851492 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
78 851292 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
79 616408 : IF (X1 <= 0.550000000000000044E+00_dp) THEN
80 221696 : IF (X1 <= 0.325000000000000011E+00_dp) THEN
81 175856 : IF (X1 <= 0.212500000000000022E+00_dp) THEN
82 161232 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
83 0 : IF (X2 <= 0.625000000000000000E-01_dp) THEN
84 0 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.888888888888888751E+01_dp
85 0 : TG2 = (2*X2 - 0.625000000000000000E-01_dp)*0.160000000000000000E+02_dp
86 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 1))
87 : ELSE
88 0 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
89 0 : IF (X2 <= 0.937500000000000000E-01_dp) THEN
90 0 : IF (X1 <= 0.128124999999999989E+00_dp) THEN
91 0 : TG1 = (2*X1 - 0.228124999999999994E+00_dp)*0.355555555555555785E+02_dp
92 0 : TG2 = (2*X2 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
93 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 2))
94 : ELSE
95 0 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
96 0 : TG2 = (2*X2 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
97 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 3))
98 : END IF
99 : ELSE
100 0 : IF (X1 <= 0.128124999999999989E+00_dp) THEN
101 0 : IF (X2 <= 0.109375000000000000E+00_dp) THEN
102 0 : TG1 = (2*X1 - 0.228124999999999994E+00_dp)*0.355555555555555785E+02_dp
103 0 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
104 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 4))
105 : ELSE
106 0 : TG1 = (2*X1 - 0.228124999999999994E+00_dp)*0.355555555555555785E+02_dp
107 0 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
108 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 5))
109 : END IF
110 : ELSE
111 0 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
112 0 : TG2 = (2*X2 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
113 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 6))
114 : END IF
115 : END IF
116 : ELSE
117 0 : TG1 = (2*X1 - 0.368750000000000022E+00_dp)*0.177777777777777715E+02_dp
118 0 : TG2 = (2*X2 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
119 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 7))
120 : END IF
121 : END IF
122 : ELSE
123 161232 : IF (X2 <= 0.187500000000000000E+00_dp) THEN
124 22184 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
125 22184 : IF (X2 <= 0.156250000000000000E+00_dp) THEN
126 576 : IF (X1 <= 0.128124999999999989E+00_dp) THEN
127 576 : TG1 = (2*X1 - 0.228124999999999994E+00_dp)*0.355555555555555785E+02_dp
128 576 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
129 576 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 8))
130 : ELSE
131 0 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
132 0 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
133 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 9))
134 : END IF
135 : ELSE
136 21608 : IF (X1 <= 0.128124999999999989E+00_dp) THEN
137 20544 : IF (X1 <= 0.114062499999999997E+00_dp) THEN
138 16576 : TG1 = (2*X1 - 0.214062499999999989E+00_dp)*0.711111111111111569E+02_dp
139 16576 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
140 16576 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 10))
141 : ELSE
142 3968 : TG1 = (2*X1 - 0.242187500000000000E+00_dp)*0.711111111111111569E+02_dp
143 3968 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
144 3968 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 11))
145 : END IF
146 : ELSE
147 1064 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
148 1064 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
149 1064 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 12))
150 : END IF
151 : END IF
152 : ELSE
153 0 : IF (X2 <= 0.156250000000000000E+00_dp) THEN
154 0 : TG1 = (2*X1 - 0.368750000000000022E+00_dp)*0.177777777777777715E+02_dp
155 0 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
156 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 13))
157 : ELSE
158 0 : TG1 = (2*X1 - 0.368750000000000022E+00_dp)*0.177777777777777715E+02_dp
159 0 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
160 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 14))
161 : END IF
162 : END IF
163 : ELSE
164 139048 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
165 92520 : IF (X1 <= 0.128124999999999989E+00_dp) THEN
166 45604 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
167 32312 : IF (X1 <= 0.114062499999999997E+00_dp) THEN
168 23440 : TG1 = (2*X1 - 0.214062499999999989E+00_dp)*0.711111111111111569E+02_dp
169 23440 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
170 23440 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 15))
171 : ELSE
172 8872 : TG1 = (2*X1 - 0.242187500000000000E+00_dp)*0.711111111111111569E+02_dp
173 8872 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
174 8872 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 16))
175 : END IF
176 : ELSE
177 13292 : TG1 = (2*X1 - 0.228124999999999994E+00_dp)*0.355555555555555785E+02_dp
178 13292 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
179 13292 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 17))
180 : END IF
181 : ELSE
182 46916 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
183 3580 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
184 3580 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
185 3580 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 18))
186 : ELSE
187 43336 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
188 43336 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
189 43336 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 19))
190 : END IF
191 : END IF
192 : ELSE
193 46528 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
194 22908 : TG1 = (2*X1 - 0.368750000000000022E+00_dp)*0.177777777777777715E+02_dp
195 22908 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
196 22908 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 20))
197 : ELSE
198 23620 : TG1 = (2*X1 - 0.368750000000000022E+00_dp)*0.177777777777777715E+02_dp
199 23620 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
200 23620 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 21))
201 : END IF
202 : END IF
203 : END IF
204 : END IF
205 : ELSE
206 14624 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
207 0 : TG1 = (2*X1 - 0.537500000000000089E+00_dp)*0.888888888888888928E+01_dp
208 0 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
209 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 22))
210 : ELSE
211 14624 : IF (X2 <= 0.187500000000000000E+00_dp) THEN
212 5864 : IF (X1 <= 0.268750000000000044E+00_dp) THEN
213 3792 : TG1 = (2*X1 - 0.481250000000000067E+00_dp)*0.177777777777777715E+02_dp
214 3792 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
215 3792 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 23))
216 : ELSE
217 2072 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.177777777777777892E+02_dp
218 2072 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
219 2072 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 24))
220 : END IF
221 : ELSE
222 8760 : IF (X1 <= 0.268750000000000044E+00_dp) THEN
223 6996 : TG1 = (2*X1 - 0.481250000000000067E+00_dp)*0.177777777777777715E+02_dp
224 6996 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
225 6996 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 25))
226 : ELSE
227 1764 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.177777777777777892E+02_dp
228 1764 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
229 1764 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 26))
230 : END IF
231 : END IF
232 : END IF
233 : END IF
234 : ELSE
235 45840 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
236 8 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.444444444444444375E+01_dp
237 8 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
238 8 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 27))
239 : ELSE
240 45832 : IF (X2 <= 0.187500000000000000E+00_dp) THEN
241 19168 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.444444444444444375E+01_dp
242 19168 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
243 19168 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 28))
244 : ELSE
245 26664 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.444444444444444375E+01_dp
246 26664 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
247 26664 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 29))
248 : END IF
249 : END IF
250 : END IF
251 : ELSE
252 394712 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
253 209024 : TG1 = (2*X1 - 0.155000000000000004E+01_dp)*0.222222222222222232E+01_dp
254 209024 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
255 209024 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 30))
256 : ELSE
257 185688 : TG1 = (2*X1 - 0.155000000000000004E+01_dp)*0.222222222222222232E+01_dp
258 185688 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
259 185688 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 31))
260 : END IF
261 : END IF
262 : ELSE
263 234884 : IF (X1 <= 0.550000000000000044E+00_dp) THEN
264 155656 : IF (X1 <= 0.325000000000000011E+00_dp) THEN
265 116804 : IF (X1 <= 0.212500000000000022E+00_dp) THEN
266 54872 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
267 9104 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
268 9104 : IF (X1 <= 0.128124999999999989E+00_dp) THEN
269 0 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
270 0 : TG1 = (2*X1 - 0.228124999999999994E+00_dp)*0.355555555555555785E+02_dp
271 0 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
272 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 32))
273 : ELSE
274 0 : TG1 = (2*X1 - 0.228124999999999994E+00_dp)*0.355555555555555785E+02_dp
275 0 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
276 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 33))
277 : END IF
278 : ELSE
279 9104 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
280 9104 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
281 9104 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
282 9104 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
283 9104 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 34))
284 : ELSE
285 0 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
286 0 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
287 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 35))
288 : END IF
289 : ELSE
290 0 : TG1 = (2*X1 - 0.284374999999999989E+00_dp)*0.355555555555555429E+02_dp
291 0 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
292 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 36))
293 : END IF
294 : END IF
295 : ELSE
296 0 : TG1 = (2*X1 - 0.256249999999999978E+00_dp)*0.177777777777777786E+02_dp
297 0 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
298 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 37))
299 : END IF
300 : ELSE
301 45768 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
302 45768 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
303 32804 : IF (X1 <= 0.184375000000000011E+00_dp) THEN
304 25304 : TG1 = (2*X1 - 0.340625000000000011E+00_dp)*0.355555555555555429E+02_dp
305 25304 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
306 25304 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 38))
307 : ELSE
308 7500 : TG1 = (2*X1 - 0.396875000000000033E+00_dp)*0.355555555555555429E+02_dp
309 7500 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
310 7500 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 39))
311 : END IF
312 : ELSE
313 12964 : IF (X1 <= 0.184375000000000011E+00_dp) THEN
314 0 : TG1 = (2*X1 - 0.340625000000000011E+00_dp)*0.355555555555555429E+02_dp
315 0 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
316 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 40))
317 : ELSE
318 12964 : TG1 = (2*X1 - 0.396875000000000033E+00_dp)*0.355555555555555429E+02_dp
319 12964 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
320 12964 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 41))
321 : END IF
322 : END IF
323 : ELSE
324 0 : IF (X1 <= 0.184375000000000011E+00_dp) THEN
325 0 : TG1 = (2*X1 - 0.340625000000000011E+00_dp)*0.355555555555555429E+02_dp
326 0 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
327 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 42))
328 : ELSE
329 0 : TG1 = (2*X1 - 0.396875000000000033E+00_dp)*0.355555555555555429E+02_dp
330 0 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
331 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 43))
332 : END IF
333 : END IF
334 : END IF
335 : ELSE
336 61932 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
337 51120 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
338 26340 : IF (X1 <= 0.268750000000000044E+00_dp) THEN
339 4672 : TG1 = (2*X1 - 0.481250000000000067E+00_dp)*0.177777777777777715E+02_dp
340 4672 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
341 4672 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 44))
342 : ELSE
343 21668 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.177777777777777892E+02_dp
344 21668 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
345 21668 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 45))
346 : END IF
347 : ELSE
348 24780 : IF (X1 <= 0.268750000000000044E+00_dp) THEN
349 20312 : TG1 = (2*X1 - 0.481250000000000067E+00_dp)*0.177777777777777715E+02_dp
350 20312 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
351 20312 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 46))
352 : ELSE
353 4468 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.177777777777777892E+02_dp
354 4468 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
355 4468 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 47))
356 : END IF
357 : END IF
358 : ELSE
359 10812 : IF (X1 <= 0.268750000000000044E+00_dp) THEN
360 1076 : TG1 = (2*X1 - 0.481250000000000067E+00_dp)*0.177777777777777715E+02_dp
361 1076 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
362 1076 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 48))
363 : ELSE
364 9736 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.177777777777777892E+02_dp
365 9736 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
366 9736 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 49))
367 : END IF
368 : END IF
369 : END IF
370 : ELSE
371 38852 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
372 25116 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
373 11400 : TG1 = (2*X1 - 0.762499999999999956E+00_dp)*0.888888888888888928E+01_dp
374 11400 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
375 11400 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 50))
376 : ELSE
377 13716 : TG1 = (2*X1 - 0.987500000000000044E+00_dp)*0.888888888888888573E+01_dp
378 13716 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
379 13716 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 51))
380 : END IF
381 : ELSE
382 13736 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
383 3896 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
384 3104 : TG1 = (2*X1 - 0.762499999999999956E+00_dp)*0.888888888888888928E+01_dp
385 3104 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
386 3104 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 52))
387 : ELSE
388 792 : TG1 = (2*X1 - 0.987500000000000044E+00_dp)*0.888888888888888573E+01_dp
389 792 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
390 792 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 53))
391 : END IF
392 : ELSE
393 9840 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
394 3136 : TG1 = (2*X1 - 0.762499999999999956E+00_dp)*0.888888888888888928E+01_dp
395 3136 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
396 3136 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 54))
397 : ELSE
398 6704 : TG1 = (2*X1 - 0.987500000000000044E+00_dp)*0.888888888888888573E+01_dp
399 6704 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
400 6704 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 55))
401 : END IF
402 : END IF
403 : END IF
404 : END IF
405 : ELSE
406 79228 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
407 66426 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
408 37350 : TG1 = (2*X1 - 0.155000000000000004E+01_dp)*0.222222222222222232E+01_dp
409 37350 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
410 37350 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 56))
411 : ELSE
412 29076 : TG1 = (2*X1 - 0.155000000000000004E+01_dp)*0.222222222222222232E+01_dp
413 29076 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
414 29076 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 57))
415 : END IF
416 : ELSE
417 12802 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
418 5278 : TG1 = (2*X1 - 0.155000000000000004E+01_dp)*0.222222222222222232E+01_dp
419 5278 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
420 5278 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 58))
421 : ELSE
422 7524 : IF (X1 <= 0.775000000000000022E+00_dp) THEN
423 1224 : TG1 = (2*X1 - 0.132500000000000018E+01_dp)*0.444444444444444464E+01_dp
424 1224 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
425 1224 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 59))
426 : ELSE
427 6300 : TG1 = (2*X1 - 0.177499999999999991E+01_dp)*0.444444444444444464E+01_dp
428 6300 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
429 6300 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 60))
430 : END IF
431 : END IF
432 : END IF
433 : END IF
434 : END IF
435 : ELSE
436 200 : IF (X1 <= 0.550000000000000044E+00_dp) THEN
437 116 : IF (X1 <= 0.325000000000000011E+00_dp) THEN
438 0 : IF (X1 <= 0.212500000000000022E+00_dp) THEN
439 0 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
440 0 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
441 0 : TG1 = (2*X1 - 0.256249999999999978E+00_dp)*0.177777777777777786E+02_dp
442 0 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
443 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 61))
444 : ELSE
445 0 : TG1 = (2*X1 - 0.368750000000000022E+00_dp)*0.177777777777777715E+02_dp
446 0 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
447 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 62))
448 : END IF
449 : ELSE
450 0 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.888888888888888751E+01_dp
451 0 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
452 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 63))
453 : END IF
454 : ELSE
455 0 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
456 0 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
457 0 : IF (X1 <= 0.268750000000000044E+00_dp) THEN
458 0 : TG1 = (2*X1 - 0.481250000000000067E+00_dp)*0.177777777777777715E+02_dp
459 0 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
460 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 64))
461 : ELSE
462 0 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.177777777777777892E+02_dp
463 0 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
464 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 65))
465 : END IF
466 : ELSE
467 0 : TG1 = (2*X1 - 0.537500000000000089E+00_dp)*0.888888888888888928E+01_dp
468 0 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
469 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 66))
470 : END IF
471 : ELSE
472 0 : TG1 = (2*X1 - 0.537500000000000089E+00_dp)*0.888888888888888928E+01_dp
473 0 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
474 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 67))
475 : END IF
476 : END IF
477 : ELSE
478 116 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
479 116 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
480 116 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
481 36 : TG1 = (2*X1 - 0.762499999999999956E+00_dp)*0.888888888888888928E+01_dp
482 36 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
483 36 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 68))
484 : ELSE
485 80 : TG1 = (2*X1 - 0.987500000000000044E+00_dp)*0.888888888888888573E+01_dp
486 80 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
487 80 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 69))
488 : END IF
489 : ELSE
490 0 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
491 0 : TG1 = (2*X1 - 0.762499999999999956E+00_dp)*0.888888888888888928E+01_dp
492 0 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
493 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 70))
494 : ELSE
495 0 : TG1 = (2*X1 - 0.987500000000000044E+00_dp)*0.888888888888888573E+01_dp
496 0 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
497 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 71))
498 : END IF
499 : END IF
500 : ELSE
501 0 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
502 0 : TG1 = (2*X1 - 0.762499999999999956E+00_dp)*0.888888888888888928E+01_dp
503 0 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
504 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 72))
505 : ELSE
506 0 : TG1 = (2*X1 - 0.987500000000000044E+00_dp)*0.888888888888888573E+01_dp
507 0 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
508 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 73))
509 : END IF
510 : END IF
511 : END IF
512 : ELSE
513 84 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
514 84 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
515 84 : IF (X1 <= 0.775000000000000022E+00_dp) THEN
516 16 : TG1 = (2*X1 - 0.132500000000000018E+01_dp)*0.444444444444444464E+01_dp
517 16 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
518 16 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 74))
519 : ELSE
520 68 : TG1 = (2*X1 - 0.177499999999999991E+01_dp)*0.444444444444444464E+01_dp
521 68 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
522 68 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 75))
523 : END IF
524 : ELSE
525 0 : IF (X1 <= 0.775000000000000022E+00_dp) THEN
526 0 : TG1 = (2*X1 - 0.132500000000000018E+01_dp)*0.444444444444444464E+01_dp
527 0 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
528 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 76))
529 : ELSE
530 0 : TG1 = (2*X1 - 0.177499999999999991E+01_dp)*0.444444444444444464E+01_dp
531 0 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
532 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 77))
533 : END IF
534 : END IF
535 : ELSE
536 0 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
537 0 : TG1 = (2*X1 - 0.155000000000000004E+01_dp)*0.222222222222222232E+01_dp
538 0 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
539 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 78))
540 : ELSE
541 0 : TG1 = (2*X1 - 0.155000000000000004E+01_dp)*0.222222222222222232E+01_dp
542 0 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
543 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 79))
544 : END IF
545 : END IF
546 : END IF
547 : END IF
548 : ELSE
549 309116 : X1 = SQRT(81.0_dp/T)
550 309116 : lower = (SQRT(T) - SQRT(eps))/n
551 309116 : IF (R >= lower) THEN
552 309116 : X2 = lower/R
553 309116 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
554 74660 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
555 0 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
556 0 : TG1 = (2*X1 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
557 0 : TG2 = (2*X2 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
558 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 80))
559 : ELSE
560 0 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
561 0 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
562 0 : TG1 = (2*X1 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
563 0 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
564 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 81))
565 : ELSE
566 0 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
567 0 : TG1 = (2*X1 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
568 0 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
569 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 82))
570 : ELSE
571 0 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
572 0 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
573 0 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
574 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 83))
575 : ELSE
576 0 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
577 0 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
578 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 84))
579 : END IF
580 : END IF
581 : END IF
582 : ELSE
583 0 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
584 0 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
585 0 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
586 0 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
587 0 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
588 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 85))
589 : ELSE
590 0 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
591 0 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
592 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 86))
593 : END IF
594 : ELSE
595 0 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
596 0 : IF (X2 <= 0.468750000000000000E+00_dp) THEN
597 0 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
598 0 : TG2 = (2*X2 - 0.906250000000000000E+00_dp)*0.320000000000000000E+02_dp
599 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 87))
600 : ELSE
601 0 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
602 0 : TG2 = (2*X2 - 0.968750000000000000E+00_dp)*0.320000000000000000E+02_dp
603 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 88))
604 : END IF
605 : ELSE
606 0 : IF (X2 <= 0.468750000000000000E+00_dp) THEN
607 0 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
608 0 : TG2 = (2*X2 - 0.906250000000000000E+00_dp)*0.320000000000000000E+02_dp
609 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 89))
610 : ELSE
611 0 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
612 0 : TG2 = (2*X2 - 0.968750000000000000E+00_dp)*0.320000000000000000E+02_dp
613 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 90))
614 : END IF
615 : END IF
616 : END IF
617 : ELSE
618 0 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
619 0 : IF (X2 <= 0.406250000000000000E+00_dp) THEN
620 0 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
621 0 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
622 0 : TG2 = (2*X2 - 0.781250000000000000E+00_dp)*0.320000000000000000E+02_dp
623 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 91))
624 : ELSE
625 0 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
626 0 : TG2 = (2*X2 - 0.781250000000000000E+00_dp)*0.320000000000000000E+02_dp
627 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 92))
628 : END IF
629 : ELSE
630 0 : TG1 = (2*X1 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
631 0 : TG2 = (2*X2 - 0.843750000000000000E+00_dp)*0.320000000000000000E+02_dp
632 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 93))
633 : END IF
634 : ELSE
635 0 : TG1 = (2*X1 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
636 0 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
637 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 94))
638 : END IF
639 : END IF
640 : END IF
641 : END IF
642 : ELSE
643 74660 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
644 0 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
645 0 : TG1 = (2*X1 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
646 0 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
647 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 95))
648 : ELSE
649 0 : IF (X2 <= 0.187500000000000000E+00_dp) THEN
650 0 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
651 0 : TG1 = (2*X1 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
652 0 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
653 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 96))
654 : ELSE
655 0 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
656 0 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
657 0 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
658 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 97))
659 : ELSE
660 0 : IF (X2 <= 0.156250000000000000E+00_dp) THEN
661 0 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
662 0 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
663 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 98))
664 : ELSE
665 0 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
666 0 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
667 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 99))
668 : END IF
669 : END IF
670 : END IF
671 : ELSE
672 0 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
673 0 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
674 0 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
675 0 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
676 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 100))
677 : ELSE
678 0 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
679 0 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
680 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 101))
681 : END IF
682 : ELSE
683 0 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
684 0 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
685 0 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
686 0 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
687 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 102))
688 : ELSE
689 0 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
690 0 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
691 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 103))
692 : END IF
693 : ELSE
694 0 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
695 0 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
696 0 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
697 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 104))
698 : ELSE
699 0 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
700 0 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
701 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 105))
702 : END IF
703 : END IF
704 : END IF
705 : END IF
706 : END IF
707 : ELSE
708 74660 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
709 13000 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
710 0 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
711 0 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
712 0 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
713 0 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
714 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 106))
715 : ELSE
716 0 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
717 0 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
718 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 107))
719 : END IF
720 : ELSE
721 0 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
722 0 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
723 0 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
724 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 108))
725 : ELSE
726 0 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
727 0 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
728 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 109))
729 : END IF
730 : END IF
731 : ELSE
732 13000 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
733 816 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
734 0 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
735 0 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
736 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 110))
737 : ELSE
738 816 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
739 816 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
740 816 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 111))
741 : END IF
742 : ELSE
743 12184 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
744 0 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
745 0 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
746 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 112))
747 : ELSE
748 12184 : IF (X1 <= 0.937500000000000000E+00_dp) THEN
749 1144 : TG1 = (2*X1 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
750 1144 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
751 1144 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 113))
752 : ELSE
753 11040 : TG1 = (2*X1 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
754 11040 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
755 11040 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 114))
756 : END IF
757 : END IF
758 : END IF
759 : END IF
760 : ELSE
761 61660 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
762 84 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
763 0 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
764 0 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
765 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 115))
766 : ELSE
767 84 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
768 0 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
769 0 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
770 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 116))
771 : ELSE
772 84 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
773 84 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
774 84 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 117))
775 : END IF
776 : END IF
777 : ELSE
778 61576 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
779 41728 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
780 8848 : IF (X1 <= 0.812500000000000000E+00_dp) THEN
781 0 : TG1 = (2*X1 - 0.156250000000000000E+01_dp)*0.160000000000000000E+02_dp
782 0 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
783 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 118))
784 : ELSE
785 8848 : TG1 = (2*X1 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
786 8848 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
787 8848 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 119))
788 : END IF
789 : ELSE
790 32880 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
791 32880 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
792 32880 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 120))
793 : END IF
794 : ELSE
795 19848 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
796 19848 : IF (X1 <= 0.937500000000000000E+00_dp) THEN
797 16364 : TG1 = (2*X1 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
798 16364 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
799 16364 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 121))
800 : ELSE
801 3484 : TG1 = (2*X1 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
802 3484 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
803 3484 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 122))
804 : END IF
805 : ELSE
806 0 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
807 0 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
808 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 123))
809 : END IF
810 : END IF
811 : END IF
812 : END IF
813 : END IF
814 : END IF
815 : ELSE
816 234456 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
817 177932 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
818 40760 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
819 0 : TG1 = (2*X1 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
820 0 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
821 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 124))
822 : ELSE
823 40760 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
824 0 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
825 0 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
826 0 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
827 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 125))
828 : ELSE
829 0 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
830 0 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
831 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 126))
832 : END IF
833 : ELSE
834 40760 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
835 120 : IF (X2 <= 0.687500000000000000E+00_dp) THEN
836 0 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
837 0 : TG2 = (2*X2 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
838 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 127))
839 : ELSE
840 120 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
841 120 : TG2 = (2*X2 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
842 120 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 128))
843 : END IF
844 : ELSE
845 40640 : IF (X2 <= 0.687500000000000000E+00_dp) THEN
846 12296 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
847 0 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
848 0 : TG2 = (2*X2 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
849 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 129))
850 : ELSE
851 12296 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
852 12296 : TG2 = (2*X2 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
853 12296 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 130))
854 : END IF
855 : ELSE
856 28344 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
857 28344 : TG2 = (2*X2 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
858 28344 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 131))
859 : END IF
860 : END IF
861 : END IF
862 : END IF
863 : ELSE
864 137172 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
865 97068 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
866 26924 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
867 0 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
868 0 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
869 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 132))
870 : ELSE
871 26924 : IF (X2 <= 0.812500000000000000E+00_dp) THEN
872 0 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
873 0 : TG2 = (2*X2 - 0.156250000000000000E+01_dp)*0.160000000000000000E+02_dp
874 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 133))
875 : ELSE
876 26924 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
877 192 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
878 192 : TG2 = (2*X2 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
879 192 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 134))
880 : ELSE
881 26732 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
882 26732 : TG2 = (2*X2 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
883 26732 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 135))
884 : END IF
885 : END IF
886 : END IF
887 : ELSE
888 70144 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
889 47684 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
890 16628 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
891 0 : TG1 = (2*X1 - 0.625000000000000000E-01_dp)*0.160000000000000000E+02_dp
892 0 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
893 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 136))
894 : ELSE
895 16628 : IF (X2 <= 0.906250000000000000E+00_dp) THEN
896 0 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
897 0 : TG2 = (2*X2 - 0.178125000000000000E+01_dp)*0.320000000000000000E+02_dp
898 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 137))
899 : ELSE
900 16628 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
901 16628 : TG2 = (2*X2 - 0.184375000000000000E+01_dp)*0.320000000000000000E+02_dp
902 16628 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 138))
903 : END IF
904 : END IF
905 : ELSE
906 31056 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
907 25884 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
908 12452 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
909 0 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
910 0 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
911 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 139))
912 : ELSE
913 12452 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
914 12452 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
915 12452 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 140))
916 : END IF
917 : ELSE
918 13432 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
919 9704 : IF (X2 <= 0.984375000000000000E+00_dp) THEN
920 6284 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
921 6284 : TG2 = (2*X2 - 0.195312500000000000E+01_dp)*0.640000000000000000E+02_dp
922 6284 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 141))
923 : ELSE
924 3420 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
925 3420 : TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
926 3420 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 142))
927 : END IF
928 : ELSE
929 3728 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
930 3728 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
931 3728 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 143))
932 : END IF
933 : END IF
934 : ELSE
935 5172 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
936 5172 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
937 4508 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
938 4508 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
939 4508 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 144))
940 : ELSE
941 664 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
942 664 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
943 664 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 145))
944 : END IF
945 : ELSE
946 0 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
947 0 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
948 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 146))
949 : END IF
950 : END IF
951 : END IF
952 : ELSE
953 22460 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
954 22460 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
955 19548 : IF (X2 <= 0.906250000000000000E+00_dp) THEN
956 17220 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
957 17220 : TG2 = (2*X2 - 0.178125000000000000E+01_dp)*0.320000000000000000E+02_dp
958 17220 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 147))
959 : ELSE
960 2328 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
961 2328 : TG2 = (2*X2 - 0.184375000000000000E+01_dp)*0.320000000000000000E+02_dp
962 2328 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 148))
963 : END IF
964 : ELSE
965 2912 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
966 2912 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
967 2912 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 149))
968 : END IF
969 : ELSE
970 0 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
971 0 : TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
972 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 150))
973 : END IF
974 : END IF
975 : END IF
976 : ELSE
977 40104 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
978 40104 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
979 39736 : IF (X2 <= 0.812500000000000000E+00_dp) THEN
980 30532 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
981 15536 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
982 15536 : TG2 = (2*X2 - 0.156250000000000000E+01_dp)*0.160000000000000000E+02_dp
983 15536 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 151))
984 : ELSE
985 14996 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
986 14996 : TG2 = (2*X2 - 0.156250000000000000E+01_dp)*0.160000000000000000E+02_dp
987 14996 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 152))
988 : END IF
989 : ELSE
990 9204 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
991 9204 : TG2 = (2*X2 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
992 9204 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 153))
993 : END IF
994 : ELSE
995 368 : IF (X2 <= 0.812500000000000000E+00_dp) THEN
996 368 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
997 368 : TG2 = (2*X2 - 0.156250000000000000E+01_dp)*0.160000000000000000E+02_dp
998 368 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 154))
999 : ELSE
1000 0 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
1001 0 : TG2 = (2*X2 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
1002 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 155))
1003 : END IF
1004 : END IF
1005 : ELSE
1006 0 : TG1 = (2*X1 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
1007 0 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1008 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 156))
1009 : END IF
1010 : END IF
1011 : END IF
1012 : ELSE
1013 56524 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
1014 56524 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
1015 52284 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
1016 34140 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
1017 9856 : IF (X2 <= 0.562500000000000000E+00_dp) THEN
1018 0 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
1019 0 : TG2 = (2*X2 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1020 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 157))
1021 : ELSE
1022 9856 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1023 1272 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1024 1272 : TG2 = (2*X2 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1025 1272 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 158))
1026 : ELSE
1027 8584 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1028 8584 : TG2 = (2*X2 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1029 8584 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 159))
1030 : END IF
1031 : END IF
1032 : ELSE
1033 24284 : IF (X2 <= 0.562500000000000000E+00_dp) THEN
1034 18040 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1035 18040 : TG2 = (2*X2 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1036 18040 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 160))
1037 : ELSE
1038 6244 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1039 6244 : TG2 = (2*X2 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1040 6244 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 161))
1041 : END IF
1042 : END IF
1043 : ELSE
1044 18144 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
1045 18144 : IF (X2 <= 0.687500000000000000E+00_dp) THEN
1046 18144 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
1047 18144 : TG2 = (2*X2 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
1048 18144 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 162))
1049 : ELSE
1050 0 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
1051 0 : TG2 = (2*X2 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
1052 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 163))
1053 : END IF
1054 : ELSE
1055 0 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1056 0 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1057 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 164))
1058 : END IF
1059 : END IF
1060 : ELSE
1061 4240 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
1062 4240 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
1063 4240 : IF (X2 <= 0.562500000000000000E+00_dp) THEN
1064 4240 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1065 4240 : TG2 = (2*X2 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1066 4240 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 165))
1067 : ELSE
1068 0 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1069 0 : TG2 = (2*X2 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1070 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 166))
1071 : END IF
1072 : ELSE
1073 0 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1074 0 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
1075 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 167))
1076 : END IF
1077 : ELSE
1078 0 : TG1 = (2*X1 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1079 0 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1080 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 168))
1081 : END IF
1082 : END IF
1083 : ELSE
1084 0 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
1085 0 : TG1 = (2*X1 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
1086 0 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1087 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 169))
1088 : ELSE
1089 0 : TG1 = (2*X1 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1090 0 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1091 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 170))
1092 : END IF
1093 : END IF
1094 : END IF
1095 : END IF
1096 : ELSE
1097 0 : DO i = 1, NDERIV + 1
1098 0 : RES(i) = 0.0_dp
1099 : END DO
1100 : END IF
1101 : END IF
1102 1160608 : END SUBROUTINE trunc_CS_poly_n20
1103 :
1104 : ! **************************************************************************************************
1105 : !> \brief ...
1106 : !> \param Nder the number of derivatives that will actually be used
1107 : !> \param iunit contains the data file to initialize the table
1108 : !> \param mepos ...
1109 : !> \param group ...
1110 : ! **************************************************************************************************
1111 2 : SUBROUTINE INIT(Nder, iunit, mepos, group)
1112 : INTEGER, INTENT(IN) :: Nder, iunit, mepos
1113 :
1114 : CLASS(mp_comm_type), INTENT(IN) :: group
1115 :
1116 : INTEGER :: I
1117 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1118 :
1119 2 : patches = 170
1120 2 : IF (Nder > nderiv_max) CPABORT("Reading data for initialization of C0 failed")
1121 2 : nderiv_init = Nder
1122 2 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1123 : ! round up to a multiple of 32 to give some generous alignment for each C0
1124 8 : ALLOCATE (C0(32*((31 + (Nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1125 : ! valgrind workaround
1126 609622 : C0 = HUGE(0.0_dp)
1127 2 : IF (mepos == 0) THEN
1128 2 : ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1129 342 : DO I = 1, patches
1130 340 : READ (iunit, *) chunk
1131 601462 : C0(1:(Nder + 1)*(degree + 1)*(degree + 2)/2, I) = chunk(1:(Nder + 1)*(degree + 1)*(degree + 2)/2)
1132 : END DO
1133 2 : DEALLOCATE (chunk)
1134 : END IF
1135 2 : CALL group%bcast(C0, 0)
1136 2 : END SUBROUTINE INIT
1137 :
1138 : ! **************************************************************************************************
1139 : !> \brief ...
1140 : ! **************************************************************************************************
1141 2 : SUBROUTINE free_C0()
1142 2 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1143 2 : nderiv_init = -1
1144 2 : END SUBROUTINE free_C0
1145 :
1146 : ! **************************************************************************************************
1147 : !> \brief ...
1148 : !> \param RES ...
1149 : !> \param NDERIV ...
1150 : !> \param TG1 ...
1151 : !> \param TG2 ...
1152 : !> \param C0 ...
1153 : ! **************************************************************************************************
1154 1160608 : SUBROUTINE PD2VAL(RES, NDERIV, TG1, TG2, C0)
1155 : REAL(KIND=dp), INTENT(OUT) :: res(*)
1156 : INTEGER, INTENT(IN) :: NDERIV
1157 : REAL(KIND=dp), INTENT(IN) :: TG1, TG2, C0(136, *)
1158 :
1159 : REAL(KIND=dp), PARAMETER :: SQRT2 = 1.4142135623730950488016887242096980785696718753_dp
1160 :
1161 : INTEGER :: K
1162 : REAL(KIND=dp) :: T1(0:15), T2(0:15)
1163 :
1164 1160608 : T1(0) = 1.0_dp
1165 1160608 : T2(0) = 1.0_dp
1166 1160608 : T1(1) = SQRT2*TG1
1167 1160608 : T2(1) = SQRT2*TG2
1168 1160608 : T1(2) = 2*TG1*T1(1) - SQRT2
1169 1160608 : T2(2) = 2*TG2*T2(1) - SQRT2
1170 1160608 : T1(3) = 2*TG1*T1(2) - T1(1)
1171 1160608 : T2(3) = 2*TG2*T2(2) - T2(1)
1172 1160608 : T1(4) = 2*TG1*T1(3) - T1(2)
1173 1160608 : T2(4) = 2*TG2*T2(3) - T2(2)
1174 1160608 : T1(5) = 2*TG1*T1(4) - T1(3)
1175 1160608 : T2(5) = 2*TG2*T2(4) - T2(3)
1176 1160608 : T1(6) = 2*TG1*T1(5) - T1(4)
1177 1160608 : T2(6) = 2*TG2*T2(5) - T2(4)
1178 1160608 : T1(7) = 2*TG1*T1(6) - T1(5)
1179 1160608 : T2(7) = 2*TG2*T2(6) - T2(5)
1180 1160608 : T1(8) = 2*TG1*T1(7) - T1(6)
1181 1160608 : T2(8) = 2*TG2*T2(7) - T2(6)
1182 1160608 : T1(9) = 2*TG1*T1(8) - T1(7)
1183 1160608 : T2(9) = 2*TG2*T2(8) - T2(7)
1184 1160608 : T1(10) = 2*TG1*T1(9) - T1(8)
1185 1160608 : T2(10) = 2*TG2*T2(9) - T2(8)
1186 1160608 : T1(11) = 2*TG1*T1(10) - T1(9)
1187 1160608 : T2(11) = 2*TG2*T2(10) - T2(9)
1188 1160608 : T1(12) = 2*TG1*T1(11) - T1(10)
1189 1160608 : T2(12) = 2*TG2*T2(11) - T2(10)
1190 1160608 : T1(13) = 2*TG1*T1(12) - T1(11)
1191 1160608 : T2(13) = 2*TG2*T2(12) - T2(11)
1192 1160608 : T1(14) = 2*TG1*T1(13) - T1(12)
1193 1160608 : T2(14) = 2*TG2*T2(13) - T2(12)
1194 1160608 : T1(15) = 2*TG1*T1(14) - T1(13)
1195 1160608 : T2(15) = 2*TG2*T2(14) - T2(13)
1196 6054088 : DO K = 1, NDERIV + 1
1197 4893480 : RES(K) = 0.0_dp
1198 83189160 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:15), C0(1:16, K))*T2(0)
1199 78295680 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:14), C0(17:31, K))*T2(1)
1200 73402200 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:13), C0(32:45, K))*T2(2)
1201 68508720 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:12), C0(46:58, K))*T2(3)
1202 63615240 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:11), C0(59:70, K))*T2(4)
1203 58721760 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:10), C0(71:81, K))*T2(5)
1204 53828280 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:9), C0(82:91, K))*T2(6)
1205 48934800 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:8), C0(92:100, K))*T2(7)
1206 44041320 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:7), C0(101:108, K))*T2(8)
1207 39147840 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:6), C0(109:115, K))*T2(9)
1208 34254360 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:5), C0(116:121, K))*T2(10)
1209 29360880 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:4), C0(122:126, K))*T2(11)
1210 24467400 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:3), C0(127:130, K))*T2(12)
1211 19573920 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:2), C0(131:133, K))*T2(13)
1212 14680440 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:1), C0(134:135, K))*T2(14)
1213 10947568 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:0), C0(136:136, K))*T2(15)
1214 : END DO
1215 1160608 : END SUBROUTINE PD2VAL
1216 :
1217 : END MODULE t_sh_p_s_c
|
