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 263544602 : 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 263544602 : use_gamma = .FALSE.
96 263544602 : upper = R**2 + 11.0_dp*R + 50.0_dp
97 263544602 : lower = R**2 - 11.0_dp*R + 0.0_dp
98 263544602 : IF (T > upper) THEN
99 8343354 : RES(1:NDERIV + 1) = 0.0_dp
100 10314230 : RETURN
101 : END IF
102 261344271 : IF (R <= 11.0_dp) THEN
103 220944892 : X2 = R/11.0_dp
104 220944892 : upper = R**2 + 11.0_dp*R + 50.0_dp
105 220944892 : lower = 0.0_dp
106 220944892 : X1 = (T - lower)/(upper - lower)
107 220944892 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
108 145921060 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
109 123032289 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
110 97109477 : IF (X2 <= 0.125000000000000000E+00_dp) THEN
111 20679533 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
112 9612581 : IF (X2 <= 0.625000000000000000E-01_dp) THEN
113 1216293 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
114 626835 : IF (X2 <= 0.312500000000000000E-01_dp) THEN
115 42453 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
116 16689 : 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 16689 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
128 8141 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
129 8141 : TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
130 8141 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 3))
131 : ELSE
132 8548 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
133 8548 : TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
134 8548 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 4))
135 : END IF
136 : END IF
137 : ELSE
138 25764 : 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 25764 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
144 25764 : TG2 = (2*X2 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
145 25764 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 6))
146 : END IF
147 : END IF
148 : ELSE
149 584382 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
150 307703 : IF (X2 <= 0.468750000000000000E-01_dp) THEN
151 135925 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
152 69268 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
153 69268 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
154 69268 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 7))
155 : ELSE
156 66657 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
157 66657 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
158 66657 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 8))
159 : END IF
160 : ELSE
161 171778 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
162 130828 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
163 130828 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
164 130828 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 9))
165 : ELSE
166 40950 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
167 40950 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
168 40950 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 10))
169 : END IF
170 : END IF
171 : ELSE
172 276679 : IF (X2 <= 0.468750000000000000E-01_dp) THEN
173 211800 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
174 211800 : TG2 = (2*X2 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
175 211800 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 11))
176 : ELSE
177 64879 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
178 64879 : TG2 = (2*X2 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
179 64879 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 12))
180 : END IF
181 : END IF
182 : END IF
183 : ELSE
184 589458 : IF (X2 <= 0.312500000000000000E-01_dp) THEN
185 51073 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
186 16232 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
187 16232 : TG2 = (2*X2 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
188 16232 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 13))
189 : ELSE
190 34841 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
191 34841 : TG2 = (2*X2 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
192 34841 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 14))
193 : END IF
194 : ELSE
195 538385 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
196 291348 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
197 291348 : TG2 = (2*X2 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
198 291348 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 15))
199 : ELSE
200 247037 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
201 247037 : TG2 = (2*X2 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
202 247037 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 16))
203 : END IF
204 : END IF
205 : END IF
206 : ELSE
207 8396288 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
208 4045577 : IF (X2 <= 0.937500000000000000E-01_dp) THEN
209 950508 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
210 608072 : IF (X2 <= 0.781250000000000000E-01_dp) THEN
211 274244 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
212 240414 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
213 240414 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
214 240414 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 17))
215 : ELSE
216 33830 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
217 33830 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
218 33830 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 18))
219 : END IF
220 : ELSE
221 333828 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
222 241530 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
223 241530 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
224 241530 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 19))
225 : ELSE
226 92298 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
227 92298 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
228 92298 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 20))
229 : END IF
230 : END IF
231 : ELSE
232 342436 : IF (X2 <= 0.781250000000000000E-01_dp) THEN
233 67001 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
234 67001 : TG2 = (2*X2 - 0.140625000000000000E+00_dp)*0.640000000000000000E+02_dp
235 67001 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 21))
236 : ELSE
237 275435 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
238 275435 : TG2 = (2*X2 - 0.171875000000000000E+00_dp)*0.640000000000000000E+02_dp
239 275435 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 22))
240 : END IF
241 : END IF
242 : ELSE
243 3095069 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
244 1585144 : IF (X2 <= 0.109375000000000000E+00_dp) THEN
245 738363 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
246 525459 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
247 525459 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
248 525459 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 23))
249 : ELSE
250 212904 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
251 212904 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
252 212904 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 24))
253 : END IF
254 : ELSE
255 846781 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
256 612959 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
257 612959 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
258 612959 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 25))
259 : ELSE
260 233822 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
261 233822 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
262 233822 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 26))
263 : END IF
264 : END IF
265 : ELSE
266 1509925 : IF (X2 <= 0.109375000000000000E+00_dp) THEN
267 633603 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
268 633603 : TG2 = (2*X2 - 0.203125000000000000E+00_dp)*0.640000000000000000E+02_dp
269 633603 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 27))
270 : ELSE
271 876322 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
272 876322 : TG2 = (2*X2 - 0.234375000000000000E+00_dp)*0.640000000000000000E+02_dp
273 876322 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 28))
274 : END IF
275 : END IF
276 : END IF
277 : ELSE
278 4350711 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
279 2098208 : IF (X2 <= 0.937500000000000000E-01_dp) THEN
280 362912 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
281 362912 : TG2 = (2*X2 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
282 362912 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 29))
283 : ELSE
284 1735296 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
285 1735296 : TG2 = (2*X2 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
286 1735296 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 30))
287 : END IF
288 : ELSE
289 2252503 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
290 2252503 : TG2 = (2*X2 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
291 2252503 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 31))
292 : END IF
293 : END IF
294 : END IF
295 : ELSE
296 11066952 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
297 5874836 : TG1 = (2*X1 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
298 5874836 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
299 5874836 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 32))
300 : ELSE
301 5192116 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
302 5192116 : TG2 = (2*X2 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
303 5192116 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 33))
304 : END IF
305 : END IF
306 : ELSE
307 76429944 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
308 33706034 : IF (X2 <= 0.187500000000000000E+00_dp) THEN
309 20907313 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
310 10204922 : IF (X2 <= 0.156250000000000000E+00_dp) THEN
311 4716541 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
312 2199675 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
313 1424540 : IF (X2 <= 0.140625000000000000E+00_dp) THEN
314 625779 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
315 625779 : TG2 = (2*X2 - 0.265625000000000000E+00_dp)*0.640000000000000000E+02_dp
316 625779 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 34))
317 : ELSE
318 798761 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
319 798761 : TG2 = (2*X2 - 0.296875000000000000E+00_dp)*0.640000000000000000E+02_dp
320 798761 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 35))
321 : END IF
322 : ELSE
323 775135 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
324 775135 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
325 775135 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 36))
326 : END IF
327 : ELSE
328 2516866 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
329 1130041 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
330 1130041 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
331 1130041 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 37))
332 : ELSE
333 1386825 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
334 1386825 : TG2 = (2*X2 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
335 1386825 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 38))
336 : END IF
337 : END IF
338 : ELSE
339 5488381 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
340 3255156 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
341 2186097 : IF (X2 <= 0.171875000000000000E+00_dp) THEN
342 903423 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
343 903423 : TG2 = (2*X2 - 0.328125000000000000E+00_dp)*0.640000000000000000E+02_dp
344 903423 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 39))
345 : ELSE
346 1282674 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
347 1282674 : TG2 = (2*X2 - 0.359375000000000000E+00_dp)*0.640000000000000000E+02_dp
348 1282674 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 40))
349 : END IF
350 : ELSE
351 1069059 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
352 1069059 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
353 1069059 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 41))
354 : END IF
355 : ELSE
356 2233225 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
357 2233225 : TG2 = (2*X2 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
358 2233225 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 42))
359 : END IF
360 : END IF
361 : ELSE
362 10702391 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
363 5133779 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
364 5133779 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
365 5133779 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 43))
366 : ELSE
367 5568612 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
368 5568612 : TG2 = (2*X2 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
369 5568612 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 44))
370 : END IF
371 : END IF
372 : ELSE
373 12798721 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
374 8381905 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
375 6175399 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
376 3152724 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
377 2173773 : IF (X2 <= 0.203125000000000000E+00_dp) THEN
378 1148165 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
379 1148165 : TG2 = (2*X2 - 0.390625000000000000E+00_dp)*0.640000000000000000E+02_dp
380 1148165 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 45))
381 : ELSE
382 1025608 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
383 1025608 : TG2 = (2*X2 - 0.421875000000000000E+00_dp)*0.640000000000000000E+02_dp
384 1025608 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 46))
385 : END IF
386 : ELSE
387 978951 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
388 978951 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
389 978951 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 47))
390 : END IF
391 : ELSE
392 3022675 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
393 2283520 : IF (X2 <= 0.234375000000000000E+00_dp) THEN
394 1088663 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
395 1088663 : TG2 = (2*X2 - 0.453125000000000000E+00_dp)*0.640000000000000000E+02_dp
396 1088663 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 48))
397 : ELSE
398 1194857 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
399 1194857 : TG2 = (2*X2 - 0.484375000000000000E+00_dp)*0.640000000000000000E+02_dp
400 1194857 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 49))
401 : END IF
402 : ELSE
403 739155 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
404 739155 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
405 739155 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 50))
406 : END IF
407 : END IF
408 : ELSE
409 2206506 : IF (X2 <= 0.218750000000000000E+00_dp) THEN
410 1322323 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
411 1322323 : TG2 = (2*X2 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
412 1322323 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 51))
413 : ELSE
414 884183 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
415 884183 : TG2 = (2*X2 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
416 884183 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 52))
417 : END IF
418 : END IF
419 : ELSE
420 4416816 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
421 2135059 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
422 2135059 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
423 2135059 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 53))
424 : ELSE
425 2281757 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
426 2281757 : TG2 = (2*X2 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
427 2281757 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 54))
428 : END IF
429 : END IF
430 : END IF
431 : ELSE
432 42723910 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
433 19189052 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
434 8781055 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
435 8781055 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
436 8781055 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 55))
437 : ELSE
438 10407997 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
439 10407997 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
440 10407997 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 56))
441 : END IF
442 : ELSE
443 23534858 : TG1 = (2*X1 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
444 23534858 : TG2 = (2*X2 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
445 23534858 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 57))
446 : END IF
447 : END IF
448 : END IF
449 : ELSE
450 25922812 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
451 19066204 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
452 16589767 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
453 14226771 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
454 9237548 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
455 5109977 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
456 4043075 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
457 2083210 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
458 2083210 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
459 2083210 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 58))
460 : ELSE
461 1959865 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
462 1959865 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
463 1959865 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 59))
464 : END IF
465 : ELSE
466 1066902 : IF (X2 <= 0.281250000000000000E+00_dp) THEN
467 573859 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
468 573859 : TG2 = (2*X2 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
469 573859 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 60))
470 : ELSE
471 493043 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
472 493043 : TG2 = (2*X2 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
473 493043 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 61))
474 : END IF
475 : END IF
476 : ELSE
477 4127571 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
478 3471192 : IF (X2 <= 0.343750000000000000E+00_dp) THEN
479 1728769 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
480 1728769 : TG2 = (2*X2 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
481 1728769 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 62))
482 : ELSE
483 1742423 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
484 1742423 : TG2 = (2*X2 - 0.718750000000000000E+00_dp)*0.320000000000000000E+02_dp
485 1742423 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 63))
486 : END IF
487 : ELSE
488 656379 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
489 656379 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
490 656379 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 64))
491 : END IF
492 : END IF
493 : ELSE
494 4989223 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
495 4398807 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
496 2729196 : 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 282140 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
502 282140 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
503 282140 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 66))
504 : END IF
505 : ELSE
506 1669611 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
507 1540952 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
508 1540952 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
509 1540952 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 67))
510 : ELSE
511 128659 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
512 128659 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
513 128659 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 68))
514 : END IF
515 : END IF
516 : ELSE
517 590416 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
518 361394 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
519 361394 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
520 361394 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 69))
521 : ELSE
522 229022 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
523 229022 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
524 229022 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 70))
525 : END IF
526 : END IF
527 : END IF
528 : ELSE
529 2362996 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
530 1815005 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
531 1180875 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
532 746918 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
533 746918 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
534 746918 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 71))
535 : ELSE
536 433957 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
537 433957 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
538 433957 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 72))
539 : END IF
540 : ELSE
541 634130 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
542 356734 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
543 356734 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
544 356734 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 73))
545 : ELSE
546 277396 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
547 277396 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
548 277396 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 74))
549 : END IF
550 : END IF
551 : ELSE
552 547991 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
553 286368 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
554 190233 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
555 190233 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
556 190233 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 75))
557 : ELSE
558 96135 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
559 96135 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
560 96135 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 76))
561 : END IF
562 : ELSE
563 261623 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
564 157275 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
565 157275 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
566 157275 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 77))
567 : ELSE
568 104348 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
569 104348 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
570 104348 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 78))
571 : END IF
572 : END IF
573 : END IF
574 : END IF
575 : ELSE
576 2476437 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
577 1964767 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
578 856696 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
579 494427 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
580 494427 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
581 494427 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 79))
582 : ELSE
583 362269 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
584 230672 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
585 230672 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
586 230672 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 80))
587 : ELSE
588 131597 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
589 131597 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
590 131597 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 81))
591 : END IF
592 : END IF
593 : ELSE
594 1108071 : IF (X2 <= 0.312500000000000000E+00_dp) THEN
595 883324 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
596 883324 : TG2 = (2*X2 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
597 883324 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 82))
598 : ELSE
599 224747 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
600 224747 : TG2 = (2*X2 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
601 224747 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 83))
602 : END IF
603 : END IF
604 : ELSE
605 511670 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
606 333602 : IF (X2 <= 0.437500000000000000E+00_dp) THEN
607 182308 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
608 182308 : TG2 = (2*X2 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
609 182308 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 84))
610 : ELSE
611 151294 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
612 87714 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
613 87714 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
614 87714 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 85))
615 : ELSE
616 63580 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
617 63580 : TG2 = (2*X2 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
618 63580 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 86))
619 : END IF
620 : END IF
621 : ELSE
622 178068 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
623 105366 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
624 105366 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
625 105366 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 87))
626 : ELSE
627 72702 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
628 72702 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
629 72702 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 88))
630 : END IF
631 : END IF
632 : END IF
633 : END IF
634 : ELSE
635 6856608 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
636 2399517 : IF (X2 <= 0.375000000000000000E+00_dp) THEN
637 1995771 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
638 953190 : TG1 = (2*X1 - 0.562500000000000000E+00_dp)*0.160000000000000000E+02_dp
639 953190 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
640 953190 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 89))
641 : ELSE
642 1042581 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
643 1042581 : TG2 = (2*X2 - 0.625000000000000000E+00_dp)*0.800000000000000000E+01_dp
644 1042581 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 90))
645 : END IF
646 : ELSE
647 403746 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
648 170781 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
649 76098 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
650 76098 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
651 76098 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 91))
652 : ELSE
653 94683 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
654 94683 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
655 94683 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 92))
656 : END IF
657 : ELSE
658 232965 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
659 232965 : TG2 = (2*X2 - 0.875000000000000000E+00_dp)*0.800000000000000000E+01_dp
660 232965 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 93))
661 : END IF
662 : END IF
663 : ELSE
664 4457091 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
665 1951091 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
666 1951091 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
667 1951091 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 94))
668 : ELSE
669 2506000 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
670 2506000 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
671 2506000 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 95))
672 : END IF
673 : END IF
674 : END IF
675 : END IF
676 : ELSE
677 22888771 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
678 13284560 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
679 9627415 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
680 8004957 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
681 7248084 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
682 6776911 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
683 6325664 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
684 4786730 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
685 4786730 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
686 4786730 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 96))
687 : ELSE
688 1538934 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
689 1538934 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
690 1538934 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 97))
691 : END IF
692 : ELSE
693 451247 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
694 346229 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
695 346229 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
696 346229 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 98))
697 : ELSE
698 105018 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
699 105018 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
700 105018 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 99))
701 : END IF
702 : END IF
703 : ELSE
704 471173 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
705 255175 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
706 255175 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
707 255175 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 100))
708 : ELSE
709 215998 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
710 215998 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
711 215998 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 101))
712 : END IF
713 : END IF
714 : ELSE
715 756873 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
716 461127 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
717 272748 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
718 272748 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
719 272748 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 102))
720 : ELSE
721 188379 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
722 188379 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
723 188379 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 103))
724 : END IF
725 : ELSE
726 295746 : IF (X1 <= 0.468750000000000000E-01_dp) THEN
727 139293 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
728 139293 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
729 139293 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 104))
730 : ELSE
731 156453 : TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
732 156453 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
733 156453 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 105))
734 : END IF
735 : END IF
736 : END IF
737 : ELSE
738 1622458 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
739 807538 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
740 343984 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
741 209729 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
742 209729 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
743 209729 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 106))
744 : ELSE
745 134255 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
746 134255 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
747 134255 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 107))
748 : END IF
749 : ELSE
750 463554 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
751 248181 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
752 248181 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
753 248181 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 108))
754 : ELSE
755 215373 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
756 215373 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
757 215373 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 109))
758 : END IF
759 : END IF
760 : ELSE
761 814920 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
762 328864 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
763 328864 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
764 328864 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 110))
765 : ELSE
766 486056 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
767 486056 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
768 486056 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 111))
769 : END IF
770 : END IF
771 : END IF
772 : ELSE
773 3657145 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
774 1415870 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
775 443403 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
776 230749 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
777 114065 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
778 114065 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
779 114065 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 112))
780 : ELSE
781 116684 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
782 116684 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
783 116684 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 113))
784 : END IF
785 : ELSE
786 212654 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
787 130698 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
788 130698 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
789 130698 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 114))
790 : ELSE
791 81956 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
792 81956 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
793 81956 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 115))
794 : END IF
795 : END IF
796 : ELSE
797 972467 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
798 424908 : IF (X1 <= 0.156250000000000000E+00_dp) THEN
799 184861 : TG1 = (2*X1 - 0.281250000000000000E+00_dp)*0.320000000000000000E+02_dp
800 184861 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
801 184861 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 116))
802 : ELSE
803 240047 : TG1 = (2*X1 - 0.343750000000000000E+00_dp)*0.320000000000000000E+02_dp
804 240047 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
805 240047 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 117))
806 : END IF
807 : ELSE
808 547559 : IF (X1 <= 0.218750000000000000E+00_dp) THEN
809 270791 : TG1 = (2*X1 - 0.406250000000000000E+00_dp)*0.320000000000000000E+02_dp
810 270791 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
811 270791 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 118))
812 : ELSE
813 276768 : TG1 = (2*X1 - 0.468750000000000000E+00_dp)*0.320000000000000000E+02_dp
814 276768 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
815 276768 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 119))
816 : END IF
817 : END IF
818 : END IF
819 : ELSE
820 2241275 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
821 1017115 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
822 535188 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
823 535188 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
824 535188 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 120))
825 : ELSE
826 481927 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
827 481927 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
828 481927 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 121))
829 : END IF
830 : ELSE
831 1224160 : 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 536702 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
843 536702 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
844 536702 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 124))
845 : END IF
846 : END IF
847 : END IF
848 : END IF
849 : ELSE
850 9604211 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
851 4345103 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
852 1572202 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
853 740780 : 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 558809 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
865 276152 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
866 276152 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
867 276152 : 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 2772901 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
893 1357658 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
894 761286 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
895 371281 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
896 371281 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
897 371281 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 132))
898 : ELSE
899 390005 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
900 390005 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
901 390005 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 133))
902 : END IF
903 : ELSE
904 596372 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
905 304946 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
906 304946 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
907 304946 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 134))
908 : ELSE
909 291426 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
910 291426 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
911 291426 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 135))
912 : END IF
913 : END IF
914 : ELSE
915 1415243 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
916 802228 : IF (X1 <= 0.343750000000000000E+00_dp) THEN
917 352134 : TG1 = (2*X1 - 0.656250000000000000E+00_dp)*0.320000000000000000E+02_dp
918 352134 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
919 352134 : 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 613015 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
927 613015 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
928 613015 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 138))
929 : END IF
930 : END IF
931 : END IF
932 : ELSE
933 5259108 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
934 1792224 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
935 858553 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
936 179716 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
937 179716 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
938 179716 : 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 933671 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
946 933671 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
947 933671 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 141))
948 : END IF
949 : ELSE
950 3466884 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
951 1690248 : 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 741250 : TG1 = (2*X1 - 0.812500000000000000E+00_dp)*0.160000000000000000E+02_dp
957 741250 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
958 741250 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 143))
959 : END IF
960 : ELSE
961 1776636 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
962 935582 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
963 935582 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
964 935582 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 144))
965 : ELSE
966 841054 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
967 841054 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
968 841054 : 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 75023832 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
977 57086634 : IF (X2 <= 0.500000000000000000E+00_dp) THEN
978 44855802 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
979 28482509 : IF (X2 <= 0.250000000000000000E+00_dp) THEN
980 22353445 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
981 22353445 : TG2 = (2*X2 - 0.250000000000000000E+00_dp)*0.400000000000000000E+01_dp
982 22353445 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 146))
983 : ELSE
984 6129064 : TG1 = (2*X1 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
985 6129064 : TG2 = (2*X2 - 0.750000000000000000E+00_dp)*0.400000000000000000E+01_dp
986 6129064 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 147))
987 : END IF
988 : ELSE
989 16373293 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
990 16373293 : TG2 = (2*X2 - 0.500000000000000000E+00_dp)*0.200000000000000000E+01_dp
991 16373293 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 148))
992 : END IF
993 : ELSE
994 12230832 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
995 5949303 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
996 2065333 : 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 1037346 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1002 1037346 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
1003 1037346 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 150))
1004 : END IF
1005 : ELSE
1006 3883970 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1007 1897806 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1008 1897806 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1009 1897806 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 151))
1010 : ELSE
1011 1986164 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1012 1986164 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1013 1986164 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 152))
1014 : END IF
1015 : END IF
1016 : ELSE
1017 6281529 : IF (X1 <= 0.687500000000000000E+00_dp) THEN
1018 3149970 : TG1 = (2*X1 - 0.131250000000000000E+01_dp)*0.160000000000000000E+02_dp
1019 3149970 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1020 3149970 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 153))
1021 : ELSE
1022 3131559 : TG1 = (2*X1 - 0.143750000000000000E+01_dp)*0.160000000000000000E+02_dp
1023 3131559 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1024 3131559 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 154))
1025 : END IF
1026 : END IF
1027 : END IF
1028 : ELSE
1029 17937198 : TG1 = (2*X1 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1030 17937198 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1031 17937198 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 155))
1032 : END IF
1033 : END IF
1034 : ELSE
1035 40399379 : IF (T < lower) THEN
1036 8113899 : use_gamma = .TRUE.
1037 8113899 : RETURN
1038 : END IF
1039 32285480 : X2 = 11.0_dp/R
1040 32285480 : X1 = (T - lower)/(upper - lower)
1041 32285480 : IF (X1 <= 0.500000000000000000E+00_dp) THEN
1042 13064869 : IF (X1 <= 0.250000000000000000E+00_dp) THEN
1043 6044567 : 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 5553259 : IF (X1 <= 0.125000000000000000E+00_dp) THEN
1061 2517458 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
1062 1192461 : IF (X2 <= 0.625000000000000000E+00_dp) THEN
1063 479206 : TG1 = (2*X1 - 0.125000000000000000E+00_dp)*0.800000000000000000E+01_dp
1064 479206 : TG2 = (2*X2 - 0.112500000000000000E+01_dp)*0.800000000000000000E+01_dp
1065 479206 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 159))
1066 : ELSE
1067 713255 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
1068 316758 : TG1 = (2*X1 - 0.625000000000000000E-01_dp)*0.160000000000000000E+02_dp
1069 316758 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1070 316758 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 160))
1071 : ELSE
1072 396497 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
1073 396497 : TG2 = (2*X2 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1074 396497 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 161))
1075 : END IF
1076 : END IF
1077 : ELSE
1078 1324997 : IF (X1 <= 0.625000000000000000E-01_dp) THEN
1079 538235 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1080 342833 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
1081 152887 : 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 46439 : TG1 = (2*X1 - 0.312500000000000000E-01_dp)*0.320000000000000000E+02_dp
1087 46439 : TG2 = (2*X2 - 0.168750000000000000E+01_dp)*0.160000000000000000E+02_dp
1088 46439 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 163))
1089 : END IF
1090 : ELSE
1091 189946 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
1092 189946 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1093 189946 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 164))
1094 : END IF
1095 : ELSE
1096 195402 : IF (X1 <= 0.312500000000000000E-01_dp) THEN
1097 91443 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1098 36408 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
1099 14212 : IF (X2 <= 0.906250000000000000E+00_dp) THEN
1100 6092 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
1101 6092 : TG2 = (2*X2 - 0.178125000000000000E+01_dp)*0.320000000000000000E+02_dp
1102 6092 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 165))
1103 : ELSE
1104 8120 : TG1 = (2*X1 - 0.156250000000000000E-01_dp)*0.640000000000000000E+02_dp
1105 8120 : TG2 = (2*X2 - 0.184375000000000000E+01_dp)*0.320000000000000000E+02_dp
1106 8120 : 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 55035 : IF (X1 <= 0.156250000000000000E-01_dp) THEN
1115 27043 : 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 10158 : IF (X1 <= 0.781250000000000000E-02_dp) THEN
1133 4048 : 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 2683 : TG1 = (2*X1 - 0.781250000000000000E-02_dp)*0.128000000000000000E+03_dp
1139 2683 : TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
1140 2683 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 172))
1141 : END IF
1142 : ELSE
1143 6110 : 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 4747 : TG1 = (2*X1 - 0.234375000000000000E-01_dp)*0.128000000000000000E+03_dp
1149 4747 : TG2 = (2*X2 - 0.198437500000000000E+01_dp)*0.640000000000000000E+02_dp
1150 4747 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 174))
1151 : END IF
1152 : END IF
1153 : END IF
1154 : ELSE
1155 27992 : IF (X2 <= 0.968750000000000000E+00_dp) THEN
1156 12699 : TG1 = (2*X1 - 0.468750000000000000E-01_dp)*0.640000000000000000E+02_dp
1157 12699 : TG2 = (2*X2 - 0.190625000000000000E+01_dp)*0.320000000000000000E+02_dp
1158 12699 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 175))
1159 : ELSE
1160 15293 : IF (X1 <= 0.234375000000000000E-01_dp) THEN
1161 7764 : TG1 = (2*X1 - 0.390625000000000000E-01_dp)*0.128000000000000000E+03_dp
1162 7764 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1163 7764 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 176))
1164 : ELSE
1165 7529 : TG1 = (2*X1 - 0.546875000000000000E-01_dp)*0.128000000000000000E+03_dp
1166 7529 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1167 7529 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 177))
1168 : END IF
1169 : END IF
1170 : END IF
1171 : END IF
1172 : ELSE
1173 103959 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1174 43513 : TG1 = (2*X1 - 0.937500000000000000E-01_dp)*0.320000000000000000E+02_dp
1175 43513 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
1176 43513 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 178))
1177 : ELSE
1178 60446 : IF (X1 <= 0.468750000000000000E-01_dp) THEN
1179 28477 : 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 7914 : TG1 = (2*X1 - 0.781250000000000000E-01_dp)*0.640000000000000000E+02_dp
1185 7914 : TG2 = (2*X2 - 0.196875000000000000E+01_dp)*0.320000000000000000E+02_dp
1186 7914 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 180))
1187 : END IF
1188 : ELSE
1189 31969 : TG1 = (2*X1 - 0.109375000000000000E+00_dp)*0.640000000000000000E+02_dp
1190 31969 : TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
1191 31969 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 181))
1192 : END IF
1193 : END IF
1194 : END IF
1195 : END IF
1196 : ELSE
1197 786762 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1198 450753 : TG1 = (2*X1 - 0.187500000000000000E+00_dp)*0.160000000000000000E+02_dp
1199 450753 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1200 450753 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 182))
1201 : ELSE
1202 336009 : IF (X1 <= 0.937500000000000000E-01_dp) THEN
1203 140837 : IF (X2 <= 0.937500000000000000E+00_dp) THEN
1204 82253 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
1205 82253 : TG2 = (2*X2 - 0.181250000000000000E+01_dp)*0.160000000000000000E+02_dp
1206 82253 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 183))
1207 : ELSE
1208 58584 : TG1 = (2*X1 - 0.156250000000000000E+00_dp)*0.320000000000000000E+02_dp
1209 58584 : TG2 = (2*X2 - 0.193750000000000000E+01_dp)*0.160000000000000000E+02_dp
1210 58584 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 184))
1211 : END IF
1212 : ELSE
1213 195172 : TG1 = (2*X1 - 0.218750000000000000E+00_dp)*0.320000000000000000E+02_dp
1214 195172 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1215 195172 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 185))
1216 : END IF
1217 : END IF
1218 : END IF
1219 : END IF
1220 : ELSE
1221 3035801 : IF (X2 <= 0.750000000000000000E+00_dp) THEN
1222 1137598 : TG1 = (2*X1 - 0.375000000000000000E+00_dp)*0.800000000000000000E+01_dp
1223 1137598 : TG2 = (2*X2 - 0.125000000000000000E+01_dp)*0.400000000000000000E+01_dp
1224 1137598 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 186))
1225 : ELSE
1226 1898203 : IF (X1 <= 0.187500000000000000E+00_dp) THEN
1227 904180 : IF (X2 <= 0.875000000000000000E+00_dp) THEN
1228 493739 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
1229 493739 : TG2 = (2*X2 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1230 493739 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 187))
1231 : ELSE
1232 410441 : TG1 = (2*X1 - 0.312500000000000000E+00_dp)*0.160000000000000000E+02_dp
1233 410441 : TG2 = (2*X2 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1234 410441 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 188))
1235 : END IF
1236 : ELSE
1237 994023 : TG1 = (2*X1 - 0.437500000000000000E+00_dp)*0.160000000000000000E+02_dp
1238 994023 : TG2 = (2*X2 - 0.175000000000000000E+01_dp)*0.400000000000000000E+01_dp
1239 994023 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 189))
1240 : END IF
1241 : END IF
1242 : END IF
1243 : END IF
1244 : ELSE
1245 7020302 : IF (X1 <= 0.375000000000000000E+00_dp) THEN
1246 2883292 : IF (X1 <= 0.312500000000000000E+00_dp) THEN
1247 1418399 : IF (X1 <= 0.281250000000000000E+00_dp) THEN
1248 711948 : TG1 = (2*X1 - 0.531250000000000000E+00_dp)*0.320000000000000000E+02_dp
1249 711948 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1250 711948 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 190))
1251 : ELSE
1252 706451 : TG1 = (2*X1 - 0.593750000000000000E+00_dp)*0.320000000000000000E+02_dp
1253 706451 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1254 706451 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 191))
1255 : END IF
1256 : ELSE
1257 1464893 : 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 1386057 : TG1 = (2*X1 - 0.687500000000000000E+00_dp)*0.160000000000000000E+02_dp
1263 1386057 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1264 1386057 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 193))
1265 : END IF
1266 : END IF
1267 : ELSE
1268 4137010 : IF (X1 <= 0.437500000000000000E+00_dp) THEN
1269 1522046 : 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 1456217 : IF (X1 <= 0.406250000000000000E+00_dp) THEN
1275 643918 : TG1 = (2*X1 - 0.781250000000000000E+00_dp)*0.320000000000000000E+02_dp
1276 643918 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1277 643918 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 195))
1278 : ELSE
1279 812299 : TG1 = (2*X1 - 0.843750000000000000E+00_dp)*0.320000000000000000E+02_dp
1280 812299 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1281 812299 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 196))
1282 : END IF
1283 : END IF
1284 : ELSE
1285 2614964 : 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 2256127 : TG1 = (2*X1 - 0.937500000000000000E+00_dp)*0.160000000000000000E+02_dp
1291 2256127 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1292 2256127 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 198))
1293 : END IF
1294 : END IF
1295 : END IF
1296 : END IF
1297 : ELSE
1298 19220611 : IF (X1 <= 0.750000000000000000E+00_dp) THEN
1299 11486132 : IF (X1 <= 0.625000000000000000E+00_dp) THEN
1300 5397403 : 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 5175827 : IF (X1 <= 0.562500000000000000E+00_dp) THEN
1312 2387475 : TG1 = (2*X1 - 0.106250000000000000E+01_dp)*0.160000000000000000E+02_dp
1313 2387475 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1314 2387475 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 201))
1315 : ELSE
1316 2788352 : TG1 = (2*X1 - 0.118750000000000000E+01_dp)*0.160000000000000000E+02_dp
1317 2788352 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1318 2788352 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 202))
1319 : END IF
1320 : END IF
1321 : ELSE
1322 6088729 : 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 5530618 : TG1 = (2*X1 - 0.137500000000000000E+01_dp)*0.800000000000000000E+01_dp
1334 5530618 : TG2 = (2*X2 - 0.150000000000000000E+01_dp)*0.200000000000000000E+01_dp
1335 5530618 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 205))
1336 : END IF
1337 : END IF
1338 : ELSE
1339 7734479 : IF (X1 <= 0.875000000000000000E+00_dp) THEN
1340 4843775 : TG1 = (2*X1 - 0.162500000000000000E+01_dp)*0.800000000000000000E+01_dp
1341 4843775 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1342 4843775 : CALL PD2VAL(RES, NDERIV, TG1, TG2, C0(1, 206))
1343 : ELSE
1344 2890704 : TG1 = (2*X1 - 0.187500000000000000E+01_dp)*0.800000000000000000E+01_dp
1345 2890704 : TG2 = (2*X2 - 0.100000000000000000E+01_dp)*0.100000000000000000E+01_dp
1346 2890704 : 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 668 : 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 668 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1367 :
1368 668 : patches = 207
1369 668 : IF (Nder > nderiv_max) CPABORT("T_C_G0 init failed")
1370 668 : nderiv_init = Nder
1371 668 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1372 : ! round up to a multiple of 32 to give some generous alignment for each C0
1373 2672 : ALLOCATE (C0(32*((31 + (Nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1374 : ! init mpi'ed buffers to silence warnings under valgrind
1375 103711808 : C0 = 1.0E99_dp
1376 668 : IF (mepos == 0) THEN
1377 334 : ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1378 69472 : DO I = 1, patches
1379 69138 : READ (iunit, *) chunk
1380 50907637 : C0(1:(Nder + 1)*(degree + 1)*(degree + 2)/2, I) = chunk(1:(Nder + 1)*(degree + 1)*(degree + 2)/2)
1381 : END DO
1382 334 : DEALLOCATE (chunk)
1383 : END IF
1384 668 : CALL group%bcast(C0, 0)
1385 :
1386 668 : END SUBROUTINE init
1387 :
1388 : ! **************************************************************************************************
1389 : !> \brief ...
1390 : ! **************************************************************************************************
1391 376 : SUBROUTINE free_C0()
1392 376 : IF (ALLOCATED(C0)) DEALLOCATE (C0)
1393 376 : nderiv_init = -1
1394 376 : END SUBROUTINE free_C0
1395 :
1396 : ! **************************************************************************************************
1397 : !> \brief ...
1398 : !> \param RES ...
1399 : !> \param NDERIV ...
1400 : !> \param TG1 ...
1401 : !> \param TG2 ...
1402 : !> \param C0 ...
1403 : ! **************************************************************************************************
1404 253230372 : 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 253230372 : T1(0) = 1.0_dp
1415 253230372 : T2(0) = 1.0_dp
1416 253230372 : T1(1) = SQRT2*TG1
1417 253230372 : T2(1) = SQRT2*TG2
1418 253230372 : T1(2) = 2*TG1*T1(1) - SQRT2
1419 253230372 : T2(2) = 2*TG2*T2(1) - SQRT2
1420 253230372 : T1(3) = 2*TG1*T1(2) - T1(1)
1421 253230372 : T2(3) = 2*TG2*T2(2) - T2(1)
1422 253230372 : T1(4) = 2*TG1*T1(3) - T1(2)
1423 253230372 : T2(4) = 2*TG2*T2(3) - T2(2)
1424 253230372 : T1(5) = 2*TG1*T1(4) - T1(3)
1425 253230372 : T2(5) = 2*TG2*T2(4) - T2(3)
1426 253230372 : T1(6) = 2*TG1*T1(5) - T1(4)
1427 253230372 : T2(6) = 2*TG2*T2(5) - T2(4)
1428 253230372 : T1(7) = 2*TG1*T1(6) - T1(5)
1429 253230372 : T2(7) = 2*TG2*T2(6) - T2(5)
1430 253230372 : T1(8) = 2*TG1*T1(7) - T1(6)
1431 253230372 : T2(8) = 2*TG2*T2(7) - T2(6)
1432 253230372 : T1(9) = 2*TG1*T1(8) - T1(7)
1433 253230372 : T2(9) = 2*TG2*T2(8) - T2(7)
1434 253230372 : T1(10) = 2*TG1*T1(9) - T1(8)
1435 253230372 : T2(10) = 2*TG2*T2(9) - T2(8)
1436 253230372 : T1(11) = 2*TG1*T1(10) - T1(9)
1437 253230372 : T2(11) = 2*TG2*T2(10) - T2(9)
1438 253230372 : T1(12) = 2*TG1*T1(11) - T1(10)
1439 253230372 : T2(12) = 2*TG2*T2(11) - T2(10)
1440 253230372 : T1(13) = 2*TG1*T1(12) - T1(11)
1441 253230372 : T2(13) = 2*TG2*T2(12) - T2(11)
1442 886821192 : DO K = 1, NDERIV + 1
1443 633590820 : RES(K) = 0.0_dp
1444 9503862300 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:13), C0(1:14, K))*T2(0)
1445 8870271480 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:12), C0(15:27, K))*T2(1)
1446 8236680660 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:11), C0(28:39, K))*T2(2)
1447 7603089840 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:10), C0(40:50, K))*T2(3)
1448 6969499020 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:9), C0(51:60, K))*T2(4)
1449 6335908200 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:8), C0(61:69, K))*T2(5)
1450 5702317380 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:7), C0(70:77, K))*T2(6)
1451 5068726560 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:6), C0(78:84, K))*T2(7)
1452 4435135740 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:5), C0(85:90, K))*T2(8)
1453 3801544920 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:4), C0(91:95, K))*T2(9)
1454 3167954100 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:3), C0(96:99, K))*T2(10)
1455 2534363280 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:2), C0(100:102, K))*T2(11)
1456 1900772460 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:1), C0(103:104, K))*T2(12)
1457 1520412012 : RES(K) = RES(K) + DOT_PRODUCT(T1(0:0), C0(105:105, K))*T2(13)
1458 : END DO
1459 253230372 : 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 20981190 : FUNCTION get_lmax_init() RESULT(lmax_init)
1468 :
1469 : INTEGER :: lmax_init
1470 :
1471 20981190 : lmax_init = nderiv_init
1472 :
1473 20981190 : END FUNCTION get_lmax_init
1474 :
1475 : END MODULE t_c_g0
|
