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 = 1.000000E-09_dp
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 252890068 : 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 252890068 : use_gamma = .FALSE.
96 252890068 : upper = R**2 + 11.0_dp*R + 50.0_dp
97 252890068 : lower = R**2 - 11.0_dp*R + 0.0_dp
98 252890068 : IF (T > upper) THEN
99 7481834 : RES(1:NDERIV + 1) = 0.0_dp
100 7398836 : RETURN
101 : END IF
102 250915274 : IF (R <= 11.0_dp) THEN
103 213191774 : X2 = R/11.0_dp
104 213191774 : upper = R**2 + 11.0_dp*R + 50.0_dp
105 213191774 : lower = 0.0_dp
106 213191774 : X1 = (T - lower)/(upper - lower)
107 213191774 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
108 141232862 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
109 118522566 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
110 95646456 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
111 20277361 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
112 9422462 : IF (X2 <= 0.625000000000000000E-01_dp) THEN
113 1186069 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
114 611570 : IF (X2 <= 0.312500000000000000E-01_dp) THEN
115 41646 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
116 16463 : 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 16463 : 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 8398 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
133 8398 : TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
134 8398 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 4))
135 : END IF
136 : END IF
137 : ELSE
138 25183 : 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 25183 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
144 25183 : TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
145 25183 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 6))
146 : END IF
147 : END IF
148 : ELSE
149 569924 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
150 301805 : IF (X2 <= 0.468750000000000000E-01_dp) THEN
151 132231 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
152 67918 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
153 67918 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
154 67918 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 7))
155 : ELSE
156 64313 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
157 64313 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
158 64313 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 8))
159 : END IF
160 : ELSE
161 169574 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
162 129882 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
163 129882 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
164 129882 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 9))
165 : ELSE
166 39692 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
167 39692 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
168 39692 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 10))
169 : END IF
170 : END IF
171 : ELSE
172 268119 : IF (X2 <= 0.468750000000000000E-01_dp) THEN
173 205090 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
174 205090 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
175 205090 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 11))
176 : ELSE
177 63029 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
178 63029 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
179 63029 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 12))
180 : END IF
181 : END IF
182 : END IF
183 : ELSE
184 574499 : IF (X2 <= 0.312500000000000000E-01_dp) THEN
185 50963 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
186 16150 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
187 16150 : TG2 = (2*X2 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
188 16150 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 13))
189 : ELSE
190 34813 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
191 34813 : TG2 = (2*X2 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
192 34813 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 14))
193 : END IF
194 : ELSE
195 523536 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
196 281909 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
197 281909 : TG2 = (2*X2 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
198 281909 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 15))
199 : ELSE
200 241627 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
201 241627 : TG2 = (2*X2 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
202 241627 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 16))
203 : END IF
204 : END IF
205 : END IF
206 : ELSE
207 8236393 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
208 3982768 : IF (X2 <= 0.937500000000000000E-01_dp) THEN
209 928904 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
210 600100 : IF (X2 <= 0.781250000000000000E-01_dp) THEN
211 272561 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
212 239466 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
213 239466 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
214 239466 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 17))
215 : ELSE
216 33095 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
217 33095 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
218 33095 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 18))
219 : END IF
220 : ELSE
221 327539 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
222 237198 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
223 237198 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
224 237198 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 19))
225 : ELSE
226 90341 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
227 90341 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
228 90341 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 20))
229 : END IF
230 : END IF
231 : ELSE
232 328804 : IF (X2 <= 0.781250000000000000E-01_dp) THEN
233 65612 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
234 65612 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
235 65612 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 21))
236 : ELSE
237 263192 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
238 263192 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
239 263192 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 22))
240 : END IF
241 : END IF
242 : ELSE
243 3053864 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
244 1566545 : IF (X2 <= 0.109375000000000000E+00_dp) THEN
245 724739 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
246 517019 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
247 517019 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
248 517019 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 23))
249 : ELSE
250 207720 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
251 207720 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
252 207720 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 24))
253 : END IF
254 : ELSE
255 841806 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
256 609086 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
257 609086 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
258 609086 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 25))
259 : ELSE
260 232720 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
261 232720 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
262 232720 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 26))
263 : END IF
264 : END IF
265 : ELSE
266 1487319 : IF (X2 <= 0.109375000000000000E+00_dp) THEN
267 616302 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
268 616302 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
269 616302 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 27))
270 : ELSE
271 871017 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
272 871017 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
273 871017 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 28))
274 : END IF
275 : END IF
276 : END IF
277 : ELSE
278 4253625 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
279 2048402 : IF (X2 <= 0.937500000000000000E-01_dp) THEN
280 342525 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
281 342525 : TG2 = (2*X2 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
282 342525 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 29))
283 : ELSE
284 1705877 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
285 1705877 : TG2 = (2*X2 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
286 1705877 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 30))
287 : END IF
288 : ELSE
289 2205223 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
290 2205223 : TG2 = (2*X2 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
291 2205223 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 31))
292 : END IF
293 : END IF
294 : END IF
295 : ELSE
296 10854899 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
297 5777634 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
298 5777634 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
299 5777634 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 32))
300 : ELSE
301 5077265 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
302 5077265 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
303 5077265 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 33))
304 : END IF
305 : END IF
306 : ELSE
307 75369095 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
308 33166329 : IF (X2 <= 0.187500000000000000E+00_dp) THEN
309 20816718 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
310 10162619 : IF (X2 <= 0.156250000000000000E+00_dp) THEN
311 4692546 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
312 2183875 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
313 1412814 : IF (X2 <= 0.140625000000000000E+00_dp) THEN
314 620750 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
315 620750 : TG2 = (2*X2 - 0.265625000000000000E+00_dp)*0.640000000000000000E+02_dp
316 620750 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 34))
317 : ELSE
318 792064 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
319 792064 : TG2 = (2*X2 - 0.296875000000000000E+00_dp)*0.640000000000000000E+02_dp
320 792064 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 35))
321 : END IF
322 : ELSE
323 771061 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
324 771061 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
325 771061 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 36))
326 : END IF
327 : ELSE
328 2508671 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
329 1126812 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
330 1126812 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
331 1126812 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 37))
332 : ELSE
333 1381859 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
334 1381859 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
335 1381859 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 38))
336 : END IF
337 : END IF
338 : ELSE
339 5470073 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
340 3238896 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
341 2171868 : IF (X2 <= 0.171875000000000000E+00_dp) THEN
342 898649 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
343 898649 : TG2 = (2*X2 - 0.328125000000000000E+00_dp)*0.640000000000000000E+02_dp
344 898649 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 39))
345 : ELSE
346 1273219 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
347 1273219 : TG2 = (2*X2 - 0.359375000000000000E+00_dp)*0.640000000000000000E+02_dp
348 1273219 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 40))
349 : END IF
350 : ELSE
351 1067028 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
352 1067028 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
353 1067028 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 41))
354 : END IF
355 : ELSE
356 2231177 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
357 2231177 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
358 2231177 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 42))
359 : END IF
360 : END IF
361 : ELSE
362 10654099 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
363 5117133 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
364 5117133 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
365 5117133 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 43))
366 : ELSE
367 5536966 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
368 5536966 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
369 5536966 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 44))
370 : END IF
371 : END IF
372 : ELSE
373 12349611 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
374 8133446 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
375 6030906 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
376 3081770 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
377 2119884 : IF (X2 <= 0.203125000000000000E+00_dp) THEN
378 1096961 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
379 1096961 : TG2 = (2*X2 - 0.390625000000000000E+00_dp)*0.640000000000000000E+02_dp
380 1096961 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 45))
381 : ELSE
382 1022923 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
383 1022923 : TG2 = (2*X2 - 0.421875000000000000E+00_dp)*0.640000000000000000E+02_dp
384 1022923 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 46))
385 : END IF
386 : ELSE
387 961886 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
388 961886 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
389 961886 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 47))
390 : END IF
391 : ELSE
392 2949136 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
393 2242122 : IF (X2 <= 0.234375000000000000E+00_dp) THEN
394 1056496 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
395 1056496 : TG2 = (2*X2 - 0.453125000000000000E+00_dp)*0.640000000000000000E+02_dp
396 1056496 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 48))
397 : ELSE
398 1185626 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
399 1185626 : TG2 = (2*X2 - 0.484375000000000000E+00_dp)*0.640000000000000000E+02_dp
400 1185626 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 49))
401 : END IF
402 : ELSE
403 707014 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
404 707014 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
405 707014 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 50))
406 : END IF
407 : END IF
408 : ELSE
409 2102540 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
410 1270466 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
411 1270466 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
412 1270466 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 51))
413 : ELSE
414 832074 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
415 832074 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
416 832074 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 52))
417 : END IF
418 : END IF
419 : ELSE
420 4216165 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
421 2032158 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
422 2032158 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
423 2032158 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 53))
424 : ELSE
425 2184007 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
426 2184007 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
427 2184007 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 54))
428 : END IF
429 : END IF
430 : END IF
431 : ELSE
432 42202766 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
433 18925153 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
434 8647601 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
435 8647601 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
436 8647601 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 55))
437 : ELSE
438 10277552 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
439 10277552 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
440 10277552 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 56))
441 : END IF
442 : ELSE
443 23277613 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
444 23277613 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
445 23277613 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 57))
446 : END IF
447 : END IF
448 : END IF
449 : ELSE
450 22876110 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
451 17464909 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
452 15741932 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
453 13738523 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
454 8953505 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
455 4963138 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
456 3958733 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
457 2034060 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
458 2034060 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
459 2034060 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 58))
460 : ELSE
461 1924673 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
462 1924673 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
463 1924673 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 59))
464 : END IF
465 : ELSE
466 1004405 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
467 549584 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
468 549584 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
469 549584 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 60))
470 : ELSE
471 454821 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
472 454821 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
473 454821 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 61))
474 : END IF
475 : END IF
476 : ELSE
477 3990367 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
478 3393113 : IF (X2 <= 0.343750000000000000E+00_dp) THEN
479 1697197 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
480 1697197 : TG2 = (2*X2 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
481 1697197 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 62))
482 : ELSE
483 1695916 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
484 1695916 : TG2 = (2*X2 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
485 1695916 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 63))
486 : END IF
487 : ELSE
488 597254 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
489 597254 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
490 597254 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 64))
491 : END IF
492 : END IF
493 : ELSE
494 4785018 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
495 4286458 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
496 2664224 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
497 2409273 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
498 2409273 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
499 2409273 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 65))
500 : ELSE
501 254951 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
502 254951 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
503 254951 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 66))
504 : END IF
505 : ELSE
506 1622234 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
507 1502526 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
508 1502526 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
509 1502526 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 67))
510 : ELSE
511 119708 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
512 119708 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
513 119708 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 68))
514 : END IF
515 : END IF
516 : ELSE
517 498560 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
518 318020 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
519 318020 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
520 318020 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 69))
521 : ELSE
522 180540 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
523 180540 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
524 180540 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 70))
525 : END IF
526 : END IF
527 : END IF
528 : ELSE
529 2003409 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
530 1567708 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
531 1056190 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
532 684350 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
533 684350 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
534 684350 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 71))
535 : ELSE
536 371840 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
537 371840 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
538 371840 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 72))
539 : END IF
540 : ELSE
541 511518 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
542 292554 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
543 292554 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
544 292554 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 73))
545 : ELSE
546 218964 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
547 218964 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
548 218964 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 74))
549 : END IF
550 : END IF
551 : ELSE
552 435701 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
553 238174 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
554 157669 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
555 157669 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
556 157669 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 75))
557 : ELSE
558 80505 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
559 80505 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
560 80505 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 76))
561 : END IF
562 : ELSE
563 197527 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
564 128367 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
565 128367 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
566 128367 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 77))
567 : ELSE
568 69160 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
569 69160 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
570 69160 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 78))
571 : END IF
572 : END IF
573 : END IF
574 : END IF
575 : ELSE
576 1722977 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
577 1473530 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
578 612763 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
579 365111 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
580 365111 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
581 365111 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 79))
582 : ELSE
583 247652 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
584 172147 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
585 172147 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
586 172147 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 80))
587 : ELSE
588 75505 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
589 75505 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
590 75505 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 81))
591 : END IF
592 : END IF
593 : ELSE
594 860767 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
595 744300 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
596 744300 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
597 744300 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 82))
598 : ELSE
599 116467 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
600 116467 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
601 116467 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 83))
602 : END IF
603 : END IF
604 : ELSE
605 249447 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
606 183691 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
607 112865 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
608 112865 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
609 112865 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 84))
610 : ELSE
611 70826 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
612 44122 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
613 44122 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
614 44122 : 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 65756 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
623 40110 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
624 40110 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
625 40110 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 87))
626 : ELSE
627 25646 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
628 25646 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
629 25646 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 88))
630 : END IF
631 : END IF
632 : END IF
633 : END IF
634 : ELSE
635 5411201 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
636 1650547 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
637 1517432 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
638 706444 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
639 706444 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
640 706444 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 89))
641 : ELSE
642 810988 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
643 810988 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
644 810988 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 90))
645 : END IF
646 : ELSE
647 133115 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
648 65444 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
649 31377 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
650 31377 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
651 31377 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 91))
652 : ELSE
653 34067 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
654 34067 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
655 34067 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 92))
656 : END IF
657 : ELSE
658 67671 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
659 67671 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
660 67671 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 93))
661 : END IF
662 : END IF
663 : ELSE
664 3760654 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
665 1598982 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
666 1598982 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
667 1598982 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 94))
668 : ELSE
669 2161672 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
670 2161672 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
671 2161672 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 95))
672 : END IF
673 : END IF
674 : END IF
675 : END IF
676 : ELSE
677 22710296 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
678 13151620 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
679 9542813 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
680 7925466 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
681 7176400 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
682 6724598 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
683 6281334 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
684 4756276 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
685 4756276 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
686 4756276 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 96))
687 : ELSE
688 1525058 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
689 1525058 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
690 1525058 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 97))
691 : END IF
692 : ELSE
693 443264 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
694 339372 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
695 339372 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
696 339372 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 98))
697 : ELSE
698 103892 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
699 103892 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
700 103892 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 99))
701 : END IF
702 : END IF
703 : ELSE
704 451802 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
705 245338 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
706 245338 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
707 245338 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 100))
708 : ELSE
709 206464 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
710 206464 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
711 206464 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 101))
712 : END IF
713 : END IF
714 : ELSE
715 749066 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
716 456028 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
717 267931 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
718 267931 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
719 267931 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 102))
720 : ELSE
721 188097 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
722 188097 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
723 188097 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 103))
724 : END IF
725 : ELSE
726 293038 : IF (X1 <= 0.468750000000000000E-01_dp) THEN
727 136699 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
728 136699 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
729 136699 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 104))
730 : ELSE
731 156339 : TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
732 156339 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
733 156339 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 105))
734 : END IF
735 : END IF
736 : END IF
737 : ELSE
738 1617347 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
739 802453 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
740 340955 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
741 208350 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
742 208350 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
743 208350 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 106))
744 : ELSE
745 132605 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
746 132605 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
747 132605 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 107))
748 : END IF
749 : ELSE
750 461498 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
751 247175 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
752 247175 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
753 247175 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 108))
754 : ELSE
755 214323 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
756 214323 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
757 214323 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 109))
758 : END IF
759 : END IF
760 : ELSE
761 814894 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
762 328816 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
763 328816 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
764 328816 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 110))
765 : ELSE
766 486078 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
767 486078 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
768 486078 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 111))
769 : END IF
770 : END IF
771 : END IF
772 : ELSE
773 3608807 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
774 1367280 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
775 395195 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
776 186187 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
777 91867 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
778 91867 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
779 91867 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 112))
780 : ELSE
781 94320 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
782 94320 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
783 94320 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 113))
784 : END IF
785 : ELSE
786 209008 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
787 128262 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
788 128262 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
789 128262 : 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 972085 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
798 424718 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
799 184713 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
800 184713 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
801 184713 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 116))
802 : ELSE
803 240005 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
804 240005 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
805 240005 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 117))
806 : END IF
807 : ELSE
808 547367 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
809 270743 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
810 270743 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
811 270743 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 118))
812 : ELSE
813 276624 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
814 276624 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
815 276624 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 119))
816 : END IF
817 : END IF
818 : END IF
819 : ELSE
820 2241527 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
821 1017271 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
822 535184 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
823 535184 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
824 535184 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 120))
825 : ELSE
826 482087 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
827 482087 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
828 482087 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 121))
829 : END IF
830 : ELSE
831 1224256 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
832 687458 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
833 343211 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
834 343211 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
835 343211 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 122))
836 : ELSE
837 344247 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
838 344247 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
839 344247 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 123))
840 : END IF
841 : ELSE
842 536798 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
843 536798 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
844 536798 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 124))
845 : END IF
846 : END IF
847 : END IF
848 : END IF
849 : ELSE
850 9558676 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
851 4342232 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
852 1570494 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
853 739744 : 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 558613 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
865 275956 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
866 275956 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
867 275956 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 127))
868 : ELSE
869 282657 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
870 282657 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
871 282657 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 128))
872 : END IF
873 : END IF
874 : ELSE
875 830750 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
876 178102 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
877 178102 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
878 178102 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 129))
879 : ELSE
880 652648 : IF (X1 <= 0.343750000000000000E+00_dp) THEN
881 318915 : TG1 = (2*X1 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
882 318915 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
883 318915 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 130))
884 : ELSE
885 333733 : TG1 = (2*X1 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
886 333733 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
887 333733 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 131))
888 : END IF
889 : END IF
890 : END IF
891 : ELSE
892 2771738 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
893 1357285 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
894 760725 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
895 371857 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
896 371857 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
897 371857 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 132))
898 : ELSE
899 388868 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
900 388868 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
901 388868 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 133))
902 : END IF
903 : ELSE
904 596560 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
905 305026 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
906 305026 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
907 305026 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 134))
908 : ELSE
909 291534 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
910 291534 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
911 291534 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 135))
912 : END IF
913 : END IF
914 : ELSE
915 1414453 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
916 801832 : IF (X1 <= 0.343750000000000000E+00_dp) THEN
917 351738 : TG1 = (2*X1 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
918 351738 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
919 351738 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 136))
920 : ELSE
921 450094 : TG1 = (2*X1 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
922 450094 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
923 450094 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 137))
924 : END IF
925 : ELSE
926 612621 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
927 612621 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
928 612621 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 138))
929 : END IF
930 : END IF
931 : END IF
932 : ELSE
933 5216444 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
934 1748824 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
935 851829 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
936 172992 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
937 172992 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
938 172992 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 139))
939 : ELSE
940 678837 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
941 678837 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
942 678837 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 140))
943 : END IF
944 : ELSE
945 896995 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
946 896995 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
947 896995 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 141))
948 : END IF
949 : ELSE
950 3467620 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
951 1690552 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
952 948998 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
953 948998 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
954 948998 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 142))
955 : ELSE
956 741554 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
957 741554 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
958 741554 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 143))
959 : END IF
960 : ELSE
961 1777068 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
962 935870 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
963 935870 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
964 935870 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 144))
965 : ELSE
966 841198 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
967 841198 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
968 841198 : 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 71958912 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
977 55172266 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
978 42995485 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
979 27446206 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
980 22036129 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
981 22036129 : TG2 = (2*X2 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
982 22036129 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 146))
983 : ELSE
984 5410077 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
985 5410077 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
986 5410077 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 147))
987 : END IF
988 : ELSE
989 15549279 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
990 15549279 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
991 15549279 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 148))
992 : END IF
993 : ELSE
994 12176781 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
995 5920135 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
996 2033925 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
997 1021423 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
998 1021423 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
999 1021423 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 149))
1000 : ELSE
1001 1012502 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1002 1012502 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
1003 1012502 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 150))
1004 : END IF
1005 : ELSE
1006 3886210 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1007 1898622 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1008 1898622 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1009 1898622 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 151))
1010 : ELSE
1011 1987588 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1012 1987588 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1013 1987588 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 152))
1014 : END IF
1015 : END IF
1016 : ELSE
1017 6256646 : IF (X1 <= 0.687500000000000000E+00_dp) THEN
1018 3130090 : TG1 = (2*X1 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
1019 3130090 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1020 3130090 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 153))
1021 : ELSE
1022 3126556 : TG1 = (2*X1 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
1023 3126556 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1024 3126556 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 154))
1025 : END IF
1026 : END IF
1027 : END IF
1028 : ELSE
1029 16786646 : TG1 = (2*X1 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1030 16786646 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1031 16786646 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 155))
1032 : END IF
1033 : END IF
1034 : ELSE
1035 37723500 : IF (T < lower) THEN
1036 5424042 : use_gamma = .TRUE.
1037 5424042 : RETURN
1038 : END IF
1039 32299458 : X2 = 11.0_dp/R
1040 32299458 : X1 = (T - lower)/(upper - lower)
1041 32299458 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
1042 13070693 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
1043 6046513 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1044 491308 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
1045 261807 : 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 259189 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
1051 259189 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
1052 259189 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 157))
1053 : END IF
1054 : ELSE
1055 229501 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
1056 229501 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1057 229501 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 158))
1058 : END IF
1059 : ELSE
1060 5555205 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
1061 2517850 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
1062 1192699 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
1063 479192 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
1064 479192 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
1065 479192 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 159))
1066 : ELSE
1067 713507 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
1068 316918 : TG1 = (2*X1 - 0.625000000000000000E-01_dp)*0.160000000000000000E+02_dp
1069 316918 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1070 316918 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 160))
1071 : ELSE
1072 396589 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
1073 396589 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1074 396589 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 161))
1075 : END IF
1076 : END IF
1077 : ELSE
1078 1325151 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
1079 537929 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1080 342973 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
1081 152927 : IF (X2 <= 0.812500000000000000E+00_dp) THEN
1082 106448 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
1083 106448 : TG2 = (2*X2 - 0.156250000000000000E+01_dp)*0.160000000000000000E+02_dp
1084 106448 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 162))
1085 : ELSE
1086 46479 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
1087 46479 : TG2 = (2*X2 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
1088 46479 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 163))
1089 : END IF
1090 : ELSE
1091 190046 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
1092 190046 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1093 190046 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 164))
1094 : END IF
1095 : ELSE
1096 194956 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
1097 90905 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1098 36516 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
1099 14320 : IF (X2 <= 0.906250000000000000E+00_dp) THEN
1100 6188 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
1101 6188 : TG2 = (2*X2 - 0.178125000000000000E+01_dp)*0.320000000000000000E+02_dp
1102 6188 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 165))
1103 : ELSE
1104 8132 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
1105 8132 : TG2 = (2*X2 - 0.184375000000000000E+01_dp)*0.320000000000000000E+02_dp
1106 8132 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 166))
1107 : END IF
1108 : ELSE
1109 22196 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
1110 22196 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
1111 22196 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 167))
1112 : END IF
1113 : ELSE
1114 54389 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
1115 26581 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
1116 16725 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
1117 8594 : IF (X2 <= 0.953125000000000000E+00_dp) THEN
1118 3300 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
1119 3300 : TG2 = (2*X2 - 0.189062500000000000E+01_dp)*0.640000000000000000E+02_dp
1120 3300 : 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 8131 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
1128 8131 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
1129 8131 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 170))
1130 : END IF
1131 : ELSE
1132 9856 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
1133 4060 : 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 2695 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
1139 2695 : TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
1140 2695 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 172))
1141 : END IF
1142 : ELSE
1143 5796 : IF (X2 <= 0.984375000000000000E+00_dp) THEN
1144 1363 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
1145 1363 : TG2 = (2*X2 - 0.195312500000000000E+01_dp)*0.640000000000000000E+02_dp
1146 1363 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 173))
1147 : ELSE
1148 4433 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
1149 4433 : TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
1150 4433 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 174))
1151 : END IF
1152 : END IF
1153 : END IF
1154 : ELSE
1155 27808 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
1156 12795 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
1157 12795 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
1158 12795 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 175))
1159 : ELSE
1160 15013 : IF (X1 <= 0.234375000000000000E-01_dp) THEN
1161 7468 : TG1 = (2*X1 - 0.390625000000000000E-01_dp)*0.128000000000000000E+03_dp
1162 7468 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1163 7468 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 176))
1164 : ELSE
1165 7545 : TG1 = (2*X1 - 0.546875000000000000E-01_dp)*0.128000000000000000E+03_dp
1166 7545 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1167 7545 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 177))
1168 : END IF
1169 : END IF
1170 : END IF
1171 : END IF
1172 : ELSE
1173 104051 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1174 43577 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
1175 43577 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
1176 43577 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 178))
1177 : ELSE
1178 60474 : IF (X1 <= 0.468750000000000000E-01_dp) THEN
1179 28489 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
1180 20563 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
1181 20563 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
1182 20563 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 179))
1183 : ELSE
1184 7926 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
1185 7926 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1186 7926 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 180))
1187 : END IF
1188 : ELSE
1189 31985 : TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
1190 31985 : TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
1191 31985 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 181))
1192 : END IF
1193 : END IF
1194 : END IF
1195 : END IF
1196 : ELSE
1197 787222 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1198 451057 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
1199 451057 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1200 451057 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 182))
1201 : ELSE
1202 336165 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
1203 140897 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1204 82285 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
1205 82285 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
1206 82285 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 183))
1207 : ELSE
1208 58612 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
1209 58612 : TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
1210 58612 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 184))
1211 : END IF
1212 : ELSE
1213 195268 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
1214 195268 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1215 195268 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 185))
1216 : END IF
1217 : END IF
1218 : END IF
1219 : END IF
1220 : ELSE
1221 3037355 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
1222 1137980 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
1223 1137980 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
1224 1137980 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 186))
1225 : ELSE
1226 1899375 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
1227 904856 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1228 493983 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
1229 493983 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1230 493983 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 187))
1231 : ELSE
1232 410873 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
1233 410873 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1234 410873 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 188))
1235 : END IF
1236 : ELSE
1237 994519 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
1238 994519 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1239 994519 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 189))
1240 : END IF
1241 : END IF
1242 : END IF
1243 : END IF
1244 : ELSE
1245 7024180 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
1246 2884810 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
1247 1419257 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
1248 712256 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
1249 712256 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1250 712256 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 190))
1251 : ELSE
1252 707001 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
1253 707001 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1254 707001 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 191))
1255 : END IF
1256 : ELSE
1257 1465553 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1258 78836 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
1259 78836 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1260 78836 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 192))
1261 : ELSE
1262 1386717 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
1263 1386717 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1264 1386717 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 193))
1265 : END IF
1266 : END IF
1267 : ELSE
1268 4139370 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
1269 1523332 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1270 65829 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
1271 65829 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1272 65829 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 194))
1273 : ELSE
1274 1457503 : IF (X1 <= 0.406250000000000000E+00_dp) THEN
1275 644584 : TG1 = (2*X1 - 0.781250000000000000E+00_dp)*0.320000000000000000E+02_dp
1276 644584 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1277 644584 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 195))
1278 : ELSE
1279 812919 : TG1 = (2*X1 - 0.843750000000000000E+00_dp)*0.320000000000000000E+02_dp
1280 812919 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1281 812919 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 196))
1282 : END IF
1283 : END IF
1284 : ELSE
1285 2616038 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1286 358837 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
1287 358837 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1288 358837 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 197))
1289 : ELSE
1290 2257201 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
1291 2257201 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1292 2257201 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 198))
1293 : END IF
1294 : END IF
1295 : END IF
1296 : END IF
1297 : ELSE
1298 19228765 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
1299 11491024 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
1300 5399639 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1301 221576 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1302 97399 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1303 97399 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1304 97399 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 199))
1305 : ELSE
1306 124177 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1307 124177 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1308 124177 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 200))
1309 : END IF
1310 : ELSE
1311 5178063 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1312 2388359 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1313 2388359 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1314 2388359 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 201))
1315 : ELSE
1316 2789704 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1317 2789704 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1318 2789704 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 202))
1319 : END IF
1320 : END IF
1321 : ELSE
1322 6091385 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
1323 558111 : IF (X1 <= 0.687500000000000000E+00_dp) THEN
1324 202501 : TG1 = (2*X1 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
1325 202501 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1326 202501 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 203))
1327 : ELSE
1328 355610 : TG1 = (2*X1 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
1329 355610 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
1330 355610 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 204))
1331 : END IF
1332 : ELSE
1333 5533274 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1334 5533274 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1335 5533274 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 205))
1336 : END IF
1337 : END IF
1338 : ELSE
1339 7737741 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
1340 4846479 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1341 4846479 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1342 4846479 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 206))
1343 : ELSE
1344 2891262 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1345 2891262 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1346 2891262 : 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 658 : 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 658 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1367 :
1368 658 : patches = 207
1369 658 : IF (Nder > nderiv_max) CPABORT("T_C_G0 init failed")
1370 658 : nderiv_init = Nder
1371 658 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1372 : ! round up to a multiple of 32 to give some generous alignment for each C0
1373 2632 : ALLOCATE (C0(32*((31 + (Nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1374 : ! init mpi'ed buffers to silence warnings under valgrind
1375 102292192 : C0 = 1.0E99_dp
1376 658 : IF (mepos == 0) THEN
1377 329 : ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1378 68432 : DO I = 1, patches
1379 68103 : READ (iunit, *) chunk
1380 50211077 : C0(1:(Nder + 1)*(degree + 1)*(degree + 2)/2, I) = chunk(1:(Nder + 1)*(degree + 1)*(degree + 2)/2)
1381 : END DO
1382 329 : DEALLOCATE (chunk)
1383 : END IF
1384 658 : CALL group%bcast(C0, 0)
1385 :
1386 658 : END SUBROUTINE init
1387 :
1388 : ! **************************************************************************************************
1389 : !> \brief ...
1390 : ! **************************************************************************************************
1391 374 : SUBROUTINE free_C0()
1392 374 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1393 374 : nderiv_init = -1
1394 374 : 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 245491232 : 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 245491232 : T1(0) = 1.0_dp
1415 245491232 : T2(0) = 1.0_dp
1416 245491232 : T1(1) = SQRT2*TG1
1417 245491232 : T2(1) = SQRT2*TG2
1418 245491232 : T1(2) = 2*TG1*T1(1) - SQRT2
1419 245491232 : T2(2) = 2*TG2*T2(1) - SQRT2
1420 245491232 : T1(3) = 2*TG1*T1(2) - T1(1)
1421 245491232 : T2(3) = 2*TG2*T2(2) - T2(1)
1422 245491232 : T1(4) = 2*TG1*T1(3) - T1(2)
1423 245491232 : T2(4) = 2*TG2*T2(3) - T2(2)
1424 245491232 : T1(5) = 2*TG1*T1(4) - T1(3)
1425 245491232 : T2(5) = 2*TG2*T2(4) - T2(3)
1426 245491232 : T1(6) = 2*TG1*T1(5) - T1(4)
1427 245491232 : T2(6) = 2*TG2*T2(5) - T2(4)
1428 245491232 : T1(7) = 2*TG1*T1(6) - T1(5)
1429 245491232 : T2(7) = 2*TG2*T2(6) - T2(5)
1430 245491232 : T1(8) = 2*TG1*T1(7) - T1(6)
1431 245491232 : T2(8) = 2*TG2*T2(7) - T2(6)
1432 245491232 : T1(9) = 2*TG1*T1(8) - T1(7)
1433 245491232 : T2(9) = 2*TG2*T2(8) - T2(7)
1434 245491232 : T1(10) = 2*TG1*T1(9) - T1(8)
1435 245491232 : T2(10) = 2*TG2*T2(9) - T2(8)
1436 245491232 : T1(11) = 2*TG1*T1(10) - T1(9)
1437 245491232 : T2(11) = 2*TG2*T2(10) - T2(9)
1438 245491232 : T1(12) = 2*TG1*T1(11) - T1(10)
1439 245491232 : T2(12) = 2*TG2*T2(11) - T2(10)
1440 245491232 : T1(13) = 2*TG1*T1(12) - T1(11)
1441 245491232 : T2(13) = 2*TG2*T2(12) - T2(11)
1442 852695279 : DO K = 1, NDERIV + 1
1443 607204047 : RES(K) = 0.0_dp
1444 9108060705 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:13), C0(1:14, K))*T2(0)
1445 8500856658 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:12), C0(15:27, K))*T2(1)
1446 7893652611 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:11), C0(28:39, K))*T2(2)
1447 7286448564 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:10), C0(40:50, K))*T2(3)
1448 6679244517 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:9), C0(51:60, K))*T2(4)
1449 6072040470 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:8), C0(61:69, K))*T2(5)
1450 5464836423 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:7), C0(70:77, K))*T2(6)
1451 4857632376 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:6), C0(78:84, K))*T2(7)
1452 4250428329 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:5), C0(85:90, K))*T2(8)
1453 3643224282 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:4), C0(91:95, K))*T2(9)
1454 3036020235 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:3), C0(96:99, K))*T2(10)
1455 2428816188 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:2), C0(100:102, K))*T2(11)
1456 1821612141 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:1), C0(103:104, K))*T2(12)
1457 1459899326 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:0), C0(105:105, K))*T2(13)
1458 : END DO
1459 245491232 : 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 10458706 : FUNCTION get_lmax_init() RESULT(lmax_init)
1468 :
1469 : INTEGER :: lmax_init
1470 :
1471 10458706 : lmax_init = nderiv_init
1472 :
1473 10458706 : END FUNCTION get_lmax_init
1474 :
1475 : END MODULE t_c_g0
|
