Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : !--------------------------------------------------------------------------------------------------!
9 : ! Copyright (c) 2008, 2009, 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 This module computes the basic integrals for the truncated coulomb operator
34 : !>
35 : !> res(1) =G_0(R,T)= ((2*erf(sqrt(t))+erf(R-sqrt(t))-erf(R+sqrt(t)))/sqrt(t))
36 : !>
37 : !> and up to 21 derivatives with respect to T
38 : !>
39 : !> res(n+1)=(-1)**n d^n/dT^n G_0(R,T)
40 : !>
41 : !> The function is only computed for values of R,T which fulfil
42 : !>
43 : !> R**2 - 11.0_dp*R + 0.0_dp < T < R**2 + 11.0_dp*R + 50.0_dp where R>=0 T>=0
44 : !>
45 : !> for T larger than the upper bound, 0 is returned
46 : !> (which is accurate at least up to 1.0E-16)
47 : !> while for T smaller than the lower bound, the caller is instructed
48 : !> to use the conventional gamma function instead
49 : !> (i.e. the limit of above expression for R to Infinity)
50 : !>
51 : !> \author Joost VandeVondele and Manuel Guidon
52 : !> \par History
53 : !> Nov 2008, 2009 Joost VandeVondele and Manuel Guidon
54 : !> May 2019 A. Bussy: Added a get_maxl_init function to get current status of nderiv_init and
55 : !> moved the file to common (made it accessible from aobasis, same place as gamma.F).
56 : ! **************************************************************************************************
57 : MODULE t_c_g0
58 : USE kinds, ONLY: dp
59 : USE message_passing, ONLY: mp_comm_type
60 : #include "../base/base_uses.f90"
61 :
62 : IMPLICIT NONE
63 : PRIVATE
64 :
65 : PUBLIC :: t_c_g0_n, init, free_C0, get_lmax_init
66 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE :: C0
67 : INTEGER, PARAMETER :: degree = 13
68 : REAL(KIND=dp), PARAMETER :: target_error = 0.100000E-08
69 : INTEGER, PARAMETER :: nderiv_max = 21
70 : INTEGER, SAVE :: nderiv_init = -1
71 : INTEGER, SAVE :: patches = -1
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief ...
77 : !> \param RES ...
78 : !> \param use_gamma ...
79 : !> \param R ...
80 : !> \param T ...
81 : !> \param NDERIV ...
82 : ! **************************************************************************************************
83 249786653 : SUBROUTINE t_c_g0_n(RES, use_gamma, R, T, NDERIV)
84 : REAL(KIND=dp), INTENT(OUT) :: RES(*)
85 : LOGICAL, INTENT(OUT) :: use_gamma
86 : REAL(KIND=dp), INTENT(IN) :: R, T
87 : INTEGER, INTENT(IN) :: NDERIV
88 :
89 : REAL(KIND=dp) :: lower, TG1, TG2, upper, X1, X2
90 :
91 249786653 : use_gamma = .FALSE.
92 249786653 : upper = R**2 + 11.0_dp*R + 50.0_dp
93 249786653 : lower = R**2 - 11.0_dp*R + 0.0_dp
94 249786653 : IF (T > upper) THEN
95 7378830 : RES(1:NDERIV + 1) = 0.0_dp
96 7357441 : RETURN
97 : END IF
98 247840700 : IF (R <= 11.0_dp) THEN
99 210080396 : X2 = R/11.0_dp
100 210080396 : upper = R**2 + 11.0_dp*R + 50.0_dp
101 210080396 : lower = 0.0_dp
102 210080396 : X1 = (T - lower)/(upper - lower)
103 210080396 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
104 138047706 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
105 115975667 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
106 95324277 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
107 20274012 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
108 9430346 : IF (X2 <= 0.625000000000000000E-01_dp) THEN
109 1190480 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
110 612287 : IF (X2 <= 0.312500000000000000E-01_dp) THEN
111 41596 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
112 16419 : IF (X2 <= 0.156250000000000000E-01_dp) THEN
113 0 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
114 0 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
115 0 : TG2 = (2*X2 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
116 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 1))
117 : ELSE
118 0 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
119 0 : TG2 = (2*X2 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
120 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 2))
121 : END IF
122 : ELSE
123 16419 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
124 8065 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
125 8065 : TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
126 8065 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 3))
127 : ELSE
128 8354 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
129 8354 : TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
130 8354 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 4))
131 : END IF
132 : END IF
133 : ELSE
134 25177 : IF (X2 <= 0.156250000000000000E-01_dp) THEN
135 0 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
136 0 : TG2 = (2*X2 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
137 0 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 5))
138 : ELSE
139 25177 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
140 25177 : TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
141 25177 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 6))
142 : END IF
143 : END IF
144 : ELSE
145 570691 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
146 302127 : IF (X2 <= 0.468750000000000000E-01_dp) THEN
147 132338 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
148 67988 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
149 67988 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
150 67988 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 7))
151 : ELSE
152 64350 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
153 64350 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
154 64350 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 8))
155 : END IF
156 : ELSE
157 169789 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
158 129960 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
159 129960 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
160 129960 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 9))
161 : ELSE
162 39829 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
163 39829 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
164 39829 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 10))
165 : END IF
166 : END IF
167 : ELSE
168 268564 : IF (X2 <= 0.468750000000000000E-01_dp) THEN
169 205389 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
170 205389 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
171 205389 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 11))
172 : ELSE
173 63175 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
174 63175 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
175 63175 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 12))
176 : END IF
177 : END IF
178 : END IF
179 : ELSE
180 578193 : IF (X2 <= 0.312500000000000000E-01_dp) THEN
181 50972 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
182 16180 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
183 16180 : TG2 = (2*X2 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
184 16180 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 13))
185 : ELSE
186 34792 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
187 34792 : TG2 = (2*X2 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
188 34792 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 14))
189 : END IF
190 : ELSE
191 527221 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
192 283055 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
193 283055 : TG2 = (2*X2 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
194 283055 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 15))
195 : ELSE
196 244166 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
197 244166 : TG2 = (2*X2 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
198 244166 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 16))
199 : END IF
200 : END IF
201 : END IF
202 : ELSE
203 8239866 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
204 3980279 : IF (X2 <= 0.937500000000000000E-01_dp) THEN
205 928803 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
206 599975 : IF (X2 <= 0.781250000000000000E-01_dp) THEN
207 272666 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
208 239540 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
209 239540 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
210 239540 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 17))
211 : ELSE
212 33126 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
213 33126 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
214 33126 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 18))
215 : END IF
216 : ELSE
217 327309 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
218 237122 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
219 237122 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
220 237122 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 19))
221 : ELSE
222 90187 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
223 90187 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
224 90187 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 20))
225 : END IF
226 : END IF
227 : ELSE
228 328828 : IF (X2 <= 0.781250000000000000E-01_dp) THEN
229 65568 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
230 65568 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
231 65568 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 21))
232 : ELSE
233 263260 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
234 263260 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
235 263260 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 22))
236 : END IF
237 : END IF
238 : ELSE
239 3051476 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
240 1566263 : IF (X2 <= 0.109375000000000000E+00_dp) THEN
241 724901 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
242 516968 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
243 516968 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
244 516968 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 23))
245 : ELSE
246 207933 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
247 207933 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
248 207933 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 24))
249 : END IF
250 : ELSE
251 841362 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
252 609209 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
253 609209 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
254 609209 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 25))
255 : ELSE
256 232153 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
257 232153 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
258 232153 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 26))
259 : END IF
260 : END IF
261 : ELSE
262 1485213 : IF (X2 <= 0.109375000000000000E+00_dp) THEN
263 614510 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
264 614510 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
265 614510 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 27))
266 : ELSE
267 870703 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
268 870703 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
269 870703 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 28))
270 : END IF
271 : END IF
272 : END IF
273 : ELSE
274 4259587 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
275 2052883 : IF (X2 <= 0.937500000000000000E-01_dp) THEN
276 342734 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
277 342734 : TG2 = (2*X2 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
278 342734 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 29))
279 : ELSE
280 1710149 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
281 1710149 : TG2 = (2*X2 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
282 1710149 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 30))
283 : END IF
284 : ELSE
285 2206704 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
286 2206704 : TG2 = (2*X2 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
287 2206704 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 31))
288 : END IF
289 : END IF
290 : END IF
291 : ELSE
292 10843666 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
293 5770828 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
294 5770828 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
295 5770828 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 32))
296 : ELSE
297 5072838 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
298 5072838 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
299 5072838 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 33))
300 : END IF
301 : END IF
302 : ELSE
303 75050265 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
304 32953185 : IF (X2 <= 0.187500000000000000E+00_dp) THEN
305 20809021 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
306 10152876 : IF (X2 <= 0.156250000000000000E+00_dp) THEN
307 4686369 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
308 2180124 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
309 1411083 : IF (X2 <= 0.140625000000000000E+00_dp) THEN
310 620094 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
311 620094 : TG2 = (2*X2 - 0.265625000000000000E+00_dp)*0.640000000000000000E+02_dp
312 620094 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 34))
313 : ELSE
314 790989 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
315 790989 : TG2 = (2*X2 - 0.296875000000000000E+00_dp)*0.640000000000000000E+02_dp
316 790989 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 35))
317 : END IF
318 : ELSE
319 769041 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
320 769041 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
321 769041 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 36))
322 : END IF
323 : ELSE
324 2506245 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
325 1124888 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
326 1124888 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
327 1124888 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 37))
328 : ELSE
329 1381357 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
330 1381357 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
331 1381357 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 38))
332 : END IF
333 : END IF
334 : ELSE
335 5466507 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
336 3238304 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
337 2172443 : IF (X2 <= 0.171875000000000000E+00_dp) THEN
338 898189 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
339 898189 : TG2 = (2*X2 - 0.328125000000000000E+00_dp)*0.640000000000000000E+02_dp
340 898189 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 39))
341 : ELSE
342 1274254 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
343 1274254 : TG2 = (2*X2 - 0.359375000000000000E+00_dp)*0.640000000000000000E+02_dp
344 1274254 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 40))
345 : END IF
346 : ELSE
347 1065861 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
348 1065861 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
349 1065861 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 41))
350 : END IF
351 : ELSE
352 2228203 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
353 2228203 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
354 2228203 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 42))
355 : END IF
356 : END IF
357 : ELSE
358 10656145 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
359 5113610 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
360 5113610 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
361 5113610 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 43))
362 : ELSE
363 5542535 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
364 5542535 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
365 5542535 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 44))
366 : END IF
367 : END IF
368 : ELSE
369 12144164 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
370 7941251 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
371 5838342 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
372 3081615 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
373 2119664 : IF (X2 <= 0.203125000000000000E+00_dp) THEN
374 1097360 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
375 1097360 : TG2 = (2*X2 - 0.390625000000000000E+00_dp)*0.640000000000000000E+02_dp
376 1097360 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 45))
377 : ELSE
378 1022304 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
379 1022304 : TG2 = (2*X2 - 0.421875000000000000E+00_dp)*0.640000000000000000E+02_dp
380 1022304 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 46))
381 : END IF
382 : ELSE
383 961951 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
384 961951 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
385 961951 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 47))
386 : END IF
387 : ELSE
388 2756727 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
389 2048130 : IF (X2 <= 0.234375000000000000E+00_dp) THEN
390 970679 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
391 970679 : TG2 = (2*X2 - 0.453125000000000000E+00_dp)*0.640000000000000000E+02_dp
392 970679 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 48))
393 : ELSE
394 1077451 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
395 1077451 : TG2 = (2*X2 - 0.484375000000000000E+00_dp)*0.640000000000000000E+02_dp
396 1077451 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 49))
397 : END IF
398 : ELSE
399 708597 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
400 708597 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
401 708597 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 50))
402 : END IF
403 : END IF
404 : ELSE
405 2102909 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
406 1270035 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
407 1270035 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
408 1270035 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 51))
409 : ELSE
410 832874 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
411 832874 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
412 832874 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 52))
413 : END IF
414 : END IF
415 : ELSE
416 4202913 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
417 2028360 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
418 2028360 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
419 2028360 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 53))
420 : ELSE
421 2174553 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
422 2174553 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
423 2174553 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 54))
424 : END IF
425 : END IF
426 : END IF
427 : ELSE
428 42097080 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
429 18879343 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
430 8623538 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
431 8623538 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
432 8623538 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 55))
433 : ELSE
434 10255805 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
435 10255805 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
436 10255805 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 56))
437 : END IF
438 : ELSE
439 23217737 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
440 23217737 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
441 23217737 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 57))
442 : END IF
443 : END IF
444 : END IF
445 : ELSE
446 20651390 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
447 15264575 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
448 13549666 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
449 11548864 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
450 7622412 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
451 4572674 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
452 3568686 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
453 1877856 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
454 1877856 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
455 1877856 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 58))
456 : ELSE
457 1690830 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
458 1690830 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
459 1690830 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 59))
460 : END IF
461 : ELSE
462 1003988 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
463 550057 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
464 550057 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
465 550057 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 60))
466 : ELSE
467 453931 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
468 453931 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
469 453931 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 61))
470 : END IF
471 : END IF
472 : ELSE
473 3049738 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
474 2453733 : IF (X2 <= 0.343750000000000000E+00_dp) THEN
475 1332666 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
476 1332666 : TG2 = (2*X2 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
477 1332666 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 62))
478 : ELSE
479 1121067 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
480 1121067 : TG2 = (2*X2 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
481 1121067 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 63))
482 : END IF
483 : ELSE
484 596005 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
485 596005 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
486 596005 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 64))
487 : END IF
488 : END IF
489 : ELSE
490 3926452 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
491 3432192 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
492 2192080 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
493 1944762 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
494 1944762 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
495 1944762 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 65))
496 : ELSE
497 247318 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
498 247318 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
499 247318 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 66))
500 : END IF
501 : ELSE
502 1240112 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
503 1124595 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
504 1124595 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
505 1124595 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 67))
506 : ELSE
507 115517 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
508 115517 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
509 115517 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 68))
510 : END IF
511 : END IF
512 : ELSE
513 494260 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
514 316190 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
515 316190 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
516 316190 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 69))
517 : ELSE
518 178070 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
519 178070 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
520 178070 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 70))
521 : END IF
522 : END IF
523 : END IF
524 : ELSE
525 2000802 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
526 1565095 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
527 1054734 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
528 683651 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
529 683651 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
530 683651 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 71))
531 : ELSE
532 371083 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
533 371083 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
534 371083 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 72))
535 : END IF
536 : ELSE
537 510361 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
538 291841 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
539 291841 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
540 291841 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 73))
541 : ELSE
542 218520 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
543 218520 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
544 218520 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 74))
545 : END IF
546 : END IF
547 : ELSE
548 435707 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
549 238150 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
550 157269 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
551 157269 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
552 157269 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 75))
553 : ELSE
554 80881 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
555 80881 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
556 80881 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 76))
557 : END IF
558 : ELSE
559 197557 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
560 128122 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
561 128122 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
562 128122 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 77))
563 : ELSE
564 69435 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
565 69435 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
566 69435 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 78))
567 : END IF
568 : END IF
569 : END IF
570 : END IF
571 : ELSE
572 1714909 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
573 1465938 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
574 608378 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
575 362285 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
576 362285 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
577 362285 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 79))
578 : ELSE
579 246093 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
580 171580 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
581 171580 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
582 171580 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 80))
583 : ELSE
584 74513 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
585 74513 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
586 74513 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 81))
587 : END IF
588 : END IF
589 : ELSE
590 857560 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
591 742012 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
592 742012 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
593 742012 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 82))
594 : ELSE
595 115548 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
596 115548 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
597 115548 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 83))
598 : END IF
599 : END IF
600 : ELSE
601 248971 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
602 183247 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
603 112073 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
604 112073 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
605 112073 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 84))
606 : ELSE
607 71174 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
608 44470 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
609 44470 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
610 44470 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 85))
611 : ELSE
612 26704 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
613 26704 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
614 26704 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 86))
615 : END IF
616 : END IF
617 : ELSE
618 65724 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
619 40500 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
620 40500 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
621 40500 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 87))
622 : ELSE
623 25224 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
624 25224 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
625 25224 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 88))
626 : END IF
627 : END IF
628 : END IF
629 : END IF
630 : ELSE
631 5386815 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
632 1640244 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
633 1507908 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
634 701089 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
635 701089 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
636 701089 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 89))
637 : ELSE
638 806819 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
639 806819 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
640 806819 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 90))
641 : END IF
642 : ELSE
643 132336 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
644 64727 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
645 30990 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
646 30990 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
647 30990 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 91))
648 : ELSE
649 33737 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
650 33737 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
651 33737 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 92))
652 : END IF
653 : ELSE
654 67609 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
655 67609 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
656 67609 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 93))
657 : END IF
658 : END IF
659 : ELSE
660 3746571 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
661 1595505 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
662 1595505 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
663 1595505 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 94))
664 : ELSE
665 2151066 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
666 2151066 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
667 2151066 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 95))
668 : END IF
669 : END IF
670 : END IF
671 : END IF
672 : ELSE
673 22072039 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
674 12461026 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
675 8833250 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
676 7211077 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
677 6469562 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
678 6076052 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
679 5854355 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
680 4375324 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
681 4375324 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
682 4375324 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 96))
683 : ELSE
684 1479031 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
685 1479031 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
686 1479031 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 97))
687 : END IF
688 : ELSE
689 221697 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
690 128355 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
691 128355 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
692 128355 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 98))
693 : ELSE
694 93342 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
695 93342 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
696 93342 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 99))
697 : END IF
698 : END IF
699 : ELSE
700 393510 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
701 203467 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
702 203467 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
703 203467 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 100))
704 : ELSE
705 190043 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
706 190043 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
707 190043 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 101))
708 : END IF
709 : END IF
710 : ELSE
711 741515 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
712 453424 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
713 266986 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
714 266986 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
715 266986 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 102))
716 : ELSE
717 186438 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
718 186438 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
719 186438 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 103))
720 : END IF
721 : ELSE
722 288091 : IF (X1 <= 0.468750000000000000E-01_dp) THEN
723 131334 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
724 131334 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
725 131334 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 104))
726 : ELSE
727 156757 : TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
728 156757 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
729 156757 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 105))
730 : END IF
731 : END IF
732 : END IF
733 : ELSE
734 1622173 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
735 803294 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
736 340882 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
737 208144 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
738 208144 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
739 208144 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 106))
740 : ELSE
741 132738 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
742 132738 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
743 132738 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 107))
744 : END IF
745 : ELSE
746 462412 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
747 247017 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
748 247017 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
749 247017 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 108))
750 : ELSE
751 215395 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
752 215395 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
753 215395 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 109))
754 : END IF
755 : END IF
756 : ELSE
757 818879 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
758 330543 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
759 330543 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
760 330543 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 110))
761 : ELSE
762 488336 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
763 488336 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
764 488336 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 111))
765 : END IF
766 : END IF
767 : END IF
768 : ELSE
769 3627776 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
770 1373185 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
771 395615 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
772 186605 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
773 92207 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
774 92207 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
775 92207 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 112))
776 : ELSE
777 94398 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
778 94398 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
779 94398 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 113))
780 : END IF
781 : ELSE
782 209010 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
783 128264 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
784 128264 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
785 128264 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 114))
786 : ELSE
787 80746 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
788 80746 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
789 80746 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 115))
790 : END IF
791 : END IF
792 : ELSE
793 977570 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
794 427212 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
795 186003 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
796 186003 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
797 186003 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 116))
798 : ELSE
799 241209 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
800 241209 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
801 241209 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 117))
802 : END IF
803 : ELSE
804 550358 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
805 272144 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
806 272144 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
807 272144 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 118))
808 : ELSE
809 278214 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
810 278214 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
811 278214 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 119))
812 : END IF
813 : END IF
814 : END IF
815 : ELSE
816 2254591 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
817 1024075 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
818 538924 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
819 538924 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
820 538924 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 120))
821 : ELSE
822 485151 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
823 485151 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
824 485151 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 121))
825 : END IF
826 : ELSE
827 1230516 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
828 690786 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
829 344811 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
830 344811 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
831 344811 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 122))
832 : ELSE
833 345975 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
834 345975 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
835 345975 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 123))
836 : END IF
837 : ELSE
838 539730 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
839 539730 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
840 539730 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 124))
841 : END IF
842 : END IF
843 : END IF
844 : END IF
845 : ELSE
846 9611013 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
847 4365958 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
848 1577614 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
849 743094 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
850 181131 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
851 89988 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
852 89988 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
853 89988 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 125))
854 : ELSE
855 91143 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
856 91143 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
857 91143 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 126))
858 : END IF
859 : ELSE
860 561963 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
861 277534 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
862 277534 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
863 277534 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 127))
864 : ELSE
865 284429 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
866 284429 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
867 284429 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 128))
868 : END IF
869 : END IF
870 : ELSE
871 834520 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
872 178082 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
873 178082 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
874 178082 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 129))
875 : ELSE
876 656438 : IF (X1 <= 0.343750000000000000E+00_dp) THEN
877 320773 : TG1 = (2*X1 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
878 320773 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
879 320773 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 130))
880 : ELSE
881 335665 : TG1 = (2*X1 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
882 335665 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
883 335665 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 131))
884 : END IF
885 : END IF
886 : END IF
887 : ELSE
888 2788344 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
889 1365059 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
890 765501 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
891 374337 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
892 374337 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
893 374337 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 132))
894 : ELSE
895 391164 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
896 391164 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
897 391164 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 133))
898 : END IF
899 : ELSE
900 599558 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
901 306444 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
902 306444 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
903 306444 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 134))
904 : ELSE
905 293114 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
906 293114 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
907 293114 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 135))
908 : END IF
909 : END IF
910 : ELSE
911 1423285 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
912 806680 : IF (X1 <= 0.343750000000000000E+00_dp) THEN
913 353698 : TG1 = (2*X1 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
914 353698 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
915 353698 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 136))
916 : ELSE
917 452982 : TG1 = (2*X1 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
918 452982 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
919 452982 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 137))
920 : END IF
921 : ELSE
922 616605 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
923 616605 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
924 616605 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 138))
925 : END IF
926 : END IF
927 : END IF
928 : ELSE
929 5245055 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
930 1757305 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
931 855629 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
932 172976 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
933 172976 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
934 172976 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 139))
935 : ELSE
936 682653 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
937 682653 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
938 682653 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 140))
939 : END IF
940 : ELSE
941 901676 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
942 901676 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
943 901676 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 141))
944 : END IF
945 : ELSE
946 3487750 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
947 1700334 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
948 954734 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
949 954734 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
950 954734 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 142))
951 : ELSE
952 745600 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
953 745600 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
954 745600 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 143))
955 : END IF
956 : ELSE
957 1787416 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
958 941470 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
959 941470 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
960 941470 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 144))
961 : ELSE
962 845946 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
963 845946 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
964 845946 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 145))
965 : END IF
966 : END IF
967 : END IF
968 : END IF
969 : END IF
970 : END IF
971 : ELSE
972 72032690 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
973 55193004 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
974 42951907 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
975 27342474 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
976 21945254 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
977 21945254 : TG2 = (2*X2 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
978 21945254 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 146))
979 : ELSE
980 5397220 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
981 5397220 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
982 5397220 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 147))
983 : END IF
984 : ELSE
985 15609433 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
986 15609433 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
987 15609433 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 148))
988 : END IF
989 : ELSE
990 12241097 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
991 5951360 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
992 2042726 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
993 1025721 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
994 1025721 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
995 1025721 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 149))
996 : ELSE
997 1017005 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
998 1017005 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
999 1017005 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 150))
1000 : END IF
1001 : ELSE
1002 3908634 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1003 1908984 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1004 1908984 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1005 1908984 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 151))
1006 : ELSE
1007 1999650 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1008 1999650 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1009 1999650 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 152))
1010 : END IF
1011 : END IF
1012 : ELSE
1013 6289737 : IF (X1 <= 0.687500000000000000E+00_dp) THEN
1014 3146367 : TG1 = (2*X1 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
1015 3146367 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1016 3146367 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 153))
1017 : ELSE
1018 3143370 : TG1 = (2*X1 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
1019 3143370 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1020 3143370 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 154))
1021 : END IF
1022 : END IF
1023 : END IF
1024 : ELSE
1025 16839686 : TG1 = (2*X1 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1026 16839686 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1027 16839686 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 155))
1028 : END IF
1029 : END IF
1030 : ELSE
1031 37760304 : IF (T < lower) THEN
1032 5411488 : use_gamma = .TRUE.
1033 5411488 : RETURN
1034 : END IF
1035 32348816 : X2 = 11.0_dp/R
1036 32348816 : X1 = (T - lower)/(upper - lower)
1037 32348816 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
1038 13092879 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
1039 6054102 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1040 489976 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
1041 260511 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
1042 2618 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
1043 2618 : TG2 = (2*X2 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
1044 2618 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 156))
1045 : ELSE
1046 257893 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
1047 257893 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
1048 257893 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 157))
1049 : END IF
1050 : ELSE
1051 229465 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
1052 229465 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1053 229465 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 158))
1054 : END IF
1055 : ELSE
1056 5564126 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
1057 2520629 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
1058 1192389 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
1059 479258 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
1060 479258 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
1061 479258 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 159))
1062 : ELSE
1063 713131 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
1064 316690 : TG1 = (2*X1 - 0.625000000000000000E-01_dp)*0.160000000000000000E+02_dp
1065 316690 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1066 316690 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 160))
1067 : ELSE
1068 396441 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
1069 396441 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1070 396441 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 161))
1071 : END IF
1072 : END IF
1073 : ELSE
1074 1328240 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
1075 538227 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1076 343190 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
1077 153059 : IF (X2 <= 0.812500000000000000E+00_dp) THEN
1078 106452 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
1079 106452 : TG2 = (2*X2 - 0.156250000000000000E+01_dp)*0.160000000000000000E+02_dp
1080 106452 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 162))
1081 : ELSE
1082 46607 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
1083 46607 : TG2 = (2*X2 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
1084 46607 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 163))
1085 : END IF
1086 : ELSE
1087 190131 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
1088 190131 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1089 190131 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 164))
1090 : END IF
1091 : ELSE
1092 195037 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
1093 90629 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1094 36488 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
1095 14181 : IF (X2 <= 0.906250000000000000E+00_dp) THEN
1096 6092 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
1097 6092 : TG2 = (2*X2 - 0.178125000000000000E+01_dp)*0.320000000000000000E+02_dp
1098 6092 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 165))
1099 : ELSE
1100 8089 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
1101 8089 : TG2 = (2*X2 - 0.184375000000000000E+01_dp)*0.320000000000000000E+02_dp
1102 8089 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 166))
1103 : END IF
1104 : ELSE
1105 22307 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
1106 22307 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
1107 22307 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 167))
1108 : END IF
1109 : ELSE
1110 54141 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
1111 26739 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
1112 16921 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
1113 8674 : IF (X2 <= 0.953125000000000000E+00_dp) THEN
1114 3380 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
1115 3380 : TG2 = (2*X2 - 0.189062500000000000E+01_dp)*0.640000000000000000E+02_dp
1116 3380 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 168))
1117 : ELSE
1118 5294 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
1119 5294 : TG2 = (2*X2 - 0.192187500000000000E+01_dp)*0.640000000000000000E+02_dp
1120 5294 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 169))
1121 : END IF
1122 : ELSE
1123 8247 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
1124 8247 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
1125 8247 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 170))
1126 : END IF
1127 : ELSE
1128 9818 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
1129 4028 : IF (X2 <= 0.984375000000000000E+00_dp) THEN
1130 1365 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
1131 1365 : TG2 = (2*X2 - 0.195312500000000000E+01_dp)*0.640000000000000000E+02_dp
1132 1365 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 171))
1133 : ELSE
1134 2663 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
1135 2663 : TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
1136 2663 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 172))
1137 : END IF
1138 : ELSE
1139 5790 : IF (X2 <= 0.984375000000000000E+00_dp) THEN
1140 1371 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
1141 1371 : TG2 = (2*X2 - 0.195312500000000000E+01_dp)*0.640000000000000000E+02_dp
1142 1371 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 173))
1143 : ELSE
1144 4419 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
1145 4419 : TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
1146 4419 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 174))
1147 : END IF
1148 : END IF
1149 : END IF
1150 : ELSE
1151 27402 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
1152 12401 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
1153 12401 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
1154 12401 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 175))
1155 : ELSE
1156 15001 : IF (X1 <= 0.234375000000000000E-01_dp) THEN
1157 7436 : TG1 = (2*X1 - 0.390625000000000000E-01_dp)*0.128000000000000000E+03_dp
1158 7436 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1159 7436 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 176))
1160 : ELSE
1161 7565 : TG1 = (2*X1 - 0.546875000000000000E-01_dp)*0.128000000000000000E+03_dp
1162 7565 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1163 7565 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 177))
1164 : END IF
1165 : END IF
1166 : END IF
1167 : END IF
1168 : ELSE
1169 104408 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1170 43814 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
1171 43814 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
1172 43814 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 178))
1173 : ELSE
1174 60594 : IF (X1 <= 0.468750000000000000E-01_dp) THEN
1175 28533 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
1176 20615 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
1177 20615 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
1178 20615 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 179))
1179 : ELSE
1180 7918 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
1181 7918 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1182 7918 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 180))
1183 : END IF
1184 : ELSE
1185 32061 : TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
1186 32061 : TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
1187 32061 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 181))
1188 : END IF
1189 : END IF
1190 : END IF
1191 : END IF
1192 : ELSE
1193 790013 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1194 451428 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
1195 451428 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1196 451428 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 182))
1197 : ELSE
1198 338585 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
1199 141792 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1200 82666 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
1201 82666 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
1202 82666 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 183))
1203 : ELSE
1204 59126 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
1205 59126 : TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
1206 59126 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 184))
1207 : END IF
1208 : ELSE
1209 196793 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
1210 196793 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1211 196793 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 185))
1212 : END IF
1213 : END IF
1214 : END IF
1215 : END IF
1216 : ELSE
1217 3043497 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
1218 1137416 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
1219 1137416 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
1220 1137416 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 186))
1221 : ELSE
1222 1906081 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
1223 907967 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1224 494185 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
1225 494185 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1226 494185 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 187))
1227 : ELSE
1228 413782 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
1229 413782 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1230 413782 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 188))
1231 : END IF
1232 : ELSE
1233 998114 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
1234 998114 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1235 998114 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 189))
1236 : END IF
1237 : END IF
1238 : END IF
1239 : END IF
1240 : ELSE
1241 7038777 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
1242 2891288 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
1243 1421999 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
1244 713643 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
1245 713643 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1246 713643 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 190))
1247 : ELSE
1248 708356 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
1249 708356 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1250 708356 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 191))
1251 : END IF
1252 : ELSE
1253 1469289 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1254 78812 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
1255 78812 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1256 78812 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 192))
1257 : ELSE
1258 1390477 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
1259 1390477 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1260 1390477 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 193))
1261 : END IF
1262 : END IF
1263 : ELSE
1264 4147489 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
1265 1526455 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1266 65553 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
1267 65553 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1268 65553 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 194))
1269 : ELSE
1270 1460902 : IF (X1 <= 0.406250000000000000E+00_dp) THEN
1271 645992 : TG1 = (2*X1 - 0.781250000000000000E+00_dp)*0.320000000000000000E+02_dp
1272 645992 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1273 645992 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 195))
1274 : ELSE
1275 814910 : TG1 = (2*X1 - 0.843750000000000000E+00_dp)*0.320000000000000000E+02_dp
1276 814910 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1277 814910 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 196))
1278 : END IF
1279 : END IF
1280 : ELSE
1281 2621034 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1282 358389 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
1283 358389 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1284 358389 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 197))
1285 : ELSE
1286 2262645 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
1287 2262645 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1288 2262645 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 198))
1289 : END IF
1290 : END IF
1291 : END IF
1292 : END IF
1293 : ELSE
1294 19255937 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
1295 11509950 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
1296 5409247 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1297 221254 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1298 97115 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1299 97115 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1300 97115 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 199))
1301 : ELSE
1302 124139 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1303 124139 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1304 124139 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 200))
1305 : END IF
1306 : ELSE
1307 5187993 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1308 2393833 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1309 2393833 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1310 2393833 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 201))
1311 : ELSE
1312 2794160 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1313 2794160 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1314 2794160 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 202))
1315 : END IF
1316 : END IF
1317 : ELSE
1318 6100703 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1319 558015 : IF (X1 <= 0.687500000000000000E+00_dp) THEN
1320 202445 : TG1 = (2*X1 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
1321 202445 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1322 202445 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 203))
1323 : ELSE
1324 355570 : TG1 = (2*X1 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
1325 355570 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1326 355570 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 204))
1327 : END IF
1328 : ELSE
1329 5542688 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1330 5542688 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1331 5542688 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 205))
1332 : END IF
1333 : END IF
1334 : ELSE
1335 7745987 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
1336 4854866 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1337 4854866 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1338 4854866 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 206))
1339 : ELSE
1340 2891121 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1341 2891121 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1342 2891121 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 207))
1343 : END IF
1344 : END IF
1345 : END IF
1346 : END IF
1347 : END SUBROUTINE t_c_g0_n
1348 :
1349 : ! **************************************************************************************************
1350 : !> \brief ...
1351 : !> \param Nder the number of derivatives that will actually be used
1352 : !> \param iunit contains the data file to initialize the table
1353 : !> \param mepos ...
1354 : !> \param group ...
1355 : ! **************************************************************************************************
1356 770 : SUBROUTINE init(Nder, iunit, mepos, group)
1357 : INTEGER, INTENT(IN) :: Nder, iunit, mepos
1358 :
1359 : CLASS(mp_comm_type), INTENT(IN) :: group
1360 :
1361 : INTEGER :: I
1362 770 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1363 :
1364 770 : patches = 207
1365 770 : IF (Nder > nderiv_max) CPABORT("T_C_G0 init failed")
1366 770 : nderiv_init = Nder
1367 770 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1368 : ! round up to a multiple of 32 to give some generous alignment for each C0
1369 3080 : ALLOCATE (C0(32*((31 + (Nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1370 : ! init mpi'ed buffers to silence warnings under valgrind
1371 117590432 : C0 = 1.0E99_dp
1372 770 : IF (mepos == 0) THEN
1373 385 : ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1374 80080 : DO I = 1, patches
1375 79695 : READ (iunit, *) chunk
1376 57721300 : C0(1:(Nder + 1)*(degree + 1)*(degree + 2)/2, I) = chunk(1:(Nder + 1)*(degree + 1)*(degree + 2)/2)
1377 : END DO
1378 385 : DEALLOCATE (chunk)
1379 : END IF
1380 770 : CALL group%bcast(C0, 0)
1381 :
1382 770 : END SUBROUTINE init
1383 :
1384 : ! **************************************************************************************************
1385 : !> \brief ...
1386 : ! **************************************************************************************************
1387 376 : SUBROUTINE free_C0()
1388 376 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1389 376 : nderiv_init = -1
1390 376 : END SUBROUTINE free_C0
1391 :
1392 : ! **************************************************************************************************
1393 : !> \brief ...
1394 : !> \param RES ...
1395 : !> \param NDERIV ...
1396 : !> \param TG1 ...
1397 : !> \param TG2 ...
1398 : !> \param C0 ...
1399 : ! **************************************************************************************************
1400 242429212 : SUBROUTINE PD2VAL(RES, NDERIV, TG1, TG2, C0)
1401 : REAL(KIND=dp), INTENT(OUT) :: res(*)
1402 : INTEGER, INTENT(IN) :: NDERIV
1403 : REAL(KIND=dp), INTENT(IN) :: TG1, TG2, C0(105, *)
1404 :
1405 : REAL(KIND=dp), PARAMETER :: SQRT2 = 1.4142135623730950488016887242096980785696718753_dp
1406 :
1407 : INTEGER :: K
1408 : REAL(KIND=dp) :: T1(0:13), T2(0:13)
1409 :
1410 242429212 : T1(0) = 1.0_dp
1411 242429212 : T2(0) = 1.0_dp
1412 242429212 : T1(1) = SQRT2*TG1
1413 242429212 : T2(1) = SQRT2*TG2
1414 242429212 : T1(2) = 2*TG1*T1(1) - SQRT2
1415 242429212 : T2(2) = 2*TG2*T2(1) - SQRT2
1416 242429212 : T1(3) = 2*TG1*T1(2) - T1(1)
1417 242429212 : T2(3) = 2*TG2*T2(2) - T2(1)
1418 242429212 : T1(4) = 2*TG1*T1(3) - T1(2)
1419 242429212 : T2(4) = 2*TG2*T2(3) - T2(2)
1420 242429212 : T1(5) = 2*TG1*T1(4) - T1(3)
1421 242429212 : T2(5) = 2*TG2*T2(4) - T2(3)
1422 242429212 : T1(6) = 2*TG1*T1(5) - T1(4)
1423 242429212 : T2(6) = 2*TG2*T2(5) - T2(4)
1424 242429212 : T1(7) = 2*TG1*T1(6) - T1(5)
1425 242429212 : T2(7) = 2*TG2*T2(6) - T2(5)
1426 242429212 : T1(8) = 2*TG1*T1(7) - T1(6)
1427 242429212 : T2(8) = 2*TG2*T2(7) - T2(6)
1428 242429212 : T1(9) = 2*TG1*T1(8) - T1(7)
1429 242429212 : T2(9) = 2*TG2*T2(8) - T2(7)
1430 242429212 : T1(10) = 2*TG1*T1(9) - T1(8)
1431 242429212 : T2(10) = 2*TG2*T2(9) - T2(8)
1432 242429212 : T1(11) = 2*TG1*T1(10) - T1(9)
1433 242429212 : T2(11) = 2*TG2*T2(10) - T2(9)
1434 242429212 : T1(12) = 2*TG1*T1(11) - T1(10)
1435 242429212 : T2(12) = 2*TG2*T2(11) - T2(10)
1436 242429212 : T1(13) = 2*TG1*T1(12) - T1(11)
1437 242429212 : T2(13) = 2*TG2*T2(12) - T2(11)
1438 842720682 : DO K = 1, NDERIV + 1
1439 600291470 : RES(K) = 0.0_dp
1440 9004372050 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:13), C0(1:14, K))*T2(0)
1441 8404080580 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:12), C0(15:27, K))*T2(1)
1442 7803789110 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:11), C0(28:39, K))*T2(2)
1443 7203497640 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:10), C0(40:50, K))*T2(3)
1444 6603206170 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:9), C0(51:60, K))*T2(4)
1445 6002914700 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:8), C0(61:69, K))*T2(5)
1446 5402623230 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:7), C0(70:77, K))*T2(6)
1447 4802331760 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:6), C0(78:84, K))*T2(7)
1448 4202040290 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:5), C0(85:90, K))*T2(8)
1449 3601748820 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:4), C0(91:95, K))*T2(9)
1450 3001457350 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:3), C0(96:99, K))*T2(10)
1451 2401165880 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:2), C0(100:102, K))*T2(11)
1452 1800874410 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:1), C0(103:104, K))*T2(12)
1453 1443012152 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:0), C0(105:105, K))*T2(13)
1454 : END DO
1455 242429212 : END SUBROUTINE PD2VAL
1456 :
1457 : ! **************************************************************************************************
1458 : !> \brief Returns the value of nderiv_init so that one can check if opening the potential file is
1459 : !> worhtwhile
1460 : !> \return ...
1461 : !> \author A. Bussy, 05.2019
1462 : ! **************************************************************************************************
1463 7171405 : FUNCTION get_lmax_init() RESULT(lmax_init)
1464 :
1465 : INTEGER :: lmax_init
1466 :
1467 7171405 : lmax_init = nderiv_init
1468 :
1469 7171405 : END FUNCTION get_lmax_init
1470 :
1471 : END MODULE t_c_g0
|
