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 262006541 : 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 262006541 : use_gamma = .FALSE.
96 262006541 : upper = R**2 + 11.0_dp*R + 50.0_dp
97 262006541 : lower = R**2 - 11.0_dp*R + 0.0_dp
98 262006541 : IF (T > upper) THEN
99 7876524 : RES(1:NDERIV + 1) = 0.0_dp
100 10165955 : RETURN
101 : END IF
102 259956418 : IF (R <= 11.0_dp) THEN
103 219540328 : X2 = R/11.0_dp
104 219540328 : upper = R**2 + 11.0_dp*R + 50.0_dp
105 219540328 : lower = 0.0_dp
106 219540328 : X1 = (T - lower)/(upper - lower)
107 219540328 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
108 145206763 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
109 122317150 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
110 96400831 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
111 20281865 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
112 9426966 : IF (X2 <= 0.625000000000000000E-01_dp) THEN
113 1186073 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
114 611574 : 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 569928 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
150 301809 : 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 169578 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
162 129886 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
163 129886 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
164 129886 : 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 8240893 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
208 3987268 : IF (X2 <= 0.937500000000000000E-01_dp) THEN
209 928916 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
210 600112 : IF (X2 <= 0.781250000000000000E-01_dp) THEN
211 272569 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
212 239474 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
213 239474 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
214 239474 : 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 327543 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
222 237202 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
223 237202 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
224 237202 : 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 3058352 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
244 1571033 : IF (X2 <= 0.109375000000000000E+00_dp) THEN
245 727939 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
246 520219 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
247 520219 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
248 520219 : 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 843094 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
256 610374 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
257 610374 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
258 610374 : 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 5777610 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
298 5777610 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
299 5777610 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 32))
300 : ELSE
301 5077289 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
302 5077289 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
303 5077289 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 33))
304 : END IF
305 : END IF
306 : ELSE
307 76118966 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
308 33594559 : IF (X2 <= 0.187500000000000000E+00_dp) THEN
309 20833442 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
310 10179296 : IF (X2 <= 0.156250000000000000E+00_dp) THEN
311 4698962 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
312 2190306 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
313 1419239 : IF (X2 <= 0.140625000000000000E+00_dp) THEN
314 622682 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
315 622682 : TG2 = (2*X2 - 0.265625000000000000E+00_dp)*0.640000000000000000E+02_dp
316 622682 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 34))
317 : ELSE
318 796557 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
319 796557 : TG2 = (2*X2 - 0.296875000000000000E+00_dp)*0.640000000000000000E+02_dp
320 796557 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 35))
321 : END IF
322 : ELSE
323 771067 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
324 771067 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
325 771067 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 36))
326 : END IF
327 : ELSE
328 2508656 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
329 1126806 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
330 1126806 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
331 1126806 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 37))
332 : ELSE
333 1381850 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
334 1381850 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
335 1381850 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 38))
336 : END IF
337 : END IF
338 : ELSE
339 5480334 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
340 3249170 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
341 2182142 : IF (X2 <= 0.171875000000000000E+00_dp) THEN
342 901219 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
343 901219 : TG2 = (2*X2 - 0.328125000000000000E+00_dp)*0.640000000000000000E+02_dp
344 901219 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 39))
345 : ELSE
346 1280923 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
347 1280923 : TG2 = (2*X2 - 0.359375000000000000E+00_dp)*0.640000000000000000E+02_dp
348 1280923 : 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 2231164 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
357 2231164 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
358 2231164 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 42))
359 : END IF
360 : END IF
361 : ELSE
362 10654146 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
363 5117139 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
364 5117139 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
365 5117139 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 43))
366 : ELSE
367 5537007 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
368 5537007 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
369 5537007 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 44))
370 : END IF
371 : END IF
372 : ELSE
373 12761117 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
374 8360486 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
375 6164450 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
376 3148173 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
377 2170545 : IF (X2 <= 0.203125000000000000E+00_dp) THEN
378 1146774 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
379 1146774 : TG2 = (2*X2 - 0.390625000000000000E+00_dp)*0.640000000000000000E+02_dp
380 1146774 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 45))
381 : ELSE
382 1023771 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
383 1023771 : TG2 = (2*X2 - 0.421875000000000000E+00_dp)*0.640000000000000000E+02_dp
384 1023771 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 46))
385 : END IF
386 : ELSE
387 977628 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
388 977628 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
389 977628 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 47))
390 : END IF
391 : ELSE
392 3016277 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
393 2281195 : IF (X2 <= 0.234375000000000000E+00_dp) THEN
394 1086763 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
395 1086763 : TG2 = (2*X2 - 0.453125000000000000E+00_dp)*0.640000000000000000E+02_dp
396 1086763 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 48))
397 : ELSE
398 1194432 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
399 1194432 : TG2 = (2*X2 - 0.484375000000000000E+00_dp)*0.640000000000000000E+02_dp
400 1194432 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 49))
401 : END IF
402 : ELSE
403 735082 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
404 735082 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
405 735082 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 50))
406 : END IF
407 : END IF
408 : ELSE
409 2196036 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
410 1318160 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
411 1318160 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
412 1318160 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 51))
413 : ELSE
414 877876 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
415 877876 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
416 877876 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 52))
417 : END IF
418 : END IF
419 : ELSE
420 4400631 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
421 2128614 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
422 2128614 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
423 2128614 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 53))
424 : ELSE
425 2272017 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
426 2272017 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
427 2272017 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 54))
428 : END IF
429 : END IF
430 : END IF
431 : ELSE
432 42524407 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
433 19091986 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
434 8733376 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
435 8733376 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
436 8733376 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 55))
437 : ELSE
438 10358610 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
439 10358610 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
440 10358610 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 56))
441 : END IF
442 : ELSE
443 23432421 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
444 23432421 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
445 23432421 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 57))
446 : END IF
447 : END IF
448 : END IF
449 : ELSE
450 25916319 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
451 19060494 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
452 16589512 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
453 14229384 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
454 9239487 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
455 5111249 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
456 4045491 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
457 2085138 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
458 2085138 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
459 2085138 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 58))
460 : ELSE
461 1960353 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
462 1960353 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
463 1960353 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 59))
464 : END IF
465 : ELSE
466 1065758 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
467 573425 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
468 573425 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
469 573425 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 60))
470 : ELSE
471 492333 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
472 492333 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
473 492333 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 61))
474 : END IF
475 : END IF
476 : ELSE
477 4128238 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
478 3472211 : IF (X2 <= 0.343750000000000000E+00_dp) THEN
479 1729021 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
480 1729021 : TG2 = (2*X2 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
481 1729021 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 62))
482 : ELSE
483 1743190 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
484 1743190 : TG2 = (2*X2 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
485 1743190 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 63))
486 : END IF
487 : ELSE
488 656027 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
489 656027 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
490 656027 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 64))
491 : END IF
492 : END IF
493 : ELSE
494 4989897 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
495 4400619 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
496 2729618 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
497 2447056 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
498 2447056 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
499 2447056 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 65))
500 : ELSE
501 282562 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
502 282562 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
503 282562 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 66))
504 : END IF
505 : ELSE
506 1671001 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
507 1542037 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
508 1542037 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
509 1542037 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 67))
510 : ELSE
511 128964 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
512 128964 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
513 128964 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 68))
514 : END IF
515 : END IF
516 : ELSE
517 589278 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
518 360306 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
519 360306 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
520 360306 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 69))
521 : ELSE
522 228972 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
523 228972 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
524 228972 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 70))
525 : END IF
526 : END IF
527 : END IF
528 : ELSE
529 2360128 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
530 1813463 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
531 1179836 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
532 746136 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
533 746136 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
534 746136 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 71))
535 : ELSE
536 433700 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
537 433700 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
538 433700 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 72))
539 : END IF
540 : ELSE
541 633627 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
542 356199 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
543 356199 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
544 356199 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 73))
545 : ELSE
546 277428 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
547 277428 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
548 277428 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 74))
549 : END IF
550 : END IF
551 : ELSE
552 546665 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
553 285802 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
554 189817 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
555 189817 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
556 189817 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 75))
557 : ELSE
558 95985 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
559 95985 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
560 95985 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 76))
561 : END IF
562 : ELSE
563 260863 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
564 156591 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
565 156591 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
566 156591 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 77))
567 : ELSE
568 104272 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
569 104272 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
570 104272 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 78))
571 : END IF
572 : END IF
573 : END IF
574 : END IF
575 : ELSE
576 2470982 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
577 1962239 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
578 855727 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
579 493211 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
580 493211 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
581 493211 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 79))
582 : ELSE
583 362516 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
584 230611 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
585 230611 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
586 230611 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 80))
587 : ELSE
588 131905 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
589 131905 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
590 131905 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 81))
591 : END IF
592 : END IF
593 : ELSE
594 1106512 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
595 881685 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
596 881685 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
597 881685 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 82))
598 : ELSE
599 224827 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
600 224827 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
601 224827 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 83))
602 : END IF
603 : END IF
604 : ELSE
605 508743 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
606 331255 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
607 180233 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
608 180233 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
609 180233 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 84))
610 : ELSE
611 151022 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
612 87610 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
613 87610 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
614 87610 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 85))
615 : ELSE
616 63412 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
617 63412 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
618 63412 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 86))
619 : END IF
620 : END IF
621 : ELSE
622 177488 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
623 105090 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
624 105090 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
625 105090 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 87))
626 : ELSE
627 72398 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
628 72398 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
629 72398 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 88))
630 : END IF
631 : END IF
632 : END IF
633 : END IF
634 : ELSE
635 6855825 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
636 2398044 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
637 1995239 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
638 953666 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
639 953666 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
640 953666 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 89))
641 : ELSE
642 1041573 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
643 1041573 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
644 1041573 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 90))
645 : END IF
646 : ELSE
647 402805 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
648 170444 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
649 75969 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
650 75969 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
651 75969 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 91))
652 : ELSE
653 94475 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
654 94475 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
655 94475 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 92))
656 : END IF
657 : ELSE
658 232361 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
659 232361 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
660 232361 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 93))
661 : END IF
662 : END IF
663 : ELSE
664 4457781 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
665 1952043 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
666 1952043 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
667 1952043 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 94))
668 : ELSE
669 2505738 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
670 2505738 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
671 2505738 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 95))
672 : END IF
673 : END IF
674 : END IF
675 : END IF
676 : ELSE
677 22889613 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
678 13286041 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
679 9630110 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
680 8009789 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
681 7251855 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
682 6778253 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
683 6326672 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
684 4786502 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
685 4786502 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
686 4786502 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 96))
687 : ELSE
688 1540170 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
689 1540170 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
690 1540170 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 97))
691 : END IF
692 : ELSE
693 451581 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
694 346417 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
695 346417 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
696 346417 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 98))
697 : ELSE
698 105164 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
699 105164 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
700 105164 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 99))
701 : END IF
702 : END IF
703 : ELSE
704 473602 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
705 256808 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
706 256808 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
707 256808 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 100))
708 : ELSE
709 216794 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
710 216794 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
711 216794 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 101))
712 : END IF
713 : END IF
714 : ELSE
715 757934 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
716 461256 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
717 272631 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
718 272631 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
719 272631 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 102))
720 : ELSE
721 188625 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
722 188625 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
723 188625 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 103))
724 : END IF
725 : ELSE
726 296678 : IF (X1 <= 0.468750000000000000E-01_dp) THEN
727 140131 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
728 140131 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
729 140131 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 104))
730 : ELSE
731 156547 : TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
732 156547 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
733 156547 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 105))
734 : END IF
735 : END IF
736 : END IF
737 : ELSE
738 1620321 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
739 805375 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
740 343355 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
741 209406 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
742 209406 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
743 209406 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 106))
744 : ELSE
745 133949 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
746 133949 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
747 133949 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 107))
748 : END IF
749 : ELSE
750 462020 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
751 247697 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
752 247697 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
753 247697 : 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 814946 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
762 328868 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
763 328868 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
764 328868 : 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 3655931 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
774 1414404 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
775 442319 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
776 230539 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
777 114043 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
778 114043 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
779 114043 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 112))
780 : ELSE
781 116496 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
782 116496 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
783 116496 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 113))
784 : END IF
785 : ELSE
786 211780 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
787 130026 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
788 130026 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
789 130026 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 114))
790 : ELSE
791 81754 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
792 81754 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
793 81754 : 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 9603572 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
851 4343736 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
852 1572006 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
853 740584 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
854 181971 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
855 90324 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
856 90324 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
857 90324 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 125))
858 : ELSE
859 91647 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
860 91647 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
861 91647 : 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 831422 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
876 178774 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
877 178774 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
878 178774 : 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 2771730 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
893 1357277 : 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 596552 : 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 291526 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
910 291526 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
911 291526 : 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 5259836 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
934 1792216 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
935 858549 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
936 179712 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
937 179712 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
938 179712 : 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 933667 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
946 933667 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
947 933667 : 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 74333565 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
977 56691674 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
978 44458166 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
979 28262411 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
980 22138428 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
981 22138428 : TG2 = (2*X2 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
982 22138428 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 146))
983 : ELSE
984 6123983 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
985 6123983 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
986 6123983 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 147))
987 : END IF
988 : ELSE
989 16195755 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
990 16195755 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
991 16195755 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 148))
992 : END IF
993 : ELSE
994 12233508 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
995 5951539 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
996 2065329 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
997 1027987 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
998 1027987 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
999 1027987 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 149))
1000 : ELSE
1001 1037342 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1002 1037342 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
1003 1037342 : 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 6281969 : IF (X1 <= 0.687500000000000000E+00_dp) THEN
1018 3150214 : TG1 = (2*X1 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
1019 3150214 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1020 3150214 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 153))
1021 : ELSE
1022 3131755 : TG1 = (2*X1 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
1023 3131755 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1024 3131755 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 154))
1025 : END IF
1026 : END IF
1027 : END IF
1028 : ELSE
1029 17641891 : TG1 = (2*X1 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1030 17641891 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1031 17641891 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 155))
1032 : END IF
1033 : END IF
1034 : ELSE
1035 40416090 : IF (T < lower) THEN
1036 8115832 : use_gamma = .TRUE.
1037 8115832 : RETURN
1038 : END IF
1039 32300258 : X2 = 11.0_dp/R
1040 32300258 : X1 = (T - lower)/(upper - lower)
1041 32300258 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
1042 13071493 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
1043 6047313 : 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 5556005 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
1061 2518650 : 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 1325951 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
1079 538729 : 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 195756 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
1097 91705 : 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 55189 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
1115 27061 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
1116 16885 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
1117 8754 : IF (X2 <= 0.953125000000000000E+00_dp) THEN
1118 3460 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
1119 3460 : TG2 = (2*X2 - 0.189062500000000000E+01_dp)*0.640000000000000000E+02_dp
1120 3460 : 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 10176 : 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 6116 : 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 4753 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
1149 4753 : TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
1150 4753 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 174))
1151 : END IF
1152 : END IF
1153 : END IF
1154 : ELSE
1155 28128 : 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 15333 : IF (X1 <= 0.234375000000000000E-01_dp) THEN
1161 7788 : TG1 = (2*X1 - 0.390625000000000000E-01_dp)*0.128000000000000000E+03_dp
1162 7788 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1163 7788 : 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 666 : 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 666 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1367 :
1368 666 : patches = 207
1369 666 : IF (Nder > nderiv_max) CPABORT("T_C_G0 init failed")
1370 666 : nderiv_init = Nder
1371 666 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1372 : ! round up to a multiple of 32 to give some generous alignment for each C0
1373 2664 : ALLOCATE (C0(32*((31 + (Nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1374 : ! init mpi'ed buffers to silence warnings under valgrind
1375 103486176 : C0 = 1.0E99_dp
1376 666 : IF (mepos == 0) THEN
1377 333 : ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1378 69264 : DO I = 1, patches
1379 68931 : READ (iunit, *) chunk
1380 50798754 : C0(1:(Nder + 1)*(degree + 1)*(degree + 2)/2, I) = chunk(1:(Nder + 1)*(degree + 1)*(degree + 2)/2)
1381 : END DO
1382 333 : DEALLOCATE (chunk)
1383 : END IF
1384 666 : CALL group%bcast(C0, 0)
1385 :
1386 666 : 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 251840586 : 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 251840586 : T1(0) = 1.0_dp
1415 251840586 : T2(0) = 1.0_dp
1416 251840586 : T1(1) = SQRT2*TG1
1417 251840586 : T2(1) = SQRT2*TG2
1418 251840586 : T1(2) = 2*TG1*T1(1) - SQRT2
1419 251840586 : T2(2) = 2*TG2*T2(1) - SQRT2
1420 251840586 : T1(3) = 2*TG1*T1(2) - T1(1)
1421 251840586 : T2(3) = 2*TG2*T2(2) - T2(1)
1422 251840586 : T1(4) = 2*TG1*T1(3) - T1(2)
1423 251840586 : T2(4) = 2*TG2*T2(3) - T2(2)
1424 251840586 : T1(5) = 2*TG1*T1(4) - T1(3)
1425 251840586 : T2(5) = 2*TG2*T2(4) - T2(3)
1426 251840586 : T1(6) = 2*TG1*T1(5) - T1(4)
1427 251840586 : T2(6) = 2*TG2*T2(5) - T2(4)
1428 251840586 : T1(7) = 2*TG1*T1(6) - T1(5)
1429 251840586 : T2(7) = 2*TG2*T2(6) - T2(5)
1430 251840586 : T1(8) = 2*TG1*T1(7) - T1(6)
1431 251840586 : T2(8) = 2*TG2*T2(7) - T2(6)
1432 251840586 : T1(9) = 2*TG1*T1(8) - T1(7)
1433 251840586 : T2(9) = 2*TG2*T2(8) - T2(7)
1434 251840586 : T1(10) = 2*TG1*T1(9) - T1(8)
1435 251840586 : T2(10) = 2*TG2*T2(9) - T2(8)
1436 251840586 : T1(11) = 2*TG1*T1(10) - T1(9)
1437 251840586 : T2(11) = 2*TG2*T2(10) - T2(9)
1438 251840586 : T1(12) = 2*TG1*T1(11) - T1(10)
1439 251840586 : T2(12) = 2*TG2*T2(11) - T2(10)
1440 251840586 : T1(13) = 2*TG1*T1(12) - T1(11)
1441 251840586 : T2(13) = 2*TG2*T2(12) - T2(11)
1442 882458912 : DO K = 1, NDERIV + 1
1443 630618326 : RES(K) = 0.0_dp
1444 9459274890 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:13), C0(1:14, K))*T2(0)
1445 8828656564 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:12), C0(15:27, K))*T2(1)
1446 8198038238 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:11), C0(28:39, K))*T2(2)
1447 7567419912 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:10), C0(40:50, K))*T2(3)
1448 6936801586 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:9), C0(51:60, K))*T2(4)
1449 6306183260 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:8), C0(61:69, K))*T2(5)
1450 5675564934 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:7), C0(70:77, K))*T2(6)
1451 5044946608 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:6), C0(78:84, K))*T2(7)
1452 4414328282 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:5), C0(85:90, K))*T2(8)
1453 3783709956 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:4), C0(91:95, K))*T2(9)
1454 3153091630 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:3), C0(96:99, K))*T2(10)
1455 2522473304 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:2), C0(100:102, K))*T2(11)
1456 1891854978 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:1), C0(103:104, K))*T2(12)
1457 1513077238 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:0), C0(105:105, K))*T2(13)
1458 : END DO
1459 251840586 : 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 19576613 : FUNCTION get_lmax_init() RESULT(lmax_init)
1468 :
1469 : INTEGER :: lmax_init
1470 :
1471 19576613 : lmax_init = nderiv_init
1472 :
1473 19576613 : END FUNCTION get_lmax_init
1474 :
1475 : END MODULE t_c_g0
|
