1 | include(dom.inc)
   2 | 
   3 |       SUBROUTINE interp3d(icount,x,y,z,F,jcount,xn,yn,zn,Fn,nmax)
   4 | 
   5 |         IMPLICIT NONE
   6 | 
   7 | !       IN
   8 |         DOM_INT                         :: icount, jcount, nmax
   9 |         DOM_REAL, DIMENSION (nmax)      :: xn,yn,zn,Fn
  10 | 
  11 | !       OUT
  12 |         DOM_REAL, DIMENSION (nmax)      :: x,y,z,F
  13 | 
  14 | !       LOCAL
  15 | 
  16 |         DOM_INT,  PARAMETER             :: degre=8
  17 | 
  18 |         DOM_REAL                        :: SUMALPHAH, SUMALPHA
  19 |         DOM_REAL                        :: ATEMP
  20 |         DOM_REAL                        :: ITEMP
  21 |         DOM_REAL                        :: xi, yi, zi
  22 |         DOM_REAL, PARAMETER             :: xmax=1000000.
  23 |         DOM_REAL, DIMENSION (nmax)      :: L
  24 |         DOM_REAL, DIMENSION (0:degre+1) :: VALPTS
  25 | 
  26 |         DOM_INT                         :: i, j, k
  27 |         DOM_INT,  DIMENSION (0:degre+1) :: PTS
  28 | 
  29 | 
  30 |         DO i=1,icount
  31 | 
  32 |         xi=x(i)
  33 |         yi=y(i)
  34 |         zi=z(i)
  35 | 
  36 |            DO k=0,degre+1
  37 | 
  38 |               PTS(k)=nmax
  39 | 
  40 |               IF ((k.eq.0).or.(k.eq.degre+1)) THEN
  41 |                  VALPTS(k)=0.
  42 |               ELSE
  43 |                  VALPTS(k)=xmax
  44 |               ENDIF
  45 | 
  46 |            ENDDO
  47 | 
  48 |            DO j=1,jcount
  49 | 
  50 |               L(j)=sqrt((xi-xn(j))**2+(yi-yn(j))**2+(zi-zn(j))**2)
  51 | 
  52 |               k=1
  53 | 
  54 |               DO WHILE (L(j)<VALPTS(k))
  55 |                  ITEMP=PTS(k)
  56 |                  ATEMP=VALPTS(k)
  57 |                  PTS(k)=j
  58 |                  VALPTS(k)=L(j)
  59 |                  PTS(k-1)=ITEMP
  60 |                  VALPTS(k-1)=ATEMP
  61 |                  k=k+1
  62 |               ENDDO
  63 | 
  64 |            ENDDO
  65 | 
  66 |            SUMALPHAH=0.
  67 |            SUMALPHA=0.
  68 |            DO j=1,degre
  69 | 
  70 |               SUMALPHAH=SUMALPHAH+Fn(PTS(j))/L(PTS(j))**2
  71 |               SUMALPHA=SUMALPHA+1/L(PTS(j))**2
  72 | 
  73 |            ENDDO
  74 | 
  75 |            F(i)=SUMALPHAH/SUMALPHA
  76 | 
  77 |         ENDDO
  78 | 
  79 |       END SUBROUTINE interp3d