interp.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE interp
!                       *****************
!
     &(xinit , yinit , ikinit , npinit , neinit ,
     & x , y , npoin , npmax , shp , elt)
!
!***********************************************************************
! PROGICIEL : STBTEL  V5.2        24/04/91    J-C GALLAND  (LNH)
!                               09/11/94    P LANG / TRIGRID (LHF)
!***********************************************************************
!
! FONCTION : INTERPOLATION DES FONDS SUR LE MAILLAGE
!
!----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________
! |      NOM       |MODE|                   ROLE
! |________________|____|______________________________________________
! |    X,Y         | -->|  COORDONNEES DES POINTS DU MAILLAGE
! |    ZF          |<-- |  COTES DU FOND
! |    XRELV,YRELV | -->|  COORDONNEES DES POINTS DE BATHY
! |    ZRELV       | -->|  COTES DES POINTS DE BATHY
! |    NBAT        | -->|  NOMBRE DE POINTS DE BATHY
! |    NBOR        | -->|  NUMEROTATION DES ELEMENTS DE BORD
! |    NPTFR       | -->|  NOMBRE DE POINTS FRONTIERE
! |    NFOND       | -->|  CANAUX DES FICHIERS DES FONDS
! |    NBFOND      | -->|  NOMBRE DE FICHIERS FONDS DONNES PAR
! |                |    |  L'UTILISATEUR (5 MAXI)
! |    FOND        | -->|  NOM DES FICHIERS DES FONDS
! |    DM          | -->|  DISTANCE MINIMALE A LA FRONTIERE
! |                |    |  POUR L'INTERPOLATION DES FONDS
! |    FONTRI      | -->|  INDICATEUR DE LECTURE DES FONDS DANS TRIGRID
! |    CORTRI      | -->|  CORRECTION DES FONDS POUR TRIGRID
! |                |    |
! | COMMON:        |    |
! |  GEO:          |    |
! |    MESH        | -->| TYPE DES ELEMENTS DU MAILLAGE
! |    NDP         | -->| NOMBRE DE NOEUDS PAR ELEMENTS
! |    NPOIN       | -->| NOMBRE TOTAL DE NOEUDS DU MAILLAGE
! |    NELEM       | -->| NOMBRE TOTAL D'ELEMENTS DU MAILLAGE
! |    NPMAX       | -->| DIMENSION EFFECTIVE DES TABLEAUX X ET Y
! |                |    | (NPMAX = NPOIN + 0.1*NELEM)
! |    NELMAX      | -->| DIMENSION EFFECTIVE DES TABLEAUX CONCERNANT
! |                |    | LES ELEMENTS (NELMAX = NELEM + 0.2*NELEM)
! |                |    |
! |________________|____|______________________________________________
! MODE : -->(DONNEE NON MODIFIEE), <--(RESULTAT), <-->(DONNEE MODIFIEE)
!----------------------------------------------------------------------
!
! APPELE PAR : STBTEL
! APPEL DE : LECFON, FASP
!
!**********************************************************************
!
      USE declarations_special
      IMPLICIT NONE
!
      INTEGER, INTENT(IN) :: NPINIT, NEINIT, NPOIN,NPMAX
      DOUBLE PRECISION, INTENT(IN) :: XINIT(npinit) , YINIT(npinit)
      INTEGER, INTENT(IN) :: IKINIT(neinit,3)
      INTEGER, INTENT(INOUT) :: ELT(npmax)
      DOUBLE PRECISION, INTENT(IN) :: X(npmax) , Y(npmax)
      DOUBLE PRECISION, INTENT(INOUT) :: SHP(npmax,3)
!
      INTEGER IELEM , JELEM , IPOIN
      DOUBLE PRECISION XP,YP,A1,A2,A3,C1,C2,X1,X2,X3,Y1,Y2,Y3
!
!=======================================================================
!
      WRITE(lu,4)
!
4     FORMAT(//,1x,'DATA INTERPOLATION',/,
     &          1x,'------------------',/)
!
      DO ipoin = 1,npoin
!
        xp = x(ipoin)
        yp = y(ipoin)
        c1 = -999999.d0
!
        DO ielem = 1,neinit
          x1 = xinit(ikinit(ielem,1))
          x2 = xinit(ikinit(ielem,2))
          x3 = xinit(ikinit(ielem,3))
          y1 = yinit(ikinit(ielem,1))
          y2 = yinit(ikinit(ielem,2))
          y3 = yinit(ikinit(ielem,3))
          a1 = (x3-x2)*(yp-y2) - (y3-y2)*(xp-x2)
          a2 = (x1-x3)*(yp-y3) - (y1-y3)*(xp-x3)
          a3 = (x2-x1)*(yp-y1) - (y2-y1)*(xp-x1)
          IF (a1.GE.0.AND.a2.GE.0.AND.a3.GE.0) GOTO 30
          c2 = min(a1,a2,a3) / ((x3-x2)*(y1-y2)-(y3-y2)*(x1-x2))
          IF (c2.GT.c1) THEN
            c1 = c2
            jelem = ielem
          ENDIF
        ENDDO
!
        WRITE(lu,*) 'EXTRAPOLATION REQUIRED FOR ',
     &              'THE NODE :',ipoin
        ielem = jelem
        x1 = xinit(ikinit(ielem,1))
        x2 = xinit(ikinit(ielem,2))
        x3 = xinit(ikinit(ielem,3))
        y1 = yinit(ikinit(ielem,1))
        y2 = yinit(ikinit(ielem,2))
        y3 = yinit(ikinit(ielem,3))
        a1 = (x3-x2)*(yp-y2) - (y3-y2)*(xp-x2)
        a2 = (x1-x3)*(yp-y3) - (y1-y3)*(xp-x3)
        a3 = (x2-x1)*(yp-y1) - (y2-y1)*(xp-x1)
!
30      CONTINUE
        c1 = (x3-x2)*(y1-y2)-(y3-y2)*(x1-x2)
        shp(ipoin,1) = a1/c1
        shp(ipoin,2) = a2/c1
        shp(ipoin,3) = a3/c1
        elt(ipoin) = ielem
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END