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 |               IF(L(j).eq.0) L(j)=1
  52 | 
  53 |               k=1
  54 | 
  55 |               DO WHILE (L(j)<VALPTS(k))
  56 |                  ITEMP=PTS(k)
  57 |                  ATEMP=VALPTS(k)
  58 |                  PTS(k)=j
  59 |                  VALPTS(k)=L(j)
  60 |                  PTS(k-1)=ITEMP
  61 |                  VALPTS(k-1)=ATEMP
  62 |                  k=k+1
  63 |               ENDDO
  64 | 
  65 |            ENDDO
  66 | 
  67 |            SUMALPHAH=0.
  68 |            SUMALPHA=0.
  69 |            DO j=1,degre
  70 | 
  71 |               SUMALPHAH=SUMALPHAH+Fn(PTS(j))/L(PTS(j))**2
  72 |               SUMALPHA=SUMALPHA+1/L(PTS(j))**2
  73 | 
  74 |            ENDDO
  75 | 
  76 |            F(i)=SUMALPHAH/SUMALPHA
  77 | 
  78 |         ENDDO
  79 | 
  80 |       END SUBROUTINE interp3d