(git:374b731)
Loading...
Searching...
No Matches
t_sh_p_s_c.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 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! **************************************************************************************************
39MODULE 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
55CONTAINS
56
57! **************************************************************************************************
58!> \brief ...
59!> \param RES ...
60!> \param R ...
61!> \param T ...
62!> \param NDERIV ...
63! **************************************************************************************************
64 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 IF (t <= 81.0_dp) THEN
75 x2 = 1/(1 + r)
76 x1 = 1/(1 + sqrt(t))
77 IF (x2 <= 0.500000000000000000e+00_dp) THEN
78 IF (x2 <= 0.250000000000000000e+00_dp) THEN
79 IF (x1 <= 0.550000000000000044e+00_dp) THEN
80 IF (x1 <= 0.325000000000000011e+00_dp) THEN
81 IF (x1 <= 0.212500000000000022e+00_dp) THEN
82 IF (x2 <= 0.125000000000000000e+00_dp) THEN
83 IF (x2 <= 0.625000000000000000e-01_dp) THEN
84 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.888888888888888751e+01_dp
85 tg2 = (2*x2 - 0.625000000000000000e-01_dp)*0.160000000000000000e+02_dp
86 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 1))
87 ELSE
88 IF (x1 <= 0.156250000000000000e+00_dp) THEN
89 IF (x2 <= 0.937500000000000000e-01_dp) THEN
90 IF (x1 <= 0.128124999999999989e+00_dp) THEN
91 tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
92 tg2 = (2*x2 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
93 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 2))
94 ELSE
95 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
96 tg2 = (2*x2 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
97 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 3))
98 END IF
99 ELSE
100 IF (x1 <= 0.128124999999999989e+00_dp) THEN
101 IF (x2 <= 0.109375000000000000e+00_dp) THEN
102 tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
103 tg2 = (2*x2 - 0.203125000000000000e+00_dp)*0.640000000000000000e+02_dp
104 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 4))
105 ELSE
106 tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
107 tg2 = (2*x2 - 0.234375000000000000e+00_dp)*0.640000000000000000e+02_dp
108 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 5))
109 END IF
110 ELSE
111 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
112 tg2 = (2*x2 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
113 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 6))
114 END IF
115 END IF
116 ELSE
117 tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
118 tg2 = (2*x2 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
119 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 7))
120 END IF
121 END IF
122 ELSE
123 IF (x2 <= 0.187500000000000000e+00_dp) THEN
124 IF (x1 <= 0.156250000000000000e+00_dp) THEN
125 IF (x2 <= 0.156250000000000000e+00_dp) THEN
126 IF (x1 <= 0.128124999999999989e+00_dp) THEN
127 tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
128 tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
129 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 8))
130 ELSE
131 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
132 tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
133 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 9))
134 END IF
135 ELSE
136 IF (x1 <= 0.128124999999999989e+00_dp) THEN
137 IF (x1 <= 0.114062499999999997e+00_dp) THEN
138 tg1 = (2*x1 - 0.214062499999999989e+00_dp)*0.711111111111111569e+02_dp
139 tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
140 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 10))
141 ELSE
142 tg1 = (2*x1 - 0.242187500000000000e+00_dp)*0.711111111111111569e+02_dp
143 tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
144 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 11))
145 END IF
146 ELSE
147 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
148 tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
149 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 12))
150 END IF
151 END IF
152 ELSE
153 IF (x2 <= 0.156250000000000000e+00_dp) THEN
154 tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
155 tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
156 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 13))
157 ELSE
158 tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
159 tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
160 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 14))
161 END IF
162 END IF
163 ELSE
164 IF (x1 <= 0.156250000000000000e+00_dp) THEN
165 IF (x1 <= 0.128124999999999989e+00_dp) THEN
166 IF (x2 <= 0.218750000000000000e+00_dp) THEN
167 IF (x1 <= 0.114062499999999997e+00_dp) THEN
168 tg1 = (2*x1 - 0.214062499999999989e+00_dp)*0.711111111111111569e+02_dp
169 tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
170 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 15))
171 ELSE
172 tg1 = (2*x1 - 0.242187500000000000e+00_dp)*0.711111111111111569e+02_dp
173 tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
174 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 16))
175 END IF
176 ELSE
177 tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
178 tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
179 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 17))
180 END IF
181 ELSE
182 IF (x2 <= 0.218750000000000000e+00_dp) THEN
183 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
184 tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
185 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 18))
186 ELSE
187 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
188 tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
189 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 19))
190 END IF
191 END IF
192 ELSE
193 IF (x2 <= 0.218750000000000000e+00_dp) THEN
194 tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
195 tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
196 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 20))
197 ELSE
198 tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
199 tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
200 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 21))
201 END IF
202 END IF
203 END IF
204 END IF
205 ELSE
206 IF (x2 <= 0.125000000000000000e+00_dp) THEN
207 tg1 = (2*x1 - 0.537500000000000089e+00_dp)*0.888888888888888928e+01_dp
208 tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
209 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 22))
210 ELSE
211 IF (x2 <= 0.187500000000000000e+00_dp) THEN
212 IF (x1 <= 0.268750000000000044e+00_dp) THEN
213 tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
214 tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
215 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 23))
216 ELSE
217 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
218 tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
219 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 24))
220 END IF
221 ELSE
222 IF (x1 <= 0.268750000000000044e+00_dp) THEN
223 tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
224 tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
225 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 25))
226 ELSE
227 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
228 tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
229 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 26))
230 END IF
231 END IF
232 END IF
233 END IF
234 ELSE
235 IF (x2 <= 0.125000000000000000e+00_dp) THEN
236 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.444444444444444375e+01_dp
237 tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
238 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 27))
239 ELSE
240 IF (x2 <= 0.187500000000000000e+00_dp) THEN
241 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.444444444444444375e+01_dp
242 tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
243 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 28))
244 ELSE
245 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.444444444444444375e+01_dp
246 tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
247 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 29))
248 END IF
249 END IF
250 END IF
251 ELSE
252 IF (x2 <= 0.125000000000000000e+00_dp) THEN
253 tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
254 tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
255 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 30))
256 ELSE
257 tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
258 tg2 = (2*x2 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
259 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 31))
260 END IF
261 END IF
262 ELSE
263 IF (x1 <= 0.550000000000000044e+00_dp) THEN
264 IF (x1 <= 0.325000000000000011e+00_dp) THEN
265 IF (x1 <= 0.212500000000000022e+00_dp) THEN
266 IF (x1 <= 0.156250000000000000e+00_dp) THEN
267 IF (x2 <= 0.375000000000000000e+00_dp) THEN
268 IF (x1 <= 0.128124999999999989e+00_dp) THEN
269 IF (x2 <= 0.312500000000000000e+00_dp) THEN
270 tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
271 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
272 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 32))
273 ELSE
274 tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
275 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
276 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 33))
277 END IF
278 ELSE
279 IF (x2 <= 0.312500000000000000e+00_dp) THEN
280 IF (x2 <= 0.281250000000000000e+00_dp) THEN
281 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
282 tg2 = (2*x2 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
283 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 34))
284 ELSE
285 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
286 tg2 = (2*x2 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
287 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 35))
288 END IF
289 ELSE
290 tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
291 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
292 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 36))
293 END IF
294 END IF
295 ELSE
296 tg1 = (2*x1 - 0.256249999999999978e+00_dp)*0.177777777777777786e+02_dp
297 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
298 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 37))
299 END IF
300 ELSE
301 IF (x2 <= 0.375000000000000000e+00_dp) THEN
302 IF (x2 <= 0.312500000000000000e+00_dp) THEN
303 IF (x1 <= 0.184375000000000011e+00_dp) THEN
304 tg1 = (2*x1 - 0.340625000000000011e+00_dp)*0.355555555555555429e+02_dp
305 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
306 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 38))
307 ELSE
308 tg1 = (2*x1 - 0.396875000000000033e+00_dp)*0.355555555555555429e+02_dp
309 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
310 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 39))
311 END IF
312 ELSE
313 IF (x1 <= 0.184375000000000011e+00_dp) THEN
314 tg1 = (2*x1 - 0.340625000000000011e+00_dp)*0.355555555555555429e+02_dp
315 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
316 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 40))
317 ELSE
318 tg1 = (2*x1 - 0.396875000000000033e+00_dp)*0.355555555555555429e+02_dp
319 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
320 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 41))
321 END IF
322 END IF
323 ELSE
324 IF (x1 <= 0.184375000000000011e+00_dp) THEN
325 tg1 = (2*x1 - 0.340625000000000011e+00_dp)*0.355555555555555429e+02_dp
326 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
327 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 42))
328 ELSE
329 tg1 = (2*x1 - 0.396875000000000033e+00_dp)*0.355555555555555429e+02_dp
330 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
331 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 43))
332 END IF
333 END IF
334 END IF
335 ELSE
336 IF (x2 <= 0.375000000000000000e+00_dp) THEN
337 IF (x2 <= 0.312500000000000000e+00_dp) THEN
338 IF (x1 <= 0.268750000000000044e+00_dp) THEN
339 tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
340 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
341 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 44))
342 ELSE
343 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
344 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
345 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 45))
346 END IF
347 ELSE
348 IF (x1 <= 0.268750000000000044e+00_dp) THEN
349 tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
350 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
351 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 46))
352 ELSE
353 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
354 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
355 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 47))
356 END IF
357 END IF
358 ELSE
359 IF (x1 <= 0.268750000000000044e+00_dp) THEN
360 tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
361 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
362 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 48))
363 ELSE
364 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
365 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
366 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 49))
367 END IF
368 END IF
369 END IF
370 ELSE
371 IF (x2 <= 0.375000000000000000e+00_dp) THEN
372 IF (x1 <= 0.437500000000000000e+00_dp) THEN
373 tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
374 tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
375 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 50))
376 ELSE
377 tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
378 tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
379 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 51))
380 END IF
381 ELSE
382 IF (x2 <= 0.437500000000000000e+00_dp) THEN
383 IF (x1 <= 0.437500000000000000e+00_dp) THEN
384 tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
385 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
386 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 52))
387 ELSE
388 tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
389 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
390 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 53))
391 END IF
392 ELSE
393 IF (x1 <= 0.437500000000000000e+00_dp) THEN
394 tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
395 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
396 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 54))
397 ELSE
398 tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
399 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
400 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 55))
401 END IF
402 END IF
403 END IF
404 END IF
405 ELSE
406 IF (x2 <= 0.375000000000000000e+00_dp) THEN
407 IF (x2 <= 0.312500000000000000e+00_dp) THEN
408 tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
409 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
410 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 56))
411 ELSE
412 tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
413 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
414 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 57))
415 END IF
416 ELSE
417 IF (x2 <= 0.437500000000000000e+00_dp) THEN
418 tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
419 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
420 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 58))
421 ELSE
422 IF (x1 <= 0.775000000000000022e+00_dp) THEN
423 tg1 = (2*x1 - 0.132500000000000018e+01_dp)*0.444444444444444464e+01_dp
424 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
425 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 59))
426 ELSE
427 tg1 = (2*x1 - 0.177499999999999991e+01_dp)*0.444444444444444464e+01_dp
428 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
429 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 IF (x1 <= 0.550000000000000044e+00_dp) THEN
437 IF (x1 <= 0.325000000000000011e+00_dp) THEN
438 IF (x1 <= 0.212500000000000022e+00_dp) THEN
439 IF (x2 <= 0.750000000000000000e+00_dp) THEN
440 IF (x1 <= 0.156250000000000000e+00_dp) THEN
441 tg1 = (2*x1 - 0.256249999999999978e+00_dp)*0.177777777777777786e+02_dp
442 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
443 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 61))
444 ELSE
445 tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
446 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
447 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 62))
448 END IF
449 ELSE
450 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.888888888888888751e+01_dp
451 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
452 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 63))
453 END IF
454 ELSE
455 IF (x2 <= 0.750000000000000000e+00_dp) THEN
456 IF (x2 <= 0.625000000000000000e+00_dp) THEN
457 IF (x1 <= 0.268750000000000044e+00_dp) THEN
458 tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
459 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
460 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 64))
461 ELSE
462 tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
463 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
464 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 65))
465 END IF
466 ELSE
467 tg1 = (2*x1 - 0.537500000000000089e+00_dp)*0.888888888888888928e+01_dp
468 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
469 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 66))
470 END IF
471 ELSE
472 tg1 = (2*x1 - 0.537500000000000089e+00_dp)*0.888888888888888928e+01_dp
473 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
474 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 67))
475 END IF
476 END IF
477 ELSE
478 IF (x2 <= 0.750000000000000000e+00_dp) THEN
479 IF (x2 <= 0.625000000000000000e+00_dp) THEN
480 IF (x1 <= 0.437500000000000000e+00_dp) THEN
481 tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
482 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
483 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 68))
484 ELSE
485 tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
486 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
487 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 69))
488 END IF
489 ELSE
490 IF (x1 <= 0.437500000000000000e+00_dp) THEN
491 tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
492 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
493 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 70))
494 ELSE
495 tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
496 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
497 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 71))
498 END IF
499 END IF
500 ELSE
501 IF (x1 <= 0.437500000000000000e+00_dp) THEN
502 tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
503 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
504 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 72))
505 ELSE
506 tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
507 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
508 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 73))
509 END IF
510 END IF
511 END IF
512 ELSE
513 IF (x2 <= 0.750000000000000000e+00_dp) THEN
514 IF (x2 <= 0.625000000000000000e+00_dp) THEN
515 IF (x1 <= 0.775000000000000022e+00_dp) THEN
516 tg1 = (2*x1 - 0.132500000000000018e+01_dp)*0.444444444444444464e+01_dp
517 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
518 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 74))
519 ELSE
520 tg1 = (2*x1 - 0.177499999999999991e+01_dp)*0.444444444444444464e+01_dp
521 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
522 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 75))
523 END IF
524 ELSE
525 IF (x1 <= 0.775000000000000022e+00_dp) THEN
526 tg1 = (2*x1 - 0.132500000000000018e+01_dp)*0.444444444444444464e+01_dp
527 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
528 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 76))
529 ELSE
530 tg1 = (2*x1 - 0.177499999999999991e+01_dp)*0.444444444444444464e+01_dp
531 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
532 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 77))
533 END IF
534 END IF
535 ELSE
536 IF (x2 <= 0.875000000000000000e+00_dp) THEN
537 tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
538 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
539 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 78))
540 ELSE
541 tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
542 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
543 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 79))
544 END IF
545 END IF
546 END IF
547 END IF
548 ELSE
549 x1 = sqrt(81.0_dp/t)
550 lower = (sqrt(t) - sqrt(eps))/n
551 IF (r >= lower) THEN
552 x2 = lower/r
553 IF (x2 <= 0.500000000000000000e+00_dp) THEN
554 IF (x1 <= 0.500000000000000000e+00_dp) THEN
555 IF (x2 <= 0.250000000000000000e+00_dp) THEN
556 tg1 = (2*x1 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
557 tg2 = (2*x2 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
558 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 80))
559 ELSE
560 IF (x2 <= 0.375000000000000000e+00_dp) THEN
561 IF (x1 <= 0.250000000000000000e+00_dp) THEN
562 tg1 = (2*x1 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
563 tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
564 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 81))
565 ELSE
566 IF (x2 <= 0.312500000000000000e+00_dp) THEN
567 tg1 = (2*x1 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
568 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
569 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 82))
570 ELSE
571 IF (x1 <= 0.375000000000000000e+00_dp) THEN
572 tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
573 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
574 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 83))
575 ELSE
576 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
577 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
578 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 84))
579 END IF
580 END IF
581 END IF
582 ELSE
583 IF (x1 <= 0.250000000000000000e+00_dp) THEN
584 IF (x2 <= 0.437500000000000000e+00_dp) THEN
585 IF (x1 <= 0.125000000000000000e+00_dp) THEN
586 tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
587 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
588 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 85))
589 ELSE
590 tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
591 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
592 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 86))
593 END IF
594 ELSE
595 IF (x1 <= 0.125000000000000000e+00_dp) THEN
596 IF (x2 <= 0.468750000000000000e+00_dp) THEN
597 tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
598 tg2 = (2*x2 - 0.906250000000000000e+00_dp)*0.320000000000000000e+02_dp
599 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 87))
600 ELSE
601 tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
602 tg2 = (2*x2 - 0.968750000000000000e+00_dp)*0.320000000000000000e+02_dp
603 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 88))
604 END IF
605 ELSE
606 IF (x2 <= 0.468750000000000000e+00_dp) THEN
607 tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
608 tg2 = (2*x2 - 0.906250000000000000e+00_dp)*0.320000000000000000e+02_dp
609 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 89))
610 ELSE
611 tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
612 tg2 = (2*x2 - 0.968750000000000000e+00_dp)*0.320000000000000000e+02_dp
613 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 90))
614 END IF
615 END IF
616 END IF
617 ELSE
618 IF (x2 <= 0.437500000000000000e+00_dp) THEN
619 IF (x2 <= 0.406250000000000000e+00_dp) THEN
620 IF (x1 <= 0.375000000000000000e+00_dp) THEN
621 tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
622 tg2 = (2*x2 - 0.781250000000000000e+00_dp)*0.320000000000000000e+02_dp
623 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 91))
624 ELSE
625 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
626 tg2 = (2*x2 - 0.781250000000000000e+00_dp)*0.320000000000000000e+02_dp
627 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 92))
628 END IF
629 ELSE
630 tg1 = (2*x1 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
631 tg2 = (2*x2 - 0.843750000000000000e+00_dp)*0.320000000000000000e+02_dp
632 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 93))
633 END IF
634 ELSE
635 tg1 = (2*x1 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
636 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
637 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 94))
638 END IF
639 END IF
640 END IF
641 END IF
642 ELSE
643 IF (x2 <= 0.250000000000000000e+00_dp) THEN
644 IF (x2 <= 0.125000000000000000e+00_dp) THEN
645 tg1 = (2*x1 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
646 tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
647 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 95))
648 ELSE
649 IF (x2 <= 0.187500000000000000e+00_dp) THEN
650 IF (x1 <= 0.750000000000000000e+00_dp) THEN
651 tg1 = (2*x1 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
652 tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
653 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 96))
654 ELSE
655 IF (x1 <= 0.875000000000000000e+00_dp) THEN
656 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
657 tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
658 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 97))
659 ELSE
660 IF (x2 <= 0.156250000000000000e+00_dp) THEN
661 tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
662 tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
663 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 98))
664 ELSE
665 tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
666 tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
667 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 99))
668 END IF
669 END IF
670 END IF
671 ELSE
672 IF (x1 <= 0.750000000000000000e+00_dp) THEN
673 IF (x1 <= 0.625000000000000000e+00_dp) THEN
674 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
675 tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
676 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 100))
677 ELSE
678 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
679 tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
680 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 101))
681 END IF
682 ELSE
683 IF (x1 <= 0.875000000000000000e+00_dp) THEN
684 IF (x2 <= 0.218750000000000000e+00_dp) THEN
685 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
686 tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
687 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 102))
688 ELSE
689 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
690 tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
691 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 103))
692 END IF
693 ELSE
694 IF (x2 <= 0.218750000000000000e+00_dp) THEN
695 tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
696 tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
697 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 104))
698 ELSE
699 tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
700 tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
701 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 IF (x2 <= 0.375000000000000000e+00_dp) THEN
709 IF (x1 <= 0.750000000000000000e+00_dp) THEN
710 IF (x2 <= 0.312500000000000000e+00_dp) THEN
711 IF (x1 <= 0.625000000000000000e+00_dp) THEN
712 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
713 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
714 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 106))
715 ELSE
716 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
717 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
718 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 107))
719 END IF
720 ELSE
721 IF (x1 <= 0.625000000000000000e+00_dp) THEN
722 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
723 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
724 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 108))
725 ELSE
726 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
727 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
728 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 109))
729 END IF
730 END IF
731 ELSE
732 IF (x2 <= 0.312500000000000000e+00_dp) THEN
733 IF (x1 <= 0.875000000000000000e+00_dp) THEN
734 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
735 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
736 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 110))
737 ELSE
738 tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
739 tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
740 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 111))
741 END IF
742 ELSE
743 IF (x1 <= 0.875000000000000000e+00_dp) THEN
744 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
745 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
746 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 112))
747 ELSE
748 IF (x1 <= 0.937500000000000000e+00_dp) THEN
749 tg1 = (2*x1 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
750 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
751 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 113))
752 ELSE
753 tg1 = (2*x1 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
754 tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
755 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 114))
756 END IF
757 END IF
758 END IF
759 END IF
760 ELSE
761 IF (x1 <= 0.750000000000000000e+00_dp) THEN
762 IF (x1 <= 0.625000000000000000e+00_dp) THEN
763 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
764 tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
765 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 115))
766 ELSE
767 IF (x2 <= 0.437500000000000000e+00_dp) THEN
768 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
769 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
770 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 116))
771 ELSE
772 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
773 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
774 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 117))
775 END IF
776 END IF
777 ELSE
778 IF (x1 <= 0.875000000000000000e+00_dp) THEN
779 IF (x2 <= 0.437500000000000000e+00_dp) THEN
780 IF (x1 <= 0.812500000000000000e+00_dp) THEN
781 tg1 = (2*x1 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
782 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
783 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 118))
784 ELSE
785 tg1 = (2*x1 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
786 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
787 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 119))
788 END IF
789 ELSE
790 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
791 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
792 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 120))
793 END IF
794 ELSE
795 IF (x2 <= 0.437500000000000000e+00_dp) THEN
796 IF (x1 <= 0.937500000000000000e+00_dp) THEN
797 tg1 = (2*x1 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
798 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
799 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 121))
800 ELSE
801 tg1 = (2*x1 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
802 tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
803 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 122))
804 END IF
805 ELSE
806 tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
807 tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
808 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 IF (x1 <= 0.500000000000000000e+00_dp) THEN
817 IF (x2 <= 0.750000000000000000e+00_dp) THEN
818 IF (x1 <= 0.250000000000000000e+00_dp) THEN
819 tg1 = (2*x1 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
820 tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
821 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 124))
822 ELSE
823 IF (x2 <= 0.625000000000000000e+00_dp) THEN
824 IF (x1 <= 0.375000000000000000e+00_dp) THEN
825 tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
826 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
827 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 125))
828 ELSE
829 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
830 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
831 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 126))
832 END IF
833 ELSE
834 IF (x1 <= 0.375000000000000000e+00_dp) THEN
835 IF (x2 <= 0.687500000000000000e+00_dp) THEN
836 tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
837 tg2 = (2*x2 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
838 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 127))
839 ELSE
840 tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
841 tg2 = (2*x2 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
842 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 128))
843 END IF
844 ELSE
845 IF (x2 <= 0.687500000000000000e+00_dp) THEN
846 IF (x1 <= 0.437500000000000000e+00_dp) THEN
847 tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
848 tg2 = (2*x2 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
849 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 129))
850 ELSE
851 tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
852 tg2 = (2*x2 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
853 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 130))
854 END IF
855 ELSE
856 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
857 tg2 = (2*x2 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
858 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 131))
859 END IF
860 END IF
861 END IF
862 END IF
863 ELSE
864 IF (x1 <= 0.250000000000000000e+00_dp) THEN
865 IF (x2 <= 0.875000000000000000e+00_dp) THEN
866 IF (x1 <= 0.125000000000000000e+00_dp) THEN
867 tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
868 tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
869 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 132))
870 ELSE
871 IF (x2 <= 0.812500000000000000e+00_dp) THEN
872 tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
873 tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
874 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 133))
875 ELSE
876 IF (x1 <= 0.187500000000000000e+00_dp) THEN
877 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
878 tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
879 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 134))
880 ELSE
881 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
882 tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
883 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 135))
884 END IF
885 END IF
886 END IF
887 ELSE
888 IF (x1 <= 0.125000000000000000e+00_dp) THEN
889 IF (x2 <= 0.937500000000000000e+00_dp) THEN
890 IF (x1 <= 0.625000000000000000e-01_dp) THEN
891 tg1 = (2*x1 - 0.625000000000000000e-01_dp)*0.160000000000000000e+02_dp
892 tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
893 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 136))
894 ELSE
895 IF (x2 <= 0.906250000000000000e+00_dp) THEN
896 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
897 tg2 = (2*x2 - 0.178125000000000000e+01_dp)*0.320000000000000000e+02_dp
898 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 137))
899 ELSE
900 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
901 tg2 = (2*x2 - 0.184375000000000000e+01_dp)*0.320000000000000000e+02_dp
902 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 138))
903 END IF
904 END IF
905 ELSE
906 IF (x1 <= 0.625000000000000000e-01_dp) THEN
907 IF (x2 <= 0.968750000000000000e+00_dp) THEN
908 IF (x1 <= 0.312500000000000000e-01_dp) THEN
909 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
910 tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
911 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 139))
912 ELSE
913 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
914 tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
915 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 140))
916 END IF
917 ELSE
918 IF (x1 <= 0.312500000000000000e-01_dp) THEN
919 IF (x2 <= 0.984375000000000000e+00_dp) THEN
920 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
921 tg2 = (2*x2 - 0.195312500000000000e+01_dp)*0.640000000000000000e+02_dp
922 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 141))
923 ELSE
924 tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
925 tg2 = (2*x2 - 0.198437500000000000e+01_dp)*0.640000000000000000e+02_dp
926 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 142))
927 END IF
928 ELSE
929 tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
930 tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
931 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 143))
932 END IF
933 END IF
934 ELSE
935 IF (x2 <= 0.968750000000000000e+00_dp) THEN
936 IF (x1 <= 0.937500000000000000e-01_dp) THEN
937 tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
938 tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
939 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 144))
940 ELSE
941 tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
942 tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
943 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 145))
944 END IF
945 ELSE
946 tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
947 tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
948 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 146))
949 END IF
950 END IF
951 END IF
952 ELSE
953 IF (x2 <= 0.937500000000000000e+00_dp) THEN
954 IF (x1 <= 0.187500000000000000e+00_dp) THEN
955 IF (x2 <= 0.906250000000000000e+00_dp) THEN
956 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
957 tg2 = (2*x2 - 0.178125000000000000e+01_dp)*0.320000000000000000e+02_dp
958 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 147))
959 ELSE
960 tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
961 tg2 = (2*x2 - 0.184375000000000000e+01_dp)*0.320000000000000000e+02_dp
962 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 148))
963 END IF
964 ELSE
965 tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
966 tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
967 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 149))
968 END IF
969 ELSE
970 tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
971 tg2 = (2*x2 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
972 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 150))
973 END IF
974 END IF
975 END IF
976 ELSE
977 IF (x2 <= 0.875000000000000000e+00_dp) THEN
978 IF (x1 <= 0.375000000000000000e+00_dp) THEN
979 IF (x2 <= 0.812500000000000000e+00_dp) THEN
980 IF (x1 <= 0.312500000000000000e+00_dp) THEN
981 tg1 = (2*x1 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
982 tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
983 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 151))
984 ELSE
985 tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
986 tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
987 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 152))
988 END IF
989 ELSE
990 tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
991 tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
992 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 153))
993 END IF
994 ELSE
995 IF (x2 <= 0.812500000000000000e+00_dp) THEN
996 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
997 tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
998 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 154))
999 ELSE
1000 tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
1001 tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
1002 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 155))
1003 END IF
1004 END IF
1005 ELSE
1006 tg1 = (2*x1 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
1007 tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1008 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 156))
1009 END IF
1010 END IF
1011 END IF
1012 ELSE
1013 IF (x2 <= 0.750000000000000000e+00_dp) THEN
1014 IF (x1 <= 0.750000000000000000e+00_dp) THEN
1015 IF (x2 <= 0.625000000000000000e+00_dp) THEN
1016 IF (x1 <= 0.625000000000000000e+00_dp) THEN
1017 IF (x2 <= 0.562500000000000000e+00_dp) THEN
1018 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1019 tg2 = (2*x2 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1020 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 157))
1021 ELSE
1022 IF (x1 <= 0.562500000000000000e+00_dp) THEN
1023 tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1024 tg2 = (2*x2 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1025 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 158))
1026 ELSE
1027 tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1028 tg2 = (2*x2 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1029 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 159))
1030 END IF
1031 END IF
1032 ELSE
1033 IF (x2 <= 0.562500000000000000e+00_dp) THEN
1034 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1035 tg2 = (2*x2 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1036 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 160))
1037 ELSE
1038 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1039 tg2 = (2*x2 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1040 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 161))
1041 END IF
1042 END IF
1043 ELSE
1044 IF (x1 <= 0.625000000000000000e+00_dp) THEN
1045 IF (x2 <= 0.687500000000000000e+00_dp) THEN
1046 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1047 tg2 = (2*x2 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
1048 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 162))
1049 ELSE
1050 tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1051 tg2 = (2*x2 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
1052 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 163))
1053 END IF
1054 ELSE
1055 tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1056 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1057 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 164))
1058 END IF
1059 END IF
1060 ELSE
1061 IF (x2 <= 0.625000000000000000e+00_dp) THEN
1062 IF (x1 <= 0.875000000000000000e+00_dp) THEN
1063 IF (x2 <= 0.562500000000000000e+00_dp) THEN
1064 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1065 tg2 = (2*x2 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1066 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 165))
1067 ELSE
1068 tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1069 tg2 = (2*x2 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1070 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 166))
1071 END IF
1072 ELSE
1073 tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1074 tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1075 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 167))
1076 END IF
1077 ELSE
1078 tg1 = (2*x1 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1079 tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1080 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 168))
1081 END IF
1082 END IF
1083 ELSE
1084 IF (x1 <= 0.750000000000000000e+00_dp) THEN
1085 tg1 = (2*x1 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
1086 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1087 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 169))
1088 ELSE
1089 tg1 = (2*x1 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1090 tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1091 CALL pd2val(res, nderiv, tg1, tg2, c0(1, 170))
1092 END IF
1093 END IF
1094 END IF
1095 END IF
1096 ELSE
1097 DO i = 1, nderiv + 1
1098 res(i) = 0.0_dp
1099 END DO
1100 END IF
1101 END IF
1102 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1118
1119 patches = 170
1120 IF (nder > nderiv_max) cpabort("Reading data for initialization of C0 failed")
1121 nderiv_init = nder
1122 IF (ALLOCATED(c0)) DEALLOCATE (c0)
1123 ! round up to a multiple of 32 to give some generous alignment for each C0
1124 ALLOCATE (c0(32*((31 + (nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1125 ! valgrind workaround
1126 c0 = huge(0.0_dp)
1127 IF (mepos == 0) THEN
1128 ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1129 DO i = 1, patches
1130 READ (iunit, *) chunk
1131 c0(1:(nder + 1)*(degree + 1)*(degree + 2)/2, i) = chunk(1:(nder + 1)*(degree + 1)*(degree + 2)/2)
1132 END DO
1133 DEALLOCATE (chunk)
1134 END IF
1135 CALL group%bcast(c0, 0)
1136 END SUBROUTINE init
1137
1138! **************************************************************************************************
1139!> \brief ...
1140! **************************************************************************************************
1141 SUBROUTINE free_c0()
1142 IF (ALLOCATED(c0)) DEALLOCATE (c0)
1143 nderiv_init = -1
1144 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 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 t1(0) = 1.0_dp
1165 t2(0) = 1.0_dp
1166 t1(1) = sqrt2*tg1
1167 t2(1) = sqrt2*tg2
1168 t1(2) = 2*tg1*t1(1) - sqrt2
1169 t2(2) = 2*tg2*t2(1) - sqrt2
1170 t1(3) = 2*tg1*t1(2) - t1(1)
1171 t2(3) = 2*tg2*t2(2) - t2(1)
1172 t1(4) = 2*tg1*t1(3) - t1(2)
1173 t2(4) = 2*tg2*t2(3) - t2(2)
1174 t1(5) = 2*tg1*t1(4) - t1(3)
1175 t2(5) = 2*tg2*t2(4) - t2(3)
1176 t1(6) = 2*tg1*t1(5) - t1(4)
1177 t2(6) = 2*tg2*t2(5) - t2(4)
1178 t1(7) = 2*tg1*t1(6) - t1(5)
1179 t2(7) = 2*tg2*t2(6) - t2(5)
1180 t1(8) = 2*tg1*t1(7) - t1(6)
1181 t2(8) = 2*tg2*t2(7) - t2(6)
1182 t1(9) = 2*tg1*t1(8) - t1(7)
1183 t2(9) = 2*tg2*t2(8) - t2(7)
1184 t1(10) = 2*tg1*t1(9) - t1(8)
1185 t2(10) = 2*tg2*t2(9) - t2(8)
1186 t1(11) = 2*tg1*t1(10) - t1(9)
1187 t2(11) = 2*tg2*t2(10) - t2(9)
1188 t1(12) = 2*tg1*t1(11) - t1(10)
1189 t2(12) = 2*tg2*t2(11) - t2(10)
1190 t1(13) = 2*tg1*t1(12) - t1(11)
1191 t2(13) = 2*tg2*t2(12) - t2(11)
1192 t1(14) = 2*tg1*t1(13) - t1(12)
1193 t2(14) = 2*tg2*t2(13) - t2(12)
1194 t1(15) = 2*tg1*t1(14) - t1(13)
1195 t2(15) = 2*tg2*t2(14) - t2(13)
1196 DO k = 1, nderiv + 1
1197 res(k) = 0.0_dp
1198 res(k) = res(k) + dot_product(t1(0:15), c0(1:16, k))*t2(0)
1199 res(k) = res(k) + dot_product(t1(0:14), c0(17:31, k))*t2(1)
1200 res(k) = res(k) + dot_product(t1(0:13), c0(32:45, k))*t2(2)
1201 res(k) = res(k) + dot_product(t1(0:12), c0(46:58, k))*t2(3)
1202 res(k) = res(k) + dot_product(t1(0:11), c0(59:70, k))*t2(4)
1203 res(k) = res(k) + dot_product(t1(0:10), c0(71:81, k))*t2(5)
1204 res(k) = res(k) + dot_product(t1(0:9), c0(82:91, k))*t2(6)
1205 res(k) = res(k) + dot_product(t1(0:8), c0(92:100, k))*t2(7)
1206 res(k) = res(k) + dot_product(t1(0:7), c0(101:108, k))*t2(8)
1207 res(k) = res(k) + dot_product(t1(0:6), c0(109:115, k))*t2(9)
1208 res(k) = res(k) + dot_product(t1(0:5), c0(116:121, k))*t2(10)
1209 res(k) = res(k) + dot_product(t1(0:4), c0(122:126, k))*t2(11)
1210 res(k) = res(k) + dot_product(t1(0:3), c0(127:130, k))*t2(12)
1211 res(k) = res(k) + dot_product(t1(0:2), c0(131:133, k))*t2(13)
1212 res(k) = res(k) + dot_product(t1(0:1), c0(134:135, k))*t2(14)
1213 res(k) = res(k) + dot_product(t1(0:0), c0(136:136, k))*t2(15)
1214 END DO
1215 END SUBROUTINE pd2val
1216
1217END MODULE t_sh_p_s_c
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
t_sh_p_s_c
...
Definition t_sh_p_s_c.F:39
t_sh_p_s_c::trunc_cs_poly_n20
subroutine, public trunc_cs_poly_n20(res, r, t, nderiv)
...
Definition t_sh_p_s_c.F:65
t_sh_p_s_c::free_c0
subroutine, public free_c0()
...
Definition t_sh_p_s_c.F:1142
t_sh_p_s_c::init
subroutine, public init(nder, iunit, mepos, group)
...
Definition t_sh_p_s_c.F:1112
message_passing::mp_comm_type
Definition message_passing.F:189