(git:b279b6b)
t_sh_p_s_c.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, 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 ...
34 !> \author Joost VandeVondele and Manuel Guidon
35 !> \par History
36 !> Nov 2008 Joost VandeVondele and Manuel Guidon
37 !> Sep 2019 A.Bussy: Moved the file to common (together with gamma.F and t_c_g0.F)
38 ! **************************************************************************************************
39 MODULE t_sh_p_s_c
40  USE kinds, ONLY: dp
41  USE message_passing, ONLY: mp_comm_type
42 #include "../base/base_uses.f90"
43 
44  IMPLICIT NONE
45  PRIVATE
46 
47  PUBLIC :: trunc_cs_poly_n20, init, free_c0
48  REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE :: c0
49  INTEGER, PARAMETER :: degree = 15
50  REAL(KIND=dp), PARAMETER :: target_error = 0.100000e-08
51  INTEGER, PARAMETER :: nderiv_max = 21
52  INTEGER, SAVE :: nderiv_init = -1
53  INTEGER, SAVE :: patches = -1
54 
55 CONTAINS
56 
57 ! **************************************************************************************************
58 !> \brief ...
59 !> \param RES ...
60 !> \param R ...
61 !> \param T ...
62 !> \param NDERIV ...
63 ! **************************************************************************************************
64  SUBROUTINE trunc_cs_poly_n20(RES, R, T, NDERIV)
65  REAL(kind=dp), INTENT(OUT) :: res(*)
66  REAL(kind=dp), INTENT(IN) :: r, t
67  INTEGER, INTENT(IN) :: nderiv
68 
69  REAL(kind=dp), PARAMETER :: eps = 34.53877639491068526026987182_dp, n = 2.0_dp
70 
71  INTEGER :: i
72  REAL(kind=dp) :: lower, tg1, tg2, x1, x2
73 
74  IF (t <= 81.0_dp) THEN
75  x2 = 1/(1 + r)
76  x1 = 1/(1 + sqrt(t))
77  IF (x2 <= 0.500000000000000000e+00_dp) THEN
78  IF (x2 <= 0.250000000000000000e+00_dp) THEN
79  IF (x1 <= 0.550000000000000044e+00_dp) THEN
80  IF (x1 <= 0.325000000000000011e+00_dp) THEN
81  IF (x1 <= 0.212500000000000022e+00_dp) THEN
82  IF (x2 <= 0.125000000000000000e+00_dp) THEN
83  IF (x2 <= 0.625000000000000000e-01_dp) THEN
84  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.888888888888888751e+01_dp
85  tg2 = (2*x2 - 0.625000000000000000e-01_dp)*0.160000000000000000e+02_dp
86  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 1))
87  ELSE
88  IF (x1 <= 0.156250000000000000e+00_dp) THEN
89  IF (x2 <= 0.937500000000000000e-01_dp) THEN
90  IF (x1 <= 0.128124999999999989e+00_dp) THEN
91  tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
92  tg2 = (2*x2 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
93  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 2))
94  ELSE
95  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
96  tg2 = (2*x2 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
97  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 3))
98  END IF
99  ELSE
100  IF (x1 <= 0.128124999999999989e+00_dp) THEN
101  IF (x2 <= 0.109375000000000000e+00_dp) THEN
102  tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
103  tg2 = (2*x2 - 0.203125000000000000e+00_dp)*0.640000000000000000e+02_dp
104  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 4))
105  ELSE
106  tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
107  tg2 = (2*x2 - 0.234375000000000000e+00_dp)*0.640000000000000000e+02_dp
108  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 5))
109  END IF
110  ELSE
111  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
112  tg2 = (2*x2 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
113  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 6))
114  END IF
115  END IF
116  ELSE
117  tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
118  tg2 = (2*x2 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
119  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 7))
120  END IF
121  END IF
122  ELSE
123  IF (x2 <= 0.187500000000000000e+00_dp) THEN
124  IF (x1 <= 0.156250000000000000e+00_dp) THEN
125  IF (x2 <= 0.156250000000000000e+00_dp) THEN
126  IF (x1 <= 0.128124999999999989e+00_dp) THEN
127  tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
128  tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
129  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 8))
130  ELSE
131  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
132  tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
133  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 9))
134  END IF
135  ELSE
136  IF (x1 <= 0.128124999999999989e+00_dp) THEN
137  IF (x1 <= 0.114062499999999997e+00_dp) THEN
138  tg1 = (2*x1 - 0.214062499999999989e+00_dp)*0.711111111111111569e+02_dp
139  tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
140  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 10))
141  ELSE
142  tg1 = (2*x1 - 0.242187500000000000e+00_dp)*0.711111111111111569e+02_dp
143  tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
144  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 11))
145  END IF
146  ELSE
147  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
148  tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
149  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 12))
150  END IF
151  END IF
152  ELSE
153  IF (x2 <= 0.156250000000000000e+00_dp) THEN
154  tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
155  tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
156  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 13))
157  ELSE
158  tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
159  tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
160  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 14))
161  END IF
162  END IF
163  ELSE
164  IF (x1 <= 0.156250000000000000e+00_dp) THEN
165  IF (x1 <= 0.128124999999999989e+00_dp) THEN
166  IF (x2 <= 0.218750000000000000e+00_dp) THEN
167  IF (x1 <= 0.114062499999999997e+00_dp) THEN
168  tg1 = (2*x1 - 0.214062499999999989e+00_dp)*0.711111111111111569e+02_dp
169  tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
170  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 15))
171  ELSE
172  tg1 = (2*x1 - 0.242187500000000000e+00_dp)*0.711111111111111569e+02_dp
173  tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
174  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 16))
175  END IF
176  ELSE
177  tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
178  tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
179  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 17))
180  END IF
181  ELSE
182  IF (x2 <= 0.218750000000000000e+00_dp) THEN
183  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
184  tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
185  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 18))
186  ELSE
187  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
188  tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
189  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 19))
190  END IF
191  END IF
192  ELSE
193  IF (x2 <= 0.218750000000000000e+00_dp) THEN
194  tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
195  tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
196  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 20))
197  ELSE
198  tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
199  tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
200  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 21))
201  END IF
202  END IF
203  END IF
204  END IF
205  ELSE
206  IF (x2 <= 0.125000000000000000e+00_dp) THEN
207  tg1 = (2*x1 - 0.537500000000000089e+00_dp)*0.888888888888888928e+01_dp
208  tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
209  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 22))
210  ELSE
211  IF (x2 <= 0.187500000000000000e+00_dp) THEN
212  IF (x1 <= 0.268750000000000044e+00_dp) THEN
213  tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
214  tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
215  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 23))
216  ELSE
217  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
218  tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
219  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 24))
220  END IF
221  ELSE
222  IF (x1 <= 0.268750000000000044e+00_dp) THEN
223  tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
224  tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
225  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 25))
226  ELSE
227  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
228  tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
229  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 26))
230  END IF
231  END IF
232  END IF
233  END IF
234  ELSE
235  IF (x2 <= 0.125000000000000000e+00_dp) THEN
236  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.444444444444444375e+01_dp
237  tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
238  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 27))
239  ELSE
240  IF (x2 <= 0.187500000000000000e+00_dp) THEN
241  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.444444444444444375e+01_dp
242  tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
243  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 28))
244  ELSE
245  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.444444444444444375e+01_dp
246  tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
247  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 29))
248  END IF
249  END IF
250  END IF
251  ELSE
252  IF (x2 <= 0.125000000000000000e+00_dp) THEN
253  tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
254  tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
255  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 30))
256  ELSE
257  tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
258  tg2 = (2*x2 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
259  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 31))
260  END IF
261  END IF
262  ELSE
263  IF (x1 <= 0.550000000000000044e+00_dp) THEN
264  IF (x1 <= 0.325000000000000011e+00_dp) THEN
265  IF (x1 <= 0.212500000000000022e+00_dp) THEN
266  IF (x1 <= 0.156250000000000000e+00_dp) THEN
267  IF (x2 <= 0.375000000000000000e+00_dp) THEN
268  IF (x1 <= 0.128124999999999989e+00_dp) THEN
269  IF (x2 <= 0.312500000000000000e+00_dp) THEN
270  tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
271  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
272  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 32))
273  ELSE
274  tg1 = (2*x1 - 0.228124999999999994e+00_dp)*0.355555555555555785e+02_dp
275  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
276  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 33))
277  END IF
278  ELSE
279  IF (x2 <= 0.312500000000000000e+00_dp) THEN
280  IF (x2 <= 0.281250000000000000e+00_dp) THEN
281  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
282  tg2 = (2*x2 - 0.531250000000000000e+00_dp)*0.320000000000000000e+02_dp
283  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 34))
284  ELSE
285  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
286  tg2 = (2*x2 - 0.593750000000000000e+00_dp)*0.320000000000000000e+02_dp
287  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 35))
288  END IF
289  ELSE
290  tg1 = (2*x1 - 0.284374999999999989e+00_dp)*0.355555555555555429e+02_dp
291  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
292  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 36))
293  END IF
294  END IF
295  ELSE
296  tg1 = (2*x1 - 0.256249999999999978e+00_dp)*0.177777777777777786e+02_dp
297  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
298  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 37))
299  END IF
300  ELSE
301  IF (x2 <= 0.375000000000000000e+00_dp) THEN
302  IF (x2 <= 0.312500000000000000e+00_dp) THEN
303  IF (x1 <= 0.184375000000000011e+00_dp) THEN
304  tg1 = (2*x1 - 0.340625000000000011e+00_dp)*0.355555555555555429e+02_dp
305  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
306  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 38))
307  ELSE
308  tg1 = (2*x1 - 0.396875000000000033e+00_dp)*0.355555555555555429e+02_dp
309  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
310  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 39))
311  END IF
312  ELSE
313  IF (x1 <= 0.184375000000000011e+00_dp) THEN
314  tg1 = (2*x1 - 0.340625000000000011e+00_dp)*0.355555555555555429e+02_dp
315  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
316  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 40))
317  ELSE
318  tg1 = (2*x1 - 0.396875000000000033e+00_dp)*0.355555555555555429e+02_dp
319  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
320  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 41))
321  END IF
322  END IF
323  ELSE
324  IF (x1 <= 0.184375000000000011e+00_dp) THEN
325  tg1 = (2*x1 - 0.340625000000000011e+00_dp)*0.355555555555555429e+02_dp
326  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
327  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 42))
328  ELSE
329  tg1 = (2*x1 - 0.396875000000000033e+00_dp)*0.355555555555555429e+02_dp
330  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
331  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 43))
332  END IF
333  END IF
334  END IF
335  ELSE
336  IF (x2 <= 0.375000000000000000e+00_dp) THEN
337  IF (x2 <= 0.312500000000000000e+00_dp) THEN
338  IF (x1 <= 0.268750000000000044e+00_dp) THEN
339  tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
340  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
341  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 44))
342  ELSE
343  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
344  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
345  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 45))
346  END IF
347  ELSE
348  IF (x1 <= 0.268750000000000044e+00_dp) THEN
349  tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
350  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
351  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 46))
352  ELSE
353  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
354  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
355  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 47))
356  END IF
357  END IF
358  ELSE
359  IF (x1 <= 0.268750000000000044e+00_dp) THEN
360  tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
361  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
362  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 48))
363  ELSE
364  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
365  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
366  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 49))
367  END IF
368  END IF
369  END IF
370  ELSE
371  IF (x2 <= 0.375000000000000000e+00_dp) THEN
372  IF (x1 <= 0.437500000000000000e+00_dp) THEN
373  tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
374  tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
375  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 50))
376  ELSE
377  tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
378  tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
379  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 51))
380  END IF
381  ELSE
382  IF (x2 <= 0.437500000000000000e+00_dp) THEN
383  IF (x1 <= 0.437500000000000000e+00_dp) THEN
384  tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
385  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
386  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 52))
387  ELSE
388  tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
389  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
390  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 53))
391  END IF
392  ELSE
393  IF (x1 <= 0.437500000000000000e+00_dp) THEN
394  tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
395  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
396  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 54))
397  ELSE
398  tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
399  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
400  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 55))
401  END IF
402  END IF
403  END IF
404  END IF
405  ELSE
406  IF (x2 <= 0.375000000000000000e+00_dp) THEN
407  IF (x2 <= 0.312500000000000000e+00_dp) THEN
408  tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
409  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
410  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 56))
411  ELSE
412  tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
413  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
414  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 57))
415  END IF
416  ELSE
417  IF (x2 <= 0.437500000000000000e+00_dp) THEN
418  tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
419  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
420  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 58))
421  ELSE
422  IF (x1 <= 0.775000000000000022e+00_dp) THEN
423  tg1 = (2*x1 - 0.132500000000000018e+01_dp)*0.444444444444444464e+01_dp
424  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
425  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 59))
426  ELSE
427  tg1 = (2*x1 - 0.177499999999999991e+01_dp)*0.444444444444444464e+01_dp
428  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
429  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 60))
430  END IF
431  END IF
432  END IF
433  END IF
434  END IF
435  ELSE
436  IF (x1 <= 0.550000000000000044e+00_dp) THEN
437  IF (x1 <= 0.325000000000000011e+00_dp) THEN
438  IF (x1 <= 0.212500000000000022e+00_dp) THEN
439  IF (x2 <= 0.750000000000000000e+00_dp) THEN
440  IF (x1 <= 0.156250000000000000e+00_dp) THEN
441  tg1 = (2*x1 - 0.256249999999999978e+00_dp)*0.177777777777777786e+02_dp
442  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
443  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 61))
444  ELSE
445  tg1 = (2*x1 - 0.368750000000000022e+00_dp)*0.177777777777777715e+02_dp
446  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
447  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 62))
448  END IF
449  ELSE
450  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.888888888888888751e+01_dp
451  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
452  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 63))
453  END IF
454  ELSE
455  IF (x2 <= 0.750000000000000000e+00_dp) THEN
456  IF (x2 <= 0.625000000000000000e+00_dp) THEN
457  IF (x1 <= 0.268750000000000044e+00_dp) THEN
458  tg1 = (2*x1 - 0.481250000000000067e+00_dp)*0.177777777777777715e+02_dp
459  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
460  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 64))
461  ELSE
462  tg1 = (2*x1 - 0.593750000000000000e+00_dp)*0.177777777777777892e+02_dp
463  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
464  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 65))
465  END IF
466  ELSE
467  tg1 = (2*x1 - 0.537500000000000089e+00_dp)*0.888888888888888928e+01_dp
468  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
469  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 66))
470  END IF
471  ELSE
472  tg1 = (2*x1 - 0.537500000000000089e+00_dp)*0.888888888888888928e+01_dp
473  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
474  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 67))
475  END IF
476  END IF
477  ELSE
478  IF (x2 <= 0.750000000000000000e+00_dp) THEN
479  IF (x2 <= 0.625000000000000000e+00_dp) THEN
480  IF (x1 <= 0.437500000000000000e+00_dp) THEN
481  tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
482  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
483  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 68))
484  ELSE
485  tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
486  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
487  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 69))
488  END IF
489  ELSE
490  IF (x1 <= 0.437500000000000000e+00_dp) THEN
491  tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
492  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
493  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 70))
494  ELSE
495  tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
496  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
497  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 71))
498  END IF
499  END IF
500  ELSE
501  IF (x1 <= 0.437500000000000000e+00_dp) THEN
502  tg1 = (2*x1 - 0.762499999999999956e+00_dp)*0.888888888888888928e+01_dp
503  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
504  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 72))
505  ELSE
506  tg1 = (2*x1 - 0.987500000000000044e+00_dp)*0.888888888888888573e+01_dp
507  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
508  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 73))
509  END IF
510  END IF
511  END IF
512  ELSE
513  IF (x2 <= 0.750000000000000000e+00_dp) THEN
514  IF (x2 <= 0.625000000000000000e+00_dp) THEN
515  IF (x1 <= 0.775000000000000022e+00_dp) THEN
516  tg1 = (2*x1 - 0.132500000000000018e+01_dp)*0.444444444444444464e+01_dp
517  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
518  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 74))
519  ELSE
520  tg1 = (2*x1 - 0.177499999999999991e+01_dp)*0.444444444444444464e+01_dp
521  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
522  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 75))
523  END IF
524  ELSE
525  IF (x1 <= 0.775000000000000022e+00_dp) THEN
526  tg1 = (2*x1 - 0.132500000000000018e+01_dp)*0.444444444444444464e+01_dp
527  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
528  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 76))
529  ELSE
530  tg1 = (2*x1 - 0.177499999999999991e+01_dp)*0.444444444444444464e+01_dp
531  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
532  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 77))
533  END IF
534  END IF
535  ELSE
536  IF (x2 <= 0.875000000000000000e+00_dp) THEN
537  tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
538  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
539  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 78))
540  ELSE
541  tg1 = (2*x1 - 0.155000000000000004e+01_dp)*0.222222222222222232e+01_dp
542  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
543  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 79))
544  END IF
545  END IF
546  END IF
547  END IF
548  ELSE
549  x1 = sqrt(81.0_dp/t)
550  lower = (sqrt(t) - sqrt(eps))/n
551  IF (r >= lower) THEN
552  x2 = lower/r
553  IF (x2 <= 0.500000000000000000e+00_dp) THEN
554  IF (x1 <= 0.500000000000000000e+00_dp) THEN
555  IF (x2 <= 0.250000000000000000e+00_dp) THEN
556  tg1 = (2*x1 - 0.500000000000000000e+00_dp)*0.200000000000000000e+01_dp
557  tg2 = (2*x2 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
558  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 80))
559  ELSE
560  IF (x2 <= 0.375000000000000000e+00_dp) THEN
561  IF (x1 <= 0.250000000000000000e+00_dp) THEN
562  tg1 = (2*x1 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
563  tg2 = (2*x2 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
564  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 81))
565  ELSE
566  IF (x2 <= 0.312500000000000000e+00_dp) THEN
567  tg1 = (2*x1 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
568  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
569  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 82))
570  ELSE
571  IF (x1 <= 0.375000000000000000e+00_dp) THEN
572  tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
573  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
574  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 83))
575  ELSE
576  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
577  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
578  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 84))
579  END IF
580  END IF
581  END IF
582  ELSE
583  IF (x1 <= 0.250000000000000000e+00_dp) THEN
584  IF (x2 <= 0.437500000000000000e+00_dp) THEN
585  IF (x1 <= 0.125000000000000000e+00_dp) THEN
586  tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
587  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
588  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 85))
589  ELSE
590  tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
591  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
592  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 86))
593  END IF
594  ELSE
595  IF (x1 <= 0.125000000000000000e+00_dp) THEN
596  IF (x2 <= 0.468750000000000000e+00_dp) THEN
597  tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
598  tg2 = (2*x2 - 0.906250000000000000e+00_dp)*0.320000000000000000e+02_dp
599  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 87))
600  ELSE
601  tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
602  tg2 = (2*x2 - 0.968750000000000000e+00_dp)*0.320000000000000000e+02_dp
603  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 88))
604  END IF
605  ELSE
606  IF (x2 <= 0.468750000000000000e+00_dp) THEN
607  tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
608  tg2 = (2*x2 - 0.906250000000000000e+00_dp)*0.320000000000000000e+02_dp
609  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 89))
610  ELSE
611  tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
612  tg2 = (2*x2 - 0.968750000000000000e+00_dp)*0.320000000000000000e+02_dp
613  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 90))
614  END IF
615  END IF
616  END IF
617  ELSE
618  IF (x2 <= 0.437500000000000000e+00_dp) THEN
619  IF (x2 <= 0.406250000000000000e+00_dp) THEN
620  IF (x1 <= 0.375000000000000000e+00_dp) THEN
621  tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
622  tg2 = (2*x2 - 0.781250000000000000e+00_dp)*0.320000000000000000e+02_dp
623  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 91))
624  ELSE
625  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
626  tg2 = (2*x2 - 0.781250000000000000e+00_dp)*0.320000000000000000e+02_dp
627  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 92))
628  END IF
629  ELSE
630  tg1 = (2*x1 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
631  tg2 = (2*x2 - 0.843750000000000000e+00_dp)*0.320000000000000000e+02_dp
632  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 93))
633  END IF
634  ELSE
635  tg1 = (2*x1 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
636  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
637  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 94))
638  END IF
639  END IF
640  END IF
641  END IF
642  ELSE
643  IF (x2 <= 0.250000000000000000e+00_dp) THEN
644  IF (x2 <= 0.125000000000000000e+00_dp) THEN
645  tg1 = (2*x1 - 0.150000000000000000e+01_dp)*0.200000000000000000e+01_dp
646  tg2 = (2*x2 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
647  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 95))
648  ELSE
649  IF (x2 <= 0.187500000000000000e+00_dp) THEN
650  IF (x1 <= 0.750000000000000000e+00_dp) THEN
651  tg1 = (2*x1 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
652  tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
653  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 96))
654  ELSE
655  IF (x1 <= 0.875000000000000000e+00_dp) THEN
656  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
657  tg2 = (2*x2 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
658  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 97))
659  ELSE
660  IF (x2 <= 0.156250000000000000e+00_dp) THEN
661  tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
662  tg2 = (2*x2 - 0.281250000000000000e+00_dp)*0.320000000000000000e+02_dp
663  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 98))
664  ELSE
665  tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
666  tg2 = (2*x2 - 0.343750000000000000e+00_dp)*0.320000000000000000e+02_dp
667  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 99))
668  END IF
669  END IF
670  END IF
671  ELSE
672  IF (x1 <= 0.750000000000000000e+00_dp) THEN
673  IF (x1 <= 0.625000000000000000e+00_dp) THEN
674  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
675  tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
676  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 100))
677  ELSE
678  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
679  tg2 = (2*x2 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
680  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 101))
681  END IF
682  ELSE
683  IF (x1 <= 0.875000000000000000e+00_dp) THEN
684  IF (x2 <= 0.218750000000000000e+00_dp) THEN
685  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
686  tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
687  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 102))
688  ELSE
689  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
690  tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
691  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 103))
692  END IF
693  ELSE
694  IF (x2 <= 0.218750000000000000e+00_dp) THEN
695  tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
696  tg2 = (2*x2 - 0.406250000000000000e+00_dp)*0.320000000000000000e+02_dp
697  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 104))
698  ELSE
699  tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
700  tg2 = (2*x2 - 0.468750000000000000e+00_dp)*0.320000000000000000e+02_dp
701  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 105))
702  END IF
703  END IF
704  END IF
705  END IF
706  END IF
707  ELSE
708  IF (x2 <= 0.375000000000000000e+00_dp) THEN
709  IF (x1 <= 0.750000000000000000e+00_dp) THEN
710  IF (x2 <= 0.312500000000000000e+00_dp) THEN
711  IF (x1 <= 0.625000000000000000e+00_dp) THEN
712  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
713  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
714  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 106))
715  ELSE
716  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
717  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
718  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 107))
719  END IF
720  ELSE
721  IF (x1 <= 0.625000000000000000e+00_dp) THEN
722  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
723  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
724  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 108))
725  ELSE
726  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
727  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
728  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 109))
729  END IF
730  END IF
731  ELSE
732  IF (x2 <= 0.312500000000000000e+00_dp) THEN
733  IF (x1 <= 0.875000000000000000e+00_dp) THEN
734  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
735  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
736  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 110))
737  ELSE
738  tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
739  tg2 = (2*x2 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
740  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 111))
741  END IF
742  ELSE
743  IF (x1 <= 0.875000000000000000e+00_dp) THEN
744  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
745  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
746  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 112))
747  ELSE
748  IF (x1 <= 0.937500000000000000e+00_dp) THEN
749  tg1 = (2*x1 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
750  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
751  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 113))
752  ELSE
753  tg1 = (2*x1 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
754  tg2 = (2*x2 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
755  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 114))
756  END IF
757  END IF
758  END IF
759  END IF
760  ELSE
761  IF (x1 <= 0.750000000000000000e+00_dp) THEN
762  IF (x1 <= 0.625000000000000000e+00_dp) THEN
763  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
764  tg2 = (2*x2 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
765  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 115))
766  ELSE
767  IF (x2 <= 0.437500000000000000e+00_dp) THEN
768  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
769  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
770  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 116))
771  ELSE
772  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
773  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
774  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 117))
775  END IF
776  END IF
777  ELSE
778  IF (x1 <= 0.875000000000000000e+00_dp) THEN
779  IF (x2 <= 0.437500000000000000e+00_dp) THEN
780  IF (x1 <= 0.812500000000000000e+00_dp) THEN
781  tg1 = (2*x1 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
782  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
783  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 118))
784  ELSE
785  tg1 = (2*x1 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
786  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
787  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 119))
788  END IF
789  ELSE
790  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
791  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
792  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 120))
793  END IF
794  ELSE
795  IF (x2 <= 0.437500000000000000e+00_dp) THEN
796  IF (x1 <= 0.937500000000000000e+00_dp) THEN
797  tg1 = (2*x1 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
798  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
799  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 121))
800  ELSE
801  tg1 = (2*x1 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
802  tg2 = (2*x2 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
803  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 122))
804  END IF
805  ELSE
806  tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
807  tg2 = (2*x2 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
808  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 123))
809  END IF
810  END IF
811  END IF
812  END IF
813  END IF
814  END IF
815  ELSE
816  IF (x1 <= 0.500000000000000000e+00_dp) THEN
817  IF (x2 <= 0.750000000000000000e+00_dp) THEN
818  IF (x1 <= 0.250000000000000000e+00_dp) THEN
819  tg1 = (2*x1 - 0.250000000000000000e+00_dp)*0.400000000000000000e+01_dp
820  tg2 = (2*x2 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
821  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 124))
822  ELSE
823  IF (x2 <= 0.625000000000000000e+00_dp) THEN
824  IF (x1 <= 0.375000000000000000e+00_dp) THEN
825  tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
826  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
827  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 125))
828  ELSE
829  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
830  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
831  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 126))
832  END IF
833  ELSE
834  IF (x1 <= 0.375000000000000000e+00_dp) THEN
835  IF (x2 <= 0.687500000000000000e+00_dp) THEN
836  tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
837  tg2 = (2*x2 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
838  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 127))
839  ELSE
840  tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
841  tg2 = (2*x2 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
842  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 128))
843  END IF
844  ELSE
845  IF (x2 <= 0.687500000000000000e+00_dp) THEN
846  IF (x1 <= 0.437500000000000000e+00_dp) THEN
847  tg1 = (2*x1 - 0.812500000000000000e+00_dp)*0.160000000000000000e+02_dp
848  tg2 = (2*x2 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
849  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 129))
850  ELSE
851  tg1 = (2*x1 - 0.937500000000000000e+00_dp)*0.160000000000000000e+02_dp
852  tg2 = (2*x2 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
853  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 130))
854  END IF
855  ELSE
856  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
857  tg2 = (2*x2 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
858  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 131))
859  END IF
860  END IF
861  END IF
862  END IF
863  ELSE
864  IF (x1 <= 0.250000000000000000e+00_dp) THEN
865  IF (x2 <= 0.875000000000000000e+00_dp) THEN
866  IF (x1 <= 0.125000000000000000e+00_dp) THEN
867  tg1 = (2*x1 - 0.125000000000000000e+00_dp)*0.800000000000000000e+01_dp
868  tg2 = (2*x2 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
869  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 132))
870  ELSE
871  IF (x2 <= 0.812500000000000000e+00_dp) THEN
872  tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
873  tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
874  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 133))
875  ELSE
876  IF (x1 <= 0.187500000000000000e+00_dp) THEN
877  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
878  tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
879  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 134))
880  ELSE
881  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
882  tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
883  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 135))
884  END IF
885  END IF
886  END IF
887  ELSE
888  IF (x1 <= 0.125000000000000000e+00_dp) THEN
889  IF (x2 <= 0.937500000000000000e+00_dp) THEN
890  IF (x1 <= 0.625000000000000000e-01_dp) THEN
891  tg1 = (2*x1 - 0.625000000000000000e-01_dp)*0.160000000000000000e+02_dp
892  tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
893  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 136))
894  ELSE
895  IF (x2 <= 0.906250000000000000e+00_dp) THEN
896  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
897  tg2 = (2*x2 - 0.178125000000000000e+01_dp)*0.320000000000000000e+02_dp
898  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 137))
899  ELSE
900  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
901  tg2 = (2*x2 - 0.184375000000000000e+01_dp)*0.320000000000000000e+02_dp
902  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 138))
903  END IF
904  END IF
905  ELSE
906  IF (x1 <= 0.625000000000000000e-01_dp) THEN
907  IF (x2 <= 0.968750000000000000e+00_dp) THEN
908  IF (x1 <= 0.312500000000000000e-01_dp) THEN
909  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
910  tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
911  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 139))
912  ELSE
913  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
914  tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
915  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 140))
916  END IF
917  ELSE
918  IF (x1 <= 0.312500000000000000e-01_dp) THEN
919  IF (x2 <= 0.984375000000000000e+00_dp) THEN
920  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
921  tg2 = (2*x2 - 0.195312500000000000e+01_dp)*0.640000000000000000e+02_dp
922  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 141))
923  ELSE
924  tg1 = (2*x1 - 0.312500000000000000e-01_dp)*0.320000000000000000e+02_dp
925  tg2 = (2*x2 - 0.198437500000000000e+01_dp)*0.640000000000000000e+02_dp
926  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 142))
927  END IF
928  ELSE
929  tg1 = (2*x1 - 0.937500000000000000e-01_dp)*0.320000000000000000e+02_dp
930  tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
931  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 143))
932  END IF
933  END IF
934  ELSE
935  IF (x2 <= 0.968750000000000000e+00_dp) THEN
936  IF (x1 <= 0.937500000000000000e-01_dp) THEN
937  tg1 = (2*x1 - 0.156250000000000000e+00_dp)*0.320000000000000000e+02_dp
938  tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
939  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 144))
940  ELSE
941  tg1 = (2*x1 - 0.218750000000000000e+00_dp)*0.320000000000000000e+02_dp
942  tg2 = (2*x2 - 0.190625000000000000e+01_dp)*0.320000000000000000e+02_dp
943  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 145))
944  END IF
945  ELSE
946  tg1 = (2*x1 - 0.187500000000000000e+00_dp)*0.160000000000000000e+02_dp
947  tg2 = (2*x2 - 0.196875000000000000e+01_dp)*0.320000000000000000e+02_dp
948  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 146))
949  END IF
950  END IF
951  END IF
952  ELSE
953  IF (x2 <= 0.937500000000000000e+00_dp) THEN
954  IF (x1 <= 0.187500000000000000e+00_dp) THEN
955  IF (x2 <= 0.906250000000000000e+00_dp) THEN
956  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
957  tg2 = (2*x2 - 0.178125000000000000e+01_dp)*0.320000000000000000e+02_dp
958  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 147))
959  ELSE
960  tg1 = (2*x1 - 0.312500000000000000e+00_dp)*0.160000000000000000e+02_dp
961  tg2 = (2*x2 - 0.184375000000000000e+01_dp)*0.320000000000000000e+02_dp
962  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 148))
963  END IF
964  ELSE
965  tg1 = (2*x1 - 0.437500000000000000e+00_dp)*0.160000000000000000e+02_dp
966  tg2 = (2*x2 - 0.181250000000000000e+01_dp)*0.160000000000000000e+02_dp
967  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 149))
968  END IF
969  ELSE
970  tg1 = (2*x1 - 0.375000000000000000e+00_dp)*0.800000000000000000e+01_dp
971  tg2 = (2*x2 - 0.193750000000000000e+01_dp)*0.160000000000000000e+02_dp
972  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 150))
973  END IF
974  END IF
975  END IF
976  ELSE
977  IF (x2 <= 0.875000000000000000e+00_dp) THEN
978  IF (x1 <= 0.375000000000000000e+00_dp) THEN
979  IF (x2 <= 0.812500000000000000e+00_dp) THEN
980  IF (x1 <= 0.312500000000000000e+00_dp) THEN
981  tg1 = (2*x1 - 0.562500000000000000e+00_dp)*0.160000000000000000e+02_dp
982  tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
983  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 151))
984  ELSE
985  tg1 = (2*x1 - 0.687500000000000000e+00_dp)*0.160000000000000000e+02_dp
986  tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
987  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 152))
988  END IF
989  ELSE
990  tg1 = (2*x1 - 0.625000000000000000e+00_dp)*0.800000000000000000e+01_dp
991  tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
992  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 153))
993  END IF
994  ELSE
995  IF (x2 <= 0.812500000000000000e+00_dp) THEN
996  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
997  tg2 = (2*x2 - 0.156250000000000000e+01_dp)*0.160000000000000000e+02_dp
998  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 154))
999  ELSE
1000  tg1 = (2*x1 - 0.875000000000000000e+00_dp)*0.800000000000000000e+01_dp
1001  tg2 = (2*x2 - 0.168750000000000000e+01_dp)*0.160000000000000000e+02_dp
1002  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 155))
1003  END IF
1004  END IF
1005  ELSE
1006  tg1 = (2*x1 - 0.750000000000000000e+00_dp)*0.400000000000000000e+01_dp
1007  tg2 = (2*x2 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1008  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 156))
1009  END IF
1010  END IF
1011  END IF
1012  ELSE
1013  IF (x2 <= 0.750000000000000000e+00_dp) THEN
1014  IF (x1 <= 0.750000000000000000e+00_dp) THEN
1015  IF (x2 <= 0.625000000000000000e+00_dp) THEN
1016  IF (x1 <= 0.625000000000000000e+00_dp) THEN
1017  IF (x2 <= 0.562500000000000000e+00_dp) THEN
1018  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1019  tg2 = (2*x2 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1020  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 157))
1021  ELSE
1022  IF (x1 <= 0.562500000000000000e+00_dp) THEN
1023  tg1 = (2*x1 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1024  tg2 = (2*x2 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1025  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 158))
1026  ELSE
1027  tg1 = (2*x1 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1028  tg2 = (2*x2 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1029  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 159))
1030  END IF
1031  END IF
1032  ELSE
1033  IF (x2 <= 0.562500000000000000e+00_dp) THEN
1034  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1035  tg2 = (2*x2 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1036  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 160))
1037  ELSE
1038  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1039  tg2 = (2*x2 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1040  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 161))
1041  END IF
1042  END IF
1043  ELSE
1044  IF (x1 <= 0.625000000000000000e+00_dp) THEN
1045  IF (x2 <= 0.687500000000000000e+00_dp) THEN
1046  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1047  tg2 = (2*x2 - 0.131250000000000000e+01_dp)*0.160000000000000000e+02_dp
1048  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 162))
1049  ELSE
1050  tg1 = (2*x1 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1051  tg2 = (2*x2 - 0.143750000000000000e+01_dp)*0.160000000000000000e+02_dp
1052  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 163))
1053  END IF
1054  ELSE
1055  tg1 = (2*x1 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1056  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1057  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 164))
1058  END IF
1059  END IF
1060  ELSE
1061  IF (x2 <= 0.625000000000000000e+00_dp) THEN
1062  IF (x1 <= 0.875000000000000000e+00_dp) THEN
1063  IF (x2 <= 0.562500000000000000e+00_dp) THEN
1064  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1065  tg2 = (2*x2 - 0.106250000000000000e+01_dp)*0.160000000000000000e+02_dp
1066  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 165))
1067  ELSE
1068  tg1 = (2*x1 - 0.162500000000000000e+01_dp)*0.800000000000000000e+01_dp
1069  tg2 = (2*x2 - 0.118750000000000000e+01_dp)*0.160000000000000000e+02_dp
1070  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 166))
1071  END IF
1072  ELSE
1073  tg1 = (2*x1 - 0.187500000000000000e+01_dp)*0.800000000000000000e+01_dp
1074  tg2 = (2*x2 - 0.112500000000000000e+01_dp)*0.800000000000000000e+01_dp
1075  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 167))
1076  END IF
1077  ELSE
1078  tg1 = (2*x1 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1079  tg2 = (2*x2 - 0.137500000000000000e+01_dp)*0.800000000000000000e+01_dp
1080  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 168))
1081  END IF
1082  END IF
1083  ELSE
1084  IF (x1 <= 0.750000000000000000e+00_dp) THEN
1085  tg1 = (2*x1 - 0.125000000000000000e+01_dp)*0.400000000000000000e+01_dp
1086  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1087  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 169))
1088  ELSE
1089  tg1 = (2*x1 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1090  tg2 = (2*x2 - 0.175000000000000000e+01_dp)*0.400000000000000000e+01_dp
1091  CALL pd2val(res, nderiv, tg1, tg2, c0(1, 170))
1092  END IF
1093  END IF
1094  END IF
1095  END IF
1096  ELSE
1097  DO i = 1, nderiv + 1
1098  res(i) = 0.0_dp
1099  END DO
1100  END IF
1101  END IF
1102  END SUBROUTINE trunc_cs_poly_n20
1103 
1104 ! **************************************************************************************************
1105 !> \brief ...
1106 !> \param Nder the number of derivatives that will actually be used
1107 !> \param iunit contains the data file to initialize the table
1108 !> \param mepos ...
1109 !> \param group ...
1110 ! **************************************************************************************************
1111  SUBROUTINE init(Nder, iunit, mepos, group)
1112  INTEGER, INTENT(IN) :: nder, iunit, mepos
1113 
1114  CLASS(mp_comm_type), INTENT(IN) :: group
1115 
1116  INTEGER :: i
1117  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: chunk
1118 
1119  patches = 170
1120  IF (nder > nderiv_max) cpabort("Reading data for initialization of C0 failed")
1121  nderiv_init = nder
1122  IF (ALLOCATED(c0)) DEALLOCATE (c0)
1123  ! round up to a multiple of 32 to give some generous alignment for each C0
1124  ALLOCATE (c0(32*((31 + (nder + 1)*(degree + 1)*(degree + 2)/2)/32), patches))
1125  ! valgrind workaround
1126  c0 = huge(0.0_dp)
1127  IF (mepos == 0) THEN
1128  ALLOCATE (chunk((nderiv_max + 1)*(degree + 1)*(degree + 2)/2))
1129  DO i = 1, patches
1130  READ (iunit, *) chunk
1131  c0(1:(nder + 1)*(degree + 1)*(degree + 2)/2, i) = chunk(1:(nder + 1)*(degree + 1)*(degree + 2)/2)
1132  END DO
1133  DEALLOCATE (chunk)
1134  END IF
1135  CALL group%bcast(c0, 0)
1136  END SUBROUTINE init
1137 
1138 ! **************************************************************************************************
1139 !> \brief ...
1140 ! **************************************************************************************************
1141  SUBROUTINE free_c0()
1142  IF (ALLOCATED(c0)) DEALLOCATE (c0)
1143  nderiv_init = -1
1144  END SUBROUTINE free_c0
1145 
1146 ! **************************************************************************************************
1147 !> \brief ...
1148 !> \param RES ...
1149 !> \param NDERIV ...
1150 !> \param TG1 ...
1151 !> \param TG2 ...
1152 !> \param C0 ...
1153 ! **************************************************************************************************
1154  SUBROUTINE pd2val(RES, NDERIV, TG1, TG2, C0)
1155  REAL(kind=dp), INTENT(OUT) :: res(*)
1156  INTEGER, INTENT(IN) :: nderiv
1157  REAL(kind=dp), INTENT(IN) :: tg1, tg2, c0(136, *)
1158 
1159  REAL(kind=dp), PARAMETER :: sqrt2 = 1.4142135623730950488016887242096980785696718753_dp
1160 
1161  INTEGER :: k
1162  REAL(kind=dp) :: t1(0:15), t2(0:15)
1163 
1164  t1(0) = 1.0_dp
1165  t2(0) = 1.0_dp
1166  t1(1) = sqrt2*tg1
1167  t2(1) = sqrt2*tg2
1168  t1(2) = 2*tg1*t1(1) - sqrt2
1169  t2(2) = 2*tg2*t2(1) - sqrt2
1170  t1(3) = 2*tg1*t1(2) - t1(1)
1171  t2(3) = 2*tg2*t2(2) - t2(1)
1172  t1(4) = 2*tg1*t1(3) - t1(2)
1173  t2(4) = 2*tg2*t2(3) - t2(2)
1174  t1(5) = 2*tg1*t1(4) - t1(3)
1175  t2(5) = 2*tg2*t2(4) - t2(3)
1176  t1(6) = 2*tg1*t1(5) - t1(4)
1177  t2(6) = 2*tg2*t2(5) - t2(4)
1178  t1(7) = 2*tg1*t1(6) - t1(5)
1179  t2(7) = 2*tg2*t2(6) - t2(5)
1180  t1(8) = 2*tg1*t1(7) - t1(6)
1181  t2(8) = 2*tg2*t2(7) - t2(6)
1182  t1(9) = 2*tg1*t1(8) - t1(7)
1183  t2(9) = 2*tg2*t2(8) - t2(7)
1184  t1(10) = 2*tg1*t1(9) - t1(8)
1185  t2(10) = 2*tg2*t2(9) - t2(8)
1186  t1(11) = 2*tg1*t1(10) - t1(9)
1187  t2(11) = 2*tg2*t2(10) - t2(9)
1188  t1(12) = 2*tg1*t1(11) - t1(10)
1189  t2(12) = 2*tg2*t2(11) - t2(10)
1190  t1(13) = 2*tg1*t1(12) - t1(11)
1191  t2(13) = 2*tg2*t2(12) - t2(11)
1192  t1(14) = 2*tg1*t1(13) - t1(12)
1193  t2(14) = 2*tg2*t2(13) - t2(12)
1194  t1(15) = 2*tg1*t1(14) - t1(13)
1195  t2(15) = 2*tg2*t2(14) - t2(13)
1196  DO k = 1, nderiv + 1
1197  res(k) = 0.0_dp
1198  res(k) = res(k) + dot_product(t1(0:15), c0(1:16, k))*t2(0)
1199  res(k) = res(k) + dot_product(t1(0:14), c0(17:31, k))*t2(1)
1200  res(k) = res(k) + dot_product(t1(0:13), c0(32:45, k))*t2(2)
1201  res(k) = res(k) + dot_product(t1(0:12), c0(46:58, k))*t2(3)
1202  res(k) = res(k) + dot_product(t1(0:11), c0(59:70, k))*t2(4)
1203  res(k) = res(k) + dot_product(t1(0:10), c0(71:81, k))*t2(5)
1204  res(k) = res(k) + dot_product(t1(0:9), c0(82:91, k))*t2(6)
1205  res(k) = res(k) + dot_product(t1(0:8), c0(92:100, k))*t2(7)
1206  res(k) = res(k) + dot_product(t1(0:7), c0(101:108, k))*t2(8)
1207  res(k) = res(k) + dot_product(t1(0:6), c0(109:115, k))*t2(9)
1208  res(k) = res(k) + dot_product(t1(0:5), c0(116:121, k))*t2(10)
1209  res(k) = res(k) + dot_product(t1(0:4), c0(122:126, k))*t2(11)
1210  res(k) = res(k) + dot_product(t1(0:3), c0(127:130, k))*t2(12)
1211  res(k) = res(k) + dot_product(t1(0:2), c0(131:133, k))*t2(13)
1212  res(k) = res(k) + dot_product(t1(0:1), c0(134:135, k))*t2(14)
1213  res(k) = res(k) + dot_product(t1(0:0), c0(136:136, k))*t2(15)
1214  END DO
1215  END SUBROUTINE pd2val
1216 
1217 END MODULE t_sh_p_s_c
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_sh_p_s_c
...
Definition: t_sh_p_s_c.F:39
t_sh_p_s_c::trunc_cs_poly_n20
subroutine, public trunc_cs_poly_n20(RES, R, T, NDERIV)
...
Definition: t_sh_p_s_c.F:65
t_sh_p_s_c::free_c0
subroutine, public free_c0()
...
Definition: t_sh_p_s_c.F:1142
t_sh_p_s_c::init
subroutine, public init(Nder, iunit, mepos, group)
...
Definition: t_sh_p_s_c.F:1112