(git:34ef472)
t_c_g0.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 !--------------------------------------------------------------------------------------------------!
9 ! Copyright (c) 2008, 2009, Joost VandeVondele and Manuel Guidon !
10 ! All rights reserved. !
11 ! !
12 ! Redistribution and use in source and binary forms, with or without !
13 ! modification, are permitted provided that the following conditions are met: !
14 ! * Redistributions of source code must retain the above copyright !
15 ! notice, this list of conditions and the following disclaimer. !
16 ! * Redistributions in binary form must reproduce the above copyright !
17 ! notice, this list of conditions and the following disclaimer in the !
18 ! documentation and/or other materials provided with the distribution. !
19 ! !
20 ! THIS SOFTWARE IS PROVIDED BY Joost VandeVondele and Manuel Guidon AS IS AND ANY !
21 ! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED !
22 ! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE !
23 ! DISCLAIMED. IN NO EVENT SHALL Joost VandeVondele or Manuel Guidon BE LIABLE FOR ANY !
24 ! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES !
25 ! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; !
26 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND !
27 ! ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT !
28 ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS !
29 ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. !
30 !--------------------------------------------------------------------------------------------------!
31 
32 ! **************************************************************************************************
33 !> \brief This module computes the basic integrals for the truncated coulomb operator
34 !>
35 !> res(1) =G_0(R,T)= ((2*erf(sqrt(t))+erf(R-sqrt(t))-erf(R+sqrt(t)))/sqrt(t))
36 !>
37 !> and up to 21 derivatives with respect to T
38 !>
39 !> res(n+1)=(-1)**n d^n/dT^n G_0(R,T)
40 !>
41 !> The function is only computed for values of R,T which fulfil
42 !>
43 !> R**2 - 11.0_dp*R + 0.0_dp < T < R**2 + 11.0_dp*R + 50.0_dp where R>=0 T>=0
44 !>
45 !> for T larger than the upper bound, 0 is returned
46 !> (which is accurate at least up to 1.0E-16)
47 !> while for T smaller than the lower bound, the caller is instructed
48 !> to use the conventional gamma function instead
49 !> (i.e. the limit of above expression for R to Infinity)
50 !>
51 !> \author Joost VandeVondele and Manuel Guidon
52 !> \par History
53 !> Nov 2008, 2009 Joost VandeVondele and Manuel Guidon
54 !> May 2019 A. Bussy: Added a get_maxl_init function to get current status of nderiv_init and
55 !> moved the file to common (made it accessible from aobasis, same place as gamma.F).
56 ! **************************************************************************************************
57 MODULE t_c_g0
58  USE kinds, ONLY: dp
59  USE message_passing, ONLY: mp_comm_type
60 #include "../base/base_uses.f90"
61 
62  IMPLICIT NONE
63  PRIVATE
64 
65  PUBLIC :: t_c_g0_n, init, free_c0, get_lmax_init
66  REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE :: c0
67  INTEGER, PARAMETER :: degree = 13
68  REAL(KIND=dp), PARAMETER :: target_error = 0.100000e-08
69  INTEGER, PARAMETER :: nderiv_max = 21
70  INTEGER, SAVE :: nderiv_init = -1
71  INTEGER, SAVE :: patches = -1
72 
73 CONTAINS
74 
75 ! **************************************************************************************************
76 !> \brief ...
77 !> \param RES ...
78 !> \param use_gamma ...
79 !> \param R ...
80 !> \param T ...
81 !> \param NDERIV ...
82 ! **************************************************************************************************
83  SUBROUTINE t_c_g0_n(RES, use_gamma, R, T, NDERIV)
84  REAL(kind=dp), INTENT(OUT) :: res(*)
85  LOGICAL, INTENT(OUT) :: use_gamma
86  REAL(kind=dp), INTENT(IN) :: r, t
87  INTEGER, INTENT(IN) :: nderiv
88 
89  REAL(kind=dp) :: lower, tg1, tg2, upper, x1, x2
90 
91  use_gamma = .false.
92  upper = r**2 + 11.0_dp*r + 50.0_dp
93  lower = r**2 - 11.0_dp*r + 0.0_dp
94  IF (t > upper) THEN
95  res(1:nderiv + 1) = 0.0_dp
96  RETURN
97  END IF
98  IF (r <= 11.0_dp) THEN
99  x2 = r/11.0_dp
100  upper = r**2 + 11.0_dp*r + 50.0_dp
101  lower = 0.0_dp
102  x1 = (t - lower)/(upper - lower)
103  IF (x1 <= 0.500000000000000000e+00_dp) THEN
104  IF (x2 <= 0.500000000000000000e+00_dp) THEN
105  IF (x2 <= 0.250000000000000000e+00_dp) THEN
106  IF (x2 <= 0.125000000000000000e+00_dp) THEN
107  IF (x1 <= 0.250000000000000000e+00_dp) THEN
108  IF (x2 <= 0.625000000000000000e-01_dp) THEN
109  IF (x1 <= 0.125000000000000000e+00_dp) THEN
110  IF (x2 <= 0.312500000000000000e-01_dp) THEN
111  IF (x1 <= 0.625000000000000000e-01_dp) THEN
112  IF (x2 <= 0.156250000000000000e-01_dp) THEN
113  IF (x1 <= 0.312500000000000000e-01_dp) THEN
114  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
115  tg2 = (2*x2 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
116  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 1))
117  ELSE
118  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
119  tg2 = (2*x2 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
120  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 2))
121  END IF
122  ELSE
123  IF (x1 <= 0.312500000000000000e-01_dp) THEN
124  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
125  tg2 = (2*x2 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
126  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 3))
127  ELSE
128  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
129  tg2 = (2*x2 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
130  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 4))
131  END IF
132  END IF
133  ELSE
134  IF (x2 <= 0.156250000000000000e-01_dp) THEN
135  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
136  tg2 = (2*x2 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
137  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 5))
138  ELSE
139  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
140  tg2 = (2*x2 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
141  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 6))
142  END IF
143  END IF
144  ELSE
145  IF (x1 <= 0.625000000000000000e-01_dp) THEN
146  IF (x2 <= 0.468750000000000000e-01_dp) THEN
147  IF (x1 <= 0.312500000000000000e-01_dp) THEN
148  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
149  tg2 = (2*x2 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
150  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 7))
151  ELSE
152  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
153  tg2 = (2*x2 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
154  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 8))
155  END IF
156  ELSE
157  IF (x1 <= 0.312500000000000000e-01_dp) THEN
158  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
159  tg2 = (2*x2 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
160  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 9))
161  ELSE
162  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
163  tg2 = (2*x2 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
164  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 10))
165  END IF
166  END IF
167  ELSE
168  IF (x2 <= 0.468750000000000000e-01_dp) THEN
169  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
170  tg2 = (2*x2 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
171  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 11))
172  ELSE
173  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
174  tg2 = (2*x2 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
175  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 12))
176  END IF
177  END IF
178  END IF
179  ELSE
180  IF (x2 <= 0.312500000000000000e-01_dp) THEN
181  IF (x1 <= 0.187500000000000000e+00_dp) THEN
182  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
183  tg2 = (2*x2 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
184  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 13))
185  ELSE
186  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
187  tg2 = (2*x2 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
188  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 14))
189  END IF
190  ELSE
191  IF (x1 <= 0.187500000000000000e+00_dp) THEN
192  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
193  tg2 = (2*x2 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
194  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 15))
195  ELSE
196  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
197  tg2 = (2*x2 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
198  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 16))
199  END IF
200  END IF
201  END IF
202  ELSE
203  IF (x1 <= 0.125000000000000000e+00_dp) THEN
204  IF (x2 <= 0.937500000000000000e-01_dp) THEN
205  IF (x1 <= 0.625000000000000000e-01_dp) THEN
206  IF (x2 <= 0.781250000000000000e-01_dp) THEN
207  IF (x1 <= 0.312500000000000000e-01_dp) THEN
208  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
209  tg2 = (2*x2 - 0.140625000000000000e+00_dp)*0.640000000000000000e+02_dp
210  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 17))
211  ELSE
212  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
213  tg2 = (2*x2 - 0.140625000000000000e+00_dp)*0.640000000000000000e+02_dp
214  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 18))
215  END IF
216  ELSE
217  IF (x1 <= 0.312500000000000000e-01_dp) THEN
218  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
219  tg2 = (2*x2 - 0.171875000000000000e+00_dp)*0.640000000000000000e+02_dp
220  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 19))
221  ELSE
222  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
223  tg2 = (2*x2 - 0.171875000000000000e+00_dp)*0.640000000000000000e+02_dp
224  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 20))
225  END IF
226  END IF
227  ELSE
228  IF (x2 <= 0.781250000000000000e-01_dp) THEN
229  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
230  tg2 = (2*x2 - 0.140625000000000000e+00_dp)*0.640000000000000000e+02_dp
231  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 21))
232  ELSE
233  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
234  tg2 = (2*x2 - 0.171875000000000000e+00_dp)*0.640000000000000000e+02_dp
235  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 22))
236  END IF
237  END IF
238  ELSE
239  IF (x1 <= 0.625000000000000000e-01_dp) THEN
240  IF (x2 <= 0.109375000000000000e+00_dp) THEN
241  IF (x1 <= 0.312500000000000000e-01_dp) THEN
242  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
243  tg2 = (2*x2 - 0.203125000000000000e+00_dp)*0.640000000000000000e+02_dp
244  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 23))
245  ELSE
246  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
247  tg2 = (2*x2 - 0.203125000000000000e+00_dp)*0.640000000000000000e+02_dp
248  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 24))
249  END IF
250  ELSE
251  IF (x1 <= 0.312500000000000000e-01_dp) THEN
252  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
253  tg2 = (2*x2 - 0.234375000000000000e+00_dp)*0.640000000000000000e+02_dp
254  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 25))
255  ELSE
256  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
257  tg2 = (2*x2 - 0.234375000000000000e+00_dp)*0.640000000000000000e+02_dp
258  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 26))
259  END IF
260  END IF
261  ELSE
262  IF (x2 <= 0.109375000000000000e+00_dp) THEN
263  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
264  tg2 = (2*x2 - 0.203125000000000000e+00_dp)*0.640000000000000000e+02_dp
265  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 27))
266  ELSE
267  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
268  tg2 = (2*x2 - 0.234375000000000000e+00_dp)*0.640000000000000000e+02_dp
269  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 28))
270  END IF
271  END IF
272  END IF
273  ELSE
274  IF (x1 <= 0.187500000000000000e+00_dp) THEN
275  IF (x2 <= 0.937500000000000000e-01_dp) THEN
276  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
277  tg2 = (2*x2 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
278  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 29))
279  ELSE
280  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
281  tg2 = (2*x2 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
282  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 30))
283  END IF
284  ELSE
285  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
286  tg2 = (2*x2 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
287  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 31))
288  END IF
289  END IF
290  END IF
291  ELSE
292  IF (x1 <= 0.375000000000000000e+00_dp) THEN
293  tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
294  tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
295  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 32))
296  ELSE
297  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
298  tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
299  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 33))
300  END IF
301  END IF
302  ELSE
303  IF (x1 <= 0.250000000000000000e+00_dp) THEN
304  IF (x2 <= 0.187500000000000000e+00_dp) THEN
305  IF (x1 <= 0.125000000000000000e+00_dp) THEN
306  IF (x2 <= 0.156250000000000000e+00_dp) THEN
307  IF (x1 <= 0.625000000000000000e-01_dp) THEN
308  IF (x1 <= 0.312500000000000000e-01_dp) THEN
309  IF (x2 <= 0.140625000000000000e+00_dp) THEN
310  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
311  tg2 = (2*x2 - 0.265625000000000000e+00_dp)*0.640000000000000000e+02_dp
312  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 34))
313  ELSE
314  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
315  tg2 = (2*x2 - 0.296875000000000000e+00_dp)*0.640000000000000000e+02_dp
316  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 35))
317  END IF
318  ELSE
319  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
320  tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
321  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 36))
322  END IF
323  ELSE
324  IF (x1 <= 0.937500000000000000e-01_dp) THEN
325  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
326  tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
327  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 37))
328  ELSE
329  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
330  tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
331  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 38))
332  END IF
333  END IF
334  ELSE
335  IF (x1 <= 0.625000000000000000e-01_dp) THEN
336  IF (x1 <= 0.312500000000000000e-01_dp) THEN
337  IF (x2 <= 0.171875000000000000e+00_dp) THEN
338  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
339  tg2 = (2*x2 - 0.328125000000000000e+00_dp)*0.640000000000000000e+02_dp
340  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 39))
341  ELSE
342  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
343  tg2 = (2*x2 - 0.359375000000000000e+00_dp)*0.640000000000000000e+02_dp
344  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 40))
345  END IF
346  ELSE
347  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
348  tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
349  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 41))
350  END IF
351  ELSE
352  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
353  tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
354  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 42))
355  END IF
356  END IF
357  ELSE
358  IF (x1 <= 0.187500000000000000e+00_dp) THEN
359  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
360  tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
361  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 43))
362  ELSE
363  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
364  tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
365  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 44))
366  END IF
367  END IF
368  ELSE
369  IF (x1 <= 0.125000000000000000e+00_dp) THEN
370  IF (x1 <= 0.625000000000000000e-01_dp) THEN
371  IF (x2 <= 0.218750000000000000e+00_dp) THEN
372  IF (x1 <= 0.312500000000000000e-01_dp) THEN
373  IF (x2 <= 0.203125000000000000e+00_dp) THEN
374  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
375  tg2 = (2*x2 - 0.390625000000000000e+00_dp)*0.640000000000000000e+02_dp
376  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 45))
377  ELSE
378  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
379  tg2 = (2*x2 - 0.421875000000000000e+00_dp)*0.640000000000000000e+02_dp
380  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 46))
381  END IF
382  ELSE
383  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
384  tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
385  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 47))
386  END IF
387  ELSE
388  IF (x1 <= 0.312500000000000000e-01_dp) THEN
389  IF (x2 <= 0.234375000000000000e+00_dp) THEN
390  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
391  tg2 = (2*x2 - 0.453125000000000000e+00_dp)*0.640000000000000000e+02_dp
392  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 48))
393  ELSE
394  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
395  tg2 = (2*x2 - 0.484375000000000000e+00_dp)*0.640000000000000000e+02_dp
396  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 49))
397  END IF
398  ELSE
399  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
400  tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
401  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 50))
402  END IF
403  END IF
404  ELSE
405  IF (x2 <= 0.218750000000000000e+00_dp) THEN
406  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
407  tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
408  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 51))
409  ELSE
410  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
411  tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
412  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 52))
413  END IF
414  END IF
415  ELSE
416  IF (x1 <= 0.187500000000000000e+00_dp) THEN
417  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
418  tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
419  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 53))
420  ELSE
421  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
422  tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
423  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 54))
424  END IF
425  END IF
426  END IF
427  ELSE
428  IF (x1 <= 0.375000000000000000e+00_dp) THEN
429  IF (x1 <= 0.312500000000000000e+00_dp) THEN
430  tg1 = (2*x1 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
431  tg2 = (2*x2 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
432  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 55))
433  ELSE
434  tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
435  tg2 = (2*x2 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
436  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 56))
437  END IF
438  ELSE
439  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
440  tg2 = (2*x2 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
441  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 57))
442  END IF
443  END IF
444  END IF
445  ELSE
446  IF (x1 <= 0.250000000000000000e+00_dp) THEN
447  IF (x1 <= 0.125000000000000000e+00_dp) THEN
448  IF (x1 <= 0.625000000000000000e-01_dp) THEN
449  IF (x2 <= 0.375000000000000000e+00_dp) THEN
450  IF (x2 <= 0.312500000000000000e+00_dp) THEN
451  IF (x1 <= 0.312500000000000000e-01_dp) THEN
452  IF (x2 <= 0.281250000000000000e+00_dp) THEN
453  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
454  tg2 = (2*x2 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
455  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 58))
456  ELSE
457  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
458  tg2 = (2*x2 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
459  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 59))
460  END IF
461  ELSE
462  IF (x2 <= 0.281250000000000000e+00_dp) THEN
463  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
464  tg2 = (2*x2 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
465  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 60))
466  ELSE
467  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
468  tg2 = (2*x2 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
469  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 61))
470  END IF
471  END IF
472  ELSE
473  IF (x1 <= 0.312500000000000000e-01_dp) THEN
474  IF (x2 <= 0.343750000000000000e+00_dp) THEN
475  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
476  tg2 = (2*x2 - 0.656250000000000000e+00_dp)*0.320000000000000000e+02_dp
477  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 62))
478  ELSE
479  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
480  tg2 = (2*x2 - 0.718750000000000000e+00_dp)*0.320000000000000000e+02_dp
481  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 63))
482  END IF
483  ELSE
484  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
485  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
486  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 64))
487  END IF
488  END IF
489  ELSE
490  IF (x1 <= 0.312500000000000000e-01_dp) THEN
491  IF (x2 <= 0.437500000000000000e+00_dp) THEN
492  IF (x1 <= 0.156250000000000000e-01_dp) THEN
493  tg1 = (2*x1 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
494  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
495  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 65))
496  ELSE
497  tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
498  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
499  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 66))
500  END IF
501  ELSE
502  IF (x1 <= 0.156250000000000000e-01_dp) THEN
503  tg1 = (2*x1 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
504  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
505  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 67))
506  ELSE
507  tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
508  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
509  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 68))
510  END IF
511  END IF
512  ELSE
513  IF (x2 <= 0.437500000000000000e+00_dp) THEN
514  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
515  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
516  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 69))
517  ELSE
518  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
519  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
520  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 70))
521  END IF
522  END IF
523  END IF
524  ELSE
525  IF (x2 <= 0.375000000000000000e+00_dp) THEN
526  IF (x2 <= 0.312500000000000000e+00_dp) THEN
527  IF (x1 <= 0.937500000000000000e-01_dp) THEN
528  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
529  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
530  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 71))
531  ELSE
532  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
533  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
534  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 72))
535  END IF
536  ELSE
537  IF (x1 <= 0.937500000000000000e-01_dp) THEN
538  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
539  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
540  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 73))
541  ELSE
542  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
543  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
544  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 74))
545  END IF
546  END IF
547  ELSE
548  IF (x1 <= 0.937500000000000000e-01_dp) THEN
549  IF (x2 <= 0.437500000000000000e+00_dp) THEN
550  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
551  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
552  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 75))
553  ELSE
554  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
555  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
556  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 76))
557  END IF
558  ELSE
559  IF (x2 <= 0.437500000000000000e+00_dp) THEN
560  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
561  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
562  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 77))
563  ELSE
564  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
565  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
566  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 78))
567  END IF
568  END IF
569  END IF
570  END IF
571  ELSE
572  IF (x2 <= 0.375000000000000000e+00_dp) THEN
573  IF (x1 <= 0.187500000000000000e+00_dp) THEN
574  IF (x2 <= 0.312500000000000000e+00_dp) THEN
575  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
576  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
577  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 79))
578  ELSE
579  IF (x1 <= 0.156250000000000000e+00_dp) THEN
580  tg1 = (2*x1 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
581  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
582  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 80))
583  ELSE
584  tg1 = (2*x1 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
585  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
586  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 81))
587  END IF
588  END IF
589  ELSE
590  IF (x2 <= 0.312500000000000000e+00_dp) THEN
591  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
592  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
593  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 82))
594  ELSE
595  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
596  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
597  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 83))
598  END IF
599  END IF
600  ELSE
601  IF (x1 <= 0.187500000000000000e+00_dp) THEN
602  IF (x2 <= 0.437500000000000000e+00_dp) THEN
603  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
604  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
605  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 84))
606  ELSE
607  IF (x1 <= 0.156250000000000000e+00_dp) THEN
608  tg1 = (2*x1 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
609  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
610  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 85))
611  ELSE
612  tg1 = (2*x1 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
613  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
614  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 86))
615  END IF
616  END IF
617  ELSE
618  IF (x1 <= 0.218750000000000000e+00_dp) THEN
619  tg1 = (2*x1 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
620  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
621  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 87))
622  ELSE
623  tg1 = (2*x1 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
624  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
625  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 88))
626  END IF
627  END IF
628  END IF
629  END IF
630  ELSE
631  IF (x1 <= 0.375000000000000000e+00_dp) THEN
632  IF (x2 <= 0.375000000000000000e+00_dp) THEN
633  IF (x1 <= 0.312500000000000000e+00_dp) THEN
634  tg1 = (2*x1 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
635  tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
636  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 89))
637  ELSE
638  tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
639  tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
640  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 90))
641  END IF
642  ELSE
643  IF (x1 <= 0.312500000000000000e+00_dp) THEN
644  IF (x1 <= 0.281250000000000000e+00_dp) THEN
645  tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
646  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
647  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 91))
648  ELSE
649  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
650  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
651  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 92))
652  END IF
653  ELSE
654  tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
655  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
656  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 93))
657  END IF
658  END IF
659  ELSE
660  IF (x1 <= 0.437500000000000000e+00_dp) THEN
661  tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
662  tg2 = (2*x2 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
663  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 94))
664  ELSE
665  tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
666  tg2 = (2*x2 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
667  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 95))
668  END IF
669  END IF
670  END IF
671  END IF
672  ELSE
673  IF (x1 <= 0.250000000000000000e+00_dp) THEN
674  IF (x1 <= 0.125000000000000000e+00_dp) THEN
675  IF (x1 <= 0.625000000000000000e-01_dp) THEN
676  IF (x1 <= 0.312500000000000000e-01_dp) THEN
677  IF (x1 <= 0.156250000000000000e-01_dp) THEN
678  IF (x1 <= 0.781250000000000000e-02_dp) THEN
679  IF (x2 <= 0.750000000000000000e+00_dp) THEN
680  tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
681  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
682  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 96))
683  ELSE
684  tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
685  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
686  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 97))
687  END IF
688  ELSE
689  IF (x2 <= 0.750000000000000000e+00_dp) THEN
690  tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
691  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
692  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 98))
693  ELSE
694  tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
695  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
696  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 99))
697  END IF
698  END IF
699  ELSE
700  IF (x2 <= 0.750000000000000000e+00_dp) THEN
701  tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
702  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
703  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 100))
704  ELSE
705  tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
706  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
707  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 101))
708  END IF
709  END IF
710  ELSE
711  IF (x2 <= 0.750000000000000000e+00_dp) THEN
712  IF (x2 <= 0.625000000000000000e+00_dp) THEN
713  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
714  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
715  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 102))
716  ELSE
717  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
718  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
719  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 103))
720  END IF
721  ELSE
722  IF (x1 <= 0.468750000000000000e-01_dp) THEN
723  tg1 = (2*x1 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
724  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
725  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 104))
726  ELSE
727  tg1 = (2*x1 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
728  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
729  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 105))
730  END IF
731  END IF
732  END IF
733  ELSE
734  IF (x2 <= 0.750000000000000000e+00_dp) THEN
735  IF (x2 <= 0.625000000000000000e+00_dp) THEN
736  IF (x1 <= 0.937500000000000000e-01_dp) THEN
737  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
738  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
739  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 106))
740  ELSE
741  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
742  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
743  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 107))
744  END IF
745  ELSE
746  IF (x1 <= 0.937500000000000000e-01_dp) THEN
747  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
748  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
749  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 108))
750  ELSE
751  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
752  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
753  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 109))
754  END IF
755  END IF
756  ELSE
757  IF (x1 <= 0.937500000000000000e-01_dp) THEN
758  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
759  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
760  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 110))
761  ELSE
762  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
763  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
764  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 111))
765  END IF
766  END IF
767  END IF
768  ELSE
769  IF (x2 <= 0.750000000000000000e+00_dp) THEN
770  IF (x2 <= 0.625000000000000000e+00_dp) THEN
771  IF (x1 <= 0.187500000000000000e+00_dp) THEN
772  IF (x1 <= 0.156250000000000000e+00_dp) THEN
773  tg1 = (2*x1 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
774  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
775  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 112))
776  ELSE
777  tg1 = (2*x1 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
778  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
779  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 113))
780  END IF
781  ELSE
782  IF (x1 <= 0.218750000000000000e+00_dp) THEN
783  tg1 = (2*x1 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
784  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
785  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 114))
786  ELSE
787  tg1 = (2*x1 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
788  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
789  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 115))
790  END IF
791  END IF
792  ELSE
793  IF (x1 <= 0.187500000000000000e+00_dp) THEN
794  IF (x1 <= 0.156250000000000000e+00_dp) THEN
795  tg1 = (2*x1 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
796  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
797  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 116))
798  ELSE
799  tg1 = (2*x1 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
800  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
801  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 117))
802  END IF
803  ELSE
804  IF (x1 <= 0.218750000000000000e+00_dp) THEN
805  tg1 = (2*x1 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
806  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
807  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 118))
808  ELSE
809  tg1 = (2*x1 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
810  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
811  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 119))
812  END IF
813  END IF
814  END IF
815  ELSE
816  IF (x1 <= 0.187500000000000000e+00_dp) THEN
817  IF (x2 <= 0.875000000000000000e+00_dp) THEN
818  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
819  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
820  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 120))
821  ELSE
822  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
823  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
824  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 121))
825  END IF
826  ELSE
827  IF (x2 <= 0.875000000000000000e+00_dp) THEN
828  IF (x1 <= 0.218750000000000000e+00_dp) THEN
829  tg1 = (2*x1 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
830  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
831  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 122))
832  ELSE
833  tg1 = (2*x1 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
834  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
835  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 123))
836  END IF
837  ELSE
838  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
839  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
840  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 124))
841  END IF
842  END IF
843  END IF
844  END IF
845  ELSE
846  IF (x1 <= 0.375000000000000000e+00_dp) THEN
847  IF (x2 <= 0.750000000000000000e+00_dp) THEN
848  IF (x1 <= 0.312500000000000000e+00_dp) THEN
849  IF (x2 <= 0.625000000000000000e+00_dp) THEN
850  IF (x1 <= 0.281250000000000000e+00_dp) THEN
851  tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
852  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
853  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 125))
854  ELSE
855  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
856  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
857  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 126))
858  END IF
859  ELSE
860  IF (x1 <= 0.281250000000000000e+00_dp) THEN
861  tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
862  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
863  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 127))
864  ELSE
865  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
866  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
867  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 128))
868  END IF
869  END IF
870  ELSE
871  IF (x2 <= 0.625000000000000000e+00_dp) THEN
872  tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
873  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
874  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 129))
875  ELSE
876  IF (x1 <= 0.343750000000000000e+00_dp) THEN
877  tg1 = (2*x1 - 0.656250000000000000e+00_dp)*0.320000000000000000e+02_dp
878  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
879  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 130))
880  ELSE
881  tg1 = (2*x1 - 0.718750000000000000e+00_dp)*0.320000000000000000e+02_dp
882  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
883  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 131))
884  END IF
885  END IF
886  END IF
887  ELSE
888  IF (x1 <= 0.312500000000000000e+00_dp) THEN
889  IF (x2 <= 0.875000000000000000e+00_dp) THEN
890  IF (x1 <= 0.281250000000000000e+00_dp) THEN
891  tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
892  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
893  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 132))
894  ELSE
895  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
896  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
897  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 133))
898  END IF
899  ELSE
900  IF (x1 <= 0.281250000000000000e+00_dp) THEN
901  tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
902  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
903  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 134))
904  ELSE
905  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
906  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
907  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 135))
908  END IF
909  END IF
910  ELSE
911  IF (x2 <= 0.875000000000000000e+00_dp) THEN
912  IF (x1 <= 0.343750000000000000e+00_dp) THEN
913  tg1 = (2*x1 - 0.656250000000000000e+00_dp)*0.320000000000000000e+02_dp
914  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
915  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 136))
916  ELSE
917  tg1 = (2*x1 - 0.718750000000000000e+00_dp)*0.320000000000000000e+02_dp
918  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
919  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 137))
920  END IF
921  ELSE
922  tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
923  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
924  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 138))
925  END IF
926  END IF
927  END IF
928  ELSE
929  IF (x2 <= 0.750000000000000000e+00_dp) THEN
930  IF (x1 <= 0.437500000000000000e+00_dp) THEN
931  IF (x2 <= 0.625000000000000000e+00_dp) THEN
932  tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
933  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
934  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 139))
935  ELSE
936  tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
937  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
938  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 140))
939  END IF
940  ELSE
941  tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
942  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
943  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 141))
944  END IF
945  ELSE
946  IF (x1 <= 0.437500000000000000e+00_dp) THEN
947  IF (x2 <= 0.875000000000000000e+00_dp) THEN
948  tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
949  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
950  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 142))
951  ELSE
952  tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
953  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
954  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 143))
955  END IF
956  ELSE
957  IF (x2 <= 0.875000000000000000e+00_dp) THEN
958  tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
959  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
960  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 144))
961  ELSE
962  tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
963  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
964  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 145))
965  END IF
966  END IF
967  END IF
968  END IF
969  END IF
970  END IF
971  ELSE
972  IF (x1 <= 0.750000000000000000e+00_dp) THEN
973  IF (x2 <= 0.500000000000000000e+00_dp) THEN
974  IF (x1 <= 0.625000000000000000e+00_dp) THEN
975  IF (x2 <= 0.250000000000000000e+00_dp) THEN
976  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
977  tg2 = (2*x2 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
978  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 146))
979  ELSE
980  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
981  tg2 = (2*x2 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
982  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 147))
983  END IF
984  ELSE
985  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
986  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
987  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 148))
988  END IF
989  ELSE
990  IF (x1 <= 0.625000000000000000e+00_dp) THEN
991  IF (x2 <= 0.750000000000000000e+00_dp) THEN
992  IF (x1 <= 0.562500000000000000e+00_dp) THEN
993  tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
994  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
995  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 149))
996  ELSE
997  tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
998  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
999  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 150))
1000  END IF
1001  ELSE
1002  IF (x1 <= 0.562500000000000000e+00_dp) THEN
1003  tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1004  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1005  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 151))
1006  ELSE
1007  tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1008  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1009  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 152))
1010  END IF
1011  END IF
1012  ELSE
1013  IF (x1 <= 0.687500000000000000e+00_dp) THEN
1014  tg1 = (2*x1 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
1015  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1016  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 153))
1017  ELSE
1018  tg1 = (2*x1 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
1019  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1020  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 154))
1021  END IF
1022  END IF
1023  END IF
1024  ELSE
1025  tg1 = (2*x1 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1026  tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1027  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 155))
1028  END IF
1029  END IF
1030  ELSE
1031  IF (t < lower) THEN
1032  use_gamma = .true.
1033  RETURN
1034  END IF
1035  x2 = 11.0_dp/r
1036  x1 = (t - lower)/(upper - lower)
1037  IF (x1 <= 0.500000000000000000e+00_dp) THEN
1038  IF (x1 <= 0.250000000000000000e+00_dp) THEN
1039  IF (x2 <= 0.500000000000000000e+00_dp) THEN
1040  IF (x1 <= 0.125000000000000000e+00_dp) THEN
1041  IF (x2 <= 0.250000000000000000e+00_dp) THEN
1042  tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
1043  tg2 = (2*x2 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
1044  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 156))
1045  ELSE
1046  tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
1047  tg2 = (2*x2 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
1048  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 157))
1049  END IF
1050  ELSE
1051  tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
1052  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1053  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 158))
1054  END IF
1055  ELSE
1056  IF (x1 <= 0.125000000000000000e+00_dp) THEN
1057  IF (x2 <= 0.750000000000000000e+00_dp) THEN
1058  IF (x2 <= 0.625000000000000000e+00_dp) THEN
1059  tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
1060  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1061  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 159))
1062  ELSE
1063  IF (x1 <= 0.625000000000000000e-01_dp) THEN
1064  tg1 = (2*x1 - 0.625000000000000000e-01_dp)*0.160000000000000000e+02_dp
1065  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1066  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 160))
1067  ELSE
1068  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
1069  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1070  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 161))
1071  END IF
1072  END IF
1073  ELSE
1074  IF (x1 <= 0.625000000000000000e-01_dp) THEN
1075  IF (x2 <= 0.875000000000000000e+00_dp) THEN
1076  IF (x1 <= 0.312500000000000000e-01_dp) THEN
1077  IF (x2 <= 0.812500000000000000e+00_dp) THEN
1078  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
1079  tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
1080  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 162))
1081  ELSE
1082  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
1083  tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
1084  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 163))
1085  END IF
1086  ELSE
1087  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
1088  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1089  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 164))
1090  END IF
1091  ELSE
1092  IF (x1 <= 0.312500000000000000e-01_dp) THEN
1093  IF (x2 <= 0.937500000000000000e+00_dp) THEN
1094  IF (x1 <= 0.156250000000000000e-01_dp) THEN
1095  IF (x2 <= 0.906250000000000000e+00_dp) THEN
1096  tg1 = (2*x1 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
1097  tg2 = (2*x2 - 0.178125000000000000e+01_dp)*0.320000000000000000e+02_dp
1098  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 165))
1099  ELSE
1100  tg1 = (2*x1 - 0.156250000000000000e-01_dp)*0.640000000000000000e+02_dp
1101  tg2 = (2*x2 - 0.184375000000000000e+01_dp)*0.320000000000000000e+02_dp
1102  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 166))
1103  END IF
1104  ELSE
1105  tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
1106  tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
1107  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 167))
1108  END IF
1109  ELSE
1110  IF (x1 <= 0.156250000000000000e-01_dp) THEN
1111  IF (x2 <= 0.968750000000000000e+00_dp) THEN
1112  IF (x1 <= 0.781250000000000000e-02_dp) THEN
1113  IF (x2 <= 0.953125000000000000e+00_dp) THEN
1114  tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
1115  tg2 = (2*x2 - 0.189062500000000000e+01_dp)*0.640000000000000000e+02_dp
1116  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 168))
1117  ELSE
1118  tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
1119  tg2 = (2*x2 - 0.192187500000000000e+01_dp)*0.640000000000000000e+02_dp
1120  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 169))
1121  END IF
1122  ELSE
1123  tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
1124  tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
1125  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 170))
1126  END IF
1127  ELSE
1128  IF (x1 <= 0.781250000000000000e-02_dp) THEN
1129  IF (x2 <= 0.984375000000000000e+00_dp) THEN
1130  tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
1131  tg2 = (2*x2 - 0.195312500000000000e+01_dp)*0.640000000000000000e+02_dp
1132  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 171))
1133  ELSE
1134  tg1 = (2*x1 - 0.781250000000000000e-02_dp)*0.128000000000000000e+03_dp
1135  tg2 = (2*x2 - 0.198437500000000000e+01_dp)*0.640000000000000000e+02_dp
1136  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 172))
1137  END IF
1138  ELSE
1139  IF (x2 <= 0.984375000000000000e+00_dp) THEN
1140  tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
1141  tg2 = (2*x2 - 0.195312500000000000e+01_dp)*0.640000000000000000e+02_dp
1142  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 173))
1143  ELSE
1144  tg1 = (2*x1 - 0.234375000000000000e-01_dp)*0.128000000000000000e+03_dp
1145  tg2 = (2*x2 - 0.198437500000000000e+01_dp)*0.640000000000000000e+02_dp
1146  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 174))
1147  END IF
1148  END IF
1149  END IF
1150  ELSE
1151  IF (x2 <= 0.968750000000000000e+00_dp) THEN
1152  tg1 = (2*x1 - 0.468750000000000000e-01_dp)*0.640000000000000000e+02_dp
1153  tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
1154  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 175))
1155  ELSE
1156  IF (x1 <= 0.234375000000000000e-01_dp) THEN
1157  tg1 = (2*x1 - 0.390625000000000000e-01_dp)*0.128000000000000000e+03_dp
1158  tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
1159  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 176))
1160  ELSE
1161  tg1 = (2*x1 - 0.546875000000000000e-01_dp)*0.128000000000000000e+03_dp
1162  tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
1163  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 177))
1164  END IF
1165  END IF
1166  END IF
1167  END IF
1168  ELSE
1169  IF (x2 <= 0.937500000000000000e+00_dp) THEN
1170  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
1171  tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
1172  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 178))
1173  ELSE
1174  IF (x1 <= 0.468750000000000000e-01_dp) THEN
1175  IF (x2 <= 0.968750000000000000e+00_dp) THEN
1176  tg1 = (2*x1 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
1177  tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
1178  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 179))
1179  ELSE
1180  tg1 = (2*x1 - 0.781250000000000000e-01_dp)*0.640000000000000000e+02_dp
1181  tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
1182  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 180))
1183  END IF
1184  ELSE
1185  tg1 = (2*x1 - 0.109375000000000000e+00_dp)*0.640000000000000000e+02_dp
1186  tg2 = (2*x2 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
1187  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 181))
1188  END IF
1189  END IF
1190  END IF
1191  END IF
1192  ELSE
1193  IF (x2 <= 0.875000000000000000e+00_dp) THEN
1194  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
1195  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1196  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 182))
1197  ELSE
1198  IF (x1 <= 0.937500000000000000e-01_dp) THEN
1199  IF (x2 <= 0.937500000000000000e+00_dp) THEN
1200  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
1201  tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
1202  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 183))
1203  ELSE
1204  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
1205  tg2 = (2*x2 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
1206  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 184))
1207  END IF
1208  ELSE
1209  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
1210  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1211  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 185))
1212  END IF
1213  END IF
1214  END IF
1215  END IF
1216  ELSE
1217  IF (x2 <= 0.750000000000000000e+00_dp) THEN
1218  tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
1219  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
1220  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 186))
1221  ELSE
1222  IF (x1 <= 0.187500000000000000e+00_dp) THEN
1223  IF (x2 <= 0.875000000000000000e+00_dp) THEN
1224  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
1225  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1226  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 187))
1227  ELSE
1228  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
1229  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1230  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 188))
1231  END IF
1232  ELSE
1233  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
1234  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1235  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 189))
1236  END IF
1237  END IF
1238  END IF
1239  END IF
1240  ELSE
1241  IF (x1 <= 0.375000000000000000e+00_dp) THEN
1242  IF (x1 <= 0.312500000000000000e+00_dp) THEN
1243  IF (x1 <= 0.281250000000000000e+00_dp) THEN
1244  tg1 = (2*x1 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
1245  tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1246  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 190))
1247  ELSE
1248  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
1249  tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1250  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 191))
1251  END IF
1252  ELSE
1253  IF (x2 <= 0.500000000000000000e+00_dp) THEN
1254  tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
1255  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1256  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 192))
1257  ELSE
1258  tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
1259  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1260  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 193))
1261  END IF
1262  END IF
1263  ELSE
1264  IF (x1 <= 0.437500000000000000e+00_dp) THEN
1265  IF (x2 <= 0.500000000000000000e+00_dp) THEN
1266  tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
1267  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1268  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 194))
1269  ELSE
1270  IF (x1 <= 0.406250000000000000e+00_dp) THEN
1271  tg1 = (2*x1 - 0.781250000000000000e+00_dp)*0.320000000000000000e+02_dp
1272  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1273  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 195))
1274  ELSE
1275  tg1 = (2*x1 - 0.843750000000000000e+00_dp)*0.320000000000000000e+02_dp
1276  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1277  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 196))
1278  END IF
1279  END IF
1280  ELSE
1281  IF (x2 <= 0.500000000000000000e+00_dp) THEN
1282  tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
1283  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1284  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 197))
1285  ELSE
1286  tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
1287  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1288  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 198))
1289  END IF
1290  END IF
1291  END IF
1292  END IF
1293  ELSE
1294  IF (x1 <= 0.750000000000000000e+00_dp) THEN
1295  IF (x1 <= 0.625000000000000000e+00_dp) THEN
1296  IF (x2 <= 0.500000000000000000e+00_dp) THEN
1297  IF (x1 <= 0.562500000000000000e+00_dp) THEN
1298  tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1299  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1300  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 199))
1301  ELSE
1302  tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1303  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1304  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 200))
1305  END IF
1306  ELSE
1307  IF (x1 <= 0.562500000000000000e+00_dp) THEN
1308  tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1309  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1310  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 201))
1311  ELSE
1312  tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1313  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1314  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 202))
1315  END IF
1316  END IF
1317  ELSE
1318  IF (x2 <= 0.500000000000000000e+00_dp) THEN
1319  IF (x1 <= 0.687500000000000000e+00_dp) THEN
1320  tg1 = (2*x1 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
1321  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1322  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 203))
1323  ELSE
1324  tg1 = (2*x1 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
1325  tg2 = (2*x2 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
1326  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 204))
1327  END IF
1328  ELSE
1329  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1330  tg2 = (2*x2 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
1331  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 205))
1332  END IF
1333  END IF
1334  ELSE
1335  IF (x1 <= 0.875000000000000000e+00_dp) THEN
1336  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1337  tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1338  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 206))
1339  ELSE
1340  tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1341  tg2 = (2*x2 - 0.100000000000000000e+01_dp)*0.100000000000000000e+01_dp
1342  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 207))
1343  END IF
1344  END IF
1345  END IF
1346  END IF
1347  END SUBROUTINE t_c_g0_n
1348 
1349 ! **************************************************************************************************
1350 !> \brief ...
1351 !> \param Nder the number of derivatives that will actually be used
1352 !> \param iunit contains the data file to initialize the table
1353 !> \param mepos ...
1354 !> \param group ...
1355 ! **************************************************************************************************
1356  SUBROUTINE init(Nder, iunit, mepos, group)
1357  INTEGER, INTENT(IN) :: nder, iunit, mepos
1358 
1359  CLASS(mp_comm_type), INTENT(IN) :: group
1360 
1361  INTEGER :: i
1362  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1363 
1364  patches = 207
1365  IF (nder > nderiv_max) cpabort("T_C_G0 init failed")
1366  nderiv_init = nder
1367  IF (ALLOCATED(c0)) DEALLOCATE (c0)
1368  ! round up to a multiple of 32 to give some generous alignment for each C0
1369  ALLOCATE (c0(32*((31 + (nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1370  ! init mpi'ed buffers to silence warnings under valgrind
1371  c0 = 1.0e99_dp
1372  IF (mepos == 0) THEN
1373  ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1374  DO i = 1, patches
1375  READ (iunit, *) chunk
1376  c0(1:(nder + 1)*(degree + 1)*(degree + 2)/2, i) = chunk(1:(nder + 1)*(degree + 1)*(degree + 2)/2)
1377  END DO
1378  DEALLOCATE (chunk)
1379  END IF
1380  CALL group%bcast(c0, 0)
1381 
1382  END SUBROUTINE init
1383 
1384 ! **************************************************************************************************
1385 !> \brief ...
1386 ! **************************************************************************************************
1387  SUBROUTINE free_c0()
1388  IF (ALLOCATED(c0)) DEALLOCATE (c0)
1389  nderiv_init = -1
1390  END SUBROUTINE free_c0
1391 
1392 ! **************************************************************************************************
1393 !> \brief ...
1394 !> \param RES ...
1395 !> \param NDERIV ...
1396 !> \param TG1 ...
1397 !> \param TG2 ...
1398 !> \param C0 ...
1399 ! **************************************************************************************************
1400  SUBROUTINE pd2val(RES, NDERIV, TG1, TG2, C0)
1401  REAL(kind=dp), INTENT(OUT) :: res(*)
1402  INTEGER, INTENT(IN) :: nderiv
1403  REAL(kind=dp), INTENT(IN) :: tg1, tg2, c0(105, *)
1404 
1405  REAL(kind=dp), PARAMETER :: sqrt2 = 1.4142135623730950488016887242096980785696718753_dp
1406 
1407  INTEGER :: k
1408  REAL(kind=dp) :: t1(0:13), t2(0:13)
1409 
1410  t1(0) = 1.0_dp
1411  t2(0) = 1.0_dp
1412  t1(1) = sqrt2*tg1
1413  t2(1) = sqrt2*tg2
1414  t1(2) = 2*tg1*t1(1) - sqrt2
1415  t2(2) = 2*tg2*t2(1) - sqrt2
1416  t1(3) = 2*tg1*t1(2) - t1(1)
1417  t2(3) = 2*tg2*t2(2) - t2(1)
1418  t1(4) = 2*tg1*t1(3) - t1(2)
1419  t2(4) = 2*tg2*t2(3) - t2(2)
1420  t1(5) = 2*tg1*t1(4) - t1(3)
1421  t2(5) = 2*tg2*t2(4) - t2(3)
1422  t1(6) = 2*tg1*t1(5) - t1(4)
1423  t2(6) = 2*tg2*t2(5) - t2(4)
1424  t1(7) = 2*tg1*t1(6) - t1(5)
1425  t2(7) = 2*tg2*t2(6) - t2(5)
1426  t1(8) = 2*tg1*t1(7) - t1(6)
1427  t2(8) = 2*tg2*t2(7) - t2(6)
1428  t1(9) = 2*tg1*t1(8) - t1(7)
1429  t2(9) = 2*tg2*t2(8) - t2(7)
1430  t1(10) = 2*tg1*t1(9) - t1(8)
1431  t2(10) = 2*tg2*t2(9) - t2(8)
1432  t1(11) = 2*tg1*t1(10) - t1(9)
1433  t2(11) = 2*tg2*t2(10) - t2(9)
1434  t1(12) = 2*tg1*t1(11) - t1(10)
1435  t2(12) = 2*tg2*t2(11) - t2(10)
1436  t1(13) = 2*tg1*t1(12) - t1(11)
1437  t2(13) = 2*tg2*t2(12) - t2(11)
1438  DO k = 1, nderiv + 1
1439  res(k) = 0.0_dp
1440  res(k) = res(k) + dot_product(t1(0:13), c0(1:14, k))*t2(0)
1441  res(k) = res(k) + dot_product(t1(0:12), c0(15:27, k))*t2(1)
1442  res(k) = res(k) + dot_product(t1(0:11), c0(28:39, k))*t2(2)
1443  res(k) = res(k) + dot_product(t1(0:10), c0(40:50, k))*t2(3)
1444  res(k) = res(k) + dot_product(t1(0:9), c0(51:60, k))*t2(4)
1445  res(k) = res(k) + dot_product(t1(0:8), c0(61:69, k))*t2(5)
1446  res(k) = res(k) + dot_product(t1(0:7), c0(70:77, k))*t2(6)
1447  res(k) = res(k) + dot_product(t1(0:6), c0(78:84, k))*t2(7)
1448  res(k) = res(k) + dot_product(t1(0:5), c0(85:90, k))*t2(8)
1449  res(k) = res(k) + dot_product(t1(0:4), c0(91:95, k))*t2(9)
1450  res(k) = res(k) + dot_product(t1(0:3), c0(96:99, k))*t2(10)
1451  res(k) = res(k) + dot_product(t1(0:2), c0(100:102, k))*t2(11)
1452  res(k) = res(k) + dot_product(t1(0:1), c0(103:104, k))*t2(12)
1453  res(k) = res(k) + dot_product(t1(0:0), c0(105:105, k))*t2(13)
1454  END DO
1455  END SUBROUTINE pd2val
1456 
1457 ! **************************************************************************************************
1458 !> \brief Returns the value of nderiv_init so that one can check if opening the potential file is
1459 !> worhtwhile
1460 !> \return ...
1461 !> \author A. Bussy, 05.2019
1462 ! **************************************************************************************************
1463  FUNCTION get_lmax_init() RESULT(lmax_init)
1464 
1465  INTEGER :: lmax_init
1466 
1467  lmax_init = nderiv_init
1468 
1469  END FUNCTION get_lmax_init
1470 
1471 END MODULE t_c_g0
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition: t_c_g0.F:57
t_c_g0::t_c_g0_n
subroutine, public t_c_g0_n(RES, use_gamma, R, T, NDERIV)
...
Definition: t_c_g0.F:84
t_c_g0::init
subroutine, public init(Nder, iunit, mepos, group)
...
Definition: t_c_g0.F:1357
t_c_g0::get_lmax_init
integer function, public get_lmax_init()
Returns the value of nderiv_init so that one can check if opening the potential file is worhtwhile.
Definition: t_c_g0.F:1464
t_c_g0::free_c0
subroutine, public free_c0()
...
Definition: t_c_g0.F:1388