initial_drogues.f

Go to the documentation of this file.

!                         **********************
                          MODULE initial_drogues
!                         **********************
!
!***********************************************************************
! TELEMAC V8P0
!***********************************************************************
!
!brief    Module containing all subroutines to deal with drogues in
!+        general including randomly sampling parcels, etc.
!+        Module variables are mainly used for generalised input and
!+        output access.
!
!history  S.E. BOURBAN (HRW)
!+        21/08/2018
!+        V8P0
!+        Initial developments - inspired from ALGAE_TRANSP
!+        by A.Joly (EDF)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief_def, ONLY : bief_obj, bief_mesh, ipid,ncsize
      USE interface_parallel, ONLY : p_max,p_sum,p_max
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     SUBROUTINES MADE AVAILABLE
!      PUBLIC :: SAMPLE_WPOIN
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     MODULE VARIABLES
!     ( MAINLY USE FOR INPUT & OUTPUT AS A GENERALISATION OF THE
!       SPECIFCS VARIABLES INCLUDED IN MODULES DROGUES_* )
!
!     NUMBER OF DIFFERENT CLASSES OF DROGUES (ALL TYPES OF CLASSES)
      INTEGER :: ndrg_clss
!
!     INCREMENTAL HIGHEST TAG NUMBER FOR DROGUES
!     TAG NUMBERS ARE UNIQUE TO ALL PARCELS
      INTEGER :: ndrg_tags = 0
!
!     CLASS DEFINED FOR EACH NODE ON THE MESH
!     THE SIZE IS NELMAX, BOTH %I AND %R ARE ALLOCATED FOR FILE
!       EXCHANGES AND THESE ARE ALLOCATED IN POINT_*
      TYPE(bief_obj), TARGET :: nodclss
!
!     CLASS DEFINED FOR EACH PARCEL
!     THE SIZE IS NFLOT_MAX, BOTH %I AND %R ARE ALLOCATED FOR FILE
!       EXCHANGES AND THESE ARE ALLOCATED IN POINT
      TYPE(bief_obj), TARGET :: parclss
!
!
!     DROGUES PROPERTIES, DIFFERENT FOR DIFFERENT TYPES OF DROGUES
!     - DRG_DENSITY: COUNTING DROGUES PER SURFACE AREA
!     - DRG_RELEASE: TIME TO RELEASE IN SECOND
      DOUBLE PRECISION, ALLOCATABLE :: drg_density(:)
      DOUBLE PRECISION, ALLOCATABLE :: drg_release(:)
!
!------------------------------------------------------------------------
!
      SAVE
!
      CONTAINS
!
!=======================================================================
!
!       1) RANDOM SAMPLING
!
!
!                   ***********************
                    SUBROUTINE sample_wpoin
!                   ***********************
!
     &( np,np_max,ncls,ntag,xp,yp,tagp,clsp,eltp,shpp,dsty,
     &  npoin,nelem,nelmax,ikle,clsn,x,y )
!
!***********************************************************************
! TELEMAC V8P0
!***********************************************************************
!
!brief    Randomly sample parcels within mesh, around a set of points.
!+    The method goes through all points in the mesh where NODCLSS
!+      is not zero, and sample parts of the surrounding triangles.
!+    In order to be "fair" and to account for variable density and
!+      triangle areas, the sampling in first carried out on a sorted
!+      list (or cumulated required number of parcels per triangle).
!
!history  S.E. BOURBAN (HRW)
!+   21/08/2016
!+   V8P0
!+   Initial implementation
!
!history  M.S.TURNBULL (HRW)
!+   04/11/2019
!+   V8P2
!+   Corrections
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CLSN    |-->| CLASS AT NODES
!| CLSP    |<->| CLASS OF EACH DROGUE
!| DSTY    |-->| DROGUE DENSITY FOR EACH CLASS
!| ELTP    |<->| ELEMENT FOR EACH DROGUE
!| IKLE    |-->| CONNECTIVITY MATRIX
!| NELEM   |-->| NUMBER OF ELEMENTS IN IKLE
!| NELMAX  |-->| MAXIMUM NUMBER OF ELEMENTS IN IKLE
!| NCLS    |-->| NUMBER OF CLASSES
!| NP      |-->| NUMBER OF DROGUES TO BE CREATED
!| NPOIN   |-->| NUMBER OF NODES IN THE MESH
!| NP_MAX  |-->| MAXIMUM NUMBER OF DROGUES TO BE CREATED
!| NTAG    |<->| NUMBER OF TAGS
!| SHPP    |<->| BARYCENTRIC COORDINATES OF DROGUES
!| TAGP    |<->| TAG OF EACH DROGUE
!| X       |-->| ABSCISSAE OF POINTS IN THE MESH
!| XP      |<->| ABSCISSAE OF DROGUES
!| Y       |-->| ORDINATES OF POINTS IN THE MESH
!| YP      |-->| ORDINATES OF DROGUES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!***********************************************************************
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(INOUT)          :: NP,NTAG
      INTEGER, INTENT(IN)             :: NP_MAX,NELEM,NELMAX,NPOIN,NCLS
      INTEGER, INTENT(IN)             :: IKLE(nelmax,3),CLSN(npoin)
      INTEGER, INTENT(INOUT)          :: TAGP(np_max),CLSP(np_max)
      INTEGER, INTENT(INOUT)          :: ELTP(np_max)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin),Y(npoin),DSTY(ncls)
      DOUBLE PRECISION, INTENT(INOUT) :: XP(np_max),YP(np_max)
      DOUBLE PRECISION, INTENT(INOUT) :: SHPP(3,np_max)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          I,J,K,IEL,NEL,IP,NJ, K1,K2,KI,N1,N2
      INTEGER          JPID,MT,NPI,NPP, ITAG,ICLS
      DOUBLE PRECISION R0,LENGTH,A
      DOUBLE PRECISION X1,Y1,X2,Y2,X3,Y3, DX1,DY1,DX2,DY2,DX3,DY3
!
      INTEGER, ALLOCATABLE          :: PPKLE(:),PELEM(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: AKLE(:)
!
!-----------------------------------------------------------------------
!
!     SET THE VARIABLE DROGUES' CLASS
      mt = 0
      DO i = 1,npoin
        mt = max( mt,clsn(i) )
      ENDDO
      IF( ncsize.GT.1 ) mt = p_max( mt )
      IF( mt.NE.ncls ) THEN
        WRITE(lu,22) mt,ncls
 22     FORMAT(1x,'CONDIN_DROGUES:',/,
     &    1x,'     NUMBER OF CLASSES OF DROGUES READ FROM',/,
     &    1x,'     THE GEOMETRY FILE ',i8,' DIFFERENT FROM THE',/,
     &    1x,'     NUMBER SET IN THE CAS FILE ',i8)
        IF( mt.GT.ncls ) THEN
          CALL plante(1)
          stop
        ENDIF
      ENDIF
!
!     LOCAL MEMORY ALLOCATION - NUMBER OF TRIANGLES CONTAINING PARCELS
!
      nel = 0
      DO j = 1,nelem
        DO k = 1,3
          i = ikle( j,k )
          IF( clsn(i).GT.0 ) nel = nel + 1
        ENDDO
      ENDDO
      ALLOCATE( pelem(2,nel) )
      ALLOCATE( ppkle(2*nel) )
      ALLOCATE( akle(2*nel) )
!
!     EVALUATE TOTAL NUMBER OF PARCELS REQUIRED FOR ONE PROCESSOR
!     BY COMPUTING THE CUMULATED SURFACE AREA WEIGHTED BY TARGET
!     PARCEL DENSITY
!
      iel = 0
      a = 0.d0
      DO j = 1,nelem
        DO k = 1,3
          i = ikle( j,k )
          IF( clsn(i).GT.0 ) THEN
            iel = iel + 1
            pelem(1,iel) = j
            pelem(2,iel) = k
!
            n1 = ikle( j,mod(k,3)+1 )
            n2 = ikle( j,mod(k+1,3)+1 )
!
            dx1 = ( x(i)-x(n1) )/2.d0
            dx2 = ( x(i)-x(n2) )/2.d0
            dy1 = ( y(i)-y(n1) )/2.d0
            dy2 = ( y(i)-y(n2) )/2.d0
            dx3 = ( 2.d0*x(i)-x(n1)-x(n2) )/3.d0
            dy3 = ( 2.d0*y(i)-y(n1)-y(n2) )/3.d0
!
            akle(2*iel-1) = a +
     &        dsty(clsn(i)) * ( dx1*dy3-dx3*dy1 )/2.d0
            a = akle(2*iel-1)
            ppkle(2*iel-1) = 0
!
            akle(2*iel) = a +
     &        dsty(clsn(i)) * ( dx3*dy2-dx2*dy3 )/2.d0
            a = akle(2*iel)
            ppkle(2*iel) = 0
!
          ENDIF
        ENDDO
      ENDDO
!
!     TOTAL NUMBER OF PARCELS FOR CURRENT PROCESSOR
      npi = int(a)
      npp = npi + np
      IF( ncsize.GT.1 ) npp = p_sum( npp )
!     CHECK AGAINST NUMBER OF PARCELS POSSIBLE
      IF( npp.GT.np_max ) THEN
        WRITE(lu,32) npp,np_max
 32     FORMAT(1x,'DROGUES::SAMPLE_WPOIN:',/,
     &    1x,'     REQUIRED NUMBER OF DROGUES (',i8,')',/,
     &    1x,'     LARGER THAN THE MAXIMUM NUMBER OF DROGUES',/,
     &    1x,'     POSSIBLE (',i8,')',/,
     &    1x,'     INCREASE THE MAXIMUM OR REDUCE YOUR DENSITY.')
        CALL plante(1)
        stop
      ENDIF
      IF (nel.GT.0) length = akle(2*nel)
!
!         ALLOCATING UNIQUE TAG NUMBERS FOR EACH PROCESSOR
!
      IF( ncsize.GT.1 ) THEN
        npp = 0
        DO jpid = 0,ncsize-1
          IF( jpid.EQ.ipid ) THEN
            npp = npp + npi
            itag = npp
          ENDIF
          npp = p_max( npp )
        ENDDO
        itag = itag - npi
      ELSE
        itag = 0
      ENDIF
!
!-----------------------------------------------------------------------
!
!     WEIGHTED SAMPLING BY POINTS (SUB-ELEMENT TRIANGLES)
!
      DO ip = 1,npi
!       RANDOM SAMPLE
        CALL random_number(r0)
        r0 = r0 * length
        IF( r0.LT.akle(1) ) THEN
!         SPECIAL CASE FOR 1ST QUADRILATERAL
          ppkle(1) = ppkle(1) + 1
        ELSE
!         DICHOTOMIE TO FIND WHERE THE PARCEL FALLS
          k1 = 1
          k2 = 2*nel
          DO WHILE( .true. )
            ki = int((k1+k2)/2)
            IF( r0.GT.akle(ki) ) THEN
              k1 = ki
            ELSE
              k2 = ki
            ENDIF
            IF( k1+1.EQ.k2 ) EXIT
          ENDDO
          ppkle(k2) = ppkle(k2) + 1
        ENDIF
!
      ENDDO
!
!-----------------------------------------------------------------------
!
!     SAMPLING AROUND WEIGHTED POINTS, ONE TRIANGLE AT A TIME
      ip = np
      DO iel = 1,nel
!
        j = pelem(1,iel)
        k = pelem(2,iel)
        i = ikle( j,k )
        icls = clsn(i)
!
        n1 = ikle( j,mod(k,3)+1 )
        n2 = ikle( j,mod(k+1,3)+1 )
!
        x1 = ( x(i)+x(n1) )/2.d0
        y1 = ( y(i)+y(n1) )/2.d0
        x2 = ( x(i)+x(n2) )/2.d0
        y2 = ( y(i)+y(n2) )/2.d0
        x3 = ( x(i)+x(n1)+x(n2) )/3.d0
        y3 = ( y(i)+y(n1)+y(n2) )/3.d0
!
!       SAMPLE FIRST TRIANGLE
        nj = ppkle(2*iel-1)
        CALL sample_triangle( ip+1,nj,np_max,j,itag,icls,
     &        tagp,clsp,eltp,shpp,
     &        xp,yp, x(i),y(i),x1,y1,x3,y3 )
        ip = ip + nj
!
!       SAMPLE SECOND TRIANGLE
        nj = ppkle(2*iel)
        CALL sample_triangle( ip+1,nj,np_max,j,itag,icls,
     &        tagp,clsp,eltp,shpp,
     &        xp,yp, x(i),y(i),x3,y3,x2,y2 )
        ip = ip + nj
!
      ENDDO
      IF( ncsize.GT.1 ) itag = p_max( itag )
      ntag = ntag + itag
      np = np + npi
!
!-----------------------------------------------------------------------
!
!     LOCAL MEMORY DE-ALLOCATION
!
      DEALLOCATE( pelem )
      DEALLOCATE( ppkle )
      DEALLOCATE( akle )
!
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE sample_wpoin
!
!                 **************************
                  SUBROUTINE sample_triangle
!                 **************************
!
     &( ip,np,np_max,ielm,itag,icls,tagp,clsp,eltp,shpp,
     &  xp,yp,x1,y1,x2,y2,x3,y3 )
!
!***********************************************************************
! TELEMAC V8P0
!***********************************************************************
!
!brief    Randomly sample parcels within a triangle.
!+    In order to be "fair" from a placement view point, R1 and R2
!+      are generated along the two edges P0P1 and P0P2 respectively,
!+      therefore covering a parallelogram with parcels (assuming
!+      the lengths of the edges P0P1 and P0P2 are not too different).
!+    Half of the parallelogram is then moved back into the triangle
!+      by symmetry with the mid-point between P1 and P2.
!
!history  S.E. BOURBAN (HRW)
!+   21/08/2016
!+   V8P0
!+   Initial implementation
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CLSP    |<->| CLASS OF EACH DROGUE
!| ELTP    |<->| ELEMENT FOR EACH DROGUE
!| ICLS    |<->| CLASS TO BE ASSIGNED
!| IELM    |-->| ELEMENT TO BE ASSIGNED
!| IP      |-->| INDEX NUMBER OF FIRST DROGUE
!| ITAG    |-->| TAGGING RANGE FOR DROGUES PER PROCESSOR
!| NP      |-->| NUMBER OF DROGUES TO BE CREATED
!| NP_MAX  |-->| MAXIMUM NUMBER OF DROGUES TO BE CREATED
!| SHPP    |<->| BARYCENTRIC COORDINATES OF DROGUES
!| TAGP    |<->| TAG OF EACH DROGUE
!| XN,YN   |-->| WHERE N = 1,2,3, COORDINATES OF THE TRIANGLE WITHIN
!|         |---|   WHICH THE NEW NP DROGUES WILL BE RANDOMLY SELECTED
!| XP,YP   |<--| COORDINATES OF THE NP NEW (RANDOMLY SAMPLED) DROGUES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!***********************************************************************
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: IP,NP,NP_MAX, IELM
      INTEGER, INTENT(INOUT)          :: TAGP(np_max),CLSP(np_max)
      INTEGER, INTENT(INOUT)          :: ELTP(np_max),ITAG,ICLS
      DOUBLE PRECISION, INTENT(INOUT) :: XP(np_max),YP(np_max)
      DOUBLE PRECISION, INTENT(IN)    :: X1,Y1,X2,Y2,X3,Y3
      DOUBLE PRECISION, INTENT(INOUT) :: SHPP(3,np_max)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          I
      DOUBLE PRECISION R1,R2, SURDET, X0,Y0
!
!-----------------------------------------------------------------------
!
      DO i = ip,ip+np-1
!
        CALL random_number(r1)
        CALL random_number(r2)
        xp(i) = x1 + r1*( x2-x1 ) + r2*( x3-x1 )
        yp(i) = y1 + r1*( y2-y1 ) + r2*( y3-y1 )
        IF( r2 .GT. ( 1.d0-r1 ) ) THEN
          xp(i) = ( x3+x2 ) - xp(i)
          yp(i) = ( y3+y2 ) - yp(i)
        ENDIF
!
        x0 = xp(i)
        y0 = yp(i)
        eltp(i) = ielm
        shpp(1,i) = ( x3-x2 )*( y0-y2 ) - ( y3-y2 )*( x0-x2 )
        shpp(2,i) = ( x1-x3 )*( y0-y3 ) - ( y1-y3 )*( x0-x3 )
        shpp(3,i) = ( x2-x1 )*( y0-y1 ) - ( y2-y1 )*( x0-x1 )
        surdet = 1.d0 / ( ( x2-x1 )*( y3-y1 ) - ( x3-x1 )*( y2-y1 ) )
        shpp(1,i) = shpp(1,i) * surdet
        shpp(2,i) = shpp(2,i) * surdet
        shpp(3,i) = shpp(3,i) * surdet
        itag = itag + 1
        tagp(i) = itag
        clsp(i) = icls
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE sample_triangle
!
!                 **************************
                  SUBROUTINE sample_polyline
!                 **************************
!
     &( np,np_max,ncls,ntag, xp,yp,tagp,clsp,eltp,shpp,dsty,
     &  ny,iy,vy,ng,xg,yg, npoin,nelem,nelmax,ikle,x,y )
!
!***********************************************************************
! TELEMAC V8P0
!***********************************************************************
!
!brief    Randomly sample parcels within a polygon defined by (XG,YG).
!+    The sampling method generate a parcel location within the bounds
!+      of the polygon, and then check if the point location is within
!+      the polygon.
!+    A better way to do it could be to use a tessalation of the
!+      polygon and then call on SAMPLE_WELEM
!
!history  S.E. BOURBAN (HRW)
!+   21/08/2016
!+   V8P0
!+   Initial implementation
!
!history  M.S.TURNBULL (HRW)
!+   21/11/2019
!+   V8P2
!+   Corrections
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CLSP    |<->| CLASS OF EACH DROGUE
!| DSTY    |-->| DROGUE DENSITY FOR EACH CLASS
!| ELTP    |<->| ELEMENT FOR EACH DROGUE
!| IKLE    |-->| CONNECTIVITY MATRIX
!| IY      |-->| ARRAYS OF NUMBER OF POINT PER POLYGONS
!| NCLS    |-->| NUMBER OF CLASSES
!| NELEM   |-->| NUMBER OF ELEMENTS IN IKLE
!| NELMAX  |-->| MAXIMUM NUMBER OF ELEMENTS IN IKLE
!| NG      |-->| NUMBER OF POLYGON VERTICES
!| NP      |-->| NUMBER OF DROGUES TO BE CREATED
!| NPOIN   |-->| NUMBER OF NODES IN THE MESH
!| NP_MAX  |-->| MAXIMUM NUMBER OF DROGUES TO BE CREATED
!| NTAG    |<->| NUMBER OF TAGS
!| NY      |-->| NUMBER OF POLYGONS
!| SHPP    |<->| BARYCENTRIC COORDINATES OF DROGUES
!| TAGP    |<->| TAG OF EACH DROGUE
!| X       |-->| ABSCISSAE OF POINTS IN THE MESH
!| XG,YG   |-->| COORDINATES OF THE VERTICES OF THE POLYGON WITHIN
!|         |---|   WHICH THE NEW NP DROGUESS WILL BE RANDOMLY SELECTED
!| XP,YP   |<--| COORDINATES OF THE NP NEW (RANDOMLY SAMPLED) DROGUES
!| VY      |-->| ARRAYS OF VALUES PER POLYGONS
!| Y       |-->| ORDINATES OF POINTS IN THE MESH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!***********************************************************************
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(INOUT)          :: NP,NTAG
      INTEGER, INTENT(IN)             :: NP_MAX,NELEM,NELMAX,NPOIN,NCLS
      INTEGER, INTENT(IN)             :: IKLE(nelmax,3)
      INTEGER, INTENT(INOUT)          :: TAGP(np_max),CLSP(np_max)
      INTEGER, INTENT(INOUT)          :: ELTP(np_max)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin),Y(npoin),DSTY(ncls)
      DOUBLE PRECISION, INTENT(INOUT) :: XP(np_max),YP(np_max)
      DOUBLE PRECISION, INTENT(INOUT) :: SHPP(3,np_max)
      INTEGER, INTENT(IN)             :: NY,IY(ny), NG
      DOUBLE PRECISION, INTENT(IN)    :: VY(ny), XG(ng),YG(ng)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          I,J,K,IEL,NEL, N1,N2,N3, NPP, JPID, MT
      INTEGER          NG1,NG2, JY
      LOGICAL          FOUND
      DOUBLE PRECISION A, R1,R2, XA,YA,XI,YI, X0,Y0,X1,Y1,X2,Y2,X3,Y3
      DOUBLE PRECISION SURDET, DET1,DET2,DET3
      DOUBLE PRECISION, PARAMETER :: CHOUIA = 1.d-10
!
      INTEGER, ALLOCATABLE        :: PPKLE(:),PELEM(:)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      LOGICAL INPOLY
      EXTERNAL inpoly
!
!-----------------------------------------------------------------------
!
!     SET THE VARIABLE DROGUES' CLASS
      mt = 0
      DO i = 1,ny
        mt = max( mt,int(vy(i)) )
      ENDDO
      IF( ncsize.GT.1 ) mt = p_max( mt )
      IF( mt.NE.ncls ) THEN
        WRITE(lu,22) mt,ncls
 22     FORMAT(1x,'CONDIN_DROGUES:',/,
     &    1x,'     NUMBER OF CLASSES OF DROGUES READ FROM',/,
     &    1x,'     THE POLYGON FILE ',i8,' DIFFERENT FROM THE',/,
     &    1x,'     NUMBER SET IN THE CAS FILE ',i8)
        IF( mt.GT.ncls ) THEN
          CALL plante(1)
          stop
        ENDIF
      ENDIF
!
      ALLOCATE( ppkle(nelem) )
!
!-----------------------------------------------------------------------
!
      ng1 = 1
      DO jy = 1,ny
!       NUMBER OF POINTS IN THE JY POLYGON
        ng2 = int( iy(jy) )
!
!-----------------------------------------------------------------------
!
!     EVALUATING THE NUMBER OF PARCELS WITHIN THE POLYGONE
!     (ACROSS ALL PROCESSORS)
!
!       SIZE AND AREA OF THE POLYGON
        a = ( xg(ng2)*yg(ng1) - xg(ng1)*yg(ng2) )
        xa = xg(ng2)
        ya = yg(ng2)
        xi = xg(ng2)
        yi = yg(ng2)
        DO i = ng1,ng2-1
          a = a + ( xg(i)*yg(i+1) - xg(i+1)*yg(i) )
          xa = max( xa,xg(i) )
          ya = max( ya,yg(i) )
          xi = min( xi,xg(i) )
          yi = min( yi,yg(i) )
        ENDDO
        a = a / 2.d0
!
!       NUMBER OF PARCELS GIVEN A SPATIAL DENSITY
        npp = int(a*dsty(int(vy(jy)))+1.d0)
!       CHECK AGAINST NUMBER OF PARCELS POSSIBLE
        IF( (np+npp).GT.np_max ) THEN
          WRITE(lu,33) (np+npp),np_max
 33       FORMAT(1x,'DROGUES::SAMPLE_POLYLINE:',/,
     &      1x,'     REQUIRED NUMBER OF DROGUES (',i8,')',/,
     &      1x,'     LARGER THAN THE MAXIMUM NUMBER OF DROGUES',/,
     &      1x,'     POSSIBLE (',i8,')',/,
     &      1x,'     INCREASE THE MAXIMUM OR REDUCE YOUR DENSITY.')
          CALL plante(1)
          stop
        ENDIF
!
!-----------------------------------------------------------------------
!
!     IDENTIFYING THOSE ELEMENTS COVERED BY THE POLYGON
!     (FOR EACH PROCESSOR)
!
!       LOCAL MEMORY ALLOCATION - NUMBER OF TRIANGLES COVERED
!
        nel = 0
        DO j = 1,nelem
          found = .false.
          ppkle(j) = 0
          DO k = 1,3
            i = ikle( j,k )
            found = found .OR.
     &        inpoly( x(i),y(i), xg(ng1:ng2),yg(ng1:ng2),ng2-ng1+1 )
          ENDDO
          IF( found ) THEN
            nel = nel + 1
            ppkle(j) = 1
          ENDIF
        ENDDO
!
        ALLOCATE( pelem(nel) )
!
        iel = 0
        DO j = 1,nelem
          IF( ppkle(j).NE.0 ) THEN
            iel = iel + 1
            pelem(iel) = j
          ENDIF
        ENDDO
!
!-----------------------------------------------------------------------
!
!     SAMPLING THE POLYGON USING A HIT AND MISS STRATEGY
!
        DO i = 1,npp
!
!       SAMPLING OVER THE ENTIRE POLYGON
          found = .false.
          DO WHILE( .NOT.found )
            IF (ncsize.GT.1) THEN
              x0 = -1.d20
              y0 = -1.d20
              IF (ipid.EQ.0) THEN
                CALL random_number(r1)
                CALL random_number(r2)
                x0 = xi + r1*( xa-xi )
                y0 = yi + r2*( ya-yi )
              ENDIF
              x0 = p_max(x0)
              y0 = p_max(y0)
            ELSE
              CALL random_number(r1)
              CALL random_number(r2)
              x0 = xi + r1*( xa-xi )
              y0 = yi + r2*( ya-yi )
            ENDIF
            found = inpoly( x0,y0, xg(ng1:ng2),yg(ng1:ng2),ng2-ng1+1 )
          ENDDO
!
!       IDENTIFYING THE TRIANGLE WHERE THE PARCELS FELL IN
!       (FOR EACH PROCESSOR, AND IT MAY BE ON THE EDGE OF TWO OR MORE)
!
          jpid = -1
          DO iel = 1,nel
!
            j = pelem(iel)
            n1 = ikle( j,1 )
            n2 = ikle( j,2 )
            n3 = ikle( j,3 )
            x1 = x(n1)
            y1 = y(n1)
            x2 = x(n2)
            y2 = y(n2)
            x3 = x(n3)
            y3 = y(n3)
!
            det1 = ( x3-x2 )*( y0-y2 ) - ( y3-y2 )*( x0-x2 )
            det2 = ( x1-x3 )*( y0-y3 ) - ( y1-y3 )*( x0-x3 )
            det3 = ( x2-x1 )*( y0-y1 ) - ( y2-y1 )*( x0-x1 )
            IF( det1.GE.-chouia .AND.
     &          det2.GE.-chouia .AND.
     &          det3.GE.-chouia ) THEN
              jpid = ipid
              EXIT
            ENDIF
!
          ENDDO
!
          IF( ncsize.GT.1 ) THEN
            jpid = p_max(jpid)
            ntag = p_max(ntag)
          ENDIF
          IF( jpid.EQ.ipid ) THEN
            np = np + 1
            xp(np) = x0
            yp(np) = y0
            eltp(np) = j
            surdet = 1.d0 / ( (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1) )
            shpp(1,np) = det1*surdet
            shpp(2,np) = det2*surdet
            shpp(3,np) = det3*surdet
            ntag = ntag + 1
            tagp(np) = ntag
            clsp(np) = int( vy(jy) )
          ENDIF
!
        ENDDO
!
        ng1 = ng2 + 1
        DEALLOCATE( pelem )
!
      ENDDO
!
      DEALLOCATE( ppkle )
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE sample_polyline
!
!                 ************************
                  SUBROUTINE sample_points
!                 ************************
!
     &( np,np_max,ncls,ntag, xp,yp,tagp,clsp,eltp,shpp, ng,xg,yg,vg,
     &  npoin,nelem,nelmax,ikle,x,y )
!
!***********************************************************************
! TELEMAC V8P0
!***********************************************************************
!
!brief    Use the X-Y-Z provided to place particles (XG,YG).
!
!history  S.E. BOURBAN (HRW)
!+   21/08/2016
!+   V8P0
!+   Initial implementation
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CLSP    |<->| CLASS OF EACH DROGUE
!| ELTP    |<->| ELEMENT FOR EACH DROGUE
!| IKLE    |-->| CONNECTIVITY MATRIX
!| NCLS    |-->| NUMBER OF CLASSES
!| NELEM   |-->| NUMBER OF ELEMENTS IN IKLE
!| NELMAX  |-->| MAXIMUM NUMBER OF ELEMENTS IN IKLE
!| NG      |-->| TOTAL NUMBER OF POINTS
!| NP      |-->| NUMBER OF DROGUES TO BE CREATED
!| NPOIN   |-->| NUMBER OF NODES IN THE MESH
!| NP_MAX  |-->| MAXIMUM NUMBER OF DROGUES TO BE CREATED
!| NTAG    |<->| NUMBER OF TAGS
!| SHPP    |<->| BARYCENTRIC COORDINATES OF DROGUES
!| TAGP    |<->| TAG OF EACH DROGUE
!| X       |-->| ABSCISSAE OF POINTS IN THE MESH
!| XG      |-->| X-COORDINATE FOR ALL POINTS
!| YG      |-->| Y-COORDINATE FOR ALL POINTS
!| XP,YP   |<--| COORDINATES OF THE NP NEW (RANDOMLY SAMPLED) DROGUES
!| VG      |-->| ARRAYS OF VALUES PER POINT
!| Y       |-->| ORDINATES OF POINTS IN THE MESH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!***********************************************************************
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(INOUT)          :: NP,NTAG
      INTEGER, INTENT(IN)             :: NP_MAX,NELEM,NELMAX,NPOIN,NCLS
      INTEGER, INTENT(IN)             :: IKLE(nelmax,3)
      INTEGER, INTENT(INOUT)          :: TAGP(np_max),CLSP(np_max)
      INTEGER, INTENT(INOUT)          :: ELTP(np_max)
      DOUBLE PRECISION, INTENT(IN)    :: X(npoin),Y(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: XP(np_max),YP(np_max)
      DOUBLE PRECISION, INTENT(INOUT) :: SHPP(3,np_max)
      INTEGER, INTENT(IN)             :: NG
      DOUBLE PRECISION, INTENT(IN)    :: XG(ng),YG(ng),VG(ng)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          I,J, N1,N2,N3, JPID, MT
      DOUBLE PRECISION X0,Y0,X1,Y1,X2,Y2,X3,Y3
      DOUBLE PRECISION SURDET, DET1,DET2,DET3
      DOUBLE PRECISION, PARAMETER :: CHOUIA = 1.d-10
!
!-----------------------------------------------------------------------
!
!     SET THE VARIABLE DROGUES' CLASS
!
      mt = 0
      DO i = 1,ng
        mt = max( mt,int(vg(i)) )
      ENDDO
      IF( ncsize.GT.1 ) mt = p_max( mt )
      IF( mt.NE.ncls ) THEN
        WRITE(lu,22) mt,ncls
 22     FORMAT(1x,'CONDIN_DROGUES:',/,
     &    1x,'     NUMBER OF CLASSES OF DROGUES READ FROM',/,
     &    1x,'     THE X-Y-Z FILE ',i8,' DIFFERENT FROM THE',/,
     &    1x,'     NUMBER SET IN THE CAS FILE ',i8)
        IF( mt.GT.ncls ) THEN
          CALL plante(1)
          stop
        ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
!
!     NUMBER OF PARCELS GIVEN A SPATIAL DENSITY
!
!     CHECK AGAINST NUMBER OF PARCELS POSSIBLE
      IF( (np+ng).GT.np_max ) THEN
        WRITE(lu,33) (np+ng),np_max
 33     FORMAT(1x,'DROGUES::SAMPLE_POLYLINE:',/,
     &    1x,'     REQUIRED NUMBER OF DROGUES (',i8,')',/,
     &    1x,'     LARGER THAN THE MAXIMUM NUMBER OF DROGUES',/,
     &    1x,'     POSSIBLE (',i8,')',/,
     &    1x,'     INCREASE THE MAXIMUM OR REDUCE YOUR DENSITY.')
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!     POSITIONING THE PROVIDED POINTS
!
      DO i = 1,ng
!
        x0 = xg(i)
        y0 = yg(i)
!
!       IDENTIFYING THE TRIANGLE WHERE THE PARCELS FELL IN
!       (FOR EACH PROCESSOR, AND IT MAY BE ON THE EDGE OF TWO OR MORE)
!
        jpid = -1
        DO j = 1,nelem
!
          n1 = ikle( j,1 )
          n2 = ikle( j,2 )
          n3 = ikle( j,3 )
          x1 = x(n1)
          y1 = y(n1)
          x2 = x(n2)
          y2 = y(n2)
          x3 = x(n3)
          y3 = y(n3)
!
          det1 = ( x3-x2 )*( y0-y2 ) - ( y3-y2 )*( x0-x2 )
          det2 = ( x1-x3 )*( y0-y3 ) - ( y1-y3 )*( x0-x3 )
          det3 = ( x2-x1 )*( y0-y1 ) - ( y2-y1 )*( x0-x1 )
          IF( det1.GE.-chouia .AND.
     &        det2.GE.-chouia .AND.
     &        det3.GE.-chouia ) THEN
            jpid = ipid
            EXIT
          ENDIF
!
        ENDDO
!
        IF( ncsize.GT.1 ) THEN
          jpid = p_max(jpid)
          ntag = p_max(ntag)
        ENDIF
        IF( jpid.EQ.ipid ) THEN
          np = np + 1
          xp(np) = x0
          yp(np) = y0
          eltp(np) = j
          surdet = 1.d0 / ( (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1) )
          shpp(1,np) = det1*surdet
          shpp(2,np) = det2*surdet
          shpp(3,np) = det3*surdet
          ntag = ntag + 1
          tagp(np) = ntag
          clsp(np) = int( vg(np) )
        ENDIF
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE sample_points
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      END MODULE initial_drogues