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