Typo fixes
[scilab.git] / scilab / modules / interpolation / src / fortran / somespline.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) Bruno Pincon
3 c
4 c This file must be used under the terms of the CeCILL.
5 c This source file is licensed as described in the file COPYING, which
6 c you should have received as part of this distribution.  The terms
7 c are also available at
8 c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
9
10       subroutine SplineCub(x, y, d, n, type, A_d, A_sd, qdy, lll)
11 *
12 *     PURPOSE
13 *        computes a cubic spline interpolation function
14 *        in Hermite form (ie computes the derivatives d(i) of the
15 *        spline in each interpolation point (x(i), y(i)))
16 *
17 *     ARGUMENTS
18 *      inputs :
19 *         n       number of interpolation points (n >= 3)
20 *         x, y    the n interpolation points, x must be in strict increasing order
21 *         type    type of the spline : currently
22 *                    type = 0 correspond to a  NOT_A_KNOT spline where it is
23 *                             imposed the conditions :
24 *                                 s'''(x(2)-) = s'''(x(2)+)
25 *                             and s'''(x(n-1)-) = s'''(x(n-1)+)
26 *                    type = 1 correspond to a NATURAL spline with the conditions :
27 *                                 s''(x1) = 0
28 *                             and s''(xn) = 0
29 *                    type = 2 correspond to a CLAMPED spline (d(1) and d(n) are given)
30 *
31 *                    type = 3 correspond to a PERIODIC spline
32 *      outputs :
33 *         d    the derivatives in each x(i) i = 1..n
34 *
35 *      work arrays :
36 *         A_d(1..n), A_sd(1..n-1), qdy(1..n-1)
37 *         lll(1..n-1) (used only in the periodic case)
38 *
39 *    NOTES
40 *         this routine requires (i)   n >= 3 (for natural) n >=4 (for not_a_knot)
41 *                               (ii)  strict increasing abscissae x(i)
42 *                               (iii) y(1) = y(n) in the periodic case
43 *         THESE CONDITIONS MUST BE TESTED IN THE CALLING CODE
44 *
45 *     AUTHOR
46 *        Bruno Pincon
47 *
48 *     July 22 2004 : correction of the case not_a_knot which worked only
49 *                    for equidistant abscissae
50 *
51       implicit none
52
53       integer n, type
54       double precision x(n), y(n), d(n), A_d(n), A_sd(n-1), qdy(n-1),
55      $                 lll(n-1)
56
57       include 'constinterp.h.f'
58       integer i
59       double precision r
60
61       if (n .eq. 2) then
62          if (type .ne. CLAMPED) then
63             d(1) = (y(2) - y(1)) / (x(2) - x(1))
64             d(2) = d(1)
65          endif
66          return
67       endif
68
69       if (n .eq. 3  .and.  type .eq. NOT_A_KNOT) then
70          call derivd(x, y, d, n, 1, FAST)
71          return
72       endif
73
74       do i = 1, n-1
75          A_sd(i) = 1.d0 / (x(i+1) - x(i))
76          qdy(i) = (y(i+1) - y(i)) * A_sd(i)**2
77       enddo
78
79 *     compute the coef matrix and r.h.s. for rows 2..n-1
80 *     (which don't relies on the type)
81       do i = 2, n-1
82          A_d(i) = 2.d0*( A_sd(i-1) + A_sd(i) )
83          d(i) = 3.d0 * ( qdy(i-1) + qdy(i) )
84       enddo
85
86 *     compute equ 1 and n in function of the type
87       if (type .eq. NATURAL) then
88          A_d(1) =  2.d0*A_sd(1)
89          d(1) = 3.d0 * qdy(1)
90          A_d(n) =  2.d0*A_sd(n-1)
91          d(n) =  3.d0 * qdy(n-1)
92          call TriDiagLDLtSolve(A_d, A_sd, d, n)
93
94       else if (type .eq. NOT_A_KNOT) then
95          !  s'''(x(2)-) = s'''(x(2)+)
96          r = A_sd(2)/A_sd(1)
97          A_d(1) = A_sd(1)/(1.d0 + r)
98          d(1) = ((3.d0*r+2.d0)*qdy(1)+r*qdy(2))/(1.d0+r)**2
99          !  s'''(x(n-1)-) = s'''(x(n-1)+)
100          r = A_sd(n-2)/A_sd(n-1)
101          A_d(n) = A_sd(n-1)/(1.d0 + r)
102          d(n) = ((3.d0*r+2.d0)*qdy(n-1)+r*qdy(n-2))/(1.d0+r)**2
103          call TriDiagLDLtSolve(A_d, A_sd, d, n)
104
105       else if (type .eq. CLAMPED) then
106          ! d(1) and d(n) are already known
107          d(2) = d(2) - d(1)*A_sd(1)
108          d(n-1) = d(n-1) - d(n)*A_sd(n-1)
109          call TriDiagLDLtSolve(A_d(2), A_sd(2), d(2), n-2)
110
111       else if (type .eq. PERIODIC) then
112          A_d(1) = 2.d0*( A_sd(1) + A_sd(n-1) )
113          d(1) = 3.d0 * ( qdy(1) + qdy(n-1) )
114          lll(1) = A_sd(n-1)
115          call dset(n-2, 0.d0, lll(2),1)  ! mise a zero
116          lll(n-2) = A_sd(n-2)
117          call CyclicTriDiagLDLtSolve(A_d, A_sd, lll, d, n-1)
118          d(n) = d(1)
119
120       endif
121
122       end ! subroutine SplineCub
123
124
125       subroutine TriDiagLDLtSolve(d, l, b, n)
126 *
127 *     PURPOSE
128 *        solve a linear system A x = b with a symmetric tridiagonal positive definite
129 *        matrix A by using an LDL^t factorization
130 *
131 *     PARAMETERS
132 *        d(1..n)   : on input the diagonal of A
133 *                    on output the diagonal of the (diagonal) matrix D
134 *        l(1..n-1) : on input the sub-diagonal of A
135 *                    on output the sub-diagonal of L
136 *        b(1..n)   : on input contains the r.h.s. b
137 *                    on output the solution x
138 *        n         : the dimension
139 *
140 *     CAUTION
141 *        no zero pivot detection
142 *
143       implicit none
144       integer n
145       integer flag
146       double precision d(n), l(n-1), b(n)
147
148       integer i
149       double precision temp
150
151       do i = 2, n
152          temp = l(i-1)
153          l(i-1) = l(i-1) / d(i-1)
154          d(i) = d(i) - temp * l(i-1)
155          b(i) = b(i) - l(i-1)*b(i-1)
156       enddo
157
158       b(n) = b(n) / d(n)
159       do i = n-1, 1, -1
160          b(i) = b(i)/d(i) - l(i)*b(i+1)
161       enddo
162
163       end ! subroutine TriDiagLDLtSolve
164
165       subroutine CyclicTriDiagLDLtSolve(d, lsd, lll, b, n)
166 *
167 *     PURPOSE
168 *        solve a linear system A x = b with a symmetric "nearly" tridiagonal
169 *        positive definite matrix A by using an LDL^t factorization,
170 *        the matrix A has the form :
171 *
172 *          |x x         x|                        |1            |
173 *          |x x x        |                        |x 1          |
174 *          |  x x x      |                        |  x 1        |
175 *          |    x x x    |  and so the L is like  |    x 1      |
176 *          |      x x x  |                        |      x 1    |
177 *          |        x x x|                        |        x 1  |
178 *          |x         x x|                        |x x x x x x 1|
179 *
180 *     PARAMETERS
181 *        d(1..n)     : on input the diagonal of A
182 *                      on output the diagonal of the (diagonal) matrix D
183 *        lsd(1..n-2) : on input the sub-diagonal of A (without  A(n,n-1))
184 *                      on output the sub-diagonal of L (without  L(n,n-1))
185 *        lll(1..n-1) : on input the last line of A (without A(n,n))
186 *                      on output the last line of L (without L(n,n))
187 *        b(1..n)     : on input contains the r.h.s. b
188 *                      on output the solution x
189 *        n           : the dimension
190 *
191 *     CAUTION
192 *        no zero pivot detection
193 *
194       implicit none
195       integer n
196       double precision d(n), lsd(n-2), lll(n-1), b(n)
197
198       integer i, j
199       double precision temp1, temp2
200
201 *     compute the LDL^t factorization
202       do i = 1, n-2
203          temp1 = lsd(i)
204          temp2 = lll(i)
205          lsd(i) = lsd(i) / d(i) ! elimination coef L(i,i-1)
206          lll(i) = lll(i) / d(i) ! elimination coef L(n,i-1)
207
208          d(i+1)   = d(i+1)   - lsd(i) * temp1 ! elimination on line i+1
209          lll(i+1) = lll(i+1) - lll(i) * temp1 ! elimination on line n
210          d(n)     = d(n)     - lll(i) * temp2 ! elimination on line n
211       enddo
212       temp2 = lll(n-1)
213       lll(n-1) = lll(n-1) / d(n-1)
214       d(n)     = d(n) - lll(n-1) * temp2
215
216 *     solve LDL^t x = b  (but use b for x and for the intermediary vectors...)
217       do i = 2, n-1
218          b(i) = b(i) - lsd(i-1)*b(i-1)
219       enddo
220       do j = 1, n-1
221          b(n) = b(n) - lll(j)*b(j)
222       enddo
223
224       do i = 1, n
225          b(i) = b(i) / d(i)
226       enddo
227
228       b(n-1) = b(n-1) - lll(n-1)*b(n)
229       do i = n-2, 1, -1
230          b(i) = b(i) - lsd(i)*b(i+1) - lll(i)*b(n)
231       enddo
232
233       end ! subroutine CyclicTriDiagLDLtSolve
234
235
236       integer function isearch(t, x, n)
237 *
238 *     PURPOSE
239 *        x(1..n) being an array (with strict increasing order and n >=2)
240 *        representing intervals, this routine return i such that :
241 *
242 *           x(i) <= t <= x(i+1)
243 *
244 *        and 0 if t is not in [x(1), x(n)]
245 *
246       implicit none
247       integer n
248       double precision t, x(n)
249
250       integer i1, i2, i
251
252       if ( x(1) .le. t  .and.  t .le. x(n)) then
253 *        dichotomic search
254          i1 = 1
255          i2 = n
256          do while ( i2 - i1 .gt. 1 )
257             i = (i1 + i2)/2
258             if ( t .le. x(i) ) then
259                i2 = i
260             else
261                i1 = i
262             endif
263          enddo
264          isearch = i1
265       else
266          isearch = 0
267       endif
268
269       end
270
271       subroutine EvalPWHermite(t, st, dst, d2st, d3st, m, x, y, d, n,
272      $                         outmode)
273 *
274 *     PURPOSE
275 *        evaluation at the abscissae t(1..m) of the piecewise hermite function
276 *        define by x(1..n), y(1..n), d(1..n) (d being the derivatives at the
277 *        x(i)) together with its derivative, second derivative and third derivative
278 *
279 *     PARAMETERS
280 *
281 *        outmode : define what return in case t(j) not in [x(1), x(n)]
282
283       implicit none
284       integer m, n, outmode
285       double precision t(m), st(m), dst(m), d2st(m), d3st(m),
286      $                 x(n), y(n), d(n)
287
288       include 'constinterp.h.f'
289       integer i, j
290       integer  isearch, isanan
291       external isearch, isanan
292       double precision tt
293       external         returnananfortran
294       logical new_call
295       common /INFO_HERMITE/new_call
296
297       new_call = .true.
298       i = 0
299       do j = 1, m
300          tt = t(j)
301          call fast_int_search(tt, x, n, i)  ! recompute i only if necessary
302
303          if ( i .ne. 0 ) then
304             call EvalHermite(tt, x(i), x(i+1), y(i), y(i+1), d(i),
305      $                       d(i+1), st(j), dst(j), d2st(j), d3st(j), i)
306
307          else   ! t(j) is outside [x(1), x(n)] evaluation depend upon outmode
308
309             if (outmode .eq. BY_NAN  .or.  isanan(tt) .eq. 1) then
310                CALL returnananfortran(st(j))
311                dst(j) = st(j)
312                d2st(j) = st(j)
313                d3st(j) = st(j)
314
315             elseif (outmode .eq. BY_ZERO) then
316                st(j) = 0.d0
317                dst(j) = 0.d0
318                d2st(j) = 0.d0
319                d3st(j) = 0.d0
320
321             elseif (outmode .eq. C0) then
322                dst(j) = 0.d0
323                d2st(j) = 0.d0
324                d3st(j) = 0.d0
325                if ( tt .lt. x(1) ) then
326                   st(j) = y(1)
327                else
328                   st(j) = y(n)
329                endif
330
331             elseif (outmode .eq. LINEAR) then
332                d2st(j) = 0.d0
333                d3st(j) = 0.d0
334                if ( tt .lt. x(1) ) then
335                   dst(j) = d(1)
336                   st(j) = y(1) + (tt - x(1))*d(1)
337                else
338                   dst(j) = d(n)
339                   st(j) = y(n) + (tt - x(n))*d(n)
340                endif
341
342             else
343                if (outmode .eq. NATURAL) then
344                   call near_interval(tt, x, n, i)
345                elseif (outmode .eq. PERIODIC) then
346                   call coord_by_periodicity(tt, x, n, i)
347                endif
348                call EvalHermite(tt, x(i), x(i+1), y(i), y(i+1), d(i),
349      $                       d(i+1), st(j), dst(j), d2st(j), d3st(j), i)
350             endif
351          endif
352       enddo
353
354       end ! subroutine EvalPWHermite
355
356       subroutine EvalHermite(t, xa, xb, ya, yb, da, db, h, dh, ddh,
357      $                       dddh, i)
358
359       implicit none
360       double precision t, xa, xb, ya, yb, da, db, h, dh, ddh, dddh
361       integer i
362
363       logical new_call
364       common /INFO_HERMITE/new_call
365
366       double precision tmxa, dx, p, c2, c3
367       integer          old_i
368       save             old_i,       c2, c3
369       data old_i/0/
370
371       if (old_i .ne. i  .or.  new_call) then
372 *        compute the following Newton form :
373 *           h(t) = ya + da*(t-xa) + c2*(t-xa)^2 + c3*(t-xa)^2*(t-xb)
374          dx = 1.d0/(xb - xa)
375          p = (yb - ya) * dx
376          c2 = (p - da) * dx
377          c3 = ( (db - p) + (da - p) ) * (dx*dx)
378          new_call = .false.
379       endif
380       old_i = i
381
382 *     eval h(t), h'(t), h"(t) and h"'(t), by a generalized Horner 's scheme
383       tmxa = t - xa
384       h = c2 + c3*(t - xb)
385       dh = h + c3*tmxa
386       ddh = 2.d0*( dh + c3*tmxa )
387       dddh = 6.d0*c3
388       h = da + h*tmxa
389       dh = h + dh*tmxa
390       h = ya + h*tmxa
391
392       end ! subroutine EvalHermite
393
394       subroutine fast_int_search(xx, x, nx, i)
395
396       implicit none
397       integer nx, i
398       double precision xx, x(nx)
399
400       integer  isearch
401       external isearch
402
403       if ( i .eq. 0 ) then
404          i = isearch(xx, x, nx)
405       elseif ( .not. (x(i) .le. xx  .and.  xx .le. x(i+1))) then
406          i = isearch(xx, x, nx)
407       endif
408
409       end
410
411       subroutine coord_by_periodicity(t, x, n, i)
412 *
413 *     PURPOSE
414 *        recompute t such that t in [x(1), x(n)] by periodicity :
415 *        and then the interval i of this new t
416 *
417       implicit none
418       integer n, i
419       double precision t, x(n)
420       integer  isearch
421       external isearch
422
423       double precision r, dx
424
425       dx = x(n) - x(1)
426       r = (t - x(1)) / dx
427
428       if (r .ge. 0.d0) then
429          t = x(1) + (r - aint(r))*dx
430       else
431          r = abs(r)
432          t = x(n) - (r - aint(r))*dx
433       endif
434
435       ! some cautions in case of roundoff errors (is necessary ?)
436       if (t .lt. x(1)) then
437          t = x(1)
438          i = 1
439       else if (t .gt. x(n)) then
440          t = x(n)
441          i  = n-1
442       else
443          i = isearch(t, x, n)
444       endif
445
446       end ! subroutine coord_by_periodicity
447
448
449       subroutine near_grid_point(xx, x, nx, i)
450 *
451 *     calcule le point de la grille le plus proche ... a detailler
452 *
453       implicit none
454       integer i, nx
455       double precision xx, x(nx)
456
457       if (xx .lt. x(1)) then
458          i = 1
459          xx = x(1)
460       else    !  xx > x(nx)
461          i = nx-1
462          xx = x(nx)
463       endif
464
465       end
466
467       subroutine near_interval(xx, x, nx, i)
468 *
469 *     idem sans modifier xx
470 *
471       implicit none
472       integer i, nx
473       double precision xx, x(nx)
474
475       if (xx .lt. x(1)) then
476          i = 1
477       else    !  xx > x(nx)
478          i = nx-1
479       endif
480
481       end
482
483       subroutine proj_by_per(t, xmin, xmax)
484 *
485 *     PURPOSE
486 *        recompute t such that t in [xmin, xmax] by periodicity.
487 *
488       implicit none
489       double precision t, xmin, xmax
490       double precision r, dx
491
492       dx = xmax - xmin
493       r = (t - xmin) / dx
494
495       if (r .ge. 0.d0) then
496          t = xmin + (r - aint(r))*dx
497       else
498          r = abs(r)
499          t = xmax - (r - aint(r))*dx
500       endif
501
502       ! some cautions in case of roundoff errors (is necessary ?)
503       if (t .lt. xmin) then
504          t = xmin
505       else if (t .gt. xmax) then
506          t = xmax
507       endif
508
509       end ! subroutine proj_by_per
510
511
512       subroutine proj_on_grid(xx, xmin, xmax)
513 *
514       implicit none
515       double precision xx, xmin, xmax
516
517       if (xx .lt. xmin) then
518          xx = xmin
519       else
520          xx = xmax
521       endif
522
523       end
524
525       subroutine BiCubicSubSpline(x, y, u, nx, ny, C, p, q, r, type)
526 *
527 *     PURPOSE
528 *        compute bicubic subsplines
529 *
530       implicit none
531       integer nx, ny, type
532       double precision x(nx), y(ny), u(nx, ny), C(4,4,nx-1,ny-1),
533      $                 p(nx, ny), q(nx, ny), r(nx, ny)
534       integer i, j
535       include 'constinterp.h.f'
536
537       if (type .eq. MONOTONE) then
538 *        approximation des derivees par SUBROUTINE DPCHIM(N,X,F,D,INCFD)
539          ! p = du/dx
540          do j = 1, ny
541             call dpchim(nx, x, u(1,j), p(1,j), 1)
542          enddo
543          ! q = du/dy
544          do i = 1, nx
545             call dpchim(ny, y, u(i,1), q(i,1), nx)
546          enddo
547          ! r = d2 u/ dx dy  approchee via  dq / dx
548          do j = 1, ny
549             call dpchim(nx, x, q(1,j), r(1,j), 1)
550          enddo
551
552       elseif (type .eq. FAST  .or.  type .eq. FAST_PERIODIC) then
553 *        approximation des derivees partielles par methode simple
554
555          ! p = du/dx
556          do j = 1, ny
557             call derivd(x, u(1,j), p(1,j), nx, 1, type)
558          enddo
559          ! q = du/dy
560          do i = 1, nx
561             call derivd(y, u(i,1), q(i,1), ny, nx, type)
562          enddo
563          ! r = d2 u/ dx dy  approchee via  dq / dx
564          do j = 1, ny
565             call derivd(x, q(1,j), r(1,j), nx, 1, type)
566          enddo
567
568       endif
569
570 *     calculs des coefficients dans les bases (x-x(i))^k (y-y(j))^l  0<= k,l <= 3
571 *     pour evaluation rapide via Horner par la suite
572       call coef_bicubic(u, p, q, r, x, y, nx, ny, C)
573
574       end ! subroutine BiCubicSubSpline
575
576       subroutine BiCubicSpline(x, y, u, nx, ny, C, p, q, r,
577      $                         A_d, A_sd, d, ll, qdu, u_temp, type)
578 *
579 *     PURPOSE
580 *        compute bicubic splines
581 *
582       implicit none
583       integer nx, ny, type
584       double precision x(nx), y(ny), u(nx, ny), C(4,4,nx-1,ny-1),
585      $                 p(nx, ny), q(nx, ny), r(nx, ny), A_d(*),
586      $                 A_sd(*), d(ny), ll(*), qdu(*), u_temp(ny)
587       include 'constinterp.h.f'
588       integer i, j
589
590       ! compute du/dx
591       do j = 1, ny
592          call SplineCub(x, u(1,j), p(1,j), nx, type, A_d, A_sd, qdu, ll)
593       enddo
594
595       ! compute du/dy
596       do i = 1, nx
597          call dcopy(ny, u(i,1), nx, u_temp, 1)
598          call SplineCub(y, u_temp, d, ny, type, A_d, A_sd, qdu, ll)
599          call dcopy(ny, d, 1, q(i,1), nx)
600       enddo
601
602       ! compute ddu/dxdy
603       call SplineCub(x, q(1,1),  r(1,1),  nx, type, A_d, A_sd, qdu, ll)
604       call SplineCub(x, q(1,ny), r(1,ny), nx, type, A_d, A_sd, qdu, ll)
605
606       do i = 1, nx
607          call dcopy(ny, p(i,1), nx, u_temp, 1)
608          d(1) = r(i,1)
609          d(ny) = r(i,ny)
610          call SplineCub(y, u_temp, d, ny, CLAMPED, A_d, A_sd, qdu, ll)
611          call dcopy(ny-2, d(2), 1, r(i,2), nx)
612       enddo
613
614 *     calculs des coefficients dans les bases (x-x(i))^k (y-y(j))^l  0<= k,l <= 3
615 *     pour evaluation rapide via Horner par la suite
616       call coef_bicubic(u, p, q, r, x, y, nx, ny, C)
617
618       end  ! subroutine BiCubicSpline
619
620
621       subroutine derivd(x, u, du, n, inc, type)
622 *
623 *     PURPOSE
624 *        given functions values u(i) at points x(i),  i = 1, ..., n
625 *        this subroutine computes approximations du(i) of the derivative
626 *        at the points x(i).
627 *
628 *     METHOD
629 *        For i in [2,n-1], the "centered" formula of order 2 is used :
630 *            d(i) = derivative at x(i) of the interpolation polynomial
631 *                   of the points {(x(j),u(j)), j in [i-1,i+1]}
632 *
633 *         For i=1 and n, if type = FAST_PERIODIC  (in which case u(n)=u(1)) then
634 *         the previus "centered" formula is also used else (type = FAST), d(1)
635 *         is the derivative at x(1) of the interpolation polynomial of
636 *         {(x(j),u(j)), j in [1,3]} and the same method is used for d(n)
637 *
638 *     ARGUMENTS
639 *      inputs :
640 *         n       integer : number of point (n >= 2)
641 *         x, u    double precision : the n points, x must be in strict increasing order
642 *         type    integer : FAST (the function is non periodic) or FAST_PERIODIC
643 *                 (the function is periodic), in this last case u(n) must be equal to u(1))
644 *         inc     integer : to deal easily with 2d applications, u(i) is in fact
645 *                 u(1,i) with u declared as u(inc,*) to avoid the direct management of
646 *                 the increment inc (the i th value given with u(1 + inc*(i-1) ...)
647 *      outputs :
648 *         d       the derivatives in each x(i) i = 1..n
649 *
650 *    NOTES
651 *         this routine requires (i)   n >= 2
652 *                               (ii)  strict increasing abscissae x(i)
653 *                               (iii) u(1)=u(n) if type = FAST_PERIODIC
654 *         ALL THESE CONDITIONS MUST BE TESTED IN THE CALLING CODE
655 *
656 *     AUTHOR
657 *        Bruno Pincon
658 *
659
660       implicit none
661       integer n, inc, type
662       double precision x(n), u(inc,*), du(inc,*)
663
664       include 'constinterp.h.f'
665       double precision dx_l, du_l, dx_r, du_r, w_l, w_r
666       integer i, k
667
668
669       if (n .eq. 2) then  ! special case used linear interp
670          du(1,1) = (u(1,2) - u(1,1))/(x(2)-x(1))
671          du(1,2) = du(1,1)
672          return
673       endif
674
675       if (type .eq. FAST_PERIODIC) then
676
677          dx_r  = x(n)-x(n-1)
678          du_r  = (u(1,1) - u(1,n-1)) / dx_r
679          do i = 1, n-1
680             dx_l = dx_r
681             du_l = du_r
682             dx_r = x(i+1) - x(i)
683             du_r = (u(1,i+1) - u(1,i))/dx_r
684             w_l  = dx_r/(dx_l + dx_r)
685             w_r = 1.d0 - w_l
686             du(1,i) = w_l*du_l + w_r*du_r
687          enddo
688          du(1,n) = du(1,1)
689
690       else if (type .eq. FAST) then
691
692          dx_l = x(2) - x(1)
693          du_l = (u(1,2) - u(1,1))/dx_l
694          dx_r = x(3) - x(2)
695          du_r = (u(1,3) - u(1,2))/dx_r
696          w_l = dx_r/(dx_l + dx_r)
697          w_r = 1.d0 - w_l
698          du(1,1) = (1.d0 + w_r)*du_l - w_r*du_r
699          du(1,2) = w_l*du_l + w_r*du_r
700          do i = 3, n-1
701             dx_l = dx_r
702             du_l  = du_r
703             dx_r = x(i+1) - x(i)
704             du_r  = (u(1,i+1) - u(1,i))/dx_r
705             w_l  = dx_r/(dx_l + dx_r)
706             w_r = 1.d0 - w_l
707             du(1,i) = w_l*du_l + w_r*du_r
708          enddo
709          du(1,n) = (1.d0 + w_l)*du_r - w_l*du_l
710
711       endif
712
713       end
714
715
716       subroutine coef_bicubic(u, p, q, r, x, y, nx, ny, C)
717 *
718 *     PURPOSE
719 *        compute for each polynomial (i,j)-patch (defined on
720 *        [x(i),x(i+1)]x[y(i),y(i+1)]) the following base
721 *        representation :
722 *           i,j        _4_  _4_   i,j       k-1       l-1
723 *          u   (x,y) = >__  >__  C   (x-x(i))  (y-y(j))
724 *                      k=1  l=1   k,l
725 *
726 *        from the "Hermite" representation (values of u, p = du/dx,
727 *        q = du/dy, r = ddu/dxdy at the 4 vertices (x(i),y(j)),
728 *        (x(i+1),y(j)), (x(i+1),y(j+1)), (x(i),y(j+1)).
729 *
730       implicit none
731
732       integer nx, ny
733       double precision u(nx, ny), p(nx, ny), q(nx, ny), r(nx, ny),
734      $                 x(nx), y(ny), C(4,4,nx-1,ny-1)
735
736       integer i, j
737       double precision a, b, cc, d, dx, dy
738
739       do j = 1, ny-1
740          dy = 1.d0/(y(j+1) - y(j))
741
742          do i = 1, nx-1
743             dx = 1.d0/(x(i+1) - x(i))
744
745             C(1,1,i,j) = u(i,j)
746             C(2,1,i,j) = p(i,j)
747             C(1,2,i,j) = q(i,j)
748             C(2,2,i,j) = r(i,j)
749
750             a = (u(i+1,j) - u(i,j))*dx
751             C(3,1,i,j) = (3.d0*a -2.d0*p(i,j) - p(i+1,j))*dx
752             C(4,1,i,j) = (p(i+1,j) + p(i,j) - 2.d0*a)*(dx*dx)
753
754             a = (u(i,j+1) - u(i,j))*dy
755             C(1,3,i,j) = (3.d0*a -2.d0*q(i,j) - q(i,j+1))*dy
756             C(1,4,i,j) = (q(i,j+1) + q(i,j) - 2.d0*a)*(dy*dy)
757
758             a = (q(i+1,j) - q(i,j))*dx
759             C(3,2,i,j) = (3.d0*a - r(i+1,j) - 2.d0*r(i,j))*dx
760             C(4,2,i,j) = (r(i+1,j) + r(i,j) - 2.d0*a)*(dx*dx)
761
762             a = (p(i,j+1) - p(i,j))*dy
763             C(2,3,i,j) = (3.d0*a - r(i,j+1) - 2.d0*r(i,j))*dy
764             C(2,4,i,j) = (r(i,j+1) + r(i,j) - 2.d0*a)*(dy*dy)
765
766             a = (u(i+1,j+1)+u(i,j)-u(i+1,j)-u(i,j+1))*(dx*dx*dy*dy)
767      $          - (p(i,j+1)-p(i,j))*(dx*dy*dy)
768      $          - (q(i+1,j)-q(i,j))*(dx*dx*dy)
769      $          + r(i,j)*(dx*dy)
770             b = (p(i+1,j+1)+p(i,j)-p(i+1,j)-p(i,j+1))*(dx*dy*dy)
771      $          - (r(i+1,j)-r(i,j))*(dx*dy)
772             cc = (q(i+1,j+1)+q(i,j)-q(i+1,j)-q(i,j+1))*(dx*dx*dy)
773      $          - (r(i,j+1)-r(i,j))*(dx*dy)
774             d = (r(i+1,j+1)+r(i,j)-r(i+1,j)-r(i,j+1))*(dx*dy)
775
776
777             C(3,3,i,j) =  9.d0*a - 3.d0*b - 3.d0*cc + d
778             C(3,4,i,j) =(-6.d0*a + 2.d0*b + 3.d0*cc - d)*dy
779             C(4,3,i,j) =(-6.d0*a + 3.d0*b + 2.d0*cc - d)*dx
780             C(4,4,i,j) =( 4.d0*a - 2.d0*b - 2.d0*cc + d)*dx*dy
781
782          enddo
783       enddo
784       end
785
786       double precision function EvalBicubic(xx, yy, xk, yk, Ck)
787
788       implicit none
789       double precision xx, yy, xk, yk, Ck(4,4)
790
791       double precision dx, dy, u
792       integer i
793
794       dx = xx - xk
795       dy = yy - yk
796
797       u = 0.d0
798       do i = 4, 1, -1
799          u = Ck(i,1) + dy*(Ck(i,2) + dy*(Ck(i,3) + dy*Ck(i,4))) + u*dx
800       enddo
801       EvalBicubic = u
802
803       end ! function EvalBicubic
804
805       subroutine EvalBicubic_with_grad(xx, yy, xk, yk, Ck,
806      $                                 u, dudx, dudy)
807
808       implicit none
809       double precision xx, yy, xk, yk, Ck(4,4), u, dudx, dudy
810
811       double precision dx, dy
812       integer i
813
814       dx = xx - xk
815       dy = yy - yk
816       u = 0.d0
817       dudx = 0.d0
818       dudy = 0.d0
819       do i = 4, 1, -1
820          u = Ck(i,1) + dy*(Ck(i,2) + dy*(Ck(i,3) + dy*Ck(i,4))) + u*dx
821          dudx = Ck(2,i)+dx*(2.d0*Ck(3,i) + dx*3.d0*Ck(4,i)) + dudx*dy
822          dudy = Ck(i,2)+dy*(2.d0*Ck(i,3) + dy*3.d0*Ck(i,4)) + dudy*dx
823       enddo
824
825       end ! subroutine EvalBicubic_with_grad
826
827
828       subroutine EvalBicubic_with_grad_and_hes(xx, yy, xk, yk, Ck,
829      $                                         u, dudx, dudy,
830      $                                         d2udx2, d2udxy, d2udy2)
831
832       implicit none
833       double precision xx, yy, xk, yk, Ck(4,4), u, dudx, dudy, d2udx2,
834      $                 d2udxy, d2udy2
835
836       double precision dx, dy
837       integer i
838
839       dx = xx - xk
840       dy = yy - yk
841       u = 0.d0
842       dudx = 0.d0
843       dudy = 0.d0
844       d2udx2 = 0.d0
845       d2udy2 = 0.d0
846       d2udxy = 0.d0
847       do i = 4, 1, -1
848          u = Ck(i,1) + dy*(Ck(i,2) + dy*(Ck(i,3) + dy*Ck(i,4))) + u*dx
849          dudx = Ck(2,i)+dx*(2.d0*Ck(3,i) + dx*3.d0*Ck(4,i)) + dudx*dy
850          dudy = Ck(i,2)+dy*(2.d0*Ck(i,3) + dy*3.d0*Ck(i,4)) + dudy*dx
851          d2udx2 = 2.d0*Ck(3,i) + dx*6.d0*Ck(4,i) + d2udx2*dy
852          d2udy2 = 2.d0*Ck(i,3) + dy*6.d0*Ck(i,4) + d2udy2*dx
853       enddo
854       d2udxy =            Ck(2,2)+dy*(2.d0*Ck(2,3)+dy*3.d0*Ck(2,4))
855      .        + dx*(2.d0*(Ck(3,2)+dy*(2.d0*Ck(3,3)+dy*3.d0*Ck(3,4)))
856      .        + dx*(3.d0*(Ck(4,2)+dy*(2.d0*Ck(4,3)+dy*3.d0*Ck(4,4)))))
857       end ! subroutine EvalBicubic_with_grad_and_hes
858
859
860       subroutine BiCubicInterp(x, y, C, nx, ny, x_eval, y_eval, z_eval,
861      $                         m, outmode)
862 *
863 *     PURPOSE
864 *        bicubic interpolation :
865 *          the grid is defined by x(1..nx), y(1..ny)
866 *          the known values are z(1..nx,1..ny), (z(i,j) being the value
867 *          at point (x(i),y(j)))
868 *          the interpolation is done on the points x_eval,y_eval(1..m)
869 *          z_eval(k) is the result of the bicubic interpolation of
870 *          (x_eval(k), y_eval(k))
871 *
872       implicit none
873       integer nx, ny, m, outmode
874       double precision x(nx), y(ny), C(4,4,nx-1,ny-1),
875      $                 x_eval(m), y_eval(m), z_eval(m)
876
877       double precision xx, yy
878       integer i, j, k
879       include 'constinterp.h.f'
880       integer  isanan
881       double precision EvalBicubic
882       external isanan, returnananfortran, EvalBicubic
883
884       i = 0
885       j = 0
886       do k = 1, m
887          xx = x_eval(k)
888          call fast_int_search(xx, x, nx, i)
889          yy = y_eval(k)
890          call fast_int_search(yy, y, ny, j)
891
892          if (i .ne. 0  .and.  j .ne. 0) then
893             z_eval(k) = EvalBicubic(xx, yy, x(i), y(j), C(1,1,i,j))
894
895          elseif (outmode .eq. BY_NAN  .or.  isanan(xx) .eq. 1
896      $                                .or.  isanan(yy) .eq. 1) then
897             CALL returnananfortran(z_eval(k))
898
899          elseif (outmode .eq. BY_ZERO) then
900             z_eval(k) = 0.d0
901
902          elseif (outmode .eq. PERIODIC) then
903             if (i .eq. 0) call coord_by_periodicity(xx, x, nx, i)
904             if (j .eq. 0) call coord_by_periodicity(yy, y, ny, j)
905             z_eval(k) = EvalBicubic(xx, yy, x(i), y(j), C(1,1,i,j))
906
907          elseif (outmode .eq. C0) then
908             if (i .eq. 0) call near_grid_point(xx, x, nx, i)
909             if (j .eq. 0) call near_grid_point(yy, y, ny, j)
910             z_eval(k) = EvalBicubic(xx, yy, x(i), y(j), C(1,1,i,j))
911
912          elseif (outmode .eq. NATURAL) then
913             if (i .eq. 0) call near_interval(xx, x, nx, i)
914             if (j .eq. 0) call near_interval(yy, y, ny, j)
915             z_eval(k) = EvalBicubic(xx, yy, x(i), y(j), C(1,1,i,j))
916          endif
917
918       enddo
919
920       end
921
922       subroutine BiCubicInterpWithGrad(x, y, C, nx, ny, x_eval, y_eval,
923      $                                 z_eval, dzdx_eval, dzdy_eval,m,
924      $                                 outmode)
925 *
926 *     PURPOSE
927 *        bicubic interpolation :
928 *          the grid is defined by x(1..nx), y(1..ny)
929 *          the known values are z(1..nx,1..ny), (z(i,j) being the value
930 *          at point (x(i),y(j)))
931 *          the interpolation is done on the points x_eval,y_eval(1..m)
932 *          z_eval(k) is the result of the bicubic interpolation of
933 *          (x_eval(k), y_eval(k)) and dzdx_eval(k), dzdy_eval(k) is the gradient
934 *
935       implicit none
936       integer nx, ny, m, outmode
937       double precision x(nx), y(ny), C(4,4,nx-1,ny-1),
938      $                 x_eval(m), y_eval(m), z_eval(m),
939      $                 dzdx_eval(m), dzdy_eval(m)
940
941       double precision xx, yy
942       integer k, i, j
943       include 'constinterp.h.f'
944       integer  isanan
945       double precision EvalBicubic
946       external isanan, returnananfortran, EvalBicubic
947       logical change_dzdx, change_dzdy
948
949       i = 0
950       j = 0
951       do k = 1, m
952          xx = x_eval(k)
953          call fast_int_search(xx, x, nx, i)
954          yy = y_eval(k)
955          call fast_int_search(yy, y, ny, j)
956
957          if (i .ne. 0  .and.  j .ne. 0) then
958             call Evalbicubic_with_grad(xx, yy, x(i), y(j), C(1,1,i,j),
959      $                                 z_eval(k),
960      $                                 dzdx_eval(k), dzdy_eval(k))
961
962          elseif ( outmode .eq. BY_NAN  .or.  isanan(xx) .eq. 1
963      $                                 .or.  isanan(yy) .eq. 1) then
964             CALL returnananfortran(z_eval(k))
965             dzdx_eval(k) = z_eval(k)
966             dzdy_eval(k) = z_eval(k)
967
968          elseif (outmode .eq. BY_ZERO ) then
969             z_eval(k) = 0.d0
970             dzdx_eval(k) = 0.d0
971             dzdy_eval(k) = 0.d0
972
973          elseif (outmode .eq. PERIODIC) then
974             if (i .eq. 0) call coord_by_periodicity(xx, x, nx, i)
975             if (j .eq. 0) call coord_by_periodicity(yy, y, ny, j)
976             call Evalbicubic_with_grad(xx, yy, x(i), y(j), C(1,1,i,j),
977      $                                 z_eval(k),
978      $                                 dzdx_eval(k), dzdy_eval(k))
979
980          elseif (outmode .eq. C0) then
981             if (i .eq. 0) then
982                call near_grid_point(xx, x, nx, i)
983                change_dzdx = .true.
984             else
985                change_dzdx = .false.
986             endif
987             if (j .eq. 0) then
988                call near_grid_point(yy, y, ny, j)
989                change_dzdy = .true.
990             else
991                change_dzdy = .false.
992             endif
993             call Evalbicubic_with_grad(xx, yy, x(i), y(j), C(1,1,i,j),
994      $                                 z_eval(k),
995      $                                 dzdx_eval(k), dzdy_eval(k))
996             if (change_dzdx) dzdx_eval(k) = 0.d0
997             if (change_dzdy) dzdy_eval(k) = 0.d0
998
999          elseif (outmode .eq. NATURAL) then
1000             if (i .eq. 0) call near_interval(xx, x, nx, i)
1001             if (j .eq. 0) call near_interval(yy, y, ny, j)
1002             call Evalbicubic_with_grad(xx, yy, x(i), y(j), C(1,1,i,j),
1003      $                                 z_eval(k),
1004      $                                 dzdx_eval(k), dzdy_eval(k))
1005
1006          endif
1007
1008       enddo
1009
1010       end
1011
1012       subroutine BiCubicInterpWithGradAndHes(x, y, C, nx, ny,
1013      $                                       x_eval, y_eval, z_eval,
1014      $                                         dzdx_eval, dzdy_eval,
1015      $                        d2zdx2_eval, d2zdxy_eval, d2zdy2_eval,
1016      $                                                    m, outmode)
1017 *
1018 *     PURPOSE
1019 *        bicubic interpolation :
1020 *          the grid is defined by x(1..nx), y(1..ny)
1021 *          the known values are z(1..nx,1..ny), (z(i,j) being the value
1022 *          at point (x(i),y(j)))
1023 *          the interpolation is done on the points x_eval,y_eval(1..m)
1024 *          z_eval(k) is the result of the bicubic interpolation of
1025 *          (x_eval(k), y_eval(k)), [dzdx_eval(k), dzdy_eval(k)] is the gradient
1026 *          and [d2zdx2(k), d2zdxy(k), d2zdy2(k)] the Hessian
1027 *
1028       implicit none
1029       integer nx, ny, m, outmode
1030       double precision x(nx), y(ny), C(4,4,nx-1,ny-1),
1031      $                 x_eval(m), y_eval(m), z_eval(m), dzdx_eval(m),
1032      $   dzdy_eval(m), d2zdx2_eval(m), d2zdxy_eval(m), d2zdy2_eval(m)
1033
1034       double precision xx, yy
1035       integer k, i, j
1036       include 'constinterp.h.f'
1037       integer  isanan
1038       double precision EvalBicubic
1039       external isanan, returnananfortran, EvalBicubic
1040       logical change_dzdx, change_dzdy
1041
1042       i = 0
1043       j = 0
1044       do k = 1, m
1045          xx = x_eval(k)
1046          call fast_int_search(xx, x, nx, i)
1047          yy = y_eval(k)
1048          call fast_int_search(yy, y, ny, j)
1049
1050          if (i .ne. 0  .and.  j .ne. 0) then
1051             call Evalbicubic_with_grad_and_hes(xx, yy, x(i), y(j),
1052      $                                      C(1,1,i,j), z_eval(k),
1053      $                                 dzdx_eval(k), dzdy_eval(k),
1054      $             d2zdx2_eval(k), d2zdxy_eval(k), d2zdy2_eval(k))
1055
1056          elseif ( outmode .eq. BY_NAN  .or.  isanan(xx) .eq. 1
1057      $                                 .or.  isanan(yy) .eq. 1) then
1058             CALL returnananfortran(z_eval(k))
1059             dzdx_eval(k) = z_eval(k)
1060             dzdy_eval(k) = z_eval(k)
1061             d2zdx2_eval(k) = z_eval(k)
1062             d2zdxy_eval(k) = z_eval(k)
1063             d2zdy2_eval(k) = z_eval(k)
1064
1065          elseif (outmode .eq. BY_ZERO ) then
1066             z_eval(k) = 0.d0
1067             dzdx_eval(k) = 0.d0
1068             dzdy_eval(k) = 0.d0
1069             d2zdx2_eval(k) = 0.d0
1070             d2zdxy_eval(k) = 0.d0
1071             d2zdy2_eval(k) = 0.d0
1072
1073          elseif (outmode .eq. PERIODIC) then
1074             if (i .eq. 0) call coord_by_periodicity(xx, x, nx, i)
1075             if (j .eq. 0) call coord_by_periodicity(yy, y, ny, j)
1076             call Evalbicubic_with_grad_and_hes(xx, yy, x(i), y(j),
1077      $                                      C(1,1,i,j), z_eval(k),
1078      $                                 dzdx_eval(k), dzdy_eval(k),
1079      $             d2zdx2_eval(k), d2zdxy_eval(k), d2zdy2_eval(k))
1080
1081          elseif (outmode .eq. C0) then
1082             if (i .eq. 0) then
1083                call near_grid_point(xx, x, nx, i)
1084                change_dzdx = .true.
1085             else
1086                change_dzdx = .false.
1087             endif
1088             if (j .eq. 0) then
1089                call near_grid_point(yy, y, ny, j)
1090                change_dzdy = .true.
1091             else
1092                change_dzdy = .false.
1093             endif
1094             call Evalbicubic_with_grad_and_hes(xx, yy, x(i), y(j),
1095      $                                      C(1,1,i,j), z_eval(k),
1096      $                                 dzdx_eval(k), dzdy_eval(k),
1097      $             d2zdx2_eval(k), d2zdxy_eval(k), d2zdy2_eval(k))
1098             if (change_dzdx) then
1099                dzdx_eval(k) = 0.d0
1100                d2zdx2_eval(k) = 0.d0
1101                d2zdxy_eval(k) = 0.d0
1102             endif
1103             if (change_dzdy) then
1104                dzdy_eval(k) = 0.d0
1105                d2zdxy_eval(k) = 0.d0
1106                d2zdy2_eval(k) = 0.d0
1107             endif
1108
1109          elseif (outmode .eq. NATURAL) then
1110             if (i .eq. 0) call near_interval(xx, x, nx, i)
1111             if (j .eq. 0) call near_interval(yy, y, ny, j)
1112             call Evalbicubic_with_grad_and_hes(xx, yy, x(i), y(j),
1113      $                                      C(1,1,i,j), z_eval(k),
1114      $                                 dzdx_eval(k), dzdy_eval(k),
1115      $             d2zdx2_eval(k), d2zdxy_eval(k), d2zdy2_eval(k))
1116          endif
1117
1118       enddo
1119
1120       end
1121
1122       subroutine driverdb3val(xp, yp, zp, fp, np, tx, ty, tz,
1123      $                        nx, ny, nz, kx, ky, kz, bcoef, work,
1124      $                        xmin, xmax, ymin, ymax, zmin, zmax,
1125      $                        outmode)
1126 *
1127 *     PURPOSE
1128 *        driver on to db3val
1129 *
1130       implicit none
1131       integer np, nx, ny, nz, kx, ky, kz, outmode
1132       double precision xp(np), yp(np), zp(np), fp(np),
1133      $                 tx(*), ty(*), tz(*), bcoef(*), work(*),
1134      $                 xmin, xmax, ymin, ymax, zmin, zmax
1135
1136       integer k
1137       logical flag_x, flag_y, flag_z
1138       double precision x, y, z
1139       include 'constinterp.h.f'
1140       integer  isanan
1141       double precision db3val
1142       external isanan, returnananfortran, db3val
1143
1144       do k = 1, np
1145          x = xp(k)
1146          if ( xmin .le. x .and. x .le. xmax ) then
1147             flag_x = .true.
1148          else
1149             flag_x = .false.
1150          endif
1151          y = yp(k)
1152          if ( ymin .le. y .and. y .le. ymax ) then
1153             flag_y = .true.
1154          else
1155             flag_y = .false.
1156          endif
1157          z = zp(k)
1158          if ( zmin .le. z .and. z .le. zmax ) then
1159             flag_z = .true.
1160          else
1161             flag_z = .false.
1162          endif
1163
1164          if ( flag_x .and. flag_y .and. flag_z ) then
1165             fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz,
1166      $                     nx, ny, nz, kx, ky, kz, bcoef, work)
1167
1168          elseif (outmode .eq. BY_NAN  .or.  isanan(x) .eq. 1
1169      $                                .or.  isanan(y) .eq. 1
1170      $                                .or.  isanan(z) .eq. 1) then
1171             CALL returnananfortran(fp(k))
1172
1173          elseif (outmode .eq. BY_ZERO) then
1174             fp(k) = 0.d0
1175
1176          else
1177             if (outmode .eq. PERIODIC) then
1178                if (.not. flag_x) call proj_by_per(x, xmin, xmax)
1179                if (.not. flag_y) call proj_by_per(y, ymin, ymax)
1180                if (.not. flag_z) call proj_by_per(z, zmin, zmax)
1181             elseif (outmode .eq. C0) then
1182                if (.not. flag_x) call proj_on_grid(x, xmin, xmax)
1183                if (.not. flag_y) call proj_on_grid(y, ymin, ymax)
1184                if (.not. flag_z) call proj_on_grid(z, zmin, zmax)
1185             endif
1186             fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz, nx, ny, nz,
1187      $                     kx, ky, kz, bcoef, work)
1188          endif
1189       enddo
1190
1191       end
1192
1193       subroutine driverdb3valwithgrad(xp, yp, zp, fp,
1194      $                        dfpdx, dfpdy, dfpdz, np, tx, ty, tz,
1195      $                        nx, ny, nz, kx, ky, kz, bcoef, work,
1196      $                        xmin, xmax, ymin, ymax, zmin, zmax,
1197      $                        outmode)
1198 *
1199 *     PURPOSE
1200 *        driver on to db3val with gradient computing
1201 *
1202       implicit none
1203       integer np, nx, ny, nz, kx, ky, kz, outmode
1204       double precision xp(np), yp(np), zp(np), fp(np),
1205      $                 dfpdx(np), dfpdy(np), dfpdz(np),
1206      $                 tx(*), ty(*), tz(*), bcoef(*), work(*),
1207      $                 xmin, xmax, ymin, ymax, zmin, zmax
1208
1209       integer k
1210       logical flag_x, flag_y, flag_z
1211       double precision x, y, z
1212       include 'constinterp.h.f'
1213       integer  isanan
1214       double precision db3val
1215       external isanan, returnananfortran, db3val
1216
1217       do k = 1, np
1218          x = xp(k)
1219          if ( xmin .le. x .and. x .le. xmax ) then
1220             flag_x = .true.
1221          else
1222             flag_x = .false.
1223          endif
1224          y = yp(k)
1225          if ( ymin .le. y .and. y .le. ymax ) then
1226             flag_y = .true.
1227          else
1228             flag_y = .false.
1229          endif
1230          z = zp(k)
1231          if ( zmin .le. z .and. z .le. zmax ) then
1232             flag_z = .true.
1233          else
1234             flag_z = .false.
1235          endif
1236
1237          if ( flag_x .and. flag_y .and. flag_z ) then
1238             fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz,
1239      $                     nx, ny, nz, kx, ky, kz, bcoef, work)
1240             dfpdx(k) = db3val(x, y, z, 1, 0, 0, tx, ty, tz,
1241      $                     nx, ny, nz, kx, ky, kz, bcoef, work)
1242             dfpdy(k) = db3val(x, y, z, 0, 1, 0, tx, ty, tz,
1243      $                     nx, ny, nz, kx, ky, kz, bcoef, work)
1244             dfpdz(k) = db3val(x, y, z, 0, 0, 1, tx, ty, tz,
1245      $                     nx, ny, nz, kx, ky, kz, bcoef, work)
1246
1247          elseif (outmode .eq. BY_NAN  .or.  isanan(x) .eq. 1
1248      $                                .or.  isanan(y) .eq. 1
1249      $                                .or.  isanan(z) .eq. 1) then
1250             CALL returnananfortran(fp(k))
1251             dfpdx(k) = fp(k)
1252             dfpdy(k) = fp(k)
1253             dfpdz(k) = fp(k)
1254
1255          elseif (outmode .eq. BY_ZERO) then
1256             fp(k) = 0.d0
1257             dfpdx(k) = 0.d0
1258             dfpdy(k) = 0.d0
1259             dfpdz(k) = 0.d0
1260
1261          elseif (outmode .eq. PERIODIC) then
1262             if (.not. flag_x) call proj_by_per(x, xmin, xmax)
1263             if (.not. flag_y) call proj_by_per(y, ymin, ymax)
1264             if (.not. flag_z) call proj_by_per(z, zmin, zmax)
1265             fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz, nx, ny, nz,
1266      $                     kx, ky, kz, bcoef, work)
1267             dfpdx(k) = db3val(x, y, z, 1, 0, 0, tx, ty, tz,
1268      $                     nx, ny, nz, kx, ky, kz, bcoef, work)
1269             dfpdy(k) = db3val(x, y, z, 0, 1, 0, tx, ty, tz,
1270      $                     nx, ny, nz, kx, ky, kz, bcoef, work)
1271             dfpdz(k) = db3val(x, y, z, 0, 0, 1, tx, ty, tz,
1272      $                     nx, ny, nz, kx, ky, kz, bcoef, work)
1273
1274          elseif (outmode .eq. C0) then
1275             if (.not. flag_x) call proj_on_grid(x, xmin, xmax)
1276             if (.not. flag_y) call proj_on_grid(y, ymin, ymax)
1277             if (.not. flag_z) call proj_on_grid(z, zmin, zmax)
1278             fp(k) = db3val(x, y, z, 0, 0, 0, tx, ty, tz, nx, ny, nz,
1279      $                     kx, ky, kz, bcoef, work)
1280             if ( flag_x ) then
1281                dfpdx(k) = 0.d0
1282             else
1283                dfpdx(k) = db3val(x, y, z, 1, 0, 0, tx, ty, tz,
1284      $                           nx, ny, nz, kx, ky, kz, bcoef, work)
1285             endif
1286             if ( flag_y ) then
1287                dfpdy(k) = 0.d0
1288             else
1289                dfpdy(k) = db3val(x, y, z, 0, 1, 0, tx, ty, tz,
1290      $                           nx, ny, nz, kx, ky, kz, bcoef, work)
1291             endif
1292             if ( flag_z ) then
1293                dfpdz(k) = 0.d0
1294             else
1295                dfpdz(k) = db3val(x, y, z, 0, 0, 1, tx, ty, tz,
1296      $                           nx, ny, nz, kx, ky, kz, bcoef, work)
1297             endif
1298
1299          endif
1300       enddo
1301
1302       end